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: Borehole accidentally drilled into 1 the repository. Radionuclides released into 2 Culebra Dolomite and transported by groundwater. Estimate travel time from release 3 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) EDZ CROWN SPACE WASTE VAULTS FAULTED GRANITE GRANITE DEEP SKIDDAW N-S SKIDDAW DEEP LATTERBARROW � Darcy’s Law: q + k ( x , ω ) ∇ p � = f ( x , ω ) N-S LATTERBARROW FAULTED TOP M-F BVG TOP M-F BVG FAULTED BLEAWATH BVG BLEAWATH BVG Incompressibility: ∇ · � = g ( x , ω ) q FAULTED F-H BVG F-H BVG FAULTED UNDIFF BVG UNDIFF BVG FAULTED N-S BVG + Boundary Conditions N-S BVG FAULTED CARB LST CARB LST FAULTED COLLYHURST COLLYHURST FAULTED BROCKRAM � Uncertainty in k = ⇒ Uncertainty in p & q BROCKRAM SHALES + EVAP FAULTED BNHM BOTTOM NHM Stochastic Modelling! FAULTED DEEP ST BEES DEEP ST BEES FAULTED N-S ST BEES N-S ST BEES FAULTED VN-S ST BEES VN-S ST BEES FAULTED DEEP CALDER DEEP CALDER FAULTED N-S CALDER N-S CALDER FAULTED VN-S CALDER VN-S CALDER MERCIA MUDSTONE QUATERNARY Geology at Sellafield (former potential UK radwaste site) � NIREX UK Ltd. c 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): � 1 � R ( x , y ) := σ 2 ρ ν λ � x − y � Usually: σ 2 > 1, λ < 1 & ν < 1 e.g. ρ 1 / 2 ( r ) = exp( − r ) or ρ ∞ ( r ) = exp( − r 2 ) Compute statistics (eg. moments) of 64 , σ 2 = 8) 1 functionals of p and � q , e.g. realisation ( λ = k ( x ) ◮ point evaluations, norms, averages k ( y ) = O (10 9 ) contrast: max x , y ◮ breakthrough time (to boundary) (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 � 2 √ ν r � ν � 2 √ ν r � σ 2 c ( x , y ) = c θ ( r ) = K ν , r = � x − y � 2 2 ν − 1 Γ( ν ) λ λ K ν : modified Bessel function of order ν σ 2 : variance θ = ( σ 2 , λ, ν ) Parameters λ : 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 � 2 √ ν r � ν � 2 √ ν r � σ 2 c ( x , y ) = c θ ( r ) = K ν , r = � x − y � 2 2 ν − 1 Γ( ν ) λ λ K ν : modified Bessel function of order ν σ 2 : variance θ = ( σ 2 , λ, ν ) Parameters λ : correlation length ν : smoothness parameter Special cases: √ c ( r ) = σ 2 exp( − ν = 1 2 : 2 r /λ ) exponential covariance c ( r ) = σ 2 � 2 r � � 2 r � ν = 1 : K 1 Bessel covariance λ λ c ( r ) = σ 2 exp( − r 2 /λ 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 ν =1/2 ν =1 ν =2 0.8 ν =3 Gauss 0.6 c(r) 0.4 0.2 0 0 1 2 3 4 5 r σ 2 = 8, λ = 1 / 64, ν = 1 / 2 σ 2 = 1, λ = 1 Smoothness: Realisations Z ( · , ω ) ∈ C η ( D ) (H¨ older), 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 : L 2 ( D ) → L 2 ( D ) , ( Cu )( x ) = u ( y ) c ( x , y ) d y D is selfadjoint, compact, nonnegative. Hence, the eigenvalues { µ m } m ∈ N form a non-increasing sequence, accumulating at most at 0. eve expansion (converges in L 2 P (Ω; L ∞ ( D ))) : Karhunen-Lo` ∞ � √ µ m φ m ( x ) Y m ( ω ) Z ( x , ω ) = Z ( x ) + m =1 where { φ m } m ∈ N are normalised eigenfunctions and Y m ∼ 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 : L 2 ( D ) → L 2 ( D ) , ( Cu )( x ) = u ( y ) c ( x , y ) d y D is selfadjoint, compact, nonnegative. Hence, the eigenvalues { µ m } m ∈ N form a non-increasing sequence, accumulating at most at 0. eve expansion (converges in L 2 P (Ω; L ∞ ( D ))) : Truncated Karhunen-Lo` s � √ µ m φ m ( x ) Y m ( ω ) Z s ( x , ω ) = Z ( x ) + m =1 where { φ m } m ∈ N are normalised eigenfunctions and Y m ∼ 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 x ∈ D ⊂ R d , ω ∈ Ω (prob. space) −∇ · ( k ( x , ω ) ∇ p ( x , ω )) = f ( x , ω ) , R. Scheichl (Heidelberg) Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 18 / 61
Computational Challenges Simulating PDEs with Highly Heterogeneous Random Coefficients x ∈ D ⊂ R d , ω ∈ Ω (prob. space) −∇ · ( k ( x , ω ) ∇ p ( x , ω )) = f ( x , ω ) , 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 x ∈ D ⊂ R d , ω ∈ Ω (prob. space) −∇ · ( k ( x , ω ) ∇ p ( x , ω )) = f ( x , ω ) , 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 x ∈ D ⊂ R d , ω ∈ Ω (prob. space) −∇ · ( k ( x , ω ) ∇ p ( x , ω )) = f ( x , ω ) , 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 k max / k min > 10 6 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 k max / k min > 10 6 Small λ = ⇒ multiscale + high stochast. dimension s > 100 ECDF of log travel time 20 000 MC samples 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 M = 30 M = 100 0.2 M = 500 M = 1000 0.1 Source: Ernst et al, 2014 ( s = M ) 0 4 4.5 5 5.5 6 log t 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 ( · , ω ) ∈ H 1 0 ( D ) s.t. � � ∀ v ∈ H 1 ∇ v · ( k ( x , ω ) ∇ p ( x , ω )) d x = f ( x , ω ) v d x , 0 ( D ) . D D ∃ ! p ( · , ω ) ∈ H 1 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 ( · , ω ) ∈ H 1 0 ( D ) s.t. � � ∀ v ∈ H 1 ∇ v · ( k ( x , ω ) ∇ p ( x , ω )) d x = f ( x , ω ) v d x , 0 ( D ) . D D ∃ ! p ( · , ω ) ∈ H 1 0 ( D ) a.s. in ω ∈ Ω (Lax-Milgram pathwise). Let V h ⊂ H 1 0 ( D ) be the space of continuous, piecewise linear FEs w.r.t. a mesh T h with mesh width h > 0. Find p h ( · , ω ) ∈ V h that satisfies weak form for all v h ∈ V h . 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 ( · , ω ) ∈ H 1 0 ( D ) s.t. � � ∀ v ∈ H 1 ∇ v · ( k ( x , ω ) ∇ p ( x , ω )) d x = f ( x , ω ) v d x , 0 ( D ) . D D ∃ ! p ( · , ω ) ∈ H 1 0 ( D ) a.s. in ω ∈ Ω (Lax-Milgram pathwise). Let V h ⊂ H 1 0 ( D ) be the space of continuous, piecewise linear FEs w.r.t. a mesh T h with mesh width h > 0. Find p h ( · , ω ) ∈ V h that satisfies weak form for all v h ∈ V h . Write p h ( x , ω ) := � M h i =1 P i ϕ i ( x ). Then this is equivalent to the random matrix system A ( ω ) P ( ω ) = F ( ω ) with � � A i , j ( ω ) := ∇ ϕ j · ( k ( x , ω ) ∇ ϕ i ) d x , F i ( ω ) := f ( x , ω ) ϕ i d x D D R. Scheichl (Heidelberg) Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 20 / 61
Computational Challenges Standard Monte Carlo Quadrature Model( h ) Output Y ( ω ) ∈ R s P ( ω ) ∈ R M h − → − → Q h , 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 Model( h ) Output Y ( ω ) ∈ R s P ( ω ) ∈ R M h − → − → Q h , s ( ω ) ∈ R random input state vector quantity of interest Real QoI Q ( ω ) inaccessible (exact PDE) , but we can assume h → 0 , s →∞ E [ Q ] and | E [ Q h , s − Q ] | = O ( h α ) + O ( s − α ′ ) E [ Q h , s ] − → Standard Monte Carlo estimator for E [ Q ]: More detail below! � N Q MC := 1 Q ( i ) ˆ h , s N i =1 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 Model( h ) Output Y ( ω ) ∈ R s P ( ω ) ∈ R M h − → − → Q h , s ( ω ) ∈ R random input state vector quantity of interest Real QoI Q ( ω ) inaccessible (exact PDE) , but we can assume h → 0 , s →∞ E [ Q ] and | E [ Q h , s − Q ] | = O ( h α ) + O ( s − α ′ ) E [ Q h , s ] − → Standard Monte Carlo estimator for E [ Q ]: More detail below! � N Q MC := 1 Q ( i ) ˆ h , s N i =1 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) : � � 2 �� ˆ � 2 � V [ Q h , s ] Q MC − E [ Q ] E = + E [ Q h , s − Q ] � �� � N � �� � � �� � =: MSE sampling error model error (“bias”) Typical: α = 1 ⇒ MSE = O ( N − 1 ) + O ( h 2 ) ≤ TOL 2 , 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) : � � 2 �� ˆ � 2 � V [ Q h , s ] Q MC − E [ Q ] E = + E [ Q h , s − Q ] � �� � N � �� � � �� � =: MSE sampling error model error (“bias”) Typical: α = 1 ⇒ MSE = O ( N − 1 ) + O ( h 2 ) ≤ TOL 2 , 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 ∼ 10 6 and Cost = O (10 12 ) 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 ∂ x 1 � L 1 ( 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 ∂ x 1 � L 1 ( D ) using mixed FEs and the AMG solver amg1r5 [Ruge, St¨ uben, 1992] Numerically observed FE-error: ≈ O ( h 3 / 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 ∂ x 1 � L 1 ( D ) using mixed FEs and the AMG solver amg1r5 [Ruge, St¨ uben, 1992] Numerically observed FE-error: ≈ O ( h 3 / 4 ) = ⇒ α ≈ 3 / 4. Numerically observed cost/sample: ≈ O ( h − 2 ) = ⇒ γ ≈ 1. ≈ O (TOL − 14 / 3 ) Total cost to get RMSE O (TOL): 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 ∂ x 1 � L 1 ( D ) using mixed FEs and the AMG solver amg1r5 [Ruge, St¨ uben, 1992] Numerically observed FE-error: ≈ O ( h 3 / 4 ) = ⇒ α ≈ 3 / 4. Numerically observed cost/sample: ≈ O ( h − 2 ) = ⇒ γ ≈ 1. ≈ O (TOL − 14 / 3 ) Total cost to get RMSE O (TOL): to get error reduction by a factor 2 → cost grows by a factor 25! Case 1: σ 2 = 1, λ = 0 . 3, ν = 0 . 5 Case 2: σ 2 = 3, λ = 0 . 1, ν = 0 . 5 h − 1 h − 1 TOL N Cost TOL N Cost 1 . 4 × 10 4 8 . 5 × 10 3 0 . 01 129 21 min 0 . 01 513 4 h 3 . 5 × 10 5 0 . 002 1025 30 days 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: N SC = 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: N SC = 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: N SC = 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 ℓ (simple geometric arguments) π d 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 ℓ (simple geometric arguments) π d 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 ℓ (simple geometric arguments) π d 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 { X k } of i.i.d. copies of a given random variable X , basic MC simulation uses the estimator E [ X ] ≈ S N N , S N = X 1 + · · · + X N . S N By the Strong Law of Large Numbers : 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 { X k } of i.i.d. copies of a given random variable X , basic MC simulation uses the estimator E [ X ] ≈ S N N , S N = X 1 + · · · + X N . S N By the Strong Law of Large Numbers : N → E [ X ] a.s. N � Also, for any measurable function f : 1 f ( X k ) → E [ f ( X )] a.s. N k =1 R. Scheichl (Heidelberg) Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 28 / 61
Basic Monte Carlo Simulation Given a sequence { X k } of i.i.d. copies of a given random variable X , basic MC simulation uses the estimator E [ X ] ≈ S N N , S N = X 1 + · · · + X N . S N By the Strong Law of Large Numbers : N → E [ X ] a.s. N � Also, for any measurable function f : 1 f ( X k ) → E [ f ( X )] a.s. N k =1 If E [ X ] = µ and Var [ X ] = σ 2 , then (via the Central Limit Theorem ) N = S N − N µ Var [ S N ] = N σ 2 and S ∗ E [ S N ] = N µ, √ → N (0 , 1) N σ 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 Since 1 �� S N � 2 � N = σ 2 = Var S N E N − µ N → 0 , we have mean square convergence of S N / N to µ . R. Scheichl (Heidelberg) Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 29 / 61
Monte Carlo Convergence Statements Since 1 �� S N � 2 � N = σ 2 = Var S N E N − µ N → 0 , we have mean square convergence of S N / N to µ . Also Chebyshev’s Inequality implies, for any ǫ > 0, 2 �� � � � � ≤ σ 2 S N � � � > N − 1 / 2+ ǫ P N − µ N 2 ǫ , � 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 Since 1 �� S N � 2 � N = σ 2 = Var S N E N − µ N → 0 , we have mean square convergence of S N / N to µ . Also Chebyshev’s Inequality implies, for any ǫ > 0, 2 �� � � � � ≤ σ 2 S N � � � > N − 1 / 2+ ǫ P N − µ N 2 ǫ , � i.e. for any ǫ > 0, the probability of the error being larger than N − 1 / 2+ ǫ converges to zero as N → ∞ . � | X − µ | 3 � If ρ := E < ∞ , then the Berry-Esseen Inequality gives 3 � � ρ � � � 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 � � 0 , 1 3 , 2 3 , 1 9 , 4 9 , 7 9 , 2 9 , 5 9 , 8 9 , 1 27 , 10 { x n } 11 n =1 = . 27 (for base 3, it is given by x n = k 3 j , where j increases monotonically and k runs through all nonnegative integers such that k / 3 j is an irreducible fraction < 1. The ordering in k is obtained by representing k in base 3 and reversing the digits) 0.1 0 − 0.1 0 0.2 0.4 0.6 0.8 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 � � 0 , 1 3 , 2 3 , 1 9 , 4 9 , 7 9 , 2 9 , 5 9 , 8 9 , 1 27 , 10 { x n } 11 n =1 = . 27 (for base 3, it is given by x n = k 3 j , where j increases monotonically and k runs through all nonnegative integers such that k / 3 j is an irreducible fraction < 1. The ordering in k is obtained by representing k in base 3 and reversing the digits) 0.1 0 − 0.1 0 0.2 0.4 0.6 0.8 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 � � 0 , 1 3 , 2 3 , 1 9 , 4 9 , 7 9 , 2 9 , 5 9 , 8 9 , 1 27 , 10 { x n } 11 n =1 = . 27 (for base 3, it is given by x n = k 3 j , where j increases monotonically and k runs through all nonnegative integers such that k / 3 j is an irreducible fraction < 1. The ordering in k is obtained by representing k in base 3 and reversing the digits) 0.1 0 − 0.1 0 0.2 0.4 0.6 0.8 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 { X k } N k =1 of X and compute the sample average S N / 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 { X k } N k =1 of X and compute the sample average S N / N . Now, let { ˜ X k } N k =1 be a second set of samples of X with sample average ˜ S N / 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 { X k } N k =1 of X and compute the sample average S N / N . Now, let { ˜ X k } N k =1 be a second set of samples of X with sample average ˜ S N / N . 2 ( S N / N + ˜ Since both averages converge to E [ X ], so does 1 S N / N ). When X k and ˜ X k are negatively correlated they are called antithetic 2 N ( S N + ˜ 1 samples, and the approximation S N ) is a more reliable 1 approximation of E [ X ] than 2 N S 2 N . 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 { X k } and { ˜ X k } be identically distributed with Cov ( X j , X k ) = Cov ( ˜ X j , ˜ X k ) = 0 for j � = k . Then � � � � � S 2 N � � S N � S N + ˜ ˜ S N + 1 S N S N Var = Var 2 Cov N , ≤ Var . 2 N 2 N N 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 { X k } and { ˜ X k } be identically distributed with Cov ( X j , X k ) = Cov ( ˜ X j , ˜ X k ) = 0 for j � = k . Then � � � � � S 2 N � � S N � S N + ˜ ˜ S N + 1 S N S N Var = Var 2 Cov N , ≤ Var . 2 N 2 N N 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 { X k } and { ˜ X k } be identically distributed with Cov ( X j , X k ) = Cov ( ˜ X j , ˜ X k ) = 0 for j � = k . Then � � � � � S 2 N � � S N � S N + ˜ ˜ S N + 1 S N S N Var = Var 2 Cov N , ≤ Var . 2 N 2 N N N Worst case: Variance of average of N samples and N antithetic samples no better than variance of N independent samples. Best case: S N / N and ˜ S N / N negatively correlated; therefore variance of N samples and N antithetic samples is smaller than variance of 2 N 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 1 (1 − u 2 ) � u 1 u = ˙ = , u (0) = u 0 . u 2 ˙ u 2 ( u 1 − 1) Let the initial condition u 0 be uncertain and modeled as a (uniform) random vector u 0 ∼ U (Γ), where Γ denotes the square Γ = u 0 + [ − ǫ, ǫ ] 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 1 (1 − u 2 ) � u 1 u = ˙ = , u (0) = u 0 . u 2 ˙ u 2 ( u 1 − 1) Let the initial condition u 0 be uncertain and modeled as a (uniform) random vector u 0 ∼ U (Γ), where Γ denotes the square Γ = u 0 + [ − ǫ, ǫ ] 2 . Goal: estimate E [ u 1 ( T )] at time T > 0. Denote by u h = u h ( ω ) the explicit Euler approximation after M time steps of length h = T M starting with initial data u 0 = u 0 ( ω ). Define the QoI Q = φ ( u ( T )) = u 1 ( T ) for u = [ u 1 , u 2 ] T and estimate E [ Q h ] using the MC method just described, where Q h = φ ( u h ). 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 1 (1 − u 2 ) � u 1 u = ˙ = , u (0) = u 0 . u 2 ˙ u 2 ( u 1 − 1) Let the initial condition u 0 be uncertain and modeled as a (uniform) random vector u 0 ∼ U (Γ), where Γ denotes the square Γ = u 0 + [ − ǫ, ǫ ] 2 . Goal: estimate E [ u 1 ( T )] at time T > 0. Denote by u h = u h ( ω ) the explicit Euler approximation after M time steps of length h = T M starting with initial data u 0 = u 0 ( ω ). Define the QoI Q = φ ( u ( T )) = u 1 ( T ) for u = [ u 1 , u 2 ] T and estimate E [ Q h ] using the MC method just described, where Q h = φ ( u h ). 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 [ Q h ] by N � Q h , N = 1 Q ( k ) Q h := � � h N k =1 (i.e the average over N samples { Q ( k ) h } N k =1 of Q h 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 [ Q h ] by N � Q h , N = 1 Q ( k ) Q h := � � h N k =1 (i.e the average over N samples { Q ( k ) h } N k =1 of Q h with step size h ) Expect better approximations for N large and h small. Error with N samples and M = T / h time steps: e h , N = | E [ Q ] − � + | E [ Q h ] − � Q h | ≤ | E [ Q ] − E [ Q h ] | Q h | � �� � � �� � discretisation error 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 [ Q h ] by N � Q h , N = 1 Q ( k ) Q h := � � h N k =1 (i.e the average over N samples { Q ( k ) h } N k =1 of Q h with step size h ) Expect better approximations for N large and h small. Error with N samples and M = T / h time steps: e h , N = | E [ Q ] − � + | E [ Q h ] − � Q h | ≤ | E [ Q ] − E [ Q h ] | Q h | � �� � � �� � discretisation error Monte Carlo error Can show also that the mean square error satisfies (with equality!) �� � 2 � � � 2 + Var [ Q h ] E [ Q ] − � E Q h = E [ Q − Q h ] 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 ) − u h � ≤ 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 ) − u h � ≤ Kh . φ Lipschitz-continuous with constant L = 1: | φ ( u ( T )) − φ ( u h ) | ≤ 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 ) − u h � ≤ Kh . φ Lipschitz-continuous with constant L = 1: | φ ( u ( T )) − φ ( u h ) | ≤ KLh . Therefore | E [ Q ] − E [ Q h ] | = | E [ Q − Q h ] | ≤ 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 [ Q h ] = σ 2 h we get (see the Problem Sheet ) �� � � � ≤ 1 . 96 σ h � � � E [ Q h ] − � > 0 . 95 + O ( N − 1 / 2 ) P Q h , N √ N 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 [ Q h ] = σ 2 h we get (see the Problem Sheet ) �� � � � ≤ 1 . 96 σ h � � � E [ Q h ] − � > 0 . 95 + O ( N − 1 / 2 ) P Q h , N √ N Combined with the discretisation error (using triangle inequality): � � e h , N ≤ 2 KLh + 1 . 96 σ h > 0 . 95 + O ( N − 1 / 2 ) . √ P N 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 [ Q h ] = σ 2 h we get (see the Problem Sheet ) �� � � � ≤ 1 . 96 σ h � � � E [ Q h ] − � > 0 . 95 + O ( N − 1 / 2 ) P Q h , N √ N Combined with the discretisation error (using triangle inequality): � � e h , N ≤ 2 KLh + 1 . 96 σ h > 0 . 95 + O ( N − 1 / 2 ) . √ P N Balancing discretization and MC errors, i.e. KLh ≈ TOL 1 . 96 σ h ≈ TOL √ and , 2 2 N leads to N ≈ 16 σ 2 h ≈ TOL h and Cost = O ( MN ) = O (TOL − 3 ) 2 KL , TOL 2 R. Scheichl (Heidelberg) Stochastic Uncertainty Modelling & MLMC RICAM Tutorial 14/12/18 36 / 61
Easy Example: Predator-Prey Dynamical System Sample trajectories 3 2.5 2 u 2 1.5 1 0.5 0 0 0.5 1 1.5 2 2.5 3 u 1 Predator-prey system integrated from 0 to T = 6 with u 0 = [0 . 5 , 2] T and ǫ = 0 . 2. Unperturbed trajectory (black) and 15 perturbed trajectories. (for the unperturbed trajectory u 1 ( 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 u 0 ∼ U (Γ), then the same holds for the random vector ˜ u 0 := 2 u 0 − u 0 . Thus, the trajectories generated by the random initial data ˜ u 0 have the same distribution as those generated by u 0 . 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 u 0 ∼ U (Γ), then the same holds for the random vector ˜ u 0 := 2 u 0 − u 0 . Thus, the trajectories generated by the random initial data ˜ u 0 have the same distribution as those generated by u 0 . Let Q h = φ ( u h ) be the basic samples and ˜ Q h = φ (˜ u h ) 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 u 0 ∼ U (Γ), then the same holds for the random vector ˜ u 0 := 2 u 0 − u 0 . Thus, the trajectories generated by the random initial data ˜ u 0 have the same distribution as those generated by u 0 . Let Q h = φ ( u h ) be the basic samples and ˜ Q h = φ (˜ u h ) the antithetic counterparts. Note that all pairs of samples are independent except each sample and its antithetic counterpart. Q h , N + � 2 ( � ˜ � 1 Use Q h , N ) instead of Q h , 2 N (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 1.52 1.52 1.51 1.51 1.5 1.5 u 1 (t end ) u 1 (t end ) 1.49 1.49 1.48 1.48 1.47 1.47 1.46 1.46 0 0.5 1 1.5 2 2.5 3 3.5 4 0 0.5 1 1.5 2 2.5 3 3.5 4 N N 4 4 x 10 x 10 Estimates of E [ u 1 ( T )] using standard MC with N samples (left) vs. antithetic sampling with N / 2 pairs of samples (right) (with 95% confidence intervals) Var [ Q h ] and Cov ( Q h , ˜ Q h ) estimated using sample variance and covariance (resp.), i.e. � N � N k =1 ( Q ( k ) Q h , N ) 2 and k =1 ( Q ( k ) Q ( k ) − � 1 − � 1 − � Q h , N )( ˜ ˜ Q h , N ) N − 1 h N − 1 h h 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 Q h ≈ Q are computable, where h > 0 denotes a discretization parameter (time step size, grid size, . . . ) and h → 0 E [ Q h ] = E [ Q ] . lim 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 Q h ≈ Q are computable, where h > 0 denotes a discretization parameter (time step size, grid size, . . . ) and h → 0 E [ Q h ] = E [ Q ] . lim More precisely, we assume error to converge in mean with rate α : | E [ Q h − 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
Recommend
More recommend
