A Multilevel Monte Carlo Asymptotic-Preserving Particle Method for Kinetic Equations in the Diffusion Limit

TL;DR

This paper introduces a multilevel Monte Carlo particle method that efficiently simulates kinetic equations in the diffusion limit, overcoming stability and accuracy issues of classical techniques.

Contribution

It combines multilevel Monte Carlo with asymptotic-preserving schemes to reduce bias and computational cost in simulating kinetic equations in the diffusive regime.

Findings

Effective bias reduction in diffusive limit simulations

Reduced computational cost compared to traditional methods

Successful numerical validation of the approach

Abstract

We propose a multilevel Monte Carlo method for a particle-based asymptotic-preserving scheme for kinetic equations. Kinetic equations model transport and collision of particles in a position-velocity phase-space. With a diffusive scaling, the kinetic equation converges to an advection-diffusion equation in the limit of zero mean free path. Classical particle-based techniques suffer from a strict time-step restriction to maintain stability in this limit. Asymptotic-preserving schemes provide a solution to this time step restriction, but introduce a first-order error in the time step size. We demonstrate how the multilevel Monte Carlo method can be used as a bias reduction technique to perform accurate simulations in the diffusive regime, while leveraging the reduced simulation cost given by the asymptotic-preserving scheme. We describe how to achieve the necessary correlation between simulation paths at different levels and demonstrate the potential of the approach via numerical experiments.

Tables

Table 1. Table 1: Results of the simulation in Section 4.2 with a geometric level sequence for E = 0.1 𝐸 0.1 E=0.1 .

Level	$\Delta t_{\ell }$	$P_{\ell }$	$𝕍\left[F_{\ell }\right]$	$𝔼\left[F_{\ell }-F_{\ell -1}\right]$	$V_{\ell }$	$𝕍\left[{\hat{Y}}_{\ell }\right]$	$C_{\ell }$	$P_{\ell }{C}_{\ell }$
0	$1.00\times {10}^{-2}$	1 393	1.32	$8.18\times {10}^{-1}$	$1.32\times {10}^0$	$9.45\times {10}^{-4}$	1	1 393
1	$5.00\times {10}^{-3}$	395	1.52	$7.91\times {10}^{-3}$	$3.58\times {10}^{-1}$	$9.07\times {10}^{-4}$	3	1 185
2	$2.50\times {10}^{-3}$	296	1.59	$2.18\times {10}^{-2}$	$4.82\times {10}^{-1}$	$1.59\times {10}^{-3}$	6	1 776
3	$1.25\times {10}^{-3}$	229	2.22	$-1.48\times {10}^{-2}$	$3.22\times {10}^{-1}$	$1.41\times {10}^{-3}$	12	2 748
4	$6.25\times {10}^{-4}$	40	1.70	$1.57\times {10}^{-3}$	$4.56\times {10}^{-2}$	$1.14\times {10}^{-3}$	24	960
$\sum$						$6.00\times {10}^{-3}$		8 062

Table 2. Table 2: Results of the simulation in Section 4.2 with a geometric level sequence for E = 0.01 𝐸 0.01 E=0.01 .

Level	$\Delta t_{\ell }$	$P_{\ell }$	$𝕍\left[F_{\ell }\right]$	$𝔼\left[F_{\ell }-F_{\ell -1}\right]$	$V_{\ell }$	$𝕍\left[{\hat{Y}}_{\ell }\right]$	$C_{\ell }$	$P_{\ell }{C}_{\ell }$
0	$1.00\times {10}^{-2}$	527 920	1.47	$8.65\times {10}^{-1}$	$1.47\times {10}^0$	$2.79\times {10}^{-6}$	1	527 920
1	$5.00\times {10}^{-3}$	165 386	1.49	$1.06\times {10}^{-2}$	$4.35\times {10}^{-1}$	$2.63\times {10}^{-6}$	3	496 158
2	$2.50\times {10}^{-3}$	112 208	1.59	$2.98\times {10}^{-2}$	$3.99\times {10}^{-1}$	$3.55\times {10}^{-6}$	6	673 248
3	$1.25\times {10}^{-3}$	69 135	1.64	$2.84\times {10}^{-2}$	$3.01\times {10}^{-1}$	$4.36\times {10}^{-6}$	12	829 620
4	$6.25\times {10}^{-4}$	39 146	1.73	$2.00\times {10}^{-2}$	$1.95\times {10}^{-1}$	$4.98\times {10}^{-6}$	24	939 504
5	$3.13\times {10}^{-4}$	20 670	1.76	$7.53\times {10}^{-3}$	$1.09\times {10}^{-1}$	$5.28\times {10}^{-6}$	48	992 160
6	$1.56\times {10}^{-4}$	10 842	1.75	$9.55\times {10}^{-3}$	$6.14\times {10}^{-2}$	$5.67\times {10}^{-6}$	96	1 040 832
7	$7.81\times {10}^{-5}$	4 894	1.91	$6.77\times {10}^{-3}$	$2.42\times {10}^{-2}$	$4.94\times {10}^{-6}$	192	939 648
8	$3.91\times {10}^{-5}$	3 937	1.77	$2.88\times {10}^{-3}$	$1.21\times {10}^{-2}$	$3.08\times {10}^{-6}$	384	1 511 808
9	$1.95\times {10}^{-5}$	2 721	1.81	$2.35\times {10}^{-3}$	$1.21\times {10}^{-2}$	$4.46\times {10}^{-6}$	768	2 089 728
10	$9.75\times {10}^{-6}$	40	1.47	$2.00\times {10}^{-3}$	$4.35\times {10}^{-4}$	$1.09\times {10}^{-5}$	1 536	61 440
$\sum$						$5.26\times {10}^{-5}$		10 102 066

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

RMSE	Classical cost	Multilevel cost	Speedup
0.1	4 544	8 062	0.56
0.01	28 627 968	10 102 066	2.83
0.001	25 167 200 256	1 742 633 519	14.4

Table 4. Table 5: Results of the simulation in Section 4.2 with an extra coarse level for E = 0.01 𝐸 0.01 E=0.01 .

Level	$\Delta t_{\ell }$	$P_{\ell }$	$𝕍\left[F_{\ell }\right]$	$𝔼\left[F_{\ell }-F_{\ell -1}\right]$	$V_{\ell }$	$𝕍\left[{\hat{Y}}_{\ell }\right]$	$C_{\ell }$	$P_{\ell }{C}_{\ell }$
0	$5.00\times {10}^{-1}$	2 978 687	1.96	$9.91\times {10}^{-1}$	$1.96\times {10}^0$	$6.60\times {10}^{-7}$	0.02	59 574
1	$1.00\times {10}^{-2}$	354 282	1.47	$-1.26\times {10}^{-1}$	$1.42\times {10}^0$	$4.01\times {10}^{-6}$	1.02	361 368
2	$5.00\times {10}^{-3}$	114 863	1.47	$1.06\times {10}^{-2}$	$4.36\times {10}^{-1}$	$3.79\times {10}^{-6}$	3	344 589
3	$2.50\times {10}^{-3}$	77 905	1.57	$3.14\times {10}^{-2}$	$4.01\times {10}^{-1}$	$5.15\times {10}^{-6}$	6	467 430
4	$1.25\times {10}^{-3}$	47 439	1.68	$2.92\times {10}^{-2}$	$3.00\times {10}^{-1}$	$6.32\times {10}^{-6}$	12	569 268
5	$6.25\times {10}^{-4}$	27 466	1.78	$2.57\times {10}^{-2}$	$2.01\times {10}^{-1}$	$7.32\times {10}^{-6}$	24	659 184
6	$3.13\times {10}^{-4}$	14 599	1.75	$9.28\times {10}^{-3}$	$1.10\times {10}^{-1}$	$7.56\times {10}^{-6}$	48	700 752
7	$1.56\times {10}^{-4}$	7 666	1.71	$3.99\times {10}^{-3}$	$6.74\times {10}^{-2}$	$8.79\times {10}^{-6}$	96	735 936
8	$7.81\times {10}^{-5}$	3 195	2.04	$4.04\times {10}^{-3}$	$2.49\times {10}^{-2}$	$7.80\times {10}^{-6}$	192	613 440
9	$3.91\times {10}^{-5}$	40	1.54	$5.51\times {10}^{-3}$	$3.65\times {10}^{-3}$	$9.12\times {10}^{-5}$	384	15 360
$\sum$						$1.42\times {10}^{-4}$		4 526 900

Equations

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

$$\rho (x,t)=\int f(x,v,t)dv.$$

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

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

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

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

$${\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+V_{p,\Delta t}^n\Delta t.$$

$$V_{p,\Delta t}^{n+1}=\{\begin{array}{ll}\pm 1/ϵ,& \text{with probability }{p}_{c,\Delta t}=\Delta t/{ϵ}^2\text{ and equal probability in the sign},\\ V_{p,\Delta t}^n,& \text{otherwise.}\end{array}$$

$$\{\begin{array}{l}{\partial }_t{f}_{+}+\frac{ϵ}{{ϵ}^2+\Delta t}{\partial }_x{f}_{+}=\frac{\Delta t}{{ϵ}^2+\Delta t}{\partial }_{xx}{f}_{+}+\frac{1}{{ϵ}^2+\Delta t}\left(\frac{\rho }{2}-f_{+}\right)\\ {\partial }_t{f}_{-}-\frac{ϵ}{{ϵ}^2+\Delta t}{\partial }_x{f}_{-}=\frac{\Delta t}{{ϵ}^2+\Delta t}{\partial }_{xx}{f}_{-}+\frac{1}{{ϵ}^2+\Delta t}\left(\frac{\rho }{2}-f_{-}\right)\end{array}.$$

$$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,\Delta t}^n=\pm \frac{ϵ}{{ϵ}^2+\Delta t},D_{\Delta t}=\frac{\Delta t}{{ϵ}^2+\Delta t}.$$

$$V_{p,\Delta t}^{n+1}=\{\begin{array}{ll}\pm \frac{ϵ}{{ϵ}^2+\Delta t},& \text{with probability }{p}_{c,\Delta t}=\frac{\Delta t}{{ϵ}^2+\Delta t}\\ & \text{blabla and equal probability in the sign},\\ V_{p,\Delta t}^n,& \text{otherwise.}\end{array}$$

$$Y(t^{\ast })=\mathbb{E}[F(X(t^{\ast }))]=\int \int F(x)f(x,v,t^{\ast })dxdv.$$

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

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

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

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

$$\{\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,{\xi }_{p,\ell -1}^n\sim \mathcal{N}(0,1),\\ X_{p,\Delta t_{\ell }}^{n+1,0}=X_{p,\Delta t_{\ell }}^{n,0}+\sum _{m=1}^M\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),{\xi }_{p,\ell }^{n,m}\sim \mathcal{N}(0,1),\end{array}$$

$${\xi }_{p,\ell -1}^n=\sum _{m=1}^M\frac{{\xi }_{p,\ell }^{n,m}}{\sqrt{M}}.$$

$${\alpha }_{p,\ell }^{n,m}\ge p_{nc,\Delta t_{\ell }}=\frac{{ϵ}^2}{{ϵ}^2+\Delta t_{\ell }}.$$

$${\alpha }_{p,\ell }^{n,\text{max}}=\underset{m}{\max }{\alpha }_{p,\ell }^{n,m},$$

$$\text{CDF}\left({\alpha }_{p,\ell }^{n,\text{max}}\right)={\left({\alpha }_{p,\ell }^{n,\text{max}}\right)}^M.$$

$${\alpha }_{p,\ell -1}^n={\left({\alpha }_{p,\ell }^{n,\text{max}}\right)}^M,$$

$$\text{sign}\left(V_{p,\Delta t_{\ell -1}}^{n+1}\right)=\text{sign}\left(V_{p,\Delta t_{\ell }}^{n,i}\right),i=\underset{1\le m\le M}{\mathrm{argmax}}\left(m|{\alpha }_{p,\ell }^{n,m}\ge p_{nc,\Delta t_{\ell }}\right).$$

$$\{\begin{array}{l}{\partial }_t{f}_{+}(x,t)={\partial }_{xx}{f}_{+}(x,t)\\ {\partial }_t{f}_{-}(x,t)={\partial }_{xx}{f}_{-}(x,t)\end{array}\Rightarrow {\partial }_t\rho (x,t)={\partial }_{xx}\rho (x,t).$$

$$\lceil 2E^{-2}\sqrt{\frac{V_{\ell }}{C_{\ell }}}\left(\sum _{\ell =0}^L\sqrt{V_{\ell }{C}_{\ell }}\right)\rceil ,$$

$$\mathbb{V}\left[{\hat{Y}}_{\ell }\right]=\frac{V_{\ell }}{P_{\ell }}.$$

$$P_C=\lceil \frac{\mathbb{V}\left[F_L\right]}{{\sum }_{\ell =0}^L\mathbb{V}\left[{\hat{Y}}_{\ell }\right]}\rceil$$

Full text

1 Introduction

Kinetic equations, modeling particle behavior in a position-velocity phase space, occur in many domains. Examples are plasma physics Birdsall2004 , bacterial chemotaxis Rousset2011c and computational fluid dynamics Pope1981 . Many of these applications exhibit a strong time-scale separation, leading to an unacceptably high simulation cost Cercignani1988 . However, one typically is only interested in computing the evolution of some macroscopic quantities of interest. These are usually some moments of the particle distribution, which can be computed as averages over velocity space. The time-scale at which these quantities of interest change is often much slower than the time-scale governing the particle dynamics. The nature of the macroscopic dynamics depends on the scaling of the problem, which can be either hyperbolic or diffusive Dimarco2014 .

The model problem in this work is a one-dimensional kinetic equation of the form

[TABLE]

where $f(x,v,t)$ represents the distribution of particles as a function of position $x\in{\mathbb{R}}$ and velocity $v\in{\mathbb{R}}$ as it evolves in time $t\in{\mathbb{R}}^{+}$. The left-hand side of (1) represents transport, while $Q(f(x,v,t))$ is a collision operator that results in discontinuous velocity changes. As the collision operator, we take the BGK model Bhatnagar1954 , which represents linear relaxation to an equilibrium distribution that only depends on the particle density

[TABLE]

We introduce a parameter $\epsilon$ that represents the mean free path. When decreasing $\epsilon$, the average time between collisions decreases. In this paper, we consider the diffusive scaling. In that case, we simultaneously increase the time scale at which we observe the evolution of the particle distribution, arriving at

[TABLE]

with $\mathcal{M}(v)$ the particles’ steady state velocity distribution. It has been shown that when taking the limit $\epsilon\to 0$, the behavior of equations of the form (3) is fully described by the diffusion equation Lapeyre2003

[TABLE]

Kinetic equations can be simulated with deterministic methods, solving the partial differential equation (PDE) that describes the evolution of the particle distribution in the position-velocity phase space. Alternatively, one can use stochastic methods that simulate a large number of particle trajectories. Deterministic methods become prohibitively expensive for higher dimensional applications. Particle-based methods do not suffer from this curse of dimensionality, at the expense of introducing a statistical error in the computed solution. The issue of time-scale separation is present in both deterministic and stochastic methods.

One way to avoid the issue of time-scale separation is through the use of asymptotic-preserving methods, which aim at reproducing a scheme for the limiting macroscopic equation in the limit of infinite time-scale separation. For deterministic discretization methods, there is a long line of such methods. We refer to Bennoune2008 ; Boscarino2013 ; Buet2007 ; Crouseilles2011 ; Dimarco2012 ; Gosse2002 ; Jin1999 ; Jin1998 ; Jin2000 ; Klar1998 ; Klar1999 ; Larsen1974 ; Lemou2008 ; Naldi2000 as a representative sample of such methods in the diffusive scaling. The recent review paper Dimarco2014 contains an overview of the state of the art on asymptotic-preserving methods for kinetic equations, and ample additional references. In the particle-based setting, only a few asymptotic-preserving methods have been developed, mostly in the hyperbolic scaling Degond2011 ; Dimarco2008 ; Dimarco2010 ; Pareschi1999 ; Pareschi2001 ; Pareschi2005 . In the diffusive scaling, there are only two works Crestetto2018 ; Dimarco2018 so far, to the best of our knowledge. Both methods avoid the time step restrictions caused by fast problem time-scales, at the expense of introducing a bias, which is of order one in the time step size.

The goal of the present paper is to combine the asymptotic-preserving scheme in Dimarco2018 with the multilevel Monte Carlo method. Given a fixed computational budget, a trade-off typically has to be made between a small bias and a low variance. The former can be obtained by reducing the time step, the latter by simulating many trajectories with large time steps. The core idea behind the multilevel Monte Carlo method Giles2008 is to reduce computational cost, by combining estimates computed with different time step sizes. The multilevel Monte Carlo method, originally developed in the context of stochastic processes, has been applied to problems across many fields, for example, finance Giles2008 and biochemistry Anderson2011 . The method has successfully been applied to simulating large PDE’s with random coefficients Cliffe2011 . Recent work has also used multilevel Monte Carlo methods in an optimization context VanBarel2019 .

The remainder of this paper is organized as follows. In Section 2, we describe the model kinetic equation on which we will demonstrate our approach, as well as the asymptotic-preserving Monte Carlo scheme that was introduced in Dimarco2018 . In Section 3, we cover the multilevel Monte Carlo method that is the core contribution of this paper. In Section 4, we present some preliminary experimental results, demonstrating the properties of the new scheme as well as its computational gain. Finally, in Section 5 we will summarize our main results and mention some possible future extensions.

2 Model problem and asymptotic-preserving scheme

2.1 Model equation in the diffusive limit

The model problem considered in this work is a one-dimensional kinetic equation in the diffusive scaling of the form (3), which we rewrite as

[TABLE]

For ease of exposition, we restrict ourselves to the case of two discrete velocities, $v=\pm 1$. Then, we can write $f_{+}(x,t)$ and $f_{-}(x,t)$ to represent the distribution of particles with, respectively, positive and negative velocities, and $\rho(x,t)=f_{+}(x,t)+f_{-}(x,t)$ represents the total density of particles. In this case, equation (5) simplifies to

[TABLE]

Equation (6) is also known as the Goldstein-Taylor model, and can be solved using a particle scheme. For this, we introduce a time step $\Delta t$ and an ensemble of $P$ particles

[TABLE]

The particle state (position and velocity) is represented as $(X,V)$, $p$ is the particle index ($1\leq p\leq P$), and $n$ represents the time index, i.e., $X_{p,\Delta t}^{n}\approx X_{p}(n\Delta t)$. Equation (6) is then solved via operator splitting as

Transport step. The position of each particle is updated based on its velocity

[TABLE] 2. 2.

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

[TABLE]

This approximation requires a time step restriction $\Delta t=\mathcal{O}(\epsilon^{2})$ as $\epsilon\rightarrow 0$, both to ensure $p_{c,\Delta t}<1$ in the collision phase, and to keep the increments in the transport phase finite. This leads to unacceptably high computational costs for small $\epsilon$.

2.2 Asymptotic-preserving Monte Carlo scheme

Recently, an asymptotic-preserving Monte Carlo scheme was proposed Dimarco2018 , based on the simulation of a modified equation

[TABLE]

In (10) we have dropped the space and time dependency of $f_{\pm}$ and $\rho$, for conciseness. The model given by (10) reduces to (6) in the limit when $\Delta t$ tends to zero and has an $\mathcal{O}(\Delta t)$ bias. In the limit when $\epsilon$ tends to zero, the equations reduce to (4).

Discretizing this equation, using operator splitting as above, again leads to a Monte Carlo scheme. For each particle $X_{p}$ and for each time step $n$, one time step now consists of a transport-diffusion and a collision step:

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 $V^{n}_{p,\Delta t}$ and diffusion coefficient $D_{\Delta t}$:

[TABLE] 2. 2.

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

[TABLE]

For more details, we refer the reader to Dimarco2018 .

3 Multilevel Monte Carlo method

3.1 Method and notation

We want to estimate some quantity of interest $Y$ that is a function of the particle distribution $f(x,v,t)$ at some specific moment $t=t^{*}$ in time, i.e., we are interested in

[TABLE]

Note that, in equation (14), the function $F$ only depends on the position $x$ and not on velocity. This is a choice we make for notational convenience and is not essential for the method we present.

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

[TABLE]

Here, $P$ denotes the number of simulated trajectories, $N$ the number of simulated time steps, $\Delta t$ the time step size, and $X^{N}_{p,\Delta t}$ is generated by the time-discretised process (11)–(13). Given a constrained computational budget, a trade-off has to be made when selecting the time step size $\Delta t$. On the one hand, a small time step reduces the bias of the simulation of each sampled trajectory, and thus of the estimated quantity of interest. On the other hand, a large time step reduces the cost per trajectory, which increases the number of trajectories that can be simulated and thus reduces the resulting variance on the estimate. The key idea behind the Multilevel Monte Carlo method Giles2008 is to generate a sequence of estimates with varying discretization accuracy and a varying number of realizations. The method achieves the bias of the finest discretization, with the variance of the coarsest discretization.

To apply the multilevel Monte Carlo method, we define a sequence of time step sizes, denoted by $\Delta t_{\ell}$ with $\ell=0\dots L$, with $\ell=L$ denoting the finest level of discretization (smallest time step), and $\ell=0$ the coarsest level. We use a fixed ratio of time steps between subsequent levels, i.e., we set $\Delta t_{\ell-1}=M\Delta t_{\ell}$ for some integer $M$. At each level, we simulate a number $P_{\ell}$ of particle trajectories. An initial coarse estimator with a large number $P_{0}$ of sample trajectories is given by

[TABLE]

This initial estimate can be improved upon by a series of difference estimators $\hat{Y}_{\ell}(t^{*})$, $\ell=1\dots L$, of the form

[TABLE]

with $N_{\ell}\Delta t_{\ell}=t^{*}$, for each value of $\ell$, and $P_{\ell}$ the number of correlated sample trajectories at each level. The estimators (17) estimate the bias induced by sampling with a simulation time step size $\Delta t_{\ell-1}$ by comparing the sample results with a simulation using a time step size $\Delta t_{\ell}$. The estimators (16)-(17) are then combined into a multilevel Monte Carlo estimator via a telescopic sum,

[TABLE]

It can easily be seen that the expected value of estimator (18) is the same as that of estimator (15) with the finest time step $\Delta t_{L}$. If the required number of particles $P_{\ell}$ at each level decreases sufficiently fast with increasing level $\ell$, the multilevel estimator will result in a reduced computational cost for a given accuracy. For more details on the multilevel Monte Carlo method, we refer to Giles2015 .

3.2 Correlating asymptotic-preserving Monte Carlo simulations

Coupled trajectories and notation

The differences in (17) will only have low variance if the simulated paths $X^{n,m}_{\Delta t_{\ell},p}$ and $X^{n}_{\Delta t_{\ell-1},p}$ are correlated. To achieve this correlation, we will couple the different sources of randomness in the simulation at consecutive levels. In each time step using the asymptotic-preserving particle scheme (11)–(13), 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 (11). On the other hand, in each collision step (13), a fraction of particles randomly get a new velocity $V_{p}^{n}$.

Particle trajectories can be coupled by separately correlating the random numbers used for the individual particles in the transport-diffusion and collision phase of each time step. To show how this is done, we introduce a 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 $m\in\{1,\dots,M\}$ and $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})$.

The key point of the algorithm is to compute the velocities $V^{n}_{p,\Delta t_{\ell-1}}$ and the Brownian increments $\xi^{n}_{p,\ell-1}$ at level $\ell-1$, based on the randomly generated values $\xi^{n,m}_{p,\ell}$ and $V^{n,m}_{p,\Delta t_{\ell}}$ at level $\ell$, instead of generating these independently. The main difficulty lies in maximizing the correlation between the velocities and Brownian increments at levels $\ell$ and $\ell-1$ while avoiding the introduction of an extra bias at level $\ell-1$. Once the coupled simulation (19) at level $\ell-1$ is performed, we can insert the results in (17) to obtain a low-variance difference estimator. In the next two subsections, we explain how we correlate the Brownian increments during the transport phase (Section 3.2) and the velocities during the collision phase (Section 3.2). We present the complete algorithm in Section 3.2.

Coupling the transport-diffusion phase

We first correlate the Brownian increments at levels $\ell$ and $\ell-1$. To this end, we first simulate the stochastic process at level $\ell$, using i.d.d. increments $\xi^{n,m}_{p,\ell}$. Then, at level $\ell-1$, we compute the Brownian increments, $\xi^{n}_{p,\ell-1}$, from those at level $\ell$, $\left\{\xi^{n,m}_{p,\ell}\right\}$, ensuring that $\xi^{n}_{p,\ell-1}\sim\mathcal{N}(0,1)$. This condition is clearly satisfied if we define $\xi^{n}_{p,\ell-1}$ as

[TABLE]

Correlating the simulations in this way means that both levels use the same Brownian path, and differences in the diffusion part of the motion only result from differences in the diffusion coefficients $D_{\ell}$ and $D_{\ell-1}$ at different levels.

In Figure 1, we show two particle trajectories, containing only diffusion behavior, i.e., (19) with $V_{p,\Delta t_{\ell-1}}^{n}=V_{p,\Delta t_{\ell}}^{n}=0$, coupled as described in (20) with $\epsilon=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 the bias caused by the paths having different diffusion coefficients.

Coupling the collision phase

While correlating the Brownian paths is relatively straightforward, the coupling of the velocities in the collision phase is more involved. Since we simulate level $\ell$ first, we have at our disposal the velocities $V^{n,m}_{p,\Delta t_{\ell}}$ at level $\ell$, which are again i.i.d. Our goal is to compute the velocities $V^{n}_{p,\Delta t_{\ell-1}}$ at level $\ell-1$ from those at level $\ell$, to maximize correlation, while ensuring that the collision probability and post-collision velocity distribution at level $\ell-1$ are satisfied. Note that, in the collision phase of the asymptotic-preserving particle scheme (13), 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 computation of $V^{n}_{p,\Delta t_{\ell-1}}$ is done in two steps. First, we will couple the occurrence of a collision at level $\ell-1$ to the occurrence of a collision in one of the $M$ sub-steps of the correlated fine simulation. If we decide to perform a collision both at level $\ell$ and $\ell-1$, we will correlate the new velocities generated in both simulations.

Let us first consider the simulation at level $\ell$. When simulating the collision step, we decide whether a collision has occurred during a time step of length $\Delta t_{\ell}$ by drawing a random number $\alpha^{n,m}_{p,\ell}\sim\mathcal{U}([0,1])$ and comparing it to the probability that no collision has occurred in the simulation, $p_{nc,\Delta t_{\ell}}=1-p_{c,\Delta t_{\ell}}$, with $p_{c,\Delta t_{\ell}}$, defined in equation (13). A collision takes place if and only if

[TABLE]

Now consider $M$ time steps of length $\Delta t_{\ell}$. At least one collision has taken place if at least one of the generated $\alpha_{p,\ell}^{n,m}$, $m\in\{1,\dots,M\}$, satisfies (21).

Deciding upon collision in the coarse simulation At level $\ell-1$, we want to use the values $\alpha_{p,\ell}^{n,m}$, $m\in\{1,\dots,M\}$ to compute a uniformly distributed number $\alpha_{p,\ell-1}^{n}$, that is correlated with the largest of the generated $\alpha_{p,\ell}^{n,m}$

[TABLE]

to compare with the collision probability $p_{nc,\Delta t_{\ell-1}}$. However, the maximum of a set of uniformly distributed random number is not uniformly distributed. The cumulative density function of $\alpha^{n,\text{max}}_{p,\ell}$ is given by

[TABLE]

Hence, by the inverse transform method, ${\left(\alpha^{n,\text{max}}_{p,\ell}\right)}^{M}\sim\mathcal{U}([0,1])$. Equation (23) implies that we can define this random number as

[TABLE]

without affecting the simulation statistics at level $\ell-1$.

It is possible to show that, given the relation in (24), a collision can occur in the fine simulation without a collision occurring in the coarse simulation. The inverse, i.e., a collision in the coarse simulation, without a fine simulation collision, is not possible.

Choosing a new velocity. If a collision takes place in both simulations in a given time step $\Delta t_{\ell-1}$, then we set the sign of the velocity of the coarse simulation, at the end of the time step to be equal in sign to the velocity of the last subdividing fine time step for which (21) holds,

[TABLE]

Because the new velocities generated in the fine simulation are i.i.d., we are free to make this selection, without altering the statistics of the coarse simulation. This approach to selecting the sign of $V_{p,\Delta t_{\ell-1}}^{n+1}$ means that the velocities going into the next time step will have the same sign.

Two particle trajectories without diffusion behavior, i.e., (19) with $D_{\Delta t_{\ell-1}}=D_{\Delta t_{\ell}}=0$ are shown in Figure 2. In this figure, a number of interesting phenomena can be observed. First of all, the fact that the particle’s characteristic velocity is dependent on the time step sizes $\Delta t_{\ell-1}$ and $\Delta t_{\ell}$ results in different slopes in the curves. This is one source of the bias that we want to estimate using the multilevel Monte Carlo method. Second of all, the collision probability between the coupled trajectories does not match precisely, as this probability also depends on $\Delta t_{\ell-1}$ and $\Delta t_{\ell}$. 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. By coincidence, the new velocity generated at $t=8$ in the fine simulation has the same sign as the coarse simulation velocity. This mismatch is also part of the bias we wish to estimate.

The complete algorithm

Combining the correlation of the Brownian increments and velocities results in Algorithm 1. The correlation of the trajectories can be seen in Figure 3 which shows the particle trajectory given by the sum of the behaviors in Figures 1 and 2.

4 Experimental Results

We will now demonstrate the viability of the suggested approach through some numerical experiments. We will simulate the model given by (10), using the multilevel Monte Carlo method to estimate a selected quantity of interest, which is the expected value of the square of the particle position, at $t^{*}$. The ensemble of particles is initialized at the origin with equal probability of having a left and right velocity. When discussing results we will replace the full expression for a sample of the quantity of interest, based on an arbitrary particle $p$, $F(X^{N,0}_{\Delta t_{\ell},p})$, with the symbol $F_{\ell}$ to simplify notation.

4.1 Model correlation behavior

In a first test, we set $t^{*}=5$ and investigate the variance of the difference estimators (17) as a function of the time step $\Delta t_{\ell}$ (or, equivalently) the level number. At level $\ell=0$, we set $\Delta t_{0}=2.5$. All finer levels ($\ell\geq 1$) are defined by setting $\Delta t_{\ell}=\Delta t_{\ell-1}/M$ with $M=2$. We fix the number of samples per difference estimator at 100 000. For a selection of values of $\epsilon$, we calculate the expected value and variance as a function of $\Delta t_{\ell}$, for $1\leq\ell$. We compute both the variance of the function samples for a given $\Delta t_{\ell}$, and the variance of the sampled differences (17), based on the coupled trajectories computed using $\Delta t_{\ell-1}$ and $\Delta t_{\ell}$. We choose $\epsilon=10$ (Figure 4), $\epsilon=1$ (Figure 5), $\epsilon=0.1$ (Figure 6) and $\epsilon=0.01$ (Figure 7).

The regime $\Delta t\ll\epsilon^{2}$. In Figures 4 through 6, we see that 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\epsilon^{2}$. This matches the weak convergence order of the Euler-Maruyama scheme, used to simulate the model (11)–(13), as well as the expected behavior from the time step dependent bias in the asymptotic-preserving model. Given this asymptotic geometric convergence, it is possible to apply the complexity theorem in Giles2008 to analyze the method’s computational cost and error bounds. This means that existing theory for multilevel Monte Carlo methods Giles2015 concerning, e.g. samples per level, convergence criteria and conditions for adding levels, can be applied in this regime.

The regime $\Delta t\gg\epsilon^{2}$. For time steps $\Delta t\gg\epsilon^{2}$, however, we see in Figures 6 and 7 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 (10) converges to the heat equation:

[TABLE]

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

The fact that the two limits approach different models can be seen most clearly in Figures 4 and 6. In the right hand panel of Figure 4 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 $\epsilon^{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 (6). For large $\Delta t$ there is convergence to (4). In practice, the size of $\Delta t$ is limited by the simulation time horizon, so it is not possible to get arbitrarily close to (4) by increasing the time step size, however. This phenomenon also has an effect on the curves in Figure 6. The curves for the mean and variance of the differences $F_{\ell}-F_{\ell-1}$ (orange lines with dots) decrease for both small and large $\Delta t$, as the model converges to the two limits.

Combining the observations from the two limits 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 $\epsilon$. The peak of the variance of the differences lies near $\Delta t\approx\epsilon^{2}$. This makes sense, as this is the region where the model parameters $D_{\Delta t}$ and $V^{n}_{p,\Delta t}$ vary the most in function of $\Delta t$. We also see a dip in the mean of the difference curves in the region of $\Delta t\approx\epsilon^{2}$. A full analysis of the behavior that occurs in the transition between the asymptotic regimes is left for future work.

4.2 Comparison with classical Monte Carlo

The analysis in Section 4.1 demonstrated a fast decay of the variance of the differences for increasingly fine levels in the region where $\Delta t\ll\epsilon^{2}$. As such, one of the necessary requirements for convergence of the multilevel Monte Carlo method is present in this region. This is, however, not the case in the regime where $\Delta t\gg\epsilon^{2}$. Here, the variance of the differences increases as the time step is refined. It is therefore highly non-trivial to perform an adequate selection of coarse levels in the regime $\Delta t\gg\epsilon^{2}$. For the fine levels, a standard multilevel Monte Carlo approach can be applied. We therefore propose two simulation strategies:

A geometric sequence of levels $\Delta t_{\ell}=\epsilon^{2}M^{-\ell}$ for $\ell>0$ starting with a coarse simulation time step $\Delta t_{0}=\epsilon^{2}$; 2. 2.

The same geometric sequence, preceded by a coarse simulation time step $t^{*}$, i.e., $\Delta t_{0}=t^{*}$, $\Delta t_{1}=\epsilon^{2}$ and $\Delta t_{\ell}=\epsilon^{2}M^{1-\ell}$ for $\ell>1$.

We compare these approaches in the following two sub-sections.

Standard MLMC refinement

We will now compute the quantity of interest described at the beginning of this section to a range of prescribed error tolerances, to verify the reduced computational cost of the multilevel Monte Carlo method. We choose to set $M=2$ and $\epsilon=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]

where $E$ is the desired root mean square error, $C_{\ell}$ is the computational cost of the estimator at level $\ell$, and $V_{\ell}$ is the estimated variance of the estimator at level $\ell$, i.e., $V_{\ell}=\mathbb{V}\left[F_{\ell}-F_{\ell-1}\right]$, where we set $F_{-1}\equiv 0$. The criterion for adding levels and determining convergence are as described in Giles2015 . The cost of a sample will be determined relative to the cost of a simulated trajectory with $\Delta t=\epsilon^{2}$. The results of the simulations for $E$ values 0.1, 0.01 and 0.001 can be found 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[F_{\ell}\right]$, expected value $\mathbb{E}\left[F_{\ell}-F_{\ell-1}\right]$ and variance $V_{\ell}$ of the differences of simulations, estimated variance of the estimator $\mathbb{V}[\hat{Y}_{\ell}]$, cost per sample $C_{\ell}$ and cost per level $P_{\ell}C_{\ell}$. The variance of the estimator at level $\ell$ is estimated as

[TABLE]

We see that the experimental results match the expected behavior of the multilevel Monte Carlo method. 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^{2}$. The cost per level $P_{\ell}C_{\ell}$ is also spread quite evenly over the levels, once the time step is a couple orders of magnitude smaller than $\epsilon^{2}$. This 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.

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 cost of each level. 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 the numbers in Table 4, $\mathbb{V}\left[F_{L}\right]$ is estimated using very few samples, so these results should not be taken to literally. They do give the correct order of magnitude of the cost of the equivalent classical Monte Carlo method, however. We now compare the cost of the classical and multilevel Monte Carlo simulations in Table 4.

As can be concluded from Table 4 the multilevel Monte Carlo method 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.

Adding a coarse level

It makes little sense to add a full sequence of levels in the regime where $\Delta t\gg\epsilon^{2}$. It could however make sense to add a single very coarse level to the simulation as the variance of $F_{\ell}$ is consistently larger than that of $F_{\ell}-F_{\ell-1}$ in Figures 4 through 7. To test this idea we repeat the experiment as before, for $E=0.01$, with a coarse level at $\Delta t_{0}=0.5$. The results of this experiment can be seen in Table 5

We see that the total cost of the simulation with the extra coarse level is lower than that of the simulation starting with $\Delta t=\epsilon^{2}$ (4 526 900 as apposed to 10 102 066). Based on this initial experiment, it makes sense to include a very coarse level when using the multilevel Monte Carlo method in this context. For a detailed analysis and more extensive numerical results, we refer to future work.

5 Conclusion

In this work, we have derived a new multilevel scheme for asymptotic-preserving particle schemes of the form given in (10). We have demonstrated that this scheme has interesting convergence behavior as the time step is refined, which is apparent in the expected value and variance of sampled differences of the quantity of interest. On the one hand, we get the expected linear convergence to the exact model in terms of $\Delta t$ for a fixed value of $\epsilon$. On the other hand, we get convergence to pure diffusion in the limit for large values of $\Delta t$. This means that we can interpret the multilevel Monte Carlo method in this setting as refining upon an initial simulation of the heat equation, gradually including transport effects until the correct regime set by $\epsilon$ has been achieved. We have shown that a significant speedup over classical Monte Carlo simulation is achieved when applying a geometric sequence of levels starting from $\Delta t=\epsilon^{2}$. We have also shown that adding an extra coarse level to the simulation further accelerates the computation in the considered test case.

The approach taken in developing the asymptotic-preserving scheme is general, and it is straightforward to apply the coupling described in Section 3.2 to other, more general, models. As such, we are confident that the ideas expressed in this paper will also be applicable to more general equations than the Goldstein-Taylor model studied here. In future work, this scheme can, for example, be extended to higher dimensional models, both in terms of position and velocity. More complicated models including, for example, absorption terms can also be studied. We intend to expand upon the results in Section 4.2, as well as considering varying $\epsilon$ together with $\Delta t$.

Acknowledgements.

We thank Pieterjan Robbe for many helpful discussions on the multilevel Monte Carlo method. We also thank the anonymous reviewers for their helpful suggestions for improving the quality of this work. 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.