Computing integrals in many dimensions Whats new? Ian H. Sloan - PowerPoint PPT Presentation

Computing integrals in many dimensions – What’s new? Ian H. Sloan i.sloan@unsw.edu.au The University of New South Wales – UNSW Australia Shanghai Jiao Tong University, 2015

Can we think of evaluating a 100 -dimensional integral, � 1 � 1 I ( F ) := . . . F ( y 1 , . . . , y 100 )d y 1 · · · d y 100 ? 0 0 Some scientists and engineers want to evaluate integrals in 100 or 1 , 000 or even 1 million dimensions.

How? We cannot use a product rule! Even a 2 -point Gauss rule in each of 100 dimensions would need 2 100 function evaluations! In general we want to approximate (with maybe large s ) � 1 � 1 I s ( F ) := . . . F ( y 1 , . . . , y s )d y 1 . . . d y s 0 0 � = [0 , 1] s F (y)dy THE WORKHORSE IS MONTE CARLO

Monte Carlo (MC) N N,s ( F ) := 1 � Q MC F (t k ) , N k =1 with t 1 , . . . , t N chosen randomly and independently from a uniform distribution on [0 , 1] s . Error: For F ∈ L 2 ([0 , 1] s ) , error MC ( F ) = σ ( F ) √ , N where σ 2 ( F ) = I s (( F − I s ( F )) 2 ) = I s ( F 2 ) − ( I s ( F )) 2 .

What sort of problems? Need to be guided by applications. Artificial problem are often too easy (or too hard)!

Uncertainty is a major source of applications PDE with random coefficients mathematical finance

Uncertainty as a main source of applications PDE with random coefficients mathematical finance In 1995 Traub and Paskov computed values of ‘mortgage-backed obligations’. In the USA people buying a house borrow money from a bank. In return, the bank holds a ’mortgage’. USA mortgages last for 30 years, and may be repaid each month, making 30 × 12 = 360 repayment possibilities. The computed quantity is a 360-dimensional expected value - hence a 360-dimensional integral. The 1995 experiments used a ‘Quasi Monte Carlo’ method. The surprising message from the finance experiments: For this problem, Quasi Monte Carlo beats Monte Carlo.

A 64 -dimensional option pricing example Here is a result of a Monte Carlo calculation in 64 dimensions for a problem of pricing a digital Asian option: 0.52 0.51 0.50 N 2 12 2 13 2 14

A 64 -dimensional option pricing example And here’s the result of a QMC calculation for the same 64 - dimensional problem of pricing a digital Asian option: 0.52 0.51 0.005 0.50 N N 2 12 2 13 2 14 2 12 2 13 2 14

Quasi-Monte Carlo (QMC) For QMC we take N I s ( F ) ≈ Q N,s ( F ) := 1 � F (t k ) , N k =1 with t 1 , . . . , t N deterministic (and cleverly chosen). How to choose t 1 , . . . , t N ? Low discrepancy points (Halton, Sobol ′ , Faure, Niederreiter, . . . ) Lattice rules

Low discrepancy points and lattice rules Algorithm: N I s ( F ) ≈ 1 � F (t k ) N k =1 t k ∈ [0 , 1] s • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • First 64 points of A lattice rule Monte Carlo method 2D Sobol ′ sequence with 64 points with 64 “random” points

Lattice rules Lattice rules were introduced by number theorists in the late 1950s and 1960s – especially Korobov, Hlawka, Zaremba, L K Hua They were designed for periodic functions. But most integrands are not periodic.

We consider only the simplest kind of lattice rule , given by N − 1 Q N,s (z; F ) = 1 k z �� �� � F , N N k =0 where z ∈ { 1 , . . . , N − 1 } s , and the braces mean that each component of the s -vector in the braces is to be replaced by its fractional part

Example of lattice & shifted lattice rules N = 34 , z = (1 , 21) N = 34 , z = (1 , 21) , ∆ = (0 . 8 , 0 . 1) 1 1 • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • 0 0 • • 0 1 0 1

We consider only the simplest kind of lattice rule, given by N − 1 Q N,s (z; F ) = 1 k z �� �� � F , N N k =0 where z ∈ { 1 , . . . , N − 1 } s , and the braces mean that each component of the s -vector in the braces is to be replaced by its fractional part. The corresponding shifted lattice rule is N − 1 ∆; F ) = 1 k z �� �� � Q N,s (z , ∆ ∆ F N + ∆ ∆ ∆ , N k =0 ∆ ∈ [0 , 1) s . ∆ where ∆

What’s new? For an important class of problems, we now know how to DESIGN lattice rules that is, we know how to find a good z in N − 1 Q N,s (z; F ) = 1 k z �� �� � F , N N k =0 to obtain a guaranteed rate of convergence, independently of dimension. The problem class is: PDE with random coefficients

PDE with random coefficients PDE with random coefficients are now attracting great interest. Example: flow through a porous medium Darcy’s law is q (x) = − a (x) ∇ p (x) , � where p = p (x) is pressure of the fluid q = � � q (x) is velocity of the fluid a = a (x) is “permeability” of the medium Incompressibility: ∇ · � q = 0 Together these give a second order elliptic PDE: ∇ · ( a (x) ∇ p (x)) = 0

Modelling the permeability Describing in all the microscopic pores and channels in a real material is commonly considered much too hard. So it is common engineering practice to model the permeability as a random field:

Example – the flow problem For each realisation of a (x) we need to solve the PDE plus boundary conditions. A 2-dimensional example is: x 2 ∂p ∂n = 0 1 p = 1 p = 0 x 1 0 1 ∂p ∂n = 0 We have Dirichlet boundary conditions on Γ D , the left and right edges, and Neumann boundary conditions on Γ N , the top and bottom edges.

Finite element method Standard piecewise-linear finite elements: x 2 ∂p ∂n = 0 1 p = 1 p = 0 x 1 0 1 ∂p ∂n = 0

Particle paths For a particular realisation of the permeability field, after we have found the pressure field p , to find the position x of a particle of the fluid solve dx d t = � q (x) = − a (x) ∇ p (x) , subject to x = (0 , 0 . 5) at t = 0 . x 2 ∂p ∂n = 0 1 p = 1 place a particle → ◦ p = 0 • x 1 0 1 ∂p ∂n = 0

Particle displacement after time t = 0.1 1 0.8 0.6 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 Note: For each realisation of the field the permeability field is different!

A simple model problem – the “uniform” case −∇ · ( a (x , y) ∇ u (x , y)) = f (x) in D , y ∈ U := [0 , 1] N , u (x , y) = 0 ∂D, on with D a bounded Lipschitz domain in R d , and ∞ � ( y j − 1 a (x , y) = a (x) + 2 ) ψ j (x) , x ∈ D, y ∈ U , j =1 where y 1 , y 2 , . . . are independent random variables uniformly with ψ j such that � distributed on [0 , 1] ; j � ψ j � ∞ < ∞ , and with a large enough and � j � ψ j � ∞ small enough to ensure a max ≥ a (x , y) ≥ a min > 0 , making the PDE strongly elliptic for every y . In practice we must truncate the sum after s terms, so the random part of the problem is

Many approaches: polynomial chaos, generalized polynomial chaos, stochastic Galerkin, stochastic collocation Monte Carlo multilevel Monte Carlo Quasi-Monte Carlo multilevel Quasi-Monte Carlo All methods face serious challenges when the effective dimensionality is high. And when all else fails, people turn to Monte Carlo methods. We aim to beat Monte Carlo.

Many contributors: Many important results using various methods by Babuska, Nobile, Tempone Harbrecht, Peters, Siebenmorgen Babuska, Tempone, Zouraris Hoang, Schwab Barth, Schwab, Zollinger Karniadakis, Xiu Charrier Kunoth, Schwab Nobile, Tempone, Webster Charrier, Scheichl, Teckentrup Cliffe, Giles, Scheichl, Teckentrup Schwab, Todor Cliffe, Graham, Scheichl, Stals Schillings, Schwab Cohen, Chkifa, Schwab Teckentrup, Scheichl, Giles, Ullmann Webster Cohen, De Vore, Schwab . . Graham, Scheichl, Ullmann . Hansen, Schwab

Application of QMC to PDE with random coefficients Graham, Kuo, Nuyens, Scheichl, Sloan (J. Comput. Physics, 2011) - Application of QMC to the lognormal case, no error analysis - Use circulant embedding to avoid truncation of KL expansion Kuo, Schwab, Sloan (SIAM J. Numer. Anal., 2012) - Application of QMC to the uniform case [cf. Cohen, De Vore, Schwab, 2010] - First complete error analysis: how to choose “weights” Kuo, Schwab, Sloan (J. FoCM, to appear) [cf. Heinrich 1998; Giles 2008] - A multi-level version of the analysis for the uniform case Schwab (Proceedings of MCQMC 2012) - Application of QMC to the uniform case for affine operator equations Le Gia (Proceedings of MCQMC 2012) - Application of QMC to the uniform case for PDE over the sphere