[PPT] - Parameter Estimation and Uncertainty Quantifjcation for Flow and PowerPoint Presentation

SLIDE 1

Parameter Estimation and Uncertainty Quantifjcation for Flow and Transport Processes

Ole Klein

WIAM 2016, Hamburg

02.09.2016

Ole Klein PE and UQ for Flow and Transport 02.09.2016 1 / 26

SLIDE 2

Introduction

Motivation

PDE-based modeling often leads to:

High-dimensional discretized parameterization
Sparse (but not necessarily low-dimensional) state
bservations

Here: Groundwater fmow equation, transport equation, Richards equation Parameter estimation should:

Produce estimates of parameters and uncertainty
Be applicable for high-resolution parameter fjelds
Remain feasible for large numbers of state
bservations

Ole Klein PE and UQ for Flow and Transport 02.09.2016 2 / 26

SLIDE 3

Introduction

Forward and Inverse Problem

S U P Z I F O G P: Parameter vectors S: Parameter fjelds U: System states Z: Measurement vectors I: Input, Interpretation F: Forward model (continuous) O: Output, Observation (sparse) G: Forward model (discrete) Task: construct mapping Z → P to estimate parameters

Ole Klein PE and UQ for Flow and Transport 02.09.2016 3 / 26

SLIDE 4

Introduction

Example: Groundwater Flow

Groundwater fmow equation: Fφ(Y, Zs, φ0; φ) := Ss(Zs)∂tφ + ∇ · jθw(Y, φ) − qθw = 0 jθw(Y, φ) := −K(Y)∇φ = − exp(Y)∇φ Convection-difgusion equation: Fc(Y, Zs, φ0, c0, φ; c) := θ∂tc + ∇ · jC(Y, φ, c) − qC = 0 jC(Y, φ, c) := −D(φ)∇c + cjθw Parameters P: log-storage term Zs, log-conductivity Y, initial conditions Measurements Z: hydraulic head φ(xi, ti), tracer concentration c(xj, tj) Parameters are underdetermined, problem is ill-posed and requires regularization

Ole Klein PE and UQ for Flow and Transport 02.09.2016 4 / 26

SLIDE 5

Bayesian Inversion

Prior Distribution

Treat parameters as random variables: P ∼ N (QPP, P∗) In more detail:      p1 p2 . . . pnP      ∼ N           Qp1p1 Qp1p2 · · · Qp1pnP Qp2p1 Qp2p2 · · · Qp2pnP . . . . . . . . . QpnPp1 QpnPp2 · · · QpnPpnP      ,      p∗

1

p∗

2

. . . p∗

nP

         

Stationary random fjelds → constituent matrices are block Toeplitz
Multiplication with QPP using circulant embedding
Multiplication with Q−1

PP is problematic

Ole Klein PE and UQ for Flow and Transport 02.09.2016 5 / 26

SLIDE 6

Bayesian Inversion

Maximum A Posteriori

Assume that measurement noise ǫ is Gaussian: ǫ = [Z − G(P)] ∼ N (QZZ, 0) Z|P ∼ N (QZZ, G(P)) Apply Bayes’ theorem fZ · fP|Z = fP · fZ|P and drop constant term: fP|Z ∝ fP · fZ|P Use maximizer Pmap of fP|Z as parameter estimate

Ole Klein PE and UQ for Flow and Transport 02.09.2016 6 / 26

SLIDE 7

Bayesian Inversion

Objective Function

For Gaussian parameters and Gaussian noise, we have: fP|Z ∝ exp ( −1 2 [ P − P∗2

Q−1

PP + Z − G(P)2

Q−1

ZZ

]) Finding Pmap is equivalent to minimizing objective function L(P) := 1 2 P − P∗2

Q−1

PP + 1

2 Z − G(P)2

Q−1

ZZ

with gradient ∇L = Q−1

PP [P − P∗] − HT ZPQ−1 ZZ [Z − G(P)]

(HZP: linearization of G around current estimate)

Ole Klein PE and UQ for Flow and Transport 02.09.2016 7 / 26

SLIDE 8

Preconditioned Conjugate Gradients

Prior Preconditioned Conjugate Gradients

Apply Conjugate Gradients to objective function expressed in Q−1/2

PP P

Equivalent to using Q−1

PP as preconditioner:

δP = −QPP∇L|Pi−1[+ conj. contrib.] = − [Pi−1 − P∗] + QPPHT

ZPQ−1 ZZ [Z − G(P)]

Allows using combined adjoint (full assembly of HZP not needed)
Completely eliminates Q−1

PP from algorithm (negative preconditioner cost)

Low-rank perturbation of identity (resulting convergence rate independent of

mesh width)

Ole Klein PE and UQ for Flow and Transport 02.09.2016 8 / 26

SLIDE 9

Preconditioned Conjugate Gradients

Prior Preconditioned Conjugate Gradients

20 40 60 80 100 10−2 10−1 100 Number of Iterations Normalized Objective Function Conjugate Gradients 20 40 60 80 100 10−2 10−1 100 Number of Iterations Normalized Objective Function Preconditioned Conjugate Gradients Convergence behavior for 2D mesh of size 64 × 64 ( ), 128 × 128 ( ), 256 × 256 ( ) and 512 × 512 ( )

Ole Klein PE and UQ for Flow and Transport 02.09.2016 9 / 26

SLIDE 10

Preconditioned Conjugate Gradients

Prior Preconditioned Conjugate Gradients

104 105 102 104 number of cells nP ttotal [s] number of measurements nZ = 16 500 500 1,000 number of measurements nZ ttotal [s] number of parameters nP = 256 × 256 Time to solution for prior preconditioned CG ( ), Gauss-Newton reference ( ) and unpreconditioned CG ( )

Ole Klein PE and UQ for Flow and Transport 02.09.2016 10 / 26

SLIDE 11

Preconditioned Conjugate Gradients

Example: 2D Groundwater Flow

Synthetic Reference Estimate for nφ = 25 Estimate for nφ = 900

Ole Klein PE and UQ for Flow and Transport 02.09.2016 11 / 26

SLIDE 12

Uncertainty Quantifjcation and Analysis

Uncertainty Quantifjcation

Provide linearized uncertainty quantifjcation through posterior covariance matrix: Qpost

PP :=

[ Q−1

PP + HT ZPQ−1 ZZHZP

]−1 Use transformation P → Q−1/2

PP P again, equivalent to split

Qpost

PP = Q1/2 PP [I + Mlike]−1 Q1/2 PP

Mlike := Q1/2

PPHT ZPQ−1 ZZHZPQ1/2 PP

Mlike has low rank, use spectral decomposition Mlike = VΛVT: Qpost

PP = Q1/2 PP

[ I + VΛVT]−1 Q1/2

PP

= Q1/2

PP

[ I − VΥVT] Q1/2

PP

(with Υ = diag(λi/(λi + 1))) = QPP − Q1/2

PPVΥVTQ1/2 PP,

Approximating Q1/2

PP through circulant embedding introduces systematic error!

Ole Klein PE and UQ for Flow and Transport 02.09.2016 12 / 26

SLIDE 13

Uncertainty Quantifjcation and Analysis

Example: 2D Groundwater Flow

Uncertainty for nφ = 25 Uncertainty for nφ = 900

Ole Klein PE and UQ for Flow and Transport 02.09.2016 13 / 26

SLIDE 14

Uncertainty Quantifjcation and Analysis

Quality Assurance / Sample Generation

Results of previous steps: Pmap and Qpost

PP

Assumption: posterior distribution can be modeled as P|Z

approx.

∼ N ( Qpost

PP, Pmap

)

Pmap is mode of posterior distribution, not its mean
Qpost

PP is linearization of true posterior variance

Higher-order terms are neglected
Construct has to be checked for consistency with data set
Realistic applications require samples from posterior distribution

Ole Klein PE and UQ for Flow and Transport 02.09.2016 14 / 26

SLIDE 15

Uncertainty Quantifjcation and Analysis

Quality Assurance / Sample Generation

Instead of spectral decomposition Mlike = VΛVT, use Qpost

PP = Q1/2 PP

[ I + LlikeLT

like

]−1 Q1/2

PP

Llike := Q1/2

PPHZPQ−1/2 ZZ

Performing singular value decomposition (SVD) Llike = VPΛ1/2VT

Z leads to

Qpost

PP = QPP − Q1/2 PPVPΥVT PQ1/2 PP

= Q1/2

PPV+ P

[ I − Υ+] [ V+

P

]T Q1/2

PP

Qpost

ZZ = QZZ − Q1/2 ZZVZΥVT ZQ1/2 ZZ

= Q1/2

ZZV+ Z

[ I − Υ+] [ V+

Z

]T Q1/2

ZZ

This can be used to check distribution of estimation error (if known) or measurement residuals, and generate conditional samples: Qpost

PP = LPLT P with LP := Q1/2 PPV+ P

[ I − Υ+]1/2 [ V+

P

]T Qpost

ZZ = LZLT Z with LZ := Q1/2 ZZV+ Z

[ I − Υ+]1/2 [ V+

Z

]T

Ole Klein PE and UQ for Flow and Transport 02.09.2016 15 / 26

SLIDE 16

Examples

Example: 2D Groundwater Flow

−20 −10 10 20 50 100 Value Count Residual of PCG, tol = 10−4 −4 −2 2 4 2,000 4,000 Value Count Error of PCG, tol = 10−4 Normalized residuals and errors for 400 ( ), 625 ( ), 900 ( ) measurements Optimization stopped early, only normalized residuals refmect this

Ole Klein PE and UQ for Flow and Transport 02.09.2016 16 / 26

SLIDE 17

Examples

Example: 2D Groundwater Flow

−4 −2 2 4 20 40 60 80 Value Count Residual of PCG, tol = 10−5 −4 −2 2 4 2,000 4,000 Value Count Error of PCG, tol = 10−5 Normalized residuals and errors for 400 ( ), 625 ( ), 900 ( ) measurements Optimization has converged, both measures indicate linearization is appropriate

Ole Klein PE and UQ for Flow and Transport 02.09.2016 17 / 26

SLIDE 18

Examples

Example: 3D Groundwater Flow

Synthetic Reference (Interior) Estimate for nφ = 225 (Interior) Estimation of 3D log-conductivity parameter fjeld 128 × 128 × 16 = 2.65 · 105 parameters, 225 measurements Moments of normalized residual: (−0.195, 0.960, 0.029, 2.99)

Ole Klein PE and UQ for Flow and Transport 02.09.2016 18 / 26

SLIDE 19

Examples

Example: 3D Groundwater Flow

Synthetic Reference (Interior) Estimate for nφ = 225 (Interior) Estimation of 3D log-conductivity parameter fjeld 128 × 128 × 16 = 2.65 · 105 parameters, 225 measurements Moments of normalized residual: (−0.195, 0.960, 0.029, 2.99)

Ole Klein PE and UQ for Flow and Transport 02.09.2016 19 / 26

SLIDE 20

Examples

Example: Convective Transport

Infmow Boundary Flow Direction Measurement Locations Outfmow Boundary

experimental setup log-conductivity Y transport of tracer c

Domain of size nΩ = 256 × 256
1.31 · 105 parameters, 1.68 · 104 transient measurements, r = 250
Initial condition of tracer is Gaussian random fjeld for simplicity

Ole Klein PE and UQ for Flow and Transport 02.09.2016 20 / 26

SLIDE 21

Examples

Example: Convective Transport

Reference for Y Estimate for Y Mean Variance Skewness Kurtosis ∆P 0.202 9.396 1.439 8.287 ∆Z −0.827 192.5 2.229 34.63

Ole Klein PE and UQ for Flow and Transport 02.09.2016 21 / 26

SLIDE 22

Examples

Example: Convective Transport

Reference for c0 Estimate for c0 Mean Variance Skewness Kurtosis ∆P 0.202 9.396 1.439 8.287 ∆Z −0.827 192.5 2.229 34.63

Ole Klein PE and UQ for Flow and Transport 02.09.2016 22 / 26

SLIDE 23

Examples

Example: Convective Transport

Uncertainty of Estimate for Y Uncertainty of Estimate for c0 Mean Variance Skewness Kurtosis ∆P 0.202 9.396 1.439 8.287 ∆Z −0.827 192.5 2.229 34.63

Ole Klein PE and UQ for Flow and Transport 02.09.2016 23 / 26

SLIDE 24

Conclusions

Parameter estimation:

Cost of PCG has same dependence on nP as standard Gauss-Newton
Cost is largely independent of number of observations nZ
Caching variant eliminates multiplication with Q−1

PP and has negative

preconditioner cost Uncertainty quantifjcation:

Randomized methods provide linearized UQ
Automatically takes correlation of state observations into account
Introduces systematic error in multiplication with Q1/2

PP

Normalized errors / residuals:

Provide consistency checks for posterior distribution
No additional costs when using randomized SVD
Normalized residuals typically more sensitive than normalized errors

Ole Klein PE and UQ for Flow and Transport 02.09.2016 24 / 26

SLIDE 25

Appendix

Caching Preconditioned CG

In addition to iteration Pi, preconditioned residual Ti and direction Di, store Ri = Q−1

PPTi,

Vi := Q−1

PPPi,

Wi := Q−1

PPDi

Allows evaluations Ri = −∇L|Pi−1 = − [Vi−1 − V∗] + HZPQ−1

ZZ [Z − G(Pi−1)]

with V∗ := Q−1

PPP∗,

Pi−1 + αDi − P∗Q−1

PP = [Pi−1 + αDi − P∗]T [Vi−1 + αWi − V∗]

= PT

i−1Vi−1 + α2DT i Wi + [P∗]T V∗

+ 2αWT

i [Pi−1 − P∗] − 2VT i−1P∗

and also holds for conjugation factors βi.

Ole Klein PE and UQ for Flow and Transport 02.09.2016 25 / 26

SLIDE 26

Appendix

Randomized SVD

Goal: generate decomposition Llike ≈ VZ,rΛ1/2

r

VT

P,r

Produce r samples of measurement noise Rj and their images Tj = LlikeRj: Rj ∼ N (I, 0) Tj := LlikeRj = Q1/2

PPHT ZPQ−1/2 ZZ Rj

Create orthonormal basis C of range space and its image B: Bj := LT

likeCj = Q−1/2 ZZ HZPQ1/2 PPCj,

B is representation of LT

like on range space of Llike

Compute SVD and recover original right singular vectors: B = VZ,rΛ1/2

r

ST VP,r := CS

Ole Klein PE and UQ for Flow and Transport 02.09.2016 26 / 26