[Help] Nested ForwardDiff.jacobian of an Integrals.IntegralProblem

Specific Domains Numerics
, ,
1

I have an integration problem (actually vector-valued, but problem here applies to scalar-valued ones as well) as a part of my model which it might be helpful to calculate its Hessian. I follow the documentation of Integrals.jl, however the nested Jacobian doesn’t seem to work for some solver algorithms. I suspect it’s due to the problem with propagating Dual types, but not sure, and hence also don’t know a priori which solvers are safe, and if there are work-around for the problematic solvers for this case. I would be happy if someone can help me to understand the issue here. Thank you!

Here is MWE, where HCubatureJL works fine but CubatureJLh doesn’t:

using Integrals, ForwardDiff
using IntegralsCubature
my_parameters = [1.0, 2.0]
my_function(x, p) = x^2 + p[1] * x + p[2]
function my_integration(p)
    my_problem = IntegralProblem(my_function, -1.0, 1.0, p)
    # return solve(my_problem, HCubatureJL(), reltol=1e-3, abstol=1e-3)  # Works
    return solve(my_problem, CubatureJLh(), reltol=1e-3, abstol=1e-3)  # Errors
end
my_solution = my_integration(my_parameters)
my_jakobian = ForwardDiff.jacobian(my_integration, my_parameters)
ForwardDiff.jacobian(x->ForwardDiff.jacobian(my_integration, x), my_parameters)

and the error:

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6})
Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  (::Type{T})(::T) where T<:Number at boot.jl:772
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
  ...
Stacktrace:
  [1] (::Cubature.var"#17#18"{Bool, Bool, Int64, Float64, Float64, Int64, Int32, Ptr{Nothing}, Cubature.IntegrandData{IntegralsCubature.var"#3#15"{IntegralProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}}, typeof(my_function), Float64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}}}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Int64})()
    @ Cubature ~/.julia/packages/Cubature/5zwuu/src/Cubature.jl:215
  [2] disable_sigint
    @ ./c.jl:473 [inlined]
  [3] cubature(xscalar::Bool, fscalar::Bool, vectorized::Bool, padaptive::Bool, fdim::Int64, f::IntegralsCubature.var"#3#15"{IntegralProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}}, typeof(my_function), Float64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}}, xmin_::Float64, xmax_::Float64, reqRelError::Float64, reqAbsError::Float64, maxEval::Int64, error_norm::Int32)
    @ Cubature ~/.julia/packages/Cubature/5zwuu/src/Cubature.jl:169
  [4] hquadrature(f::Function, xmin::Float64, xmax::Float64; reltol::Float64, abstol::Float64, maxevals::Int64)
    @ Cubature ~/.julia/packages/Cubature/5zwuu/src/Cubature.jl:230
  [5] __solvebp_call(::IntegralProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}}, typeof(my_function), Float64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::CubatureJLh, ::Integrals.ReCallVJP{Integrals.ZygoteVJP}, ::Float64, ::Float64, ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}; reltol::Float64, abstol::Float64, maxiters::Int64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ IntegralsCubature ~/.julia/packages/IntegralsCubature/adYuM/src/IntegralsCubature.jl:29
  [6] __solvebp(prob::IntegralProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}}, typeof(my_function), Float64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, alg::CubatureJLh, sensealg::Integrals.ReCallVJP{Integrals.ZygoteVJP}, lb::Float64, ub::Float64, p::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}})
    @ Integrals.IntegralsForwardDiffExt ~/.julia/packages/Integrals/9qNWp/ext/IntegralsForwardDiffExt.jl:23
  [7] solve(prob::IntegralProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}}, typeof(my_function), Float64, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, alg::CubatureJLh; sensealg::Integrals.ReCallVJP{Integrals.ZygoteVJP}, do_inf_transformation::Nothing, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}})
    @ Integrals ~/.julia/packages/Integrals/9qNWp/src/Integrals.jl:108
  [8] my_integration(p::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}})
    @ Main ~/Documents/Research_local/Projects_local/SpatialCompetition/test/main.jl:27
  [9] vector_mode_dual_eval!
    @ ~/.julia/packages/ForwardDiff/vXysl/src/apiutils.jl:24 [inlined]
 [10] vector_mode_jacobian(f::typeof(my_integration), x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:125
 [11] jacobian(f::Function, x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:21
 [12] jacobian(f::Function, x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}, 6}}}) (repeats 2 times)
    @ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:19
 [13] #11
    @ ~/Documents/Research_local/Projects_local/SpatialCompetition/test/main.jl:31 [inlined]
 [14] vector_mode_dual_eval!
    @ ~/.julia/packages/ForwardDiff/vXysl/src/apiutils.jl:24 [inlined]
 [15] vector_mode_jacobian(f::var"#11#12", x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:125
 [16] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:21
 [17] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#11#12", Float64}, Float64, 6}}}) (repeats 2 times)
    @ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:19
 [18] top-level scope
    @ ~/Documents/Research_local/Projects_local/SpatialCompetition/test/main.jl:31
2

I didn’t run the code, but does it make a difference if the bounds are given as vectors, e.g. IntegralProblem(my_function, -ones(1), ones(1), p)? I had similar looking errors once in the past, but not sure…

3

Thank you @SteffenPL, but unfortunately not. In order to have vector bounds, I vectorized the integrand my_function = @. x^2 + p[1] * x + p[2] (since otherwise 1-element vector bounds are not unwrapped automatically), but I got the same error:

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2})
Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  (::Type{T})(::T) where T<:Number at boot.jl:772
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50
  ...
Stacktrace:
  [1] (::Cubature.var"#17#18"{Bool, Bool, Int64, Float64, Float64, Int64, Int32, Ptr{Nothing}, Cubature.IntegrandData{IntegralsCubature.var"#3#15"{IntegralProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2}}, typeof(my_function), Vector{Float64}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}}}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Int64})()
    @ Cubature ~/.julia/packages/Cubature/5zwuu/src/Cubature.jl:215
  [2] disable_sigint
    @ ./c.jl:473 [inlined]
  [3] cubature(xscalar::Bool, fscalar::Bool, vectorized::Bool, padaptive::Bool, fdim::Int64, f::IntegralsCubature.var"#3#15"{IntegralProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2}}, typeof(my_function), Vector{Float64}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}}, xmin_::Vector{Float64}, xmax_::Vector{Float64}, reqRelError::Float64, reqAbsError::Float64, maxEval::Int64, error_norm::Int32)
    @ Cubature ~/.julia/packages/Cubature/5zwuu/src/Cubature.jl:169
  [4] hcubature(f::Function, xmin::Vector{Float64}, xmax::Vector{Float64}; reltol::Float64, abstol::Float64, maxevals::Int64)
    @ Cubature ~/.julia/packages/Cubature/5zwuu/src/Cubature.jl:230
  [5] __solvebp_call(::IntegralProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2}}, typeof(my_function), Vector{Float64}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::CubatureJLh, ::Integrals.ReCallVJP{Integrals.ZygoteVJP}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}; reltol::Float64, abstol::Float64, maxiters::Int64, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ IntegralsCubature ~/.julia/packages/IntegralsCubature/adYuM/src/IntegralsCubature.jl:40
  [6] __solvebp(prob::IntegralProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2}}, typeof(my_function), Vector{Float64}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, alg::CubatureJLh, sensealg::Integrals.ReCallVJP{Integrals.ZygoteVJP}, lb::Vector{Float64}, ub::Vector{Float64}, p::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2}}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}})
    @ Integrals.IntegralsForwardDiffExt ~/.julia/packages/Integrals/9qNWp/ext/IntegralsForwardDiffExt.jl:23
  [7] #solve#34
    @ ~/.julia/packages/Integrals/9qNWp/src/Integrals.jl:108 [inlined]
  [8] my_integration(p::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2}})
    @ Main ~/Documents/Research_local/Projects_local/SpatialCompetition/test/main.jl:27
  [9] vector_mode_dual_eval!
    @ ~/.julia/packages/ForwardDiff/vXysl/src/apiutils.jl:24 [inlined]
 [10] vector_mode_jacobian(f::typeof(my_integration), x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:125
 [11] jacobian(f::Function, x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2}}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:21
 [12] jacobian(f::Function, x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(my_integration), ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}, ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}, 2}}}) (repeats 2 times)
    @ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:19
 [13] #5
    @ ~/Documents/Research_local/Projects_local/SpatialCompetition/test/main.jl:31 [inlined]
 [14] vector_mode_dual_eval!
    @ ~/.julia/packages/ForwardDiff/vXysl/src/apiutils.jl:24 [inlined]
 [15] vector_mode_jacobian(f::var"#5#6", x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:125
 [16] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}}, ::Val{true})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:21
 [17] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{var"#5#6", Float64}, Float64, 2}}}) (repeats 2 times)
    @ ForwardDiff ~/.julia/packages/ForwardDiff/vXysl/src/jacobian.jl:19
 [18] top-level scope
    @ ~/Documents/Research_local/Projects_local/SpatialCompetition/test/main.jl:31
1 Like