QuadGK.jl

QuadGK.jl copied to clipboard

I have a function if I have a function f(x) which is the integral another function (using QuadGK) and x is the upper limit of integration. If I want to take an automatic derivative of that function using ForwardDiff, it results in a stack overflow. Here are instructions to reproduce the error

using QuadGK
using ForwardDiff
D(f::Function, x::Number) = ForwardDiff.derivative(f, x)

f(x) = quadgk(y -> y^2, 0, x)
D(x -> quadgk(y -> y^2, 0, x), 1)

ERROR: StackOverflowError:
Stacktrace:
 [1] #quadgk#15(::Array{Any,1}, ::Function, ::Function, ::ForwardDiff.Dual{ForwardDiff.Tag{##3#5,Int64},Float64,1}, ::ForwardDiff.Dual{ForwardDiff.Tag{##3#5,Int64},Float64,1}) at /Users/mason/.julia/v0.6/QuadGK/src/QuadGK.jl:241 (repeats 65479 times)
 [2] derivative(::##3#5, ::Int64) at /Users/mason/.julia/v0.6/ForwardDiff/src/derivative.jl:14
 [3] D(::Function, ::Int64) at ./REPL[4]:1

Thanks for the report. The quadgk endpoints only support real and complex numbers right now, and probably there is an infinite loop in the promotion rules for other number types.

It would be good to make this work, if only to make sure the methods are defined generically enough to work with any number type.

(But even if we get this working, I’m guessing that using the analytical derivative “manually” will always be much more efficient.)

@stevengj Yeah, I wasn't expecting this to work but I decided to give it a try anyways and ended up segfaulting due to the stack overflow so it was suggested I submit an issue.

Also shows up for dimensionful quantities (#17).

The case of dimensionful quantities is fixed by #36.

I'm not sure how to fix this in general — the problem is that I need a package-independent way to get the underlying real-numeric types out of things like Measurement and Dual numbers, and there doesn't seem to be an API for this right now.

I also encountered this error. I have a non-linear problem trying to find zeros. I have a working solution in NLSolve but needed additional constraints, so I turned to JuMP.

Is there a solution for Dual numbers that exists now?

If you want the derivative with respect to the endpoints of a function that involves integration, you can just take the derivative analytically and expose it to automatic differentiation via the ChainRules package. The https://github.com/SciML/Quadrature.jl package (which wraps around QuadGK) does this for you.

Leibniz integral rule strikes again! Yes I will port over the Quadrature, thanks for the feedback.

SciML actually doesn't work:

julia> f(z) = solve(IntegralProblem((x, _)->sin(x)^2, 0.0, z), QuadGKJL())[1]
f (generic function with 2 methods)

julia> f(3.14)
1.5707963254482846

julia> ForwardDiff.derivative(f, 3.14)
ERROR: StackOverflowError:
Stacktrace:
 [1] cachedrule(#unused#::Type{ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Float64}, Float64, 1}}, n::Int64) (repeats 79984 times)
   @ QuadGK ~/.julia/packages/QuadGK/ENhXl/src/gausskronrod.jl:253

The original StackOverflowError error is gone, I'm going to close this issue.

It now gives the error:

ERROR: MethodError: no method matching kronrod(::Type{ForwardDiff.Dual{ForwardDiff.Tag{var"#3#5", Int64}, Float64, 1}}, ::Int64)
Closest candidates are:
  kronrod(::Type{T}, ::Integer) where T<:AbstractFloat at ~/.julia/dev/QuadGK/src/gausskronrod.jl:150

which stems from this line. I think it could be fixed if one(Dual) returned 1 instead of a dual number, which seems better to me, but in any case it seems like a separate issue (https://github.com/JuliaDiff/ForwardDiff.jl/issues/624).