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 https://github.com/JuliaDiff/ForwardDiff.jl/issues/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.