Complex derivative of power fails

[jagot]

I'm trying to calculate complex derivatives of analytic functions, using a trick by @tpapp:

using ForwardDiff

∂(f, x::Real) = ForwardDiff.derivative(f, x)
function ∂(f, z::Complex)
    # https://discourse.julialang.org/t/automatic-differentiation-of-complex-valued-functions/30263/3
    ff = ((x,y),) -> begin
        fz = f(complex(x,y))
        [real(fz), imag(fz)]
    end
    J = ForwardDiff.jacobian(ff, [real(z), imag(z)])
    # Complex derivative for an analytic function
    #   f(x + im*y) = u(x,y) + im*v(x,y)
    # is given by
    #   f′(x + im*y) = uₓ + im*vₓ = v_y - im*u_y
    J[1,1] + im*J[2,1]
end

This works for some simple functions:

julia> ∂(sin, 0.1), ∂(sin, 0.1+0im), ∂(sin, 0.1+0.3im)
(0.9950041652780258, 0.9950041652780258 + 0.0im, 1.0401161756837587 - 0.030401301333122962im)

julia> ∂(t -> exp(-t), 0.1+0.3im)
-0.8644242021759518 + 0.26739774077289963im

But fails for functions involving powers:

julia> f = z -> z^2
#4019 (generic function with 1 method)

julia> ∂(f, 0.1)
0.2

julia> ∂(f, 0.1+0.0im)
ERROR: MethodError: no method matching Int64(::ForwardDiff.Dual{ForwardDiff.Tag{var"#4015#4016"{var"#4019#4020"}, Float64}, Float64, 2})
Closest candidates are:
  (::Type{T})(::T) where T<:Number at boot.jl:760
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
  (::Type{T})(::BigInt) where T<:Union{Int128, Int16, Int32, Int64, Int8} at gmp.jl:356
  ...
Stacktrace:
  [1] convert(#unused#::Type{Int64}, x::ForwardDiff.Dual{ForwardDiff.Tag{var"#4015#4016"{var"#4019#4020"}, Float64}, Float64, 2})
    @ Base ./number.jl:7
  [2] _cpow(z::Complex{ForwardDiff.Dual{ForwardDiff.Tag{var"#4015#4016"{var"#4019#4020"}, Float64}, Float64, 2}}, p::Complex{ForwardDiff.Dual{ForwardDiff.Tag{var"#4015#4016"{var"#4019#4020"}, Float64}, Float64, 2}})
    @ Base ./complex.jl:740
  [3] ^
    @ ./complex.jl:808 [inlined]
  [4] ^
    @ ./promotion.jl:355 [inlined]
  [5] ^
    @ ./complex.jl:813 [inlined]
  [6] macro expansion
    @ ./none:0 [inlined]
  [7] literal_pow
    @ ./none:0 [inlined]
  [8] #4019
    @ ./REPL[124]:1 [inlined]
  [9] (::var"#4015#4016"{var"#4019#4020"})(::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#4015#4016"{var"#4019#4020"}, Float64}, Float64, 2}})
    @ Main /tmp/test.jl:6
 [10] vector_mode_dual_eval
    @ ~/.julia/packages/ForwardDiff/sqhTO/src/apiutils.jl:37 [inlined]
 [11] vector_mode_jacobian(f::var"#4015#4016"{var"#4019#4020"}, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#4015#4016"{var"#4019#4020"}, Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#4015#4016"{var"#4019#4020"}, Float64}, Float64, 2}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/sqhTO/src/jacobian.jl:147
 [12] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#4015#4016"{var"#4019#4020"}, Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#4015#4016"{var"#4019#4020"}, Float64}, Float64, 2}}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/sqhTO/src/jacobian.jl:21
 [13] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#4015#4016"{var"#4019#4020"}, Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#4015#4016"{var"#4019#4020"}, Float64}, Float64, 2}}}) (repeats 2 times)
    @ ForwardDiff ~/.julia/packages/ForwardDiff/sqhTO/src/jacobian.jl:19
 [14] ∂(f::var"#4019#4020", z::ComplexF64)

Other failing examples: t -> exp(-t^2), z -> cospi(z)^2.

Unsure if related to #486

[tpapp]

It looks like Base._cpow does not handle this well, the solution could be defining a rule for the derivative.

Using z*z instead is a workaround.

[jagot]

Yes, I found out that that works. Now I only need to figure out if I want z^2 == z*z or conj(z)*z == abs2(z), but that's an issue for me 😄

[andreasnoack]

Isn't this a dup of #486?

[antoine-levitt]

This appears to be "fixed" in 1.7: it doesn't give an error, but it gives wrong results (which is worse). The issue is relatively subtle: the derivative at real numbers in complex directions is wrong:

function f(x)
    imag((3+(1+im)*x)^2)
end
println(ForwardDiff.derivative(f, 0))
println((f(0+1e-8) - f(0))/1e-8)

yields 1 and 6 respectively. Where should this be fixed? Here? In DiffRules?

[antoine-levitt]

Pinging @devmotion because of the recent complex numbers PR (which doesn't fix this issue)

[devmotion]

Where should this be fixed? Here? In DiffRules?

I didn't follow this issue bit can have a look. Does ForwardDiff use DiffRules for these calculations here? If yes, then it should be possible to obtain these incorrect derivatives without ForwardDiff by just evaluating the corresponding DiffRules definition (since we don't do anything special and just use these rules here). The recent PR just fixed a bug with rules that return Complex derivatives such that they can be used for intermediate calculations, it should be a strict improvement and finally fixes downstream tests in DiffRules (before CI always failed when adding such rules since they had to be excluded from ForwardDiff tests manually first).

[antoine-levitt]

Does ForwardDiff use DiffRules for these calculations here?

I don't think so:

@which((3 + (1 + im) * x) ^ 2) = ^(z::Complex, n::Integer) in Base at complex.jl:834

By contrast the real version (without the + im) is defined in forwarddiff. So I guess this should be added to ForwardDiff directly.

[antoine-levitt]

To reiterate, this is a pretty nasty bug, I know several people who've been bitten by this now. The simplest way forward here is to merge JuliaDiff/DiffRules.jl#37 or something like it, then define rules for f(x::Complex{Dual}) based on this.

[mcabbott]

The behaviour of the examples listed above on ForwardDiff#master is:

julia> ∂(x -> x^2, 0.1)
0.2

julia> ∂(x -> x^2, 0.1+0.0im)
0.2 + 0.0im

julia> ∂(t -> exp(-t^2), 0.1)
-0.19800996674983362

julia> ∂(t -> exp(-t^2), 0.1+0.0im)
-0.19800996674983362 + 0.0im

julia> ∂(z -> cospi(z)^2, 0.1)
-1.8465818304904567

julia> ∂(z -> cospi(z)^2, 0.1 + 0.0im)
-1.8465818304904567 + 0.0im

and

julia> function f(x)
           imag((3+(1+im)*x)^2)
       end
f (generic function with 2 methods)

julia> println(ForwardDiff.derivative(f, 0))
6.0

julia> println((f(0+1e-8) - f(0))/1e-8)
6.00000002

[antoine-levitt]

Ah I did not see this. Yes, I think we can close this now (but it's annoying that the fix is not in a release!)