Hessian throws error

[fernando-duarte]

I am trying to run the example from the Readme.md file but computing hessians instead of gradients, and getting an error for both Zygote and ForwardDiff.

using Quadrature, ForwardDiff, FiniteDiff, Zygote, Cuba
f(x,p) = sum(sin.(x .* p))
lb = ones(2)
ub = 3ones(2)
p = [1.5,2.0]

function testf(p)
    prob = QuadratureProblem(f,lb,ub,p)
    sin(solve(prob,CubaCuhre(),reltol=1e-6,abstol=1e-6)[1])
end
    
hp1 = Zygote.hessian(testf,p) # throws error
hp2 = FiniteDiff.finite_difference_hessian(testf,p) # ok
hp3 = ForwardDiff.hessian(testf,p) # throws error

Status `~/.julia/environments/v1.6/Project.toml`
[67601950] Quadrature v1.9.0
[f6369f11] ForwardDiff v0.10.18
[e88e6eb3] Zygote v0.6.14
[8a292aeb] Cuba v2.2.0
[6a86dc24] FiniteDiff v2.8.0

[fernando-duarte]

stacktrace for ForwardDiff.hessian(testf,p)

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), 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:760 (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50 ...

Stacktrace: [1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}) @ Base ./number.jl:7 [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}, i1::Int64) @ Base ./array.jl:839 [3] macro expansion @ ./broadcast.jl:984 [inlined] [4] macro expansion @ ./simdloop.jl:77 [inlined] [5] copyto! @ ./broadcast.jl:983 [inlined] [6] copyto! @ ./broadcast.jl:947 [inlined] [7] materialize! @ ./broadcast.jl:894 [inlined] [8] materialize!(dest::Vector{Float64}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}, Nothing, typeof(*), Tuple{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}, Float64}}) @ Base.Broadcast ./broadcast.jl:891 [9] (::Quadrature.var"#79#95"{QuadratureProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}, 2}}, typeof(f), Vector{Float64}, Vector{Float64}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Vector{Float64}, Vector{Float64}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}}})(x::Vector{Float64}, dx::Vector{Float64}) @ Quadrature ~/.julia/packages/Quadrature/rcCut/src/Quadrature.jl:412 [10] generic_integrand!(ndim::Int32, x_::Ptr{Float64}, ncomp::Int32, f_::Ptr{Float64}, func!::Quadrature.var"#79#95"{QuadratureProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}, 2}}, typeof(f), Vector{Float64}, Vector{Float64}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Vector{Float64}, Vector{Float64}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}}}) @ Cuba ~/.julia/packages/Cuba/KIQTB/src/Cuba.jl:92 [11] dointegrate! @ ~/.julia/packages/Cuba/KIQTB/src/cuhre.jl:40 [inlined] [12] dointegrate @ ~/.julia/packages/Cuba/KIQTB/src/Cuba.jl:195 [inlined] [13] cuhre(integrand::Quadrature.var"#79#95"{QuadratureProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}, 2}}, typeof(f), Vector{Float64}, Vector{Float64}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Vector{Float64}, Vector{Float64}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}}}, ndim::Int64, ncomp::Int64; nvec::Int64, rtol::Float64, atol::Float64, flags::Int64, minevals::Int64, maxevals::Int64, key::Int64, statefile::String, spin::Ptr{Nothing}, abstol::Missing, reltol::Missing) @ Cuba ~/.julia/packages/Cuba/KIQTB/src/cuhre.jl:97 [14] __solvebp_call(::QuadratureProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}, 2}}, typeof(f), Vector{Float64}, Vector{Float64}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::CubaCuhre, ::Quadrature.ReCallVJP{Quadrature.ZygoteVJP}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}}; reltol::Float64, abstol::Float64, maxiters::Int64, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) @ Quadrature ~/.julia/packages/Quadrature/rcCut/src/Quadrature.jl:471 [15] __solvebp(::QuadratureProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}, 2}}, typeof(f), Vector{Float64}, Vector{Float64}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::CubaCuhre, ::Quadrature.ReCallVJP{Quadrature.ZygoteVJP}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}, 2}}; kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}}) @ Quadrature ~/.julia/packages/Quadrature/rcCut/src/Quadrature.jl:583 [16] solve(::QuadratureProblem{false, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}, 2}}, typeof(f), Vector{Float64}, Vector{Float64}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::CubaCuhre; sensealg::Quadrature.ReCallVJP{Quadrature.ZygoteVJP}, kwargs::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}}) @ Quadrature ~/.julia/packages/Quadrature/rcCut/src/Quadrature.jl:153 [17] testf(p::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}, 2}}) @ Main ./In[1]:9 [18] vector_mode_dual_eval @ ~/.julia/packages/ForwardDiff/QOqCN/src/apiutils.jl:37 [inlined] [19] vector_mode_gradient(f::typeof(testf), x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(testf), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}, 2}}}) @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:106 [20] gradient @ ~/.julia/packages/ForwardDiff/QOqCN/src/gradient.jl:19 [inlined] [21] #106 @ ~/.julia/packages/ForwardDiff/QOqCN/src/hessian.jl:16 [inlined] [22] vector_mode_dual_eval @ ~/.julia/packages/ForwardDiff/QOqCN/src/apiutils.jl:37 [inlined] [23] vector_mode_jacobian(f::ForwardDiff.var"#106#107"{typeof(testf), ForwardDiff.HessianConfig{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}, 2}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}}}}, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}}}) @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/jacobian.jl:147 [24] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}}}, ::Val{false}) @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/jacobian.jl:21 [25] hessian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.HessianConfig{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}, 2}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}}}, ::Val{true}) @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/hessian.jl:17 [26] hessian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.HessianConfig{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}, 2}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(testf), Float64}, Float64, 2}}}) (repeats 2 times) @ ForwardDiff ~/.julia/packages/ForwardDiff/QOqCN/src/hessian.jl:15 [27] top-level scope @ In[5]:1 [28] eval @ ./boot.jl:360 [inlined] [29] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1094

stacktrace for Zygote.hessian(testf,p)

MethodError: no method matching Float64(::ForwardDiff.Dual{Nothing, 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:760 (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50 ...

Stacktrace: [1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{Nothing, Float64, 2}) @ Base ./number.jl:7 [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{Nothing, Float64, 2}, i1::Int64) @ Base ./array.jl:839 [3] macro expansion @ ./broadcast.jl:984 [inlined] [4] macro expansion @ ./simdloop.jl:77 [inlined] [5] copyto! @ ./broadcast.jl:983 [inlined] [6] copyto! @ ./broadcast.jl:947 [inlined] [7] materialize! @ ./broadcast.jl:894 [inlined] [8] materialize!(dest::Vector{Float64}, bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}, Nothing, typeof(*), Tuple{ForwardDiff.Dual{Nothing, Float64, 2}, Float64}}) @ Base.Broadcast ./broadcast.jl:891 [9] (::Quadrature.var"#79#95"{QuadratureProblem{false, Vector{ForwardDiff.Dual{Nothing, Float64, 2}}, typeof(f), Vector{Float64}, Vector{Float64}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Vector{Float64}, Vector{Float64}, Vector{ForwardDiff.Dual{Nothing, Float64, 2}}})(x::Vector{Float64}, dx::Vector{Float64}) @ Quadrature ~/.julia/packages/Quadrature/rcCut/src/Quadrature.jl:412 [10] generic_integrand!(ndim::Int32, x_::Ptr{Float64}, ncomp::Int32, f_::Ptr{Float64}, func!::Quadrature.var"#79#95"{QuadratureProblem{false, Vector{ForwardDiff.Dual{Nothing, Float64, 2}}, typeof(f), Vector{Float64}, Vector{Float64}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Vector{Float64}, Vector{Float64}, Vector{ForwardDiff.Dual{Nothing, Float64, 2}}}) @ Cuba ~/.julia/packages/Cuba/KIQTB/src/Cuba.jl:92 [11] dointegrate! @ ~/.julia/packages/Cuba/KIQTB/src/cuhre.jl:40 [inlined] [12] dointegrate @ ~/.julia/packages/Cuba/KIQTB/src/Cuba.jl:195 [inlined] [13] cuhre(integrand::Quadrature.var"#79#95"{QuadratureProblem{false, Vector{ForwardDiff.Dual{Nothing, Float64, 2}}, typeof(f), Vector{Float64}, Vector{Float64}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, Vector{Float64}, Vector{Float64}, Vector{ForwardDiff.Dual{Nothing, Float64, 2}}}, ndim::Int64, ncomp::Int64; nvec::Int64, rtol::Float64, atol::Float64, flags::Int64, minevals::Int64, maxevals::Int64, key::Int64, statefile::String, spin::Ptr{Nothing}, abstol::Missing, reltol::Missing) @ Cuba ~/.julia/packages/Cuba/KIQTB/src/cuhre.jl:97 [14] __solvebp_call(::QuadratureProblem{false, Vector{ForwardDiff.Dual{Nothing, Float64, 2}}, typeof(f), Vector{Float64}, Vector{Float64}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::CubaCuhre, ::Quadrature.ReCallVJP{Quadrature.ZygoteVJP}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{ForwardDiff.Dual{Nothing, Float64, 2}}; reltol::Float64, abstol::Float64, maxiters::Int64, kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}) @ Quadrature ~/.julia/packages/Quadrature/rcCut/src/Quadrature.jl:471 [15] #rrule#43 @ ~/.julia/packages/Quadrature/rcCut/src/Quadrature.jl:493 [inlined] [16] chain_rrule_kw @ ~/.julia/packages/Zygote/0da6K/src/compiler/chainrules.jl:115 [inlined] [17] macro expansion @ ~/.julia/packages/Zygote/0da6K/src/compiler/interface2.jl:0 [inlined] [18] _pullback(::Zygote.Context, ::Quadrature.var"#__solvebp##kw", ::NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}, ::typeof(Quadrature.__solvebp), ::QuadratureProblem{false, Vector{ForwardDiff.Dual{Nothing, Float64, 2}}, typeof(f), Vector{Float64}, Vector{Float64}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::CubaCuhre, ::Quadrature.ReCallVJP{Quadrature.ZygoteVJP}, ::Vector{Float64}, ::Vector{Float64}, ::Vector{ForwardDiff.Dual{Nothing, Float64, 2}}) @ Zygote ~/.julia/packages/Zygote/0da6K/src/compiler/interface2.jl:9 [19] _apply(::Function, ::Vararg{Any, N} where N) @ Core ./boot.jl:804 [20] adjoint @ ~/.julia/packages/Zygote/0da6K/src/lib/lib.jl:200 [inlined] [21] _pullback @ ~/.julia/packages/ZygoteRules/OjfTt/src/adjoint.jl:57 [inlined] [22] _pullback @ ~/.julia/packages/Quadrature/rcCut/src/Quadrature.jl:153 [inlined] [23] _pullback(::Zygote.Context, ::Quadrature.var"##solve#12", ::Quadrature.ReCallVJP{Quadrature.ZygoteVJP}, ::Base.Iterators.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}}, ::typeof(solve), ::QuadratureProblem{false, Vector{ForwardDiff.Dual{Nothing, Float64, 2}}, typeof(f), Vector{Float64}, Vector{Float64}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::CubaCuhre) @ Zygote ~/.julia/packages/Zygote/0da6K/src/compiler/interface2.jl:0 [24] _apply(::Function, ::Vararg{Any, N} where N) @ Core ./boot.jl:804 [25] adjoint @ ~/.julia/packages/Zygote/0da6K/src/lib/lib.jl:200 [inlined] [26] _pullback @ ~/.julia/packages/ZygoteRules/OjfTt/src/adjoint.jl:57 [inlined] [27] _pullback @ ~/.julia/packages/Quadrature/rcCut/src/Quadrature.jl:152 [inlined] [28] _pullback(::Zygote.Context, ::CommonSolve.var"#solve##kw", ::NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}, ::typeof(solve), ::QuadratureProblem{false, Vector{ForwardDiff.Dual{Nothing, Float64, 2}}, typeof(f), Vector{Float64}, Vector{Float64}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::CubaCuhre) @ Zygote ~/.julia/packages/Zygote/0da6K/src/compiler/interface2.jl:0 [29] _pullback @ ./In[1]:9 [inlined] [30] _pullback(ctx::Zygote.Context, f::typeof(testf), args::Vector{ForwardDiff.Dual{Nothing, Float64, 2}}) @ Zygote ~/.julia/packages/Zygote/0da6K/src/compiler/interface2.jl:0 [31] _pullback(f::Function, args::Vector{ForwardDiff.Dual{Nothing, Float64, 2}}) @ Zygote ~/.julia/packages/Zygote/0da6K/src/compiler/interface.jl:34 [32] pullback(f::Function, args::Vector{ForwardDiff.Dual{Nothing, Float64, 2}}) @ Zygote ~/.julia/packages/Zygote/0da6K/src/compiler/interface.jl:40 [33] gradient(f::Function, args::Vector{ForwardDiff.Dual{Nothing, Float64, 2}}) @ Zygote ~/.julia/packages/Zygote/0da6K/src/compiler/interface.jl:58 [34] (::Zygote.var"#77#78"{typeof(testf)})(x::Vector{ForwardDiff.Dual{Nothing, Float64, 2}}) @ Zygote ~/.julia/packages/Zygote/0da6K/src/lib/grad.jl:76 [35] forward_jacobian(f::Zygote.var"#77#78"{typeof(testf)}, x::Vector{Float64}, #unused#::Val{2}) @ Zygote ~/.julia/packages/Zygote/0da6K/src/lib/forward.jl:29 [36] forward_jacobian(f::Function, x::Vector{Float64}) @ Zygote ~/.julia/packages/Zygote/0da6K/src/lib/forward.jl:44 [37] hessian_dual @ ~/.julia/packages/Zygote/0da6K/src/lib/grad.jl:76 [inlined] [38] hessian(f::Function, x::Vector{Float64}) @ Zygote ~/.julia/packages/Zygote/0da6K/src/lib/grad.jl:74 [39] top-level scope @ In[4]:1 [40] eval @ ./boot.jl:360 [inlined] [41] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String) @ Base ./loading.jl:1094

[ChrisRackauckas]

Yes, to fix this someone would need to add a double dual dispatch which moves the Hessian calculation into the integrand. It's almost exactly the same as https://github.com/SciML/Quadrature.jl/blob/v1.9.0/src/Quadrature.jl#L581-L639

[fernando-duarte]

Thank you for your quick response, much appreciated! I think the code you link to will take me a while to understand. I will give it a shot, not sure I am the right person to tackle this.