Stochastic Modelling of Uncertainty in PDEs and Multilevel Monte - - PowerPoint PPT Presentation

Stochastic Modelling of Uncertainty in PDEs and Multilevel Monte Carlo

Robert Scheichl

Institute for Applied Mathematics & Interdisciplinary Center for Scientific Computing Universit¨ at Heidelberg Tutorial on Uncertainty Quantification – Efficient Methods for PDEs with Random Coefficients

as part of RICAM Special Semester on “Multivariate Algorithms and their Foundations in Number Theory”

Linz, December 14, 2018

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 1 / 61

Outline

Motivation: data uncertainty & stochastic modelling Model problem: subsurface flow and radwaste disposal What are the computational/mathematical challenges? The Curse of Dimensionality Monte Carlo methods Multilevel Monte Carlo Numerical Analysis for PDEs with random coefficients

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 2 / 61

UQ and the Scientific Computing Paradigm

Physical Phenomenon Mathematical Model Numerical Approximation Computer Implementation Prediction Insight Optimization Control Decision

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 3 / 61

UQ and the Scientific Computing Paradigm

Physical Phenomenon Data Quantities of Interest Mathematical Model DEs Parameters Solution Numerical Approximation Discretization Solvers Computer Implementation Software Prediction Insight Optimization Control Decision

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 3 / 61

UQ and the Scientific Computing Paradigm

Physical Phenomenon Uncertain Data Lack of Knowledge Variability Mathematical Model SDEs Random Fields Numerical Approximation ? Computer Implementation ? Prediction Insight Optimization Control Decision

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 3 / 61

UQ and the Scientific Computing Paradigm

Physical Phenomenon Uncertain Data Lack of Knowledge Variability Mathematical Model SDEs Random Fields Numerical Approximation ? Computer Implementation ? Prediction Insight Optimization Control Decision Quantified

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 3 / 61

Stochastic Modelling – Many Reasons

lack of data (e.g. data assimilation for weather prediction) data uncertainty (e.g. uncertainty quantification in subsurface flow) parameter identification (e.g. Bayesian inference in engineering) unresolvable scales (e.g. atmospheric dispersion modelling) high dimensionality (e.g. stochastic simulation in systems biology)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 4 / 61

Stochastic Modelling – Many Reasons

lack of data (e.g. data assimilation for weather prediction) data uncertainty (e.g. uncertainty quantification in subsurface flow) parameter identification (e.g. Bayesian inference in engineering) unresolvable scales (e.g. atmospheric dispersion modelling) high dimensionality (e.g. stochastic simulation in systems biology) Input: best knowledge about system (PDE), input parameter statistics, measured ouput data with error statistics,. . .

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 4 / 61

Stochastic Modelling – Many Reasons

lack of data (e.g. data assimilation for weather prediction) data uncertainty (e.g. uncertainty quantification in subsurface flow) parameter identification (e.g. Bayesian inference in engineering) unresolvable scales (e.g. atmospheric dispersion modelling) high dimensionality (e.g. stochastic simulation in systems biology) Input: best knowledge about system (PDE), input parameter statistics, measured ouput data with error statistics,. . . Output: statistics of QoIs or of entire state space, e.g. Data misfit or rainfall at some location in weather prediction Flow or ’breakthrough’ time at radioactive waste repository Production rate of an oil reservoir Amount of ash over Heathrow Certification of carbon fibre composite wing in aeronautical engin.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 4 / 61

Stochastic Modelling – Many Reasons

lack of data (e.g. data assimilation for weather prediction) data uncertainty (e.g. uncertainty quantification in subsurface flow) parameter identification (e.g. Bayesian inference in engineering) unresolvable scales (e.g. atmospheric dispersion modelling) high dimensionality (e.g. stochastic simulation in systems biology) Input: best knowledge about system (PDE), input parameter statistics, measured ouput data with error statistics,. . . Output: statistics of QoIs or of entire state space, e.g. Data misfit or rainfall at some location in weather prediction Flow or ’breakthrough’ time at radioactive waste repository Production rate of an oil reservoir Amount of ash over Heathrow Certification of carbon fibre composite wing in aeronautical engin. But: Often data very sparse/uncertain ⇒ Need good physical models!

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 4 / 61

Examples – PDEs with Random Coefficients

Navier–Stokes (e.g. flow around wing, weather forecasting): ρ ∂v ∂t + v · ∇v

• = −∇p + µ∇2v + f

in Ω subject to IC v(x, 0) = v0(x) + BCs

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 5 / 61

Examples – PDEs with Random Coefficients

Navier–Stokes (e.g. flow around wing, weather forecasting): ρ(ω) ∂v ∂t + v · ∇v

• = −∇p + µ(ω)∇2v + f

in Ω subject to IC v(x, 0) = v0(x) + BCs

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 5 / 61

Examples – PDEs with Random Coefficients

Navier–Stokes (e.g. flow around wing, weather forecasting): ρ(ω) ∂v ∂t + v · ∇v

• = −∇p + µ(ω)∇2v + f(x, ω) in Ω(ω)

subject to IC v(x, 0) = v0(x, ω) + BCs

uncertain ICs → data assimilation

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 5 / 61

Examples – PDEs with Random Coefficients

Structural Mechanics (e.g. composites, tires or bone): ∇ ·

• C : 1

2

• ∇u + ∇uT

+ F = 0 in Ω subject to BCs

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 6 / 61

Examples – PDEs with Random Coefficients

Structural Mechanics (e.g. composites, tires or bone): ∇ ·

• C(x, ω) : 1

2

• ∇u + ∇uT

+ F = 0 in Ω subject to BCs

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 6 / 61

Examples – PDEs with Random Coefficients

Structural Mechanics (e.g. composites, tires or bone): ∇ ·

• C(x, ω) : 1

2

• ∇u + ∇uT

+ F(x, ω) = 0 in Ω(ω) subject to BCs

fibre defects contact on rough surface

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 6 / 61

Predicted storm path with uncertainty cones

Source: National Hurricane Center, USA

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 7 / 61

Examples

Waves in Heterogeneous Media (e.g. seismic or tsunami):

[Behrens et al]

Subsurface Fluid Flow or mantel convection:

• ptimal well placement

[Gmeiner et al]

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 8 / 61

A Case Study: Radioactive Waste Disposal

Deep geological disposal – favoured by nearly all countries with radwaste

An area where UQ has played central role in last 30 years. Storage in containers in tunnels, several hundred meters deep, in stable geological formations. Several barriers: chemical, physical, geological. Containment must be assured for at least 10,000 years

(uncertainties due to technological complexity & long time scales).

Main escape route for radionuclides: groundwater pathway.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 9 / 61

A Case Study: Radioactive Waste Disposal

Deep geological disposal – favoured by nearly all countries with radwaste

An area where UQ has played central role in last 30 years. Storage in containers in tunnels, several hundred meters deep, in stable geological formations. Several barriers: chemical, physical, geological. Containment must be assured for at least 10,000 years

(uncertainties due to technological complexity & long time scales).

Main escape route for radionuclides: groundwater pathway. Assessing safety of potential sites of utmost importance long timescales → modelling essential! Key aspect: How to quantify uncertainties in the models?

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 9 / 61

A Case Study: Radioactive Waste Disposal

WIPP – Waste Isolation Pilot Plant

US DOE repository for radioactive waste situated near Carlsbad, NM

(fully operational since 1999).

Extensive site characterisation and performance assessment since 1976, in particular compliance certification by US EPA (every 5 years) Lots of publicly available data. http://www.wipp.energy.gov

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 10 / 61

A Case Study: Radioactive Waste Disposal

WIPP – Waste Isolation Pilot Plant

US DOE repository for radioactive waste situated near Carlsbad, NM

(fully operational since 1999).

Extensive site characterisation and performance assessment since 1976, in particular compliance certification by US EPA (every 5 years) Lots of publicly available data. http://www.wipp.energy.gov Repository at 655m depth in bedded evaporites (primarily halite, a salt) Most transmissive rock layer: Culebra Dolomite (principal pathway)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 10 / 61

A Case Study: Radioactive Waste Disposal

WIPP UQ Scenario

Important WIPP UQ Scenario:

1

Borehole accidentally drilled into the repository.

2

Radionuclides released into Culebra Dolomite and transported by groundwater.

3

Estimate travel time from release point to boundary of the region.

(to a good approximation flow is 2D)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 11 / 61

Model Problem: Uncertainty in Subsurface Flow

(e.g. risk analysis of radwaste disposal)

Darcy’s Law:

• q + k ∇p

=

• f

Incompressibility: ∇ · q = g +

Boundary Conditions

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 12 / 61

Model Problem: Uncertainty in Subsurface Flow

(e.g. risk analysis of radwaste disposal)

QUATERNARY MERCIA MUDSTONE VN-S CALDER FAULTED VN-S CALDER N-S CALDER FAULTED N-S CALDER DEEP CALDER FAULTED DEEP CALDER VN-S ST BEES FAULTED VN-S ST BEES N-S ST BEES FAULTED N-S ST BEES DEEP ST BEES FAULTED DEEP ST BEES BOTTOM NHM FAULTED BNHM SHALES + EVAP BROCKRAM FAULTED BROCKRAM COLLYHURST FAULTED COLLYHURST CARB LST FAULTED CARB LST N-S BVG FAULTED N-S BVG UNDIFF BVG FAULTED UNDIFF BVG F-H BVG FAULTED F-H BVG BLEAWATH BVG FAULTED BLEAWATH BVG TOP M-F BVG FAULTED TOP M-F BVG N-S LATTERBARROW DEEP LATTERBARROW N-S SKIDDAW DEEP SKIDDAW GRANITE FAULTED GRANITE WASTE VAULTS CROWN SPACE EDZ

Geology at Sellafield (former potential UK radwaste site) c NIREX UK Ltd.

Darcy’s Law:

• q + k(x, ω) ∇p

=

• f (x, ω)

Incompressibility: ∇ · q = g(x, ω) +

Boundary Conditions

Uncertainty in k = ⇒ Uncertainty in p &

• q

Stochastic Modelling!

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 12 / 61

Stochastic Modelling of Uncertainty

(simplified) Model uncertain conductivity tensor k as a random field k(x, ω)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 13 / 61

Stochastic Modelling of Uncertainty

(simplified) Model uncertain conductivity tensor k as a random field k(x, ω) k(x, ω) isotropic, scalar log k(x, ω) = Gaussian field

with isotropic covariance (e.g. Mat´ ern): R(x, y) := σ2ρν 1

λx − y

• Usually: σ2 > 1, λ < 1 & ν < 1

e.g. ρ1/2(r) = exp(−r) or ρ∞(r) = exp(−r2)

Compute statistics (eg. moments) of functionals of p and q, e.g.

◮ point evaluations, norms, averages ◮ breakthrough time (to boundary)

realisation (λ =

1 64, σ2 = 8)

contrast: max

x,y k(x) k(y) = O(109)

(Reasonably good fit to some field data [Gelhar, 1975], [Hoeksema et al, 1985])

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 13 / 61

A Case Study: Radioactive Waste Disposal

Example: Mat´ ern Family of Covariance Kernels

c(x, y) = cθ(r) = σ2 2ν−1 Γ(ν) 2√ν r λ ν Kν 2√ν r λ

• ,

r = x − y2 Kν : modified Bessel function of order ν Parameters θ = (σ2, λ, ν) σ2 : variance λ : correlation length ν : smoothness parameter

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 14 / 61

A Case Study: Radioactive Waste Disposal

Example: Mat´ ern Family of Covariance Kernels

c(x, y) = cθ(r) = σ2 2ν−1 Γ(ν) 2√ν r λ ν Kν 2√ν r λ

• ,

r = x − y2 Kν : modified Bessel function of order ν Parameters θ = (σ2, λ, ν) σ2 : variance λ : correlation length ν : smoothness parameter Special cases: ν = 1

2 :

c(r) = σ2 exp(− √ 2r/λ) exponential covariance ν = 1 : c(r) = σ2 2r

λ

• K1

2r

λ

• Bessel covariance

ν → ∞ : c(r) = σ2 exp(−r2/λ2) Gaussian covariance

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 14 / 61

A Case Study: Radioactive Waste Disposal

Mat´ ern Covariance Functions

1 2 3 4 5 0.2 0.4 0.6 0.8 1 r c(r) ν=1/2 ν=1 ν=2 ν=3 Gauss

σ2 = 1, λ = 1 σ2 = 8, λ = 1/64, ν = 1/2

Smoothness: Realisations Z(·, ω) ∈ C η(D) (H¨

• lder), for any η < ν.
• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 15 / 61

A Case Study: Radioactive Waste Disposal

Sampling from Z – Karhunen-Lo` eve expansion

Since c(x, y) is symmetric positive semidefinite and continuous, the covariance operator C : L2(D) → L2(D), (Cu)(x) =

• D

u(y)c(x, y) dy is selfadjoint, compact, nonnegative. Hence, the eigenvalues {µm}m∈N form a non-increasing sequence, accumulating at most at 0. Karhunen-Lo` eve expansion (converges in L2

P(Ω; L∞(D))):

Z(x, ω) = Z(x) +

∞

• m=1

√µm φm(x) Ym(ω) where {φm}m∈N are normalised eigenfunctions and Ym ∼ N(0, 1) i.i.d.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 16 / 61

A Case Study: Radioactive Waste Disposal

Sampling from Z – Karhunen-Lo` eve expansion

Since c(x, y) is symmetric positive semidefinite and continuous, the covariance operator C : L2(D) → L2(D), (Cu)(x) =

• D

u(y)c(x, y) dy is selfadjoint, compact, nonnegative. Hence, the eigenvalues {µm}m∈N form a non-increasing sequence, accumulating at most at 0. Truncated Karhunen-Lo` eve expansion (converges in L2

P(Ω; L∞(D))):

Zs(x, ω) = Z(x) +

s

• m=1

√µm φm(x) Ym(ω) where {φm}m∈N are normalised eigenfunctions and Ym ∼ N(0, 1) i.i.d. In practice, approximate Z by truncated KL expansion, using only s terms.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 16 / 61

A Case Study: Radioactive Waste Disposal

WIPP KL modes for Z = log k conditioned on 38 transmissivity observations

φm(x) unconditioned, m = 1, 2, 9, 16 φm(x) conditioned on data, m = 1, 2, 9, 16

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 17 / 61

Computational Challenges

Simulating PDEs with Highly Heterogeneous Random Coefficients

−∇ · (k(x, ω)∇p(x, ω)) = f (x, ω), x ∈ D ⊂ Rd, ω ∈ Ω (prob. space)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 18 / 61

Computational Challenges

Simulating PDEs with Highly Heterogeneous Random Coefficients

−∇ · (k(x, ω)∇p(x, ω)) = f (x, ω), x ∈ D ⊂ Rd, ω ∈ Ω (prob. space) Sampling from random field log k(x, ω) (correlated Gaussian):

◮ truncated Karhunen-Lo`

eve expansion of log k (see above)

◮ matrix factorisation, e.g. circulant embedding (FFT) ◮ via pseudodifferential “precision” operator (PDE solves)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 18 / 61

Computational Challenges

Simulating PDEs with Highly Heterogeneous Random Coefficients

−∇ · (k(x, ω)∇p(x, ω)) = f (x, ω), x ∈ D ⊂ Rd, ω ∈ Ω (prob. space) Sampling from random field log k(x, ω) (correlated Gaussian):

◮ truncated Karhunen-Lo`

eve expansion of log k (see above)

◮ matrix factorisation, e.g. circulant embedding (FFT) ◮ via pseudodifferential “precision” operator (PDE solves)

High-Dimensional Quadrature:

◮ Monte Carlo ◮ Multilevel Monte Carlo (RS), Quasi-Monte Carlo (DN) ◮ stochastic Galerkin/collocation (sparse grids, low-rank) (AN)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 18 / 61

Computational Challenges

Simulating PDEs with Highly Heterogeneous Random Coefficients

−∇ · (k(x, ω)∇p(x, ω)) = f (x, ω), x ∈ D ⊂ Rd, ω ∈ Ω (prob. space) Sampling from random field log k(x, ω) (correlated Gaussian):

◮ truncated Karhunen-Lo`

eve expansion of log k (see above)

◮ matrix factorisation, e.g. circulant embedding (FFT) ◮ via pseudodifferential “precision” operator (PDE solves)

High-Dimensional Quadrature:

◮ Monte Carlo ◮ Multilevel Monte Carlo (RS), Quasi-Monte Carlo (DN) ◮ stochastic Galerkin/collocation (sparse grids, low-rank) (AN)

Solve large number of multiscale deterministic PDEs:

◮ Efficient discretisation & FE error analysis (mesh size h) ◮ Multigrid Methods, AMG, DD Methods (robustness?)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 18 / 61

Computational Challenges

Why is it computationally so challenging?

Low regularity (global): k ∈ C η, η < ν < 1 = ⇒ fine mesh h ≪ 1 Large σ2 & exponential = ⇒ high contrast kmax/kmin > 106 Small λ = ⇒ multiscale + high stochast. dimension s > 100

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 19 / 61

Computational Challenges

Why is it computationally so challenging?

Low regularity (global): k ∈ C η, η < ν < 1 = ⇒ fine mesh h ≪ 1 Large σ2 & exponential = ⇒ high contrast kmax/kmin > 106 Small λ = ⇒ multiscale + high stochast. dimension s > 100

4 4.5 5 5.5 6 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

log t ECDF of log travel time 20 000 MC samples

M = 30 M = 100 M = 500 M = 1000

Source: Ernst et al, 2014 (s = M)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 19 / 61

Spatial Discretisation – Finite Elements

Write PDE (subject to p|∂D ≡ 0) in weak form: p(·, ω) ∈ H1

0(D) s.t.

• D

∇v · (k(x, ω)∇p(x, ω)) dx =

• D

f (x, ω)v dx, ∀v ∈ H1

0(D).

∃!p(·, ω) ∈ H1

0(D) a.s. in ω ∈ Ω (Lax-Milgram pathwise).

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 20 / 61

Spatial Discretisation – Finite Elements

Write PDE (subject to p|∂D ≡ 0) in weak form: p(·, ω) ∈ H1

0(D) s.t.

• D

∇v · (k(x, ω)∇p(x, ω)) dx =

• D

f (x, ω)v dx, ∀v ∈ H1

0(D).

∃!p(·, ω) ∈ H1

0(D) a.s. in ω ∈ Ω (Lax-Milgram pathwise).

Let Vh ⊂ H1

0(D) be the space of continuous, piecewise linear FEs

w.r.t. a mesh Th with mesh width h > 0. Find ph(·, ω) ∈ Vh that satisfies weak form for all vh ∈ Vh.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 20 / 61

Spatial Discretisation – Finite Elements

Write PDE (subject to p|∂D ≡ 0) in weak form: p(·, ω) ∈ H1

0(D) s.t.

• D

∇v · (k(x, ω)∇p(x, ω)) dx =

• D

f (x, ω)v dx, ∀v ∈ H1

0(D).

∃!p(·, ω) ∈ H1

0(D) a.s. in ω ∈ Ω (Lax-Milgram pathwise).

Let Vh ⊂ H1

0(D) be the space of continuous, piecewise linear FEs

w.r.t. a mesh Th with mesh width h > 0. Find ph(·, ω) ∈ Vh that satisfies weak form for all vh ∈ Vh. Write ph(x, ω) := Mh

i=1 Piϕi(x). Then this is equivalent to the

random matrix system A(ω)P(ω) = F(ω) with Ai,j(ω) :=

• D

∇ϕj · (k(x, ω)∇ϕi) dx, Fi(ω) :=

• D

f (x, ω)ϕi dx

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 20 / 61

Computational Challenges

Standard Monte Carlo Quadrature

Y(ω) ∈ Rs

Model(h)

− → P(ω) ∈ RMh

Output

− → Qh,s(ω) ∈ R

random input state vector quantity of interest

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 21 / 61

Computational Challenges

Standard Monte Carlo Quadrature

Y(ω) ∈ Rs

Model(h)

− → P(ω) ∈ RMh

Output

− → Qh,s(ω) ∈ R

random input state vector quantity of interest

Real QoI Q(ω) inaccessible (exact PDE), but we can assume E[Qh,s]

h→0, s→∞

− → E[Q] and |E[Qh,s − Q]| = O(hα) +O(s−α′) Standard Monte Carlo estimator for E[Q]: More detail below! ˆ QMC := 1 N

N

• i=1

Q(i)

h,s

where {Q(i)

h,s}N i=1 are i.i.d. samples computed with Model(h)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 21 / 61

Computational Challenges

Standard Monte Carlo Quadrature

Y(ω) ∈ Rs

Model(h)

− → P(ω) ∈ RMh

Output

− → Qh,s(ω) ∈ R

random input state vector quantity of interest

Real QoI Q(ω) inaccessible (exact PDE), but we can assume E[Qh,s]

h→0, s→∞

− → E[Q] and |E[Qh,s − Q]| = O(hα) +O(s−α′) Standard Monte Carlo estimator for E[Q]: More detail below! ˆ QMC := 1 N

N

• i=1

Q(i)

h,s

where {Q(i)

h,s}N i=1 are i.i.d. samples computed with Model(h)

Cost per sample is O(h−γ) (optimal: γ = d)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 21 / 61

Computational Challenges

Standard Monte Carlo Quadrature

Convergence of plain vanilla MC (mean square error): E ˆ QMC − E[Q] 2

• =: MSE

= V[Qh,s] N

sampling error

+

• E[Qh,s − Q]

2

• model error (“bias”)

Typical: α = 1 ⇒ MSE = O(N−1) + O(h2) ≤ TOL2, and so h ∼ TOL and N ∼ TOL−2.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 22 / 61

Computational Challenges

Standard Monte Carlo Quadrature

Convergence of plain vanilla MC (mean square error): E ˆ QMC − E[Q] 2

• =: MSE

= V[Qh,s] N

sampling error

+

• E[Qh,s − Q]

2

• model error (“bias”)

Typical: α = 1 ⇒ MSE = O(N−1) + O(h2) ≤ TOL2, and so h ∼ TOL and N ∼ TOL−2. Using optimal PDE solver: Cost = O(Nh−d) = O(TOL−(d+2))

(e.g. for TOL = 10−3: h ∼ 10−3, N ∼ 106 and Cost = O(1012) in 2D!!)

Quickly becomes prohibitively expensive !

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 22 / 61

Computational Challenges

Numerical Experiment with standard Monte Carlo D = (0, 1)2, unconditioned KL expansion, Q = − k ∂p

∂x1 L1(D)

using mixed FEs and the AMG solver amg1r5 [Ruge, St¨ uben, 1992]

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 23 / 61

Computational Challenges

Numerical Experiment with standard Monte Carlo D = (0, 1)2, unconditioned KL expansion, Q = − k ∂p

∂x1 L1(D)

using mixed FEs and the AMG solver amg1r5 [Ruge, St¨ uben, 1992]

Numerically observed FE-error: ≈ O(h3/4) = ⇒ α ≈ 3/4. Numerically observed cost/sample: ≈ O(h−2) = ⇒ γ ≈ 1.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 23 / 61

Computational Challenges

Numerical Experiment with standard Monte Carlo D = (0, 1)2, unconditioned KL expansion, Q = − k ∂p

∂x1 L1(D)

using mixed FEs and the AMG solver amg1r5 [Ruge, St¨ uben, 1992]

Numerically observed FE-error: ≈ O(h3/4) = ⇒ α ≈ 3/4. Numerically observed cost/sample: ≈ O(h−2) = ⇒ γ ≈ 1. Total cost to get RMSE O(TOL): ≈ O(TOL−14/3)

to get error reduction by a factor 2 → cost grows by a factor 25!

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 23 / 61

Computational Challenges

Numerical Experiment with standard Monte Carlo D = (0, 1)2, unconditioned KL expansion, Q = − k ∂p

∂x1 L1(D)

using mixed FEs and the AMG solver amg1r5 [Ruge, St¨ uben, 1992]

Numerically observed FE-error: ≈ O(h3/4) = ⇒ α ≈ 3/4. Numerically observed cost/sample: ≈ O(h−2) = ⇒ γ ≈ 1. Total cost to get RMSE O(TOL): ≈ O(TOL−14/3)

to get error reduction by a factor 2 → cost grows by a factor 25!

Case 1: σ2 = 1, λ = 0.3, ν = 0.5

TOL h−1 N Cost 0.01 129 1.4 × 104 21 min 0.002 1025 3.5 × 105 30 days

Case 2: σ2 = 3, λ = 0.1, ν = 0.5

TOL h−1 N Cost 0.01 513 8.5 × 103 4 h 0.002 Prohibitively large!! (actual numbers & CPU times on a cluster of 2GHz Intel T7300 processors)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 23 / 61

Computational Challenges

Alternatives – The Curse of Dimensionality

Stochastic Galerkin/collocation methods cost grows very fast with dimension s and polynomial order q

(faster than exponential) → #stoch. DOF: NSC = O((s + q)!/(s!q!))

lower with sparse grids (Smolyak) but growth still exponential in s! Anisotropic sparse grids or adaptive best N-term approximation can be dimension independent with sufficient smoothness (ν ≫ 1)! But best N-term theory by does not extend to nonlinear problems. Typically requires rewriting of PDE solvers and no access to good preconditioners (“intrusive”).

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 24 / 61

Computational Challenges

Alternatives – The Curse of Dimensionality

Stochastic Galerkin/collocation methods A Nouy cost grows very fast with dimension s and polynomial order q

(faster than exponential) → #stoch. DOF: NSC = O((s + q)!/(s!q!))

lower with sparse grids (Smolyak) but growth still exponential in s! Anisotropic sparse grids or adaptive best N-term approximation can be dimension independent with sufficient smoothness (ν ≫ 1)! But best N-term theory by does not extend to nonlinear problems. Typically requires rewriting of PDE solvers and no access to good preconditioners (“intrusive”).

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 24 / 61

Computational Challenges

Alternatives – The Curse of Dimensionality

Stochastic Galerkin/collocation methods A Nouy cost grows very fast with dimension s and polynomial order q

(faster than exponential) → #stoch. DOF: NSC = O((s + q)!/(s!q!))

lower with sparse grids (Smolyak) but growth still exponential in s! Anisotropic sparse grids or adaptive best N-term approximation can be dimension independent with sufficient smoothness (ν ≫ 1)! But best N-term theory by does not extend to nonlinear problems. Typically requires rewriting of PDE solvers and no access to good preconditioners (“intrusive”). Monte Carlo methods do not suffer from Curse of Dimensionality, they are “non-intrusive” and nonlinear parameter dependence is no problem, but plain vanilla version is too slow! Alternatives?

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 24 / 61

Monte Carlo Methods

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 25 / 61

The Buffon Needle Problem

In 1777, George Louis Leclerc, Comte de Buffon (1707–1788), French naturalist and mathematician, posed the following problem: Let a needle of length ℓ be thrown at random onto a horizontal plane ruled with parallel straight lines spaced by a distance d > ℓ from each other. What is the probability p that the needle will intersect one of these lines?

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 26 / 61

The Buffon Needle Problem

In 1777, George Louis Leclerc, Comte de Buffon (1707–1788), French naturalist and mathematician, posed the following problem: Let a needle of length ℓ be thrown at random onto a horizontal plane ruled with parallel straight lines spaced by a distance d > ℓ from each other. What is the probability p that the needle will intersect one of these lines? Answer: p = 2ℓ

πd

(simple geometric arguments)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 26 / 61

The Buffon Needle Problem

In 1777, George Louis Leclerc, Comte de Buffon (1707–1788), French naturalist and mathematician, posed the following problem: Let a needle of length ℓ be thrown at random onto a horizontal plane ruled with parallel straight lines spaced by a distance d > ℓ from each other. What is the probability p that the needle will intersect one of these lines? Answer: p = 2ℓ

πd

(simple geometric arguments)

Laplace later used similar randomised experiment to approximate π.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 26 / 61

The Buffon Needle Problem

In 1777, George Louis Leclerc, Comte de Buffon (1707–1788), French naturalist and mathematician, posed the following problem: Let a needle of length ℓ be thrown at random onto a horizontal plane ruled with parallel straight lines spaced by a distance d > ℓ from each other. What is the probability p that the needle will intersect one of these lines? Answer: p = 2ℓ

πd

(simple geometric arguments)

Laplace later used similar randomised experiment to approximate π. The term “Monte Carlo method” was coined by Metropolis, Ulam, von Neumann in the Manhattan project (Los Alamos, 1946).

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 26 / 61

The Buffon Needle Problem

Proceedings of the Royal Society of London, 2000

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 27 / 61

Basic Monte Carlo Simulation

Given a sequence {Xk} of i.i.d. copies of a given random variable X, basic MC simulation uses the estimator E [X] ≈ SN N , SN = X1 + · · · + XN. By the Strong Law of Large Numbers: SN N → E [X] a.s.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 28 / 61

Basic Monte Carlo Simulation

Given a sequence {Xk} of i.i.d. copies of a given random variable X, basic MC simulation uses the estimator E [X] ≈ SN N , SN = X1 + · · · + XN. By the Strong Law of Large Numbers: SN N → E [X] a.s. Also, for any measurable function f : 1 N

N

• k=1

f (Xk) → E [f (X)] a.s.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 28 / 61

Basic Monte Carlo Simulation

Given a sequence {Xk} of i.i.d. copies of a given random variable X, basic MC simulation uses the estimator E [X] ≈ SN N , SN = X1 + · · · + XN. By the Strong Law of Large Numbers: SN N → E [X] a.s. Also, for any measurable function f : 1 N

N

• k=1

f (Xk) → E [f (X)] a.s. If E [X] = µ and Var[X] = σ2, then (via the Central Limit Theorem) E [SN] = Nµ, Var[SN] = Nσ2 and S∗

N = SN − Nµ

√ Nσ → N(0, 1) i.e. the estimate is unbiased, its variance is σ2/N and the distribution of the normalised RV S∗

N becomes Gaussian as N → ∞.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 28 / 61

Monte Carlo Convergence Statements

1

Since E SN N − µ 2 = Var SN N = σ2 N → 0, we have mean square convergence of SN/N to µ.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 29 / 61

Monte Carlo Convergence Statements

1

Since E SN N − µ 2 = Var SN N = σ2 N → 0, we have mean square convergence of SN/N to µ.

2

Also Chebyshev’s Inequality implies, for any ǫ > 0, P

• SN

N − µ

• > N−1/2+ǫ
• ≤ σ2

N2ǫ ,

i.e. for any ǫ > 0, the probability of the error being larger than N−1/2+ǫ converges to zero as N → ∞.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 29 / 61

Monte Carlo Convergence Statements

1

Since E SN N − µ 2 = Var SN N = σ2 N → 0, we have mean square convergence of SN/N to µ.

2

Also Chebyshev’s Inequality implies, for any ǫ > 0, P

• SN

N − µ

• > N−1/2+ǫ
• ≤ σ2

N2ǫ ,

i.e. for any ǫ > 0, the probability of the error being larger than N−1/2+ǫ converges to zero as N → ∞.

3

If ρ := E

• |X − µ|3

< ∞, then the Berry-Esseen Inequality gives

• P{S∗

N ≤ x} − Φ(x)

• ≤

ρ 2σ3√ N , where Φ denotes cumulative density function (CDF) of N(0, 1).

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 29 / 61

Quasi-Monte Carlo Methods

In quasi-Monte Carlo (QMC) methods, samples are not chosen randomly, but special (deterministic) number sequences, known as low-discrepancy sequences, are used instead.

(discrepancy = measure of equidistribution)

−

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 30 / 61

Quasi-Monte Carlo Methods

In quasi-Monte Carlo (QMC) methods, samples are not chosen randomly, but special (deterministic) number sequences, known as low-discrepancy sequences, are used instead.

(discrepancy = measure of equidistribution)

Example: The van der Corput sequence is a low-discrepancy sequence for the unit interval. The first 11 numbers are {xn}11

n=1 =

• 0, 1

3, 2 3, 1 9, 4 9, 7 9, 2 9, 5 9, 8 9, 1 27, 10 27

• .

(for base 3, it is given by xn = k

3j , where j increases monotonically and k runs

through all nonnegative integers such that k/3j is an irreducible fraction < 1. The ordering in k is obtained by representing k in base 3 and reversing the digits)

0.2 0.4 0.6 0.8 1 −0.1 0.1

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 30 / 61

Quasi-Monte Carlo Methods

In quasi-Monte Carlo (QMC) methods, samples are not chosen randomly, but special (deterministic) number sequences, known as low-discrepancy sequences, are used instead.

(discrepancy = measure of equidistribution)

Example: The van der Corput sequence is a low-discrepancy sequence for the unit interval. The first 11 numbers are {xn}11

n=1 =

• 0, 1

3, 2 3, 1 9, 4 9, 7 9, 2 9, 5 9, 8 9, 1 27, 10 27

• .

(for base 3, it is given by xn = k

3j , where j increases monotonically and k runs

through all nonnegative integers such that k/3j is an irreducible fraction < 1. The ordering in k is obtained by representing k in base 3 and reversing the digits)

0.2 0.4 0.6 0.8 1 −0.1 0.1

Convergence rate improved from O(N−1/2) to O(N−2)!

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 30 / 61

Quasi-Monte Carlo Methods

In quasi-Monte Carlo (QMC) methods, samples are not chosen randomly, but special (deterministic) number sequences, known as low-discrepancy sequences, are used instead.

(discrepancy = measure of equidistribution)

Example: The van der Corput sequence is a low-discrepancy sequence for the unit interval. The first 11 numbers are {xn}11

n=1 =

• 0, 1

3, 2 3, 1 9, 4 9, 7 9, 2 9, 5 9, 8 9, 1 27, 10 27

• .

(for base 3, it is given by xn = k

3j , where j increases monotonically and k runs

through all nonnegative integers such that k/3j is an irreducible fraction < 1. The ordering in k is obtained by representing k in base 3 and reversing the digits)

0.2 0.4 0.6 0.8 1 −0.1 0.1

Convergence rate improved from O(N−1/2) to O(N−2)! D Nuyens

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 30 / 61

Variance Reduction: Antithetic Sampling

The constant in the MC convergence rate is the variance σ2 of the RV X. By designing an equivalent MC approximation with lower variance, we can expect faster convergence.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 31 / 61

Variance Reduction: Antithetic Sampling

The constant in the MC convergence rate is the variance σ2 of the RV X. By designing an equivalent MC approximation with lower variance, we can expect faster convergence. To approximate E [X] by standard MC, we draw independent samples {Xk}N

k=1 of X and compute the sample average SN/N.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 31 / 61

Variance Reduction: Antithetic Sampling

The constant in the MC convergence rate is the variance σ2 of the RV X. By designing an equivalent MC approximation with lower variance, we can expect faster convergence. To approximate E [X] by standard MC, we draw independent samples {Xk}N

k=1 of X and compute the sample average SN/N.

Now, let { ˜ Xk}N

k=1 be a second set of samples of X with sample

average ˜ SN/N.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 31 / 61

Variance Reduction: Antithetic Sampling

The constant in the MC convergence rate is the variance σ2 of the RV X. By designing an equivalent MC approximation with lower variance, we can expect faster convergence. To approximate E [X] by standard MC, we draw independent samples {Xk}N

k=1 of X and compute the sample average SN/N.

Now, let { ˜ Xk}N

k=1 be a second set of samples of X with sample

average ˜ SN/N. Since both averages converge to E [X], so does 1

2(SN/N + ˜

SN/N). When Xk and ˜ Xk are negatively correlated they are called antithetic samples, and the approximation

1 2N (SN + ˜

SN) is a more reliable approximation of E [X] than

1 2N S2N.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 31 / 61

Variance Reduction: Antithetic Sampling

Theorem

Let the two sequences of RVs {Xk} and { ˜ Xk} be identically distributed with Cov(Xj, Xk) = Cov( ˜ Xj, ˜ Xk) = 0 for j = k. Then Var

• SN + ˜

SN 2N

• = Var

S2N 2N

• + 1

2 Cov

• SN

N , ˜ SN N

• ≤ Var

SN N

• .
• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 32 / 61

Variance Reduction: Antithetic Sampling

Theorem

Let the two sequences of RVs {Xk} and { ˜ Xk} be identically distributed with Cov(Xj, Xk) = Cov( ˜ Xj, ˜ Xk) = 0 for j = k. Then Var

• SN + ˜

SN 2N

• = Var

S2N 2N

• + 1

2 Cov

• SN

N , ˜ SN N

• ≤ Var

SN N

• .

Worst case: Variance of average of N samples and N antithetic samples no better than variance of N independent samples.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 32 / 61

Variance Reduction: Antithetic Sampling

Theorem

Let the two sequences of RVs {Xk} and { ˜ Xk} be identically distributed with Cov(Xj, Xk) = Cov( ˜ Xj, ˜ Xk) = 0 for j = k. Then Var

• SN + ˜

SN 2N

• = Var

S2N 2N

• + 1

2 Cov

• SN

N , ˜ SN N

• ≤ Var

SN N

• .

Worst case: Variance of average of N samples and N antithetic samples no better than variance of N independent samples. Best case: SN/N and ˜ SN/N negatively correlated; therefore variance of N samples and N antithetic samples is smaller than variance of 2N indepependent samples.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 32 / 61

Easy Example: Predator-Prey Dynamical System

Consider the popular model of the dynamics of two interacting populations ˙ u = ˙ u1 ˙ u2

• =

u1(1 − u2) u2(u1 − 1)

• ,

u(0) = u0. Let the initial condition u0 be uncertain and modeled as a (uniform) random vector u0 ∼ U(Γ), where Γ denotes the square Γ = u0 + [−ǫ, ǫ]2.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 33 / 61

Easy Example: Predator-Prey Dynamical System

Consider the popular model of the dynamics of two interacting populations ˙ u = ˙ u1 ˙ u2

• =

u1(1 − u2) u2(u1 − 1)

• ,

u(0) = u0. Let the initial condition u0 be uncertain and modeled as a (uniform) random vector u0 ∼ U(Γ), where Γ denotes the square Γ = u0 + [−ǫ, ǫ]2. Goal: estimate E [u1(T)] at time T > 0. Denote by uh = uh(ω) the explicit Euler approximation after M time steps of length h = T

M starting with initial data u0 = u0(ω).

Define the QoI Q = φ(u(T)) = u1(T) for u = [u1, u2]T and estimate E [Qh] using the MC method just described, where Qh = φ(uh).

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 33 / 61

Easy Example: Predator-Prey Dynamical System

Consider the popular model of the dynamics of two interacting populations ˙ u = ˙ u1 ˙ u2

• =

u1(1 − u2) u2(u1 − 1)

• ,

u(0) = u0. Let the initial condition u0 be uncertain and modeled as a (uniform) random vector u0 ∼ U(Γ), where Γ denotes the square Γ = u0 + [−ǫ, ǫ]2. Goal: estimate E [u1(T)] at time T > 0. Denote by uh = uh(ω) the explicit Euler approximation after M time steps of length h = T

M starting with initial data u0 = u0(ω).

Define the QoI Q = φ(u(T)) = u1(T) for u = [u1, u2]T and estimate E [Qh] using the MC method just described, where Qh = φ(uh). Can easily make this problem high dimensional by adding uncertain coefficients (= environment) !

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 33 / 61

Easy Example: Predator-Prey Dynamical System

Analysis of Monte Carlo Estimator

Denote the Monte Carlo estimator for E [Qh] by

• Qh :=

Qh,N = 1 N

N

• k=1

Q(k)

h

(i.e the average over N samples {Q(k)

h }N k=1 of Qh with step size h)

Expect better approximations for N large and h small.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 34 / 61

Easy Example: Predator-Prey Dynamical System

Analysis of Monte Carlo Estimator

Denote the Monte Carlo estimator for E [Qh] by

• Qh :=

Qh,N = 1 N

N

• k=1

Q(k)

h

(i.e the average over N samples {Q(k)

h }N k=1 of Qh with step size h)

Expect better approximations for N large and h small. Error with N samples and M = T/h time steps: eh,N = |E [Q] − Qh| ≤ |E [Q] − E [Qh] |

• discretisation error

+ |E [Qh] − Qh|

• Monte Carlo error
• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 34 / 61

Easy Example: Predator-Prey Dynamical System

Analysis of Monte Carlo Estimator

Denote the Monte Carlo estimator for E [Qh] by

• Qh :=

Qh,N = 1 N

N

• k=1

Q(k)

h

(i.e the average over N samples {Q(k)

h }N k=1 of Qh with step size h)

Expect better approximations for N large and h small. Error with N samples and M = T/h time steps: eh,N = |E [Q] − Qh| ≤ |E [Q] − E [Qh] |

• discretisation error

+ |E [Qh] − Qh|

• Monte Carlo error

Can show also that the mean square error satisfies (with equality!) E

• E [Q] −

Qh 2 =

• E [Q − Qh]

2 + Var[Qh] N

Problem Sheet!

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 34 / 61

Easy Example: Predator-Prey Dynamical System

Discretisation Error (Bias)

Explicit Euler discretisation error (with some constant K > 0): u(T) − uh ≤ Kh.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 35 / 61

Easy Example: Predator-Prey Dynamical System

Discretisation Error (Bias)

Explicit Euler discretisation error (with some constant K > 0): u(T) − uh ≤ Kh. φ Lipschitz-continuous with constant L = 1: |φ(u(T)) − φ(uh)| ≤ KLh.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 35 / 61

Easy Example: Predator-Prey Dynamical System

Discretisation Error (Bias)

Explicit Euler discretisation error (with some constant K > 0): u(T) − uh ≤ Kh. φ Lipschitz-continuous with constant L = 1: |φ(u(T)) − φ(uh)| ≤ KLh. Therefore |E [Q] − E [Qh] | = |E [Q − Qh] | ≤ KLh.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 35 / 61

Easy Example: Predator-Prey Dynamical System

Balancing discretisation and MC error

Using the Berry-Esseen bound with Var[Qh] = σ2

h we get

(see the Problem Sheet)

P

• E [Qh] −

Qh,N

• ≤ 1.96σh

√ N

• > 0.95 + O(N−1/2)
• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 36 / 61

Easy Example: Predator-Prey Dynamical System

Balancing discretisation and MC error

Using the Berry-Esseen bound with Var[Qh] = σ2

h we get

(see the Problem Sheet)

P

• E [Qh] −

Qh,N

• ≤ 1.96σh

√ N

• > 0.95 + O(N−1/2)

Combined with the discretisation error (using triangle inequality): P

• eh,N ≤ 2KLh + 1.96σh

√ N

• > 0.95 + O(N−1/2).
• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 36 / 61

Easy Example: Predator-Prey Dynamical System

Balancing discretisation and MC error

Using the Berry-Esseen bound with Var[Qh] = σ2

h we get

(see the Problem Sheet)

P

• E [Qh] −

Qh,N

• ≤ 1.96σh

√ N

• > 0.95 + O(N−1/2)

Combined with the discretisation error (using triangle inequality): P

• eh,N ≤ 2KLh + 1.96σh

√ N

• > 0.95 + O(N−1/2).

Balancing discretization and MC errors, i.e. KLh ≈ TOL 2 and 1.96σh √ N ≈ TOL 2 , leads to h ≈ TOL 2KL , N ≈ 16σ2

h

TOL2 and Cost = O(MN) = O(TOL−3)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 36 / 61

Easy Example: Predator-Prey Dynamical System

Sample trajectories

0.5 1 1.5 2 2.5 3 0.5 1 1.5 2 2.5 3 u1 u2

Predator-prey system integrated from 0 to T = 6 with u0 = [0.5, 2]T and ǫ = 0.2. Unperturbed trajectory (black) and 15 perturbed trajectories.

(for the unperturbed trajectory u1(T) = 1.3942)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 37 / 61

Easy Example: Predator-Prey Dynamical System

Antithetic sampling

We can use antithetic sampling for this problem by noting that, if u0 ∼ U(Γ), then the same holds for the random vector ˜ u0 := 2u0 − u0. Thus, the trajectories generated by the random initial data ˜ u0 have the same distribution as those generated by u0.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 38 / 61

Easy Example: Predator-Prey Dynamical System

Antithetic sampling

We can use antithetic sampling for this problem by noting that, if u0 ∼ U(Γ), then the same holds for the random vector ˜ u0 := 2u0 − u0. Thus, the trajectories generated by the random initial data ˜ u0 have the same distribution as those generated by u0. Let Qh = φ(uh) be the basic samples and ˜ Qh = φ(˜ uh) the antithetic

• counterparts. Note that all pairs of samples are independent except

each sample and its antithetic counterpart.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 38 / 61

Easy Example: Predator-Prey Dynamical System

Antithetic sampling

We can use antithetic sampling for this problem by noting that, if u0 ∼ U(Γ), then the same holds for the random vector ˜ u0 := 2u0 − u0. Thus, the trajectories generated by the random initial data ˜ u0 have the same distribution as those generated by u0. Let Qh = φ(uh) be the basic samples and ˜ Qh = φ(˜ uh) the antithetic

• counterparts. Note that all pairs of samples are independent except

each sample and its antithetic counterpart. Use

1 2(

Qh,N + ˜ Qh,N) instead of

• Qh,2N (same cost).
• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 38 / 61

Easy Example: Predator-Prey Dynamical System

Numerical Experiment – Comparing standard and antithetic sampling

0.5 1 1.5 2 2.5 3 3.5 4 x 10

4

1.46 1.47 1.48 1.49 1.5 1.51 1.52 N u1(tend) 0.5 1 1.5 2 2.5 3 3.5 4 x 10

4

1.46 1.47 1.48 1.49 1.5 1.51 1.52 N u1(tend)

Estimates of E [u1(T)] using standard MC with N samples (left) vs. antithetic sampling with N/2 pairs of samples (right) (with 95% confidence intervals)

Var[Qh] and Cov(Qh, ˜ Qh) estimated using sample variance and covariance (resp.), i.e.

1 N−1

N

k=1(Q(k) h

− Qh,N)2 and

1 N−1

N

k=1(Q(k) h

− Qh,N)( ˜ Q(k)

h

− ˜ Qh,N)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 39 / 61

Multilevel Monte Carlo Methods – History

The multilevel Monte Carlo method is a powerful new variance reduction technique. First ideas for high-dimensional quadrature by Heinrich, 2000. Independently discovered and popularised by Giles, 2007 in the context of stochastic DEs in mathematical finance.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 40 / 61

Multilevel Monte Carlo Methods – History

The multilevel Monte Carlo method is a powerful new variance reduction technique. First ideas for high-dimensional quadrature by Heinrich, 2000. Independently discovered and popularised by Giles, 2007 in the context of stochastic DEs in mathematical finance. First papers in the context of UQ:

◮ [Barth, Schwab, Zollinger ’11] ◮ [Cliffe, Giles, RS, Teckentrup ’11]

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 40 / 61

Multilevel Monte Carlo Methods – History

The multilevel Monte Carlo method is a powerful new variance reduction technique. First ideas for high-dimensional quadrature by Heinrich, 2000. Independently discovered and popularised by Giles, 2007 in the context of stochastic DEs in mathematical finance. First papers in the context of UQ:

◮ [Barth, Schwab, Zollinger ’11] ◮ [Cliffe, Giles, RS, Teckentrup ’11]

Stochastic simulation of discrete state systems (biology, chemistry):

◮ [Anderson, Higham ’12], . . .

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 40 / 61

Multilevel Monte Carlo Methods – History

The multilevel Monte Carlo method is a powerful new variance reduction technique. First ideas for high-dimensional quadrature by Heinrich, 2000. Independently discovered and popularised by Giles, 2007 in the context of stochastic DEs in mathematical finance. First papers in the context of UQ:

◮ [Barth, Schwab, Zollinger ’11] ◮ [Cliffe, Giles, RS, Teckentrup ’11]

Stochastic simulation of discrete state systems (biology, chemistry):

◮ [Anderson, Higham ’12], . . .

In MCMC for Bayesian inference:

◮ [Hoang, Schwab, Stuart ’13] ◮ [Dodwell, Ketelsen, RS, Teckentrup ’15]

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 40 / 61

Standard Monte Carlo – Mean-Square Error

To estimate the expectation E [Q] of a quantity of interest Q, assume only approximations Qh ≈ Q are computable, where h > 0 denotes a discretization parameter (time step size, grid size, . . . ) and lim

h→0 E [Qh] = E [Q] .

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 41 / 61

Standard Monte Carlo – Mean-Square Error

To estimate the expectation E [Q] of a quantity of interest Q, assume only approximations Qh ≈ Q are computable, where h > 0 denotes a discretization parameter (time step size, grid size, . . . ) and lim

h→0 E [Qh] = E [Q] .

More precisely, we assume error to converge in mean with rate α: |E [Qh − Q] | hα, as h → 0, α > 0.

(in the predator-prey case α = 1)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 41 / 61

Standard Monte Carlo – Mean-Square Error

To estimate the expectation E [Q] of a quantity of interest Q, assume only approximations Qh ≈ Q are computable, where h > 0 denotes a discretization parameter (time step size, grid size, . . . ) and lim

h→0 E [Qh] = E [Q] .

More precisely, we assume error to converge in mean with rate α: |E [Qh − Q] | hα, as h → 0, α > 0.

(in the predator-prey case α = 1)

From above (and Problem Sheet) the mean square error (MSE) is E

• Qh,N − E [Q]

2 = Var[Qh] N +

• E [Qh − Q]

2.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 41 / 61

Standard Monte Carlo – Cost Scaling

Denote by C (Q(k)

h ) cost associated with computing one sample Q(k) h

(e.g. in terms of the number of floating-point operations required)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 42 / 61

Standard Monte Carlo – Cost Scaling

Denote by C (Q(k)

h ) cost associated with computing one sample Q(k) h

(e.g. in terms of the number of floating-point operations required)

Cost typcially grows linearly or with power γ ≥ 1 with M = T/h. We assume C (Q(k)

h ) h−γ,

γ ≥ 1. so that C ( Qh,N) Nh−γ (in the predator-prey case γ = 1).

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 42 / 61

Standard Monte Carlo – Cost Scaling

Denote by C (Q(k)

h ) cost associated with computing one sample Q(k) h

(e.g. in terms of the number of floating-point operations required)

Cost typcially grows linearly or with power γ ≥ 1 with M = T/h. We assume C (Q(k)

h ) h−γ,

γ ≥ 1. so that C ( Qh,N) Nh−γ (in the predator-prey case γ = 1). To balance the two MSE components, bound each of them by TOL2

2

, resulting in a total bound of TOL2 for the MSE.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 42 / 61

Standard Monte Carlo – Cost Scaling

Denote by C (Q(k)

h ) cost associated with computing one sample Q(k) h

(e.g. in terms of the number of floating-point operations required)

Cost typcially grows linearly or with power γ ≥ 1 with M = T/h. We assume C (Q(k)

h ) h−γ,

γ ≥ 1. so that C ( Qh,N) Nh−γ (in the predator-prey case γ = 1). To balance the two MSE components, bound each of them by TOL2

2

, resulting in a total bound of TOL2 for the MSE. This yields (since Qh → Q, we have Var[Qh] ≈ Var[Q] = constant) N ≥ 2 Var[Qh]TOL−2 and h TOL1/α.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 42 / 61

Standard Monte Carlo – Cost Scaling

Denote by C (Q(k)

h ) cost associated with computing one sample Q(k) h

(e.g. in terms of the number of floating-point operations required)

Cost typcially grows linearly or with power γ ≥ 1 with M = T/h. We assume C (Q(k)

h ) h−γ,

γ ≥ 1. so that C ( Qh,N) Nh−γ (in the predator-prey case γ = 1). To balance the two MSE components, bound each of them by TOL2

2

, resulting in a total bound of TOL2 for the MSE. This yields (since Qh → Q, we have Var[Qh] ≈ Var[Q] = constant) N ≥ 2 Var[Qh]TOL−2 and h TOL1/α. So total cost to have MSE < TOL2 using standard MC estimator is C ( Qh,N) TOL−2−γ/α

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 42 / 61

Multilevel Monte Carlo Methods

Main Idea

Key idea: use realisations of Qh on a hierarchy of different levels, i.e., for a sequence h0, . . . , hL of discretization parameters, and decompose E [QhL] = E [Qh0] +

L

• ℓ=1

E

• Qhℓ − Qhℓ−1
• =:

L

• ℓ=0

E [Yℓ] , where h0 > 0, hℓ = hℓ−1/s, for ℓ = 1, . . . , L, and s ∈ N \ {1}.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 43 / 61

Multilevel Monte Carlo Methods

Main Idea

Key idea: use realisations of Qh on a hierarchy of different levels, i.e., for a sequence h0, . . . , hL of discretization parameters, and decompose E [QhL] = E [Qh0] +

L

• ℓ=1

E

• Qhℓ − Qhℓ−1
• =:

L

• ℓ=0

E [Yℓ] , where h0 > 0, hℓ = hℓ−1/s, for ℓ = 1, . . . , L, and s ∈ N \ {1}. Given (unbiased) estimators { Yℓ}L

ℓ=0 for E [Yℓ], we refer to

• QML

L

:=

L

• ℓ=0
• Yℓ

as a multilevel estimator for Q

(today use standard MC on all levels).

All expectations E [Yℓ] sampled indep. ⇒ Var QML

L

= L

ℓ=0 Var

Yℓ.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 43 / 61

Multilevel Monte Carlo Methods

Multilevel Monte Carlo estimator

If each Yℓ is itself a standard Monte Carlo estimator, i.e.,

• Y0 =

Y0,N0 := 1 N0

N0

• k=0

Q(k)

h0

and

• Yℓ =

Yℓ,Nℓ := 1 Nℓ

Nℓ

• k=0
• Q(k)

hℓ − Q(k) hℓ−1

• ,

ℓ = 1, . . . , L,

• ne obtains a multilevel Monte Carlo estimator.

The associated MSE then has the standard decomposition E

• QML

L,{Nℓ} − E [Q]

2 =

L

• ℓ=0

Var Yℓ Nℓ + E [QML − Q]2 into sample variance and bias (shown as for standard MC above).

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 44 / 61

Multilevel Monte Carlo Methods

MLMC variance reduction

Choose discretisation parameters and numbers of samples again to balance the terms in the MSE. The bias term is the same as for the standard MC estimator, leading again to a choice of hL = h TOL1/α.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 45 / 61

Multilevel Monte Carlo Methods

MLMC variance reduction

Choose discretisation parameters and numbers of samples again to balance the terms in the MSE. The bias term is the same as for the standard MC estimator, leading again to a choice of hL = h TOL1/α. But why do we get variance reduction – or rather lower cost for the same variance?

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 45 / 61

Multilevel Monte Carlo Methods

MLMC variance reduction

Choose discretisation parameters and numbers of samples again to balance the terms in the MSE. The bias term is the same as for the standard MC estimator, leading again to a choice of hL = h TOL1/α. But why do we get variance reduction – or rather lower cost for the same variance? Two reasons:

1

As we coarsen the problem L → 0, the cost per sample decays rapidly from level to level, with O(sγ)

2

Since Qh → Q, then Var[Yℓ] = Var[Qhℓ − Qhℓ−1] → 0 as ℓ → ∞, allowing for smaller and smaller sample sizes Nℓ on finer and finer levels.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 45 / 61

Multilevel Monte Carlo Methods

Optimal sample sizes

The cost of the MLMC estimator is C ( QML

L,{Nℓ}) = L

• ℓ=0

NℓCℓ, Cℓ := C (Y (k)

ℓ

).

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 46 / 61

Multilevel Monte Carlo Methods

Optimal sample sizes

The cost of the MLMC estimator is C ( QML

L,{Nℓ}) = L

• ℓ=0

NℓCℓ, Cℓ := C (Y (k)

ℓ

). Treating the Nℓ as continuous variables, we can minimise the cost of the MLMC estimator for a fixed variance

L

• ℓ=0

Var Yℓ Nℓ = TOL2 2 using the Lagrange multiplier method.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 46 / 61

Multilevel Monte Carlo Methods

Optimal sample sizes

The cost of the MLMC estimator is C ( QML

L,{Nℓ}) = L

• ℓ=0

NℓCℓ, Cℓ := C (Y (k)

ℓ

). Treating the Nℓ as continuous variables, we can minimise the cost of the MLMC estimator for a fixed variance

L

• ℓ=0

Var Yℓ Nℓ = TOL2 2 using the Lagrange multiplier method. The solution to this constrained minimisation problem is Nℓ ≃

• Var[Yℓ]/Cℓ

with implied constant chosen such that the total variance is TOL2

2

(which leads to the constant

2 TOL2

• ℓ

√Cℓ Var Yℓ)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 46 / 61

Multilevel Monte Carlo Methods

MLMC cost

This results in a total cost on level ℓ proportional to √Cℓ Var Yℓ and therefore C ( QML

L,{Nℓ}) ≤

2 TOL2

• L
• ℓ=0
• Cℓ Var Yℓ

2

For comparison, the cost for standard MC is C ( QSL

hL,N) = 2 TOL2 CL Var[QhL].

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 47 / 61

Multilevel Monte Carlo Methods

MLMC cost

This results in a total cost on level ℓ proportional to √Cℓ Var Yℓ and therefore C ( QML

L,{Nℓ}) ≤

2 TOL2

• L
• ℓ=0
• Cℓ Var Yℓ

2

For comparison, the cost for standard MC is C ( QSL

hL,N) = 2 TOL2 CL Var[QhL].

If Var Yℓ decays faster than Cℓ increases, the cost on level ℓ = 0

• dominates. Since Var[Qh0] ≈ Var[QhL], the cost ratio of MLMC to

MC estimation is then approximately C0/CL s−Lγ

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 47 / 61

Multilevel Monte Carlo Methods

MLMC cost

This results in a total cost on level ℓ proportional to √Cℓ Var Yℓ and therefore C ( QML

L,{Nℓ}) ≤

2 TOL2

• L
• ℓ=0
• Cℓ Var Yℓ

2

For comparison, the cost for standard MC is C ( QSL

hL,N) = 2 TOL2 CL Var[QhL].

If Var Yℓ decays faster than Cℓ increases, the cost on level ℓ = 0

• dominates. Since Var[Qh0] ≈ Var[QhL], the cost ratio of MLMC to

MC estimation is then approximately C0/CL s−Lγ If Cℓ increases faster than Var Yℓ decays, then the cost on level ℓ = L dominates, and then the cost ratio is approximately Var[YL]/ Var[QhL] TOL2

(provided E

• (Q − QL)2

(E [Q − QL])2, which is problem dependent).

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 47 / 61

Multilevel Monte Carlo Methods

General Complexity Theorem

Theorem

Let TOL < exp(−1) and assume there are constants α, β, γ > 0 such that α ≥ min{β, γ} and, for all ℓ = 0, . . . , L, (M1) |E [Qhℓ] − E [Q] | hα

ℓ ,

(M2) Var[ Yℓ] N−1

ℓ

hβ

ℓ ,

(M3) C ( Yℓ) Nℓ h−γ

ℓ .

Then there are L and {Nℓ}L

ℓ=0 s.t. E

• QML

L,{Nℓ} − E [Q]

2 ≤ TOL2 and C ( QML

L,{Nℓ})

     TOL−2, if β > γ, TOL−2 | log TOL|2, if β = γ, TOL−2−(γ−β)/α, if β < γ.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 48 / 61

Multilevel Monte Carlo Methods

Adaptive MLMC Algorithm

The optimal values of L and N0, . . . , NL can be estimated adaptively using the sample averages ( Yℓ) and sample variances (s2

ℓ) of Yℓ.

The sample variances s2

ℓ can be used directly to estimate Var[

Yℓ].

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 49 / 61

Multilevel Monte Carlo Methods

Adaptive MLMC Algorithm

The optimal values of L and N0, . . . , NL can be estimated adaptively using the sample averages ( Yℓ) and sample variances (s2

ℓ) of Yℓ.

The sample variances s2

ℓ can be used directly to estimate Var[

Yℓ]. To bound the bias error, we assume there exists an h∗ > 0 such that for all h ≤ h∗ hα |E [Qh − Q] | hα .

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 49 / 61

Multilevel Monte Carlo Methods

Adaptive MLMC Algorithm

The optimal values of L and N0, . . . , NL can be estimated adaptively using the sample averages ( Yℓ) and sample variances (s2

ℓ) of Yℓ.

The sample variances s2

ℓ can be used directly to estimate Var[

Yℓ]. To bound the bias error, we assume there exists an h∗ > 0 such that for all h ≤ h∗ hα |E [Qh − Q] | hα . Then it follows (via the inverse triangle inequality) that |E [QhL − Q] | ≤ 1 sα − 1

• YL
• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 49 / 61

Multilevel Monte Carlo Methods

Adaptive MLMC Algorithm

The optimal values of L and N0, . . . , NL can be estimated adaptively using the sample averages ( Yℓ) and sample variances (s2

ℓ) of Yℓ.

The sample variances s2

ℓ can be used directly to estimate Var[

Yℓ]. To bound the bias error, we assume there exists an h∗ > 0 such that for all h ≤ h∗ hα |E [Qh − Q] | hα . Then it follows (via the inverse triangle inequality) that |E [QhL − Q] | ≤ 1 sα − 1

• YL

Thus, a computable bias error estimate for given hL > 0, which can be used to decide if hL is sufficiently small or L should be increased.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 49 / 61

Recall the PDE Model Problem

Case Study in Radioactive Waste Disposal

Recall

uncertain k

→ Darcy’s Law:

• q + k ∇p

= f Incompressibility: ∇ · q = →

uncertain p, q

Typical simplified model for k: log k(x, ω) = isotropic, scalar Gaussian

e.g. with exp. covariance R(x, y) = σ2 exp

• − x−y

λ

• Parametrise: log k(x, ω) ≈ J

j=1

√µjψj(x)Yj(ω)

Karhunen-Lo` eve expansion with Yj(ω) i.i.d. N(0, 1)

FE discretisation: A(ω) P(ω) = b(ω) QoI Q(ω) = φ

• p(ω)
• , e.g. particle travel time from leak to boundary
• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 50 / 61

Recall the PDE Model Problem

Numerical Experiment with standard Monte Carlo

Recall

D = (0, 1)2, unconditioned KL expansion, Q = − k ∂p

∂x1 L1(D)

using mixed FEs and AMG solver amg1r5 [Ruge, St¨ uben, 1992]

• Num. observed FE-error: ≈ O(h−3/4) ⇒ α ≈ 3/4
• Num. observed cost/sample: ≈ O(h−d) ⇒ γ ≈ 2

Total cost to get RMSE O(TOL): ≈ O(TOL−14/3)

to get error reduction by a factor 2 → cost grows by a factor 25!

Case 1: σ2 = 1, λ = 0.3, ν = 0.5

TOL h−1 N Cost 0.01 129 1.4 × 104 21 min 0.002 1025 3.5 × 105 30 days

Case 2: σ2 = 3, λ = 0.1, ν = 0.5

TOL h−1 N Cost 0.01 513 8.5 × 103 4 h 0.002 Prohibitively large!! (actual numbers & CPU times on a cluster of 2GHz Intel T7300 processors)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 51 / 61

Model Hierarchy (FE mesh & KL truncation)

L Vℓ Xℓ

Here α ≈ 1 (for smooth fctls.) and γ ≈ d (using AMG)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 52 / 61

Multilevel Monte Carlo for the PDE Model

Theoretically predicted gains

Assuming optimal AMG solver (i.e. γ ≈ d), α ≈ 1 and β ≈ 2α. Then the cost in Rd is d MC MLMC cost per sample 1 O(TOL−3) O(TOL−2) O(TOL−1) 2 O(TOL−4) O(TOL−2) O(TOL−2) 3 O(TOL−5) O(TOL−3) O(TOL−3)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 53 / 61

Multilevel Monte Carlo for the PDE Model

Theoretically predicted gains

Assuming optimal AMG solver (i.e. γ ≈ d), α ≈ 1 and β ≈ 2α. Then the cost in Rd is d MC MLMC cost per sample 1 O(TOL−3) O(TOL−2) O(TOL−1) 2 O(TOL−4) O(TOL−2) O(TOL−2) 3 O(TOL−5) O(TOL−3) O(TOL−3)

Optimality (for d = γ > β = 2α)

MLMC cost is asymptotically the same as one deterministic solve to accuracy TOL for d > 1, i.e. O(TOL−γ/α) !! (only for rough problems!)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 53 / 61

Multilevel Monte Carlo for the PDE Model

Theoretically predicted gains

Assuming optimal AMG solver (i.e. γ ≈ d), α ≈ 1 and β ≈ 2α. Then the cost in Rd is d MC MLMC cost per sample 1 O(TOL−3) O(TOL−2) O(TOL−1) 2 O(TOL−4) O(TOL−2) O(TOL−2) 3 O(TOL−5) O(TOL−3) O(TOL−3)

Optimality (for d = γ > β = 2α)

MLMC cost is asymptotically the same as one deterministic solve to accuracy TOL for d > 1, i.e. O(TOL−γ/α) !! (only for rough problems!) Can we achieve such huge gains in practice?

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 53 / 61

Multilevel Monte Carlo for the PDE Model

Numerical Experiments: D = (0, 1)2; Q = pL2(D); standard FEs

ν = 1

2, σ2 = 1, λ = 0.3, h0 = 1 8

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 54 / 61

Multilevel Monte Carlo for the PDE Model

Numerical Experiments: D = (0, 1)2; Q = pL2(D); standard FEs

hL = 1/256 (solid line is FE-error)

Matlab implementation on 3GHz Intel Core 2 Duo E8400 processor, 3.2GByte RAM, with sparse direct solver where γ ≈ 1.2

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 54 / 61

Theory: Verifying Assumptions (M1) & (M2)

For simplicity only shown for exponential covariance, i.e. ν = 1/2

• D

∇v · (k(x, ω)∇p(x, ω)) dx =

• D

f (x, ω)v dx ∀v ∈ H1

0(D)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 55 / 61

Theory: Verifying Assumptions (M1) & (M2)

For simplicity only shown for exponential covariance, i.e. ν = 1/2

• D

∇v · (k(x, ω)∇p(x, ω)) dx =

• D

f (x, ω)v dx ∀v ∈ H1

0(D)

[Barth, Schwab, Zollinger ’11]: uniformly elliptic & bounded k(·, ω) ∈ W 1,∞(D) (not satisfied here!)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 55 / 61

Theory: Verifying Assumptions (M1) & (M2)

For simplicity only shown for exponential covariance, i.e. ν = 1/2

• D

∇v · (k(x, ω)∇p(x, ω)) dx =

• D

f (x, ω)v dx ∀v ∈ H1

0(D)

[Barth, Schwab, Zollinger ’11]: uniformly elliptic & bounded k(·, ω) ∈ W 1,∞(D) (not satisfied here!) [Charrier, RS, Teckentrup ’13], [Teckentrup, RS, Giles, Ullmann ’14]: not uniformly elliptic/bounded k & low regularity

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 55 / 61

Theory: Verifying Assumptions (M1) & (M2)

For simplicity only shown for exponential covariance, i.e. ν = 1/2

• D

∇v · (k(x, ω)∇p(x, ω)) dx =

• D

f (x, ω)v dx ∀v ∈ H1

0(D)

[Barth, Schwab, Zollinger ’11]: uniformly elliptic & bounded k(·, ω) ∈ W 1,∞(D) (not satisfied here!) [Charrier, RS, Teckentrup ’13], [Teckentrup, RS, Giles, Ullmann ’14]: not uniformly elliptic/bounded k & low regularity Why more challenging?

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 55 / 61

Theory: Verifying Assumptions (M1) & (M2)

For simplicity only shown for exponential covariance, i.e. ν = 1/2

• D

∇v · (k(x, ω)∇p(x, ω)) dx =

• D

f (x, ω)v dx ∀v ∈ H1

0(D)

[Barth, Schwab, Zollinger ’11]: uniformly elliptic & bounded k(·, ω) ∈ W 1,∞(D) (not satisfied here!) [Charrier, RS, Teckentrup ’13], [Teckentrup, RS, Giles, Ullmann ’14]: not uniformly elliptic/bounded k & low regularity Why more challenging?

◮ No kmin, kmax exist such that

0 < kmin ≤ k(x, ω) ≤ kmax < ∞ a.e. in D × Ω. Cannot apply Lax-Milgram directly to ensure ∃!

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 55 / 61

Theory: Verifying Assumptions (M1) & (M2)

For simplicity only shown for exponential covariance, i.e. ν = 1/2

• D

∇v · (k(x, ω)∇p(x, ω)) dx =

• D

f (x, ω)v dx ∀v ∈ H1

0(D)

[Barth, Schwab, Zollinger ’11]: uniformly elliptic & bounded k(·, ω) ∈ W 1,∞(D) (not satisfied here!) [Charrier, RS, Teckentrup ’13], [Teckentrup, RS, Giles, Ullmann ’14]: not uniformly elliptic/bounded k & low regularity Why more challenging?

◮ No kmin, kmax exist such that

0 < kmin ≤ k(x, ω) ≤ kmax < ∞ a.e. in D × Ω. Cannot apply Lax-Milgram directly to ensure ∃!

◮ Only k(·, ω) ∈ C η(D) with η < 1

2 ⇒ p(ω, ·) /

∈ H2(D)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 55 / 61

Instead work directly (a.s. in ω) with kmin(ω) := minx k(x, ω) > 0 and kmax(ω) := maxx k(x, ω) < +∞.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 56 / 61

Instead work directly (a.s. in ω) with kmin(ω) := minx k(x, ω) > 0 and kmax(ω) := maxx k(x, ω) < +∞. Apply Lax-Milgram pathwise (a.s. in ω): ∃! p(ω) ∈ H1

0(D).

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 56 / 61

Instead work directly (a.s. in ω) with kmin(ω) := minx k(x, ω) > 0 and kmax(ω) := maxx k(x, ω) < +∞. Apply Lax-Milgram pathwise (a.s. in ω): ∃! p(ω) ∈ H1

0(D).

Use Fernique’s Theorem (for Gaussians) to show that ∀q ∈ (0, ∞) 1/kmin ∈ Lq(Ω), kmax ∈ Lq(Ω), and k ∈ Lq(Ω, C η(D)).

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 56 / 61

Instead work directly (a.s. in ω) with kmin(ω) := minx k(x, ω) > 0 and kmax(ω) := maxx k(x, ω) < +∞. Apply Lax-Milgram pathwise (a.s. in ω): ∃! p(ω) ∈ H1

0(D).

Use Fernique’s Theorem (for Gaussians) to show that ∀q ∈ (0, ∞) 1/kmin ∈ Lq(Ω), kmax ∈ Lq(Ω), and k ∈ Lq(Ω, C η(D)).

Theorem (Regularity [Charrier, RS, Teckentrup ’13])

If D is piecewise C 2, E[k] is piecewise C η, f ∈ Lq∗(Ω, Hη−1(D)), then p(ω)H1+τ(D) kmax(ω)k(·, ω)2

C η( ¯ D)

kmin(ω)4 f (·, ω)Hη−1(D), for any τ < η < 1/2. Hence p ∈ Lq(Ω, H1+τ(D)), for all q < q∗.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 56 / 61

Instead work directly (a.s. in ω) with kmin(ω) := minx k(x, ω) > 0 and kmax(ω) := maxx k(x, ω) < +∞. Apply Lax-Milgram pathwise (a.s. in ω): ∃! p(ω) ∈ H1

0(D).

Use Fernique’s Theorem (for Gaussians) to show that ∀q ∈ (0, ∞) 1/kmin ∈ Lq(Ω), kmax ∈ Lq(Ω), and k ∈ Lq(Ω, C η(D)).

Theorem (Regularity [Charrier, RS, Teckentrup ’13])

If D is piecewise C 2, E[k] is piecewise C η, f ∈ Lq∗(Ω, Hη−1(D)), then p(ω)H1+τ(D) kmax(ω)k(·, ω)2

C η( ¯ D)

kmin(ω)4 f (·, ω)Hη−1(D), for any τ < η < 1/2. Hence p ∈ Lq(Ω, H1+τ(D)), for all q < q∗.

Following “classic” proofs in [Hackbusch ’92] & [Grisvard ’85] (> 10 pages!)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 56 / 61

C´ ea’s Lemma: |p(·, ω) − ph(·, ω)|H1(D) ≤

• kmax(ω)

kmin(ω) infvh∈Vh |p(·, ω) − vh|H1(D)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 57 / 61

C´ ea’s Lemma: |p(·, ω) − ph(·, ω)|H1(D) ≤

• kmax(ω)

kmin(ω) infvh∈Vh |p(·, ω) − vh|H1(D)

Best Approximation: infvh∈Vh |p(·, ω) − vh|H1(D) ≤ C p(·, ω)H1+τ(D) hτ C = C(ω)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 57 / 61

C´ ea’s Lemma: |p(·, ω) − ph(·, ω)|H1(D) ≤

• kmax(ω)

kmin(ω) infvh∈Vh |p(·, ω) − vh|H1(D)

Best Approximation: infvh∈Vh |p(·, ω) − vh|H1(D) ≤ C p(·, ω)H1+τ(D) hτ C = C(ω) Duality Argument (for smooth, Fr´

echet-diff’ble functionals φ):

|φ (p(ω, ·)) − φ (ph(ω, ·))| ≤ kmax(ω) |p(ω, ·) − ph(ω, ·)|H1(D) |z(ω, ·) − zh(ω, ·)|H1(D)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 57 / 61

C´ ea’s Lemma: |p(·, ω) − ph(·, ω)|H1(D) ≤

• kmax(ω)

kmin(ω) infvh∈Vh |p(·, ω) − vh|H1(D)

Best Approximation: infvh∈Vh |p(·, ω) − vh|H1(D) ≤ C p(·, ω)H1+τ(D) hτ C = C(ω) Duality Argument (for smooth, Fr´

echet-diff’ble functionals φ):

|φ (p(ω, ·)) − φ (ph(ω, ·))| ≤ kmax(ω) |p(ω, ·) − ph(ω, ·)|H1(D) |z(ω, ·) − zh(ω, ·)|H1(D)

Theorem (FE Error [Charrier, RS, Teckentrup ’13])

p − phLq(Ω,H1(D)) ≤ Cτ,qf Lq∗(Ω,Hτ−1(D)) hτ, p − phLq(Ω,L2(D)) ≤ C ′

τ,qf Lq∗(Ω,Hτ−1(D)) h2τ,

   ∀ τ < η < 1

2,

∀ q < q∗.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 57 / 61

C´ ea’s Lemma: |p(·, ω) − ph(·, ω)|H1(D) ≤

• kmax(ω)

kmin(ω) infvh∈Vh |p(·, ω) − vh|H1(D)

Best Approximation: infvh∈Vh |p(·, ω) − vh|H1(D) ≤ C p(·, ω)H1+τ(D) hτ C = C(ω) Duality Argument (for smooth, Fr´

echet-diff’ble functionals φ):

|φ (p(ω, ·)) − φ (ph(ω, ·))| ≤ kmax(ω) |p(ω, ·) − ph(ω, ·)|H1(D) |z(ω, ·) − zh(ω, ·)|H1(D)

Theorem (FE Error [Teckentrup, RS, Giles, Ullmann ’14])

If there exists a Cφ ∈ Lq∗(Ω) such that φ(v) ≤ Cφ(ω)vH1−τ(D) a.s. in ω and for all v ∈ H1

0(D), then

φ(p) − φ(ph)Lq(Ω) ≤ C ′′

τ,q h2τ,

∀τ < η < 1 2 and q < q∗.

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 57 / 61

Thus, with q = 1 we get

• E
• φ(p) − φ(ph)
• = φ(p) − φ(ph)L1(Ω) = O(h2τ)

= ⇒ (M1) satisfied for any α < 1

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 58 / 61

Thus, with q = 1 we get

• E
• φ(p) − φ(ph)
• = φ(p) − φ(ph)L1(Ω) = O(h2τ)

= ⇒ (M1) satisfied for any α < 1 And with q = 2 we get V

• φ(ph) − φ(p2h)
• ≤ φ(ph) − φ(p2h)2

L2(Ω) ≤ O(h4τ)

= ⇒ (M2) satisfied for any β < 2

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 58 / 61

Thus, with q = 1 we get

• E
• φ(p) − φ(ph)
• = φ(p) − φ(ph)L1(Ω) = O(h2τ)

= ⇒ (M1) satisfied for any α < 1 And with q = 2 we get V

• φ(ph) − φ(p2h)
• ≤ φ(ph) − φ(p2h)2

L2(Ω) ≤ O(h4τ)

= ⇒ (M2) satisfied for any β < 2 Hence Cost = O(TOL−γ)

(Same as for deterministic solve!)

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 58 / 61

Thus, with q = 1 we get

• E
• φ(p) − φ(ph)
• = φ(p) − φ(ph)L1(Ω) = O(h2τ)

= ⇒ (M1) satisfied for any α < 1 And with q = 2 we get V

• φ(ph) − φ(p2h)
• ≤ φ(ph) − φ(p2h)2

L2(Ω) ≤ O(h4τ)

= ⇒ (M2) satisfied for any β < 2 Hence Cost = O(TOL−γ)

(Same as for deterministic solve!)

Hence optimal and robust deterministic solver with γ = d crucial!

This is a whole talk in itself!

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 58 / 61

Numerical Confirmation

D = (0, 1)2; covariance R(x, y) = exp

• − x−y2

0.3

• ; std. FEs, circulant embedding
• E
• φ(1)(p) − φ(1)(ph)
• where φ(1)(p) := Lω(Ψ) − bω(Ψ, v)

given Ψ(x) = x (outflow on right)

V

• φ(2)(ph) − φ(2)(p2h)
• where φ(2)(p) :=
• 1

|D∗|

• D∗ p(x) dx

2 (i.e. 2nd moment of p over small patch)

α = 1 and β = 2

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 59 / 61

Other developments in MLMC

discontinuous coefficients, point evaluations, mixed FEs, FVM, . . . Many other PDEs and applications can optimise all parameters (not just {Nℓ}) [Haji Ali, Tempone] adaptivity [Von Schwerin, Tempone et al], [Kornhuber, Youett] variance estimation [Bierig, Chernov]

• ptimal estimation of CDFs, PDFs [Giles, Nagapetyan, Ritter]

hybrid with stochastic collocation [Tesei, Nobile et al] generalisation to general multilevel quadrature [Harbrecht et al] multilevel QMC [Kuo, Schwab, Sloan] multilevel SC [Teckentrup, Jantsch, Webster, Gunzburger] MLMCMC [Hoang, Schwab, Stuart], [Dodwell, Ketelson, RS, Teckentrup] etc . . .

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 60 / 61

Some Key References

Heinrich, MC complexity of global solution of integral eq., J Complex, 1998 Giles, Multilevel Monte Carlo path simulation, Oper Res, 2008 Giles, Waterhouse, MLQMC path simulation, Radon Series, deGruyter, 2009 Cliffe, Giles, RS, Teckentrup, MLMC methods & applications to elliptic PDEs with random coefficients, Comput Visual Sci, 2011 Barth, Schwab, Zollinger, MLMC FE method for elliptic PDEs with stochastic coefficients, Numer Math, 2011 Charrier, RS, Teckentrup, FE error analysis of elliptic PDEs with random coefficients and its application to MLMC methods, SINUM, 2013 Teckentrup, RS, Giles, Ullmann, Further analysis of MLMC methods for elliptic PDEs with random coefficients, Numer Math, 2013 Kuo, Schwab, Sloan, MLQMC FE methods for a class of elliptic PDEs with random coefficients, FoCoM, 2015 Kuo, RS, Schwab, Sloan, Ullmann, MLQMC methods for lognormal diffusion problems, Math Comput, 2017 Teckentrup, Jantsch, Webster, Gunzburger, A Multilevel Stochastic Collocation Method for PDEs with Random Input Data, JUQ, 2015

• R. Scheichl (Heidelberg)

Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 61 / 61