GitHub

Advice for Performant SDE Simulations

While there exist relatively straightforward approaches to manage performance for ODE and jump simulations, this is generally not the case for SDE simulations. Below, we briefly describe some options. However, as one starts to investigate these, one quickly reaches what is (or could be) active areas of research.

SDE solver selection

We have previously described how ODE solver selection can impact simulation performance. Again, it can be worthwhile to investigate solver selection's impact on performance for SDE simulations. Throughout this documentation, we generally use the STrapezoid solver as the default choice. However, if the DifferentialEquations package is loaded

using DifferentialEquations

automatic SDE solver selection enabled (just like is the case for ODEs by default). Generally, the automatic SDE solver choice enabled by DifferentialEquations is better than just using STrapezoid. Next, if performance is critical, it can be worthwhile to check the list of available SDE solvers to find one with advantageous performance for a given problem. When doing so, it is important to pick a solver compatible with non-diagonal noise and with Ito problems.

Options for Jacobian computation

In the section on ODE simulation performance, we describe various options for computing the system Jacobian, and how these could be used to improve performance for implicit solvers. These can be used in tandem with implicit SDE solvers (such as STrapezoid). However, due to additional considerations during SDE simulations, it is much less certain whether these will actually have any impact on performance. So while these options might be worth reading about and trialling, there is no guarantee that they will be beneficial.

Parallelisation on CPUs and GPUs

We have previously described how simulation parallelisation can be used to improve performance when multiple ODE simulations are carried out. The same approaches can be used for SDE simulations. Indeed, it is often more relevant for SDEs, as these are often re-simulated using identical simulation conditions (to investigate their typical behaviour across many samples). CPU parallelisation of SDE simulations uses the same approach as ODEs. GPU parallelisation requires some additional considerations, which are described below.

GPU parallelisation of SDE simulations

GPU parallelisation of SDE simulations uses a similar approach as that for ODE simulations. The main differences are that SDE parallelisation requires a GPU SDE solver (like GPUEM) and fixed time stepping.

We will assume that we are using the CUDA GPU hardware, so we will first load the CUDA.jl backend package, as well as DiffEqGPU:

using CUDA, DiffEqGPU

Which backend package you should use depends on your available hardware, with the alternatives being listed here.

Next, we create the SDEProblem which we wish to simulate. Like for ODEs, we ensure that all vectors are static vectors and that all values are Float32s. Here we prepare the parallel simulations of a simple birth-death process.

using Catalyst, StochasticDiffEq, StaticArrays
bd_model = @reaction_network begin
    (p,d), 0 <--> X
end
@unpack X, p, d = bd_model

u0 = @SVector [X => 20.0f0]
tspan = (0.0f0, 10.0f0)
ps = @SVector [p => 10.0f0, d => 1.0f0]
sprob = SDEProblem(bd_model, u0, tspan, ps)

The SDEProblem is then used to create an EnsembleProblem.

eprob = EnsembleProblem(sprob)

Finally, we can solve our EnsembleProblem while:

  • Using a valid GPU SDE solver (either GPUEM or GPUSIEA).
  • Designating the GPU ensemble method, EnsembleGPUKernel (with the correct GPU backend as input).
  • Designating the number of trajectories we wish to simulate.
  • Designating a fixed time step size.
esol = solve(eprob, GPUEM(), EnsembleGPUKernel(CUDA.CUDABackend()); trajectories = 10000, dt = 0.01)

Above we parallelise GPU simulations with identical initial conditions and parameter values. However, varying these is also possible.

Multilevel Monte Carlo

An approach for speeding up parallel stochastic simulations is so-called multilevel Monte Carlo approaches (MLMC). These are used when a stochastic process is simulated repeatedly using identical simulation conditions. Here, instead of performing all simulations using identical tolerance, the ensemble is simulated using a range of tolerances (primarily lower ones, which yields faster simulations). Currently, StochasticDiffEq.jl do not have a native implementation for performing MLMC simulations (this will hopefully be added in the future). However, if high performance of parallel SDE simulations is required, these approaches may be worth investigating.