Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Hessian throws error #67

Closed
fernando-duarte opened this issue Jun 30, 2021 · 3 comments · Fixed by #190
Closed

Hessian throws error #67

fernando-duarte opened this issue Jun 30, 2021 · 3 comments · Fixed by #190

Comments

@fernando-duarte
Copy link

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
Copy link
Author

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
Copy link
Member

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
Copy link
Author

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

2 participants