Open dpsanders opened 7 years ago
Crossref: https://github.com/JuliaDiff/ForwardDiff.jl/issues/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.
Thanks. What's the issue with z^3
?
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.
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.
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... 😄
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))
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?
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:When I try to differentiate this, I get the following error: