1. Introduction
The efficiency of numerical methods is a very important topic for practitioners that has seen a surge of interest lately in the field of uncertainty quantification (UQ) [33, 34, 35]
. UQ seeks to combine statistical and probabilistic techniques, with traditional numerical schemes, to improve the modeling and accuracy of such numerical methods. Examples of applications include climate modeling, subsurface flow, medical imaging, more recently, deep learning
[1, 11, 30]. A particular focus has been given on the class of numerical methods known as Monte Carlo (MC) methods, which are used to solve problems incorporating elements of randomness or uncertainty [28, 32, 35], i.e. stochastic computation. One methodology which has exhibited improved efficiency and a high level of applicability, is the of multilevel Monte Carlo (MLMC).MLMC is a computational technique aimed at reducing the computational cost of Monte Carlo. Specifically it reduces the cost of attaining a mean squared error (MSE) of order , for . The methodology was first introduced by Heinrich [18] and extended and popularized by various works on diffusion processes by Giles [13, 14]
. The methodology of MLMC can be viewed as a variancereduction techniques. Since these works, MLMC has been applied in numerous areas, including stochastic filtering
[5, 20, 21, 12], Markov chain Monte Carlo (MCMC)
[6, 9] and partial differential equations with random input arising in UQ [2, 6, 17]. Despite the amount of advancements made with MLMC, the number of applications and research papers on applying the methodology to stochastic partial differential equations (SPDE) [23]. Such examples include a finitedifference solvers with applications in mathematical finance [15], and finite element methods for parabolic SPDE [3, 4]. Therefore there are open questions on how one successfully can apply MLMC to SPDE in a more general sense, which we use as our motivation for this work. Our aim from this manuscript is to present a complexity study of an alternative way to apply MLMC for SPDE, which can demonstrate computational gains. This alternative way is based on the exponential Euler method [10, 19, 24, 29] and strong pairwise coupling of solution relizations on different levels. The strong pairwise coupling approach was introduced and studied experimentally in the work on filtering methods for SPDE by Chernov et al. [7], and it is based on the exponential Euler integrator [10, 19, 26, 29] for timediscretization of reactiondiffusion type SPDE. In the work of [7] the new form of coupling is presented, however to reemphasis this has only been studied and analyzed in a filtering context.In this work our highlighting contribution is that we demonstrate the improvements of the discussed coupling approach, for numerically solving SPDE. This is presented in the standardformat of MLMC, which is a costversuserror result, of the MSE, presented in Proposition 3.2. Specifically, our findings suggest that in order to attain an MSE of in a standard setting, requires a computational cost of . This is a reduction in cost compared to other existing work, for which the cost is , cf. [4]. We also numerically verify these computational gains on a collection of numerical SPDE, which include a linear SPDE and the stochastic reaction–diffusion and Allen–Cahn equation. To demonstrate the superiority of this method over alternative ones, we compare our coupling to two others, entitled a ‘naive’ coupling and a coupling based on the Crank–Nicolson scheme.
The outline of this paper is as follows. In Section 2 we introduce and review the concepts of MLMC, while introducing our model problem which is a semilinear stochastic partial differential equation. Section 3 describes our proposed coupling method for MLMC and two alternative methods. We also summarize the theoretical properties of our main MLMC method in Proposition 3.2. Numerical experiments on various SPDE are conduced in Section 4 to demonstrate the improvement with the proposed coupling. Finally, we conclude our findings, and provide future areas of research, in Section 5. Required model assumptions are provided in the Appendix.
2. Background material
In this section we present and review the idea of (MLMC) method, which acts as the primary focus of this article. This will include a motivation behind the methodology and a review of the vanilla (MLMC) method. This will lead onto a discussion of spectral Galerkin discretizations of stochastic partial differential equations, whilst introducing the notion of Euler exponential integrators.
2.1. Notation
Let and let
be a complete probability space, where
is a probability measure on the measurable space , with a normal filtration , for . denotes a separable Hilbert space with inner product , norm and orthogonal basis . denotes the associated Bochner–Hilbert space, consisting of the set of strongly measurable maps such thatFor every , we introduce the finitedimensional subspace and the associated orthogonal projection operator for . For a (normally implicitly given) set and mappings , the notation implies there exists a such that for , and the notation means that and . For with , we set , and for we define .
2.2. Problem setup
Here, we consider a semilinear stochastic partial differential equation (SPDE) [8] of the form
(2.1) 
where is an unbounded linear operator, is a possibly randomvalued initial condition, is a potentially nonlinear reaction term and is a socalled Wiener process, cf. (A.1).
The mild solution to the (2.1) is an valued predictable processes satisfying
(2.2) 
The general form of (2.1) encapsulates numerous SPDE in practice. We will introduce and numerically study some of these in Section 4. Before introducing our underlying numerical scheme, we require a number of assumptions which we have deferred to the appendix. Suffice it to say here that we do assume that the linear operator is negativedefinite and spectrally decomposable in the considered basis:
and that the Wiener process takes the form
where is a sequence of independent scalarvalued Wiener processes. The strictly positive sequence and the nonnegative sequence are further described in the appendix.
2.3. Numerical methods
Numerical solutions of SPDE has traditionally been executed through the use of finite difference methods and finite element methods (FEM) [23, 28, 36]. For the relevance of this work, we will utilize and discuss an alternative class of Galerkinbased solvers.
2.3.1. Continuoustime Galerkin methods
2.3.2. The exponential Euler method
For a given , let and be the nodes of a uniform partition of defined by for . Then, for given , the exponential Euler approximations of are defined by
(2.4) 
where denotes the inverse operator of . Defining the components of and by , for , and recalling the spectral decomposition of the operator , we arrive at the recursive relations below
where
such that
is an eigenvalue of the operator
, presented in Assumption A.1 in the appendix. Convergence properties of the exponential Euler scheme has been studied in [24], where they demonstrate convergence in the strong sense, which highlights an improvement of the order of convergence in time against traditional numerical schemes:2.3.3. The Crank–Nicolson method
Given the discretization parameters in space and time introduced above, the Crank–Nicolson approximations of are given by and
for . On the componentlevel, the scheme can be written
where
Extrapolating from the convergence result in [36]
for the whitenoise setting and our later numerical studies, we conjecture that the Crank–Nicolson method converges strongly with the following rates:
2.4. The vanilla Multilevel Monte Carlo method
The Monte Carlo methods consist of a class of algorithms for computing/approximating expectations of random variables. They have been heavily exploited in the areas of computational chemistry, uncertainty quantification, mathematical finance and statistics, but and more recently, various extensions of the methodology has received a lot of attention in numerical analysis.
For an valued random variable
, we consider the MC estimator
where the samples are independently drawn random variables. When it is computationally costly to sample , variancereduction techniques may improve the efficiency through reducing the statistical error of the estimator. The multilevel Monte Carlo (MLMC) method is an alternative estimator based on building a hierarchy of progressively refined simulations on discretization levels . The MLMC estimator is defined by
where and the superindices in denotes the th sample on the th discretization level of . Moreover, and are pairwise coupled by sharing the same driving noise, and all draws on all levels are independent, meaning that all random variables in the sequence are independent. Pairwise coupling on this form will often reduce the variance of the estimator.
One way to assess the performance of MC methods is through the mean squared error (MSE). The following theorem describes the cost versus error of the MLMC methodology for valued random variables:
Theorem 2.3.
Suppose that there exist constants such that and

,

,

.
Then for any and , there exists a sequence such that
and
(2.6) 
The proof of this result is a straightforward extension of the original theorem presented by Giles [13] for weak approximations of stochastic differential equations.
Sketch of proof.
Let
(2.7) 
where . Then, we use that and independence, to have
Finally, one may verify that the RHS of satisfies (2.6) ∎
3. Multilevel Monte Carlo methods for SPDE
In this section we describe three different forms of pairwise coupling of levels in MLMC and present theoretical results on the performance of our coupling of choice. We will employ the following notation for the multilevel hierarchy of discretized solutions: denotes the numerical approximation of a given method on resolution level at time for , where , and , for some given and .
3.1. Correct coupling of exponential Euler
This MLMC method was first introduced and analyzed for the linear reactionterm setting in [7]. The th component of two iterations of the exponential Euler method on resolution level at time can be written
and
where
(3.1) 
for . The coupled ‘coarse solution’ is solved with timestep one iteration at time , can be written
where
The pairwise coupling over one timestep is obtained through coupling induced by the driving noise in . That is, by (3.1) we have that
which yields
(3.2) 
Remark 3.1.
We will refer to (3.2) as the ‘correct coupling’ for the exponential Euler method. To motivate the adjective ‘correct’, on top of (3.2) being an equality, note that when the socalled correct coupling ensures exact coupling between the solutions in the strongest possible sense:
But when is a nonlinear mapping it is not clear if correct coupling’ is the best pairwise coupling approach for the (MLMC) method, as it then may be that a coupling that is more sensitive to the dynamics contribution from the reaction term would lead to more variance reduction. To the best of our knowledge, it is an open problem to describe good coupling approaches for valued stochastic processes in general: couplings that satisfy the telescoping sum property
(3.3) 
and maximize the convergence rate of the strong error (see Theorem 2.3).
3.1.1. Naive coupling
For showcasing the importance of strong pairwise coupling, we also propose the alternative socalled naive coupling for the exponential Euler method:
Superficially, this may look like a natural extension of the standard approach to pairwise coupling of Wiener processes in stochastic differential equations. But this coupling does not satisfy (3.2), and our numerical tests indicate that when is nonlinear the naive coupling does not satisfy the telescoping sum property (3.3) (leading to a biased MLMC estimator) and the strong observed convergence rate is lower than when using correct coupling.
3.2. Crank–Nicolson
The last coupling we consider is for the Crank–Nicolson (CN) method. Two iterations of the CN method on resolution level at time can be written as follows for the th component:
and
where
for . One iteration of the pairwise coupled coarse solution at time , is given by
where
(3.4) 
for . We apply the coupling (3.4) for the MLMC implementation of the Crank–Nicolson method. This may be viewed as the “correct coupling” for the Crank–Nicolson method, in the sense that (3.4) is an equality for the driving noise on coupled levels. Furthermore, the resulting method satisfies the telescoping sum property (3.3).
3.3. MLMC for SPDE
In this section, we present cost versus error results for the MLMC method using the Exponential Euler method with correct coupling.
Proposition 3.2.
Proof.
A pair of correctly coupled solutions and can individually be viewed as plain vanilla exponential Euler solutions using the same driving Wiener process on different resolutions. Consequently,
This verifies the rate (ii), and the rate (i) follows by applying Jensen’s inequality to (ii). Introducing the following numberofsamplesperlevel sequence
(3.7) 
the results (3.5) and (3.6) can be proved by a similar argument as for Theorem 2.3. (The exponent in the MSE enters as is replaced by the slightly too fast decay in the formula for and is set slightly too low, in the sense that it only ensures a bias error of .) ∎
Remark 3.3.
A general framework for (MLMC) methods for reactiondiffusion type SPDE in the setting of and for numerical methods with a strong convergence rate was first developed in [4]. When and , the MSE was achieved at the computational cost for that method in [4, Theorem 4.4] compared to a cost for our method.
4. Numerical examples
In this section we numerically test our correctly coupled exponential Euler MLMC method against the two other coupling ideas described in Section 3. We study a linear parabolic equation, a sinusoidal reactionterm equation, as well as a stochastic Allen–Cahn equation. To connect our numerical results to theory, note that Proposition 3.2 applies to the linear problem, and it does not apply to the AllenCahn equation, since its reaction term does not satisfy Assumption A.3. However it is an open question whether the sinusoidal reaction term satisfies Assumption A.3, and if the aforementioned corollary applies to this equation.
The SPDE we consider here has the general form
where with Fourier basis functions for , and the final time is set to . We consider the linear operator
with eigenvalues
The Wiener process (A.1) is in the setting of
with . We consider the three different reaction terms presented in Table 1, all belonging to the class of Nemytskii operators, cf. [28].
Reaction term  

Linear  
Trigonometric  
Allen–Cahn 
Method parameters
Exponential Euler, correct coupling
Exponential Euler, naive coupling
Motivated from numerical tests of convergence rates, we set
Neglecting logarithmic factors, we observe the rates in the numerical tests (where is conservatively estimated to ). The sequence is determined by formula (2.7) with and setting . Note however that this MLMC method does not generally converge as it does not satisfy the telescoping sum property.
Crank–Nicolson
We set
Neglecting logarithmic factors we observe the rates in numerical experiments (where is conservatively estimated to ). This is consistent with Conjecture 2.2. The sequence is determined by formula (2.7) with and setting .
Remark 4.1.
When the reaction term is nonlinear and only known/easy to evaluate pointwise, e.g., , then the spectral representation of has to be approximated. For the last two of the considered SPDE we do these approximation steps using 1. IFFT to evaluate on a uniform mesh of values over , and 2. FFT of to obtain a spectral representation of the reaction term. See [7, Section 6.3.1] and [28] for further details.
4.1. Linear reaction
We first consider the SPDE with . Figure 1 presents pairwise coupled realizations for the different methods on progressively finer temporal resolutions. We observe that the “correct coupling” couples the realization pair well at all resolution levels, while the same does not hold for the other methods.
The strong convergence rate, in both time and space for each of the methods is presented in Figure 2. The convergence rate in the temporal coupling, on the left, is estimated by
with a fixed spatial resolution and independent samples of the random variable in the Monte Carlo estimator. This error describes to the convergence rate in . The correctly coupled method attains a the rate , while the other methods the rate .
The spatial coupling, on the right, is estimated by
with a fixed temporal resolution , varying the spatial resolution by the sequences of the respective methods, and using independent samples. This error describes the convergence rate in , which we observe is for all methods.
In Figure 3 we provide the approximation of by the MLMC methods for different inputs for . The correctly coupled MLMC estimator (left) has less fluctuations and converges faster than the other methods.
Figure 4 presents the MLMC approximation error versus tolerance and the computational cost versus tolerance for different input tolerances . For the former, we notice as the tolerance increase, so does the approximation error for both the naive and CN coupling. This differs to the correct coupling which seems to plateau. For the later we can see the correct coupling indeed attains the input MSE at a cost close to . Noting that naive coupling method also converges, but at a slower rate, this example illustrates the importance of strong coupling.
4.2. Trigonometric reaction
We next consider the SPDE with
A pseudoreference solution is computed by with the accuracy parameter is used in all convergence studies, as we do not know the explicit solution to this problem.
Firstly Figure 5 presents pairwise coupled realization for the different methods on progressively finer temporal resolutions. Similar to the first example, we observe that the “correct coupling” couples the realization pairs better than the other methods
As before Figure 6 shows the strong convergence rate, in terms of temporal and spatial coupling for all methods. The temporal coupling is once again for the correctly coupled method and for the other two. The spatial rate is for all methods.
Figure 7 shows the approximation of by the MLMC methods for different inputs