Method Error for ForwardDiff.hessian when ForwardDiff.Jacobian and ForwardDiff.gradient works perfectly fine. How to resolve?

New to Julia
,
1

Hello all,
This conversation is in continuation with the conversation I had here. I wanted to find the partial double and mixed derivatives of the variables c1,c2 and c3 with respect to the two variables theta and psi. The first order partial derivatives are working fine ForwardDiff.gradient or ForwardDiff.jacobian, but when I go for the ForwardDiff.hessian, I get a method error.
Below are the concerned functions:

function solve_nlprob3(thetapsi)
       th,psi = thetapsi
        
        l = 130; br = 35; bh = 18; h1 = 15; h2 = 15; br2 = br+90;
        initial_guess = [1.2566177521594696e-17; -0.48551812229559066; 0.485518122295591; -8.31462337379587e-17; -0.48551812229559155; 0.48551812229559127; -7.057740674050417e-17; -8.25000618620076e-18; -0.48551812229559127; 0.48551812229559116; 1.6511444484620922e-16; 2.3617612560724356e-17; 2.4619468236701963e-14; 7.329577692678121e-14; 170.664991614216]
        p = [th.value,psi.value,l,br,bh,h1,h2,br2]
        prob = NonlinearLeastSquaresProblem(NonlinearFunction(objective2, resid_prototype = zeros(24)), initial_guess,p)
        resu = solve(prob,reltol = 1e-12, abstol = 1e-12);
        pose2 = posematrix(th.value,psi.value,resu.u[1],resu.u[13],resu.u[14],resu.u[15],bh)
        phi,th21,th31,th41,th22,th32,th42,th52,th23,th33,th43,th53,px,py,pz=resu.u;
        PX,PY,PZ = pose2[1,4],pose2[2,4],pose2[3,4]        
        
        c1length_analytical=(((-1).*PZ+(-1).*br2.*cos(psi).*sin(th)).^2+((-1).*PX+(-1).*br2.*((-1).*cos(th).*sin(phi)+cos(phi).*sin(psi).*sin(th))).^2+(br2+(-1).*PY+(-1).*br2.*(cos(phi).*cos(th)+sin(phi).*sin(psi).*sin(th))).^2).^(1/2)
        c2length_analytical = (((-1).*PZ+(-1/2).*3 .^(1/2).*br2.*sin(psi)+(1/2).*br2.*cos(psi).*sin(th)).^2+((-1/2).*3 .^(1/2).*br2+(-1).*PX+(1/2).*3 .^(1/2).*br2.*cos(phi).*cos(psi)+(1/2).*br2.*((-1).*cos(th).*sin(phi)+cos(phi).*sin(psi).*sin(th))).^2+((-1/2).*br2+(-1).*PY+(1/2).*3 .^(1/2).*br2.*cos(psi).*sin(phi)+(1/2).*br2.*(cos(phi).*cos(th)+sin(phi).*sin(psi).*sin(th))).^2).^(1/2)
        c3length_analytical = (((-1).*PZ+(1/2).*3 .^(1/2).*br2.*sin(psi)+(1/2).*br2.*cos(psi).*sin(th)).^2+((1/2).*3 .^(1/2).*br2+(-1).*PX+(-1/2).*3 .^(1/2).*br2.*cos(phi).*cos(psi)+(1/2).*br2.*((-1).*cos(th).*sin(phi)+cos(phi).*sin(psi).*sin(th))).^2+((-1/2).*br2+(-1).*PY+(-1/2).*3 .^(1/2).*br2.*cos(psi).*sin(phi)+(1/2).*br2.*(cos(phi).*cos(th)+sin(phi).*sin(psi).*sin(th))).^2).^(1/2)
        #return [c1length_analytical, c2length_analytical, c3length_analytical]
         return c1length_analytical
   
end

gradient_result = ForwardDiff.hessian(solve_nlprob3, [0.0, 0.0])

I am not reproducing the ‘objective2’ function here, as it is too large. However, the general structure is:

function objective2(F,params,p)
    th,psi,l,br,bh,h1,h2,br2 = p
    phi,th21,th31,th41,th22,th32,th42,th52,th23,th33,th43,th53,px,py,pz = params
F[1] =   (h1+l.*cos(th21)+h2.*cos(th21+th31)).^(-1).*cos(th41).*(pz+(-1).*bh.*cos(psi).*cos(th)+(-1).*br.*cos(psi).*sin(th))+(-1).*cos(th21+th31).*(h1+l.*cos(th21)+h2.*cos(th21+th31)).^(-1).*(px+(-1).*bh.*(cos(phi).*cos(th).*sin(psi)+sin(phi).*sin(th))+(-1).*br.*((-1).*cos(th).*sin(phi)+cos(phi).*sin(psi).*sin(th))).*sin(th41) - ( (1/2).*3 .^(1/2).*((1/2).*3 .^(1/2).*(cos(th22+th32).*cos(th52)+sin(th22+th32).*sin(th42).*sin(th52))+(1/2).*((-1).*(h1+l.*cos(th22)).^(-1).*cos(th52).*((-1/2).*px+(1/2).*3 .^(1/2).*py+(1/2).*3 .^(1/2).*br.*((-1/2).*cos(phi).*cos(psi)+(1/2).*3 .^(1/2).*cos(psi).*sin(phi))+((-1).*bh+(-1).*h2).*((1/2).*3 .^(1/2).*(cos(th).*sin(phi).*sin(psi)+(-1).*cos(phi).*sin(th))+(1/2).*((-1).*cos(phi).*cos(th).*sin(psi)+(-1).*sin(phi).*sin(th)))+(1/2).*br.*((1/2).*(cos(th).*sin(phi)+(-1).*cos(phi).*sin(psi).*sin(th))+(1/2).*3 .^(1/2).*(cos(phi).*cos(th)+sin(phi).*sin(psi).*sin(th)))).*sin(th22+th32)+((-1).*(h1+l.*cos(th22)).^(-1).*cos(th42).*(pz+((-1).*bh+(-1).*h2).*cos(psi).*cos(th)+(-1/2).*3 .^(1/2).*br.*sin(psi)+(1/2).*br.*cos(psi).*sin(th))+(h1+l.*cos(th22)).^(-1).*cos(th22+th32).*((-1/2).*px+(1/2).*3 .^(1/2).*py+(1/2).*3 .^(1/2).*br.*((-1/2).*cos(phi).*cos(psi)+(1/2).*3 .^(1/2).*cos(psi).*sin(phi))+((-1).*bh+(-1).*h2).*((1/2).*3 .^(1/2).*(cos(th).*sin(phi).*sin(psi)+(-1).*cos(phi).*sin(th))+(1/2).*((-1).*cos(phi).*cos(th).*sin(psi)+(-1).*sin(phi).*sin(th)))+(1/2).*br.*((1/2).*(cos(th).*sin(phi)+(-1).*cos(phi).*sin(psi).*sin(th))+(1/2).*3 .^(1/2).*(cos(phi).*cos(th)+sin(phi).*sin(psi).*sin(th)))).*sin(th42)).*sin(th52)))+(1/2).*((1/2).*3 .^(1/2).*((-1).*cos(th52).*sin(th22+th32).*sin(th42)+cos(th22+th32).*sin(th52))+(1/2).*((h1+l.*cos(th22)).^(-1).*cos(th42).*cos(th52).*(pz+((-1).*bh+(-1).*h2).*cos(psi).*cos(th)+(-1/2).*3 .^(1/2).*br.*sin(psi)+(1/2).*br.*cos(psi).*sin(th))+(-1).*(h1+l.*cos(th22)).^(-1).*((-1/2).*px+(1/2).*3 .^(1/2).*py+(1/2).*3 .^(1/2).*br.*((-1/2).*cos(phi).*cos(psi)+(1/2).*3 .^(1/2).*cos(psi).*sin(phi))+((-1).*bh+(-1).*h2).*((1/2).*3 .^(1/2).*(cos(th).*sin(phi).*sin(psi)+(-1).*cos(phi).*sin(th))+(1/2).*((-1).*cos(phi).*cos(th).*sin(psi)+(-1).*sin(phi).*sin(th)))+(1/2).*br.*((1/2).*(cos(th).*sin(phi)+(-1).*cos(phi).*sin(psi).*sin(th))+(1/2).*3 .^(1/2).*(cos(phi).*cos(th)+sin(phi).*sin(psi).*sin(th)))).*(sin(th22).*((-1).*cos(th52).*sin(th32).*sin(th42)+cos(th32).*sin(th52))+cos(th22).*(cos(th32).*cos(th52).*sin(th42)+sin(th32).*sin(th52))))))

.......

F[24] = ......
return F
end

The error is:

MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2})

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:207
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:792
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number}
   @ Base char.jl:50
  ...


Stacktrace:
  [1] convert(#unused#::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2})
    @ Base .\number.jl:7
  [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}, i1::Int64)
    @ Base .\array.jl:969
  [3] objective2(F::Vector{Float64}, params::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, p::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}})
    @ Main .\In[93]:71
  [4] (::NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(objective2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Vector{Float64}})(::Vector{Float64}, ::Vararg{Any})
    @ SciMLBase C:\Users\user\.julia\packages\SciMLBase\2HZ5m\src\scimlfunctions.jl:2356
  [5] evaluate_f(prob::NonlinearLeastSquaresProblem{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, true, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(objective2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Vector{Float64}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, u::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}})
    @ NonlinearSolve C:\Users\user\.julia\packages\NonlinearSolve\KlGj2\src\utils.jl:151
  [6] __init(::NonlinearLeastSquaresProblem{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, true, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(objective2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Vector{Float64}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::GaussNewton{nothing, Nothing, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), LineSearch{Nothing, Nothing, Bool}, Nothing}; alias_u0::Bool, maxiters::Int64, abstol::Float64, reltol::Float64, termination_condition::Nothing, internalnorm::typeof(DiffEqBase.NONLINEARSOLVE_DEFAULT_NORM), kwargs::Base.Pairs{Symbol, Bool, Tuple{Symbol, Symbol}, NamedTuple{(:default_set, :second_time), Tuple{Bool, Bool}}})
    @ NonlinearSolve C:\Users\user\.julia\packages\NonlinearSolve\KlGj2\src\gaussnewton.jl:94
  [7] init_call(_prob::NonlinearLeastSquaresProblem{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, true, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(objective2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Vector{Float64}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, args::GaussNewton{nothing, Nothing, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), LineSearch{Nothing, Nothing, Bool}, Nothing}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::Base.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:default_set, :second_time, :reltol, :abstol), Tuple{Bool, Bool, Float64, Float64}}})
    @ DiffEqBase C:\Users\user\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:528
  [8] init_up(prob::NonlinearLeastSquaresProblem{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, true, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(objective2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Vector{Float64}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, sensealg::Nothing, u0::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, p::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, args::GaussNewton{nothing, Nothing, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), LineSearch{Nothing, Nothing, Bool}, Nothing}; kwargs::Base.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:default_set, :second_time, :reltol, :abstol), Tuple{Bool, Bool, Float64, Float64}}})
    @ DiffEqBase C:\Users\user\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:553
  [9] init(prob::NonlinearLeastSquaresProblem{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, true, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(objective2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Vector{Float64}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, args::GaussNewton{nothing, Nothing, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), LineSearch{Nothing, Nothing, Bool}, Nothing}; sensealg::Nothing, u0::Nothing, p::Nothing, kwargs::Base.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:default_set, :second_time, :reltol, :abstol), Tuple{Bool, Bool, Float64, Float64}}})
    @ DiffEqBase C:\Users\user\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:541
 [10] __solve(::NonlinearLeastSquaresProblem{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, true, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(objective2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Vector{Float64}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::GaussNewton{nothing, Nothing, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), LineSearch{Nothing, Nothing, Bool}, Nothing}; kwargs::Base.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:default_set, :second_time, :reltol, :abstol), Tuple{Bool, Bool, Float64, Float64}}})
    @ NonlinearSolve C:\Users\user\.julia\packages\NonlinearSolve\KlGj2\src\NonlinearSolve.jl:134
 [11] macro expansion
    @ C:\Users\user\.julia\packages\NonlinearSolve\KlGj2\src\default.jl:123 [inlined]
 [12] __solve(::NonlinearLeastSquaresProblem{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, true, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(objective2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Vector{Float64}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::NonlinearSolvePolyAlgorithm{:NLLS, 5, Tuple{GaussNewton{nothing, Nothing, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), LineSearch{Nothing, Nothing, Bool}, Nothing}, TrustRegion{nothing, Nothing, Rational{Int64}, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), Rational{Int64}, Nothing}, GaussNewton{nothing, Nothing, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), LineSearch{BackTracking{Float64, Int64}, Nothing, Bool}, Nothing}, TrustRegion{nothing, Nothing, Rational{Int64}, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), Rational{Int64}, Nothing}, LevenbergMarquardt{true, Nothing, Nothing, typeof(NonlinearSolve.DEFAULT_PRECS), Float64, Float64, Float64, Float64, Float64, Float64, Float64}}}; kwargs::Base.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:default_set, :second_time, :reltol, :abstol), Tuple{Bool, Bool, Float64, Float64}}})
    @ NonlinearSolve C:\Users\user\.julia\packages\NonlinearSolve\KlGj2\src\default.jl:115
 [13] __solve(::NonlinearLeastSquaresProblem{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, true, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(objective2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Vector{Float64}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::Nothing; kwargs::Base.Pairs{Symbol, Real, NTuple{4, Symbol}, NamedTuple{(:default_set, :second_time, :reltol, :abstol), Tuple{Bool, Bool, Float64, Float64}}})
    @ NonlinearSolve C:\Users\user\.julia\packages\NonlinearSolve\KlGj2\src\default.jl:398
 [14] __solve(::NonlinearLeastSquaresProblem{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, true, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(objective2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Vector{Float64}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}; default_set::Bool, second_time::Bool, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}})
    @ DiffEqBase C:\Users\user\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:1370
 [15] solve_call(::NonlinearLeastSquaresProblem{Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, true, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(objective2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Vector{Float64}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}})
    @ DiffEqBase C:\Users\user\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:608
 [16] solve_up(::NonlinearLeastSquaresProblem{Vector{Float64}, true, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(objective2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Vector{Float64}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, ::Nothing, ::Vector{Float64}, ::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}; kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}})
    @ DiffEqBase C:\Users\user\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:1049
 [17] solve(::NonlinearLeastSquaresProblem{Vector{Float64}, true, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, NonlinearFunction{true, SciMLBase.FullSpecialize, typeof(objective2), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED_NO_TIME), Nothing, SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}, Vector{Float64}}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{true}, kwargs::Base.Pairs{Symbol, Float64, Tuple{Symbol, Symbol}, NamedTuple{(:reltol, :abstol), Tuple{Float64, Float64}}})
    @ DiffEqBase C:\Users\user\.julia\packages\DiffEqBase\eTCPy\src\solve.jl:980
 [18] solve_nlprob3(thetapsi::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}, 2}})
    @ Main .\In[121]:8
 [19] vector_mode_dual_eval!
    @ C:\Users\user\.julia\packages\ForwardDiff\PcZ48\src\apiutils.jl:24 [inlined]
 [20] vector_mode_gradient(f::typeof(solve_nlprob3), x::Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}, cfg::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}, 2}}})
    @ ForwardDiff C:\Users\user\.julia\packages\ForwardDiff\PcZ48\src\gradient.jl:89
 [21] gradient
    @ C:\Users\user\.julia\packages\ForwardDiff\PcZ48\src\gradient.jl:19 [inlined]
 [22] #99
    @ C:\Users\user\.julia\packages\ForwardDiff\PcZ48\src\hessian.jl:16 [inlined]
 [23] vector_mode_dual_eval!
    @ C:\Users\user\.julia\packages\ForwardDiff\PcZ48\src\apiutils.jl:24 [inlined]
 [24] vector_mode_jacobian(f::ForwardDiff.var"#99#100"{typeof(solve_nlprob3), ForwardDiff.HessianConfig{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}, 2}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}}}, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}})
    @ ForwardDiff C:\Users\user\.julia\packages\ForwardDiff\PcZ48\src\jacobian.jl:125
 [25] jacobian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.JacobianConfig{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}}, ::Val{false})
    @ ForwardDiff C:\Users\user\.julia\packages\ForwardDiff\PcZ48\src\jacobian.jl:21
 [26] hessian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.HessianConfig{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}, 2}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}}, ::Val{true})
    @ ForwardDiff C:\Users\user\.julia\packages\ForwardDiff\PcZ48\src\hessian.jl:17
 [27] hessian(f::Function, x::Vector{Float64}, cfg::ForwardDiff.HessianConfig{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}, 2}}, Vector{ForwardDiff.Dual{ForwardDiff.Tag{typeof(solve_nlprob3), Float64}, Float64, 2}}})
    @ ForwardDiff C:\Users\user\.julia\packages\ForwardDiff\PcZ48\src\hessian.jl:15
 [28] hessian(f::Function, x::Vector{Float64})
    @ ForwardDiff C:\Users\user\.julia\packages\ForwardDiff\PcZ48\src\hessian.jl:15
 [29] top-level scope
    @ In[127]:1

I would deeply appreciate if any of you can guide me on decoding the error message.

2

Presumably, you’re forcing internal variables to have type Float64 somewhere, whereas ForwardDiff.jl needs to propagate numbers of type Dual through your code.

Can you show us this line of your code? Apparently you initialized a vector with element type Float64, which means it cannot accept Dual numbers as coefficients. This typically happens when you do A = zeros(n) instead of A = zeros(eltype(myvariable), n).

1 Like
3

Hello @gdalle , thank you for your reply. In the following line of code from the function solve_nlprob3, I think I may have initialized the output of the objective2 function as zeros(24). Could this be the problem? I haven’t initialized the variable A anywhere else.
prob = NonlinearLeastSquaresProblem(NonlinearFunction(objective2, resid_prototype = zeros(24)), initial_guess,p)

4

My bad, the variable isn’t necessarily named A, but it is a vector of Float64 so resid_prototype = zeros(24) certainly fits the bill. However, since we’re going into the internals of NonlinearSolve.jl, I’d rather let @avikpal take over.

1 Like
5

NonlinearSolve.jl specializes ForwardDiff duals, but only one level. It needs some code to handle the Hessian the same way. Open an issue.

1 Like
6

Yes, open an issue. We added the dispatches to SimpleNonlinearSolve, but the same hasn’t been done to NonlinearSolve. NonlinearSolve.jl/src/internal/forward_diff.jl at master · SciML/NonlinearSolve.jl · GitHub needs to be expanded to NLLS Problems

7

Hello @avikpal, thanks for your suggestions. I have opened an issue as recommended. Mean while, are there any good work arounds to obtain the hessian, based on what I have available from ForwardDiff.gradient? I am thinking of using Finite differences.

8

You can use the solvers from SimpleNonlinearSolve like SimpleGaussNewton, that should just work

9

Thank you @avikpal for your reply. I have tried to use SimpleGaussNewton() solver, but it is still returning a method error for ForwardDiff.hessian, while it is working fine for ForwardDiff.gradient.
I am appending a MWE for your reference, where I have used the SimpleGaussNewton() solver.

function objfn!(F, init, params)
    th1, th2 = init
    px, py, l1, l2 = params
    F[1] = l1 * cos(th1) + l2 * cos(th1 + th2) - px
    F[2] = l1 * sin(th1) + l2 * sin(th1 + th2) - py
    return F
end


function solve_nlprob(pxpy)
    px,py = pxpy
theta1 = pi/4;
theta2 = pi/4
initial_guess = [theta1;theta2];
l1 = 60; l2 = 60;
p = [px.value;py.value;l1;l2];
prob = NonlinearLeastSquaresProblem(NonlinearFunction(objfn!, resid_prototype = zeros(2)), initial_guess,p)
resu = solve(prob,GaussNewton());
th1,th2 = resu.u
    println(th1)
    println(th2)
cable1_base = [-90;0;0];
cable2_base = [-150;0;0];
cable3_base = [150;0;0];
cable1_top = [l1*cos(th1)/2;l1*sin(th1)/2;0];
cable23_top = [l1*cos(th1)+l2*cos(th1+th2)/2;l1*sin(th1)+l2*sin(th1+th2)/2;0];
c1_length=sqrt((cable1_top[1]-cable1_base[1])^2 + (cable1_top[2]-cable1_base[2])^2)
c2_length=sqrt((cable23_top[1]-cable2_base[1])^2 + (cable23_top[2]-cable2_base[2])^2)
c3_length=sqrt((cable23_top[1]-cable3_base[1])^2 + (cable23_top[2]-cable3_base[2])^2)
    return c1_length
end

ForwardDiff.hessian(solve_nlprob,[34.0,87.0])
10

Can you paste the stacktrace?

11

@gdalle, Apologies. Instead of SimpleGaussNewton(), I have used GaussNewton() as the solver. I was able to obtain the hessian while using the SimpleGaussNewton() solver, as suggested by @avikpal. My current objective is to compare the results of the ForwardDiff.hessian, with a finite difference scheme.

12

In my original post, inside the function solve_nlprob3(thetapsi), I had to extract the values from thetapsi as:

th,psi = thetapsi
th = th.value
psi = psi.value

I didn’t have to use this after using the solver SimpleNonlinearSolve() as suggested by @avikpal, as there are no more MethodErrors.
ForwardDiff.gradient is fine, but ForwardDiff.hessian returns a StackOverflowError: without a stacktrace.
Are they any way around this? I thought of investigating chunk size but I am not sure as my input vector has only two elements.

13

As mentioned earlier, it’s just not implemented. The workaround is to not use it, or help us get the SimpleNonlinearSolve.jl one updated to handle the whole NonlinearSolve.jl

14

In the meantime, perhaps ForwardDiff over Zygote would work? It’s more efficient for the Hessian anyway.
You can try it with

using DifferentiationInterface
import ForwardDiff, Zygote

hessian(f, SecondOrder(AutoForwardDiff(), AutoZygote()), x)
15

Yes, the user should probably do that, and we should add overloads. Both are true.

1 Like
16

All of you, I deeply appreciate your persistent insights. I created a new environment and installed the package DifferentiationInterface.jl. But when I tried to install the package NonlinearSolve.jl, I run into a comptability issue, which I believe is due to the FiniteDiff not being compatible with DifferentiationInterface. My current Julia version is 1.9.3. Let me check if I have any luck after I update Julia to the latest version.

ERROR: Unsatisfiable requirements detected for package FiniteDiff [6a86dc24]:
 FiniteDiff [6a86dc24] log:
 ├─possible versions are: 2.0.0-2.23.1 or uninstalled
 ├─restricted by compatibility requirements with DifferentiationInterface [a0c0ee7d] to versions: 2.23.0-2.23.1 or uninstalled
 │ └─DifferentiationInterface [a0c0ee7d] log:
 │   ├─possible versions are: 0.1.0-0.5.4 or uninstalled
 │   ├─restricted to versions * by an explicit requirement, leaving only versions: 0.1.0-0.5.4
 │   └─restricted by julia compatibility requirements to versions: 0.3.4-0.5.4 or uninstalled, leaving only versions: 0.3.4-0.5.4
 ├─restricted by julia compatibility requirements to versions: 2.0.0-2.22.0 or uninstalled, leaving only versions: uninstalled
17

You should indeed update to Julia 1.10, where all autodiff backends are tested by DifferentiationInterface.
FiniteDiff, like other parts of the (extended) SciML organization, has dropped support for Julia versions below 1.10 a while ago, so you’re stuck with a very old version if you work on 1.9, which was released before DifferentiationInterface came along.

2 Likes
18

I have updated Julia to 1.10.4. Upon trying:

using DifferentiationInterface, SciMLSensitivity
import ForwardDiff, Zygote
hessian(f, SecondOrder(AutoForwardDiff(), AutoZygote()), [34.0,87.0])

I get the following message:

Compatibility with reverse-mode automatic differentiation requires SciMLSensitivity.jl.
Please install SciMLSensitivity.jl and do `using SciMLSensitivity`/`import SciMLSensitivity`
for this functionality. For more details, see https://sensitivity.sciml.ai/dev/.

I have separately installed SciMLSensitivity package at the beginning of the project. The output of

Pkg.status()

is:

a0c0ee7d] DifferentiationInterface v0.4.2
  [f6369f11] ForwardDiff v0.10.36
  [8913a72c] NonlinearSolve v3.12.4
  [1ed8b502] SciMLSensitivity v7.61.0
  [727e6d20] SimpleNonlinearSolve v1.9.0
  [e88e6eb3] Zygote v0.6.70
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`
19

NonlinearLeastSquares is a very recent development it doesn’t have Reverse Mode bindings yet. It is on my todo list for post summer unless anyone else does it.

1 Like