Error trying to differentiate through a `Complex`

[dpsanders]

I need to treat a complex function f: C → C as a function from R^2 to R^2. I am doing this as follows:

julia 0.6> function realify(f)

               function g(x)
                   z = Complex(x[1], x[2])
                   z2 = f(z)
                   SVector(reim(z2))
               end

               return g

           end
realify (generic function with 1 method)

julia 0.6> using StaticArrays

julia 0.6> f(z) = z^3 - 1
f (generic function with 1 method)

julia 0.6> const g = realify(f)
g (generic function with 1 method)

julia 0.6> g(SVector(1, 2))
2-element SVector{2,Int64}:
 -12
  -2

When I try to differentiate this, I get the following error:

julia 0.6> using ForwardDiff

julia 0.6> ForwardDiff.jacobian(g, SVector(1, 2))
ERROR: MethodError: Cannot `convert` an object of type ForwardDiff.Dual{ForwardDiff.Tag{#g#1{#f},0x570e4f4889f9f966},Int64,2} to an object of type Signed
This may have arisen from a call to the constructor Signed(...),
since type constructors fall back to convert methods.
Stacktrace:
 [1] ^(::Complex{ForwardDiff.Dual{ForwardDiff.Tag{#g#1{#f},0x570e4f4889f9f966},Int64,2}}, ::Complex{ForwardDiff.Dual{ForwardDiff.Tag{#g#1{#f},0x570e4f4889f9f966},Int64,2}}) at ./complex.jl:625
 [2] f(::Complex{ForwardDiff.Dual{ForwardDiff.Tag{#g#1{#f},0x570e4f4889f9f966},Int64,2}}) at ./REPL[3]:1
 [3] (::#g#1{#f})(::SVector{2,ForwardDiff.Dual{ForwardDiff.Tag{#g#1{#f},0x570e4f4889f9f966},Int64,2}}) at ./REPL[1]:5
 [4] jacobian(::#g#1{#f}, ::SVector{2,Int64}) at /Users/dpsanders/.julia/v0.6/ForwardDiff/src/jacobian.jl:76

[ChrisRackauckas]

Crossref: #157

f(z) = z*z*z -1
const g = realify(f)
g(SVector(1.0, 2.0))

using ForwardDiff
ForwardDiff.jacobian(g, SVector(1.0, 2.0))

that works.

[dpsanders]

Thanks. What's the issue with z^3?

[ChrisRackauckas]

z = ForwardDiff.Dual{ForwardDiff.Tag{nothing,0xb81c941840f92fc0}}(1.0,1.0,0.0) +
    ForwardDiff.Dual{ForwardDiff.Tag{nothing,0xb81c941840f92fc0}}(2.0,0.0,1.0)*im
z^3

is an MWE, but it was surprisingly uninformative. More informative is:

z = ForwardDiff.Dual{ForwardDiff.Tag{nothing,0xb81c941840f92fc0}}(1.0,1.0,0.0) +
    ForwardDiff.Dual{ForwardDiff.Tag{nothing,0xb81c941840f92fc0}}(2.0,0.0,1.0)*im
a = 3
z^Complex(a)

is doing

zz,aa = (promote(z,a)...)
convert(Integer, aa)

which errors. So the "true MWE" is:

aa = ForwardDiff.Dual{ForwardDiff.Tag{nothing,0xb81c941840f92fc0}}(3.0,0.0,0.0) + ForwardDiff.Dual{ForwardDiff.Tag{nothing,0xb81c941840f92fc0}}(0.0,0.0,0.0)*im
convert(Integer, aa)

It's the promotion then the conversion to an integer which fails.

[ChrisRackauckas]

z = ForwardDiff.Dual{ForwardDiff.Tag{nothing,0xb81c941840f92fc0}}(1.0,1.0,0.0)
z^a
@which z^a

brings me here:

https://github.com/JuliaDiff/ForwardDiff.jl/blob/master/src/dual.jl#L100

So it looks like in the real case they had to work around this with:

https://github.com/JuliaDiff/ForwardDiff.jl/blob/master/src/dual.jl#L367

That would need to be extended to complex.

[jrevels]

That would need to be extended to complex.

The Dual type should only implement operations w.r.t. the Real interface.

This is a problem with Base's Complex exponentiation implementation, which @dpsanders and I have actually run into before (I thought we filed a Base issue/PR, but I can't find it). Here's the entry point to the problematic code.

The ^(::Complex, ::Integer) method here should just be this isinteger block by default; as it's currently defined, the code forces the exponent through a round-trip of unnecessary conversions/promotions (Integer --> Complex{Integer} --> Complex{<:Dual} --> Integer). There's no benefit to this, AFAICT; it's just a non-Julian organization of the code.

While in practice this problem can just be solved by fixing Complex's exponentiation implementation, it's is actually a symptom of a deeper assumption scattered throughout numeric Julia code: that numeric equality implies lossless convertibility to a certain set of types. For example, the problematic assumption here is that isinteger(x) implies that convert(Integer, x) is possible/lossless, where the actual documented definition of Base.isinteger(x) is "test whether x is numerically equal to some integer", not "test where x can be converted to Integer".

While such an assumption may seem naively sensible, it is untenable in practice for a language where subtype relations are purely nominal set relations (at least until we can define types like Dual{T<:Real} <: supertype(T)).

Of course, this issue is a poster child for why you shouldn't use argument-based overloading to inject/propagate metadata through native Julia code. Time to get back to work on Cassette.jl... 😄

[antoine-levitt]

Not sure if it's related or not, but I tried a complexification of the example in the manual:

using ForwardDiff

function g(x::Vector)
    sum(sin, x) + prod(tan, x) * sum(sqrt, x);
end

function f(x)
    y = g(x+im*x)
    return real(y) + imag(y)
end

srand(0)
x = rand(5) # small size for example's sake

gr = x -> ForwardDiff.gradient(f, x); # g = ∇f

display(gr(x))

That gives me a StackOverflow

julia> include("scratch.jl")
ERROR: LoadError: StackOverflowError:
Stacktrace:
 [1] sqrt(::Complex{ForwardDiff.Dual{ForwardDiff.Tag{#f,0xb046287d533b082d},Float64,5}}) at ./complex.jl:416 (repeats 74769 times)
 [2] _mapreduce(::Base.#sqrt, ::Base.#+, ::IndexLinear, ::Array{Complex{ForwardDiff.Dual{ForwardDiff.Tag{#f,0xb046287d533b082d},Float64,5}},1}) at ./reduce.jl:273
 [3] f at /home/antoine/scratch.jl:10 [inlined]
 [4] vector_mode_dual_eval(::#f, ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{#f,0xb046287d533b082d},Float64,5,Array{ForwardDiff.Dual{ForwardDiff.Tag{#f,0xb046287d533b082d},Float64,5},1}}) at /home/antoine/.julia/v0.6/ForwardDiff/src/utils.jl:36
 [5] vector_mode_gradient(::#f, ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{#f,0xb046287d533b082d},Float64,5,Array{ForwardDiff.Dual{ForwardDiff.Tag{#f,0xb046287d533b082d},Float64,5},1}}) at /home/antoine/.julia/v0.6/ForwardDiff/src/gradient.jl:94
 [6] gradient at /home/antoine/.julia/v0.6/ForwardDiff/src/gradient.jl:19 [inlined]
 [7] gradient(::#f, ::Array{Float64,1}) at /home/antoine/.julia/v0.6/ForwardDiff/src/gradient.jl:18
 [8] (::##3#4)(::Array{Float64,1}) at /home/antoine/scratch.jl:17
 [9] include_from_node1(::String) at ./loading.jl:569
 [10] include(::String) at ./sysimg.jl:14
while loading /home/antoine/scratch.jl, in expression starting on line 19

This is because

sqrt(z::Complex) = sqrt(float(z))

[jrevels]

ref JuliaLang/julia#24570

[mmikhasenko]

Do I understand correctly that this is the same issue which prevents me using JuMP minimizer involving complex functions?

In the code below, I look for zeros of the complex function minimizing abs value.
It gives the StackOverflowError because of sqrt(s). It will neither work with log, nor with ^.

exmp(s) = sqrt(s)-1-1im
exmp2(x1,x2) = begin
    v2 = abs(exmp(x1+1im*x2))
end
m = Model(solver=NLoptSolver(algorithm=:LD_MMA))
@variable(m, sp_x, start = 1.0)
@variable(m, sp_y, start = 1.0)
JuMP.register(m, :exmp2, 2, exmp2, autodiff=true)
@NLobjective(m, Min, exmp2(sp_x,sp_y))
status = solve(m)

Error:

StackOverflowError:
Stacktrace:
 [1] sqrt(::Complex{ForwardDiff.Dual{ForwardDiff.Tag{JuMP.##166#168{#exmp2},Float64},Float64,2}}) at ./complex.jl:416 (repeats 80000 times)

Is there any workaround?