The notion of "initialization" relates to the prior part of the model.

TL;DR: In the vast majority of cases, just stick to the exact Taylor-mode initialization TaylorModeInit!

Background: The prior

We model the ODE solution $y(t)$ with a Gauss–Markov prior. More precisely, let

\[\begin{aligned} Y(t) = \left[ Y^{(0)}(t), Y^{(1)}(t), \dots Y^{(q)}(t) \right], \end{aligned}\]

be the solution to the SDE

\[\begin{aligned} \text{d} Y^{(i)}(t) &= Y^{(i+1)}(t) \ \text{d}t, \qquad i = 0, \dots, q-1 \\ \text{d} Y^{(q)}(t) &= \textcolor{#389826}{A} Y(t) \ \text{d}t + \textcolor{#4063D8}{\Gamma} \ \text{d}W(t), \\ Y(0) &\sim \textcolor{#9558B2}{ \mathcal{N} \left( \mu_0, \Sigma_0 \right) }. \end{aligned}\]

Then $Y^{(i)}(t)$ models the $i$-th derivative of $y(t)$. In this section, we consider the initial distribution $\textcolor{purple}{ \mathcal{N} \left( \mu_0, \Sigma_0 \right) }$. If you're more interested in the drift matrix $\textcolor{#389826}{A}$ check out the Priors section, and for more info on the diffusion $\textcolor{#4063D8}{\Gamma}$ check out the Diffusion models and calibration section.

Setting the initial distribution

Let's assume an initial value problem of the form

\[\begin{aligned} \dot{y}(t) &= f(y(t), t), \qquad [0, T], \\ y(0) &= y_0. \end{aligned}\]

It is clear that this contains quite some information for $Y(0)$: The initial value $y_0$ and the vector field $f$ imply

\[\begin{aligned} Y^{(0)}(0) &= y_0, \\ Y^{(1)}(0) &= f(y_0, 0). \end{aligned}\]

It turns out that we can also compute higher-order derivatives of $y$ with the chain rule, and then use these to better initialize $Y^{(i)}(0)$. This, done efficiently with Taylor-mode automatic differentiation by using TaylorIntegration.jl, is what TaylorModeInit does. See also [1].

In the vast majority of cases, just stick to the exact Taylor-mode initialization TaylorModeInit!



Exact initialization via Taylor-mode automatic differentiation up to order order.

This is the recommended initialization method!

It uses TaylorIntegration.jl to efficiently compute the higher-order derivatives of the solution at the initial value, via Taylor-mode automatic differentiation.

In some special cases it can happen that TaylorIntegration.jl is incompatible with the given problem (typically because the problem definition does not allow for elements of type Taylor). If this happens, try one of SimpleInit, ForwardDiffInit (for low enough orders), ClassicSolverInit.


  • [7] Krämer et al, "Stable Implementation of Probabilistic ODE Solvers" (2020)

Simple initialization, only with the given initial value and derivative.

The remaining derivatives are set to zero with unit covariance (unless specified otherwise by setting a custom FixedDiffusion).

ClassicSolverInit(; alg=OrdinaryDiffEq.Tsit5(), init_on_ddu=false)

Initialization via regression on a few steps of a classic ODE solver.

In a nutshell, instead of specifying $\mu_0$ exactly and setting $\Sigma_0=0$ (which is what TaylorModeInit does), use a classic ODE solver to compute a few steps of the solution, and then regress on the computed values (by running a smoother) to compute $\mu_0$ and $\Sigma_0$ as the mean and covariance of the smoothing posterior at time 0. See also [2].

The initial value and derivative are set directly from the given initial value problem; optionally the second derivative can also be set via automatic differentiation by setting init_on_ddu=true.


  • alg: The solver to be used. Can be any solver from OrdinaryDiffEq.jl.
  • init_on_ddu: If true, the second derivative is also initialized exactly via automatic differentiation with ForwardDiff.jl.


  • [7] Krämer et al, "Stable Implementation of Probabilistic ODE Solvers" (2020)
  • [8] Schober et al, "A probabilistic model for the numerical solution of initial value problems", Statistics and Computing (2019)