Issue with ForwardDiff.jl

[weinbe58]

We're having difficulties using the latest patch release of ForwardDiff.jl in the following situation:

ForwardDiff.derivative(2.0) do x
    h = fill(zero(typeof(x)), (2, 2))
    F = ExponentialUtilities.exponential!(im * h, ExpMethodGeneric())
    return sum(F)
end

is returning a NaN.

[mcabbott]

On what versions does this return a NaN?

I see some errors instead:

julia> using ForwardDiff, ExponentialUtilities, FiniteDifferences

julia> FiniteDifferences.grad(central_fdm(7, 1), x -> sum(sin, exp(x)), [1 2; 3 4.0])[1]
2×2 Matrix{Float64}:
 55.4447  123.025
 82.7871  182.245

julia> FiniteDifferences.grad(central_fdm(7, 1), x -> sum(sin, ExponentialUtilities.exponential!(x, ExpMethodGeneric())), [1 2; 3 4.0])[1]
2×2 Matrix{Float64}:
 55.4447  123.025
 82.7871  182.245

julia> ForwardDiff.gradient(x -> sum(sin, ExponentialUtilities.exponential!(x, ExpMethodGeneric())), [1 2; 3 4.0])
2×2 Matrix{Float64}:
 55.4447  123.025
 82.7871  182.245

Now with x = zeros(2,2) to get errors:

julia> FiniteDifferences.grad(central_fdm(7, 1), x -> sum(sin, exp(x)), zeros(2,2))[1]
2×2 Matrix{Float64}:
 0.540302  1.0
 1.0       0.540302

julia> ForwardDiff.gradient(x -> sum(sin, ExponentialUtilities.exponential!(x, ExpMethodGeneric())), zeros(2,2))
ERROR: ArgumentError: matrix contains Infs or NaNs
Stacktrace:
  [1] chkfinite
    @ ~/.julia/dev/julia/usr/share/julia/stdlib/v1.10/LinearAlgebra/src/lapack.jl:86 [inlined]
  [2] generic_lufact!(A::Matrix{Dual{ForwardDiff.Tag{var"#63#64", Float64}, Float64, 4}}, pivot::RowMaximum; check::Bool)
    @ LinearAlgebra ~/.julia/dev/julia/usr/share/julia/stdlib/v1.10/LinearAlgebra/src/lu.jl:135
  [3] generic_lufact!
    @ ~/.julia/dev/julia/usr/share/julia/stdlib/v1.10/LinearAlgebra/src/lu.jl:133 [inlined]
  [4] #lu!#174
    @ ~/.julia/dev/julia/usr/share/julia/stdlib/v1.10/LinearAlgebra/src/lu.jl:131 [inlined]
  [5] lu! (repeats 2 times)
    @ ~/.julia/dev/julia/usr/share/julia/stdlib/v1.10/LinearAlgebra/src/lu.jl:131 [inlined]
  [6] _rdiv!
    @ ~/.julia/packages/ExponentialUtilities/Jn8Ch/src/exp_generic.jl:188 [inlined]
  [7] exp_generic_core!(y1::Matrix{Dual{ForwardDiff.Tag{var"#63#64", Float64}, Float64, 4}}, y2::Matrix{Dual{ForwardDiff.Tag{var"#63#64", Float64}, Float64, 4}}, y3::Matrix{Dual{ForwardDiff.Tag{var"#63#64", Float64}, Float64, 4}}, x::Matrix{Dual{ForwardDiff.Tag{var"#63#64", Float64}, Float64, 4}}, #unused#::Val{13})
    @ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/Jn8Ch/src/exp_generic.jl:197
  [8] exp_generic!(y1::Matrix{Dual{ForwardDiff.Tag{var"#63#64", Float64}, Float64, 4}}, y2::Matrix{Dual{ForwardDiff.Tag{var"#63#64", Float64}, Float64, 4}}, y3::Matrix{Dual{ForwardDiff.Tag{var"#63#64", Float64}, Float64, 4}}, x::Matrix{Dual{ForwardDiff.Tag{var"#63#64", Float64}, Float64, 4}}, s::Int64, #unused#::Val{13})
    @ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/Jn8Ch/src/exp_generic.jl:172
  [9] exp_generic_mutable(x::Matrix{Dual{ForwardDiff.Tag{var"#63#64", Float64}, Float64, 4}}, s::Int64, #unused#::Val{13})
    @ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/Jn8Ch/src/exp_generic.jl:168
 [10] exponential!(x::Matrix{Dual{ForwardDiff.Tag{var"#63#64", Float64}, Float64, 4}}, method::ExpMethodGeneric{Val{13}()}, cache::Nothing)
    @ ExponentialUtilities ~/.julia/packages/ExponentialUtilities/Jn8Ch/src/exp_generic.jl:138
 [11] exponential!
    @ ~/.julia/packages/ExponentialUtilities/Jn8Ch/src/exp_generic.jl:130 [inlined]
 [12] #63
    @ ./REPL[54]:1 [inlined]
 [13] vector_mode_dual_eval!
    @ ~/.julia/packages/ForwardDiff/eqMFf/src/apiutils.jl:37 [inlined]
 [14] vector_mode_gradient(f::var"#63#64", x::Matrix{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{var"#63#64", Float64}, Float64, 4, Matrix{Dual{ForwardDiff.Tag{var"#63#64", Float64}, Float64, 4}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/eqMFf/src/gradient.jl:106
 [15] gradient(f::Function, x::Matrix{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{var"#63#64", Float64}, Float64, 4, Matrix{Dual{ForwardDiff.Tag{var"#63#64", Float64}, Float64, 4}}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/eqMFf/src/gradient.jl:0
 [16] gradient(f::Function, x::Matrix{Float64}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{var"#63#64", Float64}, Float64, 4, Matrix{Dual{ForwardDiff.Tag{var"#63#64", Float64}, Float64, 4}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/eqMFf/src/gradient.jl:17
 [17] gradient(f::Function, x::Matrix{Float64})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/eqMFf/src/gradient.jl:17

# same error from example above, does not return NaN
julia> ForwardDiff.derivative(2.0) do x
           h = fill(zero(typeof(x)), (2, 2))
           F = ExponentialUtilities.exponential!(im * h, ExpMethodGeneric())
           return sum(F)  # note that this is complex, and indep of x
       end
ERROR: ArgumentError: matrix contains Infs or NaNs
Stacktrace:
  [1] chkfinite(A::Matrix{Complex{ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20", Float64}, Float64, 1}}})
    @ LinearAlgebra.LAPACK ~/.julia/dev/julia/usr/share/julia/stdlib/v1.10/LinearAlgebra/src/lapack.jl:86
  [2] generic_lufact!(A::Matrix{Complex{ForwardDiff.Dual{ForwardDiff.Tag{var"#19#20", Float64}, Float64, 1}}}, pivot::RowMaximum; check::Bool)

julia> lu([1 2; 3 4.0]) |> typeof
LU{Float64, Matrix{Float64}, Vector{Int64}}

# this does not go to /LinearAlgebra/src/lu.jl:131
julia> lu(zeros(2,2))
ERROR: SingularException(1)
Stacktrace:
 [1] checknonsingular
   @ ~/.julia/dev/julia/usr/share/julia/stdlib/v1.10/LinearAlgebra/src/factorization.jl:19 [inlined]
 [2] checknonsingular
   @ ~/.julia/dev/julia/usr/share/julia/stdlib/v1.10/LinearAlgebra/src/factorization.jl:22 [inlined]
 [3] #lu!#170
   @ ~/.julia/dev/julia/usr/share/julia/stdlib/v1.10/LinearAlgebra/src/lu.jl:82 [inlined]
 [4] lu!
   @ ~/.julia/dev/julia/usr/share/julia/stdlib/v1.10/LinearAlgebra/src/lu.jl:80 [inlined]
 [5] #lu#176
   @ ~/.julia/dev/julia/usr/share/julia/stdlib/v1.10/LinearAlgebra/src/lu.jl:299 [inlined]
 [6] lu
   @ ~/.julia/dev/julia/usr/share/julia/stdlib/v1.10/LinearAlgebra/src/lu.jl:298 [inlined]
 [7] lu(A::Matrix{Float64})
   @ LinearAlgebra ~/.julia/dev/julia/usr/share/julia/stdlib/v1.10/LinearAlgebra/src/lu.jl:298

(jl_kR2r1x) pkg> st
Status `/private/var/folders/yq/4p2zwd614y59gszh7y9ypyhh0000gn/T/jl_kR2r1x/Project.toml`
  [d4d017d3] ExponentialUtilities v1.22.0
  [26cc04aa] FiniteDifferences v0.12.25
  [f6369f11] ForwardDiff v0.10.33

With older ForwardDiff, two error cases above run:

julia> ForwardDiff.gradient(x -> sum(sin, ExponentialUtilities.exponential!(x, ExpMethodGeneric())), zeros(2,2))
2×2 Matrix{Float64}:
 0.540302  1.0
 1.0       0.540302

julia> ForwardDiff.derivative(2.0) do x  # nothing depends on x
           h = fill(zero(typeof(x)), (2, 2))
           F = ExponentialUtilities.exponential!(im * h, ExpMethodGeneric())
           return @show sum(F)  # result is complex
       end
sum(F) = Dual{ForwardDiff.Tag{var"#15#16", Float64}}(2.0,0.0) + Dual{ForwardDiff.Tag{var"#15#16", Float64}}(-0.0,0.0)*im
0.0 + 0.0im

(jl_x4ENQb) pkg> st
Status `/private/var/folders/yq/4p2zwd614y59gszh7y9ypyhh0000gn/T/jl_x4ENQb/Project.toml`
  [d4d017d3] ExponentialUtilities v1.22.0
⌃ [f6369f11] ForwardDiff v0.10.32