Error with ForwardDiff 1 and QuadGK

[sschlenkrich]

Hi ForwardDiff team,

I observe an error when updating from ForwardDiff v0.10 to v1.x. The error occurs when differentiating through quadgk integration.

A MWE to reproduce the error is as follows:

using ForwardDiff
using QuadGK

delta = [ 1.0, 10.0 ]
chi = [ 0.03, 0.30 ]
S = inv([ exp(-chi_ * delta_) for delta_ in delta, chi_ in chi ])

g(x) = quadgk(t -> S*x, 0.0, 1.0)[1]

ForwardDiff.jacobian(g, zeros(2))

With ForwardDiff 1.3, I get the following error:

ERROR: LoadError: DomainError with 0.5:
integrand produced Dual{ForwardDiff.Tag{typeof(g), Float64}}(NaN,NaN,NaN) in the interval (0.0, 1.0)

Note that the code does not throw with ForwardDiff v0.10.39.

Also note that when changing the matrx S, e.g. adding a line S .+= 0.5 then above code also works with v1.3.

Would be great if you could have a look.

[devmotion]

Can you include the full stacktrace in the issue?

[sschlenkrich]

Sure, please see below.

My guess is, it has to do with calling norm(...) on a Dual at zero.

ERROR: LoadError: DomainError with 0.5:
integrand produced Dual{ForwardDiff.Tag{typeof(g), Float64}}(NaN,NaN,NaN) in the interval (0.0, 1.0)
Stacktrace:
  [1] evalrule(f::var"#g##18#g##19"{…}, a::Float64, b::Float64, x::Vector{…}, w::Vector{…}, wg::Vector{…}, nrm::typeof(norm))
    @ QuadGK C:\Users\Sebastian\.julia\packages\QuadGK\7rND3\src\evalrule.jl:39
  [2] #do_quadgk##4
    @ C:\Users\Sebastian\.julia\packages\QuadGK\7rND3\src\adapt.jl:54 [inlined]
  [3] ntuple
    @ .\ntuple.jl:50 [inlined]
  [4] do_quadgk(f::var"#g##18#g##19"{…}, s::Tuple{…}, n::Int64, atol::Nothing, rtol::Nothing, maxevals::Int64,
 nrm::typeof(norm), _segbuf::Nothing, eval_segbuf::Nothing)
    @ QuadGK C:\Users\Sebastian\.julia\packages\QuadGK\7rND3\src\adapt.jl:52
  [5] #28
    @ C:\Users\Sebastian\.julia\packages\QuadGK\7rND3\src\api.jl:83 [inlined]
  [6] handle_infinities
    @ C:\Users\Sebastian\.julia\packages\QuadGK\7rND3\src\adapt.jl:189 [inlined]
  [7] #quadgk#26
    @ C:\Users\Sebastian\.julia\packages\QuadGK\7rND3\src\api.jl:82 [inlined]
  [8] quadgk
    @ C:\Users\Sebastian\.julia\packages\QuadGK\7rND3\src\api.jl:80 [inlined]
  [9] g
    @ C:\Users\Sebastian\01_GitHub\ufo.1\ufo.jl:9 [inlined]
 [10] vector_mode_dual_eval!
    @ C:\Users\Sebastian\.julia\packages\ForwardDiff\kQBw9\src\apiutils.jl:24 [inlined]
 [11] vector_mode_jacobian(f::typeof(g), x::Vector{…}, cfg::ForwardDiff.JacobianConfig{…})
    @ ForwardDiff C:\Users\Sebastian\.julia\packages\ForwardDiff\kQBw9\src\jacobian.jl:129
 [12] jacobian
    @ C:\Users\Sebastian\.julia\packages\ForwardDiff\kQBw9\src\jacobian.jl:22 [inlined]
 [13] jacobian(f::typeof(g), x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{…}, Float64, 2, Vector{…}})
    @ ForwardDiff C:\Users\Sebastian\.julia\packages\ForwardDiff\kQBw9\src\jacobian.jl:19
 [14] jacobian(f::typeof(g), x::Vector{Float64})
    @ ForwardDiff C:\Users\Sebastian\.julia\packages\ForwardDiff\kQBw9\src\jacobian.jl:19
 [15] top-level scope
    @ C:\Users\Sebastian\01_GitHub\ufo.1\ufo.jl:11
 [16] include(mapexpr::Function, mod::Module, _path::String)
    @ Base .\Base.jl:307
 [17] top-level scope
    @ REPL[46]:1
in expression starting at C:\Users\Sebastian\01_GitHub\ufo.1\ufo.jl:11
Some type information was truncated. Use `show(err)` to see complete types.

[sschlenkrich]

I did some further analysis. The error can be traced back to the evaluation of square root at zero.

Consider the following three calcultions of the 2-norm of x which should yield the same result:

using ForwardDiff
using LinearAlgebra

x = [
    ForwardDiff.Dual(0.0, 1.0, 0.0),
    ForwardDiff.Dual(0.0, 1.0, 1.0),
]

display( norm(x) )
display( sqrt(x[1]^2 + x[2]^2) )
display( (x[1]^2 + x[2]^2)^0.5 )

Output is

Dual{Nothing}(NaN,NaN,NaN)
Dual{Nothing}(0.0,NaN,NaN)
Dual{Nothing}(0.0,0.0,0.0)

I would say, output 2 and 3 should be equal either way. This is probbly a different discussion.

But output 1 should not have NaN as function value. This is what causes the error in quadgk.

As far as I can see, output 1 is linked to the evaluation of iszero(x). With ForwardDiff v1.x this evaluates false. And this causes division by zero.

https://github.com/JuliaLang/LinearAlgebra.jl/blob/b599095ef3da7ba7e950ee4700a3ba0fea047949/src/generic.jl#L588

I understand that iszero(x) equal to false is a deliberate choice/consequence from the change to v1.x. But in my view norm(x) should not be NaN.

Any chance this could be fixed?

[mcabbott]

Xref #547 for some earlier discussion of norm.

[devmotion]

As I mentioned in #547, sqrt behaves as expected if one enables NaN-safe mode:

julia> sqrt(x[1]^2 + x[2]^2)
Dual{Nothing}(0.0,0.0,0.0)

julia> (x[1]^2 + x[2]^2)^0.5
Dual{Nothing}(0.0,0.0,0.0)

Honestly, I still don't understand why this is not enabled by default. IMO it's the only sane way to use ForwardDiff - regardless of the performance (which was improved significantly by #777), IMO the default setting should always prioritize correctness. And I think non-NaN-safe mode is just not correct in general: If you input a number without perturbations (ie a non-Dual or a Dual with zero sensitivities), then the output should be unperturbed as well. That is, zero sensitivities should be strong zeros.

The norm case seems different and I agree with your analysis - it seems due to the (intended) behaviour in ForwardDiff@1, we miss the early exit in https://github.com/JuliaLang/LinearAlgebra.jl/blob/b599095ef3da7ba7e950ee4700a3ba0fea047949/src/generic.jl#L577 for vectors of duals with zero primals. Generally, I think norm deserves its own ForwardDiff rule. Then we'd also be able to hit the BLAS version it seems.