Numerical solution of kinetic SPDEs via stochastic Magnus expansion

In this paper, we show how the It\^o-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.


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 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 1 2 where the variables t ≥ 0, x ∈ R and v ∈ 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 where a t , b t , c t , σ t and β 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 (Ω, F, P) endowed with a filtration (F t ) t≥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 parallelin-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 Yt 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 R d×d and let (Ω, F, P, (F t ) t≥0 ) be a filtered probability space equipped with a standard Brownian motion W .For T > 0 let also X = (X t ) t∈[0,T ] be the unique strong solution to (1.1).There exists a strictly positive stopping time τ ≤ T such that: 1. Introduction (i) X t has a real logarithm Y t ∈ R d×d up to time τ , i.e.
X t = e Yt , 0 ≤ t < τ ; (ii) the following representation holds P-almost surely: t , 0 ≤ t < τ, (1.4) where Y (n) is the n-th term in the stochastic ME (see also formulas (1.6) below, in the case of constant matrices A and B); (iii) there exists a positive constant C, only dependent A, B, T and d, such that The first point of Theorem t , 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: where [A, B] := 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 erroranalysis 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.

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 (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 matrixvalued equation.We introduce the following two homogeneous grids for position and velocity respectively: X nx ax,bx with n x + 2 points on the subset [a Let us denote by vec the isomorphism of transforming a matrix to a larger column-vector by stacking each column in the matrix below each other, i.e. vec : ) 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.
The fundamental solution now can be approximated with the Magnus expansion, i.e.X t ≈ exp (Y t ).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 matrixexponential 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 tridiag nv,nv (−1, 0, 1) , Additionally, we will introduce the following matrices corresponding to the coefficient functions on the discretized spatial grid 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 for all i = 1, . . ., n x and j = 1, . . ., n v .As aforementioned, we need to extract the correct coefficient matrix for the vectorized equation (2.2).
In our notations, a derivative in x is a multiplication of the corresponding finite-difference matrix from the left to u nx,nv t , i.e.
Using the Kronecker product, it is well-known for compatible matrices In our case, this leads to Now, we need to deal with the coefficients as well.Denoting by the Hadamard, or elementwise, product, it is easy to see that Using these two observations together yields 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 Remark 2.1.We can see that the sparsity of B and A is mostly determined by the finitedifferences 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  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 where ∆t := t k+1 − t k > 0 for any k and ∆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.

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.6).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 ≥ 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 ≤ 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 ex- pensive 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 ref and by U app a benchmark and an approximate solution, respectively.
Henceforth, we will choose V nv av,bv = X nx ax,bx symmetric, centered around zero and set d := n x = n v to make our analysis a bit easier.Moreover, we choose the cut-off region of R 2 as 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 κ = 0, 1, . . .indicating the 2 −κ -th central part of the solution matrix U d t , which we will consider for our error analysis and denote the corresponding truncated matrix by to collect the corresponding indices.For our error analysis we will use the following three norms: This can be explained by the fact that the discretization is only used to compute Lebesgueintegrals of the form t 0 s p W q s 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 ∆ Leb t ≥ 10 −5 .
Since, there seems to be no computational drawback, we suggest to use a Lebesguediscretization equal to 10 −3 or 10 −4 and fix it for all tests to ∆ 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.
Throughout the experiments we will use the notation in Table 1.We used for the calcu- 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 ∆ 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.

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 In this special case, there exists an explicit fundamental solution Γ for 0 < σ < √ a (cf.[16, p. 4 Proposition 1.1.]),which is given by where ζ := (ξ, η) is the initial point and Having the fundamental solution, we can solve the Cauchy-problem by integrating against the initial datum, i.e.
To get an explicit solution (up to the stochastic integral t 0 W s ds) for the double integral we will choose φ to be Gaussian, i.e.
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 σ = 1 √ 10 , so that a − σ 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 ∆ 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 Err d, 4  T 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.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 ∆ t = 0.1 and an explosion can be seen after ∆ 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 twodimensional 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 κ = 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.
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 stepsize of the individual schemes and the other one the space discretization.Hence, we compare   T 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 ∆ t = x and "M2, x", "M3, x" denotes Magnus with step-size ∆ 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.
Let us focus on Figure 6 with d = 100.We can see that four different methods are compared: the Euler method with step-size ∆ t = 10 −3 and ∆ t = 10 −4 , as well as the Magnus method with step-size ∆ t = 0.1 of order 2 and order 3.It is notable that the Euler method with step-size ∆ t = 10 −4 and the Magnus methods perform almost the same with respect to the error.The Euler method with step-size ∆ t = 10 −3 has roughly double the error of the method with step-size ∆ 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 ∆ t = 10 −4 and has a slightly better accuracy.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 ∆ t = 0.05 for the Magnus methods.Similarly, to Figure 6, we can see that the Euler method with step-size ∆ t = 10 −3 performed worst and the Magnus methods best in terms of accuracy.However, this time the Euler method with step-size ∆ 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.
In Figure 8 with d = 300 the Euler method with step-size ∆ t = 10 −3 is exploding, since its stability criterion is violated and its errors is ∞.Therefore, we compare the Euler method  with step-size ∆ t = 10 −4 to the Magnus method with with step-size ∆ 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 ∆ 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 ∆ t = 10 −4 compared to the Magnus schemes with corresponding stepsizes, 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 ≡ 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.

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 3.2.The Magnus expansion for the stochastic Langevin equation with variable coefficients with a, σ ∈ R >0 , as above, satisfying g vv (x, v) − (σ v (x, v)) 2 ≡ 1 > 0, as in the constant coefficient case.We will also use the same initial condition as in (3.2).
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 ∆ t = 10 −4 , since an exact solution is not available, and again illustrated as orange for m2 and red for m3.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 ∆ 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 ∆ t = 10 −5 as our reference solution and compare an Euler method with ∆ t = 10 −4 , as well as m2 and m3 to it.The results are given in Table 2.
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 All in all we come to the same conclusion as in the constant coefficient case.

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 ∆ t = 10 −5 taking approximately 21 minutes, in average per trajectory, while Magnus order 3 takes only 4.62 seconds using ∆ t = 0.01, ∆ 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.

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 where the coefficients B, A ∈ R d×d are constant matrices and W is a standard one-dimensional Brownian motion.
We make the ansatz X t = exp (Y t ), where By Itô's formula (cf.[10, p. 8 Lemma 1 and p. 9 Proposition 1]) we have we have where the operators L and Q associated to the first and second order derivative of the exponential map are given by A comparison of coefficients yields 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 L, i.e.
The Bernoulli numbers are denoted by β k and we recall that β 0 = 1, β 1 = − 1 2 , β 2 = 1 6 , β 3 = 0 and β 4 = − 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 0 t = 0 R d×d and iterate for n ≥ 1 In order to derive the Magnus expansion formulas we will introduce some bookkeeping parameters , δ > 0 and substitute A by A, as well as B by δB.Henceforth, we will denote the n-th order Picard iteration with the substitution by Y n, ,δ t .
Order 1.Let us derive the first-order Magnus expansion, meaning we are interested in all terms with the first power of and δ in the first-order Picard iteration.Thus, we insert Y 0, ,δ t ≡ 0 into (A.3).Notice, that the zero matrix commutes with all matrices and therefore by definition of ad 0 Inserting this into the formula for µ yields because A commutes with itself as well.
Since, the Itô-correction term is of order 2 it will not be part of the first-order Magnus expansion and we have 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.1) and (A.2) to disregard the majority of the high-order terms.
Let us consider (A.1) first and let N = 2 denote the desired order of the Magnus expansion.
The operator ad 0 will at least result in a first order term in or δ.The commutator ad 1 Y (X) = [Y, X] = Y X − XY will at least have second order terms, because X and Y will at least have a first-order or δ coefficient.Now, the nested commutator ad 2 Y (X) = [Y, [Y, X]] 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.1) up to N − 1.Now, let us consider (A.2).We will split our consideration into two parts by linearity and focus for the moment on In the case 0 = k = n = m we will get σ 2 , which has at least order 2 terms.Hence, we will only need k = 1, . . ., 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.1).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 : (A.5) There will still be a lot of higher-order terms in (A.4) and (A.5) but it will make our derivation of the formulas easier.Now, we compute σ 2 Y 1 t : Next, we compute µ 2 Y 1 t : Counting and δ we can already see that In total, we have 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 ∈ N 0 we have t 0 s p W q s dW s = 1 q + 1 t p W q+1 t − t 0 q(q + 1) 2 s p W q−1 s + pW q+1 s s p−1 ds (A.6) s s 0 r p 2 W q 2 r drds .

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

Figure 3 .
Figure 3.: Graphical representation of U d,κ t compared to U d t for the error analysis to disregard boundary effects.

where •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 κ × I κ by taking the mean over all trajectories.The second error is the average of the first error for the corresponding grid indicated by κ.The larger κ the further away we are from the spatial boundary, as illustrated in Figure3.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 κ into account.It will serve as our main error norm in this paper.Regarding the Lebesgue discretization ∆ Leb t we have observed that the computational times are more or less constant for ∆ Leb t ∈ 10 −5 , 10 −2 and increase for lower ∆ Leb t significantly.

3. 1 .
The Magnus expansion for the stochastic Langevin equation with constant coefficients lations 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). 0

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

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

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

Figure 9 .
Figure 9.: Variable coefficients as in (3.3): Computational times and errors compared to Euler with ∆ t = 10 −4 of the Magnus expansion for varying number step-size ∆ t with fixed spatial dimension d = 200.
(B).With exactly the same arguments as for (A.1) 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.2), and split it again in two parts, i.e. we consider (σ N −k (Y t )) (n + 1)! ad m Yt (σ N −k (Y t )) (m + 1)! + ad n Yt (σ N −k (Y t )) , ad m Yt (σ N −k (Y t )) (n + m + 2)(n + 1)!m! .
Schematic of the iterative Magnus scheme 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.
Start t 5EndFigure 2.:Iterated Magnus expansionAn 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 subintervals (0, t 1 ], (t 1 , t 2 ], ..., (t n , T ].As illustrated in Figure2, 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.

Table 1 .
: Notations for the numerical experiments.Magnus order 2 with step-size ∆ t = x M3, x Magnus order 3 with step-size ∆ t = x E, x euler with step-size ∆ t = x Figure 7.: Constant coefficients as in (3.1): Computational times and errors of the Magnus expansion and Euler scheme for d = 200.