[subsection]section *[subsubsection]subsection

# Numerical solution of kinetic SPDEs via stochastic Magnus expansion

Kevin Kamm Dipartimento di Matematica, Università di Bologna, Bologna, Italy.
e-mail<EMAIL_ADDRESS>Stefano Pagliarani Dipartimento di Matematica,
Università di Bologna, Bologna, Italy. e-mail<EMAIL_ADDRESS>Andrea Pascucci Dipartimento di Matematica, Università di Bologna, Bologna,
Italy. e-mail<EMAIL_ADDRESS>
###### Abstract

In this paper, we show how the Itô-stochastic Magnus expansion can be used to
efficiently solve stochastic partial differential equations (SPDE) with two
space variables numerically. To this end, we will first discretize the SPDE in
space only by utilizing finite difference methods and vectorize the resulting
equation exploiting its sparsity.

As a benchmark, we will apply it to the case of the stochastic Langevin
equation with constant coefficients, where an explicit solution is available,
and compare the Magnus scheme with the Euler-Maruyama scheme. We will see that
the Magnus expansion is superior in terms of both accuracy and especially
computational time by using a single GPU and verify it in a variable
coefficient case. Notably, we will see speed-ups of order ranging form 20 to
200 compared to the Euler-Maruyama scheme, depending on the accuracy target
and the spatial resolution.

Keywords: Magnus Expansion, Stochastic Langevin Equation, Numerical Solutions
for SPDE, GPU Computing.
Acknowledgements: This project has received funding from the European Union’s
Horizon 2020 research and innovation programme under the Marie Sklodowska-
Curie grant agreement No 813261 and is part of the ABC-EU-XVA project.
Data Availability: All data generated or analysed during this study are
included in this published article. In particular, the code to produce the
numerical experiments is available at
https://github.com/kevinkamm/MagnusSPDE2D.

##  1 Introduction

In the recent paper [10] we have introduced a stochastic version of the Magnus
expansion [12] (hereafter abbreviated in ME), namely an exponential
representation for solutions of matrix-valued stochastic differential
equations (SDEs), in the Itô sense, of the form

$dX_{t}=B_{t}X_{t}dt+A_{t}X_{t}dW_{t},\qquad X_{0}=I_{d},$ (1.1)

for non-commuting bounded progressively measurable random matrices $A_{t}$ and
$B_{t}$. These results have recently found applications in the study of so-
called signature cumulants [7], semi-linear Itô SDEs [17], linear SDEs on
matrix Lie groups [13], stochastic modeling of motion in turbulent flows [3],
stability of multi-variate geometric Brownian motion [2] and modeling rating
transition matrices in quantitative finance [9].

In this note we investigate the application of the ME to the numerical
resolution of a class of hypoelliptic stochastic partial differential
equations (SPDEs) that naturally arises in physics and mathematical finance.
The deterministic prototype of such SPDEs is the classical Langevin equation

$\frac{1}{2}\partial_{vv}u_{t}+v\partial_{x}u-\partial_{t}u_{t}=0,$ (1.2)

where the variables $t\geq 0$, $x\in\mathbb{R}$ and $v\in\mathbb{R}$
respectively stand for time, position and velocity, and the unknown
$u_{t}=u_{t}(x,v)$ stands for the density of a particle in the phase space.
Notice that (1.2) is a degenerate, non-uniformly parabolic PDE. Perturbations
of (1.2) with variable coefficients appear in linear and non-linear form in
several applications in kinetic theory (see, for instance, [11], [5] and [4]);
also, (1.2) describes path-dependent financial derivatives such as Asian
options and volatility contracts (see, for instance, [14] and [6]).

We consider here the stochastic version of (1.2), which is the kinetic SPDE

$\displaystyle
du_{t}=\left(\frac{a_{t}}{2}\partial_{vv}+v\partial_{x}+b_{t}\partial_{v}+c_{t}\right)u_{t}dt+\left(\sigma_{t}\partial_{v}+\beta_{t}\right)u_{t}dW_{t},$
(1.3)

where $a_{t},b_{t},c_{t},\sigma_{t}$ and $\beta_{t}$ are non-constant (i.e.,
for example, $a_{t}=a_{t}(x,v)$) and possibly random coefficients. Here $W$
denotes a Wiener process defined on a complete probability space
$(\Omega,\mathcal{F},\mathbb{P})$ endowed with a filtration
$\left(\mathcal{F}_{t}\right)_{t\geq 0}$ satisfying the usual conditions. SPDE
(1.3) naturally appears in stochastic filtering theory: as shown in [15, 16],
the fundamental solution of (1.3) is the conditional transition density of a
two-dimensional stochastic process representing the position and the velocity
of a particle under partial observation.

The numerical solution of (1.3) is a challenging issue, as standard
techniques, such as Euler methods, can be cumbersome and excessively time-
consuming. In this paper, we want to demonstrate how the Itô-stochastic ME can
be used as an efficient numerical tool to solve SPDEs meanwhile exploiting
modern computer architectures such as GPUs and multiple CPUs. The theoretical
groundwork for the Itô-stochastic ME was established in [10] alongside some
numerical experiments using a GPU to fully utilize the parallel-in-time and
parallel-in-simulation features of the ME. Here, we show how the Magnus
expansion can be used iteratively to overcome the constraint of a small
convergence radius in time. Indeed, by [10, Theorem 1], the representation
$X_{t}=e^{Y_{t}}$ for the unique strong solution to (1.1) is valid only up to
a strictly positive stopping time. More precisely, we have the following:

######  Theorem 1.1 ([10]).

Let $A_{t}$ and $B_{t}$ be bounded progressively measurable matrices in
$\mathbb{R}^{d\times d}$ and let
$(\Omega,\mathcal{F},\mathbb{P},(\mathcal{F}_{t})_{t\geq 0})$ be a filtered
probability space equipped with a standard Brownian motion $W$. For $T>0$ let
also $X=(X_{t})_{t\in[0,T]}$ be the unique strong solution to (1.1). There
exists a strictly positive stopping time $\tau\leq T$ such that:

  1. (i)

$X_{t}$ has a real logarithm $Y_{t}\in\mathbb{R}^{d\times d}$ up to time
$\tau$, i.e.

$X_{t}=e^{Y_{t}},\qquad 0\leq t<\tau;$ (1.4)

  2. (ii)

the following representation holds $\mathbb{P}$-almost surely:

$Y_{t}=\sum_{n=1}^{\infty}Y^{(n)}_{t},\qquad 0\leq t<\tau,$ (1.5)

where $Y^{(n)}$ is the $n$-th term in the stochastic ME (see also formulas
(1.7) below, in the case of constant matrices $A$ and $B$);

  3. (iii)

there exists a positive constant $C$, only dependent $A$, $B$, $T$ and $d$,
such that

$\mathbb{P}(\tau\leq t)\leq Ct,\qquad t\in[0,T].$ (1.6)

The first point of Theorem 1.1 tells us that the ME only converges up to a
stopping time: to overcome this restriction, the numerical implementation of
the ME requires to apply it iteratively in time. Clearly, by (1.6) the
convergence of the ME is problem-dependent, meaning that there is no universal
best time step-size for the ME. Point (ii) of Theorem 1.1 actually yields the
numerical scheme by truncating the infinite series (1.5): we will see that in
practice it is sufficient to consider only two or three terms to obtain a good
degree of accuracy. For this reason, we will not recall the entire general
expression for the terms $Y_{t}^{(n)}$, but rather provide a user guide in
Appendix A on how to easily derive the following expansion formulas, up to
order $3$, in the case of constant $A$ and $B$:

$\begin{split}Y_{t}^{1}&=Bt+AW_{t},\\\
Y_{t}^{2}&=Y_{t}^{1}-\frac{1}{2}A^{2}t+\left[B,A\right]\int_{0}^{t}{W_{s}ds}-\frac{1}{2}\left[B,A\right]tW_{t},\\\
Y_{t}^{3}&=Y_{t}^{2}+\left[\left[B,A\right],A\right]\left(\frac{1}{2}\int_{0}^{t}{W_{s}^{2}ds}-\frac{1}{2}W_{t}\int_{0}^{t}{W_{s}ds}+\frac{1}{12}tW_{t}^{2}\right)\\\
&\quad+\left[\left[B,A\right],B\right]\left(\int_{0}^{t}{sW_{s}ds}-\frac{1}{2}t\int_{0}^{t}{W_{s}ds}-\frac{1}{12}t^{2}W_{t}\right),\end{split}$
(1.7)

where $\left[A,B\right]\coloneqq AB-BA$.

The paper is organized as follows. In Section 2 we show how to discretize a
parabolic-type SPDE by an SDE of the type (1.1). In Section 3 we show the
numerical experiments in the special case of the stochastic Langevin equation
with constant coefficients where the exact solution is available in closed
form. First, the implementation and notations for the error-analysis are
introduced and afterwards several tests concerning the parameters of the ME
are discussed. Then we compare the performance of the iterated ME with that of
standard Euler-Maruyama schemes. In Section 3.2 we consider the more general
case of SPDE (1.3) with variable coefficients. In Appendix A we provide a
heuristic derivation of the expansion formulas up to order three.

##  2 SPDEs with two space-variables

In the introduction we discussed the stochastic Langevin equation as our main
application for the ME.

In this section, we show how to derive a numerical scheme for a general
parabolic SPDE by combining space-discretization and ME. We will perform some
formal computations, which hold even for a fairly general class of SPDEs.
Since the computations are understood in a formal manner, we do not impose any
further conditions on the coefficients of the following type of SPDE

$\displaystyle\left\\{\begin{aligned} &\begin{aligned} \hskip
10.00002pt&du_{t}(x,v)=\biggl{(}h(x,v)u_{t}(x,v)+f^{x}(x,v)\partial_{x}u_{t}(x,v)+f^{v}(x,v)\partial_{v}u_{t}(x,v)\\\
&\quad+\frac{1}{2}g^{xx}(x,v)\partial_{xx}u_{t}(x,v)+g^{xv}(x,v)\partial_{xv}u_{t}(x,v)+\frac{1}{2}g^{vv}(x,v)\partial_{vv}u_{t}(x,v)\biggr{)}dt\\\
&\quad+\big{(}\sigma(x,v)u_{t}(x,v)+\sigma^{x}(x,v)\partial_{x}u_{t}(x,v)+\sigma^{v}(x,v)\partial_{v}u_{t}(x,v)\big{)}dW_{t}\end{aligned}\\\
&u_{0}(x,v)=\phi(x,v).\end{aligned}\right.$ (2.1)

For simplicity, we will only consider the case of time-independent and non-
random coefficients but it is straightforward to include them in the formal
computations.

In the following subsections, we will demonstrate how to derive an SDE of the
form (1.1) to approximate the SPDE (2.1) and apply the Magnus expansion to it.
In the end, we will recall the Euler-Maruyama scheme for the approximating
SDE.

#### Space discretization and Magnus expansion.

Following [10, Section 3.2], we will discretize the space variables but not
time. After that we will vectorize the equation to derive a matrix-valued
equation. We introduce the following two homogeneous grids for position and
velocity respectively: $\mathbb{X}^{n_{x}}_{a_{x},b_{x}}$ with $n_{x}+2$
points on the subset $[a_{x},b_{x}]\subset\mathbb{R}$,
$\mathbb{V}^{n_{v}}_{a_{v},b_{v}}$ with $n_{v}+2$ points on the subset
$[a_{v},b_{v}]\subset\mathbb{R}$:

$\displaystyle\mathbb{X}^{n_{x}}_{a_{x},b_{x}}$
$\displaystyle\coloneqq\left\\{x_{i}^{n_{x}}\in[a_{v},b_{v}]:x_{i}^{n_{v}}=a_{x}+i\Delta
x,\ i=0,\dots,n_{v}+1\right\\},$ $\displaystyle\Delta x$
$\displaystyle\coloneqq\frac{b_{x}-a_{x}}{n_{x}+1},$ (2.2)
$\displaystyle\mathbb{V}^{n_{v}}_{a_{v},b_{v}}$
$\displaystyle\coloneqq\left\\{v_{j}^{n_{v}}\in[a_{v},b_{v}]:v_{j}^{n_{v}}=a_{v}+j\Delta
v,\ j=0,\dots,n_{v}+1\right\\},$ $\displaystyle\Delta v$
$\displaystyle\coloneqq\frac{b_{v}-a_{v}}{n_{v}+1}.$ (2.3)

Let us denote by $\mathrm{vec}$ the isomorphism of transforming a matrix to a
larger column-vector by stacking each column in the matrix below each other,
i.e.

$\displaystyle\mathrm{vec}:\mathbb{R}^{n_{x}\times n_{v}}$
$\displaystyle\rightarrow\mathbb{R}^{n_{x}\cdot n_{v}\times 1},$
$\displaystyle A=[a_{ij}]$
$\displaystyle\mapsto\mathrm{vec}(A)\coloneqq[a_{1,1},\dots,a_{n_{x},1},a_{1,2},\dots,a_{n_{x},2},\dots,a_{1,n_{v}},\dots,a_{n_{x},n_{v}}]^{T}.$

Moreover, let us define

$\displaystyle
u_{t}^{n_{x},n_{v}}\coloneqq\left(u_{t}(x_{i},v_{j})\right)_{\begin{subarray}{c}i=1,\dots,n_{x}\\\
j=1,\dots,n_{v}\end{subarray}},$
$\displaystyle\phi^{n_{x},n_{v}}\coloneqq\left(\phi(x_{i},v_{j})\right)_{\begin{subarray}{c}i=1,\dots,n_{x}\\\
j=1,\dots,n_{v}\end{subarray}},$
$\displaystyle\Phi^{n_{x}n_{v}}\coloneqq\mathrm{vec}\left(\phi^{n_{x},n_{v}}\right).$

We consider the $n_{x}n_{v}$-dimensional SDE

$\displaystyle
dU^{n_{x}n_{v}}_{t}=BU^{n_{x}n_{v}}_{t}dt+AU^{n_{x}n_{v}}_{t}dW_{t},\quad
U^{n_{x}n_{v}}_{0}=\Phi^{n_{x}n_{v}},$ (2.4)

where $A$ and $B$ are $(n_{x}n_{v}\times n_{x}n_{v})$-matrices, which will be
defined via space-finite differences in the next paragraph in a way that
$U^{n_{x}n_{v}}_{t}\approx\mathrm{vec}\left(u_{t}^{n_{x},n_{v}}\right)$ with
respect to a suitable norm. This equation can be solved by first computing its
fundamental solution, then by multiplying it with the initial datum, i.e.

$\displaystyle dX_{t}$ $\displaystyle=BX_{t}dt+AX_{t}dW_{t},\quad
X_{0}=I\in\mathbb{R}^{n_{x}n_{v}\times n_{x}n_{v}},$ $\displaystyle
U^{n_{x}n_{v}}_{t}$ $\displaystyle=X_{t}\Phi^{n_{x}n_{v}}.$

The fundamental solution now can be approximated with the Magnus expansion,
i.e. $X_{t}\approx\exp\left(Y_{t}\right)$. Since $A$ and $B$ will be very
large in this case, we will utilize the sparsity of $A$ and $B$, as well as
using a special algorithm specifically designed to compute the matrix-
exponential times a vector denoted by `expmvtay2`, which does not need to
compute the whole matrix-exponential first. This is crucial for the
implementation and explained in further detail in [8]. The approximation
formulas for $X_{t}$, $Y_{t}$ and their derivation are deferred to Appendix A.

#### Computation of $A$ and $B$.

The idea is to discretize the first and second-order derivatives we will use
central differences with zero-boundary conditions. Therefore, let us introduce
the following matrices corresponding to the finite differences

$\displaystyle D^{x}$ $\displaystyle\coloneqq\frac{1}{2\Delta
x}\mathrm{tridiag}^{n_{x},n_{x}}\left(-1,0,1\right),$ $\displaystyle D^{v}$
$\displaystyle\coloneqq\frac{1}{2\Delta
v}\mathrm{tridiag}^{n_{v},n_{v}}\left(-1,0,1\right),$ $\displaystyle D^{xx}$
$\displaystyle\coloneqq\frac{1}{\left(\Delta
x\right)^{2}}\mathrm{tridiag}^{n_{x},n_{x}}\left(1,-2,1\right),$
$\displaystyle D^{vv}$ $\displaystyle\coloneqq\frac{1}{\left(\Delta
v\right)^{2}}\mathrm{tridiag}^{n_{v},n_{v}}\left(1,-2,1\right).$

Additionally, we will introduce the following matrices corresponding to the
coefficient functions on the discretized spatial grid

$\displaystyle Z^{w}$
$\displaystyle\coloneqq\left(z^{w}(x_{i},v_{j})\right)_{\begin{subarray}{c}i=1,\dots,n_{x}\\\
j=1,\dots,n_{v}\end{subarray}},$ $\displaystyle\Sigma^{w}$
$\displaystyle\coloneqq\left(\sigma^{w}(x_{i},v_{j})\right)_{\begin{subarray}{c}i=1,\dots,n_{x}\\\
j=1,\dots,n_{v}\end{subarray}}$

for $Z=F,G,H$, $z=f,g,h$, respectively, and
$w\in\left\\{x,v,xx,xv,vv\right\\}$.

Let us start with discretizing the first-order derivative with respect to $x$.
First, we replace the partial derivative by the first-order central
differences and assume zero-boundary conditions, leading to

$\displaystyle f^{x}(x_{i},v_{j})\partial_{x}u_{t}(x_{i},v_{j})$
$\displaystyle\approx
f^{x}(x_{i},v_{j})\frac{u_{t}(x_{i+1},v_{j})-u_{t}(x_{i-1},v_{j})}{2\Delta x}$

for all $i=1,\dots,n_{x}$ and $j=1,\dots,n_{v}$. As aforementioned, we need to
extract the correct coefficient matrix for the vectorized equation (2.4).

In our notations, a derivative in $x$ is a multiplication of the corresponding
finite-difference matrix from the left to $u_{t}^{n_{x},n_{v}}$, i.e.

$\displaystyle\left(\frac{u_{t}(x_{i+1},v_{j})-u_{t}(x_{i-1},v_{j})}{2\Delta
x}\right)_{\begin{subarray}{c}i=1,\dots,n_{x}\\\
j=1,\dots,n_{v}\end{subarray}}=D^{x}u_{t}^{n_{x},n_{v}}.$

Using the Kronecker product, it is well-known for compatible matrices
$D_{1}UD_{2}=C$ that

$\displaystyle\mathrm{vec}\left(C\right)=\mathrm{vec}\left(D_{1}UD_{2}\right)=\left(D_{2}^{T}\otimes
D_{1}\right)\mathrm{vec}\left(U\right).$

In our case, this leads to

$\displaystyle\mathrm{vec}\left(D^{x}u_{t}^{n_{x},n_{v}}\right)=\left(I_{n_{v}}\otimes
D^{x}\right)U_{t}^{n_{x}n_{v}}.$

Now, we need to deal with the coefficients as well. Denoting by $\odot$ the
Hadamard, or elementwise, product, it is easy to see that

$\displaystyle\mathrm{vec}\left(F^{x}\odot\left(D^{x}u_{t}^{n_{x},n_{v}}\right)\right)=\mathrm{vec}\left(F^{x}\right)\odot\mathrm{vec}\left(\left(D^{x}u_{t}^{n_{x},n_{v}}\right)\right)=\mathrm{diag}\left(\mathrm{vec}\left(F^{x}\right)\right)\mathrm{vec}\left(\left(D^{x}u_{t}^{n_{x},n_{v}}\right)\right).$

Using these two observations together yields

$\displaystyle\left(f^{x}(x_{i},v_{j})\partial_{x}u_{t}(x_{i},v_{j})\right)_{\begin{subarray}{c}i=1,\dots,n_{x}\\\
j=1,\dots,n_{v}\end{subarray}}\approx\mathrm{diag}\left(\mathrm{vec}\left(F^{x}\right)\right)\left(I_{n_{v}}\otimes
D^{x}\right)U_{t}^{n_{x}n_{v}}.$

This reasoning holds true for all other derivatives as well, i.e. an operation
in $x$ is a matrix multiplication from the left and in $v$ it is the matrix
multiplication from the right with the transposed matrix.

Conclusively, we have

$\displaystyle B$
$\displaystyle\coloneqq\mathrm{diag}\left(\mathrm{vec}\left(H\right)\right)+\mathrm{diag}\left(\mathrm{vec}\left(F^{x}\right)\right)\left(I_{n_{v}}\otimes
D^{x}\right)+\mathrm{diag}\left(\mathrm{vec}\left(F^{v}\right)\right)\left(D^{v}\otimes
I_{n_{x}}\right)$
$\displaystyle\quad+\frac{1}{2}\mathrm{diag}\left(\mathrm{vec}\left(G^{xx}\right)\right)\left(I_{n_{v}}\otimes
D^{xx}\right)+\mathrm{diag}\left(\mathrm{vec}\left(G^{xv}\right)\right)\left(D^{v}\otimes
D^{x}\right)$
$\displaystyle\quad+\frac{1}{2}\mathrm{diag}\left(\mathrm{vec}\left(G^{vv}\right)\right)\left(D^{vv}\otimes
I_{n_{x}}\right),$ $\displaystyle A$
$\displaystyle\coloneqq\mathrm{diag}\left(\mathrm{vec}\left(\Sigma\right)\right)+\mathrm{diag}\left(\mathrm{vec}\left(\Sigma^{x}\right)\right)\left(I_{n_{v}}\otimes
D^{x}\right)+\mathrm{diag}\left(\mathrm{vec}\left(\Sigma^{v}\right)\right)\left(D^{v}\otimes
I_{n_{x}}\right).$

Figure 1: Sparsity patterns of $A$, $B$ and commutators of the SPDE from
Section 3.1 with coefficients (3.4) and $d=50$. nz stands for the number of
non-zero entries.

######  Remark 2.1.

We can see that the sparsity of $B$ and $A$ is mostly determined by the
finite-differences matrices due to the elementwise product of the coefficient
functions. The coefficient functions can further reduce or increase the
sparsity by zero entries or cancellations. Additionally, this implies that the
number of non-zero diagonals is independent of the size of $B$ and $A$
respectively even though the density of non-zero elements will decline for
increasing dimensions.

An illustration of the sparsity pattern can be seen in Figure 1: In the upper
left corner is the pattern for $A$, followed by $B$, $\left[B,A\right]$,
$\left[\left[B,A\right],A\right]$ and $\left[\left[B,A\right],B\right]$. The
blue lines represent non-zero entries. The pictures are ordered by the number
of non-zero diagonals, i.e. a diagonal with at least one non-zero entry, which
are 2, 5, 5, 8 and 10 respectively. We can see that the sparsity decreases
with the order of the commutator. In Section 3.2 we will see that one
consequence of the reduced sparsity is a decrease in computational efficiency
of the Magnus expansion.

#### Iterated Magnus expansion

Figure 2: Schematic of the iterative Magnus scheme

An advantage of the Magnus expansion over pure time-iterative methods is that,
as long as it is convergent, it does not need a previous time-iteration to
compute the next time step. However, as we discussed in the introduction, the
convergence time-interval can be small for the Magnus expansion. Therefore, if
$T>0$ is larger than the convergence radius, it is necessary to split the
initial time interval $[0,T]$ into smaller sub-intervals
$(0,t_{1}],(t_{1},t_{2}],\dots,(t_{n},T]$. As illustrated in Figure 2, we will
evaluate the Magnus expansion consecutively on each of them, i.e. use the
terminal evaluation as the initial point of the next sub-interval.

We will call this method _iterated Magnus expansion_. On each sub-interval,
the Magnus expansion still has the usual parallel-in-time features.
Furthermore, due to its relatively large convergence region in time, we can
use fewer iterations compared to other iterative methods, e.g. for the Euler-
Maruyama scheme.

#### Euler-Maruyama

In this case, we do not need to vectorize the equation but we have to
discretize the time-derivative as well. With the same notation and reasoning
from above we obtain

$\displaystyle u_{t_{k+1}}^{n_{x},n_{v}}$ $\displaystyle\approx
u_{t_{k}}^{n_{x},n_{v}}+\Biggl{(}H\odot
u_{t_{k}}^{n_{x},n_{v}}+F^{x}\odot\left(D^{x}u_{t_{k}}^{n_{x},n_{v}}\right)+F^{v}\odot\left(u_{t_{k}}^{n_{x},n_{v}}\left(D^{v}\right)^{T}\right)$
(2.5)
$\displaystyle\qquad+\frac{1}{2}G^{xx}\odot\left(D^{xx}u_{t_{k}}^{n_{x},n_{v}}\right)+G^{xv}\odot\left(D^{x}u_{t_{k}}^{n_{x},n_{v}}\left(D^{v}\right)^{T}\right)+\frac{1}{2}G^{vv}\odot\left(u_{t_{k}}^{n_{x},n_{v}}\left(D^{vv}\right)^{T}\right)\Biggr{)}\Delta
t$ (2.6) $\displaystyle\quad+\Biggl{(}\Sigma\odot
u_{t_{k}}^{n_{x},n_{v}}+\Sigma^{x}\odot\left(D^{x}u_{t_{k}}^{n_{x},n_{v}}\right)+\Sigma^{v}\odot\left(u_{t_{k}}^{n_{x},n_{v}}\left(D^{v}\right)^{T}\right)\Biggr{)}\Delta
W_{t_{k}},$ (2.7)

where $\Delta t:=t_{k+1}-t_{k}>0$ for any $k$ and $\Delta
W_{t_{k}}:=W_{t_{k+1}}-W_{t_{k}}$. In the case of separable coefficients the
Hadamard product can be replaced by matrix products with diagonal matrices
from the left and right.

##  3 Numerical experiments

We present here some numerical tests that demonstrate how the Magnus expansion
can be efficiently applied to approximate the solutions of stochastic partial
differential equations (SPDEs) with two spatial variables.

In this paper, we will focus on the stochastic Langevin equation. First, we
consider the case of constant coefficients. This enables us to use an exact
solution as a benchmark. We will perform experiments to test the accuracy of
the approximate solutions, elaborate on the computational times and discuss
effects of the spatial boundaries. Afterwards, we will consider the stochastic
Langevin equation with non-constant coefficients.

To this end, we will first elaborate on the implementation of the code and
introduce some notation and norms for the error-analysis.

#### Implementation

As mentioned beforehand, we implement the algorithm iteratively as illustrated
in Figure 2. Now, we explain how to implement the algorithm for one time
interval in Matlab 2022a. For this, let us have a close look at (1.7). We can
see that, in this particular case, the computations of the Lebesgue-integrals
are decoupled from the coefficients $A$ and $B$. This makes it possible to
evaluate the Lebesgue-integrals for all trajectories in parallel by using
vectorization. In particular, we chose simple Riemann-sums to compute the
integrals.

Now, we need to compute the commutators of the scheme. We found that for
larger problems, e.g. $d\geq 100$ it is faster to copy $A$ and $B$ to the
device and compute the commutators there. As it turns out, Matlab 2022a as of
now does not support three-dimensional sparse matrices. Therefore, we decided
to loop over all desired time evaluations (for our experiments only the
terminal time) within the current sub-interval and over all simulations.
Since, the Magnus expansion is parallel-in-time and parallel-in-simulations in
each iteration, we use a threading environment to keep a single GPU busy with
computing the expansion formulas and afterwards the matrix-vector exponential.
This is also more efficient in terms of memory usage than using vectorization.
As a side-note, it is straightforward to use multiple GPUs.

For $d\leq 100$ we use CPUs only because it is faster than copying to device
in each iteration. For larger increasing spatial resolution the advantage of a
single GPU increases more and more compared to CPUs only.

The choice of algorithm for evaluating the matrix-vector exponentials turned
out to be crucial. We tested `expmv` by [1] and `expmvtay2` by [8], as well as
a Krylov-subspace implementation called `expv`. For our purposes, `expmvtay2`
was more accurate and significantly faster than the other methods. Their
algorithm makes it possible to utilize the sparsity of the problem and is GPU-
applicable without first computing the entire matrix exponential.

We tried to improve the performance of the Euler scheme as much as possible
and tested different ways to compute the matrix multiplications within each
iteration of the scheme and only save the values at the terminal points for a
further speed-up. For the matrix multiplications we tested full and sparse
matrix multiplications on GPU and CPUs. For computations using sparse
matrices, one has to loop over the simulations, since there is no analogue to
`pagemtimes` in the full case. Copying from host to device with our test cases
was more expensive than a pure CPU implementation and the low-level multi-
threading of `pagemtimes` surpassed the benefit of sparsity as well.
Therefore, we use full matrix multiplication using all available CPU cores. On
the machine we use, we could choose between 1 and 64 cores. There was no
benefit after using 12 cores, which we use for all of our experiments.

We fix in our tests the number of simulations to $M=100$. If we increase the
number of simulations, both methods will benefit from more available computer
hardware as expected. However, even for $M=100$, the Magnus method would
benefit immediately from another GPU, since the computations of matrix-vector
exponentials are costly, which is the bottleneck of our implementation.

#### Error and notations

For the numerical error analysis we will make use of the following notations.
We denote by $U^{\text{ref}}$ and by $U^{\text{app}}$ a benchmark and an
approximate solution, respectively.

Henceforth, we will choose
$\mathbb{V}^{n_{v}}_{a_{v},b_{v}}=\mathbb{X}^{n_{x}}_{a_{x},b_{x}}$ symmetric,
centered around zero and set $d\coloneqq n_{x}=n_{v}$ to make our analysis a
bit easier. Moreover, we choose the cut-off region of $\mathbb{R}^{2}$ as
$a_{x}=a_{v}=-4$ and $b_{x}=b_{v}=4$ in all experiments.

Furthermore, we want to study the impact of the zero-boundary condition and
therefore will use only a central part with varying size of the whole solution
matrix at a given time, which is illustrated in Figure 3. To vary the size we
introduce a new parameter $\kappa=0,1,\dots$ indicating the $2^{-\kappa}$-th
central part of the solution matrix $U_{t}^{d}$, which we will consider for
our error analysis and denote the corresponding truncated matrix by
$U_{t}^{\text{ref},d,\kappa},U_{t}^{\text{app},d,\kappa}\in\mathbb{R}^{\left\lfloor\frac{d}{2^{\kappa}}\right\rfloor\times\left\lfloor\frac{d}{2^{\kappa}}\right\rfloor}$.
Also we set

$\displaystyle
I^{\kappa}\coloneqq\left\\{\left\lfloor\frac{d}{2}-\frac{d}{2^{\kappa+1}}\right\rfloor,\dots,\left\lfloor\frac{d}{2}+\frac{d}{2^{\kappa+1}}\right\rfloor\right\\}$

to collect the corresponding indices.

$U_{t}^{d,\kappa}$ $i=\left\lfloor\frac{d}{2}-\frac{d}{2\cdot
2^{\kappa}}\right\rfloor$ $i=\left\lfloor\frac{d}{2}+\frac{d}{2\cdot
2^{\kappa}}\right\rfloor$ $j=\left\lfloor\frac{d}{2}-\frac{d}{2\cdot
2^{\kappa}}\right\rfloor$ $j=\left\lfloor\frac{d}{2}+\frac{d}{2\cdot
2^{\kappa}}\right\rfloor$
$\frac{d}{2^{\kappa}}$$=U_{t}^{d}$$\frac{d}{2^{\kappa}}$ Figure 3: Graphical
representation of $U_{t}^{d,\kappa}$ compared to $U_{t}^{d}$ for the error
analysis to disregard boundary effects.

For our error analysis we will use the following three norms:

$\displaystyle\mathrm{ME}_{t}^{d,\kappa}$
$\displaystyle\coloneqq\frac{1}{M}\sum_{m=1}^{M}{\left|U_{t,m}^{\text{ref},d,\kappa}-U_{t,m}^{\text{app},d,\kappa}\right|}\in\mathbb{R}^{\left\lfloor\frac{d}{2^{\kappa}}\right\rfloor{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\times}\left\lfloor\frac{d}{2^{\kappa}}\right\rfloor}$
(3.1) $\displaystyle\mathrm{AME}_{t}^{d,\kappa}$
$\displaystyle\coloneqq\frac{1}{\left|I^{\kappa}\right|^{2}}\sum_{i,j}{\big{(}\mathrm{ME}_{t}^{d,\kappa}\big{)}_{ij}}\in\mathbb{R}$
(3.2) $\displaystyle\mathrm{Err}_{t}^{d,\kappa}$
$\displaystyle\coloneqq\frac{1}{M}\sum_{m=1}^{M}{\frac{\left\|{U_{t,m}^{\text{ref},d,\kappa}-U_{t,m}^{\text{app},d,\kappa}}\right\|_{F}}{\left\|{U_{t,m}^{\text{ref},d,\kappa}}\right\|_{F}}\in\mathbb{R}},$
(3.3)

where $\left\|{\cdot}\right\|_{F}$ denotes the Frobenius norm. The first one
is a matrix consisting of a mean absolute error between the reference solution
and the approximation for each point in the grid $I^{\kappa}\times I^{\kappa}$
by taking the mean over all trajectories. The second error is the average of
the first error for the corresponding grid indicated by $\kappa$. The larger
$\kappa$ the further away we are from the spatial boundary, as illustrated in
Figure 3. The third error is the mean of a relative error between the
reference solution and the approximation taking all points in the region
corresponding to $\kappa$ into account. It will serve as our main error norm
in this paper.

Regarding the Lebesgue discretization $\Delta^{\mathrm{Leb}}_{t}$ we have
observed that the computational times are more or less constant for
$\Delta^{\mathrm{Leb}}_{t}\in\left[10^{-5},10^{-2}\right)$ and increase for
lower $\Delta^{\mathrm{Leb}}_{t}$ significantly. This can be explained by the
fact that the discretization is only used to compute Lebesgue-integrals of the
form $\int_{0}^{t}{s^{p}W_{s}^{q}ds}$. Therefore, the computation can be
vectorized, which is very fast compared to the computation of the Magnus
logarithm and the matrix-vector exponentials for
$\Delta^{\mathrm{Leb}}_{t}\geq 10^{-5}$.

Since, there seems to be no computational drawback, we suggest to use a
Lebesgue-discretization equal to $10^{-3}$ or $10^{-4}$ and fix it for all
tests to $\Delta^{\mathrm{Leb}}_{t}=10^{-4}$.

We have verified a linear behavior (with a slope less than one until the GPU
is fully saturated) in the number of simulations $M$. In our experiments, we
decided to use $M=100$ simulations and display always the average
computational times for one simulation.

Also the computational effort with respect to the finite time horizon $T$
scales linearly for all methods. Thus, we use $T=1$ as our terminal time.

Table 1: Notations for the numerical experiments. euler | Euler-Maruyama scheme (2.7)
---|---
m1 | Iterated Magnus scheme of order 1
m2 | Iterated Magnus scheme of order 2
m3 | Iterated Magnus scheme of order 3
M | Number of simulations
d | Number of grid points in $\mathbb{X}$ and $\mathbb{V}$
$\Delta_{t}$ | step-size of euler or Magnus
$\Delta^{\mathrm{Leb}}_{t}$ | discretization of the Lebesgue-integrals for Magnus
M2, x | Magnus order 2 with step-size $\Delta_{t}=x$
M3, x | Magnus order 3 with step-size $\Delta_{t}=x$
E, x |  euler with step-size $\Delta_{t}=x$

Throughout the experiments we will use the notation in Table 1. We used for
the calculations Matlab 2022a with the Parallel Computing Toolbox running on
Debian GNU/Linux 10 (buster), on a machine with the following specifications:
processor 2x AMD EPYC 7301 CPU @ 2.20 GHz, 256 GB RAM and a NVIDIA Tesla V100
PCIe (32 GB HBM2 RAM).

The next subsection is structured as follows: First, we derive the explicit
solution of the Langevin equation. Then, we discuss the impact of the number
of intervals $\Delta_{t}$ for the iterated Magnus scheme regarding
computational times and errors. Next, we look at the boundary effects over
time. This is followed by a comparison of the Magnus scheme with the Euler-
Maruyama scheme with different sizes of the space grid.

###  3.1 The Magnus expansion for the stochastic Langevin equation with
constant coefficients

In this subsection, we apply the Magnus expansion to the stochastic Langevin
equation. For further details and a solution theory in Hölder spaces under the
weak Hörmander condition we refer the reader to [16].

In the constant coefficient case, the Langevin SPDE can be recovered from
(2.1) setting

$\displaystyle h\equiv f^{v}\equiv g^{xx}\equiv
g^{xv}\equiv\sigma\equiv\sigma^{x}\equiv 0,\quad f_{x}(x,v)\coloneqq-v,\quad
g^{vv}\equiv
a\in\mathbb{R}_{>0},\quad\sigma^{v}(x,v)\equiv\sigma\in\mathbb{R}_{>0}.$ (3.4)

In this special case, there exists an explicit fundamental solution $\Gamma$
for $0<\sigma<\sqrt{a}$ (cf. [16, p. 4 Proposition 1.1.]), which is given by

$\displaystyle\Gamma\left(t,z;0,\zeta\right)$
$\displaystyle\coloneqq\Gamma_{0}\left(t,z-m_{t}(\zeta)\right),$
$\displaystyle\Gamma_{0}\big{(}t,(x,v)\big{)}$
$\displaystyle\coloneqq\frac{\sqrt{3}}{\pi
t^{2}(a-\sigma^{2})}\exp\left(-\frac{2}{a-\sigma^{2}}\left(\frac{v^{2}}{t}-\frac{3vx}{t^{2}}+\frac{3x^{2}}{t^{3}}\right)\right)$

where $\zeta\coloneqq\left(\xi,\eta\right)$ is the initial point and

$\displaystyle
m_{t}(\zeta)\coloneqq\left(\begin{array}[c]{c}\xi+t\eta-\sigma\int_{0}^{t}{W_{s}ds}\\\
\eta-\sigma W_{t}\end{array}\right).$

Having the fundamental solution, we can solve the Cauchy-problem by
integrating against the initial datum, i.e.

$\displaystyle
u_{t}(z)=\int_{\mathbb{R}^{2}}{\Gamma(t,z;0,\zeta)\phi(\zeta)d\zeta,\qquad
z=(x,v).}$

To get an explicit solution (up to the stochastic integral
$\int_{0}^{t}{W_{s}ds}$) for the double integral we will choose $\phi$ to be
Gaussian, i.e.

$\displaystyle\phi\left(\xi,\eta\right)\coloneqq\exp\left(-\frac{\left(\xi^{2}+\eta^{2}\right)}{2}\right).$
(3.5)

The formula for the exact solution is lengthy and its specific form is not
instructive for the following experiments, therefore we decided to exclude it
from this presentation. The interested reader can find it in the corresponding
Matlab 2022a code, which is publicly available.

Having an exact benchmark solution we will now perform some numerical tests to
judge the performance of the iterated Magnus scheme. Henceforth, the
parameters for the stochastic Langevin equation will be $a=1.1$ and
$\sigma=\frac{1}{\sqrt{10}}$, so that $a-\sigma^{2}=1>0$.

#### Computational effort and errors with respect to the number of iterations.

For this experiment we fix the number of grid points in each space grid to
$d=200$ but vary step-size of the Magnus scheme $\Delta_{t}$. In Figure 4 we
can see the corresponding results. The left y-axis shows the average
computational times for one simulation in a log scale and the right y-axis the
mean relative errors $\mathrm{Err}_{T}^{d,4}$ also in a log scale. The
computational times (in seconds) of m2 are depicted in light blue and of m3 in
dark blue. Moreover, the mean relative errors for m2 are orange and for m3
red.

Figure 4: Constant coefficients as in (3.4): Computational times and errors of
the Magnus expansion for varying number step-size $\Delta_{t}$ with fixed
spatial dimension $d=200$.

We can see that the mean relative errors start to fluctuate for a step-size
larger than 0.065 for m2 and see an explosion of m2 around 0.085. For m3 this
phenomenon starts outside of the picture. The fluctuations begin right after
$\Delta_{t}=0.1$ and an explosion can be seen after $\Delta_{t}=0.33$. The
explosions for large step-sizes are not surprising, since the step-size is
determined by the underlying stopping times for the convergence of the Magnus
scheme. The fluctuations beforehand are most likely due to an interplay
between error propagations due to larger step-sizes and the necessary Taylor-
terms in `expmv`, and it indicates that m3 is more stable than m2. Therefore,
experiment indicates that any step-size less than 0.05, 0.1 for $d=200$ is
well within the convergence radius of m2, m3, respectively, and yields stable
results.

For other spatial dimensions $d$ this breaking point might be different.
Moreover, we can see that the computational time increases more and more for
smaller step-size, while the error for both methods stays almost constant and
close to each other.

This suggests that one should choose the step-size as large as possible for
the iterated Magnus scheme to gain the maximal performance. However, being too
greedy will lead to blow-ups of some trajectories.

Also as a side note, usually with increasing spatial dimension $d$, one has to
choose a smaller time step-size for the Magnus methods as well: this will be
shown in Figure 6–8.

#### Mean errors and boundary effects over time.

For this experiment we fix the grid points in each space grid to $d=300$.

In Figure 5 we can see the mean absolute errors of the entire spatial grid as
a two-dimensional plot. A deep blue color indicates a small error and a bright
yellow color an error up to $10^{-2}$. The black rectangle is the
corresponding region for $\kappa=1$. The black number within the rectangle is
the average mean absolute error of the corresponding region. The picture on
the left-hand side is the area of errors at $t=0.25$ and on the right-hand
side at $t=1$.

Figure 5: Constant coefficients as in (3.4): Absolute Errors of m3 compared to
exact using $d=300$ grid points at $t=0.25$ (left) and $t=1$ (right).

We can see that the errors in the upper right and lower left corners are
significantly increasing over time. To explain this, one should note that the
Langevin equation with this specific initial datum looks like a two-
dimensional Gaussian at first and its shape changes on the diagonal from the
lower left corner to the upper right corner over time more than on the other
diagonal. Therefore, the cut-off region is getting too small for larger times
leading to boundary effects in the error plots. This also explains why the
upper left and lower right corner remain a stable small error. In the center
of the error plots we can see an increasing error over time. If this is due to
the boundary effects, error propagation due to the iterated scheme, the error
due to the order 3 truncation or the algorithm used for the matrix-vector
exponential is not apparent in this illustration, we suspect a mixture of all
of them.

#### Comparison to the Euler-Maruyama scheme

For this experiment, we will compare different choices of parameters for both
the Magnus scheme and Euler-Maruyama scheme. There are essentially two major
parameters contributing to the possible accuracy. One is the time step-size of
the individual schemes and the other one the space discretization. Hence, we
compare Euler and Magnus methods with different time step-sizes for different
space discretizations $d=100,200,300$ to increase the level of accuracy. In
the Figures 6, 7 and 8 the left y-axis shows the average computational times
in a log scale and the right y-axis the mean relative errors
$\mathrm{Err}_{T}^{d,4}$ also in a linear scale. The computational times (in
seconds) are depicted in the left blue columns and the mean relative errors in
the red right columns for each method.

As mentioned in Table 1, “E, x” denotes Euler with step-size $\Delta_{t}=x$
and “M2, x”, “M3, x” denotes Magnus with step-size $\Delta_{t}=x$ for order 2
and 3, respectively. In Figure 6 we compare the errors and computational times
of the methods with spatial dimension $d=100$, in Figure 7 with $d=200$ and in
Figure 8 with $d=300$.

Figure 6: Constant coefficients as in (3.4): Computational times and errors of
the Magnus expansion and Euler scheme for $d=100$.

Let us focus on Figure 6 with $d=100$. We can see that four different methods
are compared: the Euler method with step-size $\Delta_{t}=10^{-3}$ and
$\Delta_{t}=10^{-4}$, as well as the Magnus method with step-size
$\Delta_{t}=0.1$ of order 2 and order 3. It is notable that the Euler method
with step-size $\Delta_{t}=10^{-4}$ and the Magnus methods perform almost the
same with respect to the error. The Euler method with step-size
$\Delta_{t}=10^{-3}$ has roughly double the error of the method with step-size
$\Delta_{t}=10^{-4}$ but is ten times faster. Overall, the Magnus methods were
the fastest methods. The Magnus method of order two, three is $70$, $42$
times, respectively, faster than Euler method with step-size
$\Delta_{t}=10^{-4}$ and has a slightly better accuracy.

Figure 7: Constant coefficients as in (3.4): Computational times and errors of
the Magnus expansion and Euler scheme for $d=200$.

Now, let us consider Figure 7 with $d=200$. Again, we can see that two Euler
methods and two Magnus methods are compared but this time we have a step-size
$\Delta_{t}=0.05$ for the Magnus methods. Similarly, to Figure 6, we can see
that the Euler method with step-size $\Delta_{t}=10^{-3}$ performed worst and
the Magnus methods best in terms of accuracy. However, this time the Euler
method with step-size $\Delta_{t}=10^{-4}$ has twice the error compared to the
Magnus methods and is $25$, $23$ times slower than the Magnus method with
order 2, 3, respectively.

Figure 8: Constant coefficients as in (3.4): Computational times and errors of
the Magnus expansion and Euler scheme for $d=300$.

In Figure 8 with $d=300$ the Euler method with step-size $\Delta_{t}=10^{-3}$
is exploding, since its stability criterion is violated and its errors is
$\infty$. Therefore, we compare the Euler method with step-size
$\Delta_{t}=10^{-4}$ to the Magnus method with with step-size
$\Delta_{t}=0.025$ with order two and three. This time the Euler method is
four times worse in terms of accuracy and $30$, $27$ times slower than Magnus
with order two, three, respectively. We also performed tests with an Euler
method using step-size $\Delta_{t}=10^{-5}$. Its accuracy was still slightly
worse compared to the Magnus methods and its computational time ten times
slower than the Euler scheme in the figure. This results in a speed-up of
order 250 of the Magnus method compared to an Euler scheme with similar
accuracy.

Overall, from these observation it is clear that an Euler scheme with a fine
time-discretization is essential to make it comparable to the iterated Magnus
scheme in terms of accuracy. Moreover, increasing the number of grid points is
leading to less accurate errors using the Euler method with step-size
$\Delta_{t}=10^{-4}$ compared to the Magnus schemes with corresponding step-
sizes, while the Magnus methods remain roughly $30$ times faster in all tests.

######  Remark 3.1.

If we have a close look at all tests from above we can see that all of them
share a common feature, namely for reasonable parameters m2 and m3 were always
close. Therefore, this leads to a natural step-size control in time by
comparing the results of Magnus order 2 to order 3. If they are closer than a
given tolerance then the step-size is small enough, otherwise make it smaller
by a given factor.

With this method, the computation of the Magnus logarithms up to order 3 can
be reused for Magnus order 2. However, two matrix-vector exponentials for each
trajectory are necessary to determine if the time-step is rejected. For
implementations with a lot of trajectories, one can think about using less
randomly chosen trajectories to determine the correct step-size to increase
the overall performance.

######  Remark 3.2.

The Magnus expansion holds an advantage over all other finite-difference
method in the deterministic case. Inspecting the approximation formulas in the
case $A\equiv 0$ reveals immediately that the Magnus expansion is exact, at
order $1$, up to the initial space discretization for $x$ and $v$, meaning
that its accuracy is far more superior than e.g. explicit and implicit Euler-
schemes.

######  Remark 3.3.

We would like to point out, that it is straightforward to use the iterated
Magnus scheme in the case of deterministic piecewise constant coefficients in
time by partitioning the time interval according to the piecewise definition.
Then on each of these intervals the Magnus expansion formulas for constant
coefficients hold.

###  3.2 The Magnus expansion for the stochastic Langevin equation with
variable coefficients

In this brief subsection, we will perform some tests in the case of variable
coefficients. In particular, we choose bounded, smooth coefficients of the
form

$\displaystyle h\equiv f^{v}\equiv g^{xx}\equiv
g^{xv}\equiv\sigma\equiv\sigma^{x}\equiv 0,$ (3.6) $\displaystyle
f_{x}(x,v)\coloneqq-v,\quad
g^{vv}(x,v)=a\left(1+\frac{1}{x^{2}+1}\right),\quad\sigma^{v}(x,v)\equiv\sigma\sqrt{1+\frac{1}{x^{2}+1}},$
(3.7)

with $a,\sigma\in\mathbb{R}_{>0}$, as above, satisfying
$g^{vv}(x,v)-\left(\sigma^{v}(x,v)\right)^{2}\equiv 1>0$, as in the constant
coefficient case. We will also use the same initial condition as in (3.5).

Analog to Figure 4, we show in Figure 9 the case of varying step-sizes for the
Magnus method with fixed spatial discretization $d=200$. The average
computational times of m2 and m3 in seconds, per simulation, are again
depicted in light blue and dark blue, respectively. The errors are this time
with respect to the Euler method with $\Delta_{t}=10^{-4}$, since an exact
solution is not available, and again illustrated as orange for m2 and red for
m3.

Figure 9: Variable coefficients as in (3.7): Computational times and errors
compared to Euler with $\Delta_{t}=10^{-4}$ of the Magnus expansion for
varying number step-size $\Delta_{t}$ with fixed spatial dimension $d=200$.

The computational times look similar to the constant coefficient case and we
can notice a slight increase of the computational times for m3 at around
$0.055$. This behavior is the same in the constant coefficient case, but can
only be seen at step-sizes larger than $\Delta_{t}=0.1$. Therefore it is not
part of the previous figure. This can be explained by the necessity to use
more Taylor terms in the matrix-vector exponential to keep the errors small.

Moreover, we can see that the mean relative errors start to increase for a
step-size larger than 0.035 for m2 and 0.045 for m3. This happens earlier than
in the constant coefficient case and is the expected behavior.

In our next experiment, similarly to Figure 6–8, we show in Table 2 the mean
relative errors and computational times for $d=100,200,300$. This time we use
an Euler method with $\Delta_{t}=10^{-5}$ as our reference solution and
compare an Euler method with $\Delta_{t}=10^{-4}$, as well as m2 and m3 to it.
The results are given in Table 2.

Table 2: Variable coefficients as in (3.7): Computational times and errors compared to Euler with $\Delta_{t}=10^{-5}$ of the Magnus expansion for different spatial dimension $d=100,200,300$. Method | Mean Rel. Error (in %) | Comp. Time (in sec./simulation)
---|---|---
$d=100$
euler, $\Delta_{t}=10^{-3}$ | 0.147 % | 0.43
euler, $\Delta_{t}=10^{-4}$ | 0.047 % | 4.24
euler, $\Delta_{t}=10^{-5}$ | – | 42.24
m2, $\Delta_{t}=0.05$ | 0.015 % | 0.13
m3, $\Delta_{t}=0.05$ | 0.014 % | 0.24
$d=200$
euler, $\Delta_{t}=10^{-3}$ |  $\infty$ % | 2.03
euler, $\Delta_{t}=10^{-4}$ | 0.045 % | 21.02
euler, $\Delta_{t}=10^{-5}$ | – | 209.02
m2, $\Delta_{t}=0.025$ | 0.013 % | 1.66
m3, $\Delta_{t}=0.025$ | 0.012 % | 1.83
$d=300$
euler, $\Delta_{t}=10^{-3}$ |  $\infty$ % | 5.79
euler, $\Delta_{t}=10^{-4}$ | 0.041 % | 58.41
euler, $\Delta_{t}=10^{-5}$ | – | 582.88
m2, $\Delta_{t}=0.01$ | 0.014 % | 4.28
m3, $\Delta_{t}=0.01$ | 0.014 % | 5.0

We can see that computational times are roughly twice higher for the Magnus
scheme than in the constant coefficient case. The Euler scheme did not suffer
from a performance decrease. As aforementioned this increase of computational
effort can be attributed to the decreased sparsity. Therefore, the Magnus
scheme is only 12 times faster than the Euler scheme. Regarding the accuracy,
we can see in all cases that the Magnus expansion is more accurate than euler
with $\Delta_{t}=10^{-4}$.

All in all we come to the same conclusion as in the constant coefficient case.

##  4 Conclusions

We have seen how to derive the Itô-stochastic Magnus expansion for SDEs with
constant matrices and used it to solve two-dimensional SPDEs with a given
initial datum numerically. We derived an approximation scheme for a generic
parabolic SPDE with two space variables, then tested it on a special class of
Langevin-type SPDEs. The scheme has an excellent accuracy and its advantage in
terms of computational effort excels for higher spatial resolution. For
instance, in the case of constant coefficients to have roughly the same
accuracy for $d=400$ we need an Euler scheme with $\Delta_{t}=10^{-5}$ taking
approximately $21$ minutes, in average per trajectory, while Magnus order 3
takes only $4.62$ seconds using $\Delta_{t}=0.01$,
$\Delta^{\mathrm{Leb}}_{t}=10^{-4}$. This is a speed-up by a factor 280 just
using one GPU while sparsity ensures an almost equal memory demand. Adding
multiple GPUs will lead to further improvements with respect to computational
times.

##  Appendix A Magnus expansion formulas for constant coefficients

In order to make the paper self-contained, in this section, we demonstrate how
to heuristically derive the Itô-stochastic Magnus expansion formula for the
equation

$\displaystyle dX_{t}=BX_{t}dt+AX_{t}dW_{t},\quad X_{0}=I_{d},\quad$ (A.1)

where the coefficients $B,A\in\mathbb{R}^{d\times d}$ are constant matrices
and $W$ is a standard one-dimensional Brownian motion.

We make the ansatz $X_{t}=\exp\left(Y_{t}\right)$, where

$\displaystyle
Y_{t}=\int_{0}^{t}{\mu\left(s,Y_{s}\right)ds}+\int_{0}^{t}{\sigma\left(s,Y_{s}\right)dW_{s}},\quad
Y_{0}=0_{\mathbb{R}^{d\times d}}.$

By Itô’s formula (cf. [10, p. 8 Lemma 1 and p. 9 Proposition 1]) we have we
have

$\displaystyle dX_{t}$ $\displaystyle=BX_{t}+AX_{t}dW_{t}$
$\displaystyle=d\exp\left(Y_{t}\right)$ $\displaystyle=\begin{aligned}
&\biggl{(}\mathcal{L}_{Y_{t}}\left(\mu\left(s,Y_{s}\right)\right)+\frac{1}{2}\mathcal{Q}_{Y_{t}}\left(\sigma\left(t,Y_{t}\right),\sigma\left(t,Y_{t}\right)\right)\biggr{)}\exp\left(Y_{t}\right)dt\\\
&+\mathcal{L}_{Y_{t}}\left(\sigma\left(t,Y_{t}\right)\right)\exp\left(Y_{t}\right)dW_{t},\end{aligned}$

where the operators $\mathcal{L}$ and $\mathcal{Q}$ associated to the first
and second order derivative of the exponential map are given by

$\displaystyle\mathcal{L}_{Y}(M)$
$\displaystyle\coloneqq\sum_{n=0}^{\infty}\frac{1}{(n+1)!}\text{ad}^{n}_{Y}(M),$
$\displaystyle\mathcal{Q}_{Y}(M,N)$
$\displaystyle\coloneqq\sum_{n=0}^{\infty}{\sum_{m=0}^{\infty}{\frac{\text{ad}^{n}_{Y}(M)}{(n+1)!}\frac{\text{ad}^{m}_{Y}(N)}{(m+1)!}}}+\sum_{n=0}^{\infty}{\sum_{m=0}^{\infty}{\frac{\left[\text{ad}^{n}_{Y}(N),\text{ad}^{m}_{Y}(M)\right]}{(n+m+2)(n+1)!m!}}}.$

A comparison of coefficients yields

$\displaystyle B$
$\displaystyle\overset{!}{=}\mathcal{L}_{Y_{t}}\left(\mu\left(t,Y_{s}\right)\right)+\frac{1}{2}\mathcal{Q}_{Y_{t}}\left(\sigma\left(t,Y_{t}\right),\sigma\left(t,Y_{t}\right)\right)$
$\displaystyle A$
$\displaystyle\overset{!}{=}\mathcal{L}_{Y_{t}}\left(\sigma\left(t,Y_{t}\right)\right).$

Now, Baker’s lemma (cf. [10, p. 10 Lemma 2]) provides us with suitable
conditions and a series representation for the inverse of the first-order
derivative operator $\mathcal{L}$, i.e.

$\displaystyle\sigma\left(t,Y_{t}\right)$
$\displaystyle\equiv\sigma\left(Y_{t}\right)=\sum_{n=0}^{\infty}{\frac{\beta_{n}}{n!}\text{ad}^{n}_{Y_{t}}\left(A\right)}$
(A.2) $\displaystyle\mu\left(t,Y_{t}\right)$
$\displaystyle\equiv\mu\left(Y_{t}\right)=\begin{aligned}
\sum_{k=0}^{\infty}\frac{\beta_{k}}{k!}\text{ad}^{k}_{Y_{t}}\Biggl{(}B-\frac{1}{2}\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}&\frac{\text{ad}^{n}_{Y_{t}}\left(\sigma\left(Y_{t}\right)\right)}{(n+1)!}\frac{\text{ad}^{m}_{Y_{t}}\left(\sigma\left(Y_{t}\right)\right)}{(m+1)!}\\\
&+\frac{\left[\text{ad}^{n}_{Y_{t}}\left(\sigma\left(Y_{t}\right)\right),\text{ad}^{m}_{Y_{t}}\left(\sigma\left(Y_{t}\right)\right)\right]}{(n+m+2)(n+1)!m!}\Biggr{)},\end{aligned}$
(A.3)

The Bernoulli numbers are denoted by $\beta_{k}$ and we recall that
$\beta_{0}=1$, $\beta_{1}=-\frac{1}{2}$, $\beta_{2}=\frac{1}{6}$,
$\beta_{3}=0$ and $\beta_{4}=-\frac{1}{30}$.

Now, we know the coefficients of $Y_{t}$ and solve the Itô-SDE by means of a
stochastic Picard iteration, i.e. we start at $n=0$ with
$Y_{t}^{0}=0_{\mathbb{R}^{d\times d}}$ and iterate for $n\geq 1$

$\displaystyle
Y_{t}^{n}=\int_{0}^{t}{\mu\left(s,Y_{s}^{n-1}\right)ds}+\int_{0}^{t}{\sigma\left(s,Y_{s}^{n-1}\right)dW_{s}}.$
(A.4)

In order to derive the Magnus expansion formulas we will introduce some
bookkeeping parameters $\epsilon,\delta>0$ and substitute $A$ by $\epsilon A$,
as well as $B$ by $\delta B$. Henceforth, we will denote the $n$-th order
Picard iteration with the substitution by $Y^{n,\epsilon,\delta}_{t}$.

#### Order 1.

Let us derive the first-order Magnus expansion, meaning we are interested in
all terms with the first power of $\epsilon$ and $\delta$ in the first-order
Picard iteration. Thus, we insert $Y^{0,\epsilon,\delta}_{t}\equiv 0$ into
(A.4). Notice, that the zero matrix commutes with all matrices and therefore
by definition of $\text{ad}^{0}_{Y}(A)=A$ we have

$\displaystyle\sigma\left(Y^{0,\epsilon,\delta}_{t}\right)=\epsilon A.$

Inserting this into the formula for $\mu$ yields

$\displaystyle\mu\left(Y^{0,\epsilon,\delta}_{t}\right)=\delta
B-\frac{1}{2}\epsilon^{2}A^{2},$

because $A$ commutes with itself as well.

Since, the Itô-correction term is of order $\epsilon^{2}$ it will not be part
of the first-order Magnus expansion and we have

$\displaystyle Y^{1}_{t}=\int_{0}^{t}{Bds}+\int_{0}^{t}{AdW_{s}}=Bt+AW_{t}.$

#### Order 2.

To derive the second-order Magnus expansion let us first think about how many
terms we need at most in the infinite sums of (A.2) and (A.3) to disregard the
majority of the high-order terms.

Let us consider (A.2) first and let $N=2$ denote the desired order of the
Magnus expansion. The operator $\text{ad}^{0}$ will at least result in a first
order term in $\epsilon$ or $\delta$. The commutator
$\text{ad}^{1}_{Y}\left(X\right)=\left[Y,X\right]=YX-XY$ will at least have
second order terms, because $X$ and $Y$ will at least have a first-order
$\epsilon$ or $\delta$ coefficient. Now, the nested commutator
$\text{ad}^{2}_{Y}\left(X\right)=\left[Y,\left[Y,X\right]\right]$ will have at
least a third-order term for the same reason. Therefore, to compute the $N$-th
order of the Magnus expansion we only need to consider the infinite sum in
(A.2) up to $N-1$. Now, let us consider (A.3). We will split our consideration
into two parts by linearity and focus for the moment on

$\displaystyle\sum_{k=0}^{\infty}{\text{ad}_{Y_{t}}^{k}\left(B\right)}.$

With exactly the same arguments as for (A.2) we only need to consider $N-1$
terms of this infinite sum for the $N$-th order Magnus expansion.

Let us turn now to the rest, the Itô-correction term, of (A.3), and split it
again in two parts, i.e. we consider

$\displaystyle\sum_{k=0}^{\infty}\text{ad}^{k}_{Y_{t}}\left(\sum_{n=0}^{\infty}\sum_{m=0}^{\infty}\frac{\text{ad}_{Y_{t}}^{n}\left(\sigma\left(Y_{t}\right)\right)}{(n+1)!}\frac{\text{ad}_{Y_{t}}^{m}\left(\sigma\left(Y_{t}\right)\right)}{(m+1)!}\right).$

In the case $0=k=n=m$ we will get $\sigma^{2}$, which has at least order 2
terms. Hence, we will only need $k=1,\dots,N-2$. This also means, that $n$ and
$m$ cannot exceed $N-2$.

Increasing $k$ and keeping $n=m=0$ will also increase the order due to the
nested commutators as in the case for (A.2). Therefore, the lowest $k$ can
have the most terms resulting from $n$ and $m$. Vice versa, the higher $n$ and
$m$ the lower $k$ must be to have small enough orders.

Taking all these considerations into account, we can use the following formula
to cover all necessary terms up until order $N$:

$\displaystyle\sigma_{N}\left(Y_{t}\right)$
$\displaystyle\coloneqq\sum_{n=0}^{N-1}{\frac{\beta_{n}}{n!}\text{ad}^{n}_{Y_{t}}\left(A\right)},$
(A.5) $\displaystyle\mu_{N}\left(Y_{t}\right)$
$\displaystyle\coloneqq\sum_{k=0}^{N-1}\frac{\beta_{k}}{k!}\text{ad}^{k}_{Y_{t}}\left(B\right)$
(A.6) $\displaystyle\quad\begin{aligned}
-\frac{1}{2}\sum_{k=0}^{N-2}\frac{\beta_{k}}{k!}\text{ad}^{k}_{Y_{t}}\Biggl{(}&\sum_{n=0}^{N-2-k}\sum_{m=0}^{N-2-k}\frac{\text{ad}^{n}_{Y_{t}}\left(\sigma_{N-k}\left(Y_{t}\right)\right)}{(n+1)!}\frac{\text{ad}^{m}_{Y_{t}}\left(\sigma_{N-k}\left(Y_{t}\right)\right)}{(m+1)!}\\\
&+\frac{\left[\text{ad}^{n}_{Y_{t}}\left(\sigma_{N-k}\left(Y_{t}\right)\right),\text{ad}^{m}_{Y_{t}}\left(\sigma_{N-k}\left(Y_{t}\right)\right)\right]}{(n+m+2)(n+1)!m!}\Biggr{)}.\end{aligned}$
(A.7)

There will still be a lot of higher-order terms in (A.5) and (A.7) but it will
make our derivation of the formulas easier.

Now, we compute $\sigma_{2}\left(Y_{t}^{1}\right)$:

$\displaystyle\sigma_{2}\left(Y_{t}^{1,\epsilon,\delta}\right)$
$\displaystyle=\frac{1}{0!}\text{ad}^{0}_{Y_{t}^{1,\epsilon,\delta}}\left(\epsilon
A\right)-\frac{1}{2}\text{ad}^{1}_{Y_{t}^{1,\epsilon,\delta}}\left(\epsilon
A\right)$ $\displaystyle=\epsilon
A-\frac{1}{2}\left[Y_{t}^{1,\epsilon,\delta},\epsilon A\right]$
$\displaystyle=\epsilon A-\frac{1}{2}\left[\delta B,\epsilon
A\right]t-\frac{1}{2}\left[\epsilon A,\epsilon A\right]W_{t}$
$\displaystyle=\epsilon A-\frac{1}{2}\left[\delta B,\epsilon A\right]t.$

Next, we compute $\mu_{2}\left(Y_{t}^{1}\right)$:

$\displaystyle\mu_{2}\left(Y_{t}^{1,\epsilon,\delta}\right)=\frac{1}{0!}\text{ad}^{0}_{Y_{t}^{1,\epsilon,\delta}}\left(\delta
B\right)-\frac{1}{2}\text{ad}^{1}_{Y_{t}^{1,\epsilon,\delta}}\left(\delta
B\right)$
$\displaystyle\quad-\frac{1}{2}\Biggl{(}\frac{1}{0!}\biggl{(}\frac{\text{ad}^{0}_{Y_{t}^{1,\epsilon,\delta}}\left(\sigma_{2}\left(Y_{t}^{1,\epsilon,\delta}\right)\right)}{1!}\frac{\text{ad}^{0}_{Y_{t}^{1,\epsilon,\delta}}\left(\sigma_{2}\left(Y_{t}^{1,\epsilon,\delta}\right)\right)}{1!}$
$\displaystyle\qquad+\frac{\left[\text{ad}^{0}_{Y_{t}^{1,\epsilon,\delta}}\left(\sigma_{2}\left(Y_{t}^{1,\epsilon,\delta}\right)\right),\text{ad}^{0}_{Y_{t}^{1,\epsilon,\delta}}\left(\sigma_{2}\left(Y_{t}^{1,\epsilon,\delta}\right)\right)\right]}{2}\biggr{)}\Biggr{)}$
$\displaystyle=\delta B-\frac{1}{2}\left[Y_{t}^{1,\epsilon,\delta},\delta
B\right]$ $\displaystyle\quad-\frac{1}{2}\Biggl{(}\left(\epsilon
A-\frac{1}{2}\left[\delta B,\epsilon A\right]t\right)^{2}+\frac{\left[\epsilon
A-\frac{1}{2}\left[\delta B,\epsilon A\right]t,\epsilon
A-\frac{1}{2}\left[\delta B,\epsilon A\right]t\right]}{2}\Biggr{)}.$

Counting $\epsilon$ and $\delta$ we can already see that

$\displaystyle\left(\epsilon A-\frac{1}{2}\left[\delta B,\epsilon
A\right]t\right)^{2}$ $\displaystyle\approx\epsilon^{2}A^{2},$
$\displaystyle\left[\epsilon A-\frac{1}{2}\left[\delta B,\epsilon
A\right]t,\epsilon A-\frac{1}{2}\left[\delta B,\epsilon A\right]t\right]$
$\displaystyle\approx\left[\epsilon A,\epsilon A\right]=0.$

Thus,

$\displaystyle\mu\left(Y_{t}^{1,\epsilon,\delta};2\right)=\delta
B-\frac{1}{2}\left[\epsilon A,\delta
B\right]W_{t}-\frac{1}{2}\epsilon^{2}A^{2}.$

In total, we have

$\displaystyle
Y_{t}^{2}=\int_{0}^{t}{B-\frac{1}{2}\left[A,B\right]W_{s}-\frac{1}{2}A^{2}ds}+\int_{0}^{t}{A-\frac{1}{2}\left[B,A\right]sdW_{s}}.$

In a next step, we can apply Itô’s formula to replace the stochastic integral
by a Lebesgue integral like follows:

######  Lemma A.1.

For $p,p_{1},p_{2},q,q_{1},q_{2}\in\mathbb{N}_{0}$ we have

$\displaystyle\int_{0}^{t}{s^{p}W_{s}^{q}dW_{s}}=\frac{1}{q+1}\left(t^{p}W_{t}^{q+1}-\int_{0}^{t}{\left[\frac{q(q+1)}{2}s^{p}W_{s}^{q-1}+pW_{s}^{q+1}s^{p-1}\right]ds}\right)$
(A.8)
$\displaystyle\int_{0}^{t}{s^{p_{1}}\int_{0}^{s}{r^{p_{2}}W_{r}^{q}dr}ds}=\frac{1}{1+p_{1}}\left(t^{1+p_{1}}\int_{0}^{t}{s^{p_{2}}W_{s}^{q}ds}-\int_{0}^{t}{s^{1+p_{1}+p_{2}}W_{s}^{q}ds}\right)$
(A.9)
$\displaystyle\int_{0}^{t}{s^{p_{1}}W_{s}^{q_{1}}\int_{0}^{s}{r^{p_{2}}W_{r}^{q_{2}}dr}dW_{s}}=\begin{aligned}
\frac{1}{q_{1}+1}\Biggl{(}&t^{p_{1}}W^{q_{1}+1}_{t}\int_{0}^{t}{s^{p_{2}}W_{s}^{q_{2}}ds}\\\
&-\int_{0}^{t}{s^{p_{1}+p_{2}}W_{s}^{q_{1}+q_{2}+1}ds}\\\
&-\frac{q_{1}(q_{1}+1)}{2}\int_{0}^{t}{s^{p_{1}}W_{s}^{q_{1}-1}\int_{0}^{s}{r^{p_{2}}W_{r}^{q_{2}}dr}ds}\\\
&-p_{1}\int_{0}^{t}{s^{p_{1}-1}W^{q_{1}+1}_{s}\int_{0}^{s}{r^{p_{2}}W_{r}^{q_{2}}dr}ds}\Biggr{)}.\end{aligned}.$
(A.10)

###### Proof.

Note that (A.8) is a special case of (A.10) by setting $p=p_{1}-1$, $q=q_{1}$
and $p_{2}=q_{2}=0$.

Now, we show (A.9). With Itô’s product rule we get

$\displaystyle d\left(s^{p_{1}}\int_{0}^{s}{r^{p_{2}}W_{r}^{q}dr}\ s\right)$
$\displaystyle=s^{p_{1}}\int_{0}^{s}{r^{p_{2}}W_{r}^{q}dr}ds+s\
d\left(s^{p_{1}}\int_{0}^{s}{r^{p_{2}}W_{r}^{q}dr}\right)+0$
$\displaystyle=s^{p_{1}}\int_{0}^{s}{r^{p_{2}}W_{r}^{q}dr}ds+s\
d\left(s^{p_{1}}s^{p_{2}}W_{s}^{q}ds+\int_{0}^{s}{r^{p_{2}}W_{r}^{q}dr}p_{1}s^{p_{1}-1}ds\right)$
$\displaystyle=\left(1+p_{1}\right)s^{p_{1}}\int_{0}^{s}{r^{p_{2}}W_{r}^{q}dr}ds+s^{1+p_{1}+p_{2}}W_{s}^{q}ds.$

Rearranging the equation yields the claim.

Next, we show (A.10). With Itô’s product rule we get

$\displaystyle
d\left(s^{p_{1}}W_{s}^{q_{1}}\int_{0}^{s}{r^{p_{2}}W_{r}^{q_{2}}dr}W_{s}\right)$
$\displaystyle=s^{p_{1}}W_{s}^{q_{1}}\int_{0}^{s}{r^{p_{2}}W_{r}^{q_{2}}dr}dW_{s}+W_{s}d\left(s^{p_{1}}W_{s}^{q_{1}}\int_{0}^{s}{r^{p_{2}}W_{r}^{q_{2}}dr}\right)$
$\displaystyle\quad+d\left\langle\cdot^{p_{1}}W_{\cdot}^{q_{1}}\int_{0}^{\cdot}{r^{p_{2}}W_{r}^{q_{2}}dr},W_{\cdot}\right\rangle_{s}.$

Now, use Itô’s product rule and Itô’s formula on the the following term

$\displaystyle
d\left(s^{p_{1}}W_{s}^{q_{1}}\int_{0}^{s}{r^{p_{2}}W_{r}^{q_{2}}dr}\right)=s^{p_{1}}W_{s}^{q_{1}}d\left(\int_{0}^{s}{r^{p_{2}}W_{r}^{q_{2}}dr}\right)+\int_{0}^{s}{r^{p_{2}}W_{r}^{q_{2}}dr}d\left(s^{p_{1}}W_{s}^{q_{1}}\right)+0$
$\displaystyle=s^{p_{1}+p_{2}}W_{s}^{q_{1}+q_{2}}ds+\int_{0}^{s}{r^{p_{2}}W_{r}^{q_{2}}dr}\Biggl{(}s^{p_{1}}\biggl{(}q_{1}W_{s}^{q_{1}-1}dW_{s}+\frac{q_{1}(q_{1}-1)}{2}W_{s}^{q_{1}-2}ds\biggr{)}$
$\displaystyle\qquad+W_{s}^{q_{1}}\biggl{(}p_{1}s^{p_{1}-1}ds\biggr{)}+0\Biggr{)}$
$\displaystyle=\int_{0}^{s}{r^{p_{2}}W_{r}^{q_{2}}dr}s^{p_{1}}q_{1}W_{s}^{q_{1}-1}dW_{s}+\Biggl{[}s^{p_{1}+p_{2}}W_{s}^{q_{1}+q_{2}}+\int_{0}^{s}{r^{p_{2}}W_{r}^{q_{2}}dr}s^{p_{1}}\frac{q_{1}(q_{1}-1)}{2}W_{s}^{q_{1}-2}$
$\displaystyle\qquad+\int_{0}^{s}{r^{p_{2}}W_{r}^{q_{2}}dr}W_{s}^{q_{1}}p_{1}s^{p_{1}-1}\Biggr{]}ds$

For the quadratic variation from above we have

$\displaystyle
d\left\langle\cdot^{p_{1}}W_{\cdot}^{q_{1}}\int_{0}^{\cdot}{r^{p_{2}}W_{r}^{q_{2}}dr},W_{\cdot}\right\rangle_{s}$
$\displaystyle=\int_{0}^{s}{r^{p_{2}}W_{r}^{q_{2}}dr}s^{p_{1}}q_{1}W_{s}^{q_{1}-1}ds.$

In total, we have

$\displaystyle
d\left(s^{p_{1}}W_{s}^{q_{1}}\int_{0}^{s}{r^{p_{2}}W_{r}^{q_{2}}dr}W_{s}\right)$
$\displaystyle=(1+q_{1})s^{p_{1}}W_{s}^{q_{1}}\int_{0}^{s}{r^{p_{2}}W_{r}^{q_{2}}dr}dW_{s}+\Biggl{[}s^{p_{1}+p_{2}}W_{s}^{q_{1}+q_{2}+1}+\frac{q_{1}(q_{1}+1)}{2}s^{p_{1}}W_{s}^{q_{1}-1}\int_{0}^{s}{r^{p_{2}}W_{r}^{q_{2}}dr}$
$\displaystyle\qquad+p_{1}\int_{0}^{s}{r^{p_{2}}W_{r}^{q_{2}}dr}W_{s}^{q_{1}+1}s^{p_{1}-1}\Biggr{]}ds.$

Rearranging the equation yields the claim. ∎

Using Lemma A.1 (A.8) in the case $p=1$ and $q=0$ and the skew-symmetry of the
commutator we have finally

$\displaystyle Y_{t}^{2}$
$\displaystyle=Bt-\frac{1}{2}A^{2}t+\frac{1}{2}\left[B,A\right]\int_{0}^{t}{W_{s}ds}+AW_{t}-\frac{1}{2}\left[B,A\right]\left(tW_{t}-\int_{0}^{t}{W_{s}ds}\right)$
$\displaystyle=Y_{t}^{1}-\frac{1}{2}A^{2}t+\left[B,A\right]\int_{0}^{t}{W_{s}ds}-\frac{1}{2}\left[B,A\right]tW_{t}.$

#### Order 3.

From now on we will always repeat the steps seen from the derivation of order
2. First, identify all terms using (A.5) and (A.7) and secondly apply Lemma
A.1 to remove the iterated (stochastic) integrals by expressions with single
Lebesgue integrals.

After using (A.5) and (A.7) we get

$\displaystyle\sigma_{3}\left(Y_{t}^{2,\epsilon,\delta}\right)$
$\displaystyle=\epsilon A-\frac{1}{2}\left[\delta B,\epsilon
A\right]t-\frac{1}{12}\left[\left[\delta B,\epsilon
A\right],B\right]t^{2}-\frac{1}{2}\left[\left[\delta B,\epsilon
A\right],\epsilon A\right]\int_{0}^{t}{W_{s}ds}$
$\displaystyle\quad+\frac{1}{6}\left[\left[\delta B,\epsilon A\right],\epsilon
A\right]W_{t}$ $\displaystyle\mu_{3}\left(Y_{t}^{2,\epsilon,\delta}\right)$
$\displaystyle=\delta
B-\frac{1}{2}\epsilon^{2}A^{2}+\frac{1}{12}\left[\left[\delta B,\epsilon
A\right],\epsilon A\right]t-\frac{1}{2}\left[\left[\delta B,\epsilon
A\right],\delta B\right]\int_{0}^{t}{W_{s}ds}+\frac{1}{2}\left[\delta
B,\epsilon A\right]W_{t}$ $\displaystyle\quad+\frac{1}{3}\left[\left[\delta
B,\epsilon A\right],\delta B\right]tW_{t}+\frac{1}{12}\left[\left[\delta
B,\epsilon A\right],\epsilon A\right]W_{t}^{2}$

Applying Lemma A.1 and collecting all terms yields

$\displaystyle Y_{t}^{3}=Y_{t}^{2}$
$\displaystyle+\left[\left[B,A\right],A\right]\left(\frac{1}{2}\int_{0}^{t}{W_{s}^{2}ds}-\frac{1}{2}W_{t}\int_{0}^{t}{W_{s}ds}+\frac{1}{12}tW_{t}^{2}\right)$
$\displaystyle+\left[\left[B,A\right],B\right]\left(\int_{0}^{t}{sW_{s}ds}-\frac{1}{2}t\int_{0}^{t}{W_{s}ds}-\frac{1}{12}t^{2}W_{t}\right).$

## References

  * [1] Al-Mohy, A. H., and Higham, N. J. Computing the action of the matrix exponential, with an application to exponential integrators.  SIAM J. Sci. Comput. 33, 2 (2011), 488–511.
  * [2] Barrera, G., Högele, M., and Pardo, J. Cutoff stability of multivariate geometric Brownian motion.  Preprint arXiv: 2207.01666.
  * [3] Campana, L. Stochastic modelling of non-spherical particles in turbulence.  PhD Thesis Université Côte d’Azur et de INRIA, 2022.
  * [4] Cercignani, C. The Boltzmann equation and its applications.  Springer-Verlag, New York, 1988.
  * [5] Desvillettes, L., and Villani, C. On the trend to global equilibrium in spatially inhomogeneous entropy-dissipating systems: the linear Fokker-Planck equation.  Comm. Pure Appl. Math. 54, 1 (2001), 1–42.
  * [6] Di Francesco, M., and Pascucci, A. On the complete model with stochastic volatility by Hobson and Rogers.  Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 460, 2051 (2004), 3327–3338.
  * [7] Friz, P. K., Hager, P. P., and Tapia, N. Unified signature cumulants and generalized Magnus expansions.  Forum Math. Sigma 10 (2022), Paper No. e42, 60.
  * [8] Ibáñez, J., Alonso, J. M., Alonso-Jordá, P., Defez, E., and Sastre, J. Two Taylor algorithms for computing the action of the matrix exponential on a vector.  Algorithms 15, 2 (2022).
  * [9] Kamm, K., and Muniz, M. A novel approach to rating transition modelling via Machine Learning and SDEs on Lie groups.  Preprint arXiv: 2205.15699.
  * [10] Kamm, K., Pagliarani, S., and Pascucci, A. On the stochastic Magnus expansion and its application to SPDEs.  J. Sci. Comput. 89, 3 (2021), Paper No. 56, 31.
  * [11] Lions, P.-L. On Boltzmann and Landau equations.  Philos. Trans. Roy. Soc. London Ser. A 346, 1679 (1994), 191–204.
  * [12] Magnus, W. On the exponential solution of differential equations for a linear operator.  Comm. Pure Appl. Math. 7 (1954), 649–673.
  * [13] Muniz, M., Ehrhardt, M., Günther, M., and Winkler, R. Higher strong order methods for linear Itô SDEs on matrix Lie groups.  To appear in BIT Numer Math (2022), https://doi.org/10.1007/s10543–022–00911–5.
  * [14] Pascucci, A. PDE and martingale methods in option pricing.  Bocconi&Springer Series. Springer-Verlag, New York, 2011.
  * [15] Pascucci, A., and Pesce, A. Backward and forward filtering under the weak Hörmander condition.  To appear in Stoch. Partial Differ. Equ. Anal. Comput. (2022).
  * [16] Pascucci, A., and Pesce, A. On stochastic Langevin and Fokker-Planck equations: the two-dimensional case.  J. Differential Equations 310 (2022), 443–483.
  * [17] Yang, G., Burrage, K., Komori, Y., Burrage, P., and Ding, X. A class of new Magnus-type methods for semi-linear non-commutative Itô stochastic differential equations.  Numer. Algorithms 88, 4 (2021), 1641–1665.