Probabilistic Exponential Integrators
Authors: Nathanael Bosch
Philipp Hennig
Filip Tronarp
Venue: NeurIPS 2023 poster (2023)
Keywords: Probabilistic numerics
differential equations
exponential integrators
Kalman filters
Gaussian processes
TL;DR: We propose a new class of stable probabilistic numerical solvers for stiff ordinary differential equations.

Abstract: Probabilistic solvers provide a flexible and efficient framework for simulation, uncertainty quantification, and inference in dynamical systems. However, like standard solvers, they suffer performance penalties for certain stiff systems, where small steps are required not for reasons of numerical accuracy but for the sake of stability. This issue is greatly alleviated in semi-linear problems by the probabilistic exponential integrators developed in this paper. By including the fast, linear dynamics in the prior, we arrive at a class of probabilistic integrators with favorable properties. Namely, they are proven to be L-stable, and in a certain case reduce to a classic exponential integrator---with the added benefit of providing a probabilistic account of the numerical error. The method is also generalized to arbitrary non-linear systems by imposing piece-wise semi-linearity on the prior via Jacobians of the vector field at the previous estimates, resulting in probabilistic exponential Rosenbrock methods. We evaluate the proposed methods on multiple stiff differential equations and demonstrate their improved stability and efficiency over established probabilistic solvers. The present contribution thus expands the range of problems that can be effectively tackled within probabilistic numerics.

Probabilistic Exponential Integrators
Nathanael Bosch
Tübingen Center, University of Tübingen
nathanael www.hhumiittueiigeeeeeeeeeo
Philipp Hennig
Filipp Tronarp
Tübingen AI Center, University of Tübingen philipp Sereluummntteiinnnnceem
Lund University
filip. f htt.p..rroorrenntttat...... lu.
Abstract
Probabilistic solvers provide a flexible and efficient framework for simulation, uncertainty quantification, and inference in dynamical systems. However, like standard solvers, they suffer performance penalties for certain stiff systems, where small steps are required not for reasons of numerical accuracy but for the sake stability. This issue greatly alleviated in semi-linear problems by the prob- abilistic exponential integrators developed in this paper. By including the fast, linear dynamics in the prior, we arrive at a class of probabilistic integrators with favorable properties. Namely, they proven to be L-stable, and in certain case reduce to classic exponential integrator- with the added benefit of providing a probabilistic account of the numerical error. The method is also generalized arbitrary non-linear systems by imposing piece-wise semi-linearity on the prior via Jacobians of the vector field at the previous estimates, resulting in probabilistic exponential Rosenbrock methods. We evaluate the proposed methods on multiple stiff differential equations and demonstrate their improved stability and efficiency over established probabilistic solvers. The present contribution thus expands the range of problems that can be effectively tackled within probabilistic numerics.
1 Introduction
Dynamical systems appear throughout science and engineering, and their accurate and efficient simulation is key component in many scientific problems. There has also been a surge of interest in the intersection with machine learning, both regarding the usage of machine learning methods model and solve differential equations [36, 18, 35], and in a dynamical systems perspective on machine learning methods themselves [8, 5]. This paper focuses on the numerical simulation of dynamical systems within the framework of probabilistic numerics, which treats the numerical solvers themselves as probabilistic inference methods [11, 12, 33]. In particular, we expand the range of problems that can tackled within this framework and introduce a new class of stable probabilistic numerical methods for stiff ordinary differential equations (ODEs).
Stiff equations problems for which certain implicit methods perform much better than explicit ones [10]. But implicit methods come with increased computational complexity per step, as they typically require solving a system of equations. Exponential integrators are alternative class of methods for efficiently solving large stiff problems [48, 16, 7, 15]. They are based on the observation that, if the ODE has a semi-linear structure, the linear part can be solved exactly and only the non-linear part needs to be numerically approximated. The resulting methods are formulated in an explicit manner and do require solving system of equations, while achieving similar better stability than implicit methods. However, such methods have not yet been formulated probabilistically.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

Explicit method (not A-stable)
b. Semi-implicit method (A-stable)
Exponential integrator (L-stable)
dt= 0.01
0.25
dt= 0.25
3
Probabilistic numerical ODE solvers with different stability properties. Left: The explicit EKO solver with a 3-times integrated Wiener process prior unstable and diverges from the true solution. Center: The semi-implicit EK1 with the same prior does not diverge even though it uses larger step size, due to being Atstable, but it exhibits oscillations in the initial phase of the solution. Right: The proposed exponential integrator is L-stable and thus does not exhibit any oscillations.
Figure 1:
In this paper we develop probabilistic exponential integrators, a new class of probabilistic numerical solvers for stiff semi-linea ODEs. We build on the ODE filters which have emerged as efficient and flexible class of probabilistic numerical methods for general ODEs [40, 21, 45J. They have known convergence rates [21, 46], which have also been demonstrated empirically [2, 26, 24]. they are applicable to wide range of numerical differential equation problems [23, 25, 3], their probabilistic output can be integrated into larger inference problems [20, 39, 47], and they can be formulated parallel-in-time [4]. But while it has been shown that the choice of underlying Gauss-Markov prior does influence the resulting ODE solver [30, 45, 2]. there has not yet been strong evidence for the utility of priors other than the well-established integrated Wiener process. Probabilistic exponential integrators provide this evidence: in the probabilistic numerics framework, "solving the linear part of the ODE exactly" corresponds to an appropriate choice of prior.
Contributions Our main contribution is the development of probabilistic exponential integrators, new class of stable probabilistic solvers for stiff semi-linear ODEs. We demonstrate the close link of these methods to classic exponential integrators in Proposition 1, provide an equivalence result to classic exponential integrator in Proposition 2, and prove their L-stability in Proposition 3. To enable numerically stable implementation, we present quadrature- square-roots of the process noise covariance in Section 3.2. Finally, in Section 3.6 we also propose probabilistic exponential Rosenbrock methods for problems in which remi-linearity not known a priori. We evaluate all proposed methods on multiple stiff problems and demonstrate the improved stability and efficiency of the probabilistic exponential integrators over existing probabilistic solvers.
approach to directly compute
Numerical ODE solutions as Bayesian state estimation
Let us first consider an initial value problem with some general non-linear ODE, of the form
i(() (v((), (), y(0) Yo-
i [0,T],
(la)
(1b)
with vector field f Rd x R + Rd, initial value yo E Rd, and time span [0, Probabilistic numerical ODE solvers aim to compute posterior distribution over the ODE solution y(t) such that it satisfies the ODE on a discrete set of points T = IN C [0, TT, that is
(v(1) v(0) = yo (i((m) = ((v(1,), (.))
We call this quantity, and approximations thereof, a probabilistic numerical ODE solution. Proba- bilistic numerical ODE solvers thus compute not just single point estimate of the ODE solution, but a posterior distribution which provides a structured estimate of the numerical approximation error.
In the following, we briefly recapitulate the probabilistic ODE filter framework of Schober et al. [40] and Tronarp et al. [45] and define the prior, data model, and approximate inference scheme. In Section 3 we build on these foundations to derive the proposed probabilistic exponential integrator.

---

2.1 Gauss-Markoy prior
A priori, we model y(t) with a Gauss Markov process, defined by a stochastic differential equation (3)
(t) = (t) dt + xB (L), with state (t) E Rd(++), model matrices A E [dd B E 1xd diffusion scaling E R. and smoothness q E N. More precisely, A and B are chosen such that the state is structured as (() (t), ..... Y (l)) (t)) and then (t) models the i-th derivative of y(t). The initial value Yo E [Rd(a+1) must be chosen such that it enforces the initial condition, that is, (() yo-
Y()) = Yo:
One concrete example of such probabilistic numerical ODE solvers is the q-times Integrated Wiener process, with model matrices
process that is commonly used in the context of
0
Ia
0
AIwP(d.))
BIWP(d.))
0
La
Alternatives include the class of Matérn processes and the integrated [46]-the latter plays central role in this paper and will be discussed in detail later.
process
(t) satisfies linear Gaussian transition densities of the form [44] (t h) (t) N (O),
(5)
with transition matrices (h) and (h) given by
(h) exp (Ah), Q(h) 0
Φ(h - BBB o6 (h (h-T) -- dr.
(6)
These quantities can be computed with matrix fraction decomposition [44]. For q-times integrated Wiener process priors, closed-form expressions for the transition matrices are available [21].
2.2 Information operator
The likelihood, or data model, of a probabilistic ODE solver relates the uninformed prior to the actual ODE solution of interest with an information operator I [6], defined as
[[Y](t) = E, Y (t) f (EoY ()),
(7)
where E, E [xdx((+11) are selection matrices such that E,Y (t) = ( (t). I[Y] then captures how well Y solves the given ODE problem. In partieular, I maps the true ODE solution y the zero function, L.e. Tly] = 0. Conversely, if Z[yl (t) = 0 holds for all t E [0, T) then solves the given ODE. Unfortunately, it is in general infeasible to solve ODE exactly and enforce I[Y](t) 0 everywhere, which is why numerical ODE solvers typically discretize the time interval and take discrete steps. This leads to the data model used in most probabilistic ODE solvers [45]:
(8) Note that this specific information operator is closely linked to the IVP considered in Eq. (1a). By defining a (slightly) different data model we can also formulate probabilistic numerical IVP solvers for higher-order ODEs or differential-algebre equations, encode additional information such as conservation laws or noisy trajectory observations [3, 39].
(t.) = (t,,, f (t,Y(),)) t,,) = 0, tn E T C [0,
2.3 Approximate Gaussian inference
The resulting inference problem is described by a (1)+ Y(a)) ~ N' ( (tn),,
prior and Dirac likelihood
(9a)
Zn (t.)) ~ (E,Y - (Eo
(9b)
with क;,, := a(t.+1 - t.)) Qn = Q(tn+ttt++ttt f. initial value Y()) = Yo, discrete time grid {t,, J. and zero-valued data Z, 0 for all n. The solution of the resulting non-linear Gauss Markoy regression problem can then be efficiently approximated with Bayesian filtering and smoothing techniques [37]. Notable examples that have been used to construct probabilistic numerical ODE solvers include quadrature filters, the unscented Kalman filter, the iterated extended Kalman smoother, or particle filters [19, 45, 46). Here, we focus on the eell-estallsshed extended Kalman filter (EKF). We briefly discuss the EKF for the given state estimation problem in the following.
3

---

Prediction Given Gaussian state estimate (tx-11 ~ NN (An-1: 1n-1) and the linear conditional
distribution as Eq. (9a), the marginal distribution Y (tn) N (H4 2n also Gaussian, with
Un nnn--n1-e 1--1//-11
(10a)
E = 1 12 n- + k2 Qn--1
(10b)
Linearization To efficiently compute tractable approximation of the true posterior, the EKF linearizes the information operator I around the predicted mean 1 i.e. (4,) H, + bas
Hn = F1 - Fig Fo:
(11a)
b, = F. Poll - f(Folt, l.)
(11b)
An exact linearization with Jacobian Fy = (Fo/4 t.)) leads to a semi-implicit probabilistic ODE solver, which we call the EK1 [45]. Other choices include the zero matrix F, 0, which results in the explicit EKO solver [40, 21], or a diagonal Jacobian approximation (the DiagonalEKI which combines some stability benefits of the EK1 with the lower computational cost of the EKO [24].
Correction step In the linearized observation model, the posterior distribution of Y (tn) given the datum Zn is again Gaussian. Its posterior mean and covariance (Mm; En, are given by
Sn Hn H
(12a)
Kn Σn H ST
(12b)
In I' Kn (E1P (Eopt))
(12c)
En = (T - K, (m-K,Hn) En.
(12d)
This is also known as the update step of the EKF.
Smoothing To condition the state estimates on all data, the EKF can be followed by smoothing pass. Starting with n := JAN and E == Σ,, it consists of the following backwards recursion:
G2 Σn 1 (2
(13a)
5 = Un + Gn (1 - 1n+1):
(13b)
5 = En + Gn DTt1)
(13c)
Result
The above computations result in probabilistic numerical ODE solution with marginals
p (Y(())
- ( =) ~ N (nr 2;
(14)
which, by construction of the state Y, also contains estimates for the ODE solution as y(t) = EgY (t). Since the EKF-bbsed probabilistic solver does not compute only the marginals in Eq. (14), but full posterior distribution for the continuous object ytt), it can be evaluated for times & T (also known as "dense output" in the context of ODE solvers); it can produce joint samples from this posterior; and it can be used as
prior for subsequent inference tasks [40, 2, 47].
2.4 Practical considerations and implementation details
To improve numerical stability and preserve positive-senin semiddfieemi of the computed covariance matrices, probabilistic ODE solvers typically operate on square-roots of covariance matrices, defined
by matrix decomposition of the form M /MVM [26]. For example, the Cholesky factor is one possible square-root of a positive definite matrix. But in general, the algorithm does not require the square-roots to be upper- or lower-triangular, or even square. Additionally, we compute the exact initial state Yo from the IVP using Taylor-mode automatic differentiation [9, 26]. we compute smoothing estimates with preconditioning [26], and we calibrate uncertainties globally with quasi-man muasi-maximum likelihood approach [45, 2].

---

3 Probabilistic exponential integrators
In the remainder of the paper, unless otherwise stated, we focus IVPs with semi-linear vector-field
i(t) = (v((),, t) = Ly(t) + ())
(15)
Assuming N admits Taylor series expansion around t, the variation of constants formula provides formal expression of the solution at time t h:
(
[sp((LLC 7)) 7h dr dt (u((),
d
ytt h) exp(Lh)y(() +
Σ h++ k=0
-T)) k!
(16)
This observation is the starting point for the development of exponential integrators [31, 15]. By further defining the so-called yp-functions
px(2) /. exp(=(1-11) 7)) dr,
xk11 (k 1)!
(17)
the above identity the ODE solution simplifies to
de
ytt + h) = exp(()))))))
En
nht P (Lh) de ())()),))
(18)
In this section we develop class probabilistic exponential integrators. This is achieved by defining an appropriate class of priors that absorbs the partial linearity, which leads to the integrated Canstinn-Uhleneee processes. Proposition below directly relates this choice of prior to the classical exponential integrators. Proposition demonstrates a direct equivalence between the predictor-corrector form exponential trapezoidal rule and the once integrated Ornstein CeasstinnUhllnbeck process. Furthermore, the favorable stability properties of classical exponential integrators is retained for the probabilistic counterparts as shown Proposition 3.
3.1 The integrated Ornstein-Uhlenbers process
In Section 2.1 we highlighted the choice of the g-times integrated Wiener process prior, which essentially corresponds modeling the (q- 1)-th derivative of the right-hand side f with Wiener process. Here we follow a similar motivation, but only for the non-linear part N. Differentiating both sides of Eq. (15) q 1 times with respect to yields
d
d d(9- u(t
d N((()))
d - (t)
(19)
Then, modeling
d9-1 ()), as Wiener process and relating the result to y(() gives dea-T
dy (t) (() dl,
(20a)
( (1) Ly ( (1) di KT, dw (1).
(20b)
This process is also known as the q-times integrated Orhstein-Uhlenbeck process (IOUP), with rate parameter L and diffusion parameter k. It can be equivalently stated with the previously introduced notation (Section 2.1), by defining state (f), as the solution of a linear time- inveeiinaiiant (LTI) SDE as in Eq. (3), with system matrices
La
0
AlouP(d.g)
BioUP(d.g)
(21)
Ia L
Ia
Remark 1 (The mean of the IOUP process solves the linear part of the ODE exactly). By taking the expectation of Eq. (20b) and by linearity of integration,
see that the mean of the IOUP satisfies
( (t) Lp (t).
(0) yo-
(22)
This is line with the motivation of exponential integrators: linear part of the ODE is solved exactly, and only need to approximate the non-linear part. Figure 2 visualizes this idea.
5

---

3 Integrated Wiener process
b. Integrated
IOUP initial value
3K
Figure 2: Damped oscillator dynamics and priors with different degrees of encoded information. Left: Once-integrated Wiener process, a popular prior for probabilistic ODE solvers. Center: Once- integrated Omstein-Uhlenbeck process (IOUP) with rate parameter chosen to encode the known linearity the ODE. Right: IOUP with both the ODE information and a specified initial value and derivative. This is the kind of prior used in the probabilistic exponential integrator.
3.2
The transition parameters of the integrated Ornstein-Uhlenbed process
Since the process Y (t) is defined the solution of a linear time-invariant SDE, satisfies discrete transition densities p (1 h) Y()) . (b)) (t), =Q(h).. The following result shows that the transition parameters are intimately connected with the ş-functions defined in Eq. (17).
Proposition 1. The transition matrix of q-times integrated Orustein-Uhllnbeek process satisfies
h
Texp (exp(Awwidggn
P12(h) ] exp( (h)
with P12(h) ha P
(hh
(23)
(L)(
Proof in Appendix A. Although Proposition 1 indicates that 4 (h) may be computed more efficiently than by calling matrix eearix--xppoeential dl(g + d(y + 1) matrix, this is known be numerically less stable [41]. We therefore compute ((h) with the standard matrix-exponential formulation.
Directly computing square-roots of the process noise covariance Numerically stable probabilistic ODE solvers require square-root, /Q(h), of the process noise covariance rather than the full matrix, Q(h). For IWP priors this can be computed from the closed-form representation of Q(h) via an appropriately preconditioned Cholesky factorization [26]. However, for IOUP priors have found an analogous method that works reliably. Therefore, compute (Qhh) directly with numerical quadrature. More specifically, given a quadrature rule with nodes T, E [0, /s] and positive weights wi 0, the integral for Qh) given in Eq. (6) approximated by
X71
m
Q(h)=
wi exp(A(h BBT exp(A (h (()-T)))
2 M,
(24)
with square-roots VM, ww, exp( A((h BB of the summands, which is well-defined since w; > O. We can thus compute a square-root representation of the sum with a QR-decomposition
X. QR (IVM,
(25)
We obtain Q(h) Re R,R, and therefore approximate square-root factor is given by /Q(h) R
Similar ideas have previously been employed for time integration of Riccati equations [42, 43]. We this quadrature-trick for all IOUP methods, with Causs-Legendre quadrature on m 7 nodes.
3.3 Linearization and correction
The information operator of the probabilistic exponential integrator defined exactly as in Sec- tion 2.2. But since we now assume semi- linear vector-field / f. we have an additional option for the linearization: instead of choosing the exact F,, = J (EK1) or the zero-matrix F,, = (EKO), cheap approximate Jacobian is given by the linear part F,, L. We denote this approach by EKL. This is chosen as the default for the probabilistic exponential integrator. Note that the EKL approach can also be combined with an IWP prior, which will be serve additional baseline in the Section 4.

---

3.4 Equivalence to the classic exponential trapezoidal rule in predict-evallate-cooree mode
Now that the probabilistic exponential integrator has been defined, we can establish an equivalence result to classic exponential integrator, similarly to the closely-related equivalence statement by Schober et al. [40, Proposition 1] for the non-exponential case.
Proposition (Equivalence to the PEC exponential trapezoidal rule). The mean estimate of the probabilistic exponential integrator with once-integrated parameter I is equivalent to the classic exponential trapezoidal rule in eeeee-aeeeleeerreeeee correct mode, with the predictor being the exponential Euler method. That is, is equivalent the scheme
prior with rate
Un+1 = (o((L)) + (L))
(26a)
Vn+1 = po((h)) + (L) + (2hh) N ((n++)) -()))
(26b)
h
where Eq. (26a) corresponds to prediction step with the exponential Euler method, and Eq. (26b)
corresponds to a correction step with the exponential trapezoidal rule.
The proof is given in Appendix B. This equivalence result provides another theoretical justification for the proposed probabilistic exponential integrator. But note that the result only holds for the mean, while the probabilistic solver computes additional quantities order to track the solution uncertainty, namely covariances. These are not provided by classic exponential integrator.
3.5 L-stability of the probabilistic exponential integrator
When solving stiff ODEs, the actual efficiency of numerical method often depends on its stability, One such property A-stability: It guarantees that the numerical solution of a decaying ODE will also decay, independently of the chosen step size. In contrast, explicit methods typically only decay for sufficiently small steps. In the context of probabilistic ODE solvers, the EKO is considered to be explicit, but the EK1 with IWP prior has been shown to be A-stable [45]. Here, we show that the probabilistic exponential integrator satisfies the stronger I-stability: the numerical solution not only decays, but it decays fast, i.e. it goes different probabilistic solver stabilities. For formal definitions, see for example [27, Section 8.6].
the step size goes to infinity. Figure visualizes the
Proposition 3 (L-stability). The probabilistic exponential integrator is L-stable.
The full proof is given in Appendix C. The property essentially follows from Remark which stated that the IOUP solves linear ODEs exactly. This implies fast decay and gives L-stability.
3.6 Probabilistic exponential Rosenbrock-type methods
We conclude with a short excursion into exponential Rosenbrock methods [14, 17, 28]: Given non- linear ODE j(t) fl((,, t), exponential Rosenbrock methods perform continuous linearization of the right-hand side f around the numerical ODE solution and essentially solve a sequence of IVPs
j(t) = J.y(t) + (f(v(t).t) - J,,(t)), yt,,,) Yn
t E
(27a)
(27b)
where Jn is the Jacobian of f at the numerical solution estimate j(t,). This approach enables exponential integrators for problems where the right-hand side f is not semi-linear. Furthermore, by automatically lincarizing along the numerical solution the lincarization can be more accurate, the Lipschitz-constant of the non-linear remainder becomes smaller, and the resulting solvers can thus be more efficient than their globally lincarized counterparts [17].
This can also be done in the probabilistic setting: By linearizing the ODE right-hand side f at each step of the solver around the filtering mean Euja, we (locally) obtain a semi-linear problem. Then, updating the rate parameter of the integrated OmntteinUUlleneec process at each step of the numerical solver results in probabilistic exponential Rosenbrock -type methods. As before, the linearization of the information operator can be done with any of the EKO, EK1, or EKL. But since here the prediction relies on exact local linearization, we will by default also use an exact EK1 linearization. The resulting solver and its stability and efficiency will be evaluated in the following experiments.
7

---

10
10
10
10
10
Final error
EKO & IWP(2)
EKI & IWP(2)
10 10-10-
EKL IWP(2)
EKL & (OUP(2)
10-00
10
100
100
100
10
100
Number of steps
Number of steps
Number of steps
Number of sleps
Figure 3: The IOUP prior is more beneficial with increasing linearity of the ODE. In all three examples, the IOUP-based exponential integrator achieves lower error while requiring fewer steps than the IWP-based solvers. This effect is more pronounced for the more linear ODEs.
4 Experiments
In this section we investigate the utility and performance the proposed probabilistic exponential integrators and compare to standard non-exponential probabilistic solvers multiple ODEs. All methods are implemented in the Julia programming language [1], with special care being taken to implement the solvers in numerically stable manner, that is, with exact state initialization, preconditioned state transitions, and a square-root implementation [26]. Reference solutions are computed with the Differentiallquations.jj package [34]. All experiments level CPU. Code for the implementation and experiments is publicly available on GitHub.'
single, consumer-
4.1 Logistic equation with varying degrees of non-linearity
We start with a simple one-dimensional initial value problem: a logistic model with negative growth rate parameter and carrying capacity E R4,, of the form
(()) = v(() + K u()?, (0) = 1.
E 10).
(28a)
(28b)
The non-linearity of this problem can be directly controlled through the parameter K. Therefore, this test problem lets us investigate the IOUP's capability to leverage known linearity in the ODE.
We compare the proposed exponential integrator all introduced IWP-based solvers, with differ- ent linearization strategies: EKO approximates 0,f Rs (and is thus explicit), EKL approximates 1,
0yf= -1, and EK1 linearizes with the correct Jacobian f. The results for four different values of K are shown in Fig. 3. The explicit solver shows the largest error of all compared solvers, likely due to its lacking stability. On the other hand, the proposed exponential integrator behaves expected: the IOUP prior most beneficial for larger values of K, and as the non-linearity becomes more pronounced the performance of the IOUP approaches that of the IWP-based solver. Though for large step sizes, the IOUP outperforms the IWP prior even for the most non-linear case with K 10.
4.2 Burger's equation
Here, consider Burger's equation, which is semi linear partial differential equation (PDE)
0,((,, t) = D0Z t) - u(z, t),, t), E [0, 1].
t E [0, 1],
(29)
with diffusion coefficient D E R We transform the problem into semi-linear ODE with the method of lines [29, 38], and discretize the spatial domain on 250 equidistant points and approximate the differential operators with finite differences. The full IVP specification, including all domains, initial and boundary conditions, and additional information on the discretization, is given in Appendix D.
The results shown in Fig. 4 demonstrate the different stability properties of the solvers: The explicit EKO with IWP prior is unable to solve the IVP for any of the step sizes due its insufficient stability, and even the A-stable EKI and the more approximate EKL require small enough steps At < 10-1 On the other hand, both exponential integrators are able to compute meaningful solutions for larger range of step sizes. They both achieve lower enrors for most settings than their non-exponential
integratoor

---

ODE solution
b. Work-precision diagram
0.5
EKO IWP(2)
Final error
10
EK & TWP(2)
Time
EKL & IWP(2)
EKL & (OUP(2)
10
EK1 & IOUP(2) (RB)
0
0.5
10
10
10a
10
10
10
10a
Space
Step size
Number of evaluations
Figure 4: Benchmarking probabilistic ODE solvers on Burger's equation. Exponential and non- exponential probabilistic solvers are compared on Burger's 's equation (a) in two work-precision diagrams (b). Both exponential integrators with IOUP prior achieve lower errors than the existing IWP-based solvers, in particular for large steps. This indicates their stronger stability properties.
a. ODE solution
Work-precision diagram
10
EKO IWP(2)
Time
Final emmr
10
EKI IWP(2)
0.5
EKT. & IWP(2)
10
EKL IOUP(2)
EKI & IOUP(2) (RB)
10
10"
10
10
Space
Step size
Runtime [s]
peaction-diffusiin model. Exponential and non-exponential probabilistic solvers are compared on a -iiftion model (a) two work- precision diagrams
Figure 5: Benchmarking probabilistic ODE solvers
(b). The proposed exponential integrators with IOUP prior achieve lower errors per step size than the existing IWP-based methods. The runtime comparison shows the increased cost of the Rosenbrock-type (RB) method, while the RoseeRosnbrock probabilistic exponential integrator performs best in this comparison.
counterparts. The second diagram in Fig. 4 compares the achieved error the number of vector-field evaluations and points out a trade-off between both exponential methods: Since the Rosenbrock method additionally computes two Jacobians (with automatic differentiation) per step. needs to evaluate the vector-ficld more often than the non- non-Rosenbrock method. Thus, for expensive-to-evalute vector fields the standard probabilistic exponential integrator might be preferable.
4.3 Reaftion-diffusion model
Finally, we consider a discretized reaction-diffusion model given by semi-linear PDE
au((,, 1) = (0 + R(u(x, ()), x E [0, 1],
1 LE (0,
(30)
where E R4 is the diffusion coefficient and R(u) u(1 - u) is a logistic reaction term [22]. A fiite-difference discretization of the spatial domain transforms this PDE into an IVP with semi-linear
ODE. The full problem specification is provided in Appendix D.
Figure 5 shows the results. We again observe the improved stability of the exponential integrator variants by their lower error for large step sizes, and they outperform the IWP-based methods on all settings. The anatime-evaluation in Fig. 5 also visualizes another drawback of the Rosenbrock-type method: Since the problem is re-lincarized at each step. the IOUP also needs to re-discretized and thus matrix exponential needs to be computed. In comparison, the non- Roneooobbbock method only discretizes the IOUP prior once at the start the solve. This advantage makes the non-Rosent non-Rosenbrock probabilistic exponential integrator the most performant solver in this experiment.
Limitations
The prohabilistic exponential integrator shares many properties of both classic exponential integrators and of other iltering-basee probabilistic solvers. This also brings some challenges.

---

Cost of computing matrix exponentials The IOUP prior is more expensive to discretize than the IWP as it requires computing a matrix exponential. This trade-off well-known also in the context of classic exponential integrators. One approach to reduce computational cost is to compute the matrix exponential only approximately [32], for example with Krylov-subspace methods [13, 17]. Extending these techniques to the probabilistic solver setting thus poses an interesting direction for future work.
Cubic scaling in the ODE dimension The probabilistic exponential integrator shares the com- plexity most semi-)implicci ODE solvers: while being linear in the number of time steps, it scales cubically in the ODE dimension. By exploiting structure in the Jacobian and in the prior, some filtering-bast based ODE solvers have been formulated with linear scaling in the ODE dimension [24]. But this approach does not directly extend to the IOUP-prior Nevertheless, exploiting known structure could be particularly relevant construct solvers for specific ODEs, such as certain discretized PDEs.
Conclusion
We have presented probabilistic exponential integrators, a new class of probabilistic solvers for stiff semi-linea linear ODEs. By incorporating the fast, linear dynamics directly into the prior of the solver, the method essentially solves the linear part exactly, in similar manner as classic exponential integrators. We also extended the proposed method to general non-linear systems via iterative re-linearization and presented probabilistic exponential Rosenbrock-tye
methods. Both methods have been shown both theoretically and empirically to be more stable than their non-exponential probabilistic counterparts. This work further expands the toolbox of probabilistic numerics and opens up new possibilities for accurate and efficient probabilistic simulation and inference in stiff dynamical systems.
Acknowledgments and Disclosure of Funding
The authors gratefully acknowledge financial support by the German Federal Ministry of Education and Research (BMBF) through Project ADIMEM (FKZ LIIS18052B). and financial support by the European Research Council through ERC StG Action 757275 / PANAMA; the DFG Cluster of Excellence Machine Learning New Perspectives for Science, EXC 2064/1, project number 390727645; the German Federal Ministry of Education and Research (BMBF) through the Tübingen Al Center (FKZ: 011S 8039A); and funds from the Ministry of Science, Research and Arts of the State of Baden-Württemlle Württemberg. Filip Tronarp was partially supported by the Wallenberg AI, Autonomous Systems and Software Program (WASP) funded by the Knut and Alice Wallenberg Foundation. The authors thank the International Max Planck Research School for Intelligent Systems (IMPRS-IS) for supporting Nathanael Bosch. The authors also thank Jonathan Schmidt for many valuable discussions and for helpful feedback on the manuscript.
References
[1] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah. Julia: A fresh approach to numerical computing. SIAM review, 59(1):65 98, 2017.
[2] N. Bosch, P. Hennig, and Tronarp. Calibrated adaptive probabilistic ODE solvers. In International Conference Artificial Intelligence and Statistics. PMLR, 2021.
[3] N. Bosch, F. Tronarp, and P. Hennig. Pick-and-mix information operators for probabilistic ODE solvers. In International Conference on Artificial Intelligence and Statistics. PMLR, 2022.
[4] N. Bosch, A. Corenflos, F. Yaghoobi, F. Tronarp, P. Hennig, and S. Särkkä, Parallel-in -nn--inne probabilistic numerical ODE solvers, 2023.
[5] R. T. Q. Chen, Y. Rubanova, J. Bettencourt, and D. K. Duvenaud. Neural ordinary differential equations. In Advances in Neural Information Processing Systems. Curran Associates, Inc., 2018.
[6] Cockayne, C. Oates, T. Sullivan, and M. Girolami. Bayesian probabilistic numerical methods. SIAM Review, 61:756-789, 2019.
10

---

[7] S. M. Cox and P. C. Matthews. Exponential time differencing for stiff systems. Journal of Computational Physics, 176(2):430-455 2002.
[8] W. E. proposal on machine learning via dynamical systems. Communications in Mathematics and Statistics, 5(1):111, Mar 2017.
[9] A. Griewank and A. Walther. Evaluating Derivatives. Society for Industrial and Applied Mathematics, second edition, 2008
[10] E. Hairer and G. Wanner. Solving Ordinary Differential Equations 1I: Stiff and Differential- Algebraic Problems. Springer series in computational mathematics.
SpriinngeressSrineeee-ee Apringer-Verlag. 1991.
[11] Hennig, M. A. Osborne, and M. Girolami. Probabilistic numerics and uncertainty in computations. Proceedings. Mathematical, physical, and engineering sciences, 471(2179): Jul 2015.
[12] Hennig, M. A. Osborne, and H. Kersting. Probabilistic Numerics: Computation as Machine Leurning. Cambridge University Press, 2022.
[13] M. Hochbruck and C. Lubich. On Krylov subspace approximations to the matrix exponential operator. SIAM Journal on Numerical Analysis, 34(5):19111955, Oct 1997.
[14] M. Hochbruck and A. Ostermann. Explicit integrators of Rosenbrock-type. Oberwolfach Reports, 3(2): 100-11100111,,, 2006.
[15] M. Hochbruck and A. Ostermann. Exponential integrators. Acta Numerica, 19:209286, 2010.
[16] M. Hochbruck, C. Lubich, and H. Selhofer. Exponential integrators for large systems of differential equations. SIAM Journal on Scientific Computing, 19(5):15521554, 1998.
[17] M. Hochbruck, A. Ostermann, and J. Schweitzer. Exponential Rosenbrock-type methods. SIAM Journal on Numerical Analysis, 47(1):786 803, 2009.
[18] G. E. Karniadakis, G. Kevrekidis, I.. Lu, Perdikaris, Wang, and L. Yang. Physics- informed machine learning. Nature Reviews Physics, 3(6):422440, May 2021.
[19] H. Kersting and P. Hennig. Active uncertainty calibration in Bayesian ode solvers. In Proceed- ings of the 32nd Conference on Uncertainty in Artificial Intelligence (UAI), pages 309-318, June 2016.
[20] H. Kersting, N. Krämer, M. Schiegg, C. Daniel, M. Tiemann, and P. Hennig. Differentiable likelihoods for fast inversion of Likelihood-free dynamical systems. In International Conference on Machine Learning. PMLR, 2020.
[21] H. Kersting, T. J. Sullivan, and P. Hennig. Convergence rates of Gaussian ODE filters. Statistics and computing, 1791))1791-1818 2020.
[22] A. N. Kolmogorov. A study of the equation of diffusion with increase in the quantity of matter, and its application a biological problem. Moscow University Bulletin of Mathematics, 1:1-25, 1937.
[23] N. Krämer and P. Hennig. Linear-time probabilistic solution of boundary value problems. In Advances in Neural Information Processing Systems. Curran Associates, Inc., 2021.
[24] N. Krämer, N. Bosch, J. Schmidt, and P. Hennig. Probabilistic ODE solutions in millions of dimensions. In International Conference on Machine Leurning. PMLR, 2022.
[25] N. Krämer, J. Schmidt, and P. Hennig. Probabilistic numerical method of lines for time- dependent partial differential equations. In International Conference on Artificial Intelligence and Statisties. PMLR, 2022.
[26] N. Krämer and Hennig. Stable implementation of probabilistic ode solvers, 2020.
[27] D. Lambert. Computational Methods Ordinary Differential Equations. Introductory Mathematics for Scientists And Engineers. Wiley, 1973.
11

---

[28] V. T. Luan and A. Ostermann. Exponential Rosenbrock methods of order five construction, analysis and numerical comparisons. Journal of Computational and Applied Mathematics, 255: 417-431, 2014.
[29] N. K. Madsen. The method of lines for the numerical solution of partial differential equations. Proceedings of the SIGNUM meeting on Software for partial differential equations, 1975.
[30] E. Magnani, H. Kersting, M. Schober, and P. Hennig. Bayesian filtering for ODEs with bounded derivatives, 2017.
[31] B. V. Minchev and W. M. Wright. A review of exponential integrators for first order semi-linear problems. Technical report, Norges TekniiskNNaturvettnssaalli Universitet, 2005.
[32] C. Moler and C. Van Loan. Nineteen dubious ways to compute the exponential of matrix, twenty-five years later. SIAM Review, 45(1):349 Jan 2003.
[33] C. J. Oates and J. Sullivan. A modern retrospective on probabilistic numerics. Statistics and Computing, 29(6): 1300:13333111111 2019.
[34] C. Rackauckas and Q. Nie. DifferentialEquations for solving differential equations in julia. Journal of Open Research Software, 5(1), 2017.
performant and feature-rich ecosystem
[35] C. Rackauckas, Y. Ma, J. Martensen, C. Warner, K. Zubov, R. Supekar, Skinner, A. Ramadhan, and A. Edelman. Universal differential equations for scientific machine learning, 2021.
[36] M. Raissi, P. Perdikaris, and G. Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 378:686707, Feb 2019.
[37] S. Särkkä. Bayesian Filtering and Smoothing, volume 3 of Institute of Mathematical Statistics lextbooks. Cambridge University Press, 2013.
[38] W. E. Schiesser. The numerical method of lines: integration of partial differential equations. Elsevier, 2012.
[39] J. Schmidt, N. Krämer, and P. Hennig. A probabilistic state space model for joint inference from differential equations and data. In Advances in Neural Information Processing Systems. Curran Associates, Inc., 2021.
[40] M. Schober, S. Särkkä, and P. Hennig. A probabilistic model for the numerical solution of initial value problems. Statistics and Computing, (9(1)999-12, Jan 2019.
[41] R. B. Sidje. Expokit: A software package for computing matrix exponentials, ACM Transactions on Mathematical Software, 24(1):13015 mar 1998.
[42] Stillfjord. Low-rank second-order splitting of large-scale differential Riceati equations. IEFE Transactions on Automatic Control, 60(10))2791 1112279999777 2015.
[43] T. Stillfjord. Adaptive high-order splitting schemes for large-scale differential Riccati equations. Numerical Algorithms, 78(4): 1291151, Sep 2017.
[44] S. Särkkä and A. Solin. Applied Stochastic Differential Equations. Institute of Mathematical Statistics Textbooks. Cambridge University Press, 2019.
[45] F. Tronarp, H. Kersting, S. Särkkä, and P. Hennig. Probabilistic solutions to ordinary differential equations as nonlinear Bayesian filtering: new perspective. Statistics and Computing, 29(6): 1297-1315, 2019.
[46] F. Tronarp, S. Särkkä, and P. Hennig. Bayesian ODE solvers: The maximum a posterioni estimate. Statistics and Computing. 31(3):1-18, 2021.
[47] F. Tronarp, N. Bosch, and P. Hennig. Fenrir: Physics-enhanced regression for initial value problems. In International Conference on Machine Learning. PMLR, 2022.
[48] C. Van Loan. Computing integrals involving the matrix exponential. IEEE Transactions on Automatic Control, 23(3):395-400 1978.
12

---

Probabilistic Exponential Integrators - Appeendix
Proof of Proposition 1: Structure of the transition matrix
Proof of Proposition I. The drift-matrix ALOUP(d,g) given in Eq. (21) has block structure
AIwp(dg-1) E Eq-1 L
Asoup(d,0) =
(31)
where Eq-1== [0 0 Id] T
E Rddxdd From Van Loan [48, Theorem IJ, it follows
(h) cxp( AIwP(d.--))))
12)
(32)
0
exp(Lh)]
which is precisely Eq. (23). The same theorem also gives (h) as
12h L. exp(Awvidg-yy( r)) (d-1) exp(Lr)
(33)
Its ith d block is readily given by
/ E. exp((Amp(dg-1)( - E exp(Lr) dr
(P12(h));
1/¹ h (h - -)9-1-i exp exp(Lr) LT) dr
(4 -i!
(34)
hタ ('
79-1-i
1--1-ij exp(Lh(( - ()) T))
hh (gh))))))
where the second last equality used the change of variables T = 1-(11()) ([]-.). u). and the last line follows by definition.
B Proof of Proposition 2: Equivalence to classic exponential integrator
We first briefly recapitulate the probabilistic exponential integrator setup for the case of the once integrated OrsstennUUhe Uhlenbeck process, and then provide some auxiliary results. Then, we prove Proposition 2 in Appendix B.3.
B.1
The probabilistic exponential integrator with once-integrated Ornstein-Uhlenbedd prior
The integrated Crinstein-Uhlentee process prior with rate parameter L results in transition densities with transition matrices (from Proposition 1)
Y + h) Y(t) NV (()
T hup (I.h) s0 (T))
((h) exp(Ah)
(35)
Q(h) /" exp(Ar)BA BB T exp(AA dT
(36)
d [ (LT) Lo(LL))
0
(LT(LT)] T
po(LT)) dr
(37)
0
(LT)
dr,
(38)
Pu(Lr)po(LT))
13

---

where we assume unit diffusion o² = 1. To simplify notation, we assume equidistant time grid [tn]mg n=o with in n h for some step size and we denote the constant transition matrices simply by Φ and Q and write Yn = (4.).
Before getting to the actual proof, let also briefly recapitulate the filtering formulas that computed at each solver step. Given a Gaussian distribution Yn ~ N' (Yni Pay 2n), the prediction step computes
Pa+1 Ppn
(39)
Σ Ext1 (h) QQh).
(40)
Then, the combined lincarization and correction step compute
2n+1 E1 f( (Eulftt)>
(41)
S 1Σ
(42)
Ka+1 = 2u+1 H S
(43)
Pn+1 = Pn+1 Kn+i in+1l
(44)
Ex+1 E 54t1
Kn+1S K n+1
(45)
with observation matrix H/ = LE [-L
I, since we perform the proposed EKL linearization.
B.2 Auxiliary results
In the following, we show some properties of the transition matrices and the covariances that will be needed in the proof of Proposition 2 later.
First, note that by defining (o(2)) exp 2, the f-funcctions satisfy the following recurrence formula:
z4p(22) (z) = (z) (k 1)!
(46)
See e.g. Hochbruck and Ostermann [15]. This property will be used throughout the remainder of the section.
Lemma B.1. The transition matrices (h), Q(h) of the once integrated Ornstein-Uhlenbeek process with rate parameter L satisfy
H = [-L I),
(47)
Q(h)a = (LL))) hipi (Lh)
(48)
HQ(h) = hel,
(49)
Proof.
1
1
hupi (Lh)
HT(h) (E1 - LEo)
= 10 po([h)] - L hep(LLh) = [-L 4).
(50)
0
Po(
T (LT) (()(LT)PU(LT)) TPo(LT)q (LT) Fo(LT)PO(LT)
Q(h)H
H HTdr dr
(51)
(LT) - (LT
dr
(52)
Pu(LT) - (rFo(LT)(((LT)
(TT (Po( - -IT)
dr
(53)
(o([T) (Po(LT)T -ITI((TT
TRI (LT) (o(LT)
dr
(54)
[22 (Th) (pp)(h))
(55)
14

---

where we used LTP) pp((T) I, and 0 [*]))[]]]] = --h-1 (LT)) It follows that
y22 (Lhh) h1p1(1 (T.h))
HQ(H)
= (Ti((T))) - Ihupz(THii) = hT,
(56)
where we used ITT2((TT) (I([T) -
11
Lemma B.2. The prediction covariance Lntz satisfies
E H = Q(1)H
(57)
Proof. First, since the observation model is noiseless, the filtering covariance En satisfies
HΣ, = [0 0).
(58)
This can be shown directly from the correction step formula:
НΣn HΣ, - S,,KS,
(59)
НΣ H (Σ,, H ST²) S.K,
(60)
H H 11 (HE, S
(61)
НΣ, IS
(62)
НΣ Sn (E H 1) НΣ, ST HH
(63)
(64)
[0
(65)
Next, since the observation matrix is H = L-L I]. the filtering covariance E,, is structured as
Σn
Din [I L LJ..
(66)
L
This can be shown directly from Eq. (58):
Loo
Eo1
[0 0j -Σ= I-
[No
-20o
En L2O1]
(67)
Mio
211
and thus
Σi0 = LEpo=
(68)
Σn IE I.T. I.E
(69)
It follows
Lgo EooL T
L LE
= 1 Eoo [ L
Σ
(70)
Finally, together with Lemma B.1 we can derive the result:
E H +(h) E, (h)
H
+ Q(h)HT
(71)
(h)
Σ LT
+ Q(OHTT
(72)
+hhh
Σn 0 - Q
(73)
T
= Q(h)
(74)
15

---

B.3 Proof of Proposition
With these results, we can now prove Proposition 2.
Proof of Proposition 2. We prove the proposition by induction, showing that the filtering means are all of the form
Un Un = Lyn N (in))
(75)
where Un. Un are defined
Uo :- yo,
(76)
Un+1 P0 (Th)yn + h (Th))
(77)
Ma++ == so (h)ynn + hupi (Lh) (ūa)) - hup2(Lh) (N N ((n++)))
(78)
This result includes the statement of Proposition
Base case n 0
The initial distribution of the probabilistic solver is chosen as
Yo ]
μi0 = Lyoo + (ūo)
2o
(79)
This proves the base case 0.
Induction step n + n + 1 Now, let
Pn
Yn Lyn + N (ūn)
(80)
be the filtering mean step and En be the filtering covariance. The prediction mean is of the form
yn hup (Tu) + N(O.)) so (Th))) + N(r))) = [Tpo(T))y + hpp (L.)N().) (T. + ((,))
M = @hh)pn
(81)
The residual 2n+1 then of the form
ên+1 FI ((ooT4+1)
(82)
=o((L))(Lun + N ((n)) (po(Lh))n + hpp (LA)
(83)
=Po(L))(LLn, N(T))) L (Poo(oo))))))))) (L)) N (ūm-1)
(84)
-Po(LH) + ( - LepoLhh)u ((.)) - N (īn-1)
(85)
(po(Lh) -hipi (Lh)) (T N (2n+1) ((T))) N (fn+1),
(86)
(87)
(88)
where we used properties of the funuctionss, namely Lh (Lh) = Po(Lh) and the commutativity
(Lh(L)) Lipo(Lh) With Lemma B.2, the residual covariance S,+1 and Kalman gain Kn+1 are then of the form
Sh+1 = 11 H HQ(N)U =h1,
(89)
Kn+1 = 1nth HT HT S- = Q(h)H (h (h)2)(hh)
= (Lh)
(90)
This gives the updated mean
pn+1
#n+t - Kntt22+1
(91)
(oo(Lh) + hip (Lh) (L( + ((w)))
[L2) (N(Tm) - N(
(92)
(Lh)
[PPoLLh),, + hup (Lh)AN -hip)) (N -N
(93)
:P( (P((L) + N(jn)) - -P2( (NV -N((n+1)
16

---

This proves the first half of the mean recursion:
Bofl+11 (Lhh) + (upi(Lh)
- hhp2(Lh) (N - N ((n++1)
4n+1-
(94)
It is left to show that
E1p = Lyo+1 - N (In+1).
(95)
Starting from the right-hand side, we have
Lyn+1 + (Bn+1)
(96)
L. (Po(L))) + (
(1p2(LN) (N( ((h++1))) + N(
(97)
LP((L))L, + Lhipi (Lh) - Lhup (Lh) - (N( ((
(98)
=Po(LL) + (go(Lh)-I)N(((b) (vI1 - (N (N(n))))())))) ((
(99)
+ -N(.)) (FI(E)((((()) -N(ū++1) (Lyn
(100)
E1jm+11
(101)
This concludes the proof of the mean recursion and thus shows the equivalence of the two recursions.
C Proof of Proposition 3: L-stability
We first provide definitions L-stability and A-stability, following [27, Section 8.6].
Definition (L-stabillity) A one-step method is said to be L-stable if is A-stable and, in addition, when applied to the scalar test-equation i(() = Ay((),, X E C a complex constant with Re(X) < 0, it yields Vz+1 R(h))yn, and R(hA) - as Re(hX) - -c.
Definition (A-stability). one-step method is said be A-stable if its region of absolute stability contains the whole of the left complex half-plane. That is, when applied to the scalar test-equation complex constant with Re(A) < the method yields yn+1 = R(hX)ym.
((t) = Xy(t) with and {= C: Re(2) < 0} C {z E: R(z) < 1}.
Proof of Proposition 3. Both L-stability and A-stability directly follow from Remark 1: Since the probabilistic exponential integrator solves linear ODEs exactly its stability function is the exponential function, i.e. R(z) exp(2). A-stability and L-stability then follow: Since C- {2: |R(2) s 1} holds the method A-stable. And since [R()] - 0 as Re(2) + -o the method is L-stable. O
D Experiment details
D.1 Burger's equation
Burger's equation is semi-linear partial differential equation (PDE) of the form 0,u(z,t) t) -u(x, t)0,,(z, t) + D0Lu(x,t),
E t [0,T],
(102)
with diffusion coefficient D R4. We discretize the spatial domain Ω finite grid and approxi- mate the spatial derivatives with finite differences to obtain a semi-linear ODE the form
(t) = D L- y(t) + F(y()), t [0,T],
(103)
with N-dimensional v(1) E R.V. L E RNN the finite difference approximation of the Laplace operator 0g, and non-linear part F.
More specifically, we consider domain f2 (0, 1), which we discretize with a grid of N 250 cquidistant locations, thus we have Ar = 1/N. We consider zero-Dirichlet boundary conditions, that u(0, t) u(1,t) The discrete Laplacian is then
2
if = j.
1
LLis
if = otherwise.
(104)
A
17

---

The non non-lineat linear part of the discretized Burger's equation results from another deite-differeece approx- imation the term u 2, and chosen as
1, i d, else.
[F()): = 4
V
(105)
Vi
The initial condition is chosen
u(2, 0) = sin (3TT)3))) xz)3/²
(106)
We consider an integration time-span E [0. 1]. and choose a diffusion coefficient D 0.075.
D.2 Reaction-diffusion model
The reaction- reaction-diffusion model presented in the paper, with logistic reaction term, has been used to describe the growth and spread of biological populations [22]. is given by semi-linear PDE
0e((,, D0Z t) + R(u(x,t))
t E [0,
(107)
where D E R+ is the diffusion coefficient and R(u) u(1 u) is a logistic reaction term. We discretize the spatial domain f on finite grid and approximate the spatial derivatives with finite differences, and obtain semi-linear ODE of the form
(t) = LL)y(() + R(v(t)),
t E [0,
(108)
with N-dimensional y(1) E RR,, operator, and the reaction term R is as before but applied element-wise.
RNXN the finite difference approximation of the Laplace
We again consider domain Ω (0, 1), which we discretize on a grid of N 100 points. This time we consider zero-Neumann conditions, that is, 0zu(0, t) = a,, = 0. Including these directly into the -iffe-difference discretization, the discrete Laplacian then
if = j = orr = j= d,
2
if = j.
D
(109)
Az?
1
if = 1 otherwise.
The initial condition chosen
u(x, 30 10
(110)
The discrete ODE then solved on time-span E [0, 2], and we choose a diffusion coefficient
D 0.25.
18


Reviews:
[{"review_id": "1lw5OFqI20", "reviewer": "Reviewer_goRD", "summary": "The paper studies classical exponential integrators using the framework of probabilistic numerics, where Bayesian formalism is used to study numerical approximations of deterministic dynamical systems. The main result seems to be Proposition 1, which establishes exponential stability (called L-Stability for some reason) property for the exponential integrator.", "questions": "- Calling the Ornstein-Uhlenbeck process a Gauss-Markov in eqn 3 is overkill. \n- Notation in eqn 5 of normal distribution not clear what. Why three parameters? \n- Section 2.2 I don't see why \\mathcal I [y] = 0  when y is a solution to ODE. This says that y_t-f(y_t,t)=0 not \\frac{d}{dt} y_t -f(y_t,t)=0\n- Explain the role of Z in 9b\n- Paper seems to treat ODEs, but some numerical examples are derived for PDEs. Of course, semi-group theory for PDEs links with exponential integrators, but this seems to go beyond this paper. ", "limitations": "Methodology applies to a particular class of ODES", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The paper reasonably well introduces probabilistic numerics and studies this in the context of ODEs and tackles important notion of stability of numerical approximation. ", "weaknesses": "The paper transfers well-understood properties of numerical approximations of ODEs using the formalism of probabilistic numerics. In my view, this is not a sufficiently novel contribution. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We first want to thank the reviewer for their comments, feedback, and questions.\nIn the following we will address each open point separately.\n\n> Calling the Ornstein-Uhlenbeck process a Gauss-Markov in eqn 3 is overkill.\n\nProbabilistic numerical ODE solvers have in general been established with any\nGauss--Markov prior that allows for easy access of derivatives [1]. Since\nsection 2 is the background section of the paper, we introduced ODE filters in\ntheir generality.\n\n> Notation in eqn 5 of normal distribution not clear what. Why three parameters?\n\nThe first entry of `N(x; m, C)` is not a parameter, but the variable that is\ndescribed by the Gaussian. This is just a slighly more thorough notation than\nthe equivalent short-hand notation `N(m, C)`. To prevent confusion, we will make\nsure that the notation is consistent accross the paper.\n\n> Section 2.2 I don't see why \\mathcal I [y] = 0 when y is a solution to ODE.\nThis says that y_t-f(y_t,t)=0 not \\frac{d}{dt} y_t -f(y_t,t)=0\n\n$E_1$ is a selection matrix that selects the entry that corresponds to the first\nderivative, and similarly $E_0$ selects the zeroth derivative. In terms of the\nGauss--Markov state representation $Y$, $\\mathcal{I}$ therefore encodes exactly\nthe ODE. The statement $\\mathcal{I}\\[y\\] \\equiv 0$ comes with a slight abuse of\nnotation. We could also define $\\mathcal{I}$ not via the \"state\" $Y$ as in eq.\n(7), but in the space of the $d$-dimensional function $y$, as\n$\\mathcal{I}\\[y\\](t) = D y(t) - f(y(t), t)$,\nwhere $D$ is the derivative operator. On $Y^{(0)}$ this then gives\n$\\mathcal{I}\\[Y^{(0)\\}](t) = D Y^{(0)}(t) - f(Y^{(0)}(t), t) = Y^{(1)}(t) - f(Y^{(0)}(t), t)$.\nThis formula is equivalent to the formulation in the paper (eq. 7) and, since\nthe algorithm only operates on the state $Y$, it also closer to the actual\nimplementation of the method.\n\n> Explain the role of Z in 9b\n\n$Z_n$ formally describe the data, _before_ it is actually observed. Then,\nwhen we add the actual data to the model as defined in eq. (9), we can do\n(approximate) Bayesian inference. This notation with prior and likelihood model\nis standard in Bayesian filtering and smoothing (and Bayesian inference in\ngeneral). In the specific case of probabilistic ODE solvers the data is exactly \nzero everywhere, since we want the ODE to hold on the grid.\n\n> Paper seems to treat ODEs, but some numerical examples are derived for PDEs.\n\nThe method of lines can be used to transform semi-linear PDEs into semi-linear\nODEs, as we did in the experiment section. The semi-linear ODE is then solved\nwith the proposed method. Note that spatial discretizations of PDEs are also a\nvery common example to motivate exponential integrators; see for instance the\nintroduction of [2].\n\n> Methodology applies to a particular class of ODES\n\nThe goal of the paper was indeed to develop methods particularly for stiff\nsemi-linear ODEs. And in the context of numerical ODE solvers, developing\nspecific methods for specific problems is quite common. We still want to\nhighlight Section 3.6 of our paper: The proposed exponential integrator can also\nbe applied to any nonlinear ODE by automatically and continuously linearizing\nthe problem, in a manner that is very similar to classic exponential Rosenbrock\nmethods.\n\n\n[1] Tronarp et al, \"Bayesian ODE solvers: the maximum a posteriori estimate\",\nStatistics and Computing (2021)\n\n[2] Hochbruck et al, \"Exponential Rosenbrock-Type Methods\", SIAM Numerical\nAnalysis (2009)\n"}, {"review_id": "WSSSNNDZde", "reviewer": "Reviewer_vChc", "summary": "This paper proposed probabilistic exponential integrators, which is a new class of probabilistic solvers for stiff semi-linear ODEs. More specifically, the integrated OU process is introduced for a functional prior that directly incorporate the linear part of the dynamics. As a result, the proposed methods can be more stable than their non-exponential probabilistic counterparts. In addition, the authors also provide an extension for general non-linear systems, which can be viewed as a probablistic solver based on exponential Rosenbrock-type methods.", "questions": "1. What is the overhead computational time for computing the matrix exponentials? Is it contained in the runtime?\n2. The proposed method can be more accurate when the linear part is more dominant. Since EK1 also introduces a linear approximation, would it perform better when combined with the IOUP prior?", "limitations": "Yes, they do.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper is overall clearly written and nicely organized. \n2. The IOUP prior is neat. Also, some theoretical results of the proposed probabilistic exponential integrator have been established (e.g., equivalence to the classic exponential trapezoidal rule and L-stability).", "weaknesses": "1. The proposed methods are only for semi-linear ODEs, while other baselines can be applied to more general cases.\n2. The IOUP prior is more expensive than the IWP prior.\n3. The Kalman filter/smoother based solver has cubic scaling in the ODE dimension. \n4. It seems that the advantage of the method is more significant when the step size is large (error is large), and becomes negligible when the step size is small.\n\nNote that 2 and 3 have been adimitted by the authors.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We first want to thank the reviewer for their comments, feedback, and questions.\nIn the following we will address each open point separately.\n\n> Weakness 1. The proposed methods are only for semi-linear ODEs, while other baselines can be applied to more general cases.\n\nThe main \"probabilistic exponential integrator\" proposed in the paper is\nindeed only fully applicable to semi-linear ODEs. But, as also stated\nin the Summary section of your review, the proposed Rosenbrock-type\nmethod of Section 3.6 is applicable to general non-linear ODEs.\n\n> Weakness 4. It seems that the advantage of the method is more significant when the step size is large (error is large), and becomes negligible when the step size is small.\n\nThis is a very good point, and this is precisely why \"stability\" is so important for numerical ODE solvers. If the step size is \"small enough\", even the simplest forward Euler method eventually computes an accurate solution---but the problem is that for extremely small steps the runtime becomes prohibitively large. By using more stable solvers, such as exponential integrators, accurate solutions can be computed with larger step sizes.\n\n> Question 1. What is the overhead computational time for computing the matrix exponentials? Is it contained in the runtime?\n\nThe cost of computing the matrix exponential depends on the drift matrix and its\nsize, and thus on the specific ODE and the order of the prior. For fixed step\nsizes, the proposed probabilistic exponential integrator computes _one_ matrix\nexponential, and then re-uses the computed transition matrices for each step.\nThus, it computes exaclty one matrix exponential more than a solver with IWP\nprior. The Rosenbrock-type method however needs to compute a matrix exponential\n_at each timestep_ as the linear part changes at each step, and is thus\ncomputationally more demanding, but also more stable, as shown in Figured 4\nand 5. The reported runtimes include all parts of the solver and thus also the\ncomputation of the matrix exponential.\n\n\n> Question 2. The proposed method can be more accurate when the linear part is more dominant. Since EK1 also introduces a linear approximation, would it perform better when combined with the IOUP prior?\n\nThis is a very good question. Indeed, the proposed method can also be combined\nwith the exact first-order Taylor linearization (EK1). To keep the number of\nmethods in the experiments and figures reasonably small, we have chosen to\ncombine the probabilistic exponential integrator only with the EKL (thus not\nrequiring any additional linearization during the solve), and the\nRosenbrock-type method with the EK1 (as both the Rosenbrock-type method and the\nEK1 require local linearizations anyways); but both solvers can be combined with\nthe different linearization strategies mentioned in Section 3.3. In experiments\nnot included in the paper, we have also tested the non-Rosenbrock method with\nthe EK1, and we have observed that this can be more accurate than the EKL---but\nwith the additional cost/downside of requiring local linearizations.\n"}, {"review_id": "QAXrZEoE4P", "reviewer": "Reviewer_V9iy", "summary": "This paper presents probabilistic exponential integrators as a solver for stiff semi-linear ODEs.Their solver can also be extended to solve general non-linear ODEs with iterative re-linearization. The proposed method is shown to be L-stable in theory and empirically more stable than existing methods.", "questions": "See Weakness", "limitations": "See Weakness", "rating": 7, "confidence": 1, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper is well written. It tackles the problem of solving stiff systems, which is a great challenge in probabilistic numerics. The authors exploit the properties in semi-linear ODEs and develop a probabilistic solver that integrates the fast, linear dynamics into the prior model. The solver is also extended to general nonlinear ODEs via iterative linearization. Both theoretical and empirical analysis are sound.", "weaknesses": "The main difficulty of conventional of stiff ODEs solvers is its sensitivity to step size: small step size is required for stability. Proposition 3 shows the proposed solver is L-stable. However, the experiments mainly focus on final errors across various step sizes. A detailed investigation of stability (boundedness of solutions?) behaviors of solvers under large step size may help readers better understand the proposed approach. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We first want to thank the reviewer for their insightful comments and feedback.\n\n> the experiments mainly focus on final errors across various step\nsizes. A detailed investigation of stability (boundedness of solutions?)\nbehaviors of solvers under large step size may help readers better understand\nthe proposed approach.\n\nThank you for bringing this up. The main advantage of probabilistic exponential\nintegrators as introduced in the paper is indeed their improved stability. A\nmore detailed investigation purely focused on stability could have been helpful\nto make this message more clear; but we want to highlight that the improved\nperformance for large step sizes is demonstrated in the experiments: In both\nFigure 4 and 5, the \"error vs. step size\" plots show that as step sizes\nincrease, the non-exponential methods start to fail (either by diverging or\nachieving very large errors >>1) whereas the exponential methods still provide\nmeaningful solutions; see for example step size $10^{-1}$ in both figures. This\ndemonstrates the improved stability of the exponential integrators.\n"}, {"review_id": "vZyljVNngl", "reviewer": "Reviewer_Lu6G", "summary": "This paper proposes a probabilistic integrator for solving stiff semi-linear ODEs that decouples linear and nonlinear components for improved accuracy.  The approach uses an integrated Ornstein-Uhlenbeck process (IOUP) prior and posterior inference extends classical exponential integrator to exactly solve the linear components.  Nonlinearities are approximated via local linearization similar to the extended Kalman filter (EKF).  The authors present L-stability results and address computation through proposed approximations of the Jacobian that capture only the linear dynamics.", "questions": "* The authors suggest that probabilistic integrators may be preferable for expensive-to-evaluate functions (L259).  How?  Can the method avoid explicit evaluation via prediction?\n\n* The EKF can be quite unreliable in even simple nonlinear dynamical systems, leading to preferred approaches such as the UKF.  Would not such an approach be preferable here to avoid the local linear approximations?  What is the main challenge in a more accurate approximation?", "limitations": "Limitations are clearly described in Sec. 5", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "This is a solid paper; well-written, well-motivated, nice execution.  The results aren’t groundbreaking when compared to baseline, but there are clear benefits of the proposed method, particularly when nonlinearities are less significant.  Furthermore, results are presented honestly and limitations are clearly identified.", "weaknesses": "The paper fails to highlight any benefit of treating the initial value problem (IVP) as a Bayesian inference task.  The posterior distribution *should* characterize uncertainty over the solution, and thus it is possible to quantify and report confidence in the solution.  It is further possible to scrutinize the appropriateness of the IOUP as a prior.  The paper does not address this, but instead prioritizes numerical accuracy and execution time.  While these metrics are certainly important, they are not the primary motivation of probabilistic numerics.  This is not a critical flaw, but perhaps a missed opportunity to demonstrate the benefit of a Bayesian approach to a numerical method.\n\nUnsurprisingly, the proposed method is fairly sensitive to nonlinearity of the ODEs.  The experiment in Sec. 4.1 demonstrates that accuracy degrades as the quadratic term becomes more significant.  In fact, the proposed method seems to be more sensitive to this than baseline IWP methods, albeit with uniformly lower error.  For the simple ODE in Sec. 4.1 all methods achieve very low error (on the order of $10^{-10}$) so the observed improvement with the IOUP is marginal unless the ODE is highly linear.\n\nSome detailed comments:\n  * Fig. 5 should reference subfigures (e.g. (a) is never explicitly mentioned).\n  * L268: Typo \"Since th problem...\"", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We first want to thank the reviewer for their insightful comments, feedback, and\nquestions. In the following we will address each point separately.\n\n\n> The paper fails to highlight any benefit of treating the initial value problem\n(IVP) as a Bayesian inference task. The posterior distribution should\ncharacterize uncertainty over the solution, and thus it is possible to quantify\nand report confidence in the solution. It is further possible to scrutinize the\nappropriateness of the IOUP as a prior. The paper does not address this, but\ninstead prioritizes numerical accuracy and execution time. While these metrics\nare certainly important, they are not the primary motivation of probabilistic\nnumerics. This is not a critical flaw, but perhaps a missed opportunity to\ndemonstrate the benefit of a Bayesian approach to a numerical method.\n\n\nThe reviewer brings up a fair point. However, we would like to point out the\npaucity of papers critically examining different priors. In fact, to our\nknowledge, this is the first to give a clear argument for using anything other\nthan the canonical IWP prior. With this in mind, we hope the present\ncontribution can spark interest in taking the prior selection more seriously in\nupcoming investigations by the PN community.\n\n\n> Unsurprisingly, the proposed method is fairly sensitive to nonlinearity of the\nODEs. The experiment in Sec. 4.1 demonstrates that accuracy degrades as the\nquadratic term becomes more significant. In fact, the proposed method seems to\nbe more sensitive to this than baseline IWP methods, albeit with uniformly lower\nerror.\n\nThe performance benefits of the probabilistic exponential integrator with IOUP\nprior over a solver with IWP prior depend compeletely on the known linear part\nopf the ODE. If the linear part is set to zero, they both coincide and the\nIOUP-based solver performs exactly as the solver with IWP prior. This behaviour\nis exactly what the experiment in Section 4.1 and Figure 3 aim to visualize,\nshowing that as the linear part increasingly dominates the dynamics, the benefit\nof using an IOUP prior increases. Note also that in this experiment the IOUP\nprior never performs worse than the EKL&IWP combination, so we see this not as\n\"sensitivity to nonlinearity\" but rather as \"leveraging the semi-linearity\".\n\n\n> The authors suggest that probabilistic integrators may be preferable for\nexpensive-to-evaluate functions (L259). How? Can the method avoid explicit\nevaluation via prediction?\n\nThe comment in line 259 refers to choosing between the standard and Rosenbrock\nversion of the proposed probabilistic exponential integrator, and the suggestion\nstems purely from the required number of function evaluations per step: Since\nthe local linearization of the Rosenbrock method requires evaluating the vector\nfield, the standard probabilistic exponential integrator requires less\nevaluations and thus might be preferrable when the vector field becomes\nexceedingly expensive to evaluate. But note that both methods perform at least\none evaluation of the vector field per step.\n\n\n> The EKF can be quite unreliable in even simple nonlinear dynamical systems,\nleading to preferred approaches such as the UKF. Would not such an approach be\npreferable here to avoid the local linear approximations? What is the main\nchallenge in a more accurate approximation?\n\nThis is a very good point. It is indeed well known that the UKF or other\nfiltering/smoothing methods can be advantageous over the EKF/EKS in certain\nstate estimation problems. In the specific context of probabilistic ODE solvers,\nthe UKF has been previously suggested [1,2], but a more extensive evaluation of\nits properties and utilities would certainly be interesting in order to better\nunderstand for which types of problems it would be most beneficial.\n\nNote also that in the context of probabilistic ODE solvers, a more accurate\nlinearization not the only way to improve the accuracy of the method: we could\nalso just select a smaller step size. This is not the case for most standard\nfiltering problems. And in the regime of \"very small\" steps, a local Taylor\napproximation becomes increasingly accurate.\n\n\n[1] Kersting et al, \"Active uncertainty calibration in Bayesian ODE solvers\",\nUAI (2016)\n\n[2] Tronarp et al, \"Probabilistic solutions to ordinary differential equations\nas nonlinear Bayesian filtering: a new perspective\", Statistics and Computing\n(2019)\n"}]

Author rebuttal:

Meta-review: This is a well-written contribution that narrows the performance gap between classical and probabilistic numerical methods for ODEs and I am delighted to agree with the Reviewers that it should be accepted.
Decision: Accept (poster)