Weak adaptive solvers give rise to Out-of-memory errors

[WouterJRV]

See https://discourse.julialang.org/t/why-does-my-devectorized-code-keep-allocating-memory/62727/11

For a set of Order 1e4 coupled SDE equations with additive diagonal noise, SOSRA and similar works fine, but the DRI-methods crash with OutOfMemoryError()

Apparently, this is "because Ihat2 uses length(u)^2 memory for the pairwise 3-point random numbers"

[ChrisRackauckas]

@frankschae is there a good way to avoid allocating the matrix of all pairwise random numbers there?

[frankschae]

yeah we never use the full matrix Ihat2 but only its entries in perform_step!. These (scalar) noise increments are only used in the last stage value in all SRK schemes IIRC. So, one option is probably to compute them on the fly when they are needed. I tried to understand if we can avoid computing them for diagonal processes. Perhaps, only the diagonal of the matrix is sufficient in that case anyway.

[WouterJRV]

With StochasticDiffEq v6.36.0 it should be working, right? Unfortunately, I'm still getting the same error.

[ChrisRackauckas]

It can't be the same error because the line it was pointing to doesn't exist anymore. Maybe similar error different stack trace?

[WouterJRV]

ERROR: LoadError: OutOfMemoryError() Stacktrace: [1] Array at ./boot.jl:408 [inlined] [2] Array at ./boot.jl:416 [inlined] [3] zeros at ./array.jl:525 [inlined] [4] zeros at ./array.jl:521 [inlined] [5] alg_cache(::DRI1, ::SDEProblem{Array{Float64,3},Tuple{Float64,Float64},true,Tuple{Int64,Int64,Float64,Array{Float64,2},Float64,Array{Float64,2},Array{Float64,2},Int64,Array{Float64,2}},Nothing,SDEFunction{true,typeof(f!),typeof(g!),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(g!),Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},Nothing}, ::Array{Float64,3}, ::Array{Float64,3}, ::Array{Float64,3}, ::Tuple{Int64,Int64,Float64,Array{Float64,2},Float64,Array{Float64,2},Array{Float64,2},Int64,Array{Float64,2}}, ::Array{Float64,3}, ::Array{Float64,3}, ::Nothing, ::Type{T} where T, ::Type{T} where T, ::Type{T} where T, ::Array{Float64,3}, ::SDEFunction{true,typeof(f!),typeof(g!),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing}, ::Float64, ::Float64, ::Type{Val{true}}) at /home/wouter/.julia/packages/StochasticDiffEq/MzqBB/src/caches/srk_weak_caches.jl:598 [6] __init(::SDEProblem{Array{Float64,3},Tuple{Float64,Float64},true,Tuple{Int64,Int64,Float64,Array{Float64,2},Float64,Array{Float64,2},Array{Float64,2},Int64,Array{Float64,2}},Nothing,SDEFunction{true,typeof(f!),typeof(g!),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(g!),Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},Nothing}, ::DRI1, ::Array{Any,1}, ::Array{Any,1}, ::Type{T} where T, ::Type{Val{true}}; saveat::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_noise::Bool, save_on::Bool, save_start::Bool, save_end::Bool, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmax::Rational{Int64}, qmin::Rational{Int64}, qoldinit::Rational{Int64}, fullnormalize::Bool, failfactor::Int64, beta2::Rational{Int64}, beta1::Rational{Int64}, delta::Rational{Int64}, maxiters::Int64, dtmax::Float64, dtmin::Float64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, force_dtmin::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, initialize_integrator::Bool, seed::UInt64, alias_u0::Bool, alias_jumps::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/wouter/.julia/packages/StochasticDiffEq/MzqBB/src/solve.jl:453 [7] #__solve#99 at /home/wouter/.julia/packages/StochasticDiffEq/MzqBB/src/solve.jl:6 [inlined] [8] #solve_call#56 at /home/wouter/.julia/packages/DiffEqBase/FMY9y/src/solve.jl:61 [inlined] [9] solve_up(::SDEProblem{Array{Float64,3},Tuple{Float64,Float64},true,Tuple{Int64,Int64,Float64,Array{Float64,2},Float64,Array{Float64,2},Array{Float64,2},Int64,Array{Float64,2}},Nothing,SDEFunction{true,typeof(f!),typeof(g!),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(g!),Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},Nothing}, ::Nothing, ::Array{Float64,3}, ::Tuple{Int64,Int64,Float64,Array{Float64,2},Float64,Array{Float64,2},Array{Float64,2},Int64,Array{Float64,2}}, ::DRI1; kwargs::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}}}}}) at /home/wouter/.julia/packages/DiffEqBase/FMY9y/src/solve.jl:82 [10] #solve#57 at.../.julia/packages/DiffEqBase/FMY9y/src/solve.jl:70 [inlined] [11] top-level scope at ./timing.jl:174 in ...

[frankschae]

hmm that one looks worse. It's pointing to:

g2 = [zero(noise_rate_prototype) for k=1:m]

StochasticDiffEq.jl/src/caches/srk_weak_caches.jl

Line 598 in a4497c5

g2 = [zero(noise_rate_prototype) for k=1:m]

I don't see how one could change that to avoid the m stored copies except by tons of recomputations (because g2 and g3 are used at several stages). Do you have any simplifications of your model with respect to g? In your first code on discourse, you were using DRI1NM which doesn't have this issue because it assumes that the diagonal components do not mix.

[WouterJRV]

Yes, tried it earlier and it gives a similar error. didn't read throuth the stack trace because I assumed it was the same, but apparently that's not the case

ERROR: LoadError: OutOfMemoryError() Stacktrace: [1] Array at ./boot.jl:408 [inlined] [2] Array at ./boot.jl:416 [inlined] [3] zeros at ./array.jl:525 [inlined] [4] zeros at ./array.jl:521 [inlined] [5] alg_cache(::DRI1NM, ::SDEProblem{Array{Float64,3},Tuple{Float64,Float64},true,Tuple{Int64,Int64,Float64,Array{Float64,2},Float64,Array{Float64,2},Array{Float64,2},Int64,Array{Float64,2}},Nothing,SDEFunction{true,typeof(f!),typeof(g!),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(g!),Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},Nothing}, ::Array{Float64,3}, ::Array{Float64,3}, ::Nothing, ::Tuple{Int64,Int64,Float64,Array{Float64,2},Float64,Array{Float64,2},Array{Float64,2},Int64,Array{Float64,2}}, ::Array{Float64,3}, ::Array{Float64,3}, ::Nothing, ::Type{T} where T, ::Type{T} where T, ::Type{T} where T, ::Array{Float64,3}, ::SDEFunction{true,typeof(f!),typeof(g!),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing}, ::Float64, ::Float64, ::Type{Val{true}}) at /home/wouter/.julia/packages/StochasticDiffEq/MzqBB/src/caches/srk_weak_caches.jl:642 [6] __init(::SDEProblem{Array{Float64,3},Tuple{Float64,Float64},true,Tuple{Int64,Int64,Float64,Array{Float64,2},Float64,Array{Float64,2},Array{Float64,2},Int64,Array{Float64,2}},Nothing,SDEFunction{true,typeof(f!),typeof(g!),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(g!),Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},Nothing}, ::DRI1NM, ::Array{Any,1}, ::Array{Any,1}, ::Type{T} where T, ::Type{Val{true}}; saveat::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}, tstops::Tuple{}, d_discontinuities::Tuple{}, save_idxs::Nothing, save_everystep::Bool, save_noise::Bool, save_on::Bool, save_start::Bool, save_end::Bool, callback::Nothing, dense::Bool, calck::Bool, dt::Float64, adaptive::Bool, gamma::Rational{Int64}, abstol::Nothing, reltol::Nothing, qmax::Rational{Int64}, qmin::Rational{Int64}, qoldinit::Rational{Int64}, fullnormalize::Bool, failfactor::Int64, beta2::Rational{Int64}, beta1::Rational{Int64}, delta::Rational{Int64}, maxiters::Int64, dtmax::Float64, dtmin::Float64, internalnorm::typeof(DiffEqBase.ODE_DEFAULT_NORM), isoutofdomain::typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), unstable_check::typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), verbose::Bool, force_dtmin::Bool, timeseries_errors::Bool, dense_errors::Bool, advance_to_tstop::Bool, stop_at_next_tstop::Bool, initialize_save::Bool, progress::Bool, progress_steps::Int64, progress_name::String, progress_message::typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), userdata::Nothing, initialize_integrator::Bool, seed::UInt64, alias_u0::Bool, alias_jumps::Bool, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/wouter/.julia/packages/StochasticDiffEq/MzqBB/src/solve.jl:453 [7] #__solve#99 at /home/wouter/.julia/packages/StochasticDiffEq/MzqBB/src/solve.jl:6 [inlined] [8] #solve_call#56 at /home/wouter/.julia/packages/DiffEqBase/FMY9y/src/solve.jl:61 [inlined] [9] solve_up(::SDEProblem{Array{Float64,3},Tuple{Float64,Float64},true,Tuple{Int64,Int64,Float64,Array{Float64,2},Float64,Array{Float64,2},Array{Float64,2},Int64,Array{Float64,2}},Nothing,SDEFunction{true,typeof(f!),typeof(g!),UniformScaling{Bool},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},typeof(g!),Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},Nothing}, ::Nothing, ::Array{Float64,3}, ::Tuple{Int64,Int64,Float64,Array{Float64,2},Float64,Array{Float64,2},Array{Float64,2},Int64,Array{Float64,2}}, ::DRI1NM; kwargs::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}}}}}) at /home/wouter/.julia/packages/DiffEqBase/FMY9y/src/solve.jl:82 [10] #solve#57 at /home/wouter/.julia/packages/DiffEqBase/FMY9y/src/solve.jl:70 [inlined] [11] top-level scope at ./timing.jl:174

[frankschae]

Interesting.. and SOSRI() and SOSRA() do work?
The error is pointing to

StochasticDiffEq.jl/src/caches/srk_weak_caches.jl

Line 642 in a4497c5

g3 = zero(noise_rate_prototype)

which is just the third cached noise vector. SOSRI() should have (up to the Butcher table) exactly the same cache consumption:

StochasticDiffEq.jl/src/caches/rossler_caches.jl

Line 453 in a4497c5

g3 = zero(noise_rate_prototype); g4 = zero(noise_rate_prototype)

and SOSRA() as well (-one vector of size m for the random numbers):

StochasticDiffEq.jl/src/caches/sra_caches.jl

Line 260 in a4497c5

g3 = zero(noise_rate_prototype)

[WouterJRV]

Yes, with SOSRA I got
508.131326 seconds (18.10 M allocations: 825.235 MiB, 0.09% gc time)
and with SOSRI
762.754098 seconds (23.18 M allocations: 1.012 GiB, 0.07% gc time)