Non-Intrusive Solution of Stochastic and Parametric Equations - - PowerPoint PPT Presentation

Non-Intrusive Solution of Stochastic and Parametric Equations

Hermann G. Matthiesa

Lo¨ ıc Giraldib, Alexander Litvinenkoc, Dishi Liud, and Anthony Nouyb

aInstitute of Scientific Computing, TU Braunschweig, Brunswick, Germany b´

Ecole Centrale de Nantes, GeM, Nantes, France

cKAUST, Thuwal, Saudi Arabia dInstitute of Aerodynamics and Flow Control, DLR, Brunswick, Germany

wire@tu-bs.de http://www.wire.tu-bs.de

13 Overview-Barcelona.tex,v 5.3.1.2 2015/01/06 00:23:34 hgm Exp

2

Overview

• 1. Parametric equations
• 2. Stochastic model problem
• 3. “Plain vanilla” Galerkin
• 4. To be or no to be intrusive
• 5. Numerical Comparison
• 6. Galerkin and low-rank tensor approximation
• 7. Non-intrusive computation
• 8. Numerical examples

TU Braunschweig Institute of Scientific Computing

3

General mathematical setup

Consider operator equation, physical system modelled by A: A(p; u) = f(p) u ∈ U, f ∈ F, U — space of states, F = U∗ — dual space of actions / forcings; Operator and rhs depend on parameters p, well posed for all p ∈ P. Iterative solver — convergent for all values of p — iterates for k = 0, . . . , u(k+1)(p) = S(p; u(k)(p), R(p; u(k)(p)), with u(k)(p) → u∗(p), where S is one cycle of the solver, and the residuum: R(u(k)) := R(p; u(k)(p)) := f(p) − A(p; u(k)). When the residuum vanishes — R(p; u∗(p)) = 0 — the mapping S has a fixed point u∗(p) = S(p; u∗(p), 0).

TU Braunschweig Institute of Scientific Computing

4

Model stochastic problem

Aquifier

0.5 1 1.5 2 0.5 1 1.5 2 Geometry flow out Dirichlet b.c. flow = 0 Sources

2D Model model with stochastic data, p −∇ · (κ(x, p)∇u(x, p)) = f(x, ω) & b.c., x ∈ G ⊂ Rd (κ(x, p)∇u(x, p)) · n = g(x, p), x ∈ Γ ⊂ ∂G, p ∈ P κ stochastic conductivity, f and g stochastic sinks and sources. One p is a realisation of κ, f, g.

TU Braunschweig Institute of Scientific Computing

5

Preconditioned residual

In the iteration set u(k+1) = u(k) + ∆u(k) with ∆u(k) := S(p; u(k), R(p; u(k))) − u(k), and usually P(∆u(k)) = R(p; u(k)), so that S(p; u(k)) = u(k) + P −1(R(p; u(k))). (list of arguments shortened) Here P is some preconditioner, which may depend on p, the iteration counter k, and on the current iterate u(k); e.g. in Newton’s method P = DuA(p; u(k)).

TU Braunschweig Institute of Scientific Computing

6

Iteration

Algorithm: Start with some initial guess u(0) k ← 0 while no convergence do Compute ∆u(k) ← S(p; u(k), R(p; u(k))) − u(k) u(k+1) ← u(k) + ∆u(k) k ← k + 1 end while Uniform contraction: ∀p, u, v : S(p; u(p), R(p; u(p))) − S(p; v(p), R(p; v(p)))U ≤ ̺∗ u(p) − v(p)U. with ̺∗ < 1.

TU Braunschweig Institute of Scientific Computing

7

Discretisation I

Let S ⊂ RP be an appropriate Hilbert space of real functions on P, look for solution in tensor space U := U ⊗ S so that u(p) =

ι uι ςι(p).

Normally, discretise first U by choosing finite-dimensional UN ⊂ U, but results here independent of that. Direct Integration: To compute Quantity of Interest (QoI) Q(u) =

• P

Q(u(p), p) µ(dp) ≈

• z

wzQ(u(pz), pz) integrand and u(pz) have to be computed for all pz: expensive! But decoupled, non-intrusive solves. Want to replace it by proxy / meta model or emulator u(p) ≈ uM(p)

TU Braunschweig Institute of Scientific Computing

8

Discretisation II

(Further) discretise U ⊗ S by choosing SM = span{Ψα(p)} ⊂ S to give U ⊗ SM = UM ⊂ U ⊗ S = U . Ansatz uM(ω) =

α uαΨα(p) ∈ UM,

Often Ψα(p) = Ψα(θ(p)), where θ(p) = [. . . , θℓ(p), . . . ] are independent. If Ψα(θ) =

ℓ ψαℓ(θℓ), where α = (. . . , αℓ, . . . ),

then SM =

ℓ SM,ℓ allows for higher degree tensors.

Simplest computation for Ψα to be orthogonal (orthonormal), e.g. in inner product φ, ϕS =

• P φ(p)ϕ(p) µ(dp)

How to determine the unknown coefficients uα ?

TU Braunschweig Institute of Scientific Computing

9

Solution procedures I

Projection of the solution u(p), or of the residuum R(p; u(p)). Interpolation: Determine the uα by interpolating condition: ∀pβ : u(pβ)

!

= uM(pβ) =

• α

uαΨα(pβ). Simplest when Kronecker-δ property Ψα(pβ) = δα,β satisfied. Solve equation on interpolation points pβ — decoupled, non-intrusive solves. Pseudo-spectral projection: Simple as Ψα are orthonormal. Compute projection inner product (integral) by quadrature, i.e. uα =

• P

Ψα(p)u(p) µ(dp) ≈

• z

wzΨα(pz)u(pz), solve equation on quadrature points pz — decoupled, non-intrusive solve.

TU Braunschweig Institute of Scientific Computing

10

Solution procedures II

Mapping u(·) → uM(·) is a projection Π, and to describe a general projection, choose ˆ SM = span{Φα(p)}, projection orthogonal to ˆ SM: ∀ϕ ∈ ˆ SM : (I − Π)u, ϕ = 0, i.e. ˆ S⊥

M = im(I − Π).

Approximation properties are determined by SM, stability by ˆ SM. Collocation / Interpolation i.e. solve equation on collocation / interpolation points pβ, i.e. Φβ(p) = δ(p − pβ): R(pβ; uM(pβ)) = R(pβ;

• α

uαΨα(pβ))

!

= 0. With Kronecker-δ: R(pβ; uβ) = 0 — same as interpolation, decoupled, non-intrusive solve. We worry about the norm Π. Norm of collocation projector ΠC may grow with M.

TU Braunschweig Institute of Scientific Computing

11

Projectors

Pseudo-spectral projector ΠP is orthogonal, i.e. ΠP = 1. This means that ˆ SM = SM, normally Φα = Ψα Galerkin: Apply Galerkin weighting. ∀β : Φβ(p), R(p; uM(p)) = Φβ(p), R(p;

• α

uαΨα(p)) = 0. Coupled equations, is it intrusive? When solved in a partitioned way, residua computed by quadrature, it is non-intrusive, needs only residua on qudrature points. To have norm of projector as small as possible (Bubnov-Galerkin), choose orthogonal projection Φα = Ψα,

TU Braunschweig Institute of Scientific Computing

12

Galerkin on iteration equation

Trick: Project iteration equation. Set u(k)(p) =

α u(k) α Ψα(p)) and

u(k) = [. . . , u(k)

β , . . . ]:

u(k+1) = u(k) + ∆u(k) = S(u(k), R(u(k))) = ⇒ u(k+1) = u(k) + ∆Mu(k), with ∆Mu(k) := [. . . , Ψβ, S(p; u(k)(p), R(p; u(k)(p))), . . . ] − u(k) Define a mapping S(u): S(u) := [. . . , Ψβ, S(p;

• α

uαΨα(p), R(p;

• α

uαΨα(p))), . . . ], then ∆Mu(k) = S(u(k)) − u(k) and u(k+1) = u(k) + ∆Mu(k) = S(u(k)).

TU Braunschweig Institute of Scientific Computing

13

Convergence

Nonlinear block Jacobi algorithm: Start with some initial guess u(0) k ← 0 while no convergence do Compute ∆Mu(k) as above u(k+1) ← u(k) + ∆Mu(k) k ← k + 1 end while Theorem: The mapping S has the same contraction factor ̺∗. This means that the simple nonlinear block Jacobi algorithm converges as before.

TU Braunschweig Institute of Scientific Computing

14

The myth about intrusiveness

Folklore: Galerkin methods are intrusive. They can be, but don’t have to. Question: To be or not to be intrusive? Stochastic Galerkin conditions for iteration equation requires S(u(k)), approximated by S(u(k)) ≈ SZ(u(k)) =

• z

υz Ψβ(pz)S

• pz, u(k)(pz), R(p; u(k)(pz))
• .

to give ∆Mu(k) ≈ ∆Zu(k) = SZ(u(k)) − u(k). This requires the evaluation of preconditioned residuum — one iteration with the solver — at each pz. Theorem still holds with ∆Mu(k) replaced by ∆Zu(k) in algorithm.

TU Braunschweig Institute of Scientific Computing

15

Numerical example

2 3 4 5 1 R R R R R R R R R 6

A(p; u) := (Ku + λ1(p1)(uTu) u) = λ2(p2)f 0 =: f(p). f 0 := [1, 0, 0, 0, 0]T.

TU Braunschweig Institute of Scientific Computing

16

Numerical example spec

Case 1 Case 2 Case 3 Case 4 λ1(p1) p1 + 2 p1 + 1.1 p1 + 2 sin(4p1 + 2) λ2(p2) p2 + 25 25p2 + 0.5 10p2 + 30 10 sin(p2) + 30 c.o.v. 2.5e−2 2.9e+1 1.7e−1 2.2e−1

10

−10

10

−8

10

−6

10

−4

10

−2

10 1e−08 1e−06 1e−04 1e−02

Convergence criteria (εtol)

RMSE 2nd order polynomial 3rd order polynomial 4th order polynomial 5th order polynomial

TU Braunschweig Institute of Scientific Computing

17

Numerical results

• rder

solver calls ǫ(L2(u)) ǫ(L1(u)) ǫ(L2(Ru)) m P G P G P G P G 2 79 90 6.1e-5 6.1e-5 3.5e-5 3.5e-5 4.1e-5 4.1e-5 3 161 192 3.9e-6 3.9e-6 2.3e-6 2.3e-6 2.6e-6 2.6e-6 4 284 325 2.7e-7 2.7e-7 1.6e-7 1.6e-7 1.8e-7 1.8e-7 5 458 540 2.0e-8 2.0e-8 1.2e-8 1.2e-8 1.4e-8 1.4e-8 Low rank approximation: write u := [. . . , uα, . . . ] = (uα,n) u =

α,n uα,neα ⊗ en ≈ ur = r−1 j=1 wj ⊗ ηj.

• Use faster global methods than block Jacobi, e.g. Quasi-Newton.
• Try and keep a low-rank tensor approximation troughout,

from input fields to output solution.

TU Braunschweig Institute of Scientific Computing

18

Successive rank-one updates (SR1U)

Assume a functional J(p; u), and that A(p; u) − f(p) = δuJ(p; u) = 0, so that solution is equivalent with minimising J for each p.

• Build solution rank-one by rank-one, i.e. with already computed

ur := r−1

j=1 wj ⊗ ηj add new term wr ⊗ ηr through:

min

wr,ηr J(ur + wr ⊗ ηr) ⇔ δw,ηJ(ur + wr ⊗ ηr) = 0

⇒ successive rank-one updates (SR1U), proper generalised decomposition (PGD).

• This Galerkin procedure only solves “small” problems, good

approximations often with small r.

TU Braunschweig Institute of Scientific Computing

19

Low-rank approximation (basic PGD)

Define Jr(wr, ηr) := J(ur(p) + wr ⊗ ηr). New wr and ηr found via system δwJr(wr, ηr) = 0 ∧ δηJr(wr, ηr) = 0 ∧ wr = 1. Block-Jacobi solver: u1 ← 0; η1 ← 1; w1 ← 0; for r = 1, . . . , until ur + wr ⊗ ηr “accurate enough” do : while no convergence do ηr ← ηr/ηr; Solve δwJr(wr, ηr) = 0 for wr; wr ← wr/wr; Solve δηJr(wr, ηr) = 0 for ηr; end while ur+1 ← ur + wr ⊗ ηr; end for Output: a basic (greedy) low-rank approximation ur.

TU Braunschweig Institute of Scientific Computing

20

Non-intrusive residual for PGD

Non-intrusive approximation of first equation: δwJr(wr, ηr) = 0 ⇔ δuJ(ur + wr ⊗ ηr), ηrS = 0 in U ⇔ 0 =

• P

R(p; ur(p) + wr ⊗ ηr(p))ηr(p) dp ≈

• z

υz R(pz; ur(pz) + wr ⊗ ηr(pz)) ηr(pz), 2ndeq.: δηJr(wr, ηr) = 0 ⇔ δuJ(ur + wr ⊗ ηr), wrU = 0 in S ⇔ ∀λ ∈ S : 0 =

• P

R(p; ur(p) + wr ⊗ ηr(p)), wrUλ(p) dp ≈

• z

υz R(pz; ur(pz) + wr ⊗ ηr(pz)), wrU λ(pz)

TU Braunschweig Institute of Scientific Computing

21

Recent improvements

• Increase ur by more than one term at a time (e.g. 5–10 terms)

— larger systems to be solved.

• Use faster algorithm than block-Jacobi,

e.g. Quasi-Newton methods (here BFGS).

• Use previous iterates as control variates to have

fewer integration points per iteration.

• Increase accuracy of integration (number of integration points)

as iteration converges.

TU Braunschweig Institute of Scientific Computing

22

PGD accuracy

d = 2 d = 3 d = 4 d = 5 Block-Jacobi solver [8] 5.14 × 10−5 3.31 × 10−6 2.31 × 10−7 1.70 × 10−8 Basic PGD (Algorithm 1) r = 1 2.34 × 10−3 2.34 × 10−3 2.34 × 10−3 2.34 × 10−3 r = 2 9.67 × 10−5 8.22 × 10−5 8.22 × 10−5 8.22 × 10−5 r = 3 5.14 × 10−5 3.39 × 10−6 8.03 × 10−7 7.78 × 10−7 r = 4 5.14 × 10−5 3.31 × 10−6 2.34 × 10−7 3.63 × 10−8 r = 5 5.14 × 10−5 3.31 × 10−6 2.31 × 10−7 1.71 × 10−8 Improved PGD (Algorithm 3) r = 1 2.34 × 10−3 2.34 × 10−3 2.34 × 10−3 2.34 × 10−3 r = 2 5.14 × 10−5 3.31 × 10−6 2.85 × 10−7 1.95 × 10−7 r = 3 5.14 × 10−5 3.31 × 10−6 2.31 × 10−7 1.79 × 10−8 r = 4 5.14 × 10−5 3.31 × 10−6 2.31 × 10−7 1.79 × 10−8 r = 5 5.14 × 10−5 3.31 × 10−6 2.31 × 10−7 1.76 × 10−8

TU Braunschweig Institute of Scientific Computing

23

PGD calls

− r = 1 r = 2 r = 3 r = 4 r = 5 Basic PGD (Algorithm 1) Relative error 2.34 × 10−3 8.22 × 10−5 7.78 × 10−7 3.63 × 10−8 1.71 × 10−8 Residual calls 1044 2160 3096 3816 4464 Improved algorithm (Algorithm 3) Relative error 2.34 × 10−3 1.95 × 10−7 1.79 × 10−8 1.79 × 10−8 1.79 × 10−8 Residual calls 1044 2304 2700 2844 3024

TU Braunschweig Institute of Scientific Computing

24

Obstacle example

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 x p g(p;x)

(a) Obstacle: g(p; x).

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 x p u(p;x)

(b) Solution: u(p; x).

Figure 2: Obstacle and solution as functions of and [3].

TU Braunschweig Institute of Scientific Computing

25

Obstacle example convergence

10 20 30 40 50 60 70 80 10−9 10−8 10−7 10−6 10−5 10−4 10−3 10−2 10−1 100 r Relative error SVD of the L2-projection Algorithm 1 Algorithm 3

TU Braunschweig Institute of Scientific Computing

26

Conclusion

• Parametric problems can be emulated.
• Galerkin methods can be non-intrusive.
• Convergence can be accelerated by faster global algorithms.
• For efficiency try and use sparse representation throughout; ansatz in

low-rank tensor products, saves storage as well as computation.

• PGD / SR1U is inherently a Galerkin procedure. Can also be

non-intrusive.

• Low-rank tensor representation can be very accurate.

TU Braunschweig Institute of Scientific Computing