[PPT] - Statistical Inversion for Basal Parameters for the Antarctic Ice PowerPoint Presentation

SLIDE 1

Statistical Inversion for Basal Parameters for the Antarctic Ice Sheet

Tobin Isaac, Noémi Petra, Georg Stadler, Omar Ghattas1

Institute for Computational Engineering & Sciences The University of Texas at Austin

SIAM UQ 2014 Savannah, Georgia April 3, 2014

1Departments of Geological Sciences and Mechanical Engineering Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 1 / 28

SLIDE 2

1

Forward Model: Ice Sheet Dynamics

2

Bayesian Inversion via methods of Constrained Optimization

3

Example: Pine Island Glacier

4

Example: Antarctic Ice Sheet These slides can be found at www.ices.utexas.edu/~tisaac/slides.

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 2 / 28

SLIDE 3

Motivation

Observed surface flow velocity [Rignot et al., 2011] Antarctic ice sheet inversion for the basal friction parameter field using InSAR surface velocity measurements

IPCC Fourth Assessment Report [Meehl et al., 2007]

[T]he uncertainty in the projections of the land ice contributions [to sea level rise] is dominated by the various uncertainties in the land ice models themselves ... rather than in the temperature projections.’’

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 3 / 28

SLIDE 4

Equations

Balance of linear momentum, mass, and energy −∇ · [η(θ, u) ˙ ε − Ip] = ρg, [˙ ε = 1

2(∇u + ∇uT)]

∇ · u = 0, ρc ∂θ ∂t + u · ∇θ

− ∇ · (K∇θ) = 2 η tr(˙

ε2) Constitutive relations η(θ, u) =

A0 exp
− Q

Rθ

− 1

n

˙ ε

1−n 2n

II

[˙ εII = 1

2tr(˙

ε2)] Basal boundary conditions

T(σn) + β(x)Tu = 0

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 4 / 28

SLIDE 5

Discretization

High-order finite element discretization: Qp × Qp−2

3 1 4 2

complex geometry high-order accurate LBB stable element-wise mass conservative inf-sup constant independent of h and weakly dependent on p inf-sup constant unaffected by high aspect ratio ice sheet geometry amenable to scalable parallel AMR schemes better efficiency than low-order (increased flops per memop), if there is an efficient solver.

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 5 / 28

SLIDE 6

Solution

Nonlinear solver is inexact Newton-Krylov method: Newton’s method: linearized Stokes equations have 4th-order tensor coefficient, as do true adjoint equations. Inexact Krylov tolerances [Eisenstat and Walker, 1996]: minimizes Krylov iterations in Newton steps not believed to be in asymptotic regime. Linear solver: trilinear submesh preconditioner for high-order elements [Brown, 2010] Smoothed Aggregation AMG: tolerates heterogeneities/anisotropies of viscous operator better than geometric multigrid ASM(i)/ILU(j) smoother for high-aspect-ratio geometries

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 6 / 28

SLIDE 7

Antarctic simulation: scaling

Newton-Krylov Solver for Antarctic Ice Sheet Stokes Equations

50 100 150 200

Krylov iterations

10−13 10−10 10−7 10−4 10−1

rk/r0

p = 3 [18] p = 4 [23] p = 5 [23]

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 7 / 28

SLIDE 8

Antarctic simulation: scaling

3 6 9 12 15 10−14 10−12 10−10 10−8 10−6 10−4 10−2 100

Newton iterations (Krylov iterations labeled) ||rk||/||r0|| Newton-Krylov Solver for Antarctic Ice Sheet Stokes Equations

bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC bC

28M dofs, 1024 cores 122M dofs, 4096 cores 122M dofs, 16384 cores

bC bC bC

30 4 7 14 21 54 8 3 8 82 104 33 5 7 17 33 38 6 5 13 75 143 35 5 7 16 40 42 4 3 3 4 4 93 196 356

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 7 / 28

SLIDE 9

Inventory

We have Scalable iterative solution of nonlinear forward problem Scalable iterative solution of Newton linearization

Scalable iterative solution of adjoint equations Time to solution is more-or-less linear in the tolerance Algebraic multilevel representation of linear operator

We do not currently have Full multilevel approximation scheme of the nonlinear problem

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 8 / 28

SLIDE 10

Bayesian-inversion-in-a-slide

Let f be the map from the parameter m (controlling the basal coefficient β) to the observations through the forward problem: f (m) = B(u), A(u, m) = 0 We assume additive Gaussian noise in the measurements: d = f (m) + e, e ∼ N(0, Γ

noise)

Thus:

πlike(d|m) = exp

− 1

2(f (m) − d)TΓ−1

noise(f (m) − d)

If the prior is Gaussian with mean mprior and covariance Γ

prior, then we obtain for

the posterior pdf:

πpost(m) ∝ exp

− 1

2 f (m) − d 2

Γ−1

noise − 1

2 m − mprior 2

Γ−1

prior

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin)

Inversion for Ice Sheet Parameters SIAM UQ 2014 9 / 28

SLIDE 11

Posterior statistic 0: MAP point

The maximum a posteriori point is mMAP

def

= arg max

m

πpost(m) = arg min

m

− log πpost(m) [

def

= J(m)] = arg min

m,u 1 2 B(u) − d 2

Γ−1

noise + 1

2 m − mprior 2

Γ−1

prior

such that A(u, m) = 0. This is classical PDE-constrained optimization: 1

2Γ−1

prior characterizes the

regularization.

Γ−1

prior is (should be) given a priori by expert knowledge, not chosen

a posteriori when the solution has desirable regularity (although if ‘‘expert prior knowledge’’ is ‘‘I know it when I see it’’, there is little difference)

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 10 / 28

SLIDE 12

Scalable MAP point computation via RSINK

The KKT first-order necessary conditions (J′(m) = 0) of a solution m are

A(u, m) = 0

[state],

AT

u(u, m)p = B(u) − d

[adjoint],

AT

m(u, m)p + Γ−1

prior m = 0

[gradient], where u and p are state and adjoint variables. We solve this using Reduced Space Inexact Newton-Krylov: Reduced space: search for a minimum of J(m); given mk, uk and pk satisfy state and adjoint equations. Inexact Newton-Krylov: like inexact Newton-Krylov used to solve state equations

J′′(mk)δk = −J′(mk)

mk+1 = mk + αδk, where α is determined by line search.

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 11 / 28

SLIDE 13

Necessary components of scalable, efficient RSINK implementation

Scalable, efficient cost evaluation ⇔ Scalable, efficient nonlinear forward solver Scalable, efficient gradient computation Scalable, efficient (Gauss-)Newton Hessian application

H(m)

def

= J′′(m).

The Hessian, or an approximation, is useful for UQ: please attend Noémi’s talk in MS7 three days ago for more details.

Scalable, efficient (Gauss-)Newton Hessian preconditioner

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 12 / 28

SLIDE 14

Adjoint-methods-in-a-slide

If uk and pk solve state and adjoint equations (uk computed from mk, pk computed form uk and mk), then

J′(m) = pT

kAm(uk, mk) + mT kΓ−1

prior

All terms in equations have finite element expressions (not differentiating matrices).

Marginal cost of J′(m), given J(m): one linear PDE solve.

Linearizing the first-order KKT system of equations results in set of PDEs for computing Hessian action

Computing H(mk) ˜ m costs two (‘‘incremental’’) linear PDE solves with the same operators as Newton step / adjoint. Another opportunity to exploit inexactness: loose tolerance on incrementals approximates H(mk), but quadratic/superlinear convergence of inexact Newton only requires accurate Hessian near the solution [Hicken and Alonso, 2013].

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 13 / 28

SLIDE 15

Prior preconditioning the Hessian

We precondition the full Hessian H = H

misfit + Γ−1 prior with Γ prior, so the

preconditioned system is

H

misfit Γ prior + I ∼ Γ

1/2

prior H misfit Γ

1/2

prior + I.

The justification from PDE-constrained optimization: The precision Γ−1

prior is a differential operator.

The prior Γ

prior is compact (low effective rank), so H misfit Γ prior is compact.

Identity + compact ⇒ good (h-independent) Krylov method convergence. The implication for statistical inversion: Covariance Γ

post of Gaussianized posterior distribution is H(mMAP)−1.

The effective numerical rank of H

misfit Γ prior is the dimension of the subspace

f parameters informed by the data.

The more information the data contains, the less compact is H

misfit Γ prior Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 14 / 28

SLIDE 16

Posterior statistic 1: MAP covariance

To compute an estimate of Γ

post, we compute a low-rank approximate

eigendecomposition of Γ

1/2

prior H misfitΓ

1/2

prior . There are two main approaches:

Lanczos methods, randomized SVD [Halko et al., 2011], combinations (e.g. block Lanczos). Note that in both cases, it is not necessary to have Γ

1/2

prior , as one can use

Lanczos on H

misfit preconditioned by Γ prior, or

a randomized Γ

prior-orthogonal SVD of H misfit.

Note that randomized methods are inherently more parallel than Lanczos methods.

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 15 / 28

SLIDE 17

Sampling the Gaussianized posterior distribution

However we do our eigendecomposition, we end up with an approximation

f Γ

post of the form

Γ

post ≈ Γ prior − WrDrW T

r ,

where W T

r Γ−1

prior Wr = Ir.

How to generate a sample z ∼ N(0, Γ

post)? Given x ∼ N(0, Γ prior),

z = (I + Wr((Ir − Dr)1/2 − Ir)W T

r Γ−1

prior )x

Depending on its form, sampling Γ

prior may itself be non-trivial:

requires access to some symmetric factorization Γ

prior = GGT (such as Γ1/2 prior )

r a good approximation thereof.

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 16 / 28

SLIDE 18

Inventory

Scalable, efficient nonlinear forward and adjoint solvers buy us

cost and gradient evaluation, and Hessian application.

Assuming the information gained from the data is low rank:

a good approximation of Γ

prior buys us fast, h-independent solution of systems

involving the reduced Hessian H(m); drawing samples from N(0, Γ

post) is as easy as drawing samples from

N(0, Γ

prior). Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 17 / 28

SLIDE 19

Bayesian inference for the basal friction field: InSAR data (Pine Island region)

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 18 / 28

SLIDE 20

Bayesian inference for the basal friction field: InSAR data (Pine Island region)

Left: InSAR-based Antarctica ice surface velocity observations Right: Inferred basal friction field (MAP point)

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 19 / 28

SLIDE 21

Bayesian inference for the basal friction field: InSAR data (Pine Island region)

Left: Recovered ice surface velocity observations Right: Inferred basal friction field (MAP point)

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 19 / 28

SLIDE 22

Bayesian inference for the basal friction field: InSAR data (Pine Island region)

Left: Prior pointwise standard deviation Right: Posterior pointwise standard deviation

Isaac, Petra, Stadler, Ghattas (ICES, UT Austin) Inversion for Ice Sheet Parameters SIAM UQ 2014 20 / 28

SLIDE 23