Log In Sign Up

Improved efficiency of multilevel Monte Carlo for stochastic PDE through strong pairwise coupling

Multilevel Monte Carlo (MLMC) has become an important methodology in applied mathematics for reducing the computational cost of weak approximations. For many problems, it is well-known that strong pairwise coupling of numerical solutions in the multilevel hierarchy is needed to obtain efficiency gains. In this work, we show that strong pairwise coupling indeed is also important when (MLMC) is applied to stochastic partial differential equations (SPDE) of reaction-diffusion type, as it can improve the rate of convergence and thus improve tractability. For the (MLMC) method with strong pairwise coupling that was developed and studied numerically on filtering problems in [Chernov et al., Numer. Math., 147 (2021), 71-125], we prove that the rate of computational efficiency is higher than for existing methods. We also provide numerical comparisons with alternative coupling ideas illustrate the importance of this feature. The comparisons are conducted on a range of SPDE, which include a linear SPDE, a stochastic reaction-diffusion equation, and stochastic Allen–Cahn equation.


page 1

page 2

page 3

page 4


MG/OPT and MLMC for Robust Optimization of PDEs

An algorithm is proposed to solve robust control problems constrained by...

Multilevel Quasi-Monte Carlo for Optimization under Uncertainty

This paper considers the problem of optimizing the average tracking erro...

Pairwise coupling of convolutional neural networks for better explicability of classification systems

We examine several aspects of explicability of a classification system b...

O(N^2) Universal Antisymmetry in Fermionic Neural Networks

Fermionic neural network (FermiNet) is a recently proposed wavefunction ...

Monte Carlo Sort for unreliable human comparisons

Algorithms which sort lists of real numbers into ascending order have be...

Central Limit Theorems for Coupled Particle Filters

In this article we prove a new central limit theorem (CLT) for coupled p...

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 variance-reduction 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 finite-difference 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 time-discretization of reaction-diffusion type SPDE. In the work of [7] the new form of coupling is presented, however to re-emphasis 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 standard-format of MLMC, which is a cost-versus-error 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 that

For every , we introduce the finite-dimensional 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


where is an unbounded linear operator, is a possibly random-valued initial condition, is a potentially nonlinear reaction term and is a so-called -Wiener process, cf. (A.1).

The mild solution to the (2.1) is an -valued predictable processes satisfying


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 negative-definite and spectrally decomposable in the considered basis:

and that the -Wiener process takes the form

where is a sequence of independent scalar-valued Wiener processes. The strictly positive sequence and the non-negative 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 Galerkin-based solvers.

2.3.1. Continuous-time Galerkin methods

For , the continuous-time Galerkin solution is specified by requiring


where , and It is well-known [25] that (2.3) has a unique mild solution given by

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


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


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:

Theorem 2.1 (Jentzen and Kloeden [24]).

If the assumptions in Appendix A hold, then


where is the mild solution (2.2) of (2.1), and are the exponential Euler approximations (2.4).

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 component-level, the scheme can be written


Extrapolating from the convergence result in [36]

for the white-noise setting and our later numerical studies, we conjecture that the Crank–Nicolson method converges strongly with the following rates:

Conjecture 2.2 (Stochastic Crank–Nicolson).

Let a CFL constraint on the form hold. Then under sufficient regularity, it holds for any (cf. Assumption A.2) that

where denotes the mild solution to the SPDE (2.1) and the Crank–Nicolson approximations.

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 , variance-reduction 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 super-indices 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



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.



where . Then, we use that and independence, to have

Finally, one may verify that the RHS of satisfies (2.6) ∎

Remark 2.4.

The theorem applies in settings where one replaces in Theorem 2.3 (ii) by , and for some problems this may improve the rate . Practically, however, there is often little to gain by such an approach, as weak estimates like can be a lot more intractable than strong ones, cf. [27].

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 reaction-term setting in [7]. The -th component of two iterations of the exponential Euler method on resolution level at time can be written




for . The coupled ‘coarse solution is solved with time-step one iteration at time , can be written


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

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 so-called 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


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 so-called 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:



for . One iteration of the pairwise coupled coarse solution at time , is given by



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.

Consider the SPDE (2.1) for a linear operator with . If the assumptions in Section A hold for some , then the correctly coupled exponential Euler MLMC method with


  • .

  • .

  • .

And for any , it holds for any sufficiently small and that we can find sequence such that




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 number-of-samples-per-level sequence


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 reaction-diffusion 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.

Remark 3.4.

Under Conjecture 2.2 and the same assumptions as in Proposition 3.2, with , the computational cost of achieving an MSE of for the Crank–Nicolson based MLMC method is bounded by

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 reaction-term 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 Allen-Cahn 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
Table 1. Different reaction terms tested in the numerical experiments.

Method parameters

Exponential Euler, correct coupling

As motivated from Proposition 3.2 when , we set

Neglecting logarithmic factors, Proposition 3.2 implies the MLMC rate triplet . The sequence is determined by formula (3.7) with .

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.


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.

Figure 1. Comparison of the couplings at the finest level (blue) and the coarsest level (red) with linear reaction term . The spatial resolution is fixed to in all realizations and the temporal coupling is and from top to bottom row.

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.

Figure 2. Strong errors with respect to the temporal (left) and spatial (right) couplings, for the linear reaction term .
Figure 3. Approximation of for different values of in colored full lines and the reference solution (dashed line). Left: correct coupling. Middle: naive coupling. Right: Crank–Nicolson coupling.
Figure 4. Left: compared against the tolerance . Right: cost compared against the tolerance . Both plots are for the linear reaction term .

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 pseudo-reference 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

Figure 5. Comparison of the couplings at the finest level (blue) and the coarsest level (red) for the trigonometric reaction. The spatial resolution is fixed to and the temporal coupling is and from top to bottom row.

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 6. Strong errors with respect to temporal and spatial couplings, for the sinusoidal reaction term .

Figure 7 shows the approximation of by the MLMC methods for different inputs