A multilevel Monte Carlo method for asymptotic-preserving particle schemes in the diffusive limit

TL;DR

This paper introduces a multilevel Monte Carlo method for particle schemes in kinetic equations that effectively reduces bias and computational cost in the diffusive limit, improving long-time particle simulations.

Contribution

The paper develops a multilevel Monte Carlo approach that combines different time step estimates to reduce bias and computational effort in asymptotic-preserving particle schemes.

Findings

Significant bias reduction with larger time steps.

Reduced computational cost compared to classical Monte Carlo.

Effective correlation of trajectories across levels.

Abstract

Kinetic equations model distributions of particles in position-velocity phase space. Often, one is interested in studying the long-time behavior of particles in high-collisional regimes in which an approximate (advection)-diffusion model holds. In this paper we consider the diffusive scaling. Classical particle-based techniques suffer from a strict time-step restriction in this limit, to maintain stability. Asymptotic-preserving schemes avoid this problem, but introduce an additional time discretization error, possibly resulting in an unacceptably large bias for larger time steps. Here, we present and analyze a multilevel Monte Carlo scheme that reduces this bias by combining estimates using a hierarchy of different time step sizes. We demonstrate how to correlate trajectories from this scheme, using different time steps. We also present a strategy for selecting the levels in the multilevel scheme. Our approach significantly reduces the computation required to perform accurate simulations of the considered kinetic equations, compared to classical Monte Carlo approaches.

Tables

Table 1. Table 4: Cost comparison between classical and multilevel Monte Carlo

RMSE	Classical cost	Multilevel cost	Multilevel speedup
0.1	4 544	8 062	0.56
0.01	37 011 456	7 947 587	4.66
0.001	7 722 983 424	1 425 800 335	5.42

Table 2. Table 8: Cost comparison with and without an extra coarse level

RMSE	Without coarse level	With coarse level	Coarse level speedup
0.1	8 062	2 467	3.27
0.01	7 947 587	7 616 035	1.04
0.001	1 425 800 335	1 324 490 488	1.08

Equations

$${\partial }_{t}f(x,v,t)+v\cdot {\nabla }_{x}f(x,v,t)=Q\left(f(x,v,t)\right),$$

$$Q\left(f(x,v,t)\right)=\mathcal{M}(v)\rho (x,t)-f(x,v,t),$$

$$\rho (x,t)={\int }_{{\mathcal{D}}_v}f(x,v,t)\text{d}v.$$

$$\epsilon {\partial }_{t}f(x,v,t)+v\cdot {\nabla }_{x}f(x,v,t)=\frac{1}{\epsilon }\left(\mathcal{M}(v)\rho (x,t)-f(x,v,t)\right).$$

$${\partial }_t\rho (x,t)={\Delta }_{xx}\rho (x,t).$$

$${\partial }_{t}f(x,v,t)+\frac{v}{\epsilon }{\partial }_{x}f(x,v,t)=\frac{1}{{\epsilon }^2}\left(\mathcal{M}(v)\rho (x,t)-f(x,v,t)\right),$$

$${\left\{\left(X_{p,\Delta t}^n,V_{p,\Delta t}^n\right)\right\}}_{p=1}^P.$$

$$X_{p,\Delta t}^{n+1}=X_{p,\Delta t}^n+\Delta t{V}_{p,\Delta t}^n.$$

$$V_{p,\Delta t}^{n+1}=\{\begin{array}{ll}{V}_{p,\Delta t}^{n,\ast }\sim \frac{1}{\epsilon }\mathcal{M}(v),& \text{with probability }{p}_{c,\Delta t}=\Delta t/{\epsilon }^2,\\ V_{p,\Delta t}^n,& \text{otherwise.}\end{array}$$

$$\mathcal{M}(v)=\frac{1}{\tilde{v}}\mathcal{B}\left(\frac{v}{\tilde{v}}\right)=\frac{1}{\tilde{v}}\mathcal{B}(\bar{v}),$$

$$E\left[\bar{V}\right]=0\text{and}V\left[\bar{V}\right]=1.$$

$$\{\begin{array}{l}{\partial }_t{f}_{+}(x,t)+\frac{1}{\epsilon }{\partial }_x{f}_{+}(x,t)=\frac{1}{{\epsilon }^2}\left(\frac{\rho (x,t)}{2}-f_{+}(x,t)\right)\\ {\partial }_t{f}_{-}(x,t)-\frac{1}{\epsilon }{\partial }_x{f}_{-}(x,t)=\frac{1}{{\epsilon }^2}\left(\frac{\rho (x,t)}{2}-f_{-}(x,t)\right)\end{array},$$

$${\partial }_{t}f+\frac{\epsilon v}{{\epsilon }^2+\Delta t}{\partial }_{x}f=\frac{{\tilde{v}}^2\Delta t}{{\epsilon }^2+\Delta t}{\partial }_{xx}f+\frac{1}{{\epsilon }^2+\Delta t}\left(\mathcal{M}(v)\rho -f\right),$$

$$X_{p,\Delta t}^{n+1}=X_{p,\Delta t}^n+V_{p,\Delta t}^n\Delta t+\sqrt{2\Delta t}\sqrt{D_{\Delta t}}{\xi }_p^n,V_p^n\sim {\mathcal{M}}_{\Delta t}(v),$$

$${\mathcal{M}}_{\Delta t}(v)=\frac{\epsilon }{{\epsilon }^2+\Delta t}\mathcal{M}(v),D_{\Delta t}=\frac{{\tilde{v}}^2\Delta t}{{\epsilon }^2+\Delta t}.$$

$$V_{p,\Delta t}^{n+1}=\{\begin{array}{ll}{V}_{p,\Delta t}^{n,\ast }\sim {\mathcal{M}}_{\Delta t}(v),& \text{with probability }{p}_{c,\Delta t}=\frac{\Delta t}{{\epsilon }^2+\Delta t},\\ V_{p,\Delta t}^n,& \text{otherwise.}\end{array}$$

$$Y(t^{\ast })=E\left[F\left(X(t^{\ast }),V(t^{\ast })\right)\right]={\int }_{{\mathcal{D}}_v}{\int }_{{\mathcal{D}}_x}F(x,v)f(x,v,t^{\ast })\text{d}x\text{d}v.$$

$$\hat{Y}(t^{\ast })=\frac{1}{P}\sum _{p=1}^{P}F\left(X_{p,\Delta t}^N,V_{p,\Delta t}^N\right),t^{\ast }=N\Delta t,$$

$$V\left[\hat{Y}(t^{\ast })\right]=\frac{1}{P}V\left[F\left(X_{p,\Delta t}^N,V_{p,\Delta t}^N\right)\right],$$

$${\hat{Y}}_0(t^{\ast })=\frac{1}{P_0}\sum _{p=1}^{P_0}F\left(X_{p,\Delta t_0}^{N_0},V_{p,\Delta t_0}^{N_0}\right),t^{\ast }=N_0\Delta t_0.$$

$${\hat{Y}}_{\ell }(t^{\ast })=\frac{1}{P_{\ell }}\sum _{p=1}^{P_{\ell }}\left(F\left(X_{p,\Delta t_{\ell }}^{N_{\ell }},V_{p,\Delta t_{\ell }}^{N_{\ell }}\right)-F\left(X_{p,\Delta t_{\ell -1}}^{N_{\ell -1}},V_{p,\Delta t_{\ell -1}}^{N_{\ell -1}}\right)\right),t^{\ast }=N_{\ell }\Delta t_{\ell }.$$

$$V\left[F\left(X_{p,\Delta t_{\ell }}^{N_{\ell }},V_{p,\Delta t_{\ell }}^{N_{\ell }}\right)-F\left(X_{p,\Delta t_{\ell -1}}^{N_{\ell -1}},V_{p,\Delta t_{\ell -1}}^{N_{\ell -1}}\right)\right]$$

$$V\left[F\left(X_{p,\Delta t_{\ell }}^{N_{\ell }},V_{p,\Delta t_{\ell }}^{N_{\ell }}\right)\right]+V\left[F\left(X_{p,\Delta t_{\ell -1}}^{N_{\ell -1}},V_{p,\Delta t_{\ell -1}}^{N_{\ell -1}}\right)\right]$$

$$-2\mathrm{Cov}\left(F\left(X_{p,\Delta t_{\ell }}^{N_{\ell }},V_{p,\Delta t_{\ell }}^{N_{\ell }}\right),F\left(X_{p,\Delta t_{\ell -1}}^{N_{\ell -1}},V_{p,\Delta t_{\ell -1}}^{N_{\ell -1}}\right)\right).$$

$$V\left[F\left(X_{p,\Delta t_{\ell }}^{N_{\ell }},V_{p,\Delta t_{\ell }}^{N_{\ell }}\right)-F\left(X_{p,\Delta t_{\ell -1}}^{N_{\ell -1}},V_{p,\Delta t_{\ell -1}}^{N_{\ell -1}}\right)\right]\le V\left[F\left(X_{p,\Delta t_{\ell }}^{N_{\ell }},V_{p,\Delta t_{\ell }}^{N_{\ell }}\right)\right].$$

$$\hat{Y}(t^{\ast })=\sum _{\ell =0}^L{\hat{Y}}_{\ell }(t^{\ast }).$$

$$E\left[C\right]\le \{\begin{array}{ll}{c}_4{E}^{-2},& \beta >\gamma ,\\ c_4{E}^{-2}{\log }^{2}E,& \beta =\gamma ,\\ c_4{E}^{-2-\frac{(\gamma -\beta )}{\alpha }},& \beta <\gamma .\end{array}$$

$$\{\begin{array}{l}{X}_{p,\Delta t_{\ell -1}}^{n+1}=X_{p,\Delta t_{\ell -1}}^n+\Delta t_{\ell -1}{V}_{p,\Delta t_{\ell -1}}^n+\sqrt{2\Delta t_{\ell -1}}\sqrt{D_{\Delta t_{\ell -1}}}{\xi }_{p,\ell -1}^n\\ X_{p,\Delta t_{\ell }}^{n+1,0}=X_{p,\Delta t_{\ell }}^{n,0}+\sum _{m=0}^{M-1}\left(\Delta t_{\ell }{V}_{p,\Delta t_{\ell }}^{n,m}+\sqrt{2\Delta t_{\ell }}\sqrt{D_{\Delta t_{\ell }}}{\xi }_{p,\ell }^{n,m}\right)\end{array},$$

$$\Delta T_{p,\ell -1}^n=\Delta t_{\ell -1}{V}_{p,\Delta t_{\ell -1}}^n,\Delta T_{p,\ell }^{n,m}=\Delta t_{\ell }{V}_{p,\Delta t_{\ell }}^{n,m},$$

$$\Delta W_{p,\ell -1}^n=\sqrt{2\Delta t_{\ell -1}}\sqrt{D_{\ell -1}}{\xi }_{p,\ell -1}^n,\Delta W_{p,\ell }^{n,m}=\sqrt{2\Delta t_{\ell }}\sqrt{D_{\ell }}{\xi }_{p,\ell }^{n,m},$$

$$\{\begin{array}{l}{X}_{p,\Delta t_{\ell -1}}^{n+1}=X_{p,\Delta t_{\ell -1}}^n+\Delta T_{p,\ell -1}^n+\Delta W_{p,\ell -1}^n\\ X_{p,\Delta t_{\ell }}^{n+1,0}=X_{p,\Delta t_{\ell }}^{n,0}+\sum _{m=0}^{M-1}\left(\Delta T_{p,\ell }^{n,m}+\Delta W_{p,\ell }^{n,m}\right)\end{array},$$

Full text

1 Introduction

In many application domains, one encounters kinetic equations, modelling particle behavior in a position-velocity phase space. Examples are plasma physics Birdsall2004 , bacterial chemotaxis Rousset2011c and computational fluid dynamics Pope1981 . Many of these domains exhibit a strong time-scale separation, leading to an unacceptably high simulation cost Cercignani1988 . Typically, one is interested in some macroscopic quantities of interest, e.g., some moments of the particle distribution, which are computed as averages over velocity space. The time-scale at which these quantities of interest change is often much slower than that governing the particle dynamics, making these models very stiff problems: a naive simulation requires both small time steps to capture the fast dynamics, and long time horizons to capture the evolution of the macroscopic quantities of interest. The exact nature of the macroscopic behavior depends on the problem scaling, which can be hyperbolic or diffusive Dimarco2014 .

We are interested in a $d$-dimensional kinetic equation of the form

[TABLE]

where $f(x,v,t)$ represents the distribution of particles as a function of space $x\in\mathcal{D}_{x}\subset\mathbb{R}^{d}$, velocity $v\in\mathcal{D}_{v}\subset\mathbb{R}^{d}$ and time $t\in\mathbb{R}^{+}$ and $Q(f(x,v,t))$ is a collision operator, resulting in discontinuous velocity changes. In this paper, we consider the BGK operator Bhatnagar1954 , which linearly drives the velocity distribution to a steady state distribution $\mathcal{M}(v)$. In this case, the operator $Q\left(f(x,v,t)\right)$ is written as

[TABLE]

in which we have introduced the position density

[TABLE]

With the collision operator (2), individual particles follow a velocity-jump process.

To make the time-scale separation explicit, we consider a dimensionless, diffusively scaled version of (1). We introduce a parameter $\varepsilon$, representing the mean free path. When $\varepsilon$ decreases, the average time between collisions decreases. In the diffusive scaling, we factor out the fast collision time scale by writing the right hand side as $(1/\varepsilon)\;Q$, while simultaneously re-scaling time by $1/\varepsilon$:

[TABLE]

Taking the diffusion limit $\varepsilon\to 0$, we drive the rate of collisions to infinity, while simultaneously increasing the slow time-scale. It can be shown that, in the limit $\varepsilon\to 0$, the particle density resulting from (4) converges to the diffusion equation Lapeyre2003

[TABLE]

Equation (4) can be simulated using a broad selection of methods. Deterministic methods solve the kinetic equation (4) for $f(x,v,t)$ on a grid (using, for instance, finite differences or finite volumes), giving the particle distribution in the position-velocity phase space. This approach quickly becomes computationally infeasible as the dimension grows, as a grid must be formed over the $2d$-dimensional domain, $\mathcal{D}_{x}\times\mathcal{D}_{v}$. Stochastic methods, on the other hand, perform simulations of individual particle trajectories, with each trajectory representing a sample of the probability distribution $f(x,v,t)$. These methods do not suffer from the curse of dimensionality, but introduce a statistical error in the computed solution. When using explicit time steps, both approaches become prohibitively expensive for small values of $\varepsilon$ due to the time-scale separation.

To avoid the issues caused by time-scale separation, one can use asymptotic-preserving schemes. Such methods preserve the macroscopic limit equation, in our case (5), as $\varepsilon$ tends to zero, but do not suffer from the time step constraints caused by the time-scale separation. A large number of such methods have already been developed in literature for deterministic methods in the diffusive limit. For a non-exhaustive list, we refer to Bennoune2008 ; Boscarino2013 ; Buet2007 ; Crouseilles2011 ; Dimarco2012 ; Gosse2002 ; Jin1999 ; Jin1998 ; Jin2000 ; Klar1998 ; Klar1999 ; Larsen1974 ; Lemou2008 ; Naldi2000 and to a recent review paper Dimarco2014 , which gives an overview of the current state of the art concerning these methods. In the particle setting, only a few asymptotic-preserving methods exist, mostly in the hyperbolic scaling Degond2011 ; Dimarco2008 ; Dimarco2010 ; Pareschi1999 ; Pareschi2001 ; Pareschi2005 . In the diffusive scaling, we are only aware of three works Crestetto2018 ; Dimarco2018 ; Mortier2019 . In this paper, we make use of the scheme proposed in Dimarco2018 , where operator splitting was successfully applied to a modified kinetic equation, resulting in an unconditionally stable fixed time step particle method. This stability comes at the cost of an extra bias in the model, proportional to the size of the time step (see Section 2.2 for details).

In this paper, we present and analyse an approach to eliminate the bias introduced by using the asymptotic-preserving schemes with large time steps, through the use of the multilevel Monte Carlo method Giles2008 . Multilevel Monte Carlo methods first compute an initial estimate, using a large number of samples with a large time step (and hence a low cost per sample). This estimate has a small variance, but is expected to have a large bias. Afterwards, the bias in this initial estimate is reduced by performing corrections using a hierarchy of simulations with increasingly smaller time steps. Under correct conditions, far fewer samples with small time steps are needed, compared to a direct Monte Carlo simulation with the smallest time step, resulting in a reduced computational cost. The multilevel Monte Carlo method was first introduced in finance Giles2008 , and has since been applied to other fields such as biochemistry Anderson2011 , data science Hoel2016 and structural engineering Blondeel2020 . Recently, the method has been applied to the simulation of large PDEs with random coefficients Cliffe2011 , as well as to optimisation on these models VanBarel2019 . The method has also recently been extended to higher dimensional parameter spaces as the multi-index Monte Carlo method Haji-Ali2016 . A preliminary description of the algorithm presented here was published in Loevbak2020 , together with partial numerical results.

The rest of this paper is structured as follows: In Section 2, we present the main ideas behind kinetic equations and the asymptotic-preserving particle scheme. We also describe the model equation that will be used in the numerical experiments. In Section 3, we give an overview of the multilevel Monte Carlo method. Section 4 contains the main algorithmic contribution of this paper: we present an algorithm for generating coupled particle trajectories using the asymptotic-preserving Monte Carlo scheme at different levels of the multilevel Monte Carlo hierarchy. Section 5 contains the corresponding numerical analysis. Here, we present some lemmas on the properties of the correlation of the coupled trajectories, as a function of the time step and the time-scale separation $\varepsilon$. We then prove convergence of the scheme, making use of these lemmas. In Section 6 we present a strategy for selecting which levels to include in the scheme and motivate this strategy with numerical results. Finally, in Section 7, we summarize the main results and discuss possible extensions.

2 Kinetic equations and asymptotic-preserving particle schemes

2.1 Model equation and particle scheme

For the sake of exposition, we limit this work to one spatial dimension. The proposed method is, however, general. We will explicitly mention where care is needed when extending the method or its analysis to higher-dimensional models. We rewrite (4) as

[TABLE]

with $x\in\mathbb{R}$, $v\in\mathbb{R}$ and $t\in\mathbb{R}^{+}$.

Equation (6) can be simulated using particle schemes with a finite time step of size $\Delta t$. Each particle has a state in the position-velocity phase space $(X,V)$ at each time step $n$, i.e., $X_{p,\Delta t}^{n}\approx X_{p}(n\Delta t)$ and $V_{p,\Delta t}^{n}\approx V_{p}(n\Delta t)$. We then represent the distribution of particles by an ensemble of $P$ particles, with indices $p\in\{1,\dots,P\}$,

[TABLE]

Classically, the ensemble (7) is simulated by operator splitting, which is first order in the time step $\Delta t$ Pareschi2005a . For (6), operator splitting results in two actions for each time step:

Transport step. Each particle’s position is updated based on its velocity

[TABLE] 2. 2.

Collision step. Between transport steps, each particle’s velocity is either left unchanged (no collision) or re-sampled from $\mathcal{M}(v)$ (collision), i.e.,

[TABLE]

To simplify the analysis in Section 5, we consider a transformation of the velocity distribution $\mathcal{M}(v)$, where it is re-written as a distribution with unit variance by re-scaling both the x-axis and y-axis by its characteristic velocity $\tilde{v}$. Mathematically, we write this transformation as

[TABLE]

with $\mathcal{B}(\bar{v})$ a probability distribution for which, given $\bar{V}\sim\mathcal{B}(\bar{v})$,

[TABLE]

A visualization of the transformation shown in (10) is shown in Figure 1.

We give two examples:

• •

Two discrete velocities. To limit the cost of our simulations, we will consider $\mathcal{M}(v)=\frac{1}{2}\left(\delta_{v,-1}+\delta_{v,1}\right)$, with $\delta$ the Kronecker delta function, i.e., $v$ can take the values $\pm 1$, with equal probability. In this case, (6) becomes

[TABLE]

with $f_{+}(x,t)$ the distribution of particles with positive velocity and $f_{-}(x,t)$ that of particles with negative velocity. The total particle density is given by $\rho(x,t)=f_{+}(x,t)+f_{-}(x,t)$. Equation (12) is known as the Goldstein-Taylor model Goldstein1951 .

• •

Normal distribution. Another common choice is $\mathcal{M}(v)=\mathcal{N}(v;0,\tilde{v}^{2})$, i.e., the normally distributed with expected value 0 and variance $\tilde{v}^{2}$.

Scheme (8)–(9) has a severe time step restriction $\Delta t=\mathcal{O}(\varepsilon^{2})$ when approaching the limit $\varepsilon\to 0$. This time step restriction will often result in unacceptably high simulation costs, despite the well-defined limit Dimarco2018 .

2.2 Asymptotic-preserving Monte Carlo scheme

In Dimarco2018 , an asymptotic-preserving scheme was proposed as a solution to the high simulation cost of (8)–(9) in the limit $\varepsilon\to 0$. This asymptotic-preserving scheme works by rewriting (6) as

[TABLE]

using an approach based on the IMEX discretization. In (13), we omit the space, velocity and time dependency of $f(x,v,t)$ and $\rho(x,t)$, for conciseness. In the limit $\varepsilon\to 0$, it can be shown that the modified equation (13) converges to the diffusion limit (5). It can also be shown that, in the limit $\Delta t\to 0$, (13) converges to the original kinetic equation (6) with a rate $\mathcal{O}(\Delta t)$, see Dimarco2018 .

Particle trajectories are now simulated as follows:

Transport-diffusion step. The position of the particle is updated based on its velocity and a Brownian increment

[TABLE]

in which we have taken $\xi_{p}^{n}\sim\mathcal{N}(0,1)$ and introduced a $\Delta t$-dependent velocity distribution $\mathcal{M}_{\Delta t}(v)$ and diffusion coefficient $D_{\Delta t}$:

[TABLE]

We define the characteristic velocity of $\mathcal{M}_{\Delta t}(v)$ as $\displaystyle\tilde{v}_{\Delta t}=\frac{\varepsilon}{\varepsilon^{2}+\Delta t}\tilde{v}.$ 2. 2.

Collision step. During collisions, each particle’s velocity is updated as:

[TABLE]

For the Goldstein-Taylor model, sampling the time-step dependent velocity distribution $\mathcal{M}_{\Delta t}(v)$ means multiplying the characteristic velocity $\tilde{v}_{\Delta t}$ with $\pm 1$ with equal probability, which satisfies (10)–(11). For more details on the scheme (14)–(16) see Dimarco2018 .

3 Multilevel Monte Carlo method

We want to calculate the value of a quantity of interest (QoI) $Y(t^{*})$, which is the integral of a function $F(x,v)$ of the particle position $X(t)$ and velocity $V(t)$ at time $t=t^{*}$ with respect to $f(x,v,t)$. In practical applications, integration is commonly performed only over the velocity space, giving a quantity of interest that is time and space dependent, e.g., a density or flux. Here, we opt to also integrate over position space, i.e.,

[TABLE]

This choice simplifies the notation in the remainder of this paper, while also decreasing the simulation cost of our test problems, facilitating an in-depth numerical analysis.

A classical Monte Carlo estimator $\hat{Y}(t^{*})$ for $Y(t^{*})$ is given by

[TABLE]

with particles $\left(X^{N}_{p,\Delta t},V^{N}_{p,\Delta t}\right)$, simulated using the time discretization (14)–(16).

Given a fixed cost budget for the estimator $\hat{Y}(t^{*})$, i.e., a maximal value for the product of the number of time steps and particle simulations $N\times P$, we have to make a trade-off. On the one hand, if we choose to perform more accurate simulations by taking a small $\Delta t$, the individual trajectories will have a small bias, but the required number of time steps $N$ will be very large. As a consequence, we can only simulate a limited number of trajectories $P$. Given that the variance of (18) is given by

[TABLE]

the estimated quantity of interest $\hat{Y}(t^{*})$ will have a large variance, due to insufficient sampling of the trajectory space. On the other hand, if we choose to simulate a large number of particles $P$, reducing the number of time steps $N$ by choosing a larger $\Delta t$, the estimator variance will be smaller, but the simulation bias will be larger due to the time discretization error.

The core idea behind the multilevel Monte Carlo method (MLMC) Giles2008 ; Heinrich2001 is to avoid this trade-off by combining estimators based on trajectories with different time step sizes. This is done by starting with a coarse time step size $\Delta t_{0}$. At this coarse level, we can cheaply simulate a large number of trajectories $P_{0}$, as the number of required time steps $N_{0}$ to reach the end time $t^{*}$ is small. This estimator has a large bias, but a low variance, and is given by

[TABLE]

The estimator $\hat{Y}_{0}(t^{*})$ is then refined upon by a sequence of $L$ difference estimators at levels $\ell=1,\ldots,L$. Each difference estimator uses an ensemble of $P_{\ell}$ particle pairs

[TABLE]

Each particle pair consists of two coupled particles: a particle with a fine time step $\Delta t_{\ell}$ and a particle with a coarse time step $\Delta t_{\ell-1}=M\Delta t_{\ell}$, with $M$ a positive integer. The variance of one such particle pair

[TABLE]

can be decomposed as

[TABLE]

If we ensure that $F\left(X^{N_{\ell}}_{p,\Delta t_{\ell}},V^{N_{\ell}}_{p,\Delta t_{\ell}}\right)$ and $F\left(X^{N_{\ell-1}}_{p,\Delta t_{\ell-1}},V^{N_{\ell-1}}_{p,\Delta t_{\ell-1}}\right)$ are sufficiently correlated, then we expect that

[TABLE]

This means that fewer particle simulations are needed in order to achieve a given variance when estimating (21), as compared to a single level estimator using a time step size $\Delta t_{\ell}$. We achieve correlation in the sampled quantity of interest by making each particle pair undergo correlated simulations, which intuitively can be understood as an attempt to let two particles follow the same trajectory for two different simulation accuracies. We give detailed explanation on how this is achieved in Section 4.

One can interpret the difference estimator as using the fine simulation to estimate the bias in the coarse simulation. Given a sequence of levels $\ell\in\{0,\dots,L\}$, with decreasing step sizes, and the corresponding estimators given by (20)–(21), the multilevel Monte Carlo estimator for the quantity of interest $Y(t^{*})$ is computed by the telescopic sum

[TABLE]

It is clear that the expected value of the estimator (26) is the same as that of (18), with the finest time step $\Delta t=\Delta t_{L}$. Given a sufficiently quick reduction in the number of simulated (pairs of) trajectories $P_{\ell}$ as $\ell$ increases, it is possible to show that the multilevel Monte Carlo estimator is able to achieve the same mean square error as the classical Monte Carlo estimator at a lower computational cost.

For a detailed overview of the multilevel Monte Carlo method and its properties, we refer the reader to Giles2015 . Here, we limit ourselves to mentioning the theorem originally presented by Giles in Giles2008 and further generalized in Cliffe2011 and Giles2015 , which gives an upper bound for the computational complexity of multilevel Monte Carlo.

Remark 1 (Notation)

For conciseness, we use $\hat{F}_{\ell}$ instead of $F\left(X^{N_{\ell}}_{p,\Delta t_{\ell}},V^{N_{\ell}}_{p,\Delta t_{\ell}}\right)$ further in this work when considering approximations of the function $F\left(X(t^{*}),V(t^{*})\right)$ at level $\ell$. For $F\left(X(t^{*}),V(t^{*})\right)$ we use the short-hand notation $F$. Additionally, we introduce the requested bound on the mean square error $E$, and a series of positive constants $\alpha$, $\beta$, $\gamma$, $c_{1}$, $c_{2}$, $c_{3}$ and $c_{4}$, which are problem-dependent. In the theorem, $e$ is used to denote Euler’s constant.

Theorem 3.1

Let $F$ denote a random variable, and let $\hat{F}_{\ell}$ denote the corresponding level $\ell$ numerical approximation. If there exist independent estimators $\hat{Y}_{\ell}$ based on $P_{\ell}$ Monte Carlo samples, each with expected cost $C_{\ell}$ and variance $V_{\ell}$, and positive constants $\alpha$, $\beta$, $\gamma$, $c_{1}$, $c_{2}$, $c_{3}$ such that $\alpha\geq\frac{1}{2}\min(\beta,\gamma)$ and

$\operatorname{\mathbb{E}}\left[\hat{Y}_{\ell}\right]=\begin{cases}\operatorname{\mathbb{E}}\left[\hat{F}_{0}\right]&\ell=0,\\ \operatorname{\mathbb{E}}\left[\hat{F}_{\ell}-\hat{F}_{\ell-1}\right]&\ell>0,\end{cases}$* (Estimator notation)* 2. 2.

$\left|\operatorname{\mathbb{E}}\left[\hat{F}_{\ell}-F\right]\right|\leq c_{1}2^{-\alpha\ell}$, (Bias decreases with increasing $\ell$) 3. 3.

$V_{\ell}\leq c_{2}2^{-\beta\ell}$, (Variance decreases with increasing $\ell$) 4. 4.

$C_{\ell}\leq c_{3}2^{\gamma\ell}$, (Cost increases with increasing $\ell$)

then there exists a positive constant $c_{4}$ such that for any $E<e^{-1}$ there are values $L$ and $N_{\ell}$ for which the multilevel estimator (26) has a mean square error with bound $\text{MSE}\equiv\operatorname{\mathbb{E}}(\hat{Y}-\operatorname{\mathbb{E}}[F])^{2}<E^{2}$ with a computational complexity $C$ with bound

[TABLE]

The essential idea behind the theorem is the following. If the variance decreases faster than the cost increases, $\beta>\gamma$, then most of the work will be done at coarser levels. In this case the same complexity is achieved as a single level Monte Carlo simulation at the coarsest level. This is case in which the method is the most effective. If the variance decreases slower than the cost increases, $\beta<\gamma$, then most of the work is done at the finer levels. In this case the multilevel Monte Carlo has a much better asymptotic complexity as the finest level requires $\mathcal{O}(1)$ samples, each with cost $\mathcal{O}\left(E^{-\gamma/\alpha}\right)$. However, the constant $c_{4}$ will be large, so the multilevel speedup will be relatively small. In the intermediate case, $\beta=\gamma$, the computational cost is spread over the levels. In this case the $\log^{2}E$ factor corresponds with the number of required levels. We will refer back to this theorem in Section 5.4, where we demonstrate the convergence of our scheme.

4 Correlating particle trajectories

Recall that the time steps at levels $\ell$ and $\ell-1$ are related through $\Delta t_{\ell-1}=M\Delta t_{\ell}$. At the fine level $\ell$, we therefore define a sub-step index $m\in\{1,\dots,M\}$, i.e., we write $X^{n,m}_{p,\Delta t_{\ell}}\approx X_{p}(n\Delta t_{\ell-1}+m\Delta t_{\ell})\equiv X_{p}((nM+m)\Delta t_{\ell})$. Subsequently, we introduce a coupled pair of simulations spanning a time step with size $\Delta t_{\ell-1}$: (i) a simulation at level $\ell-1$, using a single time step of size $\Delta t_{\ell-1}$ and (ii) a simulation at level $\ell$, using $M$ time steps of size $\Delta t_{\ell}$:

[TABLE]

with $\xi^{n}_{p,\ell-1},\xi^{n,m}_{p,\ell}\sim\mathcal{N}(0,1)$, $V_{p,\Delta t_{\ell-1}}^{n}\sim\mathcal{M}_{\Delta t_{\ell-1}}(v)$ and $V_{p,\Delta t_{\ell}}^{n,m}\sim\mathcal{M}_{\Delta t_{\ell}}(v)$.

The differences in (21) will only have low variance if the simulated paths up to $X^{N_{\ell}}_{\Delta t_{\ell},p}=X^{N_{\ell-1},0}_{\Delta t_{\ell},p}$ and $X^{N_{\ell-1}}_{\Delta t_{\ell-1},p}$ are correlated. To achieve correlation, we need to perform coupled simulations of the scheme (14)–(16) with different time steps. In each time step using the asymptotic-preserving particle scheme, there are two sources of stochastic behavior. On the one hand, a new Brownian increment $\xi_{p}^{n}$ is generated for each particle in each transport-diffusion step (14). On the other hand, in each collision step (16), a fraction of particles randomly get a new velocity $V_{p}^{n,*}$ for use in the following time step. Particle trajectories can be coupled by correlating the random numbers used for the individual particles in the transport-diffusion and collision phase of each time step.

To this end, we will not draw independent samples at level $\ell-1$, but instead compute the values $\xi^{n}_{p,\ell-1}$ and $V_{p,\Delta t_{\ell-1}}^{n+1}=V_{p,\Delta t_{\ell-1}}^{n,*}$ at level $\ell-1$ based on the respective values of $\xi^{n,m}_{p,\ell}$ and $V_{p,\Delta t_{\ell}}^{n,m,*}$ in the fine simulation, while ensuring their correct statistical distribution. To achieve this, we thus first perform the $M$ fine simulation steps. Using the random numbers in these $M$ sub-steps, we then compute values with which to define the value of the random numbers in the single step of size $\Delta t_{\ell-1}$. The way in which we compute these dependent coarse random numbers at level $\ell-1$ is described in Algorithm 1. If the numbers at level $\ell-1$ have the correct statistical distribution, the coarse simulation statistics are not affected by the introduced correlation. Generalization of Algorithm 1 to higher dimensional domains for the position and velocity, can be done by simply replacing the random values $\xi$ and $\bar{V}$ with random vector quantities of equal dimension to the respective domains.

In the remainder of this section we will motivate the different steps in Algorithm 1. To achieve correlated simulations, we need to couple two sources of random behavior in (27). On the one hand, we generate a new normally distributed $\xi_{p,\ell}^{n}$ in each transport-diffusion step. On the other hand, there is a possibility that a collision occurs at each simulation step, causing the selection of a new velocity $V_{p,\ell}^{n}$. To discuss how correlation is introduced in both of these cases, we decompose each step in (27) into two parts: the transport part and the diffusion part. Defining transport increments

[TABLE]

and Brownian increments

[TABLE]

we can rewrite (27) as

[TABLE]

We first motivate the correlation of the Brownian increments in line 5 of Algorithm 1 (Section 4.1). Next, we motivate the correlation of the transport increments in lines 6–11 of Algorithm 1 (Section 4.2). After giving the motivation for the algorithm, we briefly show that this algorithm indeed achieves correlated simulations in Section 4.3.

4.1 Brownian increments

We now consider just the Brownian increments (29) for two processes spanning a time interval $\Delta t_{\ell-1}$, one spanning the interval in a single time step and the other taking $M$ time steps of size $\Delta t_{\ell}$. After $M$ fine Brownian increments $\Delta W^{n,m}_{p,\ell}$ at level $\ell$, according to (29), we get a Brownian increment spanning $\Delta t_{\ell-1}$,

[TABLE]

Given that

[TABLE]

we can compute a standard normally distributed $\xi^{n}_{p,\ell-1}$ from the $\xi^{n,m}_{p,\ell}$ as

[TABLE]

giving us the expression on line 5 of Algorithm 1.

Correlating the simulations in this way means that both simulations follow the same Brownian path, and differences in the diffusion behavior only result from differences in the diffusion coefficients $D_{\ell}$ and $D_{\ell-1}$ due to the different time steps. In Figure 2 we show two particle trajectories, containing a series of increments $\Delta W^{n,m}_{p,\ell}$ and $\Delta W^{n}_{p,\ell-1}$, using coupled normally distributed numbers as described in (33), with $\varepsilon=0.5$, $\Delta t_{\ell}=0.2$ and $M=5$. We observe that the paths have similar behavior, i.e., if the fine simulation tends towards negative values, so does the coarse simulation and vice versa. Still, there is an observable difference between them. This is due to difference in simulation variance caused by the differing diffusion coefficients.

4.2 Transport increments

While correlating Brownian increments (29) is straightforward, correlating the transport increments (28) is more involved. Since we simulate level $\ell$ first, we have the increments $\Delta T^{n,m}_{p,\ell}$, spanning a total time interval $\Delta t_{\ell}$, at our disposal. Our goal is to use the random numbers in these increments to calculate a single increment $\Delta T^{n}_{p,\ell-1}$ that spans the same time interval. Note that, in the collision phase of the asymptotic-preserving particle scheme (Section 2.2), both the value of the velocity and the probability of collision depend on the value of the time step $\Delta t$, and therefore depend on the level $\ell$. The coupling is done in two steps. First, the occurrence of a collision in each simulation step of the coarse simulation is coupled to the occurrence of a collision in at least one of the $M$ sub-steps of the coupled fine simulation. Then, if a collision occurs at both level $\ell$ and $\ell-1$, we will correlate the new velocities generated at both levels.

4.2.1 Deciding upon collision in the coarse simulation

Let us first consider the simulation at level $\ell$. In each of the $M$ fine simulation time steps of size $\Delta t_{\ell}$, we need to decide whether or not the particle collided. To this end, we draw a random number $u^{n,m}_{p,\ell}\sim\mathcal{U}([0,1])$ in the evaluation of (16). If this number is larger than the probability that no collision has occurred in the the time step $p_{nc,\Delta t_{\ell}}=1-p_{c,\Delta t_{\ell}}$

[TABLE]

then a collision is performed, i.e., a new velocity is randomly drawn from $\mathcal{M}_{\Delta t_{\ell}}$ at the end of that time step. At least one collision takes place in the interval spanning $\Delta t_{\ell-1}$ if at least one of the generated $u_{p,\ell}^{n,m}$, $m\in\{0,\dots,M-1\}$, satisfies (34), i.e.,

[TABLE]

When (35) is satisfied, we want the probability of a collision taking place in the correlated coarse simulation to be large.

We thus wish to generate a $u_{p,\ell-1}^{n}\sim\mathcal{U}\left([0,1]\right)$ based on the value of $u^{n,\text{max}}_{p,\ell}$, which we can use to test for a collision at level $\ell-1$. Since the cumulative density function of $u^{n,\text{max}}_{p,\ell}$ is given by

[TABLE]

we get, by the inverse transform method, that ${\left(u^{n,\text{max}}_{p,\ell}\right)}^{M}\sim\mathcal{U}([0,1])$. This means that we can achieve our goal by setting

[TABLE]

giving line 6 of Algorithm 1.

A collision at level $\ell-1$ then occurs when

[TABLE]

giving line 7 of Algorithm 1.

4.2.2 Choosing a new velocity

We first show that a collision at level $\ell-1$ can only occur when at least one collision has taken place in the simulation at level $\ell$.

Lemma 1

Given a simulation interval $\Delta t_{\ell-1}$ and a pair of paths $X_{p,\ell-1}$ and $X_{p,\ell}$ spanning this interval with time step sizes, respectively, $\Delta t_{\ell-1}$ and $\Delta t_{\ell}$, with collision behavior correlated as in (37). Then, the absence of collision in the $M$ sub-steps of the fine simulation $X_{p,\ell}$ guarantees the absence of collisions in the coarse simulation $X_{p,\ell-1}$.

Proof

By the definition of $p_{nc,\Delta t_{\ell-1}}$ and the fact that $\Delta t_{\ell-1}=M\Delta t_{\ell}$, checking this statement corresponds to verifying that

[TABLE]

Using (38), this is equivalent to showing that

[TABLE]

We do this by rewriting (40) as

[TABLE]

and using the binomial theorem to obtain

[TABLE]

Since $\varepsilon$, $M$ and $\Delta t_{\ell}$ are all positive, the statement is proved. ∎

Remark 2

Note that the inequality in (42) can be made strict, so it is possible that a collision is performed in the simulation at level $\ell$ without a collision occurring at level $\ell-1$. This is expected behavior, given the $\Delta t$ dependence of the collision rate in (13).

Lemma 1 implies that we never select a new velocity for the simulation at level $\ell-1$ without already selecting a new velocity for in least one of the sub-steps at level $\ell$. If a collision is performed in the simulation at level $\ell-1$, we want to correlate it with the fine simulation velocity at the end of the time interval. We consider the transformation of $\mathcal{M}(v)$ given in (10). If we use the same $\bar{V}$ to sample from $\mathcal{M}_{\Delta t_{\ell}}(v)$ and $\mathcal{M}_{\Delta t_{\ell-1}}(v)$, we can expect the resulting velocities to be correlated. In each of the $M$ fine sub-steps containing a collision we draw a $\bar{V}^{n,m,*}_{p,\ell}$ for use in the subsequent time step. Given that the $\bar{V}^{n,m,*}_{p,\ell}$ are i.i.d., we can select one freely to use as $\bar{V}^{n,*}_{p,{\ell-1}}$. To maximise the correlation of the velocities at the end of the time interval, we choose to take the last generated $\bar{V}^{n,m,*}_{p,\ell}$, i.e., we choose

[TABLE]

giving line 8 of Algorithm 1.

4.2.3 Numerical illustration

In Figure 3, we show two particle trajectories, containing a series of increments $\Delta T^{n,m}_{p,\ell}$ and $\Delta T^{n}_{p,\ell-1}$ for the two speed model (13), using coupled uniformly distributed numbers as described in (37) and (43), with $\varepsilon=0.5$, $\Delta t_{\ell}=0.2$ and $M=5$. In this model, we sample from $\mathcal{M}_{\Delta t_{\ell}}(v)$ by drawing $\bar{V}^{n,m}_{p,\ell}$ from $\{-1,1\}$ and multiplying it with $\tilde{v}_{\Delta t_{\ell}}$.

We observe that the paths have the same signs in the velocities at the end of a coarse time step whenever a collision has taken place in both simulations. Still, there is an observable difference between them. This is due to the paths having different characteristic velocities, which are a function of the time step size. There is also a probability of no collision taking place in the coarse simulation, while a collision takes place in the fine simulation, as per Remark 2. For instance, no collision occurs at $t=8$ in the coarse simulation, while a collision takes place at time $t=7.4$ and $t=8$ in the fine simulation. However, by coincidence, the new velocity generated at $t=8$ in the fine simulation has the same sign as the velocity in the coarse simulation. This mismatch is also part of the bias we want to estimate. If we add Figure 2 to Figure 3, we get the paths shown in Figure 4.

4.3 Variance of correlated paths

To conclude this section, we present a brief numerical result in which we compute the variance of 10 000 pairs trajectories as shown in Figures 2–4. The results are shown in Figure 5. We leave a detailed study of the variance structure resulting from our algorithm to the subsequent sections and limit ourselves to two observations here. We see that the variance of differences (green line with triangles) falls consistently below that of the fine simulation (blue line with squares), meaning that fewer samples are needed to estimate the difference than are needed to directly estimate the result of the fine simulation.

We also see that the variance of differences of transport processes lies above that of the coarse transport process. This is due to performing simulations with time step sizes with order of magnitude as $\varepsilon^{2}$. In Section 6, we will show that this is the parameter regime in which the asymptotic-preserving model switches between simulating an approximation of the original kinetic equation (6) and the limiting diffusion equation. As the simulation with time step $\Delta t_{\ell-1}$ contains much less transport behavior than the simulation with time step $\Delta t_{\ell}$, the variance of their difference lies between their individual variances, even though the correlation is good.

5 Analysis of the multilevel scheme

As is clear from Theorem 3.1, the behavior of the mean and variance of the estimators (21) as a function of the level $\ell$ is of key importance in the analysis of a multilevel Monte Carlo method. In this section, we present and prove lemmas concerning the mean and variance of the difference in particle position and velocity for the coupled simulations (27) as a function of the time step $\Delta t_{\ell}$ and the refinement factor $M$ at some fixed time $t^{*}$. When considering the particle position, we look separately at the Brownian increments (Section 5.1), the transport increments (Section 5.2) and particle velocities (Section 5.3). In Section 5.4, we make use of these lemmas to prove a bound on the computational complexity of our multilevel Monte Carlo scheme.

For convenience, we introduce an additional notation for the difference of an arbitrary pair of coupled increments at arbitrary levels $\ell$ and $\ell-1$ at time $t_{n}=n\Delta t_{\ell-1}$:

[TABLE]

5.1 Brownian increments

We first show that the Brownian increments have expected value 0.

Lemma 2

Given a simulation interval $\Delta t_{\ell-1}$, an independent fine simulation containing Brownian increments $\Delta W^{n,m}_{p,\ell}$ with time step sizes $\Delta t_{\ell}$ and a coarse simulation containing Brownian increments $\Delta W^{n}_{p,\ell-1}$ with time step size $\Delta t_{\ell-1}$, generated by Algorithm 1. Then,

[TABLE]

Proof

By the martingale property of the increments and the fact that the coupling preserves the statistics of the stochastic process at level $\ell-1$, we have

[TABLE]

The proof then directly follows from

[TABLE]

∎

To prove convergence of the scheme in Section 5.4, we require the following convergence result on the variance of the summed Brownian increments:

Lemma 3

Given a fixed ratio $M$ between two time step sizes $\Delta t_{\ell-1}=M\Delta t_{\ell}$, a fixed simulation horizon $t^{*}=NM\Delta t_{\ell}$, an independent fine simulation containing Brownian increments $\Delta W^{n,m}_{p,\ell}$ with time step sizes $\Delta t_{\ell}$ and a coarse simulation containing Brownian increments $\Delta W^{n}_{p,\ell-1}$, generated by Algorithm 1 with time step size $\Delta t_{\ell-1}$. Then, the variance

[TABLE]

converges $\mathcal{O}(\Delta t_{\ell})$ to 0, asymptotically for $\Delta t_{\ell}\to 0$.

Proof

To calculate the variance of the difference $\Delta_{W,n}$ for each Brownian increment at level $\ell$, we write

[TABLE]

in which we used independence of the random variables at level $\ell$.

We thus need the variances of individual increments at levels $\ell$ and $\ell-1$, which follow trivially from (29):

[TABLE]

The covariance between the increment $\Delta W^{n}_{p,\ell-1}$ at level $\ell-1$ and each sub-increment $\Delta W^{n,m}_{p,\ell}$ at level $\ell$ can be computed using (45) and (46):

[TABLE]

Using (51), (52) and (56), we can elaborate (50) to obtain

[TABLE]

Note that (58) gives the variance $\operatorname{\mathbb{V}}\left[\Delta_{W,n}\right]$ as a function of the time steps at both levels, $\Delta t_{\ell}$ and $\Delta t_{\ell-1}$, which also appear in the diffusion coefficients $D_{\Delta t_{\ell}}$ and $D_{\Delta t_{\ell-1}}$. Equivalently, using the relation $\Delta t_{\ell-1}=M\Delta t_{\ell}$, the variance can be written in terms of $\Delta t_{\ell}$ and the refinement factor $M$.

As the time steps are independent, the total variance after $N$ steps is

[TABLE]

where we omit the subscript in the right hand side, as the increments are i.i.d. We now compute the Maclaurin series in $\Delta t_{\ell}$ of (59) as

[TABLE]

We refer to the supplementary materials for a detailed computation of (60). ∎

Remark 3

In the supplementary materials we further show that

[TABLE]

This limit matches the expected convergence of (13) to (5), as Brownian motion, i.e., the Monte Carlo discretization of (5), is unbiased in the time step.

5.2 Transport increments

We now present and prove similar lemmas concerning the transport increments, i.e., the at least one collision has already taken place in both simulations, so that the initial condition on the velocity can be neglected.

Lemma 4

Given a simulation interval $\Delta t_{\ell-1}$, an independent fine simulation containing transport increments $\Delta T^{n,m}_{p,\ell}$ with time step sizes $\Delta t_{\ell}$ and a coarse simulation containing transport increments $\Delta T^{n}_{p,\ell-1}$ with time step size $\Delta t_{\ell-1}$, generated by Algorithm 1. Then,

[TABLE]

Proof

Based on (11) and the fact that the correlation preserves the coarse model statistics (as demonstrated in Section 4.2), we can show that the expected value of both paths is zero,

[TABLE]

The proof then follows from

[TABLE]

∎

Lemma 5

Given a fixed ratio $M$ between two time step sizes $\Delta t_{\ell-1}=M\Delta t_{\ell}$, a fixed simulation horizon $t^{*}=NM\Delta t_{\ell}$, an independent fine simulation containing transport increments $\Delta T^{n,m}_{p,\ell}$ with time step sizes $\Delta t_{\ell}$ and a coarse simulation containing transport increments $\Delta T^{n}_{p,\ell-1}$, generated by Algorithm 1 with time step size $\Delta t_{\ell-1}$. Then, the variance

[TABLE]

converges $\mathcal{O}(\Delta t_{\ell})$ to 0, asymptotically for $\Delta t_{\ell}\to 0$.

Proof

As there is a non-zero probability for both the simulation at level $\ell$ and at level $\ell-1$ that no collision occurs in a given time step $n$, the differences between these two simulations are themselves correlated across time steps. The variance of the difference after $N$ steps is thus given by

[TABLE]

[TABLE]

Given Lemma 4, we compute the variances of the individual increments as the expectation of their squared values:

[TABLE]

The remaining covariance sums are workout out in Appendix A.

We now plug the computed variances and covariances ((67), (99), (105) and (128)) into (65) and compute the first term of the Maclaurin series in $\Delta t_{\ell}$:

[TABLE]

We refer to the supplementary materials for a detailed computation of (68). ∎

Remark 4

In the supplementary materials we further show that

[TABLE]

This limit matches the expected convergence of (13) to (5), as (5) no longer contains transport behavior.

5.3 Velocity expectations and variances

We now present the same lemmas concerning the expectation and variance of the differences of velocities.

Lemma 6

Given a simulation interval $\Delta t_{\ell-1}$, an independent fine simulation containing velocities $V^{n,m}_{p,\ell}$ with time step size $\Delta t_{\ell}$ and a coarse simulation containing velocities $V^{n}_{p,\ell-1}$ with time step size $\Delta t_{\ell-1}$, generated by Algorithm 1. Then, given that both simulations have experienced at least one collision, $\operatorname{\mathbb{E}}\left[V_{p,\Delta t_{\ell}}^{n,m}-V_{p,\Delta t_{\ell-1}}^{n}\right]=0$.

Proof

The expected values of the individual velocities of both fine and coarse simulations are zero by (11)

[TABLE]

This means that the expected value of their difference is also zero. ∎

Lemma 7

Given a fixed ratio $M$ between two time step sizes $\Delta t_{\ell-1}=M\Delta t_{\ell}$, a fixed simulation horizon $t^{*}=NM\Delta t_{\ell}$, an independent fine simulation with time step sizes $\Delta t_{\ell}$ and a coarse simulation, generated by Algorithm 1 with time step size $\Delta t_{\ell-1}$. Then, given that both simulations have experienced at least one collision, the variance $\operatorname{\mathbb{V}}\left[V_{p,\Delta t_{\ell}}^{n,m}-V_{p,\Delta t_{\ell-1}}^{n}\right]$ converges $\mathcal{O}(\Delta t_{\ell})$ to 0, asymptotically for $\Delta t_{\ell}\to 0$.

Proof

We compute the variance of this difference as

[TABLE]

where we use (11) in the step from (72) to (73) and where the step from (73) to (74) follows the same logic as the calculation of (99) in Appendix A.

We now compute the Maclaurin series in $\Delta t_{\ell}$ of (75) as

[TABLE]

We refer to the supplementary materials for a detailed computation of (76). ∎

Remark 5

In the supplementary materials we further show that

[TABLE]

This limit matches the expected convergence of (13) to (5), as (5) no longer contains transport behavior.

5.4 Proof of convergence

Now that we have established some key properties of the correlation between the coupled trajectories at level $\ell$ and $\ell-1$, we have everything in place to derive bounds on the difference estimators (21). To this end, we assume that the quantity of interest $F(x,v)$ is Lipschitz continuous in both position and velocity, i.e., there exist constants $K_{x}$ and $K_{v}$ so that,

[TABLE]

holds, for all values in the domains $\mathcal{D}_{x}$ and $\mathcal{D}_{v}$. We can now prove the convergence of our scheme, i.e., consistency as $\ell\to\infty$. This proof is based on Lemmas 2–7 and is structured as follows: First we present three lemmas, which verify Assumptions 2–4 in Theorem 3.1. We then present a convergence theorem for our scheme in Theorem 5.1.

First, we verify the rate of decreasing bias (Theorem 3.1, Assumption 2).

Lemma 8

Given $F(x,v)$, Lipschitz in position and velocity, and a sequence of approximations $\ell=0,\dots,\infty$, coupled as described in algorithm 1 with time steps $\Delta t_{\ell-1}=M\Delta t_{\ell}$, then $\exists L_{1},0<L_{1}<\infty:\exists c_{1}:\forall\ell\geq L_{1}:\left|\operatorname{\mathbb{E}}\left[\hat{F}_{\ell}-F\right]\right|\leq c_{1}2^{-\alpha\ell}$, with $\alpha=\log_{2}(M)$.

Proof

By (78) we have

[TABLE]

with $\left(X^{*},V^{*}\right)$ the point in the velocity-phase space which produces the expected value $\hat{F}$ of $F(x,v)$, which exists by the mean value theorem. Given that both time-splitting and the IMEX-equation are linear approximations in $\Delta t$ in the limit $\Delta t\to 0$, we can observe that both $\operatorname{\mathbb{E}}\left[\left|X^{N}_{p,\Delta t}-X^{*}\right|\right]$ and $\operatorname{\mathbb{E}}\left[\left|V^{N}_{p,\Delta t}-V^{*}\right|\right]$ go to zero with the weak order of the explicit first order simulation method, i.e.,

[TABLE]

meaning there exists an upper bound $c_{1}2^{-\log_{2}(M)\ell}$, once $\ell$ is sufficiently large. ∎

Second, we verify the rate of decreasing variance (Theorem 3.1, Assumption 3).

Lemma 9

Given $F(x,v)$, Lipschitz in position and velocity, and a sequence of approximations $\ell=0,\dots,\infty$, coupled as described in algorithm 1 with time steps $\Delta t_{\ell-1}=M\Delta t_{\ell}$, then $\exists L_{2},0<L_{2}<\infty:\exists c_{2}:\forall\ell\geq L_{2}:\operatorname{\mathbb{V}}\left[\hat{F}_{\ell}-F_{\ell-1}\right]\leq c_{2}2^{-\beta\ell}$, with $\beta=\log_{2}(M)$.

Proof

By using the Lipschitz property and Lemmas 2, 4 and 6 we compute the following bound on the variance of the difference estimators:

[TABLE]

As the linear term is the first nonzero term in (60), (68) and (76) we write

[TABLE]

meaning there exists an upper bound $c_{2}2^{-\log_{2}(M)\ell}$, once $\ell$ is sufficiently large. ∎

Third, we verify the rate of increasing cost (Theorem 3.1, Assumption 4).

Lemma 10

For a sequence of difference estimators $\ell=0,\dots,\infty$, with the time step sizes following $\Delta t_{\ell-1}=M\Delta t_{\ell}$, correlated as by algorithm 1, the cost per sample $C_{\ell}$ decreases as $C_{\ell}\leq c_{3}2^{\gamma\ell}$, with $\gamma=\log_{2}(M)$ and $c_{3}$ constant, for all $\ell>0$.

Proof

The number of simulation steps $N_{\ell}$ needed for a simulation at level $\ell$ is $N_{0}M^{\ell}$, meaning that a difference estimator at level $\ell$ costs $(M+1)M^{\ell-1}$ times that of a single simulation at level 0. This means that

[TABLE]

for all $l$. ∎

Finally, we combine Theorem 3.1 and Lemmas 8–10 to prove the convergence rate of our scheme.

Theorem 5.1

Given $F(x,v)$, Lipschitz in position and velocity, and a sequence of approximations $\ell=0,\dots,\infty$, coupled as described in algorithm 1 with time steps $\Delta t_{\ell-1}=M\Delta t_{\ell}$. The multilevel Monte Carlo method, applied to this sequence of approximations algorithm, converges with a mean square error $\operatorname{\mathbb{E}}(\hat{Y}-\operatorname{\mathbb{E}}\left[F\right])^{2}<E^{2}$ and a computational complexity $C$ bounded by $\operatorname{\mathbb{E}}\left[C\right]\leq c_{4}E^{-2}(\log E)^{2}$, for a given constant $c_{4}$ and a sufficiently small $E>0$.

Proof

As $E$ constrains the bias which in turn constrains the finest level time step size $\Delta t_{L}$, a sufficiently small $E$ will enforce that the number of levels $L\geq\max(L_{1},L_{2})$. The proof then follows by insertion of Lemmas 8, 9 and 10 into Theorem 3.1. ∎

6 Selecting a level strategy

In this section, we combine analytical and numerical results to determine a strategy for selecting levels when using our asymptotic-preserving multilevel scheme. The goal is to minimize the total simulation cost. In Section 6.1 we first study the the bias and variance structure of the multilevel scheme in function of the simulation time step size $\Delta t_{\ell}$. Based on these insights, we propose two level strategies in Section 6.2, which are then compared in terms of computational cost.

6.1 Bias and variance structure

Theorem 5.1 proves that the scheme converges asymptotically as $\Delta t_{\ell}$ decreases, for general values of $\varepsilon$, $M$ and $t^{*}$, but does not give us the full picture for a larger range of values in $\Delta t_{\ell}$. To get a more complete picture, we consider the bias and variance of the difference estimators for a simple quantity of interest for different values of $\Delta t_{\ell}$. To this end we fix $t^{*}=5$ and $M=2$. At level $\ell=0$, we set $\Delta t_{0}=2.5$. At finer levels ($\ell\geq 1$) we set $\Delta t_{\ell}=\Delta t_{\ell-1}/M=\Delta t_{0}/M^{\ell}$. We fix the number of samples per difference estimator at 100 000. For a selection of values of $\varepsilon$, we calculate the expected value and variance of the individual samples $\hat{F}_{\ell}$ and difference estimators $\hat{F}_{\ell}-\hat{F}_{\ell-1}$ of the squared particle position, as a function of $\Delta t_{\ell}$, for $1\leq\ell$. We choose $\varepsilon=10$ (Figure 6), $\varepsilon=1$ (Figure 7), $\varepsilon=0.1$ (Figure 8) and $\varepsilon=0.01$ (Figure 9). We also plot the analytical bound on the variance given by (59) and (65), where we use the respective $K_{x}$-values 1.5, 5, 8 and 8, based on visually comparing results, and take $K_{v}=0$, as the QoI is independent of the velocity.

The regime $\Delta t\ll\varepsilon^{2}$. In Figures 6 through 8, we see that, as the level $\ell$ increases, the slopes of both the mean and variance curves for the differences approach an asymptotic limit $\mathcal{O}\left(\Delta t\right)$ for $\Delta t\ll\varepsilon^{2}$. This observation matches the weak convergence order of the Euler-Maruyama scheme, used to simulate the model (14)–(16), as well as the expected behavior from the time step dependent bias in the asymptotic-preserving model. This confirms the expected behavior from (60) and (68) as well as Lemmas 8 and 9. In this regime, the existing theory for multilevel Monte Carlo methods can be applied, e.g., on the required number of samples per level and conditions for adding levels Giles2015 .

The regime $\Delta t\gg\varepsilon^{2}$. For time steps $\Delta t\gg\varepsilon^{2}$, however, we see in Figures 8 and 9 that both the mean and the variance curves increase geometrically in terms of increasing level. To explain this perhaps counterintuitive result, we will look at the limit of the modified Goldstein-Taylor model when $\Delta t$ tends to infinity. In this limit, the model (13) converges to the heat equation:

[TABLE]

This means that taking increasingly larger time steps in (13) is equivalent to taking the limit $\varepsilon\to 0$.This observation is precisely the asymptotic-preserving property of the particle scheme of Section 2.2.

That the scheme approaches two different limiting models in these two limits can be seen most clearly in Figures 6 and 8. In Figure 8, the curves for the mean and variance of the differences $\hat{F}_{\ell}-\hat{F}_{\ell-1}$ (orange lines with circles) decrease for both small and large $\Delta t$, as the model converges to the two limits. In the right hand panel of Figure 6 we see that the variance of the individual simulations at level $\ell$ (blue line with squares) changes drastically as a function of $\Delta t_{\ell}$ in the region where it is of the same order of magnitude as $\varepsilon^{2}$. This is caused by the approximated models for large and small $\Delta t$ having differences in behavior, which are significant enough to be observed when plotted. The scheme thus converges to different equations for the two limits in $\Delta t$. For small $\Delta t$, there is convergence to (12). For large $\Delta t$, there is convergence to (5). In practice, the size of $\Delta t$ is limited by the simulation time horizon, so it is not possible to get arbitrarily close to (5) by increasing the time step size.

Connecting the two regimes. Combining the observations from the two limits ($\varepsilon$ tending to zero and $\Delta t_{\ell}$ tending to zero) in the time step size gives an intuitive interpretation to the multilevel Monte Carlo method in this setting: the method can be interpreted as correcting the result of a pure diffusion simulation by decreasing $\Delta t$ to get a good approximation of the transport-diffusion equation that describes the behavior for a given value of $\varepsilon$. The peak of the variance of the differences lies near $\Delta t\approx\varepsilon^{2}$. This makes sense, as this is the region where the model parameters $D_{\Delta t}$ and $\tilde{v}_{\Delta t}$ vary the most as a function of $\Delta t$. We also see a dip in the mean of the difference curves in the region of $\Delta t\approx\varepsilon^{2}$. A full analysis of the behavior that occurs in the transition between the asymptotic regimes is left for future work.

6.2 Performance and level placement strategy

In Section 6.1, we experimentally verified the asymptotic convergence rates of the bias and variance in function of increasing level number that were used in the proof of Theorem 5.1. In doing so, we observed an increasing mean and variance for the difference estimators in the region $\Delta t\gg\varepsilon^{2}$, for $M=2$ and a number of values for $\varepsilon$. In this regime, it makes little sense to include a full sequence of levels, as Theorem 3.1 only claims a speedup over classical Monte Carlo in the case of decreasing variance. Levels in this region can be interpreted as producing bias estimators that are orders of magnitude smaller than the bias in the model which they are estimating, which means wasted computation.

That a full sequence of levels makes no sense for time step sizes larger than $\varepsilon^{2}$, is therefore intuitively clear. However the question still remains as to what the best approach is to selecting levels. We consider two possible simulation strategies:

• •

Strategy 1: Geometric sequence, starting from $\varepsilon^{2}$:

We generate an initial estimate of the quantity of interest at level zero, where we simulate to $t^{*}$ using $\Delta t_{0}=\varepsilon^{2}$. 2. 2.

We continue to generate a geometric sequence of levels until an acceptably low bias has been achieved, i.e., $\Delta t_{l}=\varepsilon^{2}M^{-l}$ for $l>0$.

• •

Strategy 2: Additional inclusion of a single coarse level:

We generate an initial estimate of the quantity of interest at level zero, where we simulate to $t^{*}$ using $\Delta t_{0}=t^{*}$. 2. 2.

At level 1 we perform correlated simulations to $t^{*}$ using $\Delta t_{0}=t^{*}$ and $\Delta t_{1}=\varepsilon^{2}$. 3. 3.

We continue to generate a geometric sequence of levels until an acceptably low bias has been achieved, i.e., $\Delta t_{l}=\varepsilon^{2}M^{1-l}$ for $l>1$.

In Appendix B, we apply multilevel theory to compare both strategies for some very simple quantities of interest. In the considered cases, Strategy 2 proves more efficient. The theoretical approach does not expand to general quantities of interest, however. In the following two subsections (Section 6.2.1 for Strategy 1 and Section 6.2.2 for Strategy 2), we compare these strategies for the quantity of interest considered in Section 6.1. We choose to set $M=2$ and $\varepsilon=0.1$, and reduce the time horizon to $t^{*}=0.5$. This gives us an expensive, but computationally feasible problem. The number of samples per level is derived using the formula Giles2015

[TABLE]

with $E$ the desired bound on the root mean square error, $C_{\ell}$ the computational cost of the estimator at level $\ell$, and $V_{\ell}$ the estimated variance of the estimator at level $\ell$, i.e., $V_{\ell}=\mathbb{V}\left[\hat{F}_{\ell}-F_{\ell-1}\right]$, where we set $F_{-1}\equiv 0$. An initial estimate for $V_{\ell}$ is computed using 40, 500 and 1 000 initial samples for respective $E$-values 0.1, 0.01 and 0.001. The criteria for adding levels and determining convergence are as described in Giles2015 . The cost per sample is determined relative to the cost of a trajectory simulated with $\Delta t=\varepsilon^{2}$.

The code for performing the numerical experiments can be found at github.com/ELoevbak/APMLMC, together with the data files containing the simulation results.

6.2.1 Simulating a geometric sequence

The results of the simulations with strategy 1 for $E$ values 0.1, 0.01 and 0.001 are given in Tables 1 through 3. In these tables, we list the time step size $\Delta t_{\ell}$, number of samples $P_{\ell}$, variance of the fine simulations $\mathbb{V}\left[\hat{F}_{\ell}\right]$, expected value $\mathbb{E}\left[\hat{F}_{\ell}-F_{\ell-1}\right]$ and variance $V_{\ell}$ of the differences of simulations, estimated variance of the estimator $\mathbb{V}\left[\hat{Y}_{\ell}\right]$, cost per sample $C_{\ell}$ and level cost $P_{\ell}C_{\ell}$. The level $\ell$ estimator variance is estimated as

[TABLE]

We see that the number of samples $P_{\ell}$ needed to keep $\sum_{\ell=0}^{L}\mathbb{V}\left[\hat{Y}_{\ell}\right]<E^{2}$ decreases drastically in function of $\ell$. We also see that $\mathbb{E}\left[F_{L}-F_{L-1}\right]<E$. The cost per level $P_{\ell}C_{\ell}$ is spread quite evenly over the levels, which is to be expected as the geometric factor with which the cost increases with $\ell$ is asymptotically the same as that with which $V_{\ell}$ decreases. In short, we thus achieve the bias of the finest level, while a large amount of variance reduction is performed in the coarser levels. We can thus conclude that the experimental results match the expected behavior of the multilevel Monte Carlo method.

The total cost of each multilevel simulation, relative to the cost of a single sample at the coarsest level is computed as the sum of the costs at each level, i.e., the sum of the right most column of Tables 1 through 3. We can estimate the cost for an equivalent classical Monte Carlo simulation by considering that one needs to perform

[TABLE]

samples with the fine time step at level $L$, to achieve the same bias and variance as the multilevel estimator. The cost of each sample in the classic Monte Carlo estimator is $\frac{2}{3}C_{L}$, as we do not need to perform a correlated coarse simulation. Note that, for $E=0.1$, the variance $\mathbb{V}\left[F_{L}\right]$ is estimated using very few samples. One should thus be careful about drawing further conclusions from these tables than those made here.

We now compare the cost of the classical and multilevel Monte Carlo simulations in Table 4.

As can be concluded from the table, the multilevel Monte Carlo scheme gives a significant computational advantage when we want to compute low bias results in the setting of the modified Goldstein-Taylor model. This speedup increases as the requested accuracy of the simulation is increased and is expected to asymptotically scale with $\left(E\log^{2}E\right)^{-1}$ as the requested root mean square error is further decreased.

6.2.2 Simulating with a very coarse level

We now add an extra coarse level with $\Delta t_{0}=0.5$, and repeat the experiment as before. The results are shown in Tables 5 through 7. In these tables, we see that very little work is done on the coarsest level in comparison with the levels in the geometric sequence. The extra level thus does not have a significant cost, in comparison with the rest of the simulation. We observe that the expected behavior of the multilevel Monte Carlo method, as discussed in the previous section, is also present when including the coarser level.

We present a cost comparison with and without the coarse level in Table 8. We see that including a very coarse level consistently gives a speedup, matching expectations from multilevel theory. One observation is that the speedup becomes less significant as the requested root mean square error $E$ decreases. This makes sense as the higher the requested accuracy, the more levels are needed, and the smaller the influence of the coarse level strategy. Another thing to note is that, although $V_{1}$ is smaller than $V_{0}$, it is still relatively large, and much larger than $V_{\ell}$ in the following fine levels. We believe that it may be possible to further reduce the variance of level 1 in this strategy by using a different correlation strategy, in which we take into account that the coefficients in (13) vary strongly if $\Delta t$ takes values with different orders of magnitude.

7 Conclusion

We presented a multilevel Monte Carlo scheme for simulating a generic class of kinetic equations using asymptotic-preserving particle schemes. Although the scheme was derived in one dimension, it can be generalized to higher dimensional simulations with little extra effort. After presenting the scheme, we analyzed its convergence behavior for general functions of the particle position and velocity and provided some insights into level selection strategies.

Some analytical properties of the variance of the difference in position and velocity of two particle simulations were derived. Using these properties, we studied the behavior of the multilevel scheme as a function of the simulation time step $\Delta t_{l}$ for general Lipschitz quantities of interest. We proved that the multilevel Monte Carlo scheme converges with a computational cost that is asymptotically bounded in the root mean square error bound $E$ by $\mathcal{O}\left(E^{-2}\log^{2}(E)\right)$. The scheme’s speedup over single level Monte Carlo was confirmed through numerical experiments.

We also compared two approaches to select the levels in the multilevel scheme. After both theoretical analysis and numerical experiments, we concluded that the best option given the proposed correlation strategy is to start with a coarse level with $\Delta t_{0}=t^{*}$, followed by a geometric sequence of levels, decreasing from $\Delta t_{1}=\varepsilon^{2}$.

This work is a first step in combining the multilevel Monte Carlo method to asymptotic-preserving particle schemes for kinetic equations. In future work, we intend to look at alternative ways to increase the correlation of level 1 in the simulation strategy with a very coarse level, as this is of key importance to reduce the high computational cost of simulations for small values of $\varepsilon$. We also intend to expand our simulation code to cope with more complex models, including higher dimensional cases, absorption and position dependent model parameters, making the scheme directly applicable in relevant applications, such as fusion reactor design. An expansion to multi-index Monte Carlo, where the value of both $\Delta t$ and $\varepsilon$ are simultaneously varied, can also be considered.

Acknowledgments

We thank the anonymous reviewer for their substantial suggestions which have significantly improved this article. Emil Løvbak is funded by the Research Foundation - Flanders (FWO) under SB-Fellowship number 1SB1919N. The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and the Flemish Government – department EWI.

Conflict of interest

The authors declare that they have no conflict of interest.

Appendix A Computing the covariance sums needed to prove Lemma 5

A.1 $\displaystyle\sum_{m=0}^{M-1}\operatorname{\mathrm{Cov}}\left(\Delta T^{n,m}_{p,\ell},\Delta T^{n}_{p,\ell-1}\right)$

The covariance between the coarse increment $n$ and a fine sub-increment $(n,m)$, given Lemma 4, can be written out as

[TABLE]

To calculate the expected value in (92), we need to consider the probabilities of coupled collisions taking place in the correlated simulations. If a collision takes place in a given set of $M$ fine increments, the probability that the correlated coarse simulation time step will also simulate a collision is given by

[TABLE]

From Lemma 1, we know that it is not possible for a collision to take place in the coarse simulation, without a collision taking place in the fine simulation. This leaves three possibilities, when considering collision behavior in coarse time step $n-1$:

• •

Both at level $\ell-1$ and at level $\ell$, no collision occurred in time step $n-1$. In this case, time step $n-1$ will not affect the correlation of the velocities between the simulations. If the velocities were correlated at the beginning of time step $n-1$, they will still be so at the end of the time step, and vice versa. In this case, we thus need to look at step $n-2$, and so on, until we reach a past time step that satisfies one of the following two cases.

• •

A collision occurred at level $\ell$ in time step $n-1$, but not at level $\ell-1$. In this case, a new $\bar{V}^{n-1,m,*}_{p,\ell}$ was drawn from $\mathcal{M}(v)$, for some $m$, independently of the value of $\bar{V}^{n-1,*}_{p,\ell-1}$. Because all sampled velocities are independent, there is no correlation between $\bar{V}^{n,m-1}_{p,\ell}$ and $\bar{V}^{n-1}_{p,\ell-1}$, making the expected value of their product zero by (11).

• •

A collision occurred both at level $\ell$ and at level $\ell-1$ in time step $n-1$. We know that $\bar{V}^{n,0}_{p,\ell}=\bar{V}^{n-1,M,*}_{p,\ell}=\bar{V}^{n}_{p,\ell-1}=\bar{V}^{n-1,*}_{p,\ell-1}$ by (43). For every other fine simulation sub-step $m>1$ in time step $n$, there is a probability $1-p_{nc,\Delta t_{\ell}}^{m}$ that a collision has taken place in at least one of the steps $(n,0)$ through $(n,m)$. These collisions will not affect the coarse simulation until time step $n+1$, where $\bar{V}^{n,*}_{p,\ell}$ is used to generate a new velocity. So we have that $\bar{V}^{n,m}_{p,\ell}=\bar{V}^{n}_{p,\ell-1}$ with a probability $p_{nc,\Delta t_{\ell}}^{m}$, otherwise the two velocities are uncorrelated.

From the above list of possibilities, we conclude that the random variables $\bar{V}^{n,m}_{p,\ell}$ and $\bar{V}^{n}_{p,\ell-1}$ are equal if both of the following are true:

The last simulation step at level $\ell$ that underwent a collision is a sub-step of a simulation step $n^{\prime}$ at level $\ell-1$ which also underwent a collision. 2. 2.

The coarse simulation step $n^{\prime}$ is not the current coarse step $n$.

Otherwise, $\bar{V}^{n,m}_{p,\ell}$ and $\bar{V}^{n}_{p,\ell-1}$ are uncorrelated.

Making use of (11) and (95), we have that

[TABLE]

Plugging (96) into (92) gives

[TABLE]

where the final equality is a fortuitous coincidence that will shorten some expressions in Section 5.

A.2 $\displaystyle\sum_{m=0}^{M-2}\sum_{m^{\prime}=m+1}^{M-1}\operatorname{\mathrm{Cov}}\left(\Delta T^{n,m}_{p,\ell},\Delta T^{n,m^{\prime}}_{p,\ell}\right)$

As in Section A.1, we start by writing out the covariance between subsequent fine sub-increments $(n,m)$ and $(n,m^{\prime})$ as

[TABLE]

Calculating the sum of the covariance between subsequent fine increments is straightforward. To simplify notation, we introduce $\Delta m$ as shorthand for $m^{\prime}-m$. As the collision probability is constant across time steps, the covariance is given by

[TABLE]

Making use of (11), we get that the r.h.s. of (101) is equal to

[TABLE]

By making use of the identity

[TABLE]

and noting that $\frac{p_{nc,\Delta t_{\ell}}}{p_{c,\Delta t_{\ell}}}=\frac{\varepsilon^{2}}{\Delta t_{\ell}}$, we work out (102) as

[TABLE]

A.3 $\displaystyle\sum_{n=0}^{N-2}\sum_{n^{\prime}=n+1}^{N-1}\operatorname{\mathrm{Cov}}\left(\Delta_{T,n},\Delta_{T,n^{\prime}}\right)$

Making use of Lemma 4, we write out the covariance of differences of transport increments for a given $n$ and $n^{\prime}$ as

[TABLE]

We now calculate each term on the r.h.s. of (107) separately. We start by writing the expressions in terms of $\Delta n$, as a shorthand for $n^{\prime}-n$. In the remainder of the section we assume $n^{\prime}>n$, without loss of generality.

The increment covariance at level $\ell$ $(\textrm{I})$ is calculated from a similar starting point to (100):

[TABLE]

We now substitute the double summation in (109) with a two summations over $\mathrel{\Delta m}=|m^{\prime}-m|$:

[TABLE]

We now again make use of the identity (103), giving that (111) equals

[TABLE]

Calculating the covariance of increments at level $\ell-1$ (II) is also straightforward:

[TABLE]

The expression for the covariance between increments at level $\ell$ and level $\ell-1$ at differing time steps $n$ and $n^{\prime}$ depends on the relative position of the increments, i.e., whether the increment at level $\ell$ comes before that at level $\ell-1$, or not. If the fine increment at level $\ell$ comes first (III), we need to calculate

[TABLE]

To calculate the expected value of the r.h.s. of (115), we list the possible simulation behaviors at level $\ell-1$ in time step $n$, relative to the fine sub-step $m$:

• •

No collision occurs in the simulation at level $\ell$. This situation occurs with probability $p_{nc,\Delta t_{\ell}}^{M}.$ In this case, we have already established that no collision occurs at level $\ell-1$. The probability of the simulations being correlated at the end of time step $n$ is therefore equal to the probability of them being correlated at the start of time step $n$, given by (95).

• •

No collision occurs in sub-steps 0 through $m-1$, but at least one collision occurs in sub-steps $m$ through $M-1$. This situation occurs with probability

[TABLE]

If the simulation at level $\ell-1$ also has a collision, then $\bar{V}^{n+1}_{p,\ell-1}$ and $\bar{V}^{n,m}_{p,\ell}$ are independent by (43). If no collision occurred in the simulation at level $\ell-1$, which is the case with probability

[TABLE]

then $\bar{V}^{n}_{p,\ell-1}=\bar{V}^{n,m}_{p,\ell}$, if the trajectories were correlated at the beginning of time step $n$. This is the case with probability (95).

• •

At least one collision occurs in sub-steps 0 through $m-1$, but none occur in sub-steps $m$ through $M$. This occurs with probability

[TABLE]

By (43) we know that $\bar{V}^{n,*}_{p,\ell-1}=\bar{V}^{n,m}_{p,\ell}$, if a collision also occurs in time step $n$ of the simulation at level $\ell-1$. This happens with probability (95).

• •

Collisions happen both before sub-step $m$ and during or afterwards. In this case, there is no correlation between $\bar{V}^{n,*}_{p,\ell-1}$ and $\bar{V}^{n,m}_{p,\ell}$.

By adding the non-zero contributions from these four cases, we calculate the probability that $\bar{V}^{n,m}_{p,\ell}$ is correlated with $\bar{V}^{n,*}_{p,\ell-1}$. To get the probability of the correlation holding until $\bar{V}^{n^{\prime}}_{p,\ell-1}$, we multiply the sum by $p_{nc,\Delta t_{\ell-1}}^{\mathrel{\Delta n}-1}$. Given the properties (11), we can state that the right hand side sum in (115) is given by

[TABLE]

Making us of the fact that $\frac{p_{c,\Delta t_{\ell-1}}}{p_{nc,\Delta t_{\ell-1}}}=\frac{\Delta t_{\ell-1}}{\varepsilon^{2}}$ we can then write the r.h.s. of (117) as

[TABLE]

To calculate the covariance of coarse time steps preceding fine time steps (IV), fewer calculations are needed. The two time steps are correlated if the trajectories were correlated at the start of coarse simulation time step $n$ and this correlation was not lost due to a collision in the fine simulation between time steps $(n,0)$ through $(n^{\prime},m^{\prime}-1)$

[TABLE]

Making use of expressions (113) through (124), we get

[TABLE]

$\text{with}\qquad\sigma_{1}=\varepsilon^{2}\tilde{v}^{2}\left(p_{nc,\Delta t_{\ell}}^{1-M}+p_{nc,\Delta t_{\ell}}^{M+1}-2p_{nc,\Delta t_{\ell}}\right)-\Delta t_{\ell-1}^{2}\tilde{v}_{\Delta t_{\ell-1}}^{2}\qquad\text{and}$

[TABLE]

Appendix B Theoretical background to Section 6.2

Based on Section 2.6 in Giles2015 it can be shown that it is beneficial to leave out level 0 if

[TABLE]

with $C_{\ell}$ the cost of a sample at level $\ell$. Our analytical results from Section 6.1 can only strictly be applied in (130) if the quantity of interest is the expected particle position or velocity, as applying Lipschitz constants gives independent upper bounds for the positive and negative terms, hence we will also back up our claims with numerical experiments. Here, we only consider the particle position variance, as $\operatorname{\mathbb{V}}\left[X_{\ell}-X_{\ell-1}\right]\gg\operatorname{\mathbb{V}}\left[V_{\ell}-V_{\ell-1}\right]$ and $\operatorname{\mathbb{V}}\left[X_{\ell}\right]\gg\operatorname{\mathbb{V}}\left[V_{\ell}\right]$, for large $\Delta t_{\ell}$.

We first fix $\Delta t_{1}=\varepsilon^{2}$ and solve (130) for positive real values of $M>1$ for which the inequality becomes an equality, taking different values of $\varepsilon$ and $t^{*}$. In Figure 10 we show these unique, numerically computed $M$-values. It is clear from the range of the color bar that this $M$-value consistently exists and lies within the interval $[6,13]$, meaning that any sign change in the r.h.s. of (130) lies in this range. Next, we evaluate the left hand side of (130) for $M=6$ and $M=13$ and plot the results in Figure 11. From this figure, it is clear that (130) holds for small values of $M$, but no longer holds once the threshold value from Figure 10 is passed, meaning that (130) consistently holds once $M$ is sufficiently large. We can thus conclude that it makes sense, from a theoretical point of view, to add an extra coarse level, if it is sufficiently coarse in comparison with the level with $\Delta t=\varepsilon^{2}$.