StackOverflowError with complex `sqrt`

[fmeirinhos]

When differentiating over a complex sqrt

x0 = 5 ForwardDiff.derivative(x -> sqrt(1im * x), x0)

a StackOverflowError is encountered.

StackOverflowError: Stacktrace: [1] sqrt(::Complex{ForwardDiff.Dual{ForwardDiff.Tag{##20#21,Int64},Float16,1}}) at ./complex.jl:416 (repeats 80000 times)

The stacktrace points to the following declaration sqrt(z::Complex) = sqrt(float(z))

I have found similar issues #111 #196 #239, but I don't see how I can fix this problem.

[tkoolen]

This is basically https://github.com/JuliaLang/julia/issues/26552 I believe.

[EdoAlvarezR]

I'm running into the same issue when trying to get a derivative over a complex sqrt operation. Here I was able to reproduce the error in a simplified way:

> x = ForwardDiff.Dual{Nothing}(1,0.0,0.0,0.0)*im
> sqrt(x)
StackOverflowError:

Stacktrace:
 [1] Type at ./complex.jl:12 [inlined] (repeats 2 times)
 [2] float at ./complex.jl:967 [inlined]
 [3] sqrt(::Complex{ForwardDiff.Dual{Nothing,Float64,3}}) at ./complex.jl:461 (repeats 80000 times)

Has anyone found a solution for this?

[chriscoey]

I too am running into this

[adam-coogan]

This is a problem for me as well. Any ideas about how to deal with this?

[jfmoulin]

I run into same error with a just updated Turing/ForwardDiff.

[andreasnoack]

See the discussion in https://github.com/JuliaDiff/ForwardDiff.jl/issues/157

[baggepinnen]

I am trying this out in julia v1.6 and it's sometimes working okay now. sqrt(::Complex) contains a call to ssqs which has a branch that sometimes fails on Dual: the problematic function is Base.exponent which is not defined for Dual

 function Base.ssqs(x::T, y::T) where T <: ForwardDiff.Dual
           k::Int = 0
           ρ = x*x + y*y
           @show x
           if !isfinite(ρ) && (isinf(x) || isinf(y))
               ρ = convert(T, Inf)
           elseif isinf(ρ) || (ρ==0 && (x!=0 || y!=0)) || ρ<nextfloat(zero(T))/(2*eps(T)^2)
               m::T = max(abs(x), abs(y))
               @show m
               k = m==0 ? m : Base.exponent(m)
               xk, yk = Base.ldexp(x,-k), Base.ldexp(y,-k)
               ρ = xk*xk + yk*yk
           end
           ρ, k
       end