Multilevel and Multi-index Monte Carlo methods for the McKean–Vlasov equation

We address the approximation of functionals depending on a system of particles, described by stochastic differential equations (SDEs), in the mean-field limit when the number of particles approaches infinity. This problem is equivalent to estimating the weak solution of the limiting McKean–Vlasov SDE. To that end, our approach uses systems with finite numbers of particles and a time-stepping scheme. In this case, there are two discretization parameters: the number of time steps and the number of particles. Based on these two parameters, we consider different variants of the Monte Carlo and Multilevel Monte Carlo (MLMC) methods and show that, in the best case, the optimal work complexity of MLMC, to estimate the functional in one typical setting with an error tolerance of TOL\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {TOL}$$\end{document}, is when using the partitioning estimator and the Milstein time-stepping scheme. We also consider a method that uses the recent Multi-index Monte Carlo method and show an improved work complexity in the same typical setting of . Our numerical experiments are carried out on the so-called Kuramoto model, a system of coupled oscillators.


Introduction
In our setting, a stochastic particle system is a system of coupled d-dimensional stochastic differential equations (SDEs), each modeling the state of a "particle".Such particle systems are versatile tools that can be used to model the dynamics of various complicated phenomena using relatively simple interactions, e.g., pedestrian dynamics [22,17], collective animal behavior [10,9], interactions between cells [8] and in some numerical methods such as Ensemble Kalman filters [25].One common goal of the simulation of these particle systems is to average some quantity of interest computed on all particles, e.g., the average velocity, average exit time or average number of particles in a specific region.
Under certain conditions, most importantly the exchangeability of particles and sufficient regularity of the SDE coefficients, the stochastic particle system approaches a mean-field limit as the number of particles tends to infinity [28].Exchangeability of particles refers to the assumption that all permutations of the particles have the same joint distribution.In the mean-field limit, each particle follows a single McKean-Vlasov SDE where the advection and/or diffusion coefficients depend on the distribution of the solution to the SDE [11].In many cases, the objective is to approximate the expected value of a quantity of interest (QoI) in the mean-field limit as the number of particles tend to infinity, subject to some error tolerance, TOL.While it is possible to approximate the expectation of these QoIs by estimating the solution to a nonlinear PDE using traditional numerical methods, such methods usually suffer from the curse of dimensionality.Indeed, the cost of these method is usually of O TOL − w d for some constant w > 1 that depends on the particular numerical method.Using sparse numerical methods alleviates the curse of dimensionality but requires increasing regularity as the dimensionality of the state space increases.On the other hand, Monte Carlo methods do not suffer from this curse with respect to the dimensionality of the state space.This work explores different variants and extensions of the Monte Carlo method when the underlying stochastic particle system satisfies certain crucial assumptions.We theoretically show the validity of some of these assumptions in a somewhat general setting, while verifying the other assumptions numerically on a simple stochastic particle system, leaving further theoretical justification to a future work.
Generally, the SDEs that constitute a stochastic particle system cannot be solved exactly and their solution must instead be approximated using a time-stepping scheme with a number of time steps, N .This approximation parameter and a finite number of particles, P , are the two approximation parameters that are involved in approximating a finite average of the QoI computed for all particles in the system.Then, to approximate the expectation of this average, we use a Monte Carlo method.In such a method, multiple independent and identical stochastic particle systems, approximated with the same number of time steps, N , are simulated and the average QoI is computed from each and an overall average is then taken.Using this method, a reduction of the variance of the estimator is achieved by increasing the number of simulations of the stochastic particle system or increasing the number of particles in the system.Section 3.1 presents the Monte Carlo method more precisely in the setting of stochastic particle systems.Particle methods that are not based on Monte Carlo were also discussed in [2,3].In these methods, a single simulation of the stochastic particle system is carried out and only the number of particles is increased to reduce the variance.
As an improvement of Monte Carlo methods, the Multilevel Monte Carlo (MLMC) method was first introduced in [21] for parametric integration and in [13] for SDEs; see [14] and references therein for an overview.MLMC improves the efficiency of the Monte Carlo method when only an approximation, controlled with a single discretization parameter, of the solution to the underlying system can be computed.The basic idea is to reduce the number of required samples on the finest, most accurate but most expensive discretization, by reducing the variability of this approximation with a correlated coarser and cheaper discretization as a control variate.More details are given in Section 3.2 for the case of stochastic particle systems.The application of MLMC to particle systems has been investigated in many works [4,17,27].The same concepts have also been applied to nested expectations [14].More recently, a particle method applying the MLMC methodology to stochastic particle systems was also introduced in [26] achieving, for a linear system with a diffusion coefficient that is independent of the state variable, a work complexity of O TOL −2 (log(TOL −1 )) 5 .
Recently, the Multi-index Monte Carlo (MIMC) method [19] was introduced to tackle high dimensional problems with more than one discretization parameter.MIMC is based on the same concepts as MLMC and improves the efficiency of MLMC even further but requires mixed regularity with respect to the discretization parameters.More details are given in Section 3.3 for the case of stochastic particle systems.In that section, we demonstrate the improved work complexity of MIMC compared with the work complexity of MC and MLMC, when applied to a stochastic particle system.More specifically, we show that, when using a naive simulation method for the particle system with quadratic complexity, the optimal work complexity of MIMC is O TOL −2 log TOL −1 2 when using the Milstein time-stepping scheme and O TOL −2 log TOL −1 4 when using the Euler-Maruyama time-stepping scheme.Finally, in Section 4, we provide numerical verification for the assumptions that are made throughout the current work and the derived rates of the work complexity.
In what follows, the notation a b means that there exists a constant c that is independent of a and b such that a < cb.

Problem Setting
Consider a system of P exchangeable stochastic differential equations (SDEs) where for p = 1 . . .P , we have the following equation for X p|P (t) ∈ R d dX p|P (t) = A t, X p|P (t), λ X P (t) dt + B t, X p|P (t), λ X P (t) dW p (t) where X P (t) = {X q|P (t)} P q=1 and for some (possibly stochastic) functions, where δ is the Dirac measure, is called the empirical measure.In this setting, {W p } p≥1 are mutually independent d-dimensional Wiener processes.If, moreover, {x 0 p } p≥1 are i.i.d., then under certain conditions on the smoothness and form of A and B [28], as P → ∞ for any p ∈ N, the X p|∞ stochastic process satisfies where •) being the pdf of x 0 p which is given and is independent of p. Due to (2) and x 0 p being i.i.d, {X p|∞ } p are also i.i.d.; hence, unless we want to emphasize the particular path, we drop the p-dependence in X p|∞ and refer to the random process X ∞ instead.In any case, we are interested in computing E[ψ (X ∞ (T ))] for some given function, ψ, and some final time, T < ∞.
Kuramoto Example (Fully connected Kuramoto model for synchronized oscillators).Throughout this work, we focus on a simple, one-dimensional example of (1).For p = 1, 2, . . ., P , we seek X p|P (t) ∈ R that satisfies where σ ∈ R is a constant and {ϑ p } p are i.i.d. and independent from the set of i.i.d.random variables {x 0 p } p and the Wiener processes {W p } p .The limiting SDE as P → ∞ is Note that in terms of the generic system (1) we have with ϑ a random variable and B = σ is a constant.We are interested in a real number between zero and one that measures the level of synchronization in the system with an infinite number of oscillators [1]; with zero corresponding to total disorder.In this case, we need two estimators: one where we take ψ(•) = sin(•) and the other where we take ψ(•) = cos(•).
While it is computationally efficient to approximate E[ψ(X ∞ (T ))] by solving the McKean-Vlasov PDE, that ρ ∞ satisfies, when the state dimensionality, d, is small (cf., e.g., [17]), the cost of a standard full tensor approximation increases exponentially as the dimensionality of the state space increases.On the other hand, using sparse approximation techniques to solve the PDE requires increasing regularity assumptions as the dimensionality of the state space increases.Instead, in this work, we focus on approximating the value of E[ψ(X ∞ )] by simulating the SDE system in (1).Let us now define Here, due to exchangeability, {X p|P (T )} P p=1 are identically distributed but they are not independent since they are taken from the same realization of the particle system.Nevertheless, we have E[φ P ] = E ψ(X p|P (T )) for any p and P .In this case, with respect to the number of particles, P , the cost of a naive calculation of φ P is O P 2 due to the cost of evaluating the empirical measure in (1) for every particle in the system.It is possible to take {X p|P } P p=1 in (4) as i.i.d., i.e., for each p = 1 . . .P , X p|P is taken from a different independent realization of the system (1).In this case, the usual law of large numbers applies, but the cost of a naive calculation of φ P is O P 3 .For this reason, we focus in this work on the former method of taking identically distributed but not independent {X p|P } P p=1 .Following the setup in [7,20], our objective is to build a random estimator, A, approximating with minimal work, i.e., we wish to satisfy the constraint for a given error tolerance, TOL, and a given confidence level determined by 0 < 1.We instead impose the following, more restrictive, two constraints: Statistical constraint: for a given tolerance splitting parameter, θ ∈ (0, 1), possibly a function of TOL.To show that these bounds are sufficient note that then imposing (7) gives (5).Next, we can use Markov inequality and impose Var[A] ≤ (θTOL) 2  to satisfy (7).However, by assuming (at least asymptotic) normality of the estimator, A we can get a less stringent condition on the variance as follows: Variance constraint: Here, 0 < C is such that Φ(C ) = 1 − 2 , where Φ is the cumulative distribution function of a standard normal random variable, e.g., C ≈ 1.96 for = 0.05.The asymptotic normality of the estimator is usually shown using some form of the Central Limit Theorem (CLT) or the Lindeberg-Feller theorem (see, e.g., [7,19] for CLT results for the MLMC and MIMC estimators and Figure 3-right).
As previously mentioned, we wish to approximate the values of X ∞ by using (1) with a finite number of particles, P .For a given number of particles, P , a solution to (1) is not readily available.Instead, we have to discretize the system of SDEs using, for example, the Euler-Maruyama timestepping scheme with N time steps.For n = 0, 1, 2, . . .N − 1, .
At this point, we make the following assumptions: These assumption will be verified numerically in Section 4. In general, they translate to smoothness and boundedness assumptions on A, B and ψ.Indeed, in (P1), the weak convergence of the Euler-Maruyama method with respect to the number of time steps is a standard result shown, for example, in [23] by assuming 4-time differentiability of A, B and ψ.Showing that the constant multiplying N −1 is bounded for all P is straightforward by extending the standard proof of weak convergence the Euler-Maruyama method in [23, Chapter 14] and assuming boundedness of the derivatives A, B and ψ.On the other hand, the weak convergence with respect to the number of particles, i.e., E ψ(X p|P ) → E[ψ(X ∞ )] is a consequence of the propagation of chaos which is shown, without a convergence rate, in [28] for ψ Lipschitz, B constant and A of the the form where κ(t, •, •) is Lipschitz.On the other hand, for one-dimensional systems and using the results from [24, Theorem 3.2] we can show the weak convergence rate with respect to the number of particles and the convergence rate for the variance of φ P as the following lemma shows.Below, C(R) is the space of continuous bounded functions and C k (R) is the space of continuous bounded functions whose i'th derivative is in C(R) for i = 1, . . ., k.
Lemma 2.1 (Weak and variance convergence rates w.r.t.number of particles).Consider (1) and (2) where the norms are assumed to be uniform with respect the arguments, x and y, respectively.If, moreover, Proof.The system in this lemma is a special case of the system in [24, Theorem 3.2].From there and given the assumptions of the current lemma, (10) immediately follows.Moreover, from the same reference, we can futher conclude that for 1 ≤ p = q ≤ P .Using this we can show (11) since From here, the rate of convergence for the variance of φ N P can be shown by noting that and noting that Var[φ P ] P −1 , then showing that the first term is Var φ N P − Var[φ P ] N −1 P −1 because of the weak convergence with respect to the number of time steps.
Finally, as mentioned above, with a naive method, the total cost to compute a single sample of φ N P is O N P 2 .The quadratic power of P can be reduced by using, for example, a multipole algorithm [5,16].In general, we consider the work required to compute one sample of φ N P as O (N P γp ) for a positive constant, γ p ≥ 1.

Monte Carlo methods
In this section, we study different Monte Carlo methods that can be used to estimate the previous quantity, φ ∞ .In the following, we use the notation ω (m) p: where, for each q, ω (m) q denotes the m'th sample of the set of underlying random variables that are used in calculating X N |N q|P , i.e., the Wiener path, W q , the initial condition, x 0 q , and any random variables that are used in A or B.Moreover, we sometimes write φ N P (ω (m) 1:P ) to emphasize the dependence of the m th sample of φ N P on the underlying random variables.

Monte Carlo (MC)
The first estimator that we look at is a Monte Carlo estimator.For a given number of samples, M , number of particles, P , and number of time steps, N , we can write the MC estimator as follows: Here, Hence, due to (P1), we must have P = O TOL −1 and N = O TOL −1 to satisfy (6), and, due to (P2), we must have M = O TOL −1 to satisfy (8).Based on these choices, the total work to compute Kuramoto Example.Using a naive calculation method of φ N P (i.e., γ p = 2) gives a work complexity of O TOL −4 .See also Table 1 for the work complexities for different common values of γ p .

Multilevel Monte Carlo (MLMC)
For a given L ∈ N, define two hierarchies, {N } L =0 and {P } L =0 , satisfying P −1 ≤ P and N −1 ≤ N for all .Then, we can write the MLMC estimator as follows: where we later choose the function ϕ P L due to the telescopic sum.For MLMC to have better work complexity than that of Monte Carlo, φ N P (ω ( ,m) 1:P ) and ϕ 1:P ) must be correlated for every and m, so that their difference has a smaller variance than either φ N P (ω ( ,m) 1:P ) for all > 0. Given two discretization levels, N and N −1 , with the same number of particles, P , we can generate a sample of ϕ (ω ( ,m) 1:P ).That is, we use the same samples of the initial values, {x 0 p } p≥1 , the same Wiener paths, {W p } P p=1 , and, in case they are random as in (3), the same samples of the advection and diffusion coefficients, A and B, respectively.We can improve the correlation by using an antithetic sampler as detailed in [15] or by using a higher-order scheme like the Milstein scheme [12].In the Kuramoto example, the Euler-Maruyama and the Milstein schemes are equivalent since the diffusion coefficient is constant.
On the other hand, given two different sizes of the particle system, P and P −1 , with the same discretization level, N , we can generate a sample of ϕ N P −1 (ω 1:P ) by taking In other words, we use the same P −1 sets of random variables out of the total P sets of random variables to run an independent simulation of the stochastic system with P −1 particles.We also consider another estimator that is more correlated with φ N P (ω 1:P ).The "antithetic" estimator was first independently introduced in [17, Chapter 5] and [4] and subsequently used in other works on particle systems [27] and nested simulations [14].In this work, we call this estimator a "partitioning" estimator to clearly distinguish it from the antithetic estimator in [15].We assume that P = β p P −1 for all and some positive integer β p and take That is, we split the underlying P sets of random variables into β p identically distributed and independent groups, each of size P −1 , and independently simulate β p particle systems, each of size P −1 .Finally, for each particle system, we compute the quantity of interest and take the average of the β p quantities.
In the following subsections, we look at different settings in which either P or N depends on while the other parameter is constant for all .We begin by recalling the optimal convergence rates of MLMC when applied to a generic random variable, Y , with a trivial generalization to the case when there are two discretization parameters: one that is a function of the level, , and the other, L, that is fixed for all levels., where we assume that the samples {Y ( ,m) } ,m are mutually independent.Consider the MLMC estimator with Y ,m L,−1 = 0 and for β, w, γ, s, β, w, γ, c > 0 where s ≤ 2w, assume the following: Then, for any TOL < e −1 , there exists L, L and a sequence of {M } L =0 such that and Proof.The proof can be straightforwardly derived from the proof of [6, Theorem 1], we sketch here the main steps.First, we split the constraint (15) to a bias and variance constraints similar to ( 6) to ( 8 Finally, given L, solving for {M } L =0 to minimize the work while satisfying the variance constraint gives the desired result.

MLMC hierarchy based on the number of time steps
In this setting, we take N = (β t ) for some β t > 0 and P = P L for all , i.e., the number of particles is a constant, P L , on all levels.We make an extra assumption in this case, namely: for some constant s t > 0. The factor (β t ) −st is the usual assumption on the variance convergence of the level difference in MLMC theory [13] and is a standard result for the Euler-Maruyama scheme with s t = 1 and for the Milstein scheme with s t = 2, [23].On the other hand, the factor P −1 L can be motivated from (P2), which states that the variance of each term in the difference converges at this rate.
Due to Theorem 3.1, we can conclude that the work complexity of MLMC is Kuramoto Example.In this example, using the Milstein time-stepping scheme, we have s t = 2 (cf. Figure 1), and a naive calculation method of φ N P (γ p = 2) gives a work complexity of O TOL −3 .See also Table 1 for the work complexities for different common values of s t and γ p .Kuramoto Example.We choose β p = β t and use a naive calculation method of φ N P (yielding γ p = 2) and the partitioning sampler (yielding s p = 1).Finally, using the Milstein time-stepping scheme, we have s t = 2. Refer to Figure 1 for numerical verification.Based on these rates, we have, in (19), s = 2 log(β p ), w = log(β p ) and γ = 3 log(β p ).The MLMC work complexity in this case is O TOL −3 .See also Table 1 for the work complexities for different common values of s t and γ p .

Multi-index Monte Carlo (MIMC)
Following [19], for every multi-index α = (α 1 , α 2 ) ∈ N 2 , let P α1 = (β p ) α1 and N α2 = (β t ) α2 and define the first-order mixed-difference operator in two dimensions as The MIMC estimator is then written for a given I ⊂ N 2 as At this point, similar to the original work on MIMC [19], we make the following assumptions on the convergence of ∆φ Assumption (MIMC1) is motivated from (P1) by assuming that the mixed first order difference, ∆φ Pα 1 , gives a product of the convergence terms instead of a sum.Similarly, (MIMC2) is motivated from (MLMC1) and (MLMC2).To the best of our knowledge, there are currently no proofs of these assumptions for particle systems, but we verify them numerically for (3) in Figure 2.
Henceforth, we will assume that β t = β p for easier presentation.Following [19, Lemma 2.1] and recalling the assumption on cost per sample, Work ∆φ Nα 2 Pα 1 P γp α1 N α2 , then, for every value of L ∈ R + , the optimal set can be written as and the optimal computational complexity of MIMC is O TOL −2−2 max(0,ζ) log TOL −1 p , where Kuramoto Example.Here again, we use a naive calculation method of φ N P (yielding γ p = 2) and the partitioning sampler (yielding s p = 1).Finally, using the Milstein time-stepping scheme, we have s t = 2. Hence, ζ = 0, z = 1 and Work [A MIMC ] = O TOL −2 log TOL −1 2 .See also

Numerical Example
In this section we provide numerical evidence of the assumptions and work complexities that were made in the Section 3.This section also verifies that the constants of the work complexity (which were not tracked) are not significant for reasonable error tolerances.The results in this section were obtained using the mimclib software library [18] and GNU parallel [29].
Figure 1 shows the absolute expectation and variance of the level differences for the different MLMC settings that were outlined in Section 3.2.These figures verify Assumptions (P1), (P2) and (MLMC1)-(MLMC3) with the values s t = 2 and s p = 0 for the ϕ sampler in (13) or the value s p = 1 for the ϕ sampler in (14).For the same parameter values, Figure 2 provides numerical evidence for Assumptions (MIMC1) and (MIMC2) for the ϕ sampler (14).
We now compare the MLMC method [13] in the setting that was presented in Section 3.2.3 and the MIMC method [19] that was presented in Section 3.3.In both methods, we use the Milstein time-stepping scheme and the partitioning sampler, ϕ, in (14).Recall that in this case, we verified numerically that γ p = 2, s p = 1 and s t = 2.We also use the MLMC and MIMC algorithms that were outlined in their original work and use an initial 25 samples on each level or multi-index to compute a corresponding variance estimate that is required to compute the optimal number of samples.In the following, we refer to these methods as simply "MLMC" and "MIMC".We focus on the settings in Sections 3.2.3 and 3.3 since checking the bias of the estimator in those settings can be done straightforwardly by checking the absolute value of the level differences in MLMC or the multi-index differences in MIMC.On the other hand, checking the bias in the settings outlined in Sections 3.1, 3.2.1 and 3.2.2 is not as straightforward and determining the number of times steps and/or the number of particles to satisfy a certain error tolerance requires more sophisticated algorithms.This makes a fair numerical comparison with these later settings somewhat difficult.
Figure 3-left shows the exact errors of both MLMC and MIMC for different prescribed tolerances.This plot shows that both methods estimate the quantity of interest up to the same error tolerance; comparing their work complexity is thus fair.On the other hand, Figure 3-right is a PP plot, i.e., a plot of the cumulative distribution function (CDF) of the MLMC and MIMC estimators, normalized by their variance and shifted by their mean, versus the CDF of a standard normal distribution.This figure shows that our assumption in Section 2 of the asymptotic normality of these estimators is well founded.Figure 4 shows the maximum discretization level for both the number of time steps and the number of particles for MLMC and MIMC (cf.( 22)).Recall that, for a fixed tolerance in MIMC, 2α 2 + α 1 is bounded by a constant (cf.( 21)).Hence, Figure 4 has a direct implication on the results reported in Figure 5 where we plot the maximum cost of the samples used in both MLMC and MIMC for different tolerances.This cost represents an indivisible unit of simulation for both methods, assuming we treat the simulation of the particle system as a black box.Hence, Figure 5 shows that MIMC has better parallelization scaling, i.e., even with an infinite number of computation nodes MIMC would still be more efficient than MLMC. to (22).From the right plot, we can confirm that s t = 2 for the Milstein method.We can also deduce that using the ϕ sampler in (13) yields s p = 0 in (MIMC2) (i.e., no variance reduction compared to Var φ N P ) while using the ϕ sampler in ( 14) yields s p = 1 in (MIMC2) (i.e.O P −2 ).
Finally, we show in Figure 6 the cost estimates of MLMC and MIMC for different tolerances.This figure clearly shows the performance improvement of MIMC over MLMC and shows that the complexity rates that we derived in this work are reasonably accurate.

Conclusions
This work has shown both numerically and theoretically under certain assumptions, that could be verified numerically, the improvement of MIMC over MLMC when used to approximate a quantity of interest computed on a particle system as the number of particles goes to infinity.The application to other particle systems (or equivalently other McKean-Vlasov SDEs) is straightforward and similar improvements are expected.The same machinery was also suggested for approximating nested expectations in [14] and the analysis here applies to that setting as well.Moreover, the same machinery, i.e., multi-index structure with respect to time steps and number of particles coupled with a partitioning estimator, could be used to create control variates to reduce the computational cost of approximating quantities of interest on stochastic particle systems with a finite number of particles.
Future work includes analyzing the optimal level separation parameters, β p and β t , and the behavior of the tolerance splitting parameter, θ.Another direction could be applying the MIMC method to higher-dimensional particle systems such as the crowd model in [17].On the theoretical side, the next step is to prove the assumptions that were postulated and verified numerically in this work for certain classes of particle systems, namely: the second order convergence with respect to the number of particles of the variance of the partitioning estimator ( 14) and the convergence rates for mixed differences (MIMC1) and (MIMC2).

E
[A MC (M, P, N )] = E φ work is Work [A MC (M, P, N )] = M N P γp .

Theorem 3 . 1 (
Optimal MLMC complexity).Let Y L, be an approximation of the random variable, Y , for every ( L, ) ∈ N 2 .Denote by Y ( ,m) a sample of Y and denote its corresponding approximation by Y ( ,m) L,

1 w
), respectively.Then, since E A MLMC ( L, L) = E Y L,L , given the first assumption of the theorem and imposing the bias constraint yield L = O 1 w log( β) log(TOL −1 ) and L = O log(β) log(TOL −1 ) .The assumptions on the variance and work then give:

Table 1 :
The work complexity of the different methods presented in this work in common situations, encoded as (a, b) to represent O TOL −a (log(TOL −1 )) b .When appropriate, we use the partitioning estimator (i.e., s p = 1).In general, MIMC has always the best complexity.However, when γ p = 1 MIMC does not offer an advantage over an appropriate MLMC method.