[PPT] - Application of QMC methods to PDEs with random coefficients Frances PowerPoint Presentation

SLIDE 1

Frances Kuo

Application of QMC methods to PDEs with random coefficients

Frances Kuo

f.kuo@unsw.edu.au University of New South Wales, Sydney, Australia

Based on joint works with Ivan Graham (Bath), Rob Scheichl (Bath), Dirk Nuyens (KU Leuven), Christoph Schwab (ETH Zurich), Ian Sloan (UNSW), Josef Dick (UNSW), Thong Le Gia (UNSW), James Nichols (UNSW).

SLIDE 2

Frances Kuo

Motivating example

Uncertainty in groundwater flow

eg. risk analysis of radwaste disposal or CO2 sequestration

Darcy’s law q + a ∇p = f mass conservation law ∇ · q = 0 in D ⊂ Rd, d = 1, 2, 3 together with boundary conditions

[Cliffe, et. al. (2000)]

Uncertainty in a(x x x, ω) leads to uncertainty in q(x x x, ω) and p(x x x, ω)

SLIDE 3

Frances Kuo

Expected values of quantities of interest

To compute the expected value of some quantity of interest:

1. Generate a number of realizations of the random field

(Some approximation may be required)

2. For each realization, solve the PDE using e.g. FEM / FVM / mFEM
3. Take the average of all solutions from different realizations

This describes Monte Carlo simulation. Example: particle dispersion p = 1 p = 0

∂p ∂ n = 0 ∂p ∂ n = 0

release point →
0.2

0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1

SLIDE 4

Frances Kuo

Expected values of quantities of interest

To compute the expected value of some quantity of interest:

1. Generate a number of realizations of the random field

(Some approximation may be required)

2. For each realization, solve the PDE using e.g. FEM / FVM / mFEM
3. Take the average of all solutions from different realizations

This describes Monte Carlo simulation. NOTE: expected value = (high dimensional) integral s = stochastic dimension → use quasi-Monte Carlo methods

10

3

10

4

10

5

10

6

10

−5

10

−4

10

−3

10

−2

10

−1

error n n−1/2 n−1

r better

QMC MC

SLIDE 5

Frances Kuo

MC v.s. QMC in the unit cube

[0,1]s F (y

y y) dy y y ≈ 1 n

n

i=1

F (t t ti) Monte Carlo method Quasi-Monte Carlo methods t t ti random uniform t t ti deterministic n−1/2 convergence close to n−1 convergence (or better)

1 1 64 random points

1

1 First 64 points of a 2D Sobol′ sequence

1

1 A lattice rule with 64 points

more effective for earlier variables and lower-order projections
rder of variables irrelevant
rder of variables very important

use randomized QMC methods for error estimation

SLIDE 6

Frances Kuo

QMC

Two main families of QMC methods: (t,m,s)-nets and (t,s)-sequences lattice rules

1 1 First 64 points of a 2D Sobol′ sequence

1

1 A lattice rule with 64 points

(0,6,2)-net

Having the right number of points in various sub-cubes A group under addition modulo Z and includes the integer points

Niederreiter book (1992)

Sloan and Joe book (1994) Dick and Pillichshammer book (2010) K., Schwab, Sloan ANZIAM review (2011) Dick, K., Sloan Acta Numerica (2013)

Now also higher order digital nets

SLIDE 7

Frances Kuo

Lattice rules

Rank-1 lattice rules have points t t ti = frac i n z z z

,

i = 1, 2, . . . , n z z z ∈ Zs – the generating vector, with all components coprime to n frac(·) – means to take the fractional part of all components

1 1 A lattice rule with 64 points

SLIDE 8

Frances Kuo

Lattice rules

Rank-1 lattice rules have points t t ti = frac i n z z z

,

i = 1, 2, . . . , n z z z ∈ Zs – the generating vector, with all components coprime to n frac(·) – means to take the fractional part of all components

1 1 A lattice rule with 64 points 64 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63

∼ quality determined by the choice of z z z ∼

n = 64 z z z = (1, 19) t t ti = frac i 64 (1, 19)

“Component-by-component (CBC) construction”

[Korobov (1959); Sloan and Reztsov (2002);· · · ]

SLIDE 9

Frances Kuo

Randomly shifted lattice rules

Shifted rank-1 lattice rules have points t t ti = frac i n z z z + ∆ ∆ ∆

,

i = 1, 2, . . . , n ∆ ∆ ∆ ∈ [0, 1)s – the shift

shifted by ∆ ∆ ∆ = (0.1, 0.3) 1 1 A lattice rule with 64 points

1

1 A shifted lattice rule with 64 points

SLIDE 10

Frances Kuo

Randomly shifted lattice rules

Shifted rank-1 lattice rules have points t t ti = frac i n z z z + ∆ ∆ ∆

,

i = 1, 2, . . . , n ∆ ∆ ∆ ∈ [0, 1)s – the shift

shifted by ∆ ∆ ∆ = (0.1, 0.3) 1 1 A lattice rule with 64 points

1

1 A shifted lattice rule with 64 points

∼ use a number of random shifts for error estimation ∼

SLIDE 11

Frances Kuo

Weighted spaces

[Sloan and Wo´ zniakowski (1998);· · · ]

IDEA: the variables are not equally important (cf. effective dimension) We assume that F belongs to a weighted Sobolev space, with norm F 2

γ γ γ =

u⊆{1,...,s}

1 γu

[0,1]|u|
∂|u|F

∂y y yu (y y yu; 0)

2

dy y yu

2s subsets “anchor” at 0 “weights” Mixed first derivatives are square integrable Small weight γu means that F depends weakly on the variables y y yu

SLIDE 12

Frances Kuo

Weighted spaces

[Sloan and Wo´ zniakowski (1998);· · · ]

IDEA: the variables are not equally important (cf. effective dimension) We assume that F belongs to a weighted Sobolev space, with norm F 2

γ γ γ =

u⊆{1,...,s}

1 γu

[0,1]|u|
∂|u|F

∂y y yu (y y yu; 0)

2

dy y yu

2s subsets “anchor” at 0 “weights” Mixed first derivatives are square integrable Small weight γu means that F depends weakly on the variables y y yu

Alternatively, we may use the “unanchored” norm F 2

γ γ γ =

u⊆{1,...,s}

1 γu

[0,1]|u|
[0,1]s−|u|

∂|u|F ∂y y yu (y y yu;y y y−u)dy y y−u

2

dy y yu

SLIDE 13

Frances Kuo

Weighted spaces

[Sloan and Wo´ zniakowski (1998);· · · ]

IDEA: the variables are not equally important (cf. effective dimension) We assume that F belongs to a weighted Sobolev space, with norm F 2

γ γ γ =

u⊆{1,...,s}

1 γu

[0,1]|u|
∂|u|F

∂y y yu (y y yu; 0)

2

dy y yu

2s subsets “anchor” at 0 “weights” Mixed first derivatives are square integrable Small weight γu means that F depends weakly on the variables y y yu

Alternatively, we may use the “unanchored” norm F 2

γ γ γ =

u⊆{1,...,s}

1 γu

[0,1]|u|
[0,1]s−|u|

∂|u|F ∂y y yu (y y yu;y y y−u)dy y y−u

2

dy y yu “product” weights “order dependent” weights “POD” (product and order dependent) weights γu =

j∈u γj

γu = Γ|u| γu = Γ|u|

j∈u γj

SLIDE 14

Frances Kuo

Worst case error and CBC construction

We analyze the worst case error ewor

γ γ γ

(t t t1, . . . ,t t tn) := sup

F γ

γ γ≤1

[0,1]s F (y

y y) dy y y − 1 n

n

i=1

F (t t ti)

Then

integration error ≤ ewor

γ γ γ

(t t t1, . . . ,t t tn) F γ

γ γ

[Larcher, Leobacher, Scheicher (2003)]

Choose the weights γ γ γ = {γu} to minimize an upper bound for the RHS Choose the points t t t1, . . . ,t t tn to reduce the worst case error

SLIDE 15

Frances Kuo

Worst case error and CBC construction

We analyze the worst case error ewor

γ γ γ

(t t t1, . . . ,t t tn) := sup

F γ

γ γ≤1

[0,1]s F (y

y y) dy y y − 1 n

n

i=1

F (t t ti)

Then

integration error ≤ ewor

γ γ γ

(t t t1, . . . ,t t tn) F γ

γ γ

[Larcher, Leobacher, Scheicher (2003)]

Choose the weights γ γ γ = {γu} to minimize an upper bound for the RHS Choose the points t t t1, . . . ,t t tn to reduce the worst case error Component-by-component (CBC) construction The generating vector for a randomly shifted lattice rule “can be” obtained by a CBC algorithm to achieve close to O(n−1) convergence, with implied constant independent of s under an appropriate condition on the weights γ γ γ.

SLIDE 16

Frances Kuo

Worst case error and CBC construction

We analyze the worst case error ewor

γ γ γ

(t t t1, . . . ,t t tn) := sup

F γ

γ γ≤1

[0,1]s F (y

y y) dy y y − 1 n

n

i=1

F (t t ti)

Then

integration error ≤ ewor

γ γ γ

(t t t1, . . . ,t t tn) F γ

γ γ

[Larcher, Leobacher, Scheicher (2003)]

Choose the weights γ γ γ = {γu} to minimize an upper bound for the RHS Choose the points t t t1, . . . ,t t tn to reduce the worst case error Component-by-component (CBC) construction The generating vector for a randomly shifted lattice rule “can be” obtained by a CBC algorithm to achieve close to O(n−1) convergence, with implied constant independent of s under an appropriate condition on the weights γ γ γ.

Fast CBC with product weights: O(n log n s) operations [Nuyens and Cools, 2006] Anchored space with non-product weights: need auxiliary weights that depends on s... Fast CBC with POD weights in unanchorded space: O(n log n s + n s2) operations [K., Schwab, Sloan, 2011]

SLIDE 17

Frances Kuo

The present state of play

... in the application of QMC to PDEs with random coefficients:

1. Graham, K., Nuyens, Scheichl, Sloan (J. Comput. Physics, 2011)

Application of QMC to the lognormal case Use circulant embedding to avoid truncation of KL expansion Detailed numerical experiments, but no error analysis

2. K., Schwab, Sloan (SIAM J. Numer. Anal., 2012)

Application of QMC to the uniform case No numerical results, but we gave a complete error analysis Matches the best n term result by Cohen, De Vore, Schwab (2010) For the first time we know precisely how to choose the weights

3. K., Schwab, Sloan (submitted)

A multi-level version of the analysis for the uniform case

4. Graham, K., Nichols, Scheichl, Schwab, Sloan (submitted)

Application of QMC to the lognormal case Detailed numerical experiments as well as complete error analysis

SLIDE 18

Frances Kuo

The present state of play (continued)

5. Schwab (MCQMC 2012, to appear)

Application of QMC to the uniform case for general operator equations

6. Le Gia (MCQMC 2012, to appear)

Application of QMC to the uniform case for PDE over the sphere

7. Dick, Le Gia, K., Nuyens, Schwab (submitted)

Application of higher order QMC to the uniform case for general

perator equations

Beats the best n term result by Cohen, De Vore, Schwab (2010) New Banach (non-Hilbert) space setting New SPOD weights(smoothness-driven product and order dependent) New fast CBC construction of interlaced polynomial lattice rules

SLIDE 19

Frances Kuo

The uniform case

[K., Schwab, Sloan (2012)]

−∇ · (a(x x x,y y y) ∇u(x x x,y y y)) = f(x x x) in D, u(x x x,y y y) = 0 on ∂D a(x x x,y y y) = ¯ a(x x x) +

∞

j=1

yj ψj(x x x), y y y ∈ (− 1

2, 1 2)∞

Assumptions: There exists amin and amax such that 0 < amin ≤ a(x x x,y y y) ≤ amax There exists p ∈ (0, 1] such that ∞

j=1 ψjp L∞(D) < ∞

· · ·

SLIDE 20

Frances Kuo

The uniform case

[K., Schwab, Sloan (2012)]

−∇ · (a(x x x,y y y) ∇u(x x x,y y y)) = f(x x x) in D, u(x x x,y y y) = 0 on ∂D a(x x x,y y y) = ¯ a(x x x) +

∞

j=1

yj ψj(x x x), y y y ∈ (− 1

2, 1 2)∞

Assumptions: There exists amin and amax such that 0 < amin ≤ a(x x x,y y y) ≤ amax There exists p ∈ (0, 1] such that ∞

j=1 ψjp L∞(D) < ∞

· · · Goal: for F (y y y) = G(u(·,y y y)) a linear functional of u, we want to approximate

(− 1

2 , 1 2 )∞ F (y

y y) dy y y = lim

s→∞

(− 1

2 , 1 2 )s F (y1, . . . , ys, 0, 0, . . .) dy1 · · · dys

Our strategy: Truncate the infinite sum to s terms (truncation error) Solve the PDE using FEM (discretization error) Approximate the integral using QMC (quadrature error)

SLIDE 21

Frances Kuo

Error analysis – a sum of three terms

I(G(u)) ≈ Qs,n(G(us

h))

= 1 n

n

i=1

G(us

h(·,t

t ti))

∞

j=1 ψjp L∞(D) < ∞

f ∈ H−1+t(D) G ∈ H−1+t′(D)

Truncation error = O(s−2( 1

p −1))

use “integration”

Discretization error = O(ht+t′)

Aubin-Nitsche duality trick

SLIDE 22

Frances Kuo

Error analysis – a sum of three terms

I(G(u)) ≈ Qs,n(G(us

h))

= 1 n

n

i=1

G(us

h(·,t

t ti))

∞

j=1 ψjp L∞(D) < ∞

f ∈ H−1+t(D) G ∈ H−1+t′(D)

Truncation error = O(s−2( 1

p −1))

use “integration”

Discretization error = O(ht+t′)

Aubin-Nitsche duality trick

Quadrature error = O(n− min( 1

p − 1 2 ,1−δ))

Cf. best n term: n−( 1

p − 1 2 )

Differentiate the PDE w.r.t. y y y ∂ν

ν ν y y yu(·,y

y y)H1

0(D) ≤ C |ν

ν ν|!

j≥1 ψjνj L∞(D)

SLIDE 23

Frances Kuo

Error analysis – a sum of three terms

I(G(u)) ≈ Qs,n(G(us

h))

= 1 n

n

i=1

G(us

h(·,t

t ti))

∞

j=1 ψjp L∞(D) < ∞

f ∈ H−1+t(D) G ∈ H−1+t′(D)

Truncation error = O(s−2( 1

p −1))

use “integration”

Discretization error = O(ht+t′)

Aubin-Nitsche duality trick

Quadrature error = O(n− min( 1

p − 1 2 ,1−δ))

Cf. best n term: n−( 1

p − 1 2 )

Differentiate the PDE w.r.t. y y y ∂ν

ν ν y y yu(·,y

y y)H1

0(D) ≤ C |ν

ν ν|!

j≥1 ψjνj L∞(D)

Weighted Sobolev space with square-integrable mixed first derivatives Fast CBC construction of randomly shifted lattice rules Get POD weights (product and order dependent) γu =

c|u||u|!

j∈u ψjL∞(D)

b

SLIDE 24

Frances Kuo

POD weights arise naturally

For randomly shifted lattice rules we need to minimize, for λ ∈ (1/2, 1],   2 n

u⊆{1,...,s}

γu

λα|u| λ

 

1/λ

bound on worst case error (CBC)
u⊆{1,...,s}

(|u|!)2

j∈u ψj2 L∞(D)

γu

bound on norm

This is minimized by taking γu =

(|u|!)2

j∈u ψj2 L∞(D)

α|u|

λ

1/(1+λ) We get “POD” weights (“product and order dependent” weights).

SLIDE 25

Frances Kuo

POD weights arise naturally

For randomly shifted lattice rules we need to minimize, for λ ∈ (1/2, 1],   2 n

u⊆{1,...,s}

γu

λα|u| λ

 

1/λ

bound on worst case error (CBC)
u⊆{1,...,s}

(|u|!)2

j∈u ψj2 L∞(D)

γu

bound on norm

This is minimized by taking γu =

(|u|!)2

j∈u ψj2 L∞(D)

α|u|

λ

1/(1+λ) We get “POD” weights (“product and order dependent” weights).

Fast CBC with POD weights in unanchorded space: O(n log n s + n s2) operations Close to O(n−1) convergence when λ is close to 1/2 But require λ ≥ p 2 − p for the weights to be summable

SLIDE 26

Frances Kuo

Multi-level version

[K., Schwab, Sloan (submitted)]

I(G(u)) ≈ QL

∗ (G(u))

=

L

ℓ=0

Qsℓ,nℓ

G(usℓ

hℓ − u sℓ−1 hℓ−1)

cost = O

L

ℓ=0

sℓnℓhℓ

−d

SLIDE 27

Frances Kuo

Multi-level version

[K., Schwab, Sloan (submitted)]

I(G(u)) ≈ QL

∗ (G(u))

=

L

ℓ=0

Qsℓ,nℓ

G(usℓ

hℓ − u sℓ−1 hℓ−1)

cost = O

L

ℓ=0

sℓnℓhℓ

−d

Key estimate (challenge):

G(usℓ

hℓ − u sℓ−1 hℓ−1)γ γ γ ≤ G(usℓ hℓ − usℓ hℓ−1)γ γ γ + G(usℓ hℓ−1 − u sℓ−1 hℓ−1)γ γ γ

SLIDE 28

Frances Kuo

Multi-level version

[K., Schwab, Sloan (submitted)]

I(G(u)) ≈ QL

∗ (G(u))

=

L

ℓ=0

Qsℓ,nℓ

G(usℓ

hℓ − u sℓ−1 hℓ−1)

cost = O

L

ℓ=0

sℓnℓhℓ

−d

Key estimate (challenge):

G(usℓ

hℓ − u sℓ−1 hℓ−1)γ γ γ ≤ G(usℓ hℓ − usℓ hℓ−1)γ γ γ + G(usℓ hℓ−1 − u sℓ−1 hℓ−1)γ γ γ

RMS error with randomly shifted lattice rules: ht+t′

L

+ s−2(1/p−1)

L

+

L

ℓ=0

n− min(1/q−1/2,1−δ)

ℓ

ht+t′

ℓ−1 + s−(1/p−1) ℓ−1

New Assumption: there exists q ∈ [p, 1] such that

∞

j=1 ψjq W 1,∞(D) < ∞

differentiate PDE w.r.t. y

y y

need stronger regularity w.r.t. x

x x

use duality argument

SLIDE 29

Frances Kuo

Multi-level version

[K., Schwab, Sloan (submitted)]

I(G(u)) ≈ QL

∗ (G(u))

=

L

ℓ=0

Qsℓ,nℓ

G(usℓ

hℓ − u sℓ−1 hℓ−1)

cost = O

L

ℓ=0

sℓnℓhℓ

−d

Key estimate (challenge):

G(usℓ

hℓ − u sℓ−1 hℓ−1)γ γ γ ≤ G(usℓ hℓ − usℓ hℓ−1)γ γ γ + G(usℓ hℓ−1 − u sℓ−1 hℓ−1)γ γ γ

RMS error with randomly shifted lattice rules: ht+t′

L

+ s−2(1/p−1)

L

+

L

ℓ=0

n− min(1/q−1/2,1−δ)

ℓ

ht+t′

ℓ−1 + s−(1/p−1) ℓ−1

New Assumption: there exists q ∈ [p, 1] such that

∞

j=1 ψjq W 1,∞(D) < ∞

differentiate PDE w.r.t. y

y y

need stronger regularity w.r.t. x

x x

use duality argument

Constants depend on γu, ψjW 1,∞(D), etc. Choose γu to minimize some upper bound ⇒ get different POD weights Lagrange multipliers ⇒ choose sℓ, nℓ, hℓ

SLIDE 30

Frances Kuo

Single level cost v.s. multi-level cost

With error = O(ε) and cost =log O(ε−a) : aSL = p 2 − p + 2λp + d t + t′ aML =              max(2λq,

d t+t′ ) if “k-orthogonality” holds,

2λq if

d t+t′ ≤ 2λq − p 1−p, p 2−2p + d 2(t+t′) + max(λq, d 2(t+t′))

therwise.

where λp =     

1 2−2δ

p ∈ (0, 2

3], p 2−p

p ∈ ( 2

3, 1),

λq =     

1 2−2δ

q ∈ (0, 2

3], q 2−q

q ∈ ( 2

3, 1).

e.g. d = 2, t = t′ = 1, q = 2/3, k-orthogonality holds: error = O(M −1

hL ) and cost =log O(M 1/(1−δ) hL

)

SLIDE 31

Frances Kuo

The lognormal case

[Graham, K., Nichols, Scheichl, Schwab, Sloan (submitted)]

−∇ · (a(x x x, ω) ∇u(x x x, ω)) = f(x x x) in D, u(x x x, ω) = 0 on ∂D Karhunen-Loève expansion a(x x x, ω) = a∗(x x x) + a0(x x x) exp  

∞

j=1

√µj ξj(x x x) Yj(ω)  

(µj, ξj) - eigenpairs

D

ρ(x x x − x x x′)ξj(x x x′) dx x x′ = µj ξj(x x x) Yj(ω) - independent standard normal random numbers from R

e.g. ρ(r) = σ2 exp(−r/λ) or Gaussian or Matérn Results depend on the decay of µj

SLIDE 32

Frances Kuo

The lognormal case

[Graham, K., Nichols, Scheichl, Schwab, Sloan (submitted)]

−∇ · (a(x x x, ω) ∇u(x x x, ω)) = f(x x x) in D, u(x x x, ω) = 0 on ∂D Karhunen-Loève expansion a(x x x, ω) = a∗(x x x) + a0(x x x) exp  

∞

j=1

√µj ξj(x x x) Yj(ω)  

(µj, ξj) - eigenpairs

D

ρ(x x x − x x x′)ξj(x x x′) dx x x′ = µj ξj(x x x) Yj(ω) - independent standard normal random numbers from R

e.g. ρ(r) = σ2 exp(−r/λ) or Gaussian or Matérn Results depend on the decay of µj Dimension truncation error and FE discretization error are known

e.g. Charrier, Scheichl, Teckentrup (2011) Hoang, Schwab (2011) Teckentrup, Scheichl, Giles, Ullmann (2012)

SLIDE 33

Frances Kuo

QMC for the lognormal case

I(G(us

h)) =

Rs G(us

h(·,y

y y))

s

j=1

φ(yj) dy y y =

[0,1]s G(us

h(·, Φ−1(v

v v))) dv v v

φ – pdf of standard normal Φ−1 – icdf of standard normal

Complications: Need to map the integral into the unit cube using Φ−1

SLIDE 34

Frances Kuo

QMC for the lognormal case

I(G(us

h)) =

Rs G(us

h(·,y

y y))

s

j=1

φ(yj) dy y y =

[0,1]s G(us

h(·, Φ−1(v

v v))) dv v v

φ – pdf of standard normal Φ−1 – icdf of standard normal

Complications: Need to map the integral into the unit cube using Φ−1 Transformed integrand does not belong to the “standard” function space Need a special function space setting over Rs F 2

γ γ γ =

u⊆{1,...,s}

1 γu

R|u|
∂|u|F

∂y y yu (y y yu; 0)

2

j∈u

ψ2(yj) dy y yu

K., Sloan, Wasilkowski, Waterhouse (2010)

Need lattice convergence results for POD weights and unanchored variant

Nichols and K. (submitted)

Need to differentiate PDE w.r.t. y y y ... harder for the lognormal case

SLIDE 35

Frances Kuo

Higher order QMC rules

Classical polynomial lattice rule

[Niederreiter, 1992]

n = bm points with prime b An irreducible polynomial with degree m A generating vector of s polynomials with degree < m Interlaced polynomial lattice rule

[Goda and Dick, 2012]

Interlacing factor α An irreducible polynomial with degree m A generating vector of αs polynomials with degree < m Digit interlacing function Dα : [0, 1)α → [0, 1)

(0.x11x12x13 · · · )b, (0.x21x22x23 · · · )b, . . . , (0.xα1xα2xα3 · · · )b becomes (0.x11x21 · · · xα1x12x22 · · · xα2x13x23 · · · xα3 · · · )b

SLIDE 36

Frances Kuo

Higher order QMC rules

Classical polynomial lattice rule

[Niederreiter, 1992]

n = bm points with prime b An irreducible polynomial with degree m A generating vector of s polynomials with degree < m Interlaced polynomial lattice rule

[Goda and Dick, 2012]

Interlacing factor α An irreducible polynomial with degree m A generating vector of αs polynomials with degree < m Digit interlacing function Dα : [0, 1)α → [0, 1)

(0.x11x12x13 · · · )b, (0.x21x22x23 · · · )b, . . . , (0.xα1xα2xα3 · · · )b becomes (0.x11x21 · · · xα1x12x22 · · · xα2x13x23 · · · xα3 · · · )b

Different: higher order polynomial lattice rule

[Baldeaux and Dick, 2009]

An irreducible polynomial with degree αm A generating vector of s polynomials with degree < αm

SLIDE 37

Frances Kuo

Weighted Banach space norm

[Dick, Le Gia, K., Nuyens, Schwab (submitted)]

F q,r =

u⊆{1:s}
1

γuq

v⊆u
τ

τ τ u\v∈ {1:α}|u\v|

[0,1]|v|
[0,1]s−|v| (∂

(α α αv,τ τ τ u\v,0 0) y y y

F )(y y y) dy y y{1:s}\v

q

dy y yv r/q1/r ↑ as q ↑ ↑ as r ↓

Combines Lq norm for functions and ℓr norm for vectors Hilbert space setting with q = r = 2

SLIDE 38

Frances Kuo

Weighted Banach space norm

[Dick, Le Gia, K., Nuyens, Schwab (submitted)]

F q,r =

u⊆{1:s}
1

γuq

v⊆u
τ

τ τ u\v∈ {1:α}|u\v|

[0,1]|v|
[0,1]s−|v| (∂

(α α αv,τ τ τ u\v,0 0) y y y

F )(y y y) dy y y{1:s}\v

q

dy y yv r/q1/r ↑ as q ↑ ↑ as r ↓

Combines Lq norm for functions and ℓr norm for vectors Hilbert space setting with q = r = 2 Worst case error bound for digital nets, with 1

r + 1 r′ = 1

sup

F q,r≤1

|Is(F ) − Qs,n(F )| ≤

∅=u⊆{1:s}
C|u| γu
k

k ku∈D∗

u

b−µα(k

k ku)

r′1/r′ =: er′(S)

SLIDE 39

Frances Kuo

Weighted Banach space norm

[Dick, Le Gia, K., Nuyens, Schwab (submitted)]

F q,r =

u⊆{1:s}
1

γuq

v⊆u
τ

τ τ u\v∈ {1:α}|u\v|

[0,1]|v|
[0,1]s−|v| (∂

(α α αv,τ τ τ u\v,0 0) y y y

F )(y y y) dy y y{1:s}\v

q

dy y yv r/q1/r ↑ as q ↑ ↑ as r ↓

Combines Lq norm for functions and ℓr norm for vectors Hilbert space setting with q = r = 2 Worst case error bound for digital nets, with 1

r + 1 r′ = 1

sup

F q,r≤1

|Is(F ) − Qs,n(F )| ≤

∅=u⊆{1:s}
C|u| γu
k

k ku∈D∗

u

b−µα(k

k ku)

r′1/r′ =: er′(S)

Take r′ = 1 and r = ∞, with any q. Choose γu to make F q,∞ ≤ 1. Get SPOD weights (smoothness-driven product and order dependent)

γu =

ν

ν νu∈{1:α}|u|

|ν ν νu|!

j∈u
c ψj

νj L∞(D)

SLIDE 40

Frances Kuo

Fast CBC with SPOD weights

Error criterion with r′ = 1

e1(S) =

∅=u⊆{1:s}

C|u| γu

k

k ku∈D∗

u

b−µα(k

k ku)

=

∅=v⊆{1:αs}

C|u(v)| γu(v)

ℓ

ℓ ℓv∈D∗

v

b−µα(Eα(ℓ

ℓ ℓv,0 0))

≤

∅=v⊆{1:αs}

C|u(v)| γu(v)bα(α−1)|u(v)|/2

ℓ

ℓ ℓv∈D∗

v

b−αµ1(ℓ

ℓ ℓv)

Et(q1, . . . , qt) :=

∅=v⊆{1:t}
γv
ℓ

ℓ ℓv∈D∗

v

b−αµ1(ℓ

ℓ ℓv)

SLIDE 41

Frances Kuo

Fast CBC with SPOD weights

Error criterion with r′ = 1

e1(S) =

∅=u⊆{1:s}

C|u| γu

k

k ku∈D∗

u

b−µα(k

k ku)

=

∅=v⊆{1:αs}

C|u(v)| γu(v)

ℓ

ℓ ℓv∈D∗

v

b−µα(Eα(ℓ

ℓ ℓv,0 0))

≤

∅=v⊆{1:αs}

C|u(v)| γu(v)bα(α−1)|u(v)|/2

ℓ

ℓ ℓv∈D∗

v

b−αµ1(ℓ

ℓ ℓv)

Et(q1, . . . , qt) :=

∅=v⊆{1:t}
γv
ℓ

ℓ ℓv∈D∗

v

b−αµ1(ℓ

ℓ ℓv)

Fast (using FFT) CBC construction with SPOD weights Work in blocks of α dimensions Store intermediate sums and products Different update procedures for “within block” and for “switching block” Cost O(αsn log n + α2s2n) operations

SLIDE 42

Frances Kuo

Higher order convergence

CBC error bound for PDE example

e1(S) ≤

2

n − 1

0=ν

ν ν∈{0:α}s

(|ν ν ν|!)λ

s

j=1

νj>0

c ψj

νj L∞(D)

λ 1/λ , 1 α < λ ≤ 1 <

2

n − 1

u⊂N

|u|<∞

|u|!
j∈u

βj λ1/λ ,

∞

j=1

βp

j < ∞ ⇐

⇒

∞

j=1

ψjp

L∞(D) < ∞

SLIDE 43

Frances Kuo

Higher order convergence

CBC error bound for PDE example

e1(S) ≤

2

n − 1

0=ν

ν ν∈{0:α}s

(|ν ν ν|!)λ

s

j=1

νj>0

c ψj

νj L∞(D)

λ 1/λ , 1 α < λ ≤ 1 <

2

n − 1

u⊂N

|u|<∞

|u|!
j∈u

βj λ1/λ ,

∞

j=1

βp

j < ∞ ⇐

⇒

∞

j=1

ψjp

L∞(D) < ∞

Choose λ = p and α = ⌊1/p⌋ + 1 Convergence rate O(n− 1

p ), with implied constant independent of s

Beats best n term result O(n−( 1

p − 1 2 ))

SLIDE 44

Frances Kuo

Higher order convergence

CBC error bound for PDE example

e1(S) ≤

2

n − 1

0=ν

ν ν∈{0:α}s

(|ν ν ν|!)λ

s

j=1

νj>0

c ψj

νj L∞(D)

λ 1/λ , 1 α < λ ≤ 1 <

2

n − 1

u⊂N

|u|<∞

|u|!
j∈u

βj λ1/λ ,

∞

j=1

βp

j < ∞ ⇐

⇒

∞

j=1

ψjp

L∞(D) < ∞

Choose λ = p and α = ⌊1/p⌋ + 1 Convergence rate O(n− 1

p ), with implied constant independent of s

Beats best n term result O(n−( 1

p − 1 2 ))

Interlacing factor α is at least 2 Fully deterministic

SLIDE 45

Frances Kuo

The present state of play

... in the application of QMC to PDEs with random coefficients:

1. Graham, K., Nuyens, Scheichl, Sloan (J. Comput. Physics, 2011)

Application of QMC to the lognormal case Use circulant embedding to avoid truncation of KL expansion Detailed numerical experiments, but no error analysis

2. K., Schwab, Sloan (SIAM J. Numer. Anal., 2012)

Application of QMC to the uniform case No numerical results, but we gave a complete error analysis Matches the best n term result by Cohen, De Vore, Schwab (2010) For the first time we know precisely how to choose the weights

3. K., Schwab, Sloan (submitted)

A multi-level version of the analysis for the uniform case

4. Graham, K., Nichols, Scheichl, Schwab, Sloan (submitted)

Application of QMC to the lognormal case Detailed numerical experiments as well as complete error analysis

SLIDE 46

Frances Kuo

Summary and extensions

5. Schwab (MCQMC 2012, to appear)

Application of QMC to the uniform case for general operator equations

6. Le Gia (MCQMC 2012, to appear)

Application of QMC to the uniform case for PDE over the sphere

7. Dick, Le Gia, K., Nuyens, Schwab (submitted)

Application of higher order QMC to the uniform case for general

perator equations

Beats the best n term result by Cohen, De Vore, Schwab (2010) New Banach (non-Hilbert) space setting New SPOD weights(smoothness-driven product and order dependent) New fast CBC construction of interlaced polynomial lattice rules

SLIDE 47

Frances Kuo

Summary and extensions

5. Schwab (MCQMC 2012, to appear)

Application of QMC to the uniform case for general operator equations

6. Le Gia (MCQMC 2012, to appear)

Application of QMC to the uniform case for PDE over the sphere

7. Dick, Le Gia, K., Nuyens, Schwab (submitted)

Application of higher order QMC to the uniform case for general

perator equations

Beats the best n term result by Cohen, De Vore, Schwab (2010) New Banach (non-Hilbert) space setting New SPOD weights(smoothness-driven product and order dependent) New fast CBC construction of interlaced polynomial lattice rules Extensions Galerkin, mixed FEM, higher order FEM Non-linear functional, higher moments Multi-level, lognormal, · · ·