SLIDE 1 Mathematical Institute University of Oxford MLQMC Methods for Elliptic PDEs Driven by White Noise
- M. Croci (Oxford), M. B. Giles (Oxford), M. E. Rognes (Simula), P. E. Farrell (Oxford)
SIAM CSE and ICIAM 2019
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 2
Mathematical Institute University of Oxford Overview
Introduction and Background Multilevel Monte Carlo Multilevel Quasi Monte Carlo Conclusions
SLIDE 3 Mathematical Institute University of Oxford Introduction
The motivation of our research is the solution of partial differential equations with random coefficients via the multilevel Quasi-Monte Carlo method. Today we focus on the ubiquitous model problem, −∇ · (eu(x,ω)∇p(x, ω)) = 1, x ∈ G ⊂ Rd, ω ∈ Ω, p(x, ω) = 0, x ∈ ∂G, ω ∈ Ω. Here u(x, ω) ∼ N(0, C(x, y)) is a (Mat´ ern) Gaussian random field and we are interested in computing E[P], where P(ω) = ||p||2
L2(G)(ω).
Simple approach: standard Monte Carlo (MC) method, E[P] ≈ 1 N
N
P(ωn). Convergence rate O(N−1/2). O(ε−q) cost per sample ⇒ O(ε−2−q) complexity for ε tol.
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 4 Mathematical Institute University of Oxford From standard Monte Carlo to Quasi Monte Carlo
Approximate E[P] with an s-dimensional integral over [0, 1]s, E[P] ≈
N
N
Y (xn),
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
Pseudo-random points. Low-discrepancy point sequence (Sobol’).
QMC convergence rate up to O(N−1) ⇒ up to O(ε−1−q) complexity.
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 5 Mathematical Institute University of Oxford Randomised Quasi Monte Carlo
E[P] ≈
M
M
N
N
Y (ˆ xn,m)
where {ˆ xn,m}N
n=1 is the m-th randomisation of a low-discrepancy point set {xn}N n=1.
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
First N = 256 Sobol’ points before (left) and after (right) scrambling.
QMC convergence rate up to O(N−1) ⇒ up to O(ε−1−q) complexity.
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 6 Mathematical Institute University of Oxford Multilevel Monte Carlo [Giles 2008]
Solve for P using a hierarchy of L + 1 (possibly non-nested) meshes to obtain the approximations Pℓ of different accuracy for ℓ = 0, . . . , L, then E[P] ≈ E[PL] = E[P0] +
L
E[Pℓ − Pℓ−1].
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 7 Mathematical Institute University of Oxford Multilevel Monte Carlo [Giles 2008]
Solve for P using a hierarchy of L + 1 (possibly non-nested) meshes to obtain the approximations Pℓ of different accuracy for ℓ = 0, . . . , L, then E[P] ≈ E[PL] = E[P0] +
L
E[Pℓ − Pℓ−1]. Apply standard MC to each term on the RHS: E[P] ≈ 1 N0
N0
P0(ωn
0) + L
1 Nℓ
Nℓ
[Pℓ(ωn
ℓ ) − Pℓ−1(ωn ℓ )].
Under suitable conditions = ⇒ optimal Nℓ known and O(ε−2) complexity, O(ε−q) better than MC.
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 8 Mathematical Institute University of Oxford Multilevel Quasi Monte Carlo [Giles and Waterhouse 2009]
We now approximate each expectation in the telescoping sum with randomised quasi Monte Carlo (QMC). Rewrite E[Pℓ − Pℓ−1] as the sℓ-dimensional integral over [0, 1]sℓ, E[Pℓ − Pℓ−1] =
Yℓ(x)dx ≈ 1 M
M
Nℓ
Nℓ
Yℓ(ˆ xℓ
n,m)
Yℓ, where {ˆ xℓ
n,m}Nℓ n=1 is the m-th randomisation of a low-discrepancy point set {xℓ n}Nℓ n=1. We
use M = 32 and random digital shifted Sobol’ sequences1. QMC integration convergence2 up to O(N−1
ℓ
) = ⇒ complexity up to O(ε−1).
1Python-wrapped Intel
R
MKL library Sobol’ sequence implementation augmented with Joe and
Kuo’s primitive polynomials and direction numbers (smax = 21201) [Joe and Kuo 2008].
2Lots of recent work with randomly shifted lattice rules [Kuo, Schwab, Sloan 2015, Kuo, Scheichl,
Schwab, Sloan, Ullmann 2017, Herrmann and Schwab 2017,...].
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 9
Mathematical Institute University of Oxford Gaussian field sampling
Sampling the Gaussian field u(x, ω) ∼ N(0, C) is hard!
Lots of recent developments related to ML(Q)MC: [Kuo et al. 2015, Kuo et al. 2018, Graham et al. 2018, Drzisga et al. 2017, Herrmann and Schwab 2017, Osborn et al. 2018, C. et al. 2018, ...]
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 10
Mathematical Institute University of Oxford Gaussian field sampling
Sampling the Gaussian field u(x, ω) ∼ N(0, C) is hard!
Na¨ ıve approach
Discretise u and compute a Cholesky factorization of the covariance matrix.
Better approaches
Karhunen-Lo` eve. FFT + circulant embeddings. [Dietrich and Newsam 1997] PDE approach [Lindgren et al. 2009] (see next).
Lots of recent developments related to ML(Q)MC: [Kuo et al. 2015, Kuo et al. 2018, Graham et al. 2018, Drzisga et al. 2017, Herrmann and Schwab 2017, Osborn et al. 2018, C. et al. 2018, ...]
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 11 Mathematical Institute University of Oxford PDE approach to Mat´ ern field sampling [Lindgren et al. 2009]
If the covariance of u is of the Mat´ ern class with smoothness parameter ν, then u approximately satisfies the linear elliptic SPDE, Lu(x, ω) =
k u(x, ω) = η ˙ W , x ∈ D ⊂ Rd, ν = 2k − d/2 > 0, where ˙ W is spatial white noise, G ⊂ D, η ∈ R and for today k ∈ N, η = 1.
Spatial white noise in [0, 1]2.
Refs: [Abrahamsen 1997, Scheuerer 2010, Lindgren et al. 2009, Bolin et al. 2017, Khristenko et al. 2018]
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 12 Mathematical Institute University of Oxford PDE approach to Mat´ ern field sampling [Lindgren et al. 2009] Definition (Spatial White Noise ˙ W (Hida et al. 1993))
For any φ ∈ L2(D), define ˙ W , φ :=
W . For any φi, φj ∈ L2(D), bi = ˙ W , φi, bj = ˙ W , φj are zero-mean Gaussian random variables, with, E[bibj] =
φiφj dx =: Mij, b ∼ N(0, M). (1)
Spatial white noise in [0, 1]2.
Refs: [Abrahamsen 1997, Scheuerer 2010, Lindgren et al. 2009, Bolin et al. 2017, Khristenko et al. 2018]
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 13
Mathematical Institute University of Oxford Overview
Introduction and Background Multilevel Monte Carlo Multilevel Quasi Monte Carlo Conclusions
SLIDE 14
Mathematical Institute University of Oxford Efficient white noise sampling for MLMC [C. et al. 2018]
Let {Vℓ}L
ℓ=0 be a hierarchy of (possibly non-nested) FEM approximation subspaces with
Vℓ = span({φℓ
i } nℓ
dofs
i=1 ) ⊆ H1 0(D). On each MLMC level, we need to solve for uℓ and uℓ−1,
Luℓ = ˙ W , and Luℓ−1 = ˙ W , if ℓ > 0, where we use the same white noise sample on both levels to enforce the MLMC coupling.
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 15 Mathematical Institute University of Oxford Efficient white noise sampling for MLMC [C. et al. 2018]
Let {Vℓ}L
ℓ=0 be a hierarchy of (possibly non-nested) FEM approximation subspaces with
Vℓ = span({φℓ
i } nℓ
dofs
i=1 ) ⊆ H1 0(D). On each MLMC level, we need to solve for uℓ and uℓ−1,
Luℓ = ˙ W , and Luℓ−1 = ˙ W , if ℓ > 0, where we use the same white noise sample on both levels to enforce the MLMC coupling. After discretisation, the above yields the linear system,
Aℓ−1 uℓ uℓ−1
bℓ−1
where bℓ
i = ˙
W , φℓ
i .
Hence b ∼ N(0, M), where M is the mass matrix over Vℓ + Vℓ−1, (set V−1 = ∅), i.e. M =
Mℓ,ℓ−1 (Mℓ,ℓ−1)T Mℓ−1
Mℓ,ℓ−1
ij
=
i φℓ−1 j
dx, if ℓ > 0. NOTE: we do not require the FEM approximation subspaces to be nested!
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 16
Mathematical Institute University of Oxford How to sample b?
Sampling b is hard!
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 17 Mathematical Institute University of Oxford How to sample b?
Sampling b is hard!
Na¨ ıve approach
Factorise the covariance M using Cholesky (cubic complexity!) Works well if M diagonal. Previous work under this assumption [Lindgren et
- al. 2009, Osborn et al. 2017, Drzisga, et al. 2017, Du and Zhang 2002].
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 18 Mathematical Institute University of Oxford How to sample b?
Sampling b is hard!
Na¨ ıve approach
Factorise the covariance M using Cholesky (cubic complexity!) Works well if M diagonal. Previous work under this assumption [Lindgren et
- al. 2009, Osborn et al. 2017, Drzisga, et al. 2017, Du and Zhang 2002].
Our Work [C. et al. 2018]
We do not require M to be diagonal. We can sample b with linear complexity.
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 19
Mathematical Institute University of Oxford Sampling b in linear complexity (Sketch) [C. et al. 2018]
Two ingredients: supermeshing [Farrell 2009] and local factorisation [Wathen 1987].
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 20 Mathematical Institute University of Oxford Sampling b in linear complexity (Sketch) [C. et al. 2018]
Two ingredients: supermeshing [Farrell 2009] and local factorisation [Wathen 1987].
- 1. Supermeshing: construct a FEM subspace Sℓ such that Vℓ and Vℓ−1 are both
nested within Sℓ. This requires a supermesh construction: Sample bℓ
S ∼ N(0, Mℓ S) where Mℓ S is the mass matrix over Sℓ and get bℓ and bℓ−1
by transferring bℓ
S onto Vℓ and Vℓ−1 using nested interpolation. EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 21 Mathematical Institute University of Oxford Sampling b in linear complexity (Sketch) [C. et al. 2018]
Two ingredients: supermeshing [Farrell 2009] and local factorisation [Wathen 1987].
- 1. Supermeshing: construct a FEM subspace Sℓ such that Vℓ and Vℓ−1 are both
nested within Sℓ. This requires a supermesh construction: Sample bℓ
S ∼ N(0, Mℓ S) where Mℓ S is the mass matrix over Sℓ and get bℓ and bℓ−1
by transferring bℓ
S onto Vℓ and Vℓ−1 using nested interpolation.
- 2. Local factorisation: spatially disjoint pieces of white noise are independent,
we can sample small independent local white noise vectors bℓ
e ∼ N(0, Mℓ e) on
each supermesh cell e and assemble the contributions together to obtain bℓ
S.
RESULT: Linear cost complexity and trivially parallelisable!
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 22
Mathematical Institute University of Oxford Overview
Introduction and Background Multilevel Monte Carlo Multilevel Quasi Monte Carlo Conclusions
SLIDE 23
Mathematical Institute University of Oxford Related work
Herrmann and Schwab (2017 pre-print) Same application problem. Some aspects of the theory are very general. Algorithm is detailed only for 1D, nested, structured grids. Truncated QMC. Our work Focussed on algorithm development for 2D and 3D, non-nested, unstructured grids. Non-truncated hybrid QMC/MC.
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 24
Mathematical Institute University of Oxford Extending work on MLMC to MLQMC: the challenge
Low-discrepancy sequences are extremely uniform in the first few dimensions and in low-dimensional projections, but less so across the whole hypercube. Therefore QMC works best when the integrand has low effective dimension (adapted from [Joe and Kuo 2008]).
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 25 Mathematical Institute University of Oxford Extending work on MLMC to MLQMC: the challenge
Low-discrepancy sequences are extremely uniform in the first few dimensions and in low-dimensional projections, but less so across the whole hypercube. Therefore QMC works best when the integrand has low effective dimension (adapted from [Joe and Kuo 2008]). For good QMC convergence we need to order the dimensions in our QMC integrands in order of decaying importance so that the largest error components are
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 26 Mathematical Institute University of Oxford Haar wavelet expansion of white noise
From now on, D = [0, 1]d. We expand ˙ W into a Haar wavelet series. Let l ∈ ({−1} ∩ N)d, n ∈ Nd, |l| = maxi(li), x+ = max(x, 0), H−1,0(x) = 1D(x),
1 x −1 1 H0,0(x) 1 x −21/2 21/2 H1,0(x) 1 x −21/2 21/2 H1,1(x) 1 x −2 2 H2,0(x) 1 x −2 2 H2,1(x) 1 x −2 2 H2,2(x) 1 x −2 2 H2,3(x)
˙ W =
|l|=∞
(2l −1)+
zl,n(ω)Hl,n(x), zl,n ∼ N(0, 1) i.i.d.
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 27
Mathematical Institute University of Oxford Multilevel Quasi Monte Carlo white noise sampling
The H0,0 = H0,0(x)H0,0(y) 2D Haar wavelet and the level L = |l| = 1 Haar mesh.
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 28
Mathematical Institute University of Oxford Multilevel Quasi Monte Carlo white noise sampling
The H0,0 = H0,0(x)H0,0(y) 2D Haar wavelet and the level L = |l| = 1 Haar mesh.
Truncate the series of ˙ W at Haar level L . Let ˙ WL := truncation, ˙ WR := remainder. ˙ W = ˙ WL + ˙ WR. 1) Order all the NL = 2d(L +1) coefficients in ˙ WL according to |l|1. 2) Use hybrid MC/QMC sampling (figure).
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 29 Mathematical Institute University of Oxford Multilevel Quasi Monte Carlo white noise sampling
˙ WR is sampled with pseudo-random points so no ordering needed. The covariance
˙ WR is known and a similar technique as in MLMC can be used for the sampling.
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 30 Mathematical Institute University of Oxford Multilevel Quasi Monte Carlo white noise sampling
˙ WR is sampled with pseudo-random points so no ordering needed. The covariance
˙ WR is known and a similar technique as in MLMC can be used for the sampling. This time a 3-way supermesh is required. Let NS = # supermesh cells. The Haar mesh is “nice” so we expect NS ≤ C(d) (# cells of the finest parent mesh). Numerical experiments suggest C(d) < 3 in 1D and 2D and C(d) ≈ O(50) in 3D. Overall Mat´ ern field sampling cost is O(NL log(NL ) + NS), where the log term can be dropped if using locally supported Haar wavelets or if d = 1. L ≤ Lmax.
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 31
Mathematical Institute University of Oxford MC vs QMC vs MLQMC (2D, ν = 1, M = 256)
1 2 3 4 5 MLQMC level ℓ 10−12 10−11 10−10 10−9 10−8 log2(NℓVar)
MLQMC approx ν = 1.0 (2D)
Nℓ = 1 Nℓ = 16 Nℓ = 256 Nℓ = 4096
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 32
Mathematical Institute University of Oxford MC vs QMC vs MLQMC (2D, ν = 1, M = 256)
1 2 3 4 5 MLQMC level ℓ 10−12 10−11 10−10 10−9 10−8 log2(NℓVar)
MLQMC approx ν = 1.0 (2D)
Nℓ = 1 Nℓ = 16 Nℓ = 256 Nℓ = 4096
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 33
Mathematical Institute University of Oxford MLQMC convergence and cost (2D, ν = 1)
1 2 3 4 5
MLQMC level ℓ
100 101 102 103 104
Nℓ ε = 5.000e − 05 ε = 2.500e − 05 ε = 1.250e − 05 ε = 6.250e − 06 ε = 3.125e − 06 ε = 1.563e − 06 ε = 7.813e − 07
10−6 10−5
ε
10−4 10−3
ε2Cost Std QMC MLMC MLQMC
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 34
Mathematical Institute University of Oxford Overview
Introduction and Background Multilevel Monte Carlo Multilevel Quasi Monte Carlo Conclusions
SLIDE 35 Mathematical Institute University of Oxford Conclusions Outlook
Ordering the variables is needed for good QMC convergence. Asymptotic QMC convergence rate is still O(N−1/2). However, large gains to be found in the pre-asymptotic regime while keeping the QMC dimension contained. We can couple white noise between different FEM approximation subspaces. A supermesh construction is not needed in the nested case. 2-way supermesh needed for MLMC and 3-way supermesh needed for MLQMC. The overall order of complexity is linear in the number of elements of the supermesh and it can be trivially parallelised. Standard techniques usually have cubic complexity. Software used: FEniCS, libsupermesh, Intel R
- MKL. Fast FEniCS-based Mat´
ern field and white noise sampling software to appear on Bitbucket soon! JUQ paper (MLMC only) and slides: https://croci.github.io MLQMC paper and PhD thesis to be submitted soon, please read both!
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
SLIDE 36 Mathematical Institute University of Oxford Selected references - Thank you for listening!
JUQ paper [1] and slides: https://croci.github.io
[1] M. Croci, M. B. Giles, M. E. Rognes, and P. E. Farrell. Efficient white noise sampling and coupling for multilevel Monte Carlo with non-nested meshes. SIAM/ASA Journal on Uncertainty Quantification, 6(4):1630–1655, 2018. doi: 10.1137/18M1175239. [2] P. E. Farrell, M. D. Piggott, C. C. Pain, G. J. Gorman, and C. R. Wilson. Conservative interpolation between unstructured meshes via supermesh construction. Computer Methods in Applied Mechanics and Engineering, 198:2632–2642, 2009. [3] M. B. Giles. Multilevel Monte Carlo path simulation. Operations Research, 56(3):607–617, 2008. doi: 10.1287/opre.1070.0496. [4] M. B. Giles and B. J. Waterhouse. Multilevel quasi-Monte Carlo path simulation. Advanced Financial Modelling, Radon Series on Computational and Applied Mathematics, (8):165–181, 2009. [5] L. Herrmann and C. Schwab. Multilevel quasi-Monte Carlo integration with product weights for elliptic PDEs with lognormal coefficients Multilevel quasi-Monte Carlo integration with product weights for elliptic PDEs with lognormal coefficients. Technical report, ETH Zurich, Zurich, 2017. [6] F. Y. Kuo, R. Scheichl, C. Schwab, I. H. Sloan, and E. Ullmann. Multilevel Quasi-Monte Carlo Methods for Lognormal Diffusion Problems. Mathematics of Computation, 86(308):2827–2860, 2017. [7] F. Lindgren, H. Rue, and J. Lindstr¨
- m. An explicit link between Gaussian fields and Gaussian
Markov random fields: the stochastic partial differential equation approach. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(4):423–498, 2009. [8] A. J. Wathen. Realistic eigenvalue bounds for the Galerkin mass matrix. IMA Journal of Numerical Analysis, 7(August 1985):449–457, 1987. doi: 10.1093/imanum/7.4.449.
EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling
