Help me beat Lsoda


I want to simulate a Jump process for biological applications quite precisely and fast. I developped a package PiecewiseDeterministicMarkovProcesses.jl based on DifferentialEquations.jl to simulate these jump processes.

This package has been optimized a lot to reduce allocations and be complient with AD. Note that for my application where the dimension is > 35 elements, I dont want to use StaticArrays.

In my application, I find lsoda to be the fastest and the most reliable despite the fact that I cannot use an ODE iterator with it (and hence it allocates). Thus, I want to be able to be as close as possible performance wise, to the lsoda solver but I dont manage to tweak the options to achieve this. If anybody has experience with this, I would not mind a little push.

Thank you.

PS: the MWE below gives the following timings. Note that we have an analytical expression of the jump process so we can attest the error. The results are

Solvers comparison

--> norm difference = 0.007230211093428807  - solver = lsoda
  0.000935 seconds (9.61 k allocations: 502.719 KiB)

--> norm difference = 0.00043447196503620944  - solver = rodas5P
  0.001447 seconds (2.08 k allocations: 160.797 KiB)

--> norm difference = 12.221574659298312  - solver = TRBDF2
  0.002779 seconds (2.07 k allocations: 158.578 KiB)

--> norm difference = 0.010638187872245908  - solver = CVODEAdams
  0.001777 seconds (17.38 k allocations: 480.062 KiB)

--> norm difference = 0.005070017916750658  - solver = CVODEBDF
  0.002406 seconds (18.05 k allocations: 521.891 KiB)

--> norm difference = 0.0001112569666474883  - solver = rodas4P
  0.002659 seconds (4.47 k allocations: 347.438 KiB)

--> norm difference = 0.005069408256076713  - solver = cvode
  0.003115 seconds (24.58 k allocations: 820.281 KiB)

--> norm difference = 0.13379368165760752  - solver = Rosenbrock23
  0.005789 seconds (14.07 k allocations: 1.071 MiB)

--> norm difference = 0.01244811034121085  - solver = AutoTsit5-RS23
  0.000439 seconds (329 allocations: 24.422 KiB)
Code on PiecewiseDeterministicMarkovProcesses#master
# using Revise
using PiecewiseDeterministicMarkovProcesses, LinearAlgebra, Random, DifferentialEquations, Sundials
	const PDMP = PiecewiseDeterministicMarkovProcesses

function AnalyticalSampleCHV(xc0, xd0, ti, nj::Int64)
	xch = [xc0[1]]
	xdh = [xd0[1]]
	th  = [ti]
	t = ti
	while length(th)<nj
		xc = xch[end]
		xd = xdh[end]
		S = -log(rand())
		if mod(xd,2) == 0
			t += 1/10*log(1+10*S/xc)
			push!(xch,xc + 10 * S )
			t += 1/(3xc)*(exp(3S)-1)
			push!(xch,xc * exp(-3S) )
		push!(xdh,xd + 1 )

		S = -log(rand())
	return th, xch, xdh

function F!(ẋ, xc, xd, parms, t)
	if mod(xd[1], 2)==0
		ẋ[1] = 10xc[1]
		ẋ[1] = -3xc[1]^2

R(x) = x

function R!(rate, xc, xd, parms, t, issum::Bool)
	# rate function
	if issum == false
		rate[1] = R(xc[1])
		rate[2] = parms[1]
		return 0., parms[1] + 50.
		return R(xc[1]) + parms[1], parms[1] + 50.

xc0 = [1.0]
	xd0 = [0, 0]

	nu = [1 0;0 -1]
	parms = [.0]
	ti = 0.322156
	tf = 100000.
	nj = 50

errors = Float64[]

	res_a_chv = AnalyticalSampleCHV(xc0,xd0,ti,nj)

problem = PDMP.PDMPProblem(F!, R!, nu, xc0, xd0, parms, (ti, tf))
println("\n\nSolvers comparison")
	for ode in [	(:lsoda,"lsoda"),
		res = PDMP.solve(problem, CHV(ode[1]); n_jumps = nj, abstol = 1e-9, reltol = 1e-7)
		printstyled(color=:green, "\n--> norm difference = ", norm(res.time - res_a_chv[1], Inf64), "  - solver = ",ode[2],"\n")
		res = @time PDMP.solve(problem, CHV(ode[1]); n_jumps = nj, abstol = 1e-9, reltol = 1e-7)
		push!(errors,norm(res.time - res_a_chv[1], Inf64))

Have you tried FBDF (QNDF would also be good to try)? In my experience, that is often very good for stiff systems.

1 Like

None of them return the result, it errors because of dtmin.

1 Like

That’s slightly odd. Another to try is rodas5P which isn’t documented yet, but should be better than Rodas4P.

1 Like

Yes it comes close (see above) and it much more precise. I tried to change the tolerance to decrease the computation time of Rodas5 while having similar error to LSODA. I arrive at

--> norm difference = 0.008140039428781165  - solver = rodas5P
  0.000922 seconds (1.40 k allocations: 107.453 KiB)

so roughly similar to LSODA but with much smaller allocations.

However in my application, it is 10x slower maybe because of jacobian computation


You should read through Solving Large Stiff Equations · DifferentialEquations.jl and in particular, try jacobian free methods (Newton krylov).

1 Like

Are you sure it’s stiff?

1 Like

We solve either

dx/dt = G(x) = 10x / x (= 10)


dx/dt = G(x) = -3x^2 / x (= x)

which is indeed not stiff but the last simplification is not done in the code and x can be small to zero.

Thank you for your remark, my example is maybe to simple for a MWE

LSODA switches between a stiff and non-stiff ODE solver. If the problem is non-stiff, then LSODA will switch over to the non-stiff method which will have a much lower overhead than the stuff for stiff equations, and it’ll beat every single one handedly. But of course, that’s just because stiff ODE solvers are not optimized for non-stiff ODEs!

I just had a collaborator make this mistake the other day. There’s a few benchmarks which directly show this behavior:


The point is, if LSODA is lower than everything you’re trying, it may just be in the other class of methods. Or, auto-switch may just be better.

Note you probably want to try AutoVern8(Rodas5P()) for this case. Tsit5 + Rosenbrock23 is only for very high tolerances.


Thanks a lot for your input! In my use case, I was able to beat lsoda with CVODE_BDF(linear_solver=:GMRES). I tried AutoVern8(Rodas5P()) but it could not compute the first jump time.

┌ Warning: Instability detected. Aborting
└ @ SciMLBase ~/.julia/packages/SciMLBase/xWByK/src/integrator_interface.jl:523
ERROR: AssertionError: Could not compute next jump time 1.
Return code = Unstable
 0.0 < 0.0,

“0.0 < 0.0” est un assert de PDMP

1 Like

Try tweaking the switching parameters? If we can get this case for benchmarks that would be helpful

1 Like

I tried with AutoVern8(Rodas5P(), nonstifftol = 1//100) but it did not succeed

Can we turn this into a SciMLBenchmarks work-precision diagram? That will make it more concrete for performance tracking.

1 Like

Will do once the paper is out.

Just to mention, there is a part similar to this in the code,

dxdt = F(x, y)
dydt = G(x, y)
dzdt = [x, y] ∈ fixed_region

Hence, the vector field is only piecewise continuous although it should not have a big effect because z does not feedback on x,y

The paper is accepted! Hence, I can put the code online here. Any improvement would be a major win for us

1 Like

Can we get the ODEProblem definition out?

Tough, it is a PDMP sampled with the CHV method. Maybe it is a benchamark for JumpProcesses?

It would be good to test the new Coevolve against that. I don’t know if you’ve seen that paper but it tests against the CHV in a few problems. Can you share the MWE and we can translate it?

I don’t know if you’ve seen that paper

Yes I have seen this paper but can’t find it again. I had a few comments about it but I forgot

Can you share the MWE and we can translate it?

It is written for PiecewiseDeterministicMarkovProcesses.jl, not for JumpProcesses. The drift is here, the total rate is here and the solve.

It would benefit from PeriodicCallback but this is not in PiecewiseDeterministicMarkovProcesses.jl

I’m one of the authors of the Coevolve paper, Let me see if I can map this problem to Coevolve tomorrow and see how it does.

1 Like