Multilevel and Multi-index Monte Carlo methods for the McKean-Vlasov - - PowerPoint PPT Presentation

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

Abdul-Lateef Haji-Ali∗ Ra´ ul Tempone†

∗Heriot-Watt University, United Kingdom † Alexander von Humboldt Professor, RWTH Aachen, Germany,

KAUST, Saudi Arabia. DCSE Fall School on ROM and UQ

November, 2019

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Motivation

Particle Systems in the Mean-field

Particle systems are a collection of coupled, usually identical and simple, models that can be used to model complicated phenomena.

Molecular dynamics, Crowd simulation, Oscillators, . . .

Certain stochastic particles systems have a mean-field limit when the number of particles increase. Such limits can be useful to understand their complicated phenomena.

1/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Motivation

Main reference

This presentation is based on the manuscript

• ”Multilevel and Multi-index Monte Carlo methods for the

McKean-Vlasov equation” by A. L. Haji Ali and R. T. Statistics and Computing, Vol. 28(4), pp. 923–935, 2018. A McKean-Vlasov process is a stochastic process described by a SDE whose coefficients depend on the distribution of the solution

• itself. They relate to the Vlasov model for plasma evolution and

were first studied by Henry McKean in 1966. For 0 < t ≤ T the process X(t) solves dX(t) =a(X(t), µ(t))dt + b(X(t), µ(t))dW (t), µ(t) =L(X(t)) and µ(0) given. Goal: approximate E[g(X(T))] for some given g.

2/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Motivation

Convergence to the mean-field

For particles Xp|P, p = 1, . . . , P, (evolving in a system of size P) define ”shadow” particles Xp|∞, (evolving in a system of ∞ size) Xp|P(t) = x0

p +

t  a(Xp|P(t)) + 1 P

P

• q=1

A(Xp|P(t), Xq|P(t))   dt + σWp(t) Xp|∞(t) = x0

p +

t

• a(Xp|∞(t)) +
• A(Xp|∞(t), y)µ∞(t)(dy)
• dt

+ σWp(t), with µ∞(t) the marginal distribution for Xp|∞(t) for any p. Consistency: The initial values x0

p, p = 1, . . . , P, are i.i.d. from

µ∞(0).

3/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Motivation

For t > 0, and all x, the pdf of the marginal distribution µ∞(t) of the infinite size system satisfies a nonlinear Fokker Planck equation ∂tρ∞(t, x)+div (ρ∞(t, x)(a(x) + ρ∞(t, ·) ∗ A(x, ·))) =

• i

σ2

i

2 ∂2

i ρ∞(t, x)

with a given initial condition ρ∞(0, ·) and suitable b.c. Question: What about the rate of weak convergence? E

• g(Xp|P(T)) − g(Xp|∞(T))
• O(1/P)

4/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Simple Example

Kuramoto oscillator model †

For p = 1, 2, . . . , P consider equally coupled oscillators with intrinsic natural frequencies ϑp that follow a system of Itˆ

• SDEs

dXp|P(t) =

• ϑp + 1

P

P

• q=1

sin(Xp|P(t) − Xq|P(t))

• dt + σdWp(t)

Xp|P(0) = x0

p

where we are interested in Total order =

• 1

P

P

• p=1

cos

• Xp|P(T)
• 2

+

• 1

P

P

• p=1

sin

• Xp|P(T)
• 2

; a real number between zero and one that measures the level of synchronization of the coupled oscillators.

†Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence, Springer, Berlin, 1984. 5/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Simple Example

Kuramoto oscillator model †

For p = 1, 2, . . . , P consider equally coupled oscillators with intrinsic natural frequencies ϑp that follow a system of Itˆ

• SDEs

dXp|P(t) =

• ϑp + 1

P

P

• q=1

sin(Xp|P(t) − Xq|P(t))

• dt + σdWp(t)

Xp|P(0) = x0

p

where we are interested in: φP = 1 P

P

• p=1

cos

• Xp|P(T)
• ,

Mean-field limit: φP → φ∞ = E

• cos(Xp|∞(T))
• as

P ↑ ∞ dXp|∞ =

• ϑp +
• R

sin(Xp|∞(t) − y)µ∞(t, dy)

• dt + σdWp(t)

†Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence, Springer, Berlin, 1984. 5/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Simple Example

Kuramoto oscillator model †, Euler-Maruyama

For p = 1, 2, . . . , P consider equally coupled oscillators with intrinsic natural frequencies ϑp that follow a system of Itˆ

• SDEs

X n|N

p|P − X n−1|N p|P

=

• ϑp + 1

P

P

• q=1

sin(X n|N

p|P − X n|N q|P )

• T

N + σ∆W n|N

p

X 0|N

p|P = x0 p

where we are interested in: φN

P = 1

P

P

• p=1

cos

• X N|N

p|P

• ,

Mean-field limit: φP → φ∞ = E

• cos(Xp|∞(T))
• as

P ↑ ∞ dXp|∞ =

• ϑp +
• R

sin(Xp|∞(t) − y)µ∞(t, dy)

• dt + σdWp(t)

†Y. Kuramoto, Chemical Oscillations, Waves, and Turbulence, Springer, Berlin, 1984. 5/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Simple Example

Objective

Our objective is to build an estimator A ≈ φ∞ with minimal work where P(|A − φ∞| ≤ TOL) ≥ 1 − ǫ for a given accuracy TOL and a given confidence level determined by 0 < ǫ ≪ 1.

6/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Simple Example

Objective

Our objective is to build an estimator A ≈ φ∞ with minimal work where P(|A − φ∞| ≤ TOL) ≥ 1 − ǫ for a given accuracy TOL and a given confidence level determined by 0 < ǫ ≪ 1. We instead impose the following, more restrictive, two constraints: Bias constraint: |E[A] − φ∞| ≤ TOL/3, Statistical constraint: P (|A − E[A]| ≥ 2TOL/3)

6/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Simple Example

Objective

Our objective is to build an estimator A ≈ φ∞ with minimal work where P(|A − φ∞| ≤ TOL) ≥ 1 − ǫ for a given accuracy TOL and a given confidence level determined by 0 < ǫ ≪ 1. We instead impose the following, more restrictive, two constraints: Bias constraint: |E[A] − φ∞| ≤ TOL/3, Variance constraint: Var[A] ≤ (2TOL/3Cǫ)2 . assuming (at least asymptotic a) normality of the estimator, A. Here, 0 < Cǫ is such that Φ(Cǫ) = 1 − ǫ

2, where Φ is the c.d.f. of a

standard normal random variable.

• aN. Collier, A.-L. Haji-Ali, E. von Schwerin, F. Nobile, and R. Tempone. “A continuation

multilevel Monte Carlo algorithm”. BIT Numerical Mathematics, 55(2):399-432, (2015).

6/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Monte Carlo (MC)

Monte Carlo

The simplest estimator is the Monte Carlo estimator AMC = 1 M

M

• m=1

φN

P (ωm P ).

For a given P, N and M that we can choose to satisfy the error constraints and minimize the work. Here ωm

P =

• ωm

p

P

p=1 and for

each particle, ωm

p denotes the independent, identically distributed

(i.i.d.) samples of the set of underlying random variables that are used in calculating X N|N

p|P , 1 ≤ p ≤ P.

7/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Monte Carlo (MC)

Monte Carlo work complexity

In our 1D example, we can check (at least numerically) that Minimize total work: Work(AMC), such that: Bias(AMC) =

• φ∞ − E
• φN

P

• ≤ TOL

3 and: Var[AMC] = Var

• φN

P

• M

≤ 2TOL 3Cǫ 2

8/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Monte Carlo (MC)

Monte Carlo work complexity

In our 1D example, we can check (at least numerically) that Minimize total work: Work(AMC) = O

• MNP2

such that: Bias(AMC) = O

• N−1

+ O

• P−1

≤ TOL 3 and: Var[AMC] = O

• P−1

M ≤ 2TOL 3Cǫ 2

8/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Monte Carlo (MC)

Monte Carlo work complexity

In our 1D example, we can check (at least numerically) that Minimize total work: Work(AMC) = O

• MNP2

such that: Bias(AMC) = O

• N−1

+ O

• P−1

≤ TOL 3 and: Var[AMC] = O

• P−1

M ≤ 2TOL 3Cǫ 2 In this case, we choose P = O

• TOL−1

, N = O

• TOL−1

, M = O

• TOL−1

and the total cost of a naive Monte Carlo is O

• TOL−4

. Observe: The cost of a “single cloud” naive method with M = 1 is O

• TOL−5

8/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – Introduction

Following (Heinrich, 2001) and (Giles, 2008), For a given L ∈ N, define two hierarchies {Nℓ}L

ℓ=1 and {Pℓ}L ℓ=1 satisfying Pℓ−1 ≤ Pℓ

and Nℓ−1 ≤ Nℓ for all ℓ. Recall the telescopic decomposition φ∞ ≈ E

• φNL

PL

• = E
• φN0

P0

• +

L

• ℓ=1

E

• φNℓ

Pℓ − ϕNℓ−1 Pℓ−1

• 9/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – Introduction

Following (Heinrich, 2001) and (Giles, 2008), For a given L ∈ N, define two hierarchies {Nℓ}L

ℓ=1 and {Pℓ}L ℓ=1 satisfying Pℓ−1 ≤ Pℓ

and Nℓ−1 ≤ Nℓ for all ℓ. Recall the telescopic decomposition φ∞ ≈ E

• φNL

PL

• = E
• φN0

P0

• +

L

• ℓ=1

E

• φNℓ

Pℓ − ϕNℓ−1 Pℓ−1

• =

L

• ℓ=0

E[∆ℓφ]. where ∆ℓφ =

• φN0

P0

if ℓ = 0, φNℓ

Pℓ − ϕNℓ−1 Pℓ−1

if ℓ > 0. Here, we assume that the auxiliary estimator ϕ satisfies E

• ϕNℓ−1

Pℓ−1

• = E
• φNℓ−1

Pℓ−1

• 9/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – Introduction

Then, using Monte Carlo to approximate each level independently, the MLMC estimator can be written as AMLMC =

L

• ℓ=0

1 Mℓ

Mℓ

• m=1

∆ℓφ(ωℓ,m

Pℓ ).

where Mℓ is optimally chosen. High correlation is crucial (between the pairs (Nℓ, Nℓ−1)? (Pℓ, Pℓ−1)?) to ensure that Var[∆ℓφ] goes to zero sufficiently fast.

10/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – Introduction

Recall: MLMC optimal work complexity ‡ §

Bias: |E[∆ℓφ]| = O (exp(−wℓ)), Variance: Var[∆ℓφ] = O (exp(−sℓ)), Work: Work[∆ℓφ] = O (exp(γℓ)).

‡Cliffe, K.A. and Giles, M.B. and Scheichl, R. and Teckentrup, A. Computing and

Visualization in Science, “Multilevel Monte Carlo methods and applications to elliptic PDEs with random coefficients” (2011).

§Giles, Acta Numerica 2015. 11/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – Introduction

Recall: MLMC optimal work complexity ‡ §

Bias: |E[∆ℓφ]| = O (exp(−wℓ)), Variance: Var[∆ℓφ] = O (exp(−sℓ)), Work: Work[∆ℓφ] = O (exp(γℓ)). The optimal work of MLMC is        O

• TOL−2

s > γ O

• TOL−2

log

• TOL−12

s = γ O

• TOL−2− γ−s

w

• s < γ

Recall the total cost of Monte Carlo is O

• TOL−2− γ

w

• ‡Cliffe, K.A. and Giles, M.B. and Scheichl, R. and Teckentrup, A. Computing and

Visualization in Science, “Multilevel Monte Carlo methods and applications to elliptic PDEs with random coefficients” (2011).

§Giles, Acta Numerica 2015. 11/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – Introduction

In the following, we look at different settings in which either Pℓ or Nℓ depends on ℓ while the other parameter is constant for all ℓ.

12/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles

MLMC in number of particles, P

Pℓ = 2ℓ and Nℓ = NL. Build correlated samples by sampling 2ℓ and sub-sampling 2ℓ−1 particles out of them (e.g. the first 2ℓ−1). Use the same initial conditions, Brownian paths or any other random variables associated to a particle. ϕNL

Pℓ−1(ωℓ,m Pℓ ) = φNL Pℓ−1

• ωℓ,m

p

Pℓ−1

p=1

• 14/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles

MLMC in number of particles, P

Pℓ = 2ℓ and Nℓ = NL.

1 2 3 4 5 6 7 8 9 10 ℓ 10−6 10−5 10−4 10−3 10−2 10−1 100 101

Expectation

1 2 3 4 5 6 7 8 9 10 ℓ 10−10 10−9 10−8 10−7 10−6 10−5 10−4 10−3 10−2

Variance

φNℓ

P0

φN0

Pℓ

φNℓ

P0 − φNℓ−1 P0

φN0

Pℓ − ϕN0 Pℓ−1

O

• 2−ℓ

O

• 2−2ℓ

s = w = 1

14/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles

MLMC in number of particles, P

Pℓ = 2ℓ and Nℓ = NL.

2 4 6 8 10 12 14 16 ℓ 10−6 10−5 10−4 10−3 10−2 10−1 100 Seconds φNℓ

P0 − φNℓ−1 P0

φN0

Pℓ − ϕN0 Pℓ−1

O

• 2ℓ

O

• 22ℓ

O

• 23ℓ

Work, γ = 2

14/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles

MLMC in number of particles, P

Pℓ = 2ℓ and Nℓ = NL. Summary: w = 1, s = 1, γ = 2 Fixing NL, the optimal work of biased MLMC is O

• TOL−3

.

14/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles

MLMC in number of particles, P

Pℓ = 2ℓ and Nℓ = NL. Summary: w = 1, s = 1, γ = 2 Fixing NL, the optimal work of biased MLMC is O

• TOL−3

. To control the bias O

• N−1

L

• , choose NL = O
• TOL−1

Cost per sample: O (NLPℓ) Variance: O

• P−1

ℓ

• then the total cost becomes

O

• TOL−4

. No advantage of MLMC over MC? Need a more correlated estimator!

14/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles

Summary

Method Work complexity Monte Carlo O

• TOL−4

MLMC in N O

• TOL−3

MLMC in P O

• TOL−4

15/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles

Reducing the variance w.r.t P; partitioning (antithetic) sampler

The crucial element is how fast Var[∆φℓ] goes to zero compared to how much it costs to compute ∆φℓ. Better choice of ϕ(·) to reduce Var[∆ℓφ].

16/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles

Reducing the variance w.r.t P; partitioning (antithetic) sampler

The crucial element is how fast Var[∆φℓ] goes to zero compared to how much it costs to compute ∆φℓ. Better choice of ϕ(·) to reduce Var[∆ℓφ]. If Pℓ−1 = Pℓ/2, then given a group of Pℓ particles, we can subsample two identically-distributed, independent groups that have Pℓ−1 particles and average the QoI

ϕNL

Pℓ−1(ωℓ,m Pℓ ) = 1

2

• φNL

Pℓ−1

• ωℓ,m

p

Pℓ−1

p=1

• + φNL

Pℓ−1

• ωℓ,m

p

Pℓ

p=1+Pℓ−1

• 16/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles, partitioning

MLMC in number of particles P, with partitioning sampler

Pℓ = 2ℓ and Nℓ = NL. Build correlated samples by sampling 2ℓ and sub-sampling two identically-distributed, independent groups of 2ℓ−1 particles

• ut of them. Use the same initial conditions, Brownian paths
• r any other random variables associated to a particle.

ϕNL

Pℓ−1(ωℓ,m Pℓ ) = 1

2

• φNL

Pℓ−1

• ωℓ,m

p

Pℓ−1

p=1

• + φNL

Pℓ−1

• ωℓ,m

p

Pℓ

p=1+Pℓ−1

• 17/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles, partitioning

MLMC in number of particles P, with partitioning sampler

Pℓ = 2ℓ and Nℓ = NL.

1 2 3 4 5 6 7 8 9 10 ℓ 10−6 10−5 10−4 10−3 10−2 10−1 100 101

Expectation

1 2 3 4 5 6 7 8 9 10 ℓ 10−10 10−9 10−8 10−7 10−6 10−5 10−4 10−3 10−2

Variance

φNℓ

P0

φN0

Pℓ

φNℓ

P0 − φNℓ−1 P0

φN0

Pℓ − ϕN0 Pℓ−1

φN0

Pℓ −

ϕN0

Pℓ−1

O

• 2−ℓ

O

• 2−2ℓ

Bias, w = 1, s = 2

17/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles, partitioning

MLMC in number of particles P, with partitioning sampler

Pℓ = 2ℓ and Nℓ = NL.

2 4 6 8 10 12 14 16 ℓ 10−6 10−5 10−4 10−3 10−2 10−1 100 Seconds φNℓ

P0 − φNℓ−1 P0

φN0

Pℓ − ϕN0 Pℓ−1

O

• 2ℓ

O

• 22ℓ

O

• 23ℓ

Work, γ = 2

17/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles, partitioning

MLMC in number of particles P, with partitioning sampler

Pℓ = 2ℓ and Nℓ = NL. Summary: w = 1, s = 2, γ = 2 Fixing NL, the optimal work of biased MLMC is O

• TOL−2 (log TOL)2

.

17/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles, partitioning

MLMC in number of particles P, with partitioning sampler

Pℓ = 2ℓ and Nℓ = NL. Summary: w = 1, s = 2, γ = 2 Fixing NL, the optimal work of biased MLMC is O

• TOL−2 (log TOL)2

. To control bias O

• N−1

L

• , choose NL = O
• TOL−1

Cost per sample: O (NLPℓ) Variance: O

• P−1

ℓ

• then the total cost becomes

O

• TOL−3 (log TOL)2

. Can we do better?

17/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles, partitioning

Summary

Method Work complexity Monte Carlo O

• TOL−4

MLMC in N O

• TOL−3

MLMC in P O

• TOL−4

MLMC in P, partitioning O

• TOL−3 log(TOL−1)2

18/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles and time-steps

MLMC in P and N, with partitioning samplers

Pℓ = 2ℓ and Nℓ = 2ℓ. Build correlated samples by

Sampling 2ℓ and sub-sampling two identically-distributed, independent groups of 2ℓ−1 particles out of them. Use the same initial conditions, Brownian paths or any other random variables associated to a particle At the same time, by using the same Brownian paths discretized with different meshes 2ℓ and 2ℓ−1. ϕNℓ−1

Pℓ−1(ωℓ,m Pℓ ) = 1

2

• φNℓ−1

Pℓ−1

• ωℓ,m

p

Pℓ−1

p=1

• + φNℓ−1

Pℓ−1

• ωℓ,m

p

Pℓ

p=1+Pℓ−1

• 19/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles and time-steps

MLMC in P and N, with partitioning samplers

Pℓ = 2ℓ and Nℓ = 2ℓ.

1 2 3 4 5 6 7 8 9 10 ℓ 10−6 10−5 10−4 10−3 10−2 10−1 100 101

Expectation

1 2 3 4 5 6 7 8 9 10 ℓ 10−10 10−9 10−8 10−7 10−6 10−5 10−4 10−3 10−2

Variance

φNℓ

P0

φN0

Pℓ

φNℓ

P0 − φNℓ−1 P0

φN0

Pℓ − ϕN0 Pℓ−1

φN0

Pℓ −

ϕN0

Pℓ−1

φNℓ

Pℓ −

ϕNℓ−1

Pℓ−1

O

• 2−ℓ

O

• 2−2ℓ

w = 1, s = 2

19/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles and time-steps

MLMC in P and N, with partitioning samplers

Pℓ = 2ℓ and Nℓ = 2ℓ.

2 4 6 8 10 12 14 16 ℓ 10−6 10−5 10−4 10−3 10−2 10−1 100 Seconds φNℓ

P0 − φNℓ−1 P0

φN0

Pℓ − ϕN0 Pℓ−1

φNℓ

Pℓ − ϕNℓ−1 Pℓ−1

O

• 2ℓ

O

• 22ℓ

O

• 23ℓ

Work, γ = 3

19/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles and time-steps

MLMC in P and N, with partitioning samplers

Pℓ = 2ℓ and Nℓ = 2ℓ. Summary: w = 1, s = 2, γ = 3 The optimal work of asymptotically unbiased MLMC is O

• TOL−3

= O

• TOL−(2+ γ−s

w )

. The number of particles is not helping. Can we do better?

19/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multilevel Monte Carlo (MLMC) – In number of particles and time-steps

Summary

Method Work complexity Monte Carlo O

• TOL−4

MLMC in N O

• TOL−3

MLMC in P O

• TOL−4

MLMC in P, partitioning O

• TOL−3 log(TOL−1)2

MLMC in P and N O

• TOL−4

MLMC in P and N, partitioning O

• TOL−3

20/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multi-Index Monte Carlo

Variance reduction: MLMC

P N

21/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multi-Index Monte Carlo

Variance reduction: Further potential

P N

21/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multi-Index Monte Carlo – General Framework

MIMC Estimator¶

For α = (α1, α2) ∈ N2, let Pα1 = 2α1 and Nα2 = 2α2. Define the first order mixed difference ∆αφ = ∆1,α(∆2,αφ) = (φ

Nα2 Pα1 − φ Nα2−1 Pα1

) − (ϕ

Nα2 Pα1−1 − ϕ Nα2−1 Pα1−1).

¶A.-L. Haji-Ali, F. Nobile, and R. Tempone. “Multi-Index Monte Carlo: When Sparsity Meets

Sampling”. Numerische Mathematik, 1-40, (2015)

22/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multi-Index Monte Carlo – General Framework

MIMC Estimator¶

For α = (α1, α2) ∈ N2, let Pα1 = 2α1 and Nα2 = 2α2. Define the first order mixed difference ∆αφ = ∆1,α(∆2,αφ) = (φ

Nα2 Pα1 − φ Nα2−1 Pα1

) − (ϕ

Nα2 Pα1−1 − ϕ Nα2−1 Pα1−1).

Then the MIMC estimator can be written as AMIMC =

• α∈I

1 Mα

Mα

• m=1

∆αφ(ωα,m) for some properly chosen index set I ⊂ N2 and optimal number of samples Mα for every α ∈ I.

¶A.-L. Haji-Ali, F. Nobile, and R. Tempone. “Multi-Index Monte Carlo: When Sparsity Meets

Sampling”. Numerische Mathematik, 1-40, (2015)

22/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multi-Index Monte Carlo – General Framework

MIMC optimal work complexity

Bias: |E[∆αφ]| = O (exp(−w1α1 − w2α2)), Variance: Var[∆αφ] = O (exp(−s1α1 − s2α2)), Work: Work[∆αφ] = O (exp(γ1α1 + γ2α2)).

23/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multi-Index Monte Carlo – General Framework

MIMC optimal work complexity

Bias: |E[∆αφ]| = O (exp(−w1α1 − w2α2)), Variance: Var[∆αφ] = O (exp(−s1α1 − s2α2)), Work: Work[∆αφ] = O (exp(γ1α1 + γ2α2)). The optimal set I I(L) =

• (α1, α2) ∈ N2 :

(2w1 + γ1 − s1) α1 + (2w2 + γ2 − s2)α2 ≤ L

• 23/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multi-Index Monte Carlo – General Framework

MIMC optimal work complexity

Bias: |E[∆αφ]| = O (exp(−w1α1 − w2α2)), Variance: Var[∆αφ] = O (exp(−s1α1 − s2α2)), Work: Work[∆αφ] = O (exp(γ1α1 + γ2α2)). The optimal set I I(L) =

• (α1, α2) ∈ N2 :

(2w1 + γ1 − s1) α1 + (2w2 + γ2 − s2)α2 ≤ L

• Letting ζ = max
• γ1−s1

2w1 , γ2−s2 2w2

• , then the optimal work of MIMC is

O

• TOL−2(1+max(0,ζ)) log
• TOL−1p

for p ≥ 0.

23/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multi-Index Monte Carlo – General Framework

Optimal choice of I: Total degree set

L α1 α2 Iδ(L) =

• α ∈ Nd : δ · α ≤ L
• 24/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multi-Index Monte Carlo – Application

MIMC, with partitioning samplers

Let Pα1 = 2α1 and Nα2 = 2α2. Build correlated samples by

Sampling 2α1 and sub-sampling two identically-distributed, independent groups of 2α1−1 particles out of them. At the same time, by using the same Brownian paths discretized with different meshes 2α2 and 2α2−1. Use MIMC levels: Mixed differences!

25/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multi-Index Monte Carlo – Application

MIMC, with partitioning samplers

Let Pα1 = 2α1 and Nα2 = 2α2.

3 6 9 12 15 18 21 24 i 10−10 10−9 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 100

Expectation

3 6 9 12 15 18 21 24 i 10−19 10−17 10−15 10−13 10−11 10−9 10−7 10−5 10−3

Variance

α = (i, 0) α = (0, i) α = (i, i) 2−i 2−2i 2−4i

w1, w2 = 1, s1 = s2 = 2 Notice higher rates for mixed difference.

25/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multi-Index Monte Carlo – Application

MIMC, with partitioning samplers

Let Pα1 = 2α1 and Nα2 = 2α2.

5 10 15 20 25 i 10−6 10−5 10−4 10−3 10−2 10−1 100 101 102

Time

α = (i, 0) α = (0, i) α = (i, i) 2i 22i 23i

γ1 = 2, γ2 = 1

25/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multi-Index Monte Carlo – Application

MIMC, with partitioning samplers

Let Pα1 = 2α1 and Nα2 = 2α2. Summary: w1 = w2 = 1 s1 = s2 = 2 γ1 = 2γ2 = 2    = ⇒ ζ = max γ1 − s1 2w1 , γ2 − s2 2w2

• = 0

The optimal set I(L) =

• (α1, α2) ∈ N2 : 2α1 + 3α2 ≤ L
• The optimal work of the asymptotically unbiased MIMC is

O

• TOL−2 log
• TOL−12

25/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multi-Index Monte Carlo – Application

Summary

Method Work complexity Monte Carlo O

• TOL−4

MLMC in N O

• TOL−3

MLMC in P O

• TOL−4

MLMC in P, partitioning O

• TOL−3 log(TOL−1)2

MLMC in P and N O

• TOL−4

MLMC in P and N, partitioning O

• TOL−3

MIMC O

• TOL−2 log
• TOL−12

26/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multi-Index Monte Carlo – Results

Numerical Example: MIMC vs. MLMC

X n|N

p|P − X n−1|N p|P

=  ϑp + 0.4 P

P

• q=1

sin(X n|N

p|P − X n|N q|P )

  T N + 0.4∆W n|N

p|P

X 0|N

p|P ∼ N(0, 0.2)

where ϑp ∼ U(−0.2, 0.2). The quantity of interest is φN

P = 1

P

P

• p=1

cos

• X N|N

p|P

• .

for T = 1.

27/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multi-Index Monte Carlo – Results

Numerical Example: MIMC vs. MLMC

for T = 1.

10−6 10−5 10−4 10−3 10−2 10−1 TOL 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 Error TOL MIMC MLMC

27/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multi-Index Monte Carlo – Results

Numerical Example: MIMC vs. MLMC

for T = 1.

10−5 10−4 10−3 10−2 10−1 TOL 102 103 104 105 106 107 108 109 1010 1011 Work Estimate TOL−2 log(TOL−1)2 TOL−3 MIMC MLMC

27/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Multi-Index Monte Carlo – Results

Numerical Example: MIMC vs. MLMC

for T = 1.

10−5 10−4 10−3 10−2 10−1 TOL 10−2 10−1 100 101 102 103 104 105 Running time (sec) TOL−2 log(TOL−1)2 TOL−3 MIMC MLMC

27/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Conclusions

Conclusions

MLMC applied to particle systems in the mean-field limit is not trivial since the variance of the quantity of interest depends on the number of particles, P. The partitioning (antithetic) estimator achieves much better L2 convergence rate w.r.t. P. Considerable computational saving by using MIMC, from TOL−4 to TOL−2 log2(TOL). Other applications: higher dimension particle systems (e.g. crowd flow). MLMC and MIMC Theory: working on estimates on mixed differences.

28/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Conclusions

References

Sznitman, A.S., 1991. Topics in propagation of chaos. In Ecole d’t de probabilits de Saint-Flour XIX1989 (pp. 165-251). Springer, Berlin, Heidelberg. Talay, D. and Tubaro, L., 1990. Expansion of the global error for numerical schemes solving stochastic differential equations. Stochastic analysis and applications, 8(4), pp.483-509. Bally, V. and Talay, D., 1996. The law of the Euler scheme for stochastic differential equations. Probability theory and related fields, 104(1), pp.43-60. Szepessy, A., Tempone, R. and Zouraris, G.E., 2001. Adaptive weak approximation of stochastic differential equations. Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 54(10), pp.1169-1214. Katsoulakis, M.A. and Szepessy, A., 2006. Stochastic hydrodynamical limits of particle systems. Communications in Mathematical Sciences, 4(3), pp.513-549. Haji-Ali, A.L., Nobile, F. and Tempone, R., 2016. Multi-index Monte Carlo: when sparsity meets sampling. Numerische Mathematik, 132(4), pp.767-806. Haji-Ali, A.L. and Tempone, R., 2018. Multilevel and Multi-index Monte Carlo methods for the McKeanVlasov equation. Statistics and Computing, 28(4), pp.923-935.

29/0

Monte Carlo Methods for McKean-Vlasov [Haji-Ali, Tempone] Appendix – Multi-pole algorithms

A multi-pole algorithm can be used to reduce the computation cost with respect to the number of particles, P. We can possibly achieve γP ≈ 1 In that case, increasing the number of particles, P, is “better” than sampling M stochastic systems each of size P. We can build a MLMC estimator using dependent samples. Use mean-field theory instead of Monte Carlo theory. We can also use the interaction radius as a discretization parameter and apply MIMC.

30/0