SDE Solvers (seemingly without reason) hitting maxiters

[TorkelE]

Originally posted as https://discourse.julialang.org/t/strange-maxiters-numeric-instability-occured-when-solving-certain-sde/49392 however, after looking closer at it I think it might be a bug in the package, so it probably fits better here.

Short Summary

• I have an SDEProblem which when I try to solve it sometimes, strangely, hits maxiters.
• The cause seems to be dt suddenly becoming very small.
• Solving for the same parameter set, the error only occurs sometimes. By fixing the seed it becomes consistent.
• Some parameter values seem likelier to produce the error than others.
• The problem does not seem to relate to the amount of noise in the system.
• The problem only seems to occur for Ito type methods, but not for Stratonovich type methods.

Full Story

Basically, I am investigating a simple SDE system, performing a large number of simulation over a large number of parameters. However, I have recently started getting strange

┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase /home/torkelloman/.julia/packages/DiffEqBase/V7P18/src/integrator_interface.jl:329

This is rather rare (been when solving 10s of thousands of times it comes up quite a lot). It does not come up every time I solve a certain SDEProblem (but I can make it reproducible by specifying a seed). There seem to be some patterns though, where for some parameter selections the solver is much more likely to hit the maxiters.

This is a basic reproducible example:

using StochasticDiffEq

function sys_f(du,u,p,t)
    σ,A = u
    S,D,τ,v0,n,η = p
    du[1] = v0+(S*σ)^n/((S*σ)^n+(D*A)^n+1) - σ
    du[2] = τ*(σ-A)
end
function sys_g(du,u,p,t)
    σ,A = u
    S,D,τ,v0,n,η = p
    du[1,1] = η*sqrt(max(v0+(S*σ)^n/((S*σ)^n+(D*A)^n+1),0.))
    du[1,2] = -η*sqrt(max(σ,0.))
    du[1,3] = 0
    du[1,4] = 0
    du[2,1] = 0
    du[2,2] = 0
    du[2,3] = η*sqrt(max(τ*σ,0.))
    du[2,4] = -η*sqrt(max(τ*A,0.))
end


u0 = [0.025, 0.025]
tspan = (0.,200.)
p = [28.48035868435802, 0.2527461586781728, 0.05, 0.025, 3.0, 0.1]

sde_prob = SDEProblem(sys_f,sys_g,u0,tspan,p,noise_rate_prototype=zeros(2,4));

sol = solve(sde_prob,ImplicitEM(),seed=4963697075065472807);

Solving without this seed, things are mostly fine:

for repeat = 1:1000
    sol = solve(sde_prob,ImplicitEM())
end

yields maybe only 1 such error per run.

Even increasing maxiters a bit, I receive no benefit:

sol = solve(sde_prob,ImplicitEM(),seed=4963697075065472807,maxiters=10000000);

still errors.

To try ad figure out what is happening I plot the solution, and the time step length, for an attempt where I ran into the problem (and one where I didn't).

Here's the solutions using something like

sol = solve(sde_prob,ImplicitEM(),seed=4963697075065472807);
dts = diff(sol.t)
plot(sol.t[2:end],dts,framestyle=:box)

I can get the dt:

In these, I only plot the successful solution to the same length as the time point where the failed one failed. Here's a comparison where I plot the full length of the successful simulation.

Looking at the (failed) solution directly, the last 20 values of sol.t becomes

 3.0696546597633163
 3.0696553650908487
 3.0696561585843227
 3.0696568640351365
 3.069657657667302
 3.069658363134364
 3.069659156784809
 3.0696598622394653
 3.069660655875954
 3.0696613614302923
 3.069662155178923
 3.0696628609243257
 3.0696636548879033
 3.0696643602892197
 3.0696651538657007
 3.0696658593007053
 3.0696666529150853
 3.0696673582528358
 3.069668151757805
 3.069668856962127
 3.0696696503169894

while the last 20 values of dts are

 7.934824788335959e-7
 7.053275323798402e-7
 7.934934740383426e-7
 7.054508137649407e-7
 7.936321653190248e-7
 7.054670621009507e-7
 7.936504449190807e-7
 7.054546564688735e-7
 7.936364885274827e-7
 7.055543385092733e-7
 7.937486308229325e-7
 7.057454025627408e-7
 7.939635775500165e-7
 7.054013164697892e-7
 7.935764809730017e-7
 7.054350046331592e-7
 7.936143799902595e-7
 7.053377504284697e-7
 7.935049692875396e-7
 7.052043220490134e-7
 7.933548622496289e-7

The problem does not seem to depend on the amount of noise in the system. I can check the outputs of sys_f and sys_g at the final value (where it hits maxiters):

u_in = zeros(2)
sys_f(u_in,sol.u[end],p,0.)
u_in

gives

2-element Array{Float64,1}:
 0.5475974457463257
 0.023042909797453985

and

u_in = zeros(2,4)
sys_g(u_in,sol.u[end],p,0.)
u_in

gives

2×4 Array{Float64,2}:
 0.101223  -0.0690655  0.0         0.0
 0.0        0.0        0.0154435  -0.00284127

Further proving this point, I look at another parameter set. Here we note that the last parameter scales the noise, by making this very small, the noise become negligible. This error happens even when this parameter is 0.001 (very small).

using StochasticDiffEq

function sys_f(du,u,p,t)
    σ,A = u
    S,D,τ,v0,n,η = p
    du[1] = v0+(S*σ)^n/((S*σ)^n+(D*A)^n+1) - σ
    du[2] = τ*(σ-A)
end
function sys_g(du,u,p,t)
    σ,A = u
    S,D,τ,v0,n,η = p
    du[1,1] = η*sqrt(max(v0+(S*σ)^n/((S*σ)^n+(D*A)^n+1),0.))
    du[1,2] = -η*sqrt(max(σ,0.))
    du[1,3] = 0
    du[1,4] = 0
    du[2,1] = 0
    du[2,2] = 0
    du[2,3] = η*sqrt(max(τ*σ,0.))
    du[2,4] = -η*sqrt(max(τ*A,0.))
end

u0 = [0.025, 0.025]
tspan = (0.,200.)
p = [6.7341506577508214, 1.2224984640507832, 0.05, 0.025, 3.0, 0.001]

sde_prob = SDEProblem(sys_f,sys_g,u0,tspan,p,noise_rate_prototype=zeros(2,4));

sol1 = solve(sde_prob,ImplicitEM(),seed=5105556673070507935); sol = sol1;
sol2 = solve(sde_prob,ImplicitEM());

┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase /home/torkelloman/.julia/packages/DiffEqBase/V7P18/src/integrator_interface.jl:329

using Plots
gr(); default(fmt = :png);
p1 = plot(sol1,framestyle=:box,title="Bad stuff",xlimit=(0.,sol1.t[end]))
p2 = plot(sol2,framestyle=:box,title="Everything is fine",xlimit=(0.,sol1.t[end]))
p3 = plot(sol2,framestyle=:box,title="Everything is fine",xlimit=(0.,sol2.t[end]))
plot_trajectories = plot(p1,p2,p3,size=(1200,400),layout=(1,3))

p1 = plot(sol1.t[2:end],diff(sol1.t),framestyle=:box,xlimit=(0.,sol1.t[end]),title="Bad stuff")
p2 = plot(sol2.t[2:end],diff(sol2.t),framestyle=:box,xlimit=(0.,sol1.t[end]),title="Everything is fine")
p3 = plot(sol2.t[2:end],diff(sol2.t),framestyle=:box,xlimit=(0.,sol2.t[end]),title="Everything is fine")
plot_dts = plot(p1,p2,p3,size=(1200,400),layout=(1,3))

Finally, I have tried investigating this over several solvers. Looking at the initial problem I try:

sol_ImplicitEM = solve(sde_prob,ImplicitEM(),seed=4963697075065472807);
sol_STrapezoid = solve(sde_prob,STrapezoid(),seed=4963697075065472807);
sol_SImplicitMidpoint = solve(sde_prob,SImplicitMidpoint(),seed=4963697075065472807);
sol_ImplicitEulerHeun = solve(sde_prob,ImplicitEulerHeun(),seed=4963697075065472807);
sol_ISSEM = solve(sde_prob,ISSEM(),seed=4963697075065472807);
sol_ISSEulerHeun = solve(sde_prob,ISSEulerHeun(),seed=4963697075065472807);

Of these, sol_ImplicitEM, sol_STrapezoid , sol_SImplicitMidpoint all fails, while ImplicitEulerHeun, ISSEM, ISSEulerHeun succeeds. See plots (know it is not super helpful, but best I got to show how they fail. Solutions: dts:

Finally, I try scanning for other seeds which causes failure in the remaining three. Did not find any for ImplicitEulerHeun or ISSEulerHeun, however, ISSEM fails for seed=15088718199444158177

sol_ISSEM = solve(sde_prob,ISSEM(),seed=15088718199444158177);

Looking at this seed specifically, ImplicitEM and SImplicitMidpoint also fails, while the remaining three succeeds.

Trying to find a pattern, it seems like I can find errors for the Ito methods, but not for the Stratonovich methods.

This is the output of Pkg.status()

  [717857b8] DSP v0.6.8
  [0c46a032] DifferentialEquations v6.15.0
  [5ab0869b] KernelDensity v0.6.1
  [91a5bcdd] Plots v1.6.8
  [f27b6e38] Polynomials v1.1.10
  [10745b16] Statistics

[ChrisRackauckas]

I thought I mentioned that it's not odd or without a reason, just random. I show in https://arxiv.org/abs/1804.04344 that the stability of a method can be heavily dependent on the noise terms that you get. This is showing all of the signs of instability though, hence it's definitely what I termed pathwise stiffness and the standard drift-implicit methods do nothing here (since it's "stiffness in the noise"). This will be a good test for some of the new methods we're cooking up though, like a TruncatedFullyImplicitEM and the like, but that will take some time to develop.

Note that looking at a single seed doesn't make much sense, because with adaptivity two methods on the same seed will still get a different Brownian path. So it's really just, for a given method on a given seed, the Brownian path that it sees is an unstable path.

[ChrisRackauckas]

Though @YingboMa if you could look at the nonlinear solvers in this case to have fresh eyes and maybe see if I missed anything there, that would be great.

[TorkelE]

Seems that (yet again) when I start to believe that I'm getting a vague grip on how to solve things numerically, something will pop up to rudely prove me the difference!

[ChrisRackauckas]

SDEs are hard. Stiffness is SDEs is even weirder than ODEs 🤷

[TorkelE]

On a general level, it is kinda cool that there's this kind of strange behaviours in numerical systems, makes the field seem quite cool. On a local level, it kind of sucks.

Practically, would it be appropriate to use Stratonovich type methods to solve this kind of chemical Langevin problems?

[TorkelE]

Actually, I think I have found similar errors, but for other parameter sets, for the Stratonovich methods as well.