Forward differentiation and differential equations with solvers like LSODA, Sundials etc

Specific Domains Modelling & Simulations
, , ,
1

I have similar questions to this thread

I do explicit type conversion.
When I use “standard” solver like Tsit5() everything is OK. But with solvers like lsoda() or CVODE_BDF() I’m getting an error:

Example:

using DifferentialEquations, Optim, LinearAlgebra, Sundials

function f(dxdt,x,p,t)
    dxdt[1] = -p[1]*x[1]
    dxdt[2] = p[1]*x[1]
end

tspan = (0.,8.)
tdata = range(0.,stop=8.,length=10)
p0 = [1.]
x0 = [10.,0.]

prob = ODEProblem(f,x0,tspan,p0)

sol = DifferentialEquations.solve(prob; saveat=tdata)

xdata = sol[2,:] + randn(length(sol[2,:])) # in-silico data

function loss_function(x)
    _prob = remake(prob; u0=convert.(eltype(x),prob.u0), p=x)
#    sol_new = DifferentialEquations.solve(_prob, Tsit5(); saveat=tdata)  # OK
    sol_new = DifferentialEquations.solve(_prob, CVODE_BDF(); saveat=tdata)  # Error
    norm(sol_new[2,:] - xdata)^2
end
optimize(loss_function, [2.0], LBFGS(),autodiff=:forward)

Error:

ERROR: LoadError: MethodError: no method matching Sundials.NVector(::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1})
Closest candidates are:
  Sundials.NVector(::Array{Float64,1}) at C:\Users\mazhenir\.julia\packages\Sundials\eAJ5O\src\nvector_wrapper.jl:16
  Sundials.NVector(::Ptr{Sundials._generic_N_Vector}) at C:\Users\mazhenir\.julia\packages\Sundials\eAJ5O\src\nvector_wrapper.jl:24
Stacktrace:
 [1] #__init#20(::Bool, ::Nothing, ::Float64, ::Float64, ::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}, ::Array{Float64,1}, ::Int64, ::Nothing, ::Float64, ::Float64, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Bool, ::Nothing, ::Bool, ::Bool, ::Nothing, ::Bool, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::typeof(DiffEqBase.__init), ::ODEProblem{Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CVODE_BDF{:Newton,:Dense}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}) at C:\Users\mazhenir\.julia\packages\Sundials\eAJ5O\src\common_interface\solve.jl:121
 [2] (::getfield(DiffEqBase, Symbol("#kw##__init")))(::NamedTuple{(:saveat,),Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}, ::typeof(DiffEqBase.__init), ::ODEProblem{Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CVODE_BDF{:Newton,:Dense}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}) at .\none:0
 [3] #init#428(::Base.Iterators.Pairs{Symbol,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Tuple{Symbol},NamedTuple{(:saveat,),Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}}, ::Function, ::ODEProblem{Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CVODE_BDF{:Newton,:Dense}, ::Vararg{Any,N} where N) at C:\Users\mazhenir\.julia\packages\DiffEqBase\ZQVwI\src\solve.jl:20
 [4] (::getfield(DiffEqBase, Symbol("#kw##init")))(::NamedTuple{(:saveat,),Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}, ::typeof(init), ::ODEProblem{Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CVODE_BDF{:Newton,:Dense}, ::Array{Any,1}, ::Vararg{Array{Any,1},N} where N) at .\none:0
 [5] #__solve#19(::Base.Iterators.Pairs{Symbol,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Tuple{Symbol},NamedTuple{(:saveat,),Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}}, ::Function, ::ODEProblem{Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CVODE_BDF{:Newton,:Dense}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at C:\Users\mazhenir\.julia\packages\Sundials\eAJ5O\src\common_interface\solve.jl:10
 [6] (::getfield(DiffEqBase, Symbol("#kw##__solve")))(::NamedTuple{(:saveat,),Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}, ::typeof(DiffEqBase.__solve), ::ODEProblem{Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CVODE_BDF{:Newton,:Dense}, ::Array{Any,1}, ::Array{Any,1}, ::Array{Any,1}, ::Type{Val{true}}) at .\none:0 (repeats 5 times)
 [7] #solve#429(::Base.Iterators.Pairs{Symbol,StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},Tuple{Symbol},NamedTuple{(:saveat,),Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}}, ::Function, ::ODEProblem{Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CVODE_BDF{:Newton,:Dense}) at C:\Users\mazhenir\.julia\packages\DiffEqBase\ZQVwI\src\solve.jl:39
 [8] (::getfield(DiffEqBase, Symbol("#kw##solve")))(::NamedTuple{(:saveat,),Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}, ::typeof(solve), ::ODEProblem{Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},Tuple{Float64,Float64},true,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1},ODEFunction{true,typeof(f),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Nothing,DiffEqBase.StandardODEProblem}, ::CVODE_BDF{:Newton,:Dense}) at .\none:0
 [9] loss_function(::Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1}) at d:\Playground\julia-experiment\rail-switcher\gradient_check.jl:22
 [10] vector_mode_gradient!(::DiffResults.MutableDiffResult{1,Float64,Tuple{Array{Float64,1}}}, ::typeof(loss_function), ::Array{Float64,1}, ::ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1}}) at C:\Users\mazhenir\.julia\packages\ForwardDiff\N0wMF\src\apiutils.jl:37
 [11] gradient! at C:\Users\mazhenir\.julia\packages\ForwardDiff\N0wMF\src\gradient.jl:35 [inlined]
 [12] gradient! at C:\Users\mazhenir\.julia\packages\ForwardDiff\N0wMF\src\gradient.jl:33 [inlined]
 [13] (::getfield(NLSolversBase, Symbol("##14#18")){Float64,typeof(loss_function),ForwardDiff.GradientConfig{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1,Array{ForwardDiff.Dual{ForwardDiff.Tag{typeof(loss_function),Float64},Float64,1},1}}})(::Array{Float64,1}, ::Array{Float64,1}) at C:\Users\mazhenir\.julia\packages\NLSolversBase\KG9Ie\src\objective_types\oncedifferentiable.jl:69
 [14] value_gradient!!(::OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}) at C:\Users\mazhenir\.julia\packages\NLSolversBase\KG9Ie\src\interface.jl:82
 [15] initial_state(::LBFGS{Nothing,LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},getfield(Optim, Symbol("##22#24"))}, ::Optim.Options{Float64,Nothing}, ::OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}) at C:\Users\mazhenir\.julia\packages\Optim\Agd3B\src\multivariate\solvers\first_order\l_bfgs.jl:158
 [16] optimize(::OnceDifferentiable{Float64,Array{Float64,1},Array{Float64,1}}, ::Array{Float64,1}, ::LBFGS{Nothing,LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},getfield(Optim, Symbol("##22#24"))}, ::Optim.Options{Float64,Nothing}) at C:\Users\mazhenir\.julia\packages\Optim\Agd3B\src\multivariate\optimize\optimize.jl:33
 [17] #optimize#87 at C:\Users\mazhenir\.julia\packages\Optim\Agd3B\src\multivariate\optimize\interface.jl:116 [inlined]
 [18] (::getfield(Optim, Symbol("#kw##optimize")))(::NamedTuple{(:autodiff,),Tuple{Symbol}}, ::typeof(optimize), ::Function, ::Array{Float64,1}, ::LBFGS{Nothing,LineSearches.InitialStatic{Float64},LineSearches.HagerZhang{Float64,Base.RefValue{Bool}},getfield(Optim, Symbol("##22#24"))}, ::Optim.Options{Float64,Nothing}) at .\none:0 (repeats 2 times)
 [19] top-level scope at none:0
 [20] include at .\boot.jl:326 [inlined]
 [21] include_relative(::Module, ::String) at .\loading.jl:1038
 [22] include(::Module, ::String) at .\sysimg.jl:29
 [23] include(::String) at .\client.jl:403
 [24] top-level scope at none:0
in expression starting at d:\Playground\julia-experiment\rail-switcher\gradient_check.jl:25
2

This will never happen. It’s a fundamental engineering problem that precompiled binaries from C++/Fortran etc. cannot be autodifferentiated because autodifferentiation always works by intercepting the code itself and creating a new code which calculates the derivative. Cross language automatic differentiation just really isn’t a thing, and in many ways it’s not really possible with JITing C++/Fortran and a lot of other technical hurdles.

(It’s possible to do with code that bridges to R, and theoretically with Python, because both of those languages are interpreted, though any compiled code (like NumPy, SciPy, and Numba/Cython) would cause issues, but then the code likely isn’t too useful since it would be the slowest parts of these languages…)

As an aside, this technical problem is one of the core reasons for building native Julia solvers BTW.

Other forms of differentiation are definitely possible here though. Sensitivity analysis is explained in the docs here:

http://docs.juliadiffeq.org/latest/analysis/sensitivity.html

and a fallback to forward differencing is always a possibility.

1 Like
3

Thank you, @ChrisRackauckas. I will read about “sensitivity analysis”.