# Why the separation of \`ODEProblem\` and \`solve\` in DifferentialEquations.jl?

**URL:** <https://discourse.julialang.org/t/why-the-separation-of-odeproblem-and-solve-in-differentialequations-jl/43737>\
**Category:** General Usage\
**Tags:** diffeq\
**Created:** [July 27, 2020, 2:20am UTC](https://discourse.julialang.org/t/why-the-separation-of-odeproblem-and-solve-in-differentialequations-jl/43737 "2020-07-27T02:20:19Z")\
**Posts on this page:** 14\
**Page:** 1

<div class="post-metadata">

**Author:** ![Samuel\_Ainsworth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samuel_ainsworth/32/15644_2.png) [@Samuel\_Ainsworth](https://discourse.julialang.org/u/Samuel_Ainsworth)\
**Post date:** [July 27, 2020, 2:20am UTC](https://discourse.julialang.org/t/why-the-separation-of-odeproblem-and-solve-in-differentialequations-jl/43737/1 "2020-07-27T02:20:19Z")

</div>

Is there a performance benefit to reusing `ODEProblem`s? I’m curious why the `ODEProblem` struct/function if the equivalent info could also just be handed off to `solve`?

I have a problem where I repeatedly need to solve the same large ODE with slightly different conditions repeatedly, and I’d like to avoid reuse work as much as possible if possible. I’m currently putting a lot of pressure on the GC by repeatedly creating `ODEProblem`s and then `solve`ing them. Is there a perf trick that I’m missing?

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 27, 2020, 3:58am UTC](https://discourse.julialang.org/t/why-the-separation-of-odeproblem-and-solve-in-differentialequations-jl/43737/2 "2020-07-27T03:58:52Z")

</div>

> [@Samuel\_Ainsworth](#):
>
> Is there a performance benefit to reusing `ODEProblem` s? I’m curious why the `ODEProblem` struct/function if the equivalent info could also just be handed off to `solve` ?

We need this information in order to know how to dispatch between ODEs, SDEs, DAEs, DDEs, jump equations, etc. In some cases, a function, initial condition, and timespan are not enough to distinguish between these cases. It also let’s us to a lot of engineering and introspection at a higher level. See the coming JuliaCon talk on automated optimization and parallelism which relies on introspection of the problem types pre-solve. Would could have people paste the same 20 arguments in multiple points of their code, but seems decidedly against DRY.

> [@Samuel\_Ainsworth](#):
>
> I have a problem where I repeatedly need to solve the same large ODE with slightly different conditions repeatedly, and I’d like to avoid reuse work as much as possible if possible. I’m currently putting a lot of pressure on the GC by repeatedly creating `ODEProblem` s and then `solve` ing them. Is there a perf trick that I’m missing?

The first thing to note is that if you are creating thousands of really small ODEProblems and want to avoid the dynamic checking going on, then you should do `ODEProblem{false}(...)` or `ODEProblem{true}(...)`, i.e. directly declare whether it’s in in-place or out-of-place form. Normally this isn’t a huge deal so it’s not mentioned very often in the documentation, but it’s in there since there are scenarios where this helps.

Secondly, note that there are a lot of optimizations in Julia v1.5 (coming out next week) and Julia v1.6 which are specifically designed to optimize this case. Most notably,

> <https://github.com/JuliaLang/julia/blob/v1.5.0-rc1/NEWS.md#compilerruntime-improvements>

> - Immutable structs (including tuples) that contain references can now be allocated on the stack, and allocated inline within arrays and other structs ([#33886](https://github.com/JuliaLang/julia/issues/33886)). This significantly reduces the number of heap allocations in some workloads. Code that requires assumptions about object layout and addresses (usually for interoperability with C or other languages) might need to be updated; for example any object that needs a stable address should be a `mutable struct` . As a result, Array `view` s no longer allocate ([#34126]).

On v1.5 the `ODEProblem`, since it is [a fully typed immutable struct](https://github.com/SciML/DiffEqBase.jl/blob/v6.40.9/src/problems/ode_problems.jl#L16-L40), should be completely elided on Julia v1.5. So that issue should completely disappear.

That’s the issue that you know about. Let me just quickly mention though that’s not the real issue: this is only 80 bytes and while annoying it’s not the allocation you’re actually worried about. The one you’re actually worried about a little bit more insidious is:

> <https://github.com/SciML/DifferentialEquations.jl/issues/603>
>
> Hi,
> 
> could you please try to help understand why this minimal working example …does not return an object of concrete type?
> \`\`\`julia
> function f(du,u,p,t)  
> du\[1\] = p.a
> du\[2\] = p.b
> end
> 
> function solve\_ode(f::F, p::P) where {F,P}
> 
> tspan = (0., 1.0)
> Δt = tspan\[2\] - tspan\[1\]
> dt = 1/252
> nodes = Int(ceil(Δt / dt) + 1)
> t = T = \[tspan\[1\] + (i - 1) \* dt for i = 1:nodes\]
> 
> # if I do not set {true}, prob type Any...
> prob = ODEProblem{true}(f, \[0., 0.\], tspan, p)
> # prob = ODEProblem(f, \[0., 0.\], tspan, p)
> 
> prob\_func = (prob, i, repeat) -\> begin
> remake(prob, tspan = (T\[i + 1\], t\[1\]))
> end
> 
> # ensemble problem
> odes = EnsembleProblem(prob, prob\_func = prob\_func)
> 
> sol = DifferentialEquations.solve(
> odes, Tsit5(), EnsembleThreads(), trajectories = nodes - 1, saveat = -dt
> )
> 
> return sol
> end
> \`\`\`
> When running the following:
> \`\`\`julia
> @code\_warntype solve\_ode(f, (a = 1, b = 1)
> Variables
> #self#::Core.Compiler.Const(solve\_ode, false)
> f::Core.Compiler.Const(f, false)
> p::NamedTuple{(:a, :b),Tuple{Int64,Int64}}
> #25::var"#25#27"{Tuple{Float64,Float64},Float64}
> #26::var"#26#28"{Array{Float64,1},Array{Float64,1}}
> tspan::Tuple{Float64,Float64}
> Δt::Float64
> dt::Float64
> nodes::Int64
> T::Array{Float64,1}
> t::Array{Float64,1}
> prob::ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,NamedTuple{(:a, :b),Tuple{Int64,Int64}},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}
> prob\_func::var"#26#28"{Array{Float64,1},Array{Float64,1}}
> odes::EnsembleProblem{ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,NamedTuple{(:a, :b),Tuple{Int64,Int64}},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},var"#26#28"{Array{Float64,1},Array{Float64,1}},DiffEqBase.var"#342#348",DiffEqBase.var"#344#350",Array{Any,1}}
> sol::EnsembleSolution{\_A,\_B,\_C} where \_C where \_B where \_A
> 
> Body::EnsembleSolution{\_A,\_B,\_C} where \_C where \_B where \_A
> 1 ─ (tspan = Core.tuple(0.0, 1.0))
> │ %2 = Base.getindex(tspan::Core.Compiler.Const((0.0, 1.0), false), 2)::Core.Compiler.Const(1.0, false)
> │ %3 = Base.getindex(tspan::Core.Compiler.Const((0.0, 1.0), false), 1)::Core.Compiler.Const(0.0, false)
> │ (Δt = %2 - %3)
> │ (dt = 1 / 252)
> │ %6 = (Δt::Core.Compiler.Const(1.0, false) / dt::Core.Compiler.Const(0.003968253968253968, false))::Core.Compiler.Const(252.0, false)
> │ %7 = Main.ceil(%6)::Core.Compiler.Const(252.0, false)
> │ %8 = (%7 + 1)::Core.Compiler.Const(253.0, false)
> │ (nodes = Main.Int(%8))
> │ %10 = Main.:(var"#25#27")::Core.Compiler.Const(var"#25#27", false)
> │ %11 = Core.typeof(tspan::Core.Compiler.Const((0.0, 1.0), false))::Core.Compiler.Const(Tuple{Float64,Float64}, false)
> │ %12 = Core.typeof(dt::Core.Compiler.Const(0.003968253968253968, false))::Core.Compiler.Const(Float64, false)
> │ %13 = Core.apply\_type(%10, %11, %12)::Core.Compiler.Const(var"#25#27"{Tuple{Float64,Float64},Float64}, false)
> │ %14 = tspan::Core.Compiler.Const((0.0, 1.0), false)::Core.Compiler.Const((0.0, 1.0), false)
> │ (#25 = %new(%13, %14, dt::Core.Compiler.Const(0.003968253968253968, false)))
> │ %16 = #25::Core.Compiler.Const(var"#25#27"{Tuple{Float64,Float64},Float64}((0.0, 1.0), 0.003968253968253968), false)::Core.Compiler.Const(var"#25#27"{Tuple{Float64,Float64},Float64}((0.0, 1.0), 0.003968253968253968), false)
> │ %17 = (1:nodes::Core.Compiler.Const(253, false))::Core.Compiler.Const(1:253, false)
> │ %18 = Base.Generator(%16, %17)::Core.Compiler.Const(Base.Generator{UnitRange{Int64},var"#25#27"{Tuple{Float64,Float64},Float64}}(var"#25#27"{Tuple{Float64,Float64},Float64}((0.0, 1.0), 0.003968253968253968), 1:253), false)
> │ %19 = Base.collect(%18)::Array{Float64,1}
> │ (T = %19)
> │ (t = %19)
> │ %22 = Core.apply\_type(Main.ODEProblem, true)::Core.Compiler.Const(ODEProblem{true,tType,isinplace,P,F,K,PT} where PT where K where F where P where isinplace where tType, false)
> │ %23 = Base.vect(0.0, 0.0)::Array{Float64,1}
> │ %24 = tspan::Core.Compiler.Const((0.0, 1.0), false)::Core.Compiler.Const((0.0, 1.0), false)
> │ (prob = (%22)(f, %23, %24, p))
> │ %26 = Main.:(var"#26#28")::Core.Compiler.Const(var"#26#28", false)
> │ %27 = Core.typeof(T)::Core.Compiler.Const(Array{Float64,1}, false)
> │ %28 = Core.typeof(t)::Core.Compiler.Const(Array{Float64,1}, false)
> │ %29 = Core.apply\_type(%26, %27, %28)::Core.Compiler.Const(var"#26#28"{Array{Float64,1},Array{Float64,1}}, false)
> │ %30 = T::Array{Float64,1}
> │ (#26 = %new(%29, %30, t))
> │ (prob\_func = #26)
> │ %33 = (:prob\_func,)::Core.Compiler.Const((:prob\_func,), false)
> │ %34 = Core.apply\_type(Core.NamedTuple, %33)::Core.Compiler.Const(NamedTuple{(:prob\_func,),T} where T\<:Tuple, false)
> │ %35 = Core.tuple(prob\_func)::Tuple{var"#26#28"{Array{Float64,1},Array{Float64,1}}}
> │ %36 = (%34)(%35)::NamedTuple{(:prob\_func,),Tuple{var"#26#28"{Array{Float64,1},Array{Float64,1}}}}
> │ %37 = Core.kwfunc(Main.EnsembleProblem)::Core.Compiler.Const(Core.var"#Type##kw"(), false)
> │ (odes = (%37)(%36, Main.EnsembleProblem, prob::Core.Compiler.PartialStruct(ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,NamedTuple{(:a, :b),Tuple{Int64,Int64}},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, Any\[Core.Compiler.Const(ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing}(f, LinearAlgebra.UniformScaling{Bool}(true), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing), false), Array{Float64,1}, Core.Compiler.Const((0.0, 1.0), false), NamedTuple{(:a, :b),Tuple{Int64,Int64}}, Core.Compiler.Const(Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}(), false), Core.Compiler.Const(DiffEqBase.StandardODEProblem(), false)\])))
> │ %39 = Main.Tsit5()::Core.Compiler.Const(Tsit5(), false)
> │ %40 = Main.EnsembleThreads()::Core.Compiler.Const(EnsembleThreads(), false)
> │ %41 = (nodes::Core.Compiler.Const(253, false) - 1)::Core.Compiler.Const(252, false)
> │ %42 = -dt::Core.Compiler.Const(0.003968253968253968, false)::Core.Compiler.Const(-0.003968253968253968, false)
> │ %43 = (:trajectories, :saveat)::Core.Compiler.Const((:trajectories, :saveat), false)
> │ %44 = Core.apply\_type(Core.NamedTuple, %43)::Core.Compiler.Const(NamedTuple{(:trajectories, :saveat),T} where T\<:Tuple, false)
> │ %45 = Core.tuple(%41, %42)::Core.Compiler.Const((252, -0.003968253968253968), false)
> │ %46 = (%44)(%45)::NamedTuple{(:trajectories, :saveat),Tuple{Int64,Float64}}
> │ %47 = DifferentialEquations.solve::Core.Compiler.Const(DiffEqBase.solve, false)
> │ %48 = Core.kwfunc(%47)::Core.Compiler.Const(DiffEqBase.var"#solve##kw"(), false)
> │ %49 = DifferentialEquations.solve::Core.Compiler.Const(DiffEqBase.solve, false)
> │ %50 = odes::Core.Compiler.PartialStruct(EnsembleProblem{ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,NamedTuple{(:a, :b),Tuple{Int64,Int64}},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},var"#26#28"{Array{Float64,1},Array{Float64,1}},DiffEqBase.var"#342#348",DiffEqBase.var"#344#350",Array{Any,1}}, Any\[Core.Compiler.PartialStruct(ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,NamedTuple{(:a, :b),Tuple{Int64,Int64}},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, Any\[Core.Compiler.Const(ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing}(f, LinearAlgebra.UniformScaling{Bool}(true), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing), false), Array{Float64,1}, Core.Compiler.Const((0.0, 1.0), false), NamedTuple{(:a, :b),Tuple{Int64,Int64}}, Core.Compiler.Const(Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}(), false), Core.Compiler.Const(DiffEqBase.StandardODEProblem(), false)\]), var"#26#28"{Array{Float64,1},Array{Float64,1}}, Core.Compiler.Const(DiffEqBase.var"#342#348"(), false), Core.Compiler.Const(DiffEqBase.var"#344#350"(), false), Array{Any,1}\])::Core.Compiler.PartialStruct(EnsembleProblem{ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,NamedTuple{(:a, :b),Tuple{Int64,Int64}},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},var"#26#28"{Array{Float64,1},Array{Float64,1}},DiffEqBase.var"#342#348",DiffEqBase.var"#344#350",Array{Any,1}}, Any\[Core.Compiler.PartialStruct(ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,NamedTuple{(:a, :b),Tuple{Int64,Int64}},ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem}, Any\[Core.Compiler.Const(ODEFunction{true,typeof(f),LinearAlgebra.UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing}(f, LinearAlgebra.UniformScaling{Bool}(true), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing), false), Array{Float64,1}, Core.Compiler.Const((0.0, 1.0), false), NamedTuple{(:a, :b),Tuple{Int64,Int64}}, Core.Compiler.Const(Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}(), false), Core.Compiler.Const(DiffEqBase.StandardODEProblem(), false)\]), var"#26#28"{Array{Float64,1},Array{Float64,1}}, Core.Compiler.Const(DiffEqBase.var"#342#348"(), false), Core.Compiler.Const(DiffEqBase.var"#344#350"(), false), Array{Any,1}\])
> │ (sol = (%48)(%46, %49, %50, %39, %40))
> └── return sol
> \`\`\`
> the returned value type is not concrete. Or there is something I am not understanding?
> 
> Thanks!

Essentially previous versions of Julia were unable to rely on constant values inside of keyword arguments, and so because you could do things like `save_idxs=1` to change the output from an array to a scalar of just the first value, the existence of these features caused some instabilities which you have to be careful about. We were able to prove this was the case by isolated every issue and get it statically compiling on Julia v1.4, but it required removing a few keyword arguments which we noted were all linked to this idea.

This was fixed in [more precise inference of `splatnew` by JeffBezanson · Pull Request #35976 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/pull/35976) and it was the last thing added to the v1.5 backports: [Backports for Julia 1.5-RC1 (or beta2) by KristofferC · Pull Request #36098 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/pull/36098) so it should be fixed in the RC1. So please test the RC and see if that’s all fixed up. If that doesn’t reduce the last small allocation, I have a hook:

> <https://github.com/SciML/DiffEqBase.jl/blob/master/src/solve.jl#L75-L82>

which would allow for specializing the compiled output, but it requires these changes in v1.5 to handle all of the keyword arguments (otherwise it’ll just error outputs which don’t match the default).

That said, one person did find a case which wasn’t specialized by this: [Add keyword argument constant propogation to News by ChrisRackauckas · Pull Request #36292 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/pull/36292#issuecomment-652055307) but we can make sure to open up an issue so this gets fixed in Julia v1.6 if that effects DiffEq users.

This sounds like a good JuliaCon time issue though: let’s get this all cleaned up.

---

<div class="post-metadata">

**Author:** ![Samuel\_Ainsworth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samuel_ainsworth/32/15644_2.png) [@Samuel\_Ainsworth](https://discourse.julialang.org/u/Samuel_Ainsworth)\
**Post date:** [July 27, 2020, 5:16am UTC](https://discourse.julialang.org/t/why-the-separation-of-odeproblem-and-solve-in-differentialequations-jl/43737/3 "2020-07-27T05:16:12Z")

</div>

Thank you for such a thorough answer, Chris! I’ll start digging through each of these references further. Looking forward to v1.5!

After some more detailed profiling, I found that in fact performance of the `solve` on the problem is fast; it’s actually the `solve` on the `ODEAdjointProblem` that is much more problematic. Here are some numbers on a very small, toy example:

```julia
forward: 0.000738 seconds (4.65 k allocations: 958.266 KiB)
adjoint: 0.007064 seconds (98.62 k allocations: 5.837 MiB)

```

This is just running the adjoint process on the “u” variables, no parameters yet.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 27, 2020, 5:32am UTC](https://discourse.julialang.org/t/why-the-separation-of-odeproblem-and-solve-in-differentialequations-jl/43737/4 "2020-07-27T05:32:22Z")

</div>

Can I get that example? Reverse mode will always need to allocate unless it’s [a reversible AD](https://github.com/GiggleLiu/NiLang.jl) simply because it needs to cache values for the reverse pass, but I don’t know if it should be that high. I can take a look and profile it some and see if there’s any major offender to fix.

Note that if calculations are small enough to consider the allocations of the ODEProblem itself, maybe you should be using forward sensitivities?

---

<div class="post-metadata">

**Author:** ![Samuel\_Ainsworth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samuel_ainsworth/32/15644_2.png) [@Samuel\_Ainsworth](https://discourse.julialang.org/u/Samuel_Ainsworth)\
**Post date:** [July 27, 2020, 9:54am UTC](https://discourse.julialang.org/t/why-the-separation-of-odeproblem-and-solve-in-differentialequations-jl/43737/5 "2020-07-27T09:54:33Z")

</div>

Yeah, absolutely. I’m super curious to figure out what’s going on! Here’s a small MWE: [https://gist.github.com/samuela/8d7cf55cf921decfffc7559691bfad12](https://gist.github.com/samuela/8d7cf55cf921decfffc7559691bfad12).

On this toy problem obv everything is already very fast, but in more realistic scenarios I’m seeing 10-100x slowdowns. Forward solves that take less than 100ms, and corresponding adjoint solves that take 10s. That’s on a 12-dimensional quadrotor system defined in much the same way. Happy to point you towards some example code for that as well if you’re interested!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 27, 2020, 10:14am UTC](https://discourse.julialang.org/t/why-the-separation-of-odeproblem-and-solve-in-differentialequations-jl/43737/6 "2020-07-27T10:14:41Z")

</div>

[https://diffeq.sciml.ai/latest/analysis/sensitivity/](https://diffeq.sciml.ai/latest/analysis/sensitivity/)

> `TrackerAdjoint` is able to use a `TrackedArray` form with out-of-place functions `du = f(u,p,t)` but requires an `Array{TrackedReal}` form for `f(du,u,p,t)` mutating `du` . The latter has much more overhead, and should be avoided if possible. Thus if solving non-ODEs with lots of parameters, using `TrackerAdjoint` with an out-of-place definition may be the current best option.

The same is true about the vjps, though if you don’t have branching you can `ReverseDiffVJP(true)` to compile the backpass. But I think your problem is probably much better off doing

```julia
function aug_dynamics!(z, policy_params, t)
    x = @view z[2:end]
    u = policy(x, policy_params)
    [x' * x + u' * u;u]
end

```

to accommodate for those factors in reverse mode AD (or use the AoS-\>SoA conversion and convert back).

---

<div class="post-metadata">

**Author:** ![Samuel\_Ainsworth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samuel_ainsworth/32/15644_2.png) [@Samuel\_Ainsworth](https://discourse.julialang.org/u/Samuel_Ainsworth)\
**Post date:** [July 27, 2020, 7:50pm UTC](https://discourse.julialang.org/t/why-the-separation-of-odeproblem-and-solve-in-differentialequations-jl/43737/7 "2020-07-27T19:50:21Z")

</div>

Gotcha, I misunderstood in-place dynamics as always being faster than out-of-place ones. I’m getting ~2x speedup in the adjoint passes now with the OOP version:

```julia
forward: 0.001867 seconds (9.42 k allocations: 7.273 MiB)
adjoint: 0.020280 seconds (45.02 k allocations: 43.124 MiB, 38.89% gc time)

```

(x\_dim = 64 this time.) But it’s still ~10x slower. Is that normal? It still seems a bit high since the adjoint process is the same size as the forward with nearly identical dynamics.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 27, 2020, 7:57pm UTC](https://discourse.julialang.org/t/why-the-separation-of-odeproblem-and-solve-in-differentialequations-jl/43737/8 "2020-07-27T19:57:20Z")

</div>

> [@Samuel\_Ainsworth](#):
>
> I misunderstood in-place dynamics as always being faster than out-of-place ones

That’s always true for forward passes, but with reverse passes it can sometimes mix better with the vjp calculation to do it OOP. That should go away with the new AD framework that’s being worked on though IIRC.

> [@Samuel\_Ainsworth](#):
>
> But it’s still ~10x slower. Is that normal? It still seems a bit high since the adjoint process is the same size as the forward with nearly identical dynamics.

That does still seem high. I’ll want to take a deep look at that example.

---

<div class="post-metadata">

**Author:** ![Samuel\_Ainsworth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samuel_ainsworth/32/15644_2.png) [@Samuel\_Ainsworth](https://discourse.julialang.org/u/Samuel_Ainsworth)\
**Post date:** [July 27, 2020, 8:00pm UTC](https://discourse.julialang.org/t/why-the-separation-of-odeproblem-and-solve-in-differentialequations-jl/43737/9 "2020-07-27T20:00:52Z")

</div>

Ok, I just update [the gist](https://gist.github.com/samuela/8d7cf55cf921decfffc7559691bfad12) to OOP dynamics and for clarity!

Here’s the `@btime` results:

```julia
[ Info: forward
  1.443 ms (9420 allocations: 7.27 MiB)
[ Info: BacksolveAdjoint
  50.525 ms (115033 allocations: 113.32 MiB)
[ Info: InterpolatingAdjoint
  30.859 ms (81025 allocations: 75.65 MiB)
[ Info: QuadratureAdjoint
  21.833 ms (50868 allocations: 48.77 MiB)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 27, 2020, 8:08pm UTC](https://discourse.julialang.org/t/why-the-separation-of-odeproblem-and-solve-in-differentialequations-jl/43737/10 "2020-07-27T20:08:34Z")

</div>

Perfect, thanks.

---

<div class="post-metadata">

**Author:** ![Samuel\_Ainsworth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/samuel_ainsworth/32/15644_2.png) [@Samuel\_Ainsworth](https://discourse.julialang.org/u/Samuel_Ainsworth)\
**Post date:** [July 27, 2020, 8:45pm UTC](https://discourse.julialang.org/t/why-the-separation-of-odeproblem-and-solve-in-differentialequations-jl/43737/11 "2020-07-27T20:45:04Z")

</div>

In case you’re curious, here’s the actual problem I’m trying to speed up atm: [https://gist.github.com/samuela/9a1daca41fd46ce5e67f5df150933373](https://gist.github.com/samuela/9a1daca41fd46ce5e67f5df150933373).

Curiously I’ve found that in-place is actually quite a bit faster than out-of-place for this system:

```julia
In-place:
[ Info: forward
  20.009 ms (23413 allocations: 58.50 MiB)
[ Info: BacksolveAdjoint
  5.056 s (28059882 allocations: 5.93 GiB)
[ Info: InterpolatingAdjoint
  1.596 s (9013271 allocations: 1.89 GiB)
[ Info: QuadratureAdjoint
  473.475 ms (2782989 allocations: 599.98 MiB)
  
Out-of-place:
[ Info: forward
  22.901 ms (45777 allocations: 60.26 MiB)
[ Info: BacksolveAdjoint
  7.333 s (25147876 allocations: 3.42 GiB)
[ Info: InterpolatingAdjoint
  5.475 s (18755710 allocations: 2.54 GiB)
[ Info: QuadratureAdjoint
  1.659 s (5810587 allocations: 807.58 MiB)

```

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 27, 2020, 9:03pm UTC](https://discourse.julialang.org/t/why-the-separation-of-odeproblem-and-solve-in-differentialequations-jl/43737/12 "2020-07-27T21:03:52Z")

</div>

For that problem I would suggest use `QuadratureAdjoint(autojacvec=ReverseDiffVJP(true))`

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 28, 2020, 12:54am UTC](https://discourse.julialang.org/t/why-the-separation-of-odeproblem-and-solve-in-differentialequations-jl/43737/13 "2020-07-28T00:54:30Z")

</div>

For reference, before:

```julia
15.938 ms (50871 allocations: 48.77 MiB)

```

after:

```julia
6.676 ms (65565 allocations: 1.66 MiB)

```

with discussion here: [https://github.com/SciML/DiffEqSensitivity.jl/issues/316](https://github.com/SciML/DiffEqSensitivity.jl/issues/316)

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [July 28, 2020, 1:10am UTC](https://discourse.julialang.org/t/why-the-separation-of-odeproblem-and-solve-in-differentialequations-jl/43737/14 "2020-07-28T01:10:54Z")

</div>

We can probably make this faster still with a MTK generated VJP.