Multi-Index Monte Carlo and Multi-Index Stochastic Collocation Ra - - PowerPoint PPT Presentation

Multi-Index Monte Carlo and Multi-Index Stochastic Collocation

Ra´ ul Tempone

Alexander von Humboldt Professor, RWTH Aachen, Germany and KAUST, Saudi Arabia. Course material developed together with Fabio Nobile (CSQI-MATHICSE, ´ Ecole Polytechnique F´ ed´ erale de Lausanne, Switzerland)

DCSE Fall School on ROM and UQ, November 4–8, 2019

1 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone

Part I Monte Carlo, Multilevel Monte Carlo, Multi-index Monte Carlo methods

2 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Variance Reduction techniques – Control variates

Outline

Variance Reduction techniques – Control variates Multi Level Monte Carlo (MLMC) Main idea and complexity analysis Further optimization of MLMC Normality of MLMC estimator CMLMC Optimal hierarchies Conclusions Multi-Index Monte Carlo General Framework Choosing the Multi-Index Set in MIMC Main Theorem Numerical Results Conclusions Multi-index Stochastic Collocation Complexity analysis Infinite dimensional case (parametric regularity)

3 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Variance Reduction techniques – Control variates

Variance reduction

The Monte Carlo Error is E[Y ] − 1 M

M

• m=1

Y (ωm) ≈ co

• Var[Y ]

M Idea: introduce techniques to ”reduce” Var[Y ] while keeping E[Y ] unchanged. In practice, we look for a new random variable Z such that E[Z] = E[Y ] Var[Z] ≪ Var[Y ] Then, we apply Monte Carlo to the variable Z instead of Y E[Y ] − 1 M

M

• m=1

Z(ωm) ≈ co

• Var [Z]

M ≪ co

• Var[Y ]

M Remark: When adopting these techniques the Monte Carlo rate O(1/ √ M) is unchanged, but we improve the multiplicative constant. For these techniques to be efficient, we need to use particular features

• f Y .

4 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Variance Reduction techniques – Control variates

Variance reduction

The Monte Carlo Error is E[Y ] − 1 M

M

• m=1

Y (ωm) ≈ co

• Var[Y ]

M Idea: introduce techniques to ”reduce” Var[Y ] while keeping E[Y ] unchanged. In practice, we look for a new random variable Z such that E[Z] = E[Y ] Var[Z] ≪ Var[Y ] Then, we apply Monte Carlo to the variable Z instead of Y E[Y ] − 1 M

M

• m=1

Z(ωm) ≈ co

• Var [Z]

M ≪ co

• Var[Y ]

M Remark: When adopting these techniques the Monte Carlo rate O(1/ √ M) is unchanged, but we improve the multiplicative constant. For these techniques to be efficient, we need to use particular features

• f Y .

4 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Variance Reduction techniques – Control variates

Variance reduction

The Monte Carlo Error is E[Y ] − 1 M

M

• m=1

Y (ωm) ≈ co

• Var[Y ]

M Idea: introduce techniques to ”reduce” Var[Y ] while keeping E[Y ] unchanged. In practice, we look for a new random variable Z such that E[Z] = E[Y ] Var[Z] ≪ Var[Y ] Then, we apply Monte Carlo to the variable Z instead of Y E[Y ] − 1 M

M

• m=1

Z(ωm) ≈ co

• Var [Z]

M ≪ co

• Var[Y ]

M Remark: When adopting these techniques the Monte Carlo rate O(1/ √ M) is unchanged, but we improve the multiplicative constant. For these techniques to be efficient, we need to use particular features

• f Y .

4 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Variance Reduction techniques – Control variates

Variance reduction

The Monte Carlo Error is E[Y ] − 1 M

M

• m=1

Y (ωm) ≈ co

• Var[Y ]

M Idea: introduce techniques to ”reduce” Var[Y ] while keeping E[Y ] unchanged. In practice, we look for a new random variable Z such that E[Z] = E[Y ] Var[Z] ≪ Var[Y ] Then, we apply Monte Carlo to the variable Z instead of Y E[Y ] − 1 M

M

• m=1

Z(ωm) ≈ co

• Var [Z]

M ≪ co

• Var[Y ]

M Remark: When adopting these techniques the Monte Carlo rate O(1/ √ M) is unchanged, but we improve the multiplicative constant. For these techniques to be efficient, we need to use particular features

• f Y .

4 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Variance Reduction techniques – Control variates

Control Variates

Take any auxiliary random variable X for which we know the exact value

• f E[X], called Control Variate.

Then, we can define the new random variable Z(β) = Y − β (X − E[X]) Clearly E[Z(β)] = E[Y ] − βE[X − E[X]] = E[Y ], so applying Monte Carlo to Z produces an unbiased estimator of E[Y ]. What can we say about Var [Z(β)]? Var [Z(β)] = E[((Y − E[Y ]) − β(X − E[X]))2] = Var[Y ] + β2Var[X] − 2β Cov[Y , X] which is a quadratic function of β. Minimizing over beta yields β∗ = Cov[Y , X] Var[X] and Var[Z](β∗) = Var[Y ]

• 1 − (Cov[Y , X])2

Var[X]Var[Y ]

• < Var[Y ]

5 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Variance Reduction techniques – Control variates

Control Variates

Take any auxiliary random variable X for which we know the exact value

• f E[X], called Control Variate.

Then, we can define the new random variable Z(β) = Y − β (X − E[X]) Clearly E[Z(β)] = E[Y ] − βE[X − E[X]] = E[Y ], so applying Monte Carlo to Z produces an unbiased estimator of E[Y ]. What can we say about Var [Z(β)]? Var [Z(β)] = E[((Y − E[Y ]) − β(X − E[X]))2] = Var[Y ] + β2Var[X] − 2β Cov[Y , X] which is a quadratic function of β. Minimizing over beta yields β∗ = Cov[Y , X] Var[X] and Var[Z](β∗) = Var[Y ]

• 1 − (Cov[Y , X])2

Var[X]Var[Y ]

• < Var[Y ]

5 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Variance Reduction techniques – Control variates

Control Variates

Take any auxiliary random variable X for which we know the exact value

• f E[X], called Control Variate.

Then, we can define the new random variable Z(β) = Y − β (X − E[X]) Clearly E[Z(β)] = E[Y ] − βE[X − E[X]] = E[Y ], so applying Monte Carlo to Z produces an unbiased estimator of E[Y ]. What can we say about Var [Z(β)]? Var [Z(β)] = E[((Y − E[Y ]) − β(X − E[X]))2] = Var[Y ] + β2Var[X] − 2β Cov[Y , X] which is a quadratic function of β. Minimizing over beta yields β∗ = Cov[Y , X] Var[X] and Var[Z](β∗) = Var[Y ]

• 1 − (Cov[Y , X])2

Var[X]Var[Y ]

• < Var[Y ]

5 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Variance Reduction techniques – Control variates

Control Variates

Take any auxiliary random variable X for which we know the exact value

• f E[X], called Control Variate.

Then, we can define the new random variable Z(β) = Y − β (X − E[X]) Clearly E[Z(β)] = E[Y ] − βE[X − E[X]] = E[Y ], so applying Monte Carlo to Z produces an unbiased estimator of E[Y ]. What can we say about Var [Z(β)]? Var [Z(β)] = E[((Y − E[Y ]) − β(X − E[X]))2] = Var[Y ] + β2Var[X] − 2β Cov[Y , X] which is a quadratic function of β. Minimizing over beta yields β∗ = Cov[Y , X] Var[X] and Var[Z](β∗) = Var[Y ]

• 1 − (Cov[Y , X])2

Var[X]Var[Y ]

• < Var[Y ]

5 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Variance Reduction techniques – Control variates

Remarks: We always obtain a variance reduction as long as X is correlated with Y In practice, we can approximate β∗ by using sample covariances and variances (at least for large M) In some cases (see e.g. Multi Level Monte Carlo), if one knows a priori that X is strongly positively (resp. negatively) correlated with Y , then the choice β = 1 (resp. β = −1) works well.

6 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Variance Reduction techniques – Control variates

Remarks: We always obtain a variance reduction as long as X is correlated with Y In practice, we can approximate β∗ by using sample covariances and variances (at least for large M) In some cases (see e.g. Multi Level Monte Carlo), if one knows a priori that X is strongly positively (resp. negatively) correlated with Y , then the choice β = 1 (resp. β = −1) works well.

6 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Variance Reduction techniques – Control variates

Remarks: We always obtain a variance reduction as long as X is correlated with Y In practice, we can approximate β∗ by using sample covariances and variances (at least for large M) In some cases (see e.g. Multi Level Monte Carlo), if one knows a priori that X is strongly positively (resp. negatively) correlated with Y , then the choice β = 1 (resp. β = −1) works well.

6 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Variance Reduction techniques – Control variates

Does control variates pay in terms of computational work?

Assume that the work to generate the pair (Xj, Yj) is a factor (1 + θ) times the work to generate Yj. Method # samples

• Computat. Work

MC

• C

TOL

2 Var[Y ]

• C

TOL

2 Var[Y ] MC + Cont.Var.

• C

TOL

2 Var[Z(β∗)] (1 + θ)

• C

TOL

2 Var[Z(β∗)] We see that the strategy only pays if Var[Y ] > (1 + θ)Var[Z(β∗)] i.e. 1 > (1 + θ)(1 − (ρXY )2) iff (ρXY )2 > θ 1 + θ. Obs: As said before, we need to use some knowledge of Y to find a sufficiently correlated X, i.e. that it has a relatively small θ and a large ρXY . . . Question: Can we use the structure of SDEs and random PDEs for this purpose?

7 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC)

Outline

Variance Reduction techniques – Control variates Multi Level Monte Carlo (MLMC) Main idea and complexity analysis Further optimization of MLMC Normality of MLMC estimator CMLMC Optimal hierarchies Conclusions Multi-Index Monte Carlo General Framework Choosing the Multi-Index Set in MIMC Main Theorem Numerical Results Conclusions Multi-index Stochastic Collocation Complexity analysis Infinite dimensional case (parametric regularity)

8 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

The Multi Level Monte Carlo idea

Consider again a random vector y(ω) = (y1(ω), . . . , yN(ω)) with density ρ(y) : Γ → R+, a Hilbert-space valued function u(y) : Γ → V , u ∈ L2

ρ(Γ; V ), solution

• f a differential problem

a Lipshitz functional Q : V → R, with Q(0) = 0. Goal: Compute E[Q(u(y(ω)))] In practice the solution u(y(ω)) can not be computed exactly and we will approximate the differential problem by some discretization scheme. Let h be the discretization parameter and uh(y(ω)) the approximate solution such that uh(y(ω)) → u(y(ω)) (in the proper norm) as h → 0.

9 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

The Multi Level Monte Carlo idea

Consider again a random vector y(ω) = (y1(ω), . . . , yN(ω)) with density ρ(y) : Γ → R+, a Hilbert-space valued function u(y) : Γ → V , u ∈ L2

ρ(Γ; V ), solution

• f a differential problem

a Lipshitz functional Q : V → R, with Q(0) = 0. Goal: Compute E[Q(u(y(ω)))] In practice the solution u(y(ω)) can not be computed exactly and we will approximate the differential problem by some discretization scheme. Let h be the discretization parameter and uh(y(ω)) the approximate solution such that uh(y(ω)) → u(y(ω)) (in the proper norm) as h → 0.

9 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

The Multi Level Monte Carlo idea

Introduce a sequence of finer and finer discretizations h0 > h1 > . . . , hL and assume that mesh size hL achieves the desired accuracy

h0 h1 hL

.....

Notation: g(ω) = Q(u(y(ω))) (output quantity) gℓ(ω) = Q(uhℓ(y(ω))), ℓ = 0, . . . , L. Goal: compute E[gL] Idea: write the above expectation as a telescopic sum E[gL] = E[g0] +

L

• ℓ=1

E[gℓ − gℓ−1]

10 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

The Multi Level Monte Carlo idea

Introduce a sequence of finer and finer discretizations h0 > h1 > . . . , hL and assume that mesh size hL achieves the desired accuracy

h0 h1 hL

.....

Notation: g(ω) = Q(u(y(ω))) (output quantity) gℓ(ω) = Q(uhℓ(y(ω))), ℓ = 0, . . . , L. Goal: compute E[gL] Idea: write the above expectation as a telescopic sum E[gL] = E[g0] +

L

• ℓ=1

E[gℓ − gℓ−1]

10 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Multi Level Monte Carlo estimator

MLMC estimator A := 1 M0

M0

• j=1

g0(ω0

j ) + L

• ℓ=1

1 Mℓ

Mℓ

• j=1

(gℓ − gℓ−1)(ωℓ

j )

with {ωℓ

j , j = 1, . . . , Mℓ} independent samples of iid replica on each

• level. Notice that E[A] = E[gL].

Mean Square error (MSE) E[(A − E[g])2] = E[(A − E[gL])2]

• Variance (statistical error)

+ (E[gL − g])2

• Bias (discret. error)

Variance of the estimator. Thanks to the independent samples among levels: Var [A] = E[(A − E[gL])2] = Var [g0] M0 +

L

• ℓ=1

Var [gℓ − gℓ−1] Mℓ Key point: Since Var [gℓ − gℓ−1] gets smaller and smaller for ℓ large, one can take Mℓ smaller and smaller Only few samples on the fine grid hL.

11 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Multi Level Monte Carlo estimator

MLMC estimator A := 1 M0

M0

• j=1

g0(ω0

j ) + L

• ℓ=1

1 Mℓ

Mℓ

• j=1

(gℓ − gℓ−1)(ωℓ

j )

with {ωℓ

j , j = 1, . . . , Mℓ} independent samples of iid replica on each

• level. Notice that E[A] = E[gL].

Mean Square error (MSE) E[(A − E[g])2] = E[(A − E[gL])2]

• Variance (statistical error)

+ (E[gL − g])2

• Bias (discret. error)

Variance of the estimator. Thanks to the independent samples among levels: Var [A] = E[(A − E[gL])2] = Var [g0] M0 +

L

• ℓ=1

Var [gℓ − gℓ−1] Mℓ Key point: Since Var [gℓ − gℓ−1] gets smaller and smaller for ℓ large, one can take Mℓ smaller and smaller Only few samples on the fine grid hL.

11 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Multi Level Monte Carlo estimator

MLMC estimator A := 1 M0

M0

• j=1

g0(ω0

j ) + L

• ℓ=1

1 Mℓ

Mℓ

• j=1

(gℓ − gℓ−1)(ωℓ

j )

with {ωℓ

j , j = 1, . . . , Mℓ} independent samples of iid replica on each

• level. Notice that E[A] = E[gL].

Mean Square error (MSE) E[(A − E[g])2] = E[(A − E[gL])2]

• Variance (statistical error)

+ (E[gL − g])2

• Bias (discret. error)

Variance of the estimator. Thanks to the independent samples among levels: Var [A] = E[(A − E[gL])2] = Var [g0] M0 +

L

• ℓ=1

Var [gℓ − gℓ−1] Mℓ Key point: Since Var [gℓ − gℓ−1] gets smaller and smaller for ℓ large, one can take Mℓ smaller and smaller Only few samples on the fine grid hL.

11 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Multi Level Monte Carlo estimator

MLMC estimator A := 1 M0

M0

• j=1

g0(ω0

j ) + L

• ℓ=1

1 Mℓ

Mℓ

• j=1

(gℓ − gℓ−1)(ωℓ

j )

with {ωℓ

j , j = 1, . . . , Mℓ} independent samples of iid replica on each

• level. Notice that E[A] = E[gL].

Mean Square error (MSE) E[(A − E[g])2] = E[(A − E[gL])2]

• Variance (statistical error)

+ (E[gL − g])2

• Bias (discret. error)

Variance of the estimator. Thanks to the independent samples among levels: Var [A] = E[(A − E[gL])2] = Var [g0] M0 +

L

• ℓ=1

Var [gℓ − gℓ−1] Mℓ Key point: Since Var [gℓ − gℓ−1] gets smaller and smaller for ℓ large, one can take Mℓ smaller and smaller Only few samples on the fine grid hL.

11 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Variance reduction: MLMC

Recall: With Monte Carlo we have to satisfy Var[AMC] = 1 ML Var[gL] ≈ 1 ML Var[g] ≤ TOL2. Main point: MLMC reduces the variance of the deepest level using samples on coarser ( less expensive) levels! Var[AMLMC] = 1 M0 Var[g0] +

L

• ℓ=1

1 Mℓ Var[∆gℓ] ≤ TOL2. Observe: Level 0 in MLMC is usually determined by both stability and accuracy, i.e. Var[∆g1] << Var[g0] ≈ Var[g] < ∞.

12 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Optimal choice of Mℓ

Let C0: cost of generating one realization of g0 Cℓ: cost of generating one realization of gℓ − gℓ−1, ℓ > 0 V0 = Var [g0] Vℓ = Var [gℓ − gℓ−1], ℓ > 0 Then Total work: WMLMC =

L

• ℓ=0

MℓCℓ, Total variance: Var [A] =

L

• ℓ=0

Vℓ Mℓ . Find optimal {Mℓ} to minimize the cost at a fixed variance level min

{Mℓ} L

• ℓ=0

MℓCℓ subject to

L

• ℓ=0

M−1

ℓ Vℓ ≤ TOL2

If we replace (relaxation) Mℓ by real variables, the optimal solution is Mℓ = TOL−2

• Vℓ

Cℓ

L

• j=0
• VjCj

13 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Optimal choice of Mℓ

Let C0: cost of generating one realization of g0 Cℓ: cost of generating one realization of gℓ − gℓ−1, ℓ > 0 V0 = Var [g0] Vℓ = Var [gℓ − gℓ−1], ℓ > 0 Then Total work: WMLMC =

L

• ℓ=0

MℓCℓ, Total variance: Var [A] =

L

• ℓ=0

Vℓ Mℓ . Find optimal {Mℓ} to minimize the cost at a fixed variance level min

{Mℓ} L

• ℓ=0

MℓCℓ subject to

L

• ℓ=0

M−1

ℓ Vℓ ≤ TOL2

If we replace (relaxation) Mℓ by real variables, the optimal solution is Mℓ = TOL−2

• Vℓ

Cℓ

L

• j=0
• VjCj

13 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Optimal choice of Mℓ

Let C0: cost of generating one realization of g0 Cℓ: cost of generating one realization of gℓ − gℓ−1, ℓ > 0 V0 = Var [g0] Vℓ = Var [gℓ − gℓ−1], ℓ > 0 Then Total work: WMLMC =

L

• ℓ=0

MℓCℓ, Total variance: Var [A] =

L

• ℓ=0

Vℓ Mℓ . Find optimal {Mℓ} to minimize the cost at a fixed variance level min

{Mℓ} L

• ℓ=0

MℓCℓ subject to

L

• ℓ=0

M−1

ℓ Vℓ ≤ TOL2

If we replace (relaxation) Mℓ by real variables, the optimal solution is Mℓ = TOL−2

• Vℓ

Cℓ

L

• j=0
• VjCj

13 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Proof: write the Lagrangian function L(M0, . . . , ML, λ) =

L

• ℓ=0

MℓCℓ − λ

• TOL2 −

L

• j=0

Vj Mj

• Then

∂L ∂Mℓ = Cℓ − λ Vℓ M2

ℓ

= 0, hence Mℓ =

• λ Vℓ

Cℓ .

Replacing in the constraint gives

L

• j=0
• VjCj

λ = TOL2, = ⇒ √ λ = TOL−2

L

• j=0
• VjCj

In practice, one should take the ceiling of the real value Mℓ (important if Mℓ < 1). Optimal sample sizes: Mℓ =     TOL−2

• Vℓ

Cℓ

L

• j=0
• VjCj

    Optimal work: WMLMC ≤ TOL−2  

L

• j=0
• VjCj

 

2

+

L

• ℓ=0

Cℓ

14 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Proof: write the Lagrangian function L(M0, . . . , ML, λ) =

L

• ℓ=0

MℓCℓ − λ

• TOL2 −

L

• j=0

Vj Mj

• Then

∂L ∂Mℓ = Cℓ − λ Vℓ M2

ℓ

= 0, hence Mℓ =

• λ Vℓ

Cℓ .

Replacing in the constraint gives

L

• j=0
• VjCj

λ = TOL2, = ⇒ √ λ = TOL−2

L

• j=0
• VjCj

In practice, one should take the ceiling of the real value Mℓ (important if Mℓ < 1). Optimal sample sizes: Mℓ =     TOL−2

• Vℓ

Cℓ

L

• j=0
• VjCj

    Optimal work: WMLMC ≤ TOL−2  

L

• j=0
• VjCj

 

2

+

L

• ℓ=0

Cℓ

14 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Proof: write the Lagrangian function L(M0, . . . , ML, λ) =

L

• ℓ=0

MℓCℓ − λ

• TOL2 −

L

• j=0

Vj Mj

• Then

∂L ∂Mℓ = Cℓ − λ Vℓ M2

ℓ

= 0, hence Mℓ =

• λ Vℓ

Cℓ .

Replacing in the constraint gives

L

• j=0
• VjCj

λ = TOL2, = ⇒ √ λ = TOL−2

L

• j=0
• VjCj

In practice, one should take the ceiling of the real value Mℓ (important if Mℓ < 1). Optimal sample sizes: Mℓ =     TOL−2

• Vℓ

Cℓ

L

• j=0
• VjCj

    Optimal work: WMLMC ≤ TOL−2  

L

• j=0
• VjCj

 

2

+

L

• ℓ=0

Cℓ

14 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Proof: write the Lagrangian function L(M0, . . . , ML, λ) =

L

• ℓ=0

MℓCℓ − λ

• TOL2 −

L

• j=0

Vj Mj

• Then

∂L ∂Mℓ = Cℓ − λ Vℓ M2

ℓ

= 0, hence Mℓ =

• λ Vℓ

Cℓ .

Replacing in the constraint gives

L

• j=0
• VjCj

λ = TOL2, = ⇒ √ λ = TOL−2

L

• j=0
• VjCj

In practice, one should take the ceiling of the real value Mℓ (important if Mℓ < 1). Optimal sample sizes: Mℓ =     TOL−2

• Vℓ

Cℓ

L

• j=0
• VjCj

    Optimal work: WMLMC ≤ TOL−2  

L

• j=0
• VjCj

 

2

+

L

• ℓ=0

Cℓ

14 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Complexity analysis (error versus cost)

[Giles 2008], [Cliffe, Giles, Scheichl and Teckentrup 2011]

For a problem in Rd (d-dimensional), assume

• 1. hℓ = h0δℓ,

0 < δ < 1 (geometric meshes)

• 2. |E[g − gℓ]| ≤ Cwhα

ℓ (weak convergence rate)

• 3. E[(g − gℓ)2] ≤ Cshβ

ℓ (strong convergence rate)

• 4. Cℓ = Cch−dγ

ℓ

(γ = 3 for direct solver and full matrix; γ ≈ 1 for

• ptimal solver with linear complexity)

Notice that from assumption 3. we have

• 5. Vℓ ≤ Cvhβ

ℓ−1, with Cv = 2Cs(δβ + 1)

Indeed Vℓ = Var [gℓ − gℓ−1] ≤ E[(gℓ−gℓ−1)2] ≤ 2E[(g−gℓ)2]+2E[(g−gℓ−1)2] ≤ 2Cs(δβ+1)hβ

ℓ−1

Moreover one always has β ≤ 2α (typically β = 2α for PDEs with smooth random coefficients)

Indeed by Cauchy-Schwarz |E[g − gℓ]| ≤ E[(g − gℓ)2]

1 2 ≤ Csh β 2

ℓ ,

hence α ≥ β 2

15 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Complexity analysis (error versus cost)

[Giles 2008], [Cliffe, Giles, Scheichl and Teckentrup 2011]

For a problem in Rd (d-dimensional), assume

• 1. hℓ = h0δℓ,

0 < δ < 1 (geometric meshes)

• 2. |E[g − gℓ]| ≤ Cwhα

ℓ (weak convergence rate)

• 3. E[(g − gℓ)2] ≤ Cshβ

ℓ (strong convergence rate)

• 4. Cℓ = Cch−dγ

ℓ

(γ = 3 for direct solver and full matrix; γ ≈ 1 for

• ptimal solver with linear complexity)

Notice that from assumption 3. we have

• 5. Vℓ ≤ Cvhβ

ℓ−1, with Cv = 2Cs(δβ + 1)

Indeed Vℓ = Var [gℓ − gℓ−1] ≤ E[(gℓ−gℓ−1)2] ≤ 2E[(g−gℓ)2]+2E[(g−gℓ−1)2] ≤ 2Cs(δβ+1)hβ

ℓ−1

Moreover one always has β ≤ 2α (typically β = 2α for PDEs with smooth random coefficients)

Indeed by Cauchy-Schwarz |E[g − gℓ]| ≤ E[(g − gℓ)2]

1 2 ≤ Csh β 2

ℓ ,

hence α ≥ β 2

15 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Complexity analysis (error versus cost)

[Giles 2008], [Cliffe, Giles, Scheichl and Teckentrup 2011]

For a problem in Rd (d-dimensional), assume

• 1. hℓ = h0δℓ,

0 < δ < 1 (geometric meshes)

• 2. |E[g − gℓ]| ≤ Cwhα

ℓ (weak convergence rate)

• 3. E[(g − gℓ)2] ≤ Cshβ

ℓ (strong convergence rate)

• 4. Cℓ = Cch−dγ

ℓ

(γ = 3 for direct solver and full matrix; γ ≈ 1 for

• ptimal solver with linear complexity)

Notice that from assumption 3. we have

• 5. Vℓ ≤ Cvhβ

ℓ−1, with Cv = 2Cs(δβ + 1)

Indeed Vℓ = Var [gℓ − gℓ−1] ≤ E[(gℓ−gℓ−1)2] ≤ 2E[(g−gℓ)2]+2E[(g−gℓ−1)2] ≤ 2Cs(δβ+1)hβ

ℓ−1

Moreover one always has β ≤ 2α (typically β = 2α for PDEs with smooth random coefficients)

Indeed by Cauchy-Schwarz |E[g − gℓ]| ≤ E[(g − gℓ)2]

1 2 ≤ Csh β 2

ℓ ,

hence α ≥ β 2

15 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Complexity analysis (error versus cost)

Theorem [Cliffe-Giles-Scheichl-Teckentrup 2011] Under the assumptions 1-4 above, if 2α ≥ min(β, dγ), the computational work required to approximate E[g] with MLMC with accuracy 0 < TOL < 1/e in mean square sense, that is E[(A − E[g])2] ≤ TOL2 is bounded as follows: WMLMC ≤ C      TOL−2, for β > dγ, TOL−2 log2(TOL), for β = dγ, TOL−2−(dγ−β)/α, for β < dγ, Remark: standard MC has corresponding complexity WMC ∝ TOL−2−dγ/α.

16 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Complexity analysis (error versus cost)

Theorem [Cliffe-Giles-Scheichl-Teckentrup 2011] Under the assumptions 1-4 above, if 2α ≥ min(β, dγ), the computational work required to approximate E[g] with MLMC with accuracy 0 < TOL < 1/e in mean square sense, that is E[(A − E[g])2] ≤ TOL2 is bounded as follows: WMLMC ≤ C      TOL−2, for β > dγ, TOL−2 log2(TOL), for β = dγ, TOL−2−(dγ−β)/α, for β < dγ, Remark: standard MC has corresponding complexity WMC ∝ TOL−2−dγ/α.

16 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Proof: Enforce a MSE ≤ TOL2 as Bias constraint: |E[g − gL]|2 ≤ 1 2 TOL2, Variance constraint: Var [A] ≤ 1 2 TOL2 From the Bias constraint we get L(TOL) =

• log(

√ 2Cwhα

0 TOL−1)

α log δ−1

• ∼ logδ TOL

1 α

and setting ˜ Cv = Cvhβ

0 and ˜

Cc = Cch−dγ , the total work is (disregarding the term L

j=0 Cj)

WMLMC ≤ TOL−2  

L

• j=1
• ˜

Cv ˜ Ccδj β−dγ

2

 

2

Case β > dγ: WMLMC ≤ TOL−2 ∞

j=1

˜ Cv ˜ Ccδj β−dγ

2

2 ≤ C TOL−2 Case β = dγ: WMLMC ≤ C TOL−2L ≤ C TOL−2(log TOL−1)2 Case β < dγ: WMLMC ≤ CTOL−2δL β−dγ

2

≤ C TOL−2− dγ−β

α

If, moreover, α ≥ min(β, dγ) then the one sample work/level, L

j=0 Cj, is negligible

with respect to the TOL dependent one.

17 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Proof: Enforce a MSE ≤ TOL2 as Bias constraint: |E[g − gL]|2 ≤ 1 2 TOL2, Variance constraint: Var [A] ≤ 1 2 TOL2 From the Bias constraint we get L(TOL) =

• log(

√ 2Cwhα

0 TOL−1)

α log δ−1

• ∼ logδ TOL

1 α

and setting ˜ Cv = Cvhβ

0 and ˜

Cc = Cch−dγ , the total work is (disregarding the term L

j=0 Cj)

WMLMC ≤ TOL−2  

L

• j=1
• ˜

Cv ˜ Ccδj β−dγ

2

 

2

Case β > dγ: WMLMC ≤ TOL−2 ∞

j=1

˜ Cv ˜ Ccδj β−dγ

2

2 ≤ C TOL−2 Case β = dγ: WMLMC ≤ C TOL−2L ≤ C TOL−2(log TOL−1)2 Case β < dγ: WMLMC ≤ CTOL−2δL β−dγ

2

≤ C TOL−2− dγ−β

α

If, moreover, α ≥ min(β, dγ) then the one sample work/level, L

j=0 Cj, is negligible

with respect to the TOL dependent one.

17 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

Proof: Enforce a MSE ≤ TOL2 as Bias constraint: |E[g − gL]|2 ≤ 1 2 TOL2, Variance constraint: Var [A] ≤ 1 2 TOL2 From the Bias constraint we get L(TOL) =

• log(

√ 2Cwhα

0 TOL−1)

α log δ−1

• ∼ logδ TOL

1 α

and setting ˜ Cv = Cvhβ

0 and ˜

Cc = Cch−dγ , the total work is (disregarding the term L

j=0 Cj)

WMLMC ≤ TOL−2  

L

• j=1
• ˜

Cv ˜ Ccδj β−dγ

2

 

2

Case β > dγ: WMLMC ≤ TOL−2 ∞

j=1

˜ Cv ˜ Ccδj β−dγ

2

2 ≤ C TOL−2 Case β = dγ: WMLMC ≤ C TOL−2L ≤ C TOL−2(log TOL−1)2 Case β < dγ: WMLMC ≤ CTOL−2δL β−dγ

2

≤ C TOL−2− dγ−β

α

If, moreover, α ≥ min(β, dγ) then the one sample work/level, L

j=0 Cj, is negligible

with respect to the TOL dependent one.

17 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Main idea and complexity analysis

(Shocking) comments on the MLMC rates

Let us focus on two particular cases: Fastest convergence rate, β > dγ. Here the convergence rate is TOL−2, which is the same of Monte Carlo sampling when the cost to sample each realization is fixed. This means that we do not feel the effect of the discretization in the rates! Smooth noise, β = 2α and β < dγ. Here the resulting convergence rate is TOL−dγ/α, which is the complexity of solving just one realization in the deepest level!

18 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Further optimization of MLMC

MLMC: Beyond geometric hierarchies

The previous result is based on geometric discretizations, i.e. hℓ = h0δℓ. Questions: Can we construct non-geometric discretization MLMC hierarchies? Do they reduce the computational work of MLMC even further? [HvSST] “Implementation and analysis of an adaptive multilevel Monte Carlo algorithm” by H. Hoel, E. von Schwerin, A. Szepessy and R. Tempone, Monte Carlo Methods and Applications. Volume 20, Issue 1,

• pp. 1–41, November 2013.

[ANvST00] “Optimization of mesh hierarchies for Multilevel Monte Carlo”, by A.-L Haji-Ali, F. Nobile, E. von Schwerin and R. Tempone. Stochastic Partial Differential Equations: Analysis and Computations,

• Vol. 4, Issue 1, Pages 76–112, 2016

[NANvST01] “A Continuation Multilevel Monte Carlo”, by N. Collier, A.-L Haji-Ali, F. Nobile, E. von Schwerin and R. Tempone. BIT Numer. Math., 55, pp. 399–432, 2015.

19 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Further optimization of MLMC

Our goal is to compute E [g] where g = g(u). Here g is either a bounded linear functional or a Lips- chitz functional with respect to u, and u solves a stochastic equation. Example: −∇ · (a(x; ω)∇u(x; ω)) = f (x; ω) for x ∈ D = [0, 1]d, u(x; ω) = 0 for x ∈ ∂D, and g(u) =

• D

κ(x)u(x)dx, for sufficiently smooth a, f , κ.

20 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Further optimization of MLMC

Following the standard MLMC approach, we introduce a hierarchy of L + 1 meshes defined by decreasing mesh sizes {hℓ}L

ℓ=0 and we denote

the approximation of g using mesh size hℓ by gℓ. We then write the MLMC estimator as A = 1 M0

M0

• m=1

g0(ω0,m) +

L

• ℓ=1

1 Mℓ

Mℓ

• m=1

(gℓ(ωℓ,m) − gℓ−1(ωℓ,m)) . (2) We assume the following models:

• E [gℓ − g]
• ≈ QW hq1

ℓ ,

(3a) Var[gℓ − gℓ−1] := Vℓ ≈ QShq2

ℓ−1,

(3b) Work per sample of level ℓ := Wℓ ≈ h−dγ

ℓ

. (3c) for some positive constants QW , QS, q1, q2, d and γ.

21 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Further optimization of MLMC

Examples for q1, q2: q1 = q2 = 1 for an SDE with Euler-Maruyama approximation. In our example: a PDE with smooth random coefficients and for piecewise linear or piecewise bilinear continuous finite element approximations we have q1 = 2 and q2 = 4. Examples for γ: γ = 1 for an SDE with Euler-Maruyama approximation. In our PDE example: d = 3 and γ = 3 for a naive Gaussian elimination implementation. Moreover, Using an Iterative solver has γ ≈ 1 while using Direct solver has γ ≈ 1.5. We define: χ = q2 dγ and η = q1 dγ . In our PDE example: χ ≈ 1.34 for iterative solver and χ ≈ 0.89 for direct solver.

22 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Further optimization of MLMC

Examples for q1, q2: q1 = q2 = 1 for an SDE with Euler-Maruyama approximation. In our example: a PDE with smooth random coefficients and for piecewise linear or piecewise bilinear continuous finite element approximations we have q1 = 2 and q2 = 4. Examples for γ: γ = 1 for an SDE with Euler-Maruyama approximation. In our PDE example: d = 3 and γ = 3 for a naive Gaussian elimination implementation. Moreover, Using an Iterative solver has γ ≈ 1 while using Direct solver has γ ≈ 1.5. We define: χ = q2 dγ and η = q1 dγ . In our PDE example: χ ≈ 1.34 for iterative solver and χ ≈ 0.89 for direct solver.

22 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Further optimization of MLMC

Problem (Optimization of computational work)

Given L ∈ N and θ ∈ (0, 1), find H = ({hℓ}L

ℓ=0, {Mℓ}L ℓ=0) ∈ RL+1 +

× RL+1

+

such that W (H) =

L

• ℓ=0

Mℓ hdγ

ℓ

, (4a) is minimized while satisfying the constraints QW hq1

L ≤ (1 − θ)TOL,

(4b) V0 M0 + QS

L

• ℓ=1

hq2

ℓ−1

Mℓ ≤ θTOL Cα 2 , (4c) for a confidence parameter, Cα, such that Φ(Cα) = 1 − α

2 ; here,

0 < α ≪ 1 and Φ is the cumulative distribution function of a standard normal random variable.

23 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Normality of MLMC estimator

Generalized CLT

Theorem (Lindeberg-Feller CLT)

For each n, let Xn,m, for 1 ≤ n ≤ m, be independent random variables (not necessarily identical). Denote an =

n

• m=1

Xn,m, Yn,m = Xn,m − E [Xn,m] , s2

n = n

• m=1

E

• Y 2

n,m

• .

Suppose the following Lindeberg condition is satisfied for all ǫ > 0: lim

n→∞ s−2 n n

• m=1

E

• Y 2

n,m1|Yn,m|>ǫsn

• = 0.

(5) Then, lim

n→∞ P

an − E [an] sn ≤ z

• = Φ(z).

24 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Normality of MLMC estimator

On normality of MLMC estimator

As presented in [NANvST01], [HvSST] has similar results for SDE.

Lemma ([Collier, Haji-Ali, N., von Schwerin, Tempone, 2014])

Consider the MLMC estimator A given by A =

L

• ℓ=0

Mℓ

• m=1

Gℓ(ωℓ,m) Mℓ , where Gℓ(ωℓ,m) denote as usual i.i.d. samples of the random variable Gℓ = gℓ − gℓ−1. The family of random variables, (Gℓ)ℓ≥0, is also assumed

• independent. Denote Yℓ = |Gℓ − E[Gℓ]| and assume the following

C1β−q3ℓ ≤ E[Y 2

ℓ ]

for all ℓ ≥ 0, (6a) E[Y 2+δ

ℓ

] ≤ C2β−τℓ for all ℓ ≥ 0, (6b) for some β > 1 and strictly positive constants C1, C2, q3, δ and τ.

25 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Normality of MLMC estimator

Lemma (cont.)

Choose the number of samples on each level Mℓ to satisfy, for q2 > 0 and a strictly positive sequence {Hℓ}ℓ≥0 Mℓ ≥ β−q2ℓTOL−2H−1

ℓ

L

• ℓ=0

Hℓ

• for all ℓ ≥ 0.

(7) Moreover, choose the number of levels L to satisfy L ≤ max

• 0, c log
• TOL−1

log β + C

• (8)

for some constants C, and c > 0. Finally, denoting p = (1 + δ/2)q3 + (δ/2)q2 − τ, if we have that either p > 0 or c < δ/p, then lim

TOL→0 P

• A − E [A]
• Var[A]

≤ z

• = Φ (z) .

26 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Normality of MLMC estimator

QQ Plot: Experiment verification of normality of MLMC estimator

0.0 0.2 0.4 0.6 0.8 1.0 Empirical CDF 0.0 0.2 0.4 0.6 0.8 1.0 Normal CDF

Direct

0.0 0.2 0.4 0.6 0.8 1.0 Empirical CDF 0.0 0.2 0.4 0.6 0.8 1.0 Normal CDF

Iterative

TOL= 0.005 TOL= 0.00707107

This plot shows that the density of the normalized statistical error is well approximated by a standard normal density.

27 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – CMLMC

Error splitting

Using the asymptotic normality of the MLMC estimator, we can require separately Bias[A] ≈ Cwhα

L ≤ (1 − θ)TOL,

Var [A] ≈ V0 M0 +

L

• ℓ=1

Vℓ Mℓ ≤ θTOL cδ 2 , to guarantee a total error smaller than TOL with probability ≥ 1 − δ. θ ∈ (0, 1) is a free parameter chosen to minimize the computational work. The Bias constraint determines the maximum level L to use Minimization of work under variance constraint leads to optimal sample sizes {Mℓ}L

ℓ=0

To find optimal L and {Mℓ}ℓ one needs good estimates of the rates (α, β, γ), constants (Cw, Cv, Cc) and variances {Vℓ}L

l=0.

28 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – CMLMC

Error splitting

Using the asymptotic normality of the MLMC estimator, we can require separately Bias[A] ≈ Cwhα

L ≤ (1 − θ)TOL,

Var [A] ≈ V0 M0 +

L

• ℓ=1

Vℓ Mℓ ≤ θTOL cδ 2 , to guarantee a total error smaller than TOL with probability ≥ 1 − δ. θ ∈ (0, 1) is a free parameter chosen to minimize the computational work. The Bias constraint determines the maximum level L to use Minimization of work under variance constraint leads to optimal sample sizes {Mℓ}L

ℓ=0

To find optimal L and {Mℓ}ℓ one needs good estimates of the rates (α, β, γ), constants (Cw, Cv, Cc) and variances {Vℓ}L

l=0.

28 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – CMLMC

Error splitting

Using the asymptotic normality of the MLMC estimator, we can require separately Bias[A] ≈ Cwhα

L ≤ (1 − θ)TOL,

Var [A] ≈ V0 M0 +

L

• ℓ=1

Vℓ Mℓ ≤ θTOL cδ 2 , to guarantee a total error smaller than TOL with probability ≥ 1 − δ. θ ∈ (0, 1) is a free parameter chosen to minimize the computational work. The Bias constraint determines the maximum level L to use Minimization of work under variance constraint leads to optimal sample sizes {Mℓ}L

ℓ=0

To find optimal L and {Mℓ}ℓ one needs good estimates of the rates (α, β, γ), constants (Cw, Cv, Cc) and variances {Vℓ}L

l=0.

28 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – CMLMC

A “continuation” adaptive algorithm (CMLMC)

[Haji-Ali, Nobile, Tempone 2014]

Choose a sequence of decreasing tolerances: TOL0 > TOL1 > . . . > TOLK = TOL and an initial guess of the rates (α(0), β(0), γ(0)), constants (C (0)

w , C (0) v , C (0) c ) and variances {Vℓ}L(0) ℓ=0,

FOR k = 1, . . . , K Based on rates (α(k−1), β(k−1), γ(k−1)), constants (C (k−1)

w

, C (k−1)

v

, C (k−1)

c

) and variances {Vℓ}L(k−1)

ℓ=0

compute optimal (L(k), θ(k)) = arg min

θ∈(0,1) L(k−1)≤L≤Lmax

Work(L, θ), s.t. Cwhα

L ≤ (1 − θ)TOLk

compute optimal {M(k)

ℓ }L(k) ℓ=0 to satisfy TOLk.

run MLMC with L(k), {M(k)

ℓ }L(k) ℓ=0

update rates (α(k), β(k), γ(k)), constants (C (k)

w , C (k) v

, C (k)

c

) and variances {Vℓ}L(k)

ℓ=0 based on the new simulations performed

k = k + 1

end FOR

29 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – CMLMC

A “continuation” adaptive algorithm (CMLMC)

[Haji-Ali, Nobile, Tempone 2014]

Choose a sequence of decreasing tolerances: TOL0 > TOL1 > . . . > TOLK = TOL and an initial guess of the rates (α(0), β(0), γ(0)), constants (C (0)

w , C (0) v , C (0) c ) and variances {Vℓ}L(0) ℓ=0,

FOR k = 1, . . . , K Based on rates (α(k−1), β(k−1), γ(k−1)), constants (C (k−1)

w

, C (k−1)

v

, C (k−1)

c

) and variances {Vℓ}L(k−1)

ℓ=0

compute optimal (L(k), θ(k)) = arg min

θ∈(0,1) L(k−1)≤L≤Lmax

Work(L, θ), s.t. Cwhα

L ≤ (1 − θ)TOLk

compute optimal {M(k)

ℓ }L(k) ℓ=0 to satisfy TOLk.

run MLMC with L(k), {M(k)

ℓ }L(k) ℓ=0

update rates (α(k), β(k), γ(k)), constants (C (k)

w , C (k) v

, C (k)

c

) and variances {Vℓ}L(k)

ℓ=0 based on the new simulations performed

k = k + 1

end FOR

29 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – CMLMC

A “continuation” adaptive algorithm (CMLMC)

[Haji-Ali, Nobile, Tempone 2014]

Choose a sequence of decreasing tolerances: TOL0 > TOL1 > . . . > TOLK = TOL and an initial guess of the rates (α(0), β(0), γ(0)), constants (C (0)

w , C (0) v , C (0) c ) and variances {Vℓ}L(0) ℓ=0,

FOR k = 1, . . . , K Based on rates (α(k−1), β(k−1), γ(k−1)), constants (C (k−1)

w

, C (k−1)

v

, C (k−1)

c

) and variances {Vℓ}L(k−1)

ℓ=0

compute optimal (L(k), θ(k)) = arg min

θ∈(0,1) L(k−1)≤L≤Lmax

Work(L, θ), s.t. Cwhα

L ≤ (1 − θ)TOLk

compute optimal {M(k)

ℓ }L(k) ℓ=0 to satisfy TOLk.

run MLMC with L(k), {M(k)

ℓ }L(k) ℓ=0

update rates (α(k), β(k), γ(k)), constants (C (k)

w , C (k) v

, C (k)

c

) and variances {Vℓ}L(k)

ℓ=0 based on the new simulations performed

k = k + 1

end FOR

29 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – CMLMC

A “continuation” adaptive algorithm (CMLMC)

[Haji-Ali, Nobile, Tempone 2014]

Choose a sequence of decreasing tolerances: TOL0 > TOL1 > . . . > TOLK = TOL and an initial guess of the rates (α(0), β(0), γ(0)), constants (C (0)

w , C (0) v , C (0) c ) and variances {Vℓ}L(0) ℓ=0,

FOR k = 1, . . . , K Based on rates (α(k−1), β(k−1), γ(k−1)), constants (C (k−1)

w

, C (k−1)

v

, C (k−1)

c

) and variances {Vℓ}L(k−1)

ℓ=0

compute optimal (L(k), θ(k)) = arg min

θ∈(0,1) L(k−1)≤L≤Lmax

Work(L, θ), s.t. Cwhα

L ≤ (1 − θ)TOLk

compute optimal {M(k)

ℓ }L(k) ℓ=0 to satisfy TOLk.

run MLMC with L(k), {M(k)

ℓ }L(k) ℓ=0

update rates (α(k), β(k), γ(k)), constants (C (k)

w , C (k) v

, C (k)

c

) and variances {Vℓ}L(k)

ℓ=0 based on the new simulations performed

k = k + 1

end FOR

29 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – CMLMC

A “continuation” adaptive algorithm (CMLMC)

[Haji-Ali, Nobile, Tempone 2014]

Choose a sequence of decreasing tolerances: TOL0 > TOL1 > . . . > TOLK = TOL and an initial guess of the rates (α(0), β(0), γ(0)), constants (C (0)

w , C (0) v , C (0) c ) and variances {Vℓ}L(0) ℓ=0,

FOR k = 1, . . . , K Based on rates (α(k−1), β(k−1), γ(k−1)), constants (C (k−1)

w

, C (k−1)

v

, C (k−1)

c

) and variances {Vℓ}L(k−1)

ℓ=0

compute optimal (L(k), θ(k)) = arg min

θ∈(0,1) L(k−1)≤L≤Lmax

Work(L, θ), s.t. Cwhα

L ≤ (1 − θ)TOLk

compute optimal {M(k)

ℓ }L(k) ℓ=0 to satisfy TOLk.

run MLMC with L(k), {M(k)

ℓ }L(k) ℓ=0

update rates (α(k), β(k), γ(k)), constants (C (k)

w , C (k) v

, C (k)

c

) and variances {Vℓ}L(k)

ℓ=0 based on the new simulations performed

k = k + 1

end FOR

29 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – CMLMC

An adaptive algorithm (Continuation-MLMC)

A critical issue is how to estimate variances Vℓ on fine levels where very few simulations are run. In [Haji-Ali, N., Tempone 2014] we have proposed a Bayesian update: prior model Vℓ ≈ Cvhβ

ℓ ,

30 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – CMLMC

An adaptive algorithm (Continuation-MLMC)

A critical issue is how to estimate variances Vℓ on fine levels where very few simulations are run. In [Haji-Ali, N., Tempone 2014] we have proposed a Bayesian update: prior model Vℓ ≈ Cvhβ

ℓ ,

30 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – CMLMC

[NANvST01], Given a hierarchy {hℓ}L

ℓ=0 and the previous models with

certain estimates, it is easy to find some optimal parameters Splitting parameter is given in terms of hL as θ = 1 − QW hq1

L

TOL .

Optimal number of samples in R M∗

ℓ =

• Cα

θTOL 2 Vℓ Wℓ L

• ℓ=0
• VℓWℓ
• .

(9) Of course, Mℓ ∈ N. We can take Mℓ = ⌈M∗

ℓ ⌉ ≤ M∗ ℓ + 1.

Total work =

L

• ℓ=0

MℓWℓ ≤

L

• ℓ=0

M∗

ℓ Wℓ + L

• ℓ=0

Wℓ. However, under certain conditions (Recall 2q1 ≥ min(q2, dγ)) the first term dominates the second one as TOL → 0. Hence, we can consider Mℓ ≈ M∗

ℓ to get the same work estimate, approximately.

Optimal L ∈ N can be found with brute-force optimization in some limited range.

31 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – CMLMC

We propose a Continuation MLMC algorithm that, given a hierarchy, solves the given approximation problem for a sequence of decreasing tolerances, ending with the desired one. The sequence is chosen such that the total work is dominated by the last iteration.

1: Compute with an initial hierarchy. 2: Estimate problem parameters {Vℓ}L

ℓ=0 , QS, QW , q1 and q2.

3: Set i = 0. 4: repeat 5:

Find optimal integer L.

6:

Generate hierarchy {hℓ}L

ℓ=0 for TOLi.

7:

Using TOLi and the variance estimates {Vℓ}L

ℓ=0 and opti-

mal θ, compute {Mℓ}L

ℓ=0 according to (9).

8:

Compute with the resulting MLMC hierarchy.

9:

Estimate problem parameters, {Vℓ}L

ℓ=0 , QS, QW , q1 and q2.

10:

Estimate the total error.

11:

Set i = i + 1

12: until Total error estimate is less than TOL

32 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – CMLMC

Estimating q1, q2, QS and QW is done using a Bayesian approach that allows assuming a prior on q1 and q2. Moreover, samples from multiple levels are used to estimate these parameters. Estimating Vℓ for ℓ > 0 is also done using a Bayesian approach with the models (3) as priors. For r1 > 1 and r2 > 1 (e.g. r1 = 2 and r2 = 1.1) we choose TOLi =

• r iE −i

1

r −1

2 TOL

i < iE, r iE −i

2

r −1

2 TOL

i ≥ iE, Here iE =

• − log(TOL)+log(r2)+log(TOLmax)

log(r1)

• , imposes

TOL0 = TOLmax for some maximum tolerance. Moreover TOLiE = r−1

2 TOL.

33 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – CMLMC

Error plot of CMLMC

10−3 10−2 10−1 TOL 10−6 10−5 10−4 10−3 10−2 10−1 Error

Direct

2 3 2 1 10−3 10−2 10−1 TOL 10−5 10−4 10−3 10−2 10−1 Error

Iterative

3 3 1 1 1 2 TOL

The algorithm was run with Cα = 2 so that the bound holds with 95% confidence.

34 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – CMLMC

Total work of CMLMC

10−3 10−2 10−1 TOL 10−1 100 101 102 103 Time (sec.)

Direct

10−3 10−2 10−1 TOL 10−1 100 101 102 103 Time (sec.)

Iterative

Reference Algorithm

Reference lines are TOL−2.25 and TOL−2, respectively. This is consistent with Corollary 7.

35 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – CMLMC

10−3 10−2 10−1 TOL 1 2 3 4 5 6 7 8 Normalized running time

Direct

10−3 10−2 10−1 TOL 1 2 3 4 5 6 7 Normalized running time

Iterative

CMLMC CMLMC (θ = 0.5) SMLMC( M = 3, θ = 0.5) SMLMC( M = 25, θ = 0.5)

Improvement in running time due to better choice of splitting parameter, θ.

36 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – CMLMC

10−3 10−2 10−1 TOL 1 2 3 4 5 6 7 Normalized running time

Direct

10−3 10−2 10−1 TOL 1 2 3 4 5 6 7 8 Normalized running time

Iterative

CMLMC with reuse CMLMC without reuse SMLMC( M = 3, θ = 0.5) with reuse SMLMC( M = 3, θ = 0.5) without reuse

Reusing samples in CMLMC does not significantly improve running, since the work is dominated by the work of the last iteration.

37 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Optimal hierarchies

Optimal hierarchies

The optimal hierarchies are presented in [ANvST00].

Theorem (On the optimal hierarchies when χ = 1)

For any fixed L ∈ N, with χ = 1, the optimal sequences {hℓ}L

ℓ=0 and {Mℓ}L ℓ=0 in

Problem 1 are given by hℓ = β(ℓ−L) (1 − θ)TOL QW 1

q1 ,

for l = 0, 1, 2, . . . , L, and the optimal choice of the splitting parameter is θ(1, η, L) =

• 1 + 1

2η 1 L + 1 −1 → 1 as L → ∞. For this case the optimal number of levels, L, satisfies asymptotically lim

TOL→0

L + 1 log TOL−1 = 1 2η .

38 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Optimal hierarchies

Theorem (On the optimal hierarchies when χ = 1)

For any fixed L ∈ N, with χ = 1, the optimal sequences, {hℓ}L

ℓ=0 in Problem 1 are

given by hℓ(θ, L) = (1 − θ) TOL QW 1

q1

1−χℓ+1 1−χL+1 V0

QS 1

dγ

χℓ−χL 1−χL+1

· χ

− 1

dγ

2 1−χ

• χL+1−χℓ+1

1−χL+1

+

L(1−χℓ+1)−ℓ(1−χL+1)

1−χL+1

• ,

where the optimal choice of the splitting parameter is θ(χ, η, L) =

• 1 + 1

2η 1 − χ 1 − χL+1 −1 →

• 1 + 1 − min(χ, 1)

2η −1 as L → ∞. For this case the optimal number of levels, L, satisfies asymptotically 1 2η χ − 1 log χ ≤ lim inf

TOL→0

L + 1 log (TOL−1) ≤ lim sup

TOL→0

L + 1 log (TOL−1) ≤ max {1, χ} 2η χ − 1 log χ .

39 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Optimal hierarchies

Corollary (On the asymptotic work with optimal hierarchies)

For the these optimal hierarchies and using asymptotic upper bounds on L, the total computational complexity W (H) TOL−2

• 1+ 1−χ

2η

→ C0,

as TOL ց 0 for 0 < χ < 1 W (H) TOL−2(log(TOL))2 → C1, as TOL ց 0 for χ = 1, W (H) TOL−2 → C2, as TOL ց 0 for χ > 1, with known constants of proportionality, C0 = C 2

αQSQ 1−χ

η

• W

χ

• − 2χ

1−χ

1

2η 2 1 + 2η 1 − χ 2

• 1+ 1−χ

2η

• ,

C1 = C 2

αQS exp(2)

1 2η 2 , C2 = C 2

αV0 χ−1

χ

• Q
• 1

χ

• S

χ2

• χ

χ−1

• (χ − 1)−2 .

40 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Optimal hierarchies

For χ = 1, the hierarchies are not geometric in general hℓ+1 hℓ = V0 QS −

(χ−1)χℓ dγ(χ1+L−1)

χ

2 dγ ( 1 1−χ + (L+1)χℓ+1 χL+1−1 )

(1 − θ)TOL QW (χ−1)χℓ+1

q1(χL+1−1)

. However, we argue that for a nearby L ∈ R (not necessarily optimal,

• r applicable) that these hierarchies become geometric.

Moreover, for geometric hierarchies hℓ = h0βℓ (not necessarily nested) with β =

• χ

2 dγ(1−χ) ,

if χ ∈ R+ \ {1}, exp

• − 2

q2

• ,

if χ = 1, and h0 =

• V0

QS

1

q2 χ 1 dγ(1−χ) , if χ > 1, the asymptotic computational

complexity is the same as the computational complexity of the

• ptimized hierarchies.

41 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Optimal hierarchies

10−7 10−6 10−5 10−4 10−3 10−2 10−1 TOL 0.9 1.0 1.1 1.2 1.3 Normalized Work

χ = 1.33333333333

10−7 10−6 10−5 10−4 10−3 10−2 10−1 TOL 0.9 1.0 1.1 1.2 1.3 1.4 Normalized Work

χ = 0.888888888889

Real-valued optimized hierarchies Integer-valued hierarchies Local integer optimization Real-valued geometric hierarchies integer-valued geometric hierarchies

The number of levels in each hierarchy was chosen using brute-force

• ptimization.

42 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi Level Monte Carlo (MLMC) – Conclusions

Conclusions

Showed normality of MLMC estimator under certain conditions through the use of Lindeberg central limit theorem. We use this in the formulation of our MLMC algorithm and the work optimization problem to control the error in weak sense. Computational saving through better tolerance splitting between bias and statistical error contributions. A more stable continuation MLMC algorithm to learn problem parameters with a small overhead. In CMLMC, reusing samples does not introduce significant computational savings. We show that geometric hierarchies, which are not integer subdivisions, are near-optimal. We derive the MLMC computational complexity with known rates and multiplicative constants. The knowledge of these constants may be used to improve the algorithm further.

43 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo

Outline

Variance Reduction techniques – Control variates Multi Level Monte Carlo (MLMC) Main idea and complexity analysis Further optimization of MLMC Normality of MLMC estimator CMLMC Optimal hierarchies Conclusions Multi-Index Monte Carlo General Framework Choosing the Multi-Index Set in MIMC Main Theorem Numerical Results Conclusions Multi-index Stochastic Collocation Complexity analysis Infinite dimensional case (parametric regularity)

44 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo

Variance reduction: MLMC

45 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo

Variance reduction: further potential

45 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – General Framework

MIMC Estimator

Consider discretization parameters possibly different in each direction hi,αi = hi,0 β−αi

i

with βi > 1. For a multi-index α = (αi)d

i=1 ∈ Nd, we denote by g α the

approximation of g calculated using a discretization defined by α. For i = 1, . . . , d, define the first order difference operators ∆igα =

• gα

if αi = 0, gα − gα−ei if αi > 0, and construct the first order mixed difference ∆gα =

• ⊗d

i=1∆i

• gα =
• j∈{0,1}d

(−1)|j|gα−j with |j| = d

i=1 ji. Requires 2d evaluations of g on different grids.

46 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – General Framework

MIMC Estimator

Consider discretization parameters possibly different in each direction hi,αi = hi,0 β−αi

i

with βi > 1. For a multi-index α = (αi)d

i=1 ∈ Nd, we denote by g α the

approximation of g calculated using a discretization defined by α. For i = 1, . . . , d, define the first order difference operators ∆igα =

• gα

if αi = 0, gα − gα−ei if αi > 0, and construct the first order mixed difference ∆gα =

• ⊗d

i=1∆i

• gα =
• j∈{0,1}d

(−1)|j|gα−j with |j| = d

i=1 ji. Requires 2d evaluations of g on different grids.

46 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – General Framework

MIMC Estimator

Consider discretization parameters possibly different in each direction hi,αi = hi,0 β−αi

i

with βi > 1. For a multi-index α = (αi)d

i=1 ∈ Nd, we denote by g α the

approximation of g calculated using a discretization defined by α. For i = 1, . . . , d, define the first order difference operators ∆igα =

• gα

if αi = 0, gα − gα−ei if αi > 0, and construct the first order mixed difference ∆gα =

• ⊗d

i=1∆i

• gα =
• j∈{0,1}d

(−1)|j|gα−j with |j| = d

i=1 ji. Requires 2d evaluations of g on different grids.

46 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – General Framework

Left: Tensor domain, cylinder. Center: Non-tensor domain immersed in a tensor box. Right: Non-tensor domain with a structured mesh.

47 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – General Framework

Example: Computing gα in d = 2

For α = (α1, α2), we have ∆g(α1,α2) = ∆2(∆1g(α1,α2)) = ∆2 (gα1,α2 − gα1−1,α2) = gα1,α2 − gα1−1,α2 − gα1,α2−1 + gα1−1,α2−1. Notice that in general, ∆gα requires 2d eval- uations of g at different discretization pa- rameters, the largest work of which corre- sponds precisely to the index appearing in ∆gα, namely α. α1 α2

48 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – General Framework

MIMC Estimator

Then, assuming that E [gα] → E [g] as αi → ∞ for all i = 1, . . . , d, it is not difficult to see that E [g] =

• α∈Nd

E [∆gα] ≈

• α∈I

E [∆gα] where I ⊂ Nd is a properly chosen index set. As in MLMC, appoximating each term by independent MC samplers, the MIMC estimator can be written as AMIMC =

• α∈I

1 Mα

Mα

• m=1

∆gα(ωα,m) with properly chosen sample sizes (Mα)α∈I.

49 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – General Framework

MIMC Estimator

Then, assuming that E [gα] → E [g] as αi → ∞ for all i = 1, . . . , d, it is not difficult to see that E [g] =

• α∈Nd

E [∆gα] ≈

• α∈I

E [∆gα] where I ⊂ Nd is a properly chosen index set. As in MLMC, appoximating each term by independent MC samplers, the MIMC estimator can be written as AMIMC =

• α∈I

1 Mα

Mα

• m=1

∆gα(ωα,m) with properly chosen sample sizes (Mα)α∈I.

49 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – General Framework

MIMC Estimator

Then, assuming that E [gα] → E [g] as αi → ∞ for all i = 1, . . . , d, it is not difficult to see that E [g] =

• α∈Nd

E [∆gα] ≈

• α∈I

E [∆gα] where I ⊂ Nd is a properly chosen index set. As in MLMC, appoximating each term by independent MC samplers, the MIMC estimator can be written as AMIMC =

• α∈I

1 Mα

Mα

• m=1

∆gα(ωα,m) with properly chosen sample sizes (Mα)α∈I.

49 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – General Framework

Our objective is to build an estimator A = AMIMC where P(|A − E [g] | ≤ TOL) ≥ 1 − ǫ (10) 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 − g] | ≤ (1 − θ)TOL, (11) Statistical constraint: P (|A − E [A] | ≤ θTOL) ≥ 1 − ǫ. (12) For a given fixed θ ∈ (0, 1). Moreover, motivated by the asymptotic normality of the estimator, A, we approximate (12) by Var[A] ≤ θTOL Cǫ 2 . (13) Here, 0 < Cǫ is such that Φ(Cǫ) = 1 − ǫ

2, where Φ is the cumulative

distribution function of a standard normal random variable.

50 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Choosing the Multi-Index Set in MIMC

Given variance and work estimates we can already optimize for the

• ptimal number of samples M∗

α ∈ R to satisfy the variance constraint

(13) M∗

α =

• Cǫ

θTOL 2 Vα Wα

• α∈I
• VαWα
• .

Taking M∗

α ≤ Mα ≤ M∗ α + 1 such that Mα ∈ N and substituting in the

total work gives Work(I) ≤

• Cǫ

θTOL 2

α∈I

• VαWα

2 +

• α∈I

Wα

• Min. cost of I

. Observe:The work now depends on the index set, I, only.

51 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Choosing the Multi-Index Set in MIMC

Choosing the Index set I – Tensor Product sets

An obvious choice of I is the Full Tensor index-set I(L) = {α ∈ Nd : αi ≤ Li for i ∈ {1 · · · d}} for some L ∈ Rd. It turns out, unsurprisingly, that Full Tensor (FT) index-sets impose restrictive conditions

• n the weak rates wi and yield sub-optimal

complexity rates. L1 L2 α1 α2 Remark: The Bias = E

• g −

α∈IL gα

• = E [g − gL] corresponds to

the discretization error on a (possibly anisotropic) full tensor grid of level L = (L1, . . . , Ld).

52 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Choosing the Multi-Index Set in MIMC

Choosing the Index set I – Tensor Product sets

An obvious choice of I is the Full Tensor index-set I(L) = {α ∈ Nd : αi ≤ Li for i ∈ {1 · · · d}} for some L ∈ Rd. It turns out, unsurprisingly, that Full Tensor (FT) index-sets impose restrictive conditions

• n the weak rates wi and yield sub-optimal

complexity rates. L1 L2 α1 α2 Remark: The Bias = E

• g −

α∈IL gα

• = E [g − gL] corresponds to

the discretization error on a (possibly anisotropic) full tensor grid of level L = (L1, . . . , Ld).

52 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Choosing the Multi-Index Set in MIMC

Choosing the Index set I – Tensor Product sets

An obvious choice of I is the Full Tensor index-set I(L) = {α ∈ Nd : αi ≤ Li for i ∈ {1 · · · d}} for some L ∈ Rd. It turns out, unsurprisingly, that Full Tensor (FT) index-sets impose restrictive conditions

• n the weak rates wi and yield sub-optimal

complexity rates. L1 L2 α1 α2 Remark: The Bias = E

• g −

α∈IL gα

• = E [g − gL] corresponds to

the discretization error on a (possibly anisotropic) full tensor grid of level L = (L1, . . . , Ld).

52 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Choosing the Multi-Index Set in MIMC

MIMC general analysis framework

Question: How do we find optimal index set I for MIMC? min

I⊂Nd Work(I)

such that Bias =

• α/

∈I

Eα ≤ (1 − θ)TOL, Assumption: MIMC work is not dominated by the work to compute a single sample corresponding to each α. Then, minimizing equivalently

• Work(I), the previous min problem can

be recast into a knapsack problem with profits defined for each multi-index α. The corresponding α profit is Pα = Bias contribution Work contribution = Eα √VαWα

53 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Choosing the Multi-Index Set in MIMC

Bias for MIMC

What is the Bias in a MIMC approximation? Bias = E

• g −
• α∈I

gα

• =
• α/

∈I

E [gα] It corresponds to the discretization error on a sparse grid by the so called Combination Technique [Griebel-Schneider-Zenger ’91]

(from [Burgartz-Griebel, Acta Num. ’04])

54 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Choosing the Multi-Index Set in MIMC

Bias for MIMC

What is the Bias in a MIMC approximation? Bias = E

• g −
• α∈I

gα

• =
• α/

∈I

E [gα] It corresponds to the discretization error on a sparse grid by the so called Combination Technique [Griebel-Schneider-Zenger ’91]

(from [Burgartz-Griebel, Acta Num. ’04])

54 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Choosing the Multi-Index Set in MIMC

MIMC general analysis framework

Define the total error associated with an index-set I as E(I) =

• α/

∈I

Eα and the corresponding total work estimate as W(I) =

• α∈I
• VαWα.

Then we can show the following optimality result with respect to E(I) and W(I), namely:

Lemma (Optimal profit sets)

The index-set I(ν) = {α ∈ Nd : Pα ≥ ν} for Pα =

Eα √VαWα is optimal in the sense that any other index-set, ˜

I, with smaller work, W(˜ I) < W(I(ν)), leads to a larger error, E(˜ I) > E(I(ν)).

55 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Choosing the Multi-Index Set in MIMC

MIMC general analysis framework

Once the shape of I is determined, we find I(TOL) to be the minimum set of the family that satisfies the bias constraint E(I(TOL)) =

• α/

∈I(TOL)

Eα ≤ (1 − θ)TOL yielding the corresponding computational work

• Cǫ

θTOL 2  

• α∈I(TOL)
• VαWα

 

2

TOL−(2+p) with p ≥ 0 and possibly some multiplicative log factors in the above

• estimate. To get sharper complexity results we need particular, problem

dependent, assumptions, as we see next.

56 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Choosing the Multi-Index Set in MIMC

Particular Assumptions for MIMC

For every α, we assume the following Assumption 1 (Bias) : Eα = |E [∆gα]| d

i=1 β−αiwi i

Assumption 2 (Variance) : Vα = Var[∆gα] d

i=1 β−αisi i

, Assumption 3 (Work) : Wα = Work(∆gα) d

i=1 βαiγi i

, For positive constants γi, wi, si ≤ 2wi and for i = 1 . . . d. Work(MIMC) =

• α∈I

MαWα

• α∈I

Mα d

i=1 βαiγi i

• .

57 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Choosing the Multi-Index Set in MIMC

Remark on product rate particular assumptions

The Assumptions 1 & 2 in MIMC that the rates for each α are products of 1D rates Eα ∝ d

i=1 β−αiwi i

, Vα ∝ d

i=1 β−αisi i

are stronger than the corresponding assumptions in MLMC. They imply existence of mixed derivatives of the solution of the PDE (and possibly the solution of the adjoint problem associated to the functional Ψ), as opposed to standard derivatives for MLMC.

58 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Choosing the Multi-Index Set in MIMC

Remark on product rate particular assumptions

The Assumptions 1 & 2 in MIMC that the rates for each α are products of 1D rates Eα ∝ d

i=1 β−αiwi i

, Vα ∝ d

i=1 β−αisi i

are stronger than the corresponding assumptions in MLMC. They imply existence of mixed derivatives of the solution of the PDE (and possibly the solution of the adjoint problem associated to the functional Ψ), as opposed to standard derivatives for MLMC.

58 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Choosing the Multi-Index Set in MIMC

Defining the optimal index-set for MIMC

Under Assumptions 1-3 (assuming that they hold as sharp estimates and not just as upper bounds) we have Pα = d

i=1 β −(wi+ γi −si

2

)αi i

= e− d

i=1 log(βi)(wi+ γi −si 2

)αi = e−Cδ d

i=1 δiαi

with δi = log(βi)(wi + γi−si

2

) Cδ , and Cδ =

d

• j=1

log(βj)(wj + γj − sj 2 ). Observe that 0 < δi ≤ 1, since si ≤ 2wi and γi > 0. Moreover,

d

• i=1

δi = 1. Then, for any L ∈ R, the optimal index-set can be written as Iδ(L) = {α ∈ Nd : α · δ =

d

• i=1

αiδi ≤ L}. (14)

59 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Choosing the Multi-Index Set in MIMC

Defining the optimal index-set for MIMC

Under Assumptions 1-3 (assuming that they hold as sharp estimates and not just as upper bounds) we have Pα = d

i=1 β −(wi+ γi −si

2

)αi i

= e− d

i=1 log(βi)(wi+ γi −si 2

)αi = e−Cδ d

i=1 δiαi

with δi = log(βi)(wi + γi−si

2

) Cδ , and Cδ =

d

• j=1

log(βj)(wj + γj − sj 2 ). Observe that 0 < δi ≤ 1, since si ≤ 2wi and γi > 0. Moreover,

d

• i=1

δi = 1. Then, for any L ∈ R, the optimal index-set can be written as Iδ(L) = {α ∈ Nd : α · δ =

d

• i=1

αiδi ≤ L}. (14)

59 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Choosing the Multi-Index Set in MIMC

s L α1 α2

60 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Main Theorem

MIMC work estimate

η = min

i∈{1···d}

log(βi)wi δi , ζ = max

i∈{1···d}

γi − si 2wi , z = #{i ∈ {1 · · · d} : γi − si 2wi = ζ}.

Theorem (Work estimate with optimal weights)

Let the total-degree index set Iδ(L) be given by (14) and (??), taking L = 1 η

• log(TOL−1) + (z − 1) log

1 η log(TOL−1)

• + C
• .

Under Assumptions 1-3, the bias constraint in (11) is satisfied asymptotically and the total work, W (Iδ), of the MIMC estimator, A, subject to the variance constraint (13) satisfies: lim sup

TOL↓0

W (Iδ) TOL−2−2 max(0,ζ) (log (TOL−1))p < ∞, where 0 ≤ p ≤ 3d + 2(d − 1)ζ is known and depends on d, γ, w, s and β.

61 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Main Theorem

Powers of the logarithmic term

ξ = min

i∈{1···d}

2wi − si γi , d2 = #{i ∈ {1 · · · d} : γi = si}, ζ = max

i∈{1···d}

γi − si 2wi , z = #{i ∈ {1 · · · d} : γi − si 2wi = ζ}. Cases for p: A) if ζ ≤ 0 and ζ < ξ,

• r ζ = ξ = 0

then p = 2d2. B) if ζ > 0 and ξ > 0 then p = 2(z − 1)(ζ + 1). C-D) if ζ ≥ 0 and ξ = 0 then p = d − 1 + 2(z − 1)(ζ + 1).

62 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Main Theorem

Fully Isotropic Case: Rough noise case

Assume wi = w, si = s < 2w, βi = β and γi = γ for all i ∈ {1 · · · d}. Then the optimal work is Work(MC) = O

• TOL−2− dγ

w

• .

Work(MLMC) =        O

• TOL−2

, s > dγ, O

• TOL−2

log

• TOL−12

, s = dγ, O

• TOL−(2+ (dγ−s)

w

) , s < dγ. Work(MIMC) =          O

• TOL−2

, s > γ, O

• TOL−2

log

• TOL−12d

, s = γ, O

• TOL−(2+ γ−s

w ) log

• TOL−1(d−1) γ−s

w

• ,

s < γ.

63 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Main Theorem

Fully Isotropic Case: Smooth noise case

Assume wi = w, si = 2w, βi = β and γi = γ for all i ∈ {1 · · · d} and d ≥ 3. Then the optimal work is Work(MC) = O

• TOL−2− dγ

w

• .

Work(MLMC) =        O

• TOL−2

, 2w > dγ, O

• TOL−2

log

• TOL−12

, 2w = dγ, O

• TOL− dγ

w

• ,

2w < dγ. Work(MIMC) =        O

• TOL−2

, 2w > γ, O

• TOL−2

log

• TOL−13(d−1)

, 2w = γ, O

• TOL− γ

w

log

• TOL−1(d−1)(1+γ/w)

, 2w < γ, Up to a multiplicative logarithmic term, Work(MIMC) is the same as solving just a

• ne dimensional deterministic problem.

64 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Main Theorem

MIMC: Case with a single worst direction

Recall ζ = maxi∈{1···d}

γi−si 2wi

and z = #{i ∈ {1 · · · d} : γi−si

2wi

= ζ}. In the special case when ζ > 0 and z = 1, i.e. when the directions are dominated by a single “worst” direction with the maximum difference between the work rate and the rate of variance convergence. In this case, the value of L becomes L = 1 η

• log(TOL−1) + log(C)
• and MIMC with a TD index-set achieves a better rate for the

computational complexity, namely O

• TOL2−2ζ

. In other words, the logarithmic term disappears from the computational complexity. Observe: TD-MIMC with a single worst direction has the same rate of computational complexity as a

• ne-dimensional MLMC along that

single direction.

65 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Numerical Results

Problem description

We test our methods on a three-dimensional, linear elliptic PDE with variable, smooth, stochastic coefficients. The problem is isotropic and we have γi = 2, wi = 2, and si = 4 as TOL → 0.

66 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Numerical Results

Problem description

−∇ · (a(x; ω)∇u(x; ω)) = 1 for x ∈ (0, 1)3, u(x; ω) = 0 for x ∈ ∂(0, 1)3, where a(x; ω) = 1 + exp

• 2Y1Φ121(x) + 2Y2Φ877(x)
• .

Here, Y1 and Y2 are i.i.d. uniform random variables in the range [−1, 1]. We also take Φijk(x) = φi(x1)φj(x2)φk(x3), and φi(x) =

• cos

i

2πx

• i is even,

sin i+1

2 πx

• i is odd,

Finally, the quantity of interest, g, is g = 100

• 2πσ2 −3

2

D

exp

• −x − x02

2

2σ2

• u(x)dx,

and the selected parameters are σ = 0.04 and x0 = [0.5, 0.2, 0.6].

67 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Numerical Results

Numerical test: Computational Errors

10−4 10−3 10−2 10−1 100 TOL 10−9 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 100 Error TOL MIMC TD MIMC FT MLMC

Several runs for different TOL values. Error is satisfied in probability but not

• ver-killed.

68 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Numerical Results

Numerical test: Maximum degrees of freedom

10−4 10−3 10−2 10−1 100 TOL 102 103 104 105 106 107 Maximum DoF TOL−0.5 TOL−1.5 MIMC TD MIMC FT MLMC

Maximum number of degrees of freedom of a sample PDE solve for different TOL

• values. This is

an indication

• f required

memory.

69 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Numerical Results

Numerical test: Running time, 3D problem

10−4 10−3 10−2 10−1 100 TOL 10−2 10−1 100 101 102 103 104 105 106 Running time TOL−2 TOL−3 MIMC TD MIMC FT MLMC

Recall that the work complexity of MC is O

• TOL−5

70 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Numerical Results

Numerical test: Running time, 4D problem

10−4 10−3 10−2 10−1 100 TOL 10−1 100 101 102 103 104 Running time TOL−2 TOL−4 MIMC TD MLMC

A similar PDE problem with d=4 . The work complexity of MC is O

• TOL−6

71 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Numerical Results

Numerical test: QQ-plot

0.0 0.2 0.4 0.6 0.8 1.0 Empirical CDF 0.0 0.2 0.4 0.6 0.8 1.0 Normal CDF TOL = 1.6 × 10−3 TOL = 8 × 10−4

Numerical verification of asymptotic normality of the MIMC estimator. A corresponding statement and proof of the normality of an MIMC estimator can be found in (Haji-Ali et al. 2014).

72 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-Index Monte Carlo – Conclusions

MIMC Conclusions and Extra Points

MIMC is a generalization of MLMC and performs better, especially in higher dimensions, provided mixed regularity between discretization parameters. MIMC general analysis framework, identifying optimal index-set through profit thresholding. Each particular set of regularity assumptions yield its optimal index-set and related complexity.

A MIMC direction does not have to be a spatial dimension. It can represent any form of discretization parameter! Example: 1-DIM Stochastic Particle Systems, MIMC brings complexity down from O(TOL−4) to O(TOL−2 log

• TOL−12).

”A study of Monte Carlo methods for weak approximations of stochastic particle systems in the mean-field”, by A. L. Haji Ali and R. T. May 2016.

Observe, connection to Ensemble Kalman Filter (EnKF): ML-MIMC can compute other statistics, for instance the covariance.

73 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation

Outline

Variance Reduction techniques – Control variates Multi Level Monte Carlo (MLMC) Main idea and complexity analysis Further optimization of MLMC Normality of MLMC estimator CMLMC Optimal hierarchies Conclusions Multi-Index Monte Carlo General Framework Choosing the Multi-Index Set in MIMC Main Theorem Numerical Results Conclusions Multi-index Stochastic Collocation Complexity analysis Infinite dimensional case (parametric regularity)

74 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation

Stochastic and deterministic discretizations

Consider now a quantity of interest g = Q(u(y)) Stochastic approximation by interpolation on a tensor grid: For an interpolation level β = (β1, . . . , βN), denote gβ = E

• U m(β)[g]
• and Hierarchical surplus

∆[gβ] =

N

• n=1

∆jgβ, ∆jgβ =

• gβ − gβ−ej,

if βj > 1, gβ if βj = 1. Deterministic discretization: for an approximation level α = (α1, . . . , αd), corresponding to discretization parameters hi,αi (e.g. hi,αi = hi,02−αi), denote g α(y) = Q(uα(y)) the approximate QoI and Hierarchical surplus ∆[g α] =

d

• i=1

∆ig α, ∆ig α =

• g α − g α−ei,

if αi > 1, g α if αi = 1,

75 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation

Multi Index Stochastic Collocation (MISC)

Idea: build a generalized sparse construction in both approximation levels β and α.

[MISC1]: A.-L. Haji-Ali, F. Nobile, L. Tamellini and R. Tempone. “Multi-Index Stochastic Collocation for random PDEs”, Computer Methods in Applied Mechanics and Engineering, Vol. 306, pp. 95–122, July 2016. [MISC2]: A.-L. Haji-Ali, F. Nobile, L. Tamellini and R. Tempone. “Multi-Index Stochastic Collocation convergence rates for random PDEs with parametric regularity”, arXiv:1511.05393. FoCM, Vol. 16(6), Pages 1555-1605, 2016.

76 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation

MISC main operator

Denote gα,β = E

• U m(β)[g α]
• = E
• U m(β)[Q(uα)]
• We can define Delta operators along the stochastic and deterministic

dimensions ∆d

i gα,β =

• gα,β − gα−ei,β,

if αi > 1, gα,β if αi = 1, ∆s

j gα,β =

• gα,β − gα,β−ej,

if βj > 1, gα,β if βj = 1, Mixed difference operators: deterministic difference : ∆d[gα,β] =

d

• i=1

∆d

i gα,β

stochastic difference : ∆s[gα,β] =

N

• j=1

∆s

j gα,β

77 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation

MISC main operator

Denote gα,β = E

• U m(β)[g α]
• = E
• U m(β)[Q(uα)]
• We can define Delta operators along the stochastic and deterministic

dimensions ∆d

i gα,β =

• gα,β − gα−ei,β,

if αi > 1, gα,β if αi = 1, ∆s

j gα,β =

• gα,β − gα,β−ej,

if βj > 1, gα,β if βj = 1, Mixed difference operators: deterministic difference : ∆d[gα,β] =

d

• i=1

∆d

i gα,β

stochastic difference : ∆s[gα,β] =

N

• j=1

∆s

j gα,β

77 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation

MISC main operator

Denote gα,β = E

• U m(β)[g α]
• = E
• U m(β)[Q(uα)]
• We can define Delta operators along the stochastic and deterministic

dimensions ∆d

i gα,β =

• gα,β − gα−ei,β,

if αi > 1, gα,β if αi = 1, ∆s

j gα,β =

• gα,β − gα,β−ej,

if βj > 1, gα,β if βj = 1, Mixed difference operators: deterministic difference : ∆d[gα,β] =

d

• i=1

∆d

i gα,β

stochastic difference : ∆s[gα,β] =

N

• j=1

∆s

j gα,β

77 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation

MISC estimator

MISC estimator for E [g] AMISC(I) =

• (α,β)∈I

∆s ∆d[gα,β]

• for some index set I ∈ Nd+N.

78 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation

Quasi-optimal index sets

The “best” index sets are given by a knapsack problem [Griebel-Knapek ’09,

Gerstner-Griebel ’03, Bungartz-Griebel ’04]

Error representation |E [g] − AMISC(I)| =

• (α,β)/

∈I

∆s ∆d[gα,β]

• ≤
• (α,β)/

∈I

|∆s∆d[gα,β]| Define Error contribution: ∆Eα,β ≈ |∆s(∆d[gα,β])| Work contribution: ∆Wα,β Knapsack problem max

xα,β

• (α,β)∈Nd+N

xα,β∆Eα,β s.t.

• (α,β)∈Nd+N

xα,β∆Wα,β ≤ Wmax, xα,β ∈ {0, 1}

79 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation

Quasi-optimal index sets

The “best” index sets are given by a knapsack problem [Griebel-Knapek ’09,

Gerstner-Griebel ’03, Bungartz-Griebel ’04]

Error representation |E [g] − AMISC(I)| =

• (α,β)/

∈I

∆s ∆d[gα,β]

• ≤
• (α,β)/

∈I

|∆s∆d[gα,β]| Define Error contribution: ∆Eα,β ≈ |∆s(∆d[gα,β])| Work contribution: ∆Wα,β Knapsack problem max

xα,β

• (α,β)∈Nd+N

xα,β∆Eα,β s.t.

• (α,β)∈Nd+N

xα,β∆Wα,β ≤ Wmax, xα,β ∈ {0, 1}

79 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation

Quasi-optimal sparse grids

The Knapsack problem is, in general, computationally intractable, but its relaxation xα,β ∈ [0, 1] has a simple solution (Dantzig algorithm) Define Profit: Pα,β = ∆Eα,β

∆Wα,β

Sort Pα,β by decreasing profit Add (α, β) to I according to such order until Wmax is reached Quasi-optimal index sets I = I(ǫ) = {(α, β) ∈ Nd+N : Pα,β ≥ ǫ} This construction is valid also for N = ∞ random parameters (distributed fields) as long as Pα,β ≥ ǫ impies that ∃jǫ ∈ N such that βk = 1, ∀k > jǫ. Here it is crucial that m(1) = 1 so that ∆Wα,β ∝ ∞

n=1 m(βn) < ∞.

80 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation

Quasi-optimal sparse grids

The Knapsack problem is, in general, computationally intractable, but its relaxation xα,β ∈ [0, 1] has a simple solution (Dantzig algorithm) Define Profit: Pα,β = ∆Eα,β

∆Wα,β

Sort Pα,β by decreasing profit Add (α, β) to I according to such order until Wmax is reached Quasi-optimal index sets I = I(ǫ) = {(α, β) ∈ Nd+N : Pα,β ≥ ǫ} This construction is valid also for N = ∞ random parameters (distributed fields) as long as Pα,β ≥ ǫ impies that ∃jǫ ∈ N such that βk = 1, ∀k > jǫ. Here it is crucial that m(1) = 1 so that ∆Wα,β ∝ ∞

n=1 m(βn) < ∞.

80 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation

Quasi-optimal sparse grids

The Knapsack problem is, in general, computationally intractable, but its relaxation xα,β ∈ [0, 1] has a simple solution (Dantzig algorithm) Define Profit: Pα,β = ∆Eα,β

∆Wα,β

Sort Pα,β by decreasing profit Add (α, β) to I according to such order until Wmax is reached Quasi-optimal index sets I = I(ǫ) = {(α, β) ∈ Nd+N : Pα,β ≥ ǫ} This construction is valid also for N = ∞ random parameters (distributed fields) as long as Pα,β ≥ ǫ impies that ∃jǫ ∈ N such that βk = 1, ∀k > jǫ. Here it is crucial that m(1) = 1 so that ∆Wα,β ∝ ∞

n=1 m(βn) < ∞.

80 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Complexity analysis

Complexity analysis – N finite

Refinement schemes: stochastic Collocation: m(βj) ∼ 2βj−1. deterministic directions: hi,αi = hi,02−αi Assumptions There exist strictly positive constants QW , Cwork, and rates wi, γi for i = 1 . . . d and gj for j = 1 . . . N, such that ∆Eα,β ≤ QW  

N

• j=1

e−gjm(βj)   d

• i=1

e−wiαi

• ∆Wα,β ≤ Cwork

 

N

• j=1

2βj   d

• i=1

eγiαi

• .

Exponential convergence in the stochastic variables Algebraic convergence in the deterministic variables

81 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Complexity analysis

Complexity analysis – N finite

Refinement schemes: stochastic Collocation: m(βj) ∼ 2βj−1. deterministic directions: hi,αi = hi,02−αi Assumptions There exist strictly positive constants QW , Cwork, and rates wi, γi for i = 1 . . . d and gj for j = 1 . . . N, such that ∆Eα,β ≤ QW  

N

• j=1

e−gjm(βj)   d

• i=1

(hi,αi) ˜

wi

• ∆Wα,β ≤ Cwork

 

N

• j=1

2βj   d

• i=1

(hi,αi)−˜

γi

• .

Exponential convergence in the stochastic variables Algebraic convergence in the deterministic variables

81 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Complexity analysis

Complexity analysis – N finite

Refinement schemes: stochastic Collocation: m(βj) ∼ 2βj−1. deterministic directions: hi,αi = hi,02−αi Assumptions There exist strictly positive constants QW , Cwork, and rates wi, γi for i = 1 . . . d and gj for j = 1 . . . N, such that ∆Eα,β ≤ QW  

N

• j=1

e−gjm(βj)   d

• i=1

e−wiαi

• ∆Wα,β ≤ Cwork

 

N

• j=1

2βj   d

• i=1

eγiαi

• .

Exponential convergence in the stochastic variables Algebraic convergence in the deterministic variables

81 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Complexity analysis

Implication of previous assumptions

The assumption on product rates ∆Eα,β ≤ QW

N

• j=1

e−gjm(βj)

d

• i=1

e−wiαi implies some mixed Sobolev regularity in the deterministic variables “mixed analytic” regularity (Analyticity in poly-ellipses) in the stochastic variables mixed deterministic / stochastic regualrity These assumptions are fulfilled e.g. for the elliptic problem −∇ · (a(x, y) ∇u(x, y)) = f (x) with finite number of random variables a(x, y) = exp{N

j=1 ψj(x)yj}.

82 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Complexity analysis

Implication of previous assumptions

The assumption on product rates ∆Eα,β ≤ QW

N

• j=1

e−gjm(βj)

d

• i=1

e−wiαi implies some mixed Sobolev regularity in the deterministic variables “mixed analytic” regularity (Analyticity in poly-ellipses) in the stochastic variables mixed deterministic / stochastic regualrity These assumptions are fulfilled e.g. for the elliptic problem −∇ · (a(x, y) ∇u(x, y)) = f (x) with finite number of random variables a(x, y) = exp{N

j=1 ψj(x)yj}.

82 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Complexity analysis

Optimal index sets

Under the previous assumptions, the optimal index-sets have the form I(L) = {(α, β) ∈ Nd+N :

d

• i=1

(wi + γi)αi +

N

• j=1

(log(2)βj + gi2βj) ≤ L}

5 10 15 20 25 1 2 3 4 5 6 7 8 α β L=1 L=2 L=3 L=4

index sets for diffusion problem −∇ · (a(y)∇u) = f in 1D (d=1) and with 1 random variable (N=1)

83 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Complexity analysis

MISC Error estimate

Theorem (Error estimate with optimal sets)

[MISC1] Under previous assumptions on the error and work contributions, for any W large enough there exists an index-set I(W ) such that W AMISC(I(W )) ≤ W , Error[AMISC(I(W ))] ≤ CE(N)W −ζ (log (W ))(ζ+1)(z−1). where ζ = mini=1,...,d

wi γi and z = #{i = 1, . . . d : wi γi = ζ}.

Up to logarithmic terms, asymptotically Error ∼ W −ζ: complexity of a 1D deterministic problem (the worst one). The complexity of the Stochastic Collocation is not seen. The asymptotic rate is independent of the number N of random variables (but the constant CE is not!)

84 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Complexity analysis

MISC Error estimate

Theorem (Error estimate with optimal sets)

[MISC1] Under previous assumptions on the error and work contributions, for any W large enough there exists an index-set I(W ) such that W AMISC(I(W )) ≤ W , Error[AMISC(I(W ))] ≤ CE(N)W −ζ (log (W ))(ζ+1)(z−1). where ζ = mini=1,...,d

wi γi and z = #{i = 1, . . . d : wi γi = ζ}.

Up to logarithmic terms, asymptotically Error ∼ W −ζ: complexity of a 1D deterministic problem (the worst one). The complexity of the Stochastic Collocation is not seen. The asymptotic rate is independent of the number N of random variables (but the constant CE is not!)

84 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Complexity analysis

MISC Error estimate

Theorem (Error estimate with optimal sets)

[MISC1] Under previous assumptions on the error and work contributions, for any W large enough there exists an index-set I(W ) such that W AMISC(I(W )) ≤ W , Error[AMISC(I(W ))] ≤ CE(N)W −ζ (log (W ))(ζ+1)(z−1). where ζ = mini=1,...,d

wi γi and z = #{i = 1, . . . d : wi γi = ζ}.

Up to logarithmic terms, asymptotically Error ∼ W −ζ: complexity of a 1D deterministic problem (the worst one). The complexity of the Stochastic Collocation is not seen. The asymptotic rate is independent of the number N of random variables (but the constant CE is not!)

84 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Complexity analysis

Numerical test case and comparison

Problem: −∇ · (a(y)∇u) = 1, in B = [0, 1]d, u = 0, on ∂B a(y) = e

N

n=1

√λnψnyn,

ψn: Fourier modes, yn

iid

∼ U([−1, 1]), λn ∼ e−n Clenshaw-Curtis collocation points; finite difference discretization in space. Comparison MISC with a-priori and a-posteriori construction of optimal set MLSC (uniform refinement in space; only one discretization param.) SGSC: Single level SC (with optimal balance of spatial and stochastic errors) Stochastic Composit Collocation Method (SCC) [Bieri 2011]. Corresponds to the choice I = {(α, β) : α +

N

• n=1

βn ≤ L} Multi-Index Monte Carlo (MIMC)

85 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Complexity analysis

Numerical test case and comparison

Problem: −∇ · (a(y)∇u) = 1, in B = [0, 1]d, u = 0, on ∂B a(y) = e

N

n=1

√λnψnyn,

ψn: Fourier modes, yn

iid

∼ U([−1, 1]), λn ∼ e−n Clenshaw-Curtis collocation points; finite difference discretization in space. Comparison MISC with a-priori and a-posteriori construction of optimal set MLSC (uniform refinement in space; only one discretization param.) SGSC: Single level SC (with optimal balance of spatial and stochastic errors) Stochastic Composit Collocation Method (SCC) [Bieri 2011]. Corresponds to the choice I = {(α, β) : α +

N

• n=1

βn ≤ L} Multi-Index Monte Carlo (MIMC)

85 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Complexity analysis

Numerical comparison

Test case in dimension d = 3 with N = 10 random variables

100 101 102 103 104 105 106 107 108 Work 10−9 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 Error

a-prior MISC a-post MISC SCC MIMC SGSC-QO a-prior MLSC-QO a-post MLSC-QO E = W −2 log(W)6 E = W −0.5

86 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Infinite dimensional case (parametric regularity)

Infinite dimensional case (N = ∞) – Model problem

We apply MISC to the following problem on a hypercube domain B ⊂ Rd −∇ · (a(x, y) ∇u(x, y)) = f (x) in B u(x, y) = 0

• n

∂B, where a(x, y) = eκ(x,y), with κ(x, y) =

• j∈N+

ψj(x)yj. Here, yj

iid

∼ U([−1, 1]). QoI: g = Q(u) Goal: compute E [g] = E [Q(u)] Stochastic collocation on Clenshaw-Curtis points (sparse) piecewise liner finite elements in space.

87 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Infinite dimensional case (parametric regularity)

Infinite dimensional case (N = ∞) – Model problem

We apply MISC to the following problem on a hypercube domain B ⊂ Rd −∇ · (a(x, y) ∇u(x, y)) = f (x) in B u(x, y) = 0

• n

∂B, where a(x, y) = eκ(x,y), with κ(x, y) =

• j∈N+

ψj(x)yj. Here, yj

iid

∼ U([−1, 1]). QoI: g = Q(u) Goal: compute E [g] = E [Q(u)] Stochastic collocation on Clenshaw-Curtis points (sparse) piecewise liner finite elements in space.

87 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Infinite dimensional case (parametric regularity)

Assumptions

There exist smax ∈ N+ and an increasing sequence p0 < p1 < . . . < psmax with psmax < 1

2 such that

• j≥1

ψjps

C s(B) ≤ κps s < ∞,

s = 0, 1, . . . , smax. with κs sufficiently small. B and f are such that u(y) ∈ H1+smax(B) for any y ∈ [−1, 1]N. Then, there exist {gj(s)}j∈N with

j≥1 e−psgj(s) < ∞ for s = 0, . . . , smax

such that ∆Eα,β ≤ QW min

s=0,...,smax

 

∞

• j=1

e−gj(s)m(βj)   d

• i=1

e− min{1, s

d }αi

• 88 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Infinite dimensional case (parametric regularity)

MISC convergence result

Theorem

[MISC2] Under technical assumptions the profit-based MISC estimator built using Stochastic Collocation over Clenshaw-Curtis points and piecewise multilinear finite elements for solving the deterministic problems, we have, for any δ > 0, Error[AMISC(I)] ≤ ˜ CP(δ) W AMISC(I)−rMISC+δ. The rate rMISC is as follows: Case 1 if min{ 1

γ , smax γd } ≤ 1 psmax − 2, then rMISC = min{ 1 γ , smax γd },

Case 2 if min{ 1

γ , smax γd } ≥ 1 psmax − 2, then

rMISC = 1 p0 − 2

• max

s=1,...,smax

• 1

min{ 1

γ , s γd }

1 p0 − 1 ps

• + 1

−1 .

89 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Infinite dimensional case (parametric regularity)

Example: log uniform field with parametric regularity

Here, the regularity of κ = log(a) is determined through ν > 0 κ(x, y) =

• k∈Nd

Ak

• ℓ∈{0,1}d

yk,ℓ

d

• j=1
• cos

π L kjxj ℓj sin π L kjxj 1−ℓj , where the coefficients Ak are taken as Ak= √ 3

• 2

|k|0 2 (1 + |k|2)− ν+d/2 2

. (Same decay as Mat´ ern covariance of parameter ν). We have p0 > ν

d + 1 2

−1 and ps > ν−s

d

+ 1

2

−1. Applying the theorem with optimal s we obtain rMISC = min 1 γ , ν d − 3 2

• 1

1 + γ

• 90 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Infinite dimensional case (parametric regularity)

Example: log uniform field with parametric regularity

Here, the regularity of κ = log(a) is determined through ν > 0 κ(x, y) =

• k∈Nd

Ak

• ℓ∈{0,1}d

yk,ℓ

d

• j=1
• cos

π L kjxj ℓj sin π L kjxj 1−ℓj , where the coefficients Ak are taken as Ak= √ 3

• 2

|k|0 2 (1 + |k|2)− ν+d/2 2

. (Same decay as Mat´ ern covariance of parameter ν). We have p0 > ν

d + 1 2

−1 and ps > ν−s

d

+ 1

2

−1. Applying the theorem with optimal s we obtain rMISC = min 1 γ , ν d − 3 2

• 1

1 + γ

• 90 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Infinite dimensional case (parametric regularity)

Application of main theorem [MISC2]

1 2 3 4 5 6 7 8 9 ν 0.0 0.5 1.0 1.5 2.0 rMISC Deterministic problem

Theory Improved Square summ. Square summ., improved Observed for d = 1

5 10 15 20 25 30 35 ν 0.0 0.5 1.0 1.5 2.0 rMISC Deterministic problem

Theory Improved Square summ. Square summ., improved Observed for d = 3

Error ∝ Work−rMISC (ν,d) A similar analysis shows the corresponding ν-dependent convergence rates

• f MIMC but based on ℓ2 summability of bs and Fernique type of results.

91 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Infinite dimensional case (parametric regularity)

MISC numerical results [MISC2]

100 101 102 103 104 105 106 107 108 Wmax 10−9 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 100 Relative Error

Simplified a-posteriori MIMC E = W −1.61 E = W −0.5

100 101 102 103 104 105 106 Wmax 10−7 10−6 10−5 10−4 10−3 10−2 10−1 100 Relative Error

Simplified a-posteriori MIMC E = W −1.11 E = W −0.5

Left: d = 1, ν = 2.5. Right: d = 3, ν = 4.5. Error ∝ Work−rMISC (ν,d)

92 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Infinite dimensional case (parametric regularity)

MISC numerical results – 1D

100 101 102 103 104 105 106 107 108 Wmax 2 4 6 8 10 12 max (α, β)

Simplified a-posteriori MIMC

100 101 102 103 104 105 106 107 108 Wmax 100 101 102 103

Simplified a-posteriori MIMC

d = 1 and ν = 2.5. Extreme values of α and β included in the MISC set I. Specifically, left-solid is the maximum space discretization level, max(α,β)∈I (max (α)), left-dashed is the maximum quadrature level, max(α,β)∈I (max (β)), right-solid is the index of the last activated random variable, max(α,β)∈I

• maxβj >1 j
• , and right-dashed is the maximum number of jointly

activated variables, max(α,β)∈I (|β − 1|0).

93 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Infinite dimensional case (parametric regularity)

MISC numerical results [MISC2]

100 101 102 103 104 105 106 Wmax 1 2 3 4 5 6 7 8 max (α, β)

Simplified a-posteriori MIMC

100 101 102 103 104 105 106 Wmax 100 101 102 103 104

Simplified a-posteriori MIMC

d = 3 and ν = 4.5. Extreme values of α and β included in the MISC set I. Specifically, left-solid is the maximum space discretization level, max(α,β)∈I (max (α)), left-dashed is the maximum quadrature level, max(α,β)∈I (max (β)), right-solid is the index of the last activated random variable, max(α,β)∈I

• maxβj >1 j
• , and right-dashed is the maximum number of jointly

activated variables, max(α,β)∈I (|β − 1|0).

94 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Infinite dimensional case (parametric regularity)

Deterministic runs, numerical results [MISC2]

These plots shows the non-asymptotic effect of the logarithmic factor for d > 1 (as discussed in [Thm. MISC1]) on the linear convergence fit in log-log scale.

100 101 102 103 104 Wmax 10−9 10−8 10−7 10−6 10−5 10−4 10−3 10−2 Relative Error

Simplified a-posteriori E = W −1.96

100 101 102 103 104 105 106 Wmax 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 Relative Error

Simplified a-posteriori E = W −1.40 E = W −2 log(W)2

Left: d = 1. Right: d = 3.

95 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Infinite dimensional case (parametric regularity)

Acknowledgements

Ivo Babuska Joakim Beck H˙ akon Hoel Abdul Lateef Haji-Ali Quan Long Giovanni Migliorati Mohammad Motamed Marco Scavino Erik von Schwerin Anders Szepessy Lorenzo Tamellini Francesco Tesei Suojin Wang Clayton Webster Georgios Zouraris

Italian project FIRB-IDEAS (’09) Advanced Numerical Techniques for Uncertainty Quantification in Engineering and Life Science Problems Academic Excellency Alliance KAUST–UT Austin project “Predictability and uncertainty quantification for models of porous media” KAUST Strategic Research Initiative (SRI) Center for Uncertainty Quantification in Computational Science and Engineering. Swiss National Science Foundation Project No. 140574 “Efficient numerical methods for flow and transport phenomena in heterogeneous random porous media” 96 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Infinite dimensional case (parametric regularity)

References

Polynomial approximation and sparse grids

• I. Babuˇ

ska, F. Nobile and R. Tempone. A stochastic collocation method for elliptic PDEs with random input data, SIAM Review, 52(2):317–355, 2010

• J. B¨

ack, F. Nobile, L. Tamellini and R. Tempone Stochastic spectral Galerkin and collocation methods for PDEs with random coefficients: a numerical comparison, LNCSE Springer, 76:43-62, 2011

• F. Nobile, L. Tamellini, R. Tempone,

Convergence of quasi-optimal sparse grid approximation of Hilbert-valued functions: application to random elliptic PDEs, Numer. Math. 2015. Published online.

• J. Beck, F. Nobile, L. Tamellini, and R. Tempone,

Convergence of quasi-optimal Stochastic Galerkin methods for a class of PDES with random coefficients, Computers & Mathematics with Applications. Vol 67, No. 4, 732–751, 2014.

• J. Beck, F. Nobile, L. Tamellini, and R. Tempone,

A quasi-optimal sparse grids procedure for groundwater flows, Springer LNCSE, Vol. 95, 1-16, 2014.

• J. B¨

ack and F. Nobile and L. Tamellini and R. Tempone On the optimal polynomial approximation of stochastic PDEs by Galerkin and Collocation methods, Math. Mod. Methods Appli. Sci. (M3AS), 22(9), 2012.

• I. Babuska, R. Tempone and G. Zouraris,

Galerkin finite element approximations of stochastic elliptic partial differential equations. SIAM J. Numer. Anal. 42 (2004), no. 2, 800–825.

97 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Infinite dimensional case (parametric regularity)

References

• F. Nobile, R. Tempone and C. Webster

An anisotropic sparse grid stochastic collocation method for PDEs with random input data, SIAM J. Numer. Anal., 46(5):2411–2442, 2008

• F. Nobile, R. Tempone and C. Webster

A sparse grid stochastic collocation method for partial differential equations with random input data. SIAM J. Numer. Anal. 46(5):2309–2345, 2008.

• F. Nobile and R. Tempone

Analysis and implementation issues for the numerical approximation of parabolic equations with random coefficients, IJNME, 80:979–1006, 2009

• L. Tamellini,

Polynomial approximation of PDEs with stochastic coefficients, PhD Thesis, Politecnico di

• Milano. March 2012.

98 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Infinite dimensional case (parametric regularity)

References

Multi Level/Index Monte Carlo

A.-L. Haji-Ali, F. Nobile, L. Tamellini, and R. Tempone, Multi-index stochastic collocation convergence rates for random PDEs with parametric regularity, MATHICSE Technical Report 29.2015. To appear in FoCM, 2016. A.-L. Haji-Ali, F. Nobile, L. Tamellini, and R. Tempone, Multi-Index Stochastic Collocation for random PDEs, Comput. Methods Appl. Mech. Engrg., vol. 306, pp. 95–122, 2016.

• F. Nobile, F. Tesei,

A multi level Monte Carlo method with control variate for elliptic PDEs with log-normal coefficients, Stoch. PDEs: Anal. & Comp., vol. 3, no. 3, pp. 398–444, 2015. A-L. Haji-Ali, F. Nobile, R. Tempone, Multi index Monte Carlo: when sparsity meets sampling, Numer. Math. vol. 132, no. 4, pp. 767–806, 2016. A-L. Haji-Ali, F. Nobile, E. von Schwerin, R. Tempone, Optimization of mesh hierarchies in multilevel Monte Carlo samplers, SPDEs: theory and applications, vol. 4, no. 1, pp.76–112, 2016.

• N. Collier, A-L. Haji-Ali, F. Nobile, E. von Schwerin, R. Tempone,

A continuation multilevel Monte Carlo algorithm, BIT Numerical Mathematics, 2015, vol 55/2, pp. 399-432.

• C. Bayer, H. Hoel, E. Von Schwerin, and R. Tempone.

On non-asymptotic optimal stopping criteria in Monte Carlo Simulations. Siam J. Sci. Comp., vol. 36, num. 2, pp. A869-A885, 2014.

99 / 101

Monte Carlo, Multilevel Monte Carlo, Multi-Index Monte Carlo Ra´ ul Tempone Multi-index Stochastic Collocation – Infinite dimensional case (parametric regularity)

References

• H. Hoel, E. von Schwerin, A. Szepessy, and R. Tempone.

Implementation and Analysis of an Adaptive Multi Level Monte Carlo Algorithm. Monte Carlo Methods and Applications. Vol. 20, Issue 1, Pages 1–41, March 2014.

• H. Hoel, E. Von Schwerin, A. Szepessy and R. Tempone.

Adaptive Multi Level Monte Carlo Simulation, LNCSE Springer, Vol. 82, pp. 217-234, 2012

100 / 101