Incorrect Hessian by `exp`

function · Issue #677 · JuliaDiff/ForwardDiff.jl (github.com)

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`

.