University of Oxford Efficient white noise sampling and coupling for - - PowerPoint PPT Presentation

Mathematical Institute University of Oxford Efficient white noise sampling and coupling for multilevel Monte Carlo

• M. Croci (Oxford), M. B. Giles (Oxford), P. E. Farrell (Oxford), M. E. Rognes (Simula)

MCQMC2018 - July 4, 2018

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford Overview

Introduction White noise sampling Numerical results Conclusions and further work

Mathematical Institute University of Oxford Overview

Introduction White noise sampling Numerical results Conclusions and further work

Mathematical Institute University of Oxford Motivation

The motivation of our research is the sampling of lognormal Gaussian fields. A Mat´ ern Gaussian field (approximately) satisfies a linear elliptic SPDE of the form Lu = ˙ W , x ∈ D, ω ∈ Ω + BCs, where u = u(x, ω) and ˙ W is spatial white noise. Other approaches can be used (with pros and cons), but we will not discuss them here.

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford Motivation

The motivation of our research is the sampling of lognormal Gaussian fields. A Mat´ ern Gaussian field (approximately) satisfies a linear elliptic SPDE of the form Lu = ˙ W , x ∈ D, ω ∈ Ω + BCs, where u = u(x, ω) and ˙ W is spatial white noise. Other approaches can be used (with pros and cons), but we will not discuss them here. The same techniques can be used to solve a more general class of SPDEs, e.g. N(u) + Lu = ˙ W , x ∈ D, ω ∈ Ω + BCs. In this case solving means to compute E[P(u)] for some functional P of the solution.

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford Motivation

The motivation of our research is the sampling of lognormal Gaussian fields. A Mat´ ern Gaussian field (approximately) satisfies a linear elliptic SPDE of the form Lu = ˙ W , x ∈ D, ω ∈ Ω + BCs, where u = u(x, ω) and ˙ W is spatial white noise. Other approaches can be used (with pros and cons), but we will not discuss them here. The same techniques can be used to solve a more general class of SPDEs, e.g. N(u) + Lu = ˙ W , x ∈ D, ω ∈ Ω + BCs. In this case solving means to compute E[P(u)] for some functional P of the solution. Common applications: finance, geology, meteorology, biology. . .

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford Motivation

The motivation of our research is the sampling of lognormal Gaussian fields. A Mat´ ern Gaussian field (approximately) satisfies a linear elliptic SPDE of the form Lu = ˙ W , x ∈ D, ω ∈ Ω + BCs, where u = u(x, ω) and ˙ W is spatial white noise. Other approaches can be used (with pros and cons), but we will not discuss them here. The same techniques can be used to solve a more general class of SPDEs, e.g. N(u) + Lu = ˙ W , x ∈ D, ω ∈ Ω + BCs. In this case solving means to compute E[P(u)] for some functional P of the solution. Common applications: finance, geology, meteorology, biology. . . Main issue: sampling ˙ W is hard!

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford Motivation

The motivation of our research is the sampling of lognormal Gaussian fields. A Mat´ ern Gaussian field (approximately) satisfies a linear elliptic SPDE of the form Lu = ˙ W , x ∈ D, ω ∈ Ω + BCs, where u = u(x, ω) and ˙ W is spatial white noise. Other approaches can be used (with pros and cons), but we will not discuss them here. The same techniques can be used to solve a more general class of SPDEs, e.g. N(u) + Lu = ˙ W , x ∈ D, ω ∈ Ω + BCs. In this case solving means to compute E[P(u)] for some functional P of the solution. Common applications: finance, geology, meteorology, biology. . . Main issue: sampling ˙ W is hard! The efficient sampling of ˙ W is the focus of this talk.

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford White noise (1D)

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

• 150
• 100
• 50

50 100 150

WARNING! Point evaluation not defined!

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford White noise (2D)

WARNING! Point evaluation not defined!

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford White noise (2D)

WARNING! Point evaluation not defined! IDEA! Avoid point evaluation by integrating ˙ W .

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford White Noise (practical definition) Definition (Spatial White Noise ˙ W )

For any φ ∈ L2(D), define ˙ W , φ :=

• D ˙

W φ dx. For any φi, φj ∈ L2(D), bi = ˙ W , φi, bj = ˙ W , φj are zero-mean Gaussian random variables, with, E[bibj] =

• D

φiφj dx =: Mij, b ∼ N(0, M). (1)

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford White Noise (practical definition)

IMPORTANT NOTE: Generalised random fields of type I and of type II.

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford Finite element (FEM) framework

When solving SPDEs (see 1st slide) with FEM, we get (for linear problems) Discrete weak form: find uh ∈ Vh s.t. for all vh ∈ Vh, a(uh, vh) = ˙ W , vh, (2) Where Vh = span({φi}n

i=0), (e.g. with Lagrange elements). EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford Finite element (FEM) framework

When solving SPDEs (see 1st slide) with FEM, we get (for linear problems) Discrete weak form: find uh ∈ Vh s.t. for all vh ∈ Vh, a(uh, vh) = ˙ W , vh, (2) Where Vh = span({φi}n

i=0), (e.g. with Lagrange elements).

FEM linear system: uh =

i uiφi, u = [u0, . . . , un]T,

Au = b(ω), (3) where the entries of b are given by, ˙ W , φi(ω) = bi(ω), (4) with b ∼ N(0, M) as before. M is the mass matrix of Vh.

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford Multilevel Monte Carlo [Giles 2008] + FEM framework

For MLMC, we have two approximation levels ℓ and ℓ − 1. For any particular ω ∈ Ω, we need to solve: find uℓ

h ∈ V ℓ h, uℓ−1 h

∈ V ℓ−1

h

s.t. for all vℓ

h ∈ V ℓ h, vℓ−1 h

∈ V ℓ−1

h

, a(uℓ

h, vℓ h) = ˙

W , vℓ

h(ω),

(5) a(uℓ−1

h

, vℓ−1

h

) = ˙ W , vℓ−1

h

(ω). (6)

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford Multilevel Monte Carlo [Giles 2008] + FEM framework

For MLMC, we have two approximation levels ℓ and ℓ − 1. For any particular ω ∈ Ω, we need to solve: find uℓ

h ∈ V ℓ h, uℓ−1 h

∈ V ℓ−1

h

s.t. for all vℓ

h ∈ V ℓ h, vℓ−1 h

∈ V ℓ−1

h

, a(uℓ

h, vℓ h) = ˙

W , vℓ

h(ω),

(5) a(uℓ−1

h

, vℓ−1

h

) = ˙ W , vℓ−1

h

(ω). (6) This yields the linear system

• Aℓ

Aℓ−1 uℓ uℓ−1

• =
• bℓ

bℓ−1

• = b,

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford Multilevel Monte Carlo [Giles 2008] + FEM framework

For MLMC, we have two approximation levels ℓ and ℓ − 1. For any particular ω ∈ Ω, we need to solve: find uℓ

h ∈ V ℓ h, uℓ−1 h

∈ V ℓ−1

h

s.t. for all vℓ

h ∈ V ℓ h, vℓ−1 h

∈ V ℓ−1

h

, a(uℓ

h, vℓ h) = ˙

W , vℓ

h(ω),

(5) a(uℓ−1

h

, vℓ−1

h

) = ˙ W , vℓ−1

h

(ω). (6) This yields the linear system

• Aℓ

Aℓ−1 uℓ uℓ−1

• =
• bℓ

bℓ−1

• = b,

where b ∼ N(0, M). Let V ℓ

h = span({φℓ i }nℓ i=0) and V ℓ−1 h

= span({φℓ−1

i

}nℓ−1

i=0 ), then

M =

• Mℓ

Mℓ,ℓ−1 (Mℓ,ℓ−1)T Mℓ−1

• ,

Mℓ,ℓ−1

ij

=

• φℓ

i φℓ−1 j

dx.

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford Multilevel Monte Carlo [Giles 2008] + FEM framework

For MLMC, we have two approximation levels ℓ and ℓ − 1. For any particular ω ∈ Ω, we need to solve: find uℓ

h ∈ V ℓ h, uℓ−1 h

∈ V ℓ−1

h

s.t. for all vℓ

h ∈ V ℓ h, vℓ−1 h

∈ V ℓ−1

h

, a(uℓ

h, vℓ h) = ˙

W , vℓ

h(ω),

(5) a(uℓ−1

h

, vℓ−1

h

) = ˙ W , vℓ−1

h

(ω). (6) This yields the linear system

• Aℓ

Aℓ−1 uℓ uℓ−1

• =
• bℓ

bℓ−1

• = b,

where b ∼ N(0, M). Let V ℓ

h = span({φℓ i }nℓ i=0) and V ℓ−1 h

= span({φℓ−1

i

}nℓ−1

i=0 ), then

M =

• Mℓ

Mℓ,ℓ−1 (Mℓ,ℓ−1)T Mℓ−1

• ,

Mℓ,ℓ−1

ij

=

• φℓ

i φℓ−1 j

dx. NOTE: we do not require the FEM approximation subspaces to be nested!

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford Overview

Introduction White noise sampling Numerical results Conclusions and further work

Mathematical Institute University of Oxford The challenge

SAMPLING PROBLEM 1: single level realisations: sample b ∼ N(0, M), where M is the mass matrix of Vh. SAMPLING PROBLEM 2: coupled realisations: sample b ∼ N(0, M), where M is the block mass matrix given by V ℓ

h and V ℓ−1 h

.

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford How to sample b?

Sampling b is hard!

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford How to sample b?

Sampling b is hard!

Na¨ ıve approach

• Factorise M = HHT (cubic complexity!) and set b = Hz, with z ∼ N(0, I).

⇒ E[bbT] = E[Hz(Hz)T] = HE[zzT]HT = HIHT = M.

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford How to sample b?

Sampling b is hard!

Na¨ ıve approach

• Factorise M = HHT (cubic complexity!) and set b = Hz, with z ∼ N(0, I).

⇒ E[bbT] = E[Hz(Hz)T] = HE[zzT]HT = HIHT = M.

• Works well if M diagonal: previous work used either mass-lumping [Lindgren, Rue

and Lindstr¨

• m 2009], piecewise constant elements [Osborn, Vassilevski and Villa

2017] or a piecewise constant approximation of white noise [Drzisga, et al. 2017, Du and Zhang 2002].

• We do not require M to be diagonal (and we do not approximate white noise).
• We can sample b with linear complexity.

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford How to sample b?

Sampling b is hard!

Na¨ ıve approach

• Factorise M = HHT (cubic complexity!) and set b = Hz, with z ∼ N(0, I).

⇒ E[bbT] = E[Hz(Hz)T] = HE[zzT]HT = HIHT = M.

• Works well if M diagonal: previous work used either mass-lumping [Lindgren, Rue

and Lindstr¨

• m 2009], piecewise constant elements [Osborn, Vassilevski and Villa

2017] or a piecewise constant approximation of white noise [Drzisga, et al. 2017, Du and Zhang 2002].

• We do not require M to be diagonal (and we do not approximate white noise).
• We can sample b with linear complexity.

IDEA! H does not need to be square, maybe we can find a more efficient factorisation!

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford White noise sampling: single level realisations

SAMPLING PROBLEM 1: need to sample b ∼ N(0, M). Exploit the FEM assembly M = LT     M1 · · · M2 ... . . . ... ...     L = LTdiage(Me)L. (7)

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford White noise sampling: single level realisations

SAMPLING PROBLEM 1: need to sample b ∼ N(0, M). Exploit the FEM assembly N(0, M) ∼ b = LT    b1 b2 . . .    = LTvstacke(be) (8)

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford White noise sampling: single level realisations

SAMPLING PROBLEM 1: need to sample b ∼ N(0, M). Exploit the FEM assembly

• Each be can be sampled as be = Heze with ze ∼ N(0, I) and HeHT

e = Me.

• b = LTvstacke(be) is N(0, M) since

E[bbT] = LTE[vstacke(be)vstacke(be)T]L = LTdiage(He)diage(HT

e )L = LTdiage(Me)L = M.

• If the mapping to the FEM reference element is affine (e.g. Lagrange elements on

simplices) we have that Me/|e| = const on each element and only one local factorisation is needed. This approach is trivially parallelisable!

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford White noise sampling: coupled realisations

SAMPLING PROBLEM 2: need to sample b ∼ N(0, M), where M is now the block mixed mass matrix.

Definition (Supermesh, [Farrell 2009])

Let A and B be two (possibly non-nested) meshes. Their supermesh S is one of their common refinements. A and B are both nested within S.

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford White noise sampling: coupled realisations

SAMPLING PROBLEM 2: need to sample b ∼ N(0, M), where M is now the block mixed mass matrix.

• Factorise locally, this time on each supermesh element.
• Sample b on S, then interpolate/project the result onto A and B (this step can be

performed locally).

• Since A and B are nested within S, this operation is exact. Note that A and B

need not be nested. Previous work on white noise coupling for MLMC used either a nested hierarchy [Drzisga et al. 2017, Osborn et al. 2017] or an algebraically constructed hierarchy of agglomerated meshes [Osborn, Vassilevski and Villa 2017].

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford Complexity overview

• ffline cost
• nline cost (per sample)

memory storage single level 0 (or O(m3N)) O(m3N) (or O(m2N)) O(m2) (or O(m2N)) single l. (affine) O(m3) O(m2N) O(m2) coupled 0 (or O(m3NS)) O(m3NS) (or O(m2NS)) O(m2) (or O(m2NS)) coupled (affine) O(m3) O(m2NS) O(m2)

Table: Memory and cost complexity of our white noise sampling strategy. In the non-affine case the cost per sample can be lowered by precomputing and storing the local factorisations (see entries in blue). NS is the number of supermesh elements. In our experience with MLMC, NS ≤ cdNℓ and cd = 2 (1D), cd = 2.5 (2D), cd = 45 (3D).

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford Overview

Introduction White noise sampling Numerical results Conclusions and further work

Mathematical Institute University of Oxford Numerical results: convergence of P(u) = u2

L2(D)

Consider the linear elliptic SPDE [Lindgren, Rue and Lindstr¨

• m 2009], [Bolin, Kirchner

and Kov´ acs 2017],

• I − κ−2∆

k u(x, ω) = η ˙ W , x ∈ D ⊆ Rd, ω ∈ Ω, ν = 2k − d/2 > 0. We compute FEM solutions {uℓ

h}ℓ=8 ℓ=1 with a non-nested hierarchy of subspaces {V ℓ h}ℓ=8 ℓ=1.

1 2 3 4 5 6 7 8

• 30
• 25
• 20
• 15
• 10
• 5

5 1 2 3 4 5 6 7 8

• 70
• 60
• 50
• 40
• 30
• 20
• 10

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford Numerical results: covariance convergence

C(r) = E[u(x)u(y)] = 1 2ν−1Γ(ν)(κr)νKν(κr), r = x − y2, κ = √ 8ν λ , x, y ∈ D,

0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford Overview

Introduction White noise sampling Numerical results Conclusions and further work

Mathematical Institute University of Oxford Conclusions and further work Outlook

• White noise is an extremely non-smooth object and is defined through its integral.
• We can sample single level white noise realisations efficiently.
• We can couple white noise between different FEM approximation subspaces. A

supermesh construction is not needed in the nested case.

• 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. Further work: extensions to QMC and MLQMC. Paper: https://arxiv.org/abs/1803.04857

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling

Mathematical Institute University of Oxford References - Thank you for listening!

[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. Preprint, 2018. URL https://arxiv.org/abs/1803.04857. [2] D. Bolin, K. Kirchner, and M. Kov´

• acs. Numerical solutions of fractional elliptic stochastic PDEs

with spatial white noise. Preprint, 2017. [3] D. Drzisga, B. Gmeiner, U. Ude, R. Scheichl, and B. Wohlmuth. Scheduling massively parallel multigrid for multilevel Monte Carlo methods. SIAM Journal of Scientific Computing, 39(5): 873–897, 2017. [4] Q. Du and T. Zhang. Numerical approximation of some linear stochastic partial differential equations driven by special additive noises. SIAM Journal on Numerical Analysis, 40(4):1421–1445, 2002. [5] M. B. Giles. Multilevel Monte Carlo path simulation. Operations Research, 56(3):607–617, 2008. [6] 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. [7] S. Osborn, P. S. Vassilevski, and U. Villa. A multilevel, hierarchical sampling technique for spatially correlated random fields. Preprint, 2017. [8] S. Osborn, P. Zulian, T. Benson, U. Villa, R. Krause, and P. S. Vassilevski. Scalable hierarchical PDE sampler for generating spatially correlated random fields using non-matching meshes. Preprint, 2017. [9] A. J. Wathen. Realistic eigenvalue bounds for the Galerkin mass matrix. IMA Journal of Numerical Analysis, 7:449–457, 1987.

EPSRC Centre for Doctoral Training in Industrially Focused Mathematical Modelling