Incorrect Hessian by exp function

[rzyu45]

I was trying to calculate the Hessian of a simple function such as

using ForwardDiff
function busPQ(Ybus, Vθ)
    n = length(Vθ)
    V=Vθ[1:Int(n/2)]
    θ=Vθ[Int(n/2)+1:n]
    Vbus=V.*exp.(im .*θ)
    # Vbus=V.*(cos.(θ)+im .*sin.(θ))
    # print(Vbus-Vbus)
    I = Ybus.*Vbus
    return vcat(reim(Vbus.*conj(I))...)
end

Ybus=[1 - 3im]
Vm=[1]
Va=[0]
Vθ = vcat(Vm,Va)
function Hessian(f, x)
    n = length(x)
    out = ForwardDiff.jacobian(x -> ForwardDiff.jacobian(f, x), x)
    return reshape(out, n, n, n)
end
H(y) = Hessian(x-> busPQ(Ybus, x), y)

At point [1,0], the above codes produce the wrong Hessian. However, if I set very small perturbation to Vθ, then the results are correct.

julia> H([1,0])
2×2×2 Array{Float64, 3}:
[:, :, 1] =
 2.0  0.0
 6.0  0.0

[:, :, 2] =
 0.0  2.0
 0.0  6.0

julia> H([1,1e-10])
2×2×2 Array{Float64, 3}:
[:, :, 1] =
 2.0  0.0
 6.0  0.0

[:, :, 2] =
 0.0  0.0
 0.0  0.0

julia> H([1,-1e-10])
2×2×2 Array{Float64, 3}:
[:, :, 1] =
 2.0  0.0
 6.0  0.0

[:, :, 2] =
 0.0  0.0
 0.0  0.0

It is found to be more robust to use cos.(θ)+im .*sin.(θ), which functioned correctly at point [1,0], instead of exp.(im .*θ).
I think the problem has something in common with #653.

The difference between cos.(θ)+im .*sin.(θ) and exp.(im .*θ) is printed as follows

Complex{ForwardDiff.Dual{ForwardDiff.Tag{var"#21#22", ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20"{var"#21#22"}, Int64}, Int64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20"{var"#21#22"}, Int64}, Float64, 2}, 2}}[Dual{ForwardDiff.Tag{var"#21#22", ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20"{var"#21#22"}, Int64}, Int64, 2}}}(Dual{ForwardDiff.Tag{var"#19#20"{var"#21#22"}, Int64}}(0.0,0.0,0.0),Dual{ForwardDiff.Tag{var"#19#20"{var"#21#22"}, Int64}}(0.0,0.0,0.0),Dual{ForwardDiff.Tag{var"#19#20"{var"#21#22"}, Int64}}(0.0,0.0,1.0)) + Dual{ForwardDiff.Tag{var"#21#22", ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20"{var"#21#22"}, Int64}, Int64, 2}}}(Dual{ForwardDiff.Tag{var"#19#20"{var"#21#22"}, Int64}}(0.0,0.0,0.0),Dual{ForwardDiff.Tag{var"#19#20"{var"#21#22"}, Int64}}(0.0,0.0,0.0),Dual{ForwardDiff.Tag{var"#19#20"{var"#21#22"}, Int64}}(0.0,0.0,0.0))*im]2×2×2 Array{Float64, 3}:

It is clear that the above Dual is not all zero.

The dependencies are ForwardDiff v0.10.36 and Julia V1.9.3.

[mcabbott]

This is probably fixed by #481, as master gives the right answer:

julia> H([1,0]) ≈ H([1,1e-10])
true

(jl_JRlWH7) pkg> st ForwardDiff
Status `/private/var/folders/yq/4p2zwd614y59gszh7y9ypyhh0000gn/T/jl_JRlWH7/Project.toml`
  [f6369f11] ForwardDiff v0.11.0-DEV `https://github.com/JuliaDiff/ForwardDiff.jl.git#master`

(#481 was pulled from v0.10 as it seemed that some people relied on the old behaviour, and we haven't tagged v0.11 yet... because making a breaking release of such a widely used package creates a lot of work downstream, and people have ideas for other things they would quite like to roll into a breaking 1.0 release.)