Error with ForwardDiff 1 and QuadGK

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.

Can you include the full stacktrace in the issue?

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.

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?

Xref #547 for some earlier discussion of norm.

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.

Navigation Menu

Repository navigation

Error with ForwardDiff 1 and QuadGK #785

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Error with ForwardDiff 1 and QuadGK #785

Description

Activity

devmotion commented on Nov 17, 2025

sschlenkrich commented on Nov 17, 2025

sschlenkrich commented on Nov 17, 2025

mcabbott commented on Nov 17, 2025

devmotion commented on Nov 17, 2025

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions