Optimal Experimental Design for Large-Scale Bayesian Inverse - - PowerPoint PPT Presentation

Optimal Experimental Design for Large-Scale Bayesian Inverse Problems via Multi-PDE-Constrained Optimization

Umberto Villa1 and Omar Ghattas2

1Washington University, Dept of Electrical & Systems Engineering 2The Oden Institute for Computational Engineering and Sciences 2Departments of Geological Sciences & Mechanical Engineering 2The University of Texas at Austin

Motivation

• PDE models (and other) often depend on some uncertain input

parameters (fields), whose values cannot be directly measured

• Bayesian inversion seeks to estimate the parameters and quantify the

uncertainty in the inversion, given observations of the outputs

• But for many problems, data acquisition is very expensive; how do we

choose optimal experiments?

– what is measured, where, and how often; which experiments to run, etc.

• Optimal experimental design (OED): find optimal data acquisition

strategy to minimize the uncertainty in the inferred parameters or maximize the information gained about the model

• Significant challenges:

– the parameter to be inferred is often a spatially correlated field – leads to high dimensional parameter space after discretization – the forward model is often expensive to solve (PDEs) – OED problem may have combinatorial complexity

Bayesian inversion for basal friction field in Antarctica

InSAR-based ice surface velocity observations Inferred mean of basal friction field

Reconstructed ice surface velocity field (based on inferred mean of basal friction field) Inferred uncertainty in basal friction field (standard deviation of Laplace approx of posterior of log basal friction)

Joint work with Tobin Isaac (Georgia Tech), Noemi Petra (UC-Merced), Georg Stadler (NYU)

Bayesian global seismic inversion

Joint work with Tan Bui-Thanh (UT Austin), Carsten Burstedde (Bonn), Georg Stadler (NYU), Lucas Wilcox (Naval Postgraduate School)

Bayesian poroelastic inversion for management of subsurface fluid injection

Use observations (InSAR, GPS) of surface deformation induced by either wastewater injection or production of reservoir or aquifer, in addition to well pressure measurements, to infer subsurface permeability and elastic properties. Forward predict and then ultimately optimize injection or extraction processes to avoid induced seismicity or preserve reservoir integrity.

Joint work with Amal Alghamdi, Ann Chen, Marc Hesse, Georg Stadler

Calibration of reduced models using DNS “data”: Which DNS experiments to run?

Turbulence/reduced chemical mechanism model for turbulent combustion, calibrated using DNS data (w/ Bob Moser, Todd Oliver, Umberto Villa, Peng Chen, UT Austin) Inference of a mesoscopic nonlocal Cahn-Hilliard model parameters of block copolymers from ”DNS” molecular dynamics simulations and X-ray scattering experiments (w/Frank Alexander (Brookhaven), Tinsley Oden, Karen Willcox (UT Austin), Danial Faghihi (Buffalo))

Goals

Devise scalable & efficient algorithms for OED problems governed by Bayesian inverse problems that are in turn governed by forward PDE models with uncertain parameter fields [an outer loop around an outer loop]

• “Scalable” means the cost, measured in # of PDE solves,

is independent of the parameter (m), data (y), and experimental design (w) dimensions

• Efficient means the constant is “small” (difficult to

achieve)

• A. Alexanderian and A. K. Saibaba, Efficient D-optimal design of experiments for infinite-

dimensional Bayesian linear inverse problems, ArXiv, 2017

• A. Alexanderian, N. Petra, G. Stadler, and O. Ghattas, A fast and scalable method for A-optimal

design of experiments for infinite-dimensional Bayesian nonlinear inverse problems, SISC, 2016

Optimality Criteria

– Bayes risk seeks to minimize expected error between true & inferred parameter – A-optimal seeks to minimize the expectation of the trace of posterior covariance (i.e. the average variance of m) – Expected information Gain (EIG) seeks to maximize the expected Kullback-Leibler divergence from the prior to the posterior : Note: EIG reduces to D-optimal for linear inverse problems

Z

Y

1 2kmMAP(y, w) myk2π(y)dy 7! min

Z

Y

trace [Cpost(y, w)] π(y)dy 7! min

Z

Y

DKL(µpost(·|y; w)kµpr) π(y) dy 7! max

• A. Alexanderian, P. Gloor, and O. Ghattas, On Bayesian A- and D-optimal experimental designs in

infinite dimensions, Bayesian Analysis, 2016.

• EIG is the expected value of the KL divergence from prior to posterior
• ver all possible realizations of the data

where

• Sample average approximation (SAA) gives

where

Expected Information Gain (EIG)

EIG := EY [DKL(µpost(·|y)kµpr)] = Z

Y

DKL(µpost(·|y)kµpr) π(y) dy = Z

M

Z

Y

DKL(µpost(·|y)kµpr)πlike(y|m)dy dµpr(m)

DKL(µpost(·|y)kµpr) := Z

M

ln ✓dµpost(m|y) dµpr(m) ◆ dµpost(m|y)

EIG ⇡ 1 Ns

Ns

X

i=1

DKL(µpost(·|yi)kµpr)

yi = Bu(mi) + εi with mi ∼ µpr, εi ∼ µn i.i.d.

Evaluation of EIG

• Double Loop Monte Carlo (DLMC):

– Numerical underflow issues for concentrated posteriors – Extremely expensive to compute

• Monte Carlo with Laplace Approximation (MCLA):

– Laplace approx introduces bias, which decays with the number of

• bservations

– It exploits closed form expression for KL distance between Gaussian distributions – Efficient evaluation of the function & its gradient

• Double Loop Monte Carlo with Importance Sampling (DLMC-IS):

– Bias can be made arbitrarily small by increasing number of inner samples – Uses Laplace approx as biasing distribution to evaluate evidence – Extremely expensive to compute

• J. Beck, B. Mansour Dia, L. F. Espath, Q. Long, and R. Tempone, Fast Bayesian experimental

design: Laplace-based importance sampling for the expected information gain, CMAME, 2018

Laplace Approximation of the Posterior

• Computing the KL divergence is intractable for large-scale problems
• Laplace approximation is a fast and scalable approach for:

– Approximates posterior as Gaussian with mean given by maximum a posteriori (MAP) point, and covariance by inverse of Hessian of negative log posterior – Combined with low-rank approximation of Hessian of negative log likelihood, renders posterior estimate for large-scale inverse problems tractable – Permits tractable computation of KL divergence (since between two Gaussians)

• Our approach exploits three approximations:

– Sample average approximation (SAA) of expectation over the data – Laplace approximation to the true posterior – Spectrum truncation & randomized eigensolver to efficiently represent Hessian

• f negative log likelihood
• U. Villa, O. Ghattas. Scalable Laplace approximation-based optimal experimental

design for infinite dimensional Bayesian inverse problems, in preparation, 2019

Bayesian Inverse Problem

• Forward problem

where: – is state – is parameter to be inferred (random field)

• Observation process ( is observation operator)
• Bayes theorem

where

B r(u, m) = 0

u

m ∼ N(mpr, Cpr)

y = Bu(m) + ε, ε ∼ N(0, Γn)

dµpost(m|y) ∝ πlike(y|m) dµpr(m)

πlike(y|m) / exp ✓ 1 2kBu(m) yk2

Γ−1

n

◆

Laplace Approximation

• For nonlinear inverse problems, the posterior is non-Gaussian
• Laplace approximation of posterior reads

where mMAP is the maximum a posteriori point (minimizes negative log posterior) and Cpost is the inverse of the (Gauss Newton) Hessian at mMAP

• Closed form expression for KL distance between Gaussian distributions

m ∼ ν := N (mMAP, Cpost)

mMAP := argmin

m∈M

1 2kBu(m) yk2

Γ−1

n + 1

2km mprk2

C−1

pr

Cpost =

• Hm + C−1

pr

−1

DKL(νkµpr) = 1 2 h logdet ⇣ C

1 2

prC−1 postC

1 2

pr

⌘ trace ⇣ C

1 2

post(C−1 post C−1 pr )C

1 2

post

⌘ +kmMAP mprk2

C−1

pr

i

• For nonlinear inverse problems, the posterior is non-Gaussian
• Laplace approximation of posterior reads

where mMAP is the maximum a posteriori point (maximizes negative log posterior) and Cpost is the inverse of the (Gauss Newton) Hessian at mMAP

• Closed form expression for KL distance between Gaussian distributions

where is prior-preconditioned data misfit Hessian

Laplace Approximation

m ∼ ν := N (mMAP, Cpost)

mMAP := argmin

m∈M

1 2kBu(m) yk2

Γ−1

n + 1

2km mprk2

C−1

pr

Cpost =

• Hm + C−1

pr

−1

e Hm = C

1 2

prHmC

1 2

pr

DKL(νkµpr) = 1 2 h logdet ⇣ I + e Hm ⌘ trace ⇣ C

1 2

postHmC

1 2

post

⌘ + kmMAP mprk2

C−1

pr

i

Finding the MAP point scalably

mMAP := argmin

m∈M

1 2kBu(m) yk2

Γ−1

n + 1

2km mprk2

C−1

pr

• Equivalent to solving a deterministic inverse problem:
• Minimize a regularized data misfit
• Can be done at cost (measured in forw/adj PDE solves) that is independent of

the parameter & data dimensions; depends only on “information dimension” (how many modes in parameter space are informed by the data)

• But must be done properly:
• Hessian-free inexact Newton-conjugate gradient solver (terminate CG when

negative curvature direction encountered & avoid oversolving)

• Hessian action in CG directions formed by solving a pair of second-order

linearized forward/adjoint PDEs

• Precondition by prior covariance to yield compact perturbation of identity

à CG converges in dimension-independent steps

• Globalize via continuation and trust region or line search

Example: Scalability of MAP Point Computation for Ice Sheet Inverse Problem

• T. Isaac, N. Petra, G. Stadler, O. Ghattas, Scalable and efficient algorithms for the propagation of

uncertainty from data through inference to prediction for large-scale problems, with application to flow of the Antarctic ice sheet, Journal of Computational Physics, 2015.

Spectrum Truncation for DKL Computation

• Computing logdet and trace of large implicitly defined operators is

intractable using naïve algorithms

• Exploit fact that prior-preconditioned Hessian misfit is compact:
• Prior covariance operator is compact
• Hessian of data misfit is usually compact (observables are

sensitive to a limited number of modes in parameter space)

• Make low rank approximation of
• Let be the r dominant eigenvalues of , then
• Randomized eigensolvers can extract accurate approximation of

dominant eigenvalues of at a cost (measured in PDE solves)

• that is independent of the parameter and data dimensions
• and depends only on the numerical rank r
• Randomized eigensolver requires only Hessian actions; can be formed
• perator-free at cost of one incremental forward/adjoint PDE solve

e Hm

Cpr

Hm

e Hm

λj (j = 1, . . . , r)

e Hm

logdet ⇣ I + e Hm ⌘ ≈

r

X

j=1

ln(1 + λj), trace ⇣ C

1 2

postHmC

1 2

post

⌘ ≈

r

X

j=1

λj 1 + λj

e Hm

Randomized Eigensolver

• 1. Generate i.i.d. Gaussian matrix Ω ∈ Rn×r
• 2. Form Y = e

HΩ,

• 3. Compute Q = orthonormal basis for Y
• 4. Compute T = QT e

HQ

• 5. Decompose T = ZΛZT .
• 6. Low rank approximation e

H = VΛVT , where V ∈ Rn×r := QZ e

• N. Halko, P. G. Martinsson, and J. A. Tropp. Finding structure with randomness: Probabilistic

algorithms for constructing approximate matrix decompositions. SIAM Review, 2011.

• A. K. Saibaba, A. Alexanderian, and I. C. F. Ipsen, Randomized matrix-free trace and log-

determinant estimators, Numerische Mathematik, 2017.

Main cost: 2r incremental forward solves, 2r incremental adjoint solves, 4r prior square root applications Double pass algorithm:

Example: Inversion for basal friction field in Antarctica

InSAR-based ice surface velocity observations Inferred mean of basal friction field

Reconstructed ice surface velocity field (based on inferred mean of basal friction field) Inferred uncertainty in basal friction field (standard deviation of Laplace approx of posterior of log basal friction)

Joint work with Tobin Isaac (Georgia Tech), Noemi Petra (UC-Merced), Georg Stadler (NYU)

Nonlinear Stokes Ice Sheet Flow

Forward & Adjoint Stokes Equations

Incremental Forward Stokes Equation

Incremental Adjoint Stokes Equation

27

Spectrum of the prior-preconditioned data misfit Hessian

e Hm = C

1 2

prHmC

1 2

pr

• Spectrum of for Antarctica inverse problem with

410K and 1.19M basal sliding parameters (observed to decay like i-3)

• ~5000 dominant modes, indep. of parameter & data dimensions
• intrinsic problem dimension depends on information content of data

Design Space & Weighted Likelihood

design := ⇢ x1, . . . , xnd w1, . . . , wnd

• <latexit sha1_base64="3WlDCJsQEo1P5s4Zjr1c9RiegPg=">AIpXichVtj9w0EM6WQpcU6JXyrV8sThUFnU4bVKkICam9o+ohcWVRe3dFlyVyEu+eOcOtnPXleX+T7zQxg7SbvrpCLSaifzPDPzZMYvec2o0rPZP5MbH938+JNb0/j259/sWdnbtfnirRyIKcFIJ+TrHijDKyYmpHXtS4yhk5y8PHX52RaSigr/S65osKrzidEkLrMGV7bxNXmjTUkUXGLfvwJpYwsdWriNCcryg2WEq+tYTY2ab5Eb2yW7AGnFrtod5leFbaNI1TCok2Ca3Dw3FKeNmlQ6mkqwud2mxnd7Y/8w8aGkln7EbdM8/u3rqTlqJoKsJ1wbBS58ms1gvIq2nBCFRpFKlxcYlX5BxMjiuiFsY3yqIH4CnRUkj4cY28dzPC4EqpdZUDs8L6QoWYc45h541e/rAwlNeNJrxoCy0bhrRAruopJIUmq3BwIWkoBUVF1jiQsNs4jl5LoQVYWhQzAQpe15sgCLspKY3cTaOH6ATrGkGOaqtumwRrAm1jR2zP8LqDFXo9BT21t/fjdO8MHvEw1oLx3wEnoA5VPXlAIzc2LDZDW1lQjXp8fjwBeGvwPhTmwk9UmGFDmWIaSjgeScPmXNfWI1+f+ewTwkuB/KMmBnaQ2wZACm8Ca68D7M1Gh0LNQKErwoHAYEeXYac86HR5Y1CVixyq8sykskLO9qtoLqmQe2gulN6mV0tnZSmDYAXOyBA0CbFvQakQ1FvfM/hAG5rvIM/VOuwLxUyR0rWsiKYdyvig+JxvcU4fjp3/djivBCywmyj6osB5aiak3lR3ZLFza9XvWNYy+3cUbAT4CRvErvYRz+0KIMmw0F1SRjU7zwTdTBsHWGS3Dvw5IfZ6eNcjzvGc8x+DJeABfB/hZgB9sTfsgzO4X4Dpw/pYPlv0fI4FuoFdYklpRJkJdtd8Ufi8XdnH8ErAB3cbQ6ImEmshg8NS1fZH649w90QpvX7q+8hcidud0N9GxRmYgUDHUnQAf+fActVRXmboeElHP9EO3qQr6OFjcFXfd+0Ns/cd8dwhSbhTk0Tr/fT2b7ye+Pdp8cdJfpNLofR09jJLocfQkOorm0UlURP9O4sm9yVfTb6bH01fT05Z6Y9LF3Iu2nmn2H4GRMFY=</latexit><latexit sha1_base64="3WlDCJsQEo1P5s4Zjr1c9RiegPg=">AIpXichVtj9w0EM6WQpcU6JXyrV8sThUFnU4bVKkICam9o+ohcWVRe3dFlyVyEu+eOcOtnPXleX+T7zQxg7SbvrpCLSaifzPDPzZMYvec2o0rPZP5MbH938+JNb0/j259/sWdnbtfnirRyIKcFIJ+TrHijDKyYmpHXtS4yhk5y8PHX52RaSigr/S65osKrzidEkLrMGV7bxNXmjTUkUXGLfvwJpYwsdWriNCcryg2WEq+tYTY2ab5Eb2yW7AGnFrtod5leFbaNI1TCok2Ca3Dw3FKeNmlQ6mkqwud2mxnd7Y/8w8aGkln7EbdM8/u3rqTlqJoKsJ1wbBS58ms1gvIq2nBCFRpFKlxcYlX5BxMjiuiFsY3yqIH4CnRUkj4cY28dzPC4EqpdZUDs8L6QoWYc45h541e/rAwlNeNJrxoCy0bhrRAruopJIUmq3BwIWkoBUVF1jiQsNs4jl5LoQVYWhQzAQpe15sgCLspKY3cTaOH6ATrGkGOaqtumwRrAm1jR2zP8LqDFXo9BT21t/fjdO8MHvEw1oLx3wEnoA5VPXlAIzc2LDZDW1lQjXp8fjwBeGvwPhTmwk9UmGFDmWIaSjgeScPmXNfWI1+f+ewTwkuB/KMmBnaQ2wZACm8Ca68D7M1Gh0LNQKErwoHAYEeXYac86HR5Y1CVixyq8sykskLO9qtoLqmQe2gulN6mV0tnZSmDYAXOyBA0CbFvQakQ1FvfM/hAG5rvIM/VOuwLxUyR0rWsiKYdyvig+JxvcU4fjp3/djivBCywmyj6osB5aiak3lR3ZLFza9XvWNYy+3cUbAT4CRvErvYRz+0KIMmw0F1SRjU7zwTdTBsHWGS3Dvw5IfZ6eNcjzvGc8x+DJeABfB/hZgB9sTfsgzO4X4Dpw/pYPlv0fI4FuoFdYklpRJkJdtd8Ufi8XdnH8ErAB3cbQ6ImEmshg8NS1fZH649w90QpvX7q+8hcidud0N9GxRmYgUDHUnQAf+fActVRXmboeElHP9EO3qQr6OFjcFXfd+0Ns/cd8dwhSbhTk0Tr/fT2b7ye+Pdp8cdJfpNLofR09jJLocfQkOorm0UlURP9O4sm9yVfTb6bH01fT05Z6Y9LF3Iu2nmn2H4GRMFY=</latexit><latexit sha1_base64="3WlDCJsQEo1P5s4Zjr1c9RiegPg=">AIpXichVtj9w0EM6WQpcU6JXyrV8sThUFnU4bVKkICam9o+ohcWVRe3dFlyVyEu+eOcOtnPXleX+T7zQxg7SbvrpCLSaifzPDPzZMYvec2o0rPZP5MbH938+JNb0/j259/sWdnbtfnirRyIKcFIJ+TrHijDKyYmpHXtS4yhk5y8PHX52RaSigr/S65osKrzidEkLrMGV7bxNXmjTUkUXGLfvwJpYwsdWriNCcryg2WEq+tYTY2ab5Eb2yW7AGnFrtod5leFbaNI1TCok2Ca3Dw3FKeNmlQ6mkqwud2mxnd7Y/8w8aGkln7EbdM8/u3rqTlqJoKsJ1wbBS58ms1gvIq2nBCFRpFKlxcYlX5BxMjiuiFsY3yqIH4CnRUkj4cY28dzPC4EqpdZUDs8L6QoWYc45h541e/rAwlNeNJrxoCy0bhrRAruopJIUmq3BwIWkoBUVF1jiQsNs4jl5LoQVYWhQzAQpe15sgCLspKY3cTaOH6ATrGkGOaqtumwRrAm1jR2zP8LqDFXo9BT21t/fjdO8MHvEw1oLx3wEnoA5VPXlAIzc2LDZDW1lQjXp8fjwBeGvwPhTmwk9UmGFDmWIaSjgeScPmXNfWI1+f+ewTwkuB/KMmBnaQ2wZACm8Ca68D7M1Gh0LNQKErwoHAYEeXYac86HR5Y1CVixyq8sykskLO9qtoLqmQe2gulN6mV0tnZSmDYAXOyBA0CbFvQakQ1FvfM/hAG5rvIM/VOuwLxUyR0rWsiKYdyvig+JxvcU4fjp3/djivBCywmyj6osB5aiak3lR3ZLFza9XvWNYy+3cUbAT4CRvErvYRz+0KIMmw0F1SRjU7zwTdTBsHWGS3Dvw5IfZ6eNcjzvGc8x+DJeABfB/hZgB9sTfsgzO4X4Dpw/pYPlv0fI4FuoFdYklpRJkJdtd8Ufi8XdnH8ErAB3cbQ6ImEmshg8NS1fZH649w90QpvX7q+8hcidud0N9GxRmYgUDHUnQAf+fActVRXmboeElHP9EO3qQr6OFjcFXfd+0Ns/cd8dwhSbhTk0Tr/fT2b7ye+Pdp8cdJfpNLofR09jJLocfQkOorm0UlURP9O4sm9yVfTb6bH01fT05Z6Y9LF3Iu2nmn2H4GRMFY=</latexit><latexit sha1_base64="3WlDCJsQEo1P5s4Zjr1c9RiegPg=">AIpXichVtj9w0EM6WQpcU6JXyrV8sThUFnU4bVKkICam9o+ohcWVRe3dFlyVyEu+eOcOtnPXleX+T7zQxg7SbvrpCLSaifzPDPzZMYvec2o0rPZP5MbH938+JNb0/j259/sWdnbtfnirRyIKcFIJ+TrHijDKyYmpHXtS4yhk5y8PHX52RaSigr/S65osKrzidEkLrMGV7bxNXmjTUkUXGLfvwJpYwsdWriNCcryg2WEq+tYTY2ab5Eb2yW7AGnFrtod5leFbaNI1TCok2Ca3Dw3FKeNmlQ6mkqwud2mxnd7Y/8w8aGkln7EbdM8/u3rqTlqJoKsJ1wbBS58ms1gvIq2nBCFRpFKlxcYlX5BxMjiuiFsY3yqIH4CnRUkj4cY28dzPC4EqpdZUDs8L6QoWYc45h541e/rAwlNeNJrxoCy0bhrRAruopJIUmq3BwIWkoBUVF1jiQsNs4jl5LoQVYWhQzAQpe15sgCLspKY3cTaOH6ATrGkGOaqtumwRrAm1jR2zP8LqDFXo9BT21t/fjdO8MHvEw1oLx3wEnoA5VPXlAIzc2LDZDW1lQjXp8fjwBeGvwPhTmwk9UmGFDmWIaSjgeScPmXNfWI1+f+ewTwkuB/KMmBnaQ2wZACm8Ca68D7M1Gh0LNQKErwoHAYEeXYac86HR5Y1CVixyq8sykskLO9qtoLqmQe2gulN6mV0tnZSmDYAXOyBA0CbFvQakQ1FvfM/hAG5rvIM/VOuwLxUyR0rWsiKYdyvig+JxvcU4fjp3/djivBCywmyj6osB5aiak3lR3ZLFza9XvWNYy+3cUbAT4CRvErvYRz+0KIMmw0F1SRjU7zwTdTBsHWGS3Dvw5IfZ6eNcjzvGc8x+DJeABfB/hZgB9sTfsgzO4X4Dpw/pYPlv0fI4FuoFdYklpRJkJdtd8Ufi8XdnH8ErAB3cbQ6ImEmshg8NS1fZH649w90QpvX7q+8hcidud0N9GxRmYgUDHUnQAf+fActVRXmboeElHP9EO3qQr6OFjcFXfd+0Ns/cd8dwhSbhTk0Tr/fT2b7ye+Pdp8cdJfpNLofR09jJLocfQkOorm0UlURP9O4sm9yVfTb6bH01fT05Z6Y9LF3Iu2nmn2H4GRMFY=</latexit>
• xi: candidate sensor locations
• wi: weights
• Ideally:

Combinatorial complexity

• Relax:
• Sparsity-inducing penalization:

0 ≤ wi ≤ 1

wi ∈ {0, 1}

P✏(w) ⇠ ⇢ kwk1 ✏ = 1 kwk0 ✏ ! 0

Weighted likelihood function πlike(y|m; w) / exp ✓ 1 2kBu(m) yk2

Γ−1

W

◆ , ΓW = diag ⇣ w− 1

2

⌘ Γn diag ⇣ w− 1

2

⌘

• A. Alexanderian, N. Petra, G. Stadler, and O. Ghattas, A-optimal design of experiments for infinite

dimensional Bayesian linear inverse problems with regularized l0-sparsification, SISC, 2014

OED Problem with EIG-MCLA

yi = Bu(mi) + εi with mi ∼ µpr, εi ∼ µn i.i.d.

min

w∈W

1 2 Ns

Ns

X

i=1

2 4

ri

X

j=1

✓ ln(1 + λij) λij 1 + λij ◆ + kmi mprk2

C−1

pr

3 5 + βP(w)

     mi := argmin

m∈M 1 2kB ui yik2 Γ−1

W + 1

2km mprk2 C−1

pr

subject to: hr(ui, m), ˜ pi = 0 8˜ p 2 U i = 1, . . . , Ns

( h˜ aij, Hm(mi, w) aiji = λijh˜ aij, C−1

pr aiji

˜ aij 2 U kaijk2

C−1

pr = 1

EIG-MCLA Sparsity-inducing penalization Subject to:

• Inverse problem

with

• Eigenproblem

j = 1, . . . , ri; i = 1, . . . , Ns

OED Problem with EIG-MCLA

min

w∈W

1 2 Ns

Ns

X

i=1

2 4

ri

X

j=1

✓ ln(1 + λij) λij 1 + λij ◆ + kmi mprk2

C−1

pr

3 5 + βP(w)

Subject to:

• Inverse problem
• Eigenproblem

   hr(ui, mi), ˜ pii = 0 8˜ pi 2 U h∂ur|?[˜ ui], pii = hB˜ ui, Γ−1

W (B ui yi)i

8˜ ui 2 U h ˜ mi, C−1

pr (mi mpr)i = h∂mr|?[ ˜

mi], pii 8 ˜ mi 2 M i = 1, . . . , Ns

             h∂ur|?[vij], ˜ qiji = h∂mr|?[aij], ˜ qiji 8˜ qij 2 U h∂ur|?[˜ vij], qiji = hB˜ vij, Γ−1

W Biviji

8˜ vij 2 U h∂mr|?[ ˜ mij], qiji = λijh ˜ mij, C−1

pr aiji

8 ˜ mij 2 M haij, C−1

pr aiji˜

λij = ˜ λij 8˜ λij 2 R i = 1, . . . , Ns, j = 1, . . . , ri.

Cost Per Optimization Iteration (Measured in Forward/Adjoint PDE Solves)

• Objective functional:
• MAP Point: Solve an inverse problem using inexact Newton-CG

algorithm for each SAA sample.

• Eigenvalues of Hessian at MAP point: Compute the low-rank

factorization of the prior-preconditioned Hessian of the data misfit using randomized eigensolver for each SAA sample

• Gradient (needed for quasi-Newton optimization):
• Adjoint of Hessian eigenproblem: Adjoint eigenvalues are ~ free
• Adjoint of MAP point necessary conditions: Solve a linear system

involving the Hessian of the inverse problem for each SAA sample; spectral decomposition of Hessian is available, so cost is ~ free Newton method, variational adjoint-based methods, & randomized eigensolvers allow for efficient computation of objective functional and its gradient at a cost (measured in number of PDE solves) that is asymptotically independent of parameter dimension and the data/sensor dimension.

• Prior distribution:
• Forward model
• Design space

9 equispaced candidate sensor locations Find 4 optimal locations

• Small problem: Only 126 different designs to compare using brute force

approach

1D Elliptic PDE Coefficient Inversion

49

• U. Villa, N. Petra, and O. Ghattas, hIPPYlib: An Extensible Software Framework for

Large-Scale Deterministic and Linearized Bayesian Inverse Problems http://hippylib.github.io m ⇠ N(0, Cpr),

Cpr = σ2 κ−2I − ∆ −1

( − d

dx

• em du

dx

• = 0,

x ∈ (0, 1) u(0) = 0, u(1) = 1

Blue squares: All 126 possible designs (9 choose 4) EIG: Double loop Monte Carlo (w/IS) with 1000 inner + 1000 outer samples EIG-LA: Monte Carlo-Laplace Approximation with 1000 outer samples Error bars: 3 standard deviations (99.7% CI) Correlation coefficient: 0.98

Accuracy of Laplace Approximation to EIG

Sparsifying Penalty

Penalty

• A. Alexanderian, N. Petra, G. Stadler, O. Ghattas,

A-optimal design of experiments for infinite- dimensional Bayesian linear inverse problems with regularized l0-sparsification, SISC, 2014

Pε(wi) =   

wi ε ,

0 ≤ wi ≤ 1

2ε

pε(wi),

1 2ε ≤ wi ≤ 2ε

1, 2ε ≤ wi ≤ 1 pε : cubic polynomial s.t. Pε ∈ C1

Sparsifying Penalty

Optimal designs (Ns = 512) Penalty

Effect of SAA Sample Size

Based on all 126 possible sensor locations L-BFGS optimizer with bound constraints

2D Elliptic PDE Coefficient Inversion

• Prior distribution:
• Forward model:
• Design space:

10-by-10 grid of candidate sensor locations Select 15 locations

• Brute force approach is not feasible:
• Solve continuous relaxation using our
• ptimization method

   div (emru) = 0 in Ω u = ud

• n ΓD ⇢ ∂Ω

emru · n = 0

• n ΓN ⇢ ∂Ω

m ⇠ N(0, Cpr),

), Cpr = σ2 [I div(Θr)]−2

✓ 100 15 ◆ ≈ 1017

combinations m: random realization u(m): forward solution

Blue squares: 123 random designs tested (from 1017 combinations) + 5 optimal in red EIG: Double loop Monte Carlo (w/IS) with 1000 inner + 1000 outer samples EIG-LA: Monte Carlo-Laplace Approximation with 1000 outer samples Error bars: 3 standard deviations (99.7% CI) Correlation coefficient: 0.99

Accuracy of Laplace Approximation

Sparsifying Penalty (continuation toward l0)

Ns = 512. L-BFGS optimizer with bound constraints

Effect of SAA Sample Size

100-choose-15 sensor locations L-BFGS optimizer with bound constraints

Computational Cost: Design-Dim

Outer BFGS Iterations PDE solves Only a moderate increase in computational cost w.r.t. the dimension of the design space Number of candidate sensor locations

Computational Cost: Param-Dim

Outer BFGS Iterations PDE solves Computational cost (measured in terms of PDE solves) is independent of parameter dimension Number of degrees of freedom in parameter space

• Repeated evaluation of EIG using double Loop Monte Carlo (w/IS) is

prohibitive for large-scale problems

• MC-Laplace approximation allows for efficient & accurate estimation of EIG:

– SAA over possible realization of the data – Closed form expression for KL distance between Gaussian involves “only” log-determinant & trace of (preconditioned) posterior covariance – Randomized eigensolvers permit efficient evaluation of these invariants

• Gradient computation using variational adjoint methods
• Cost measured in forward PDE solves is independent of parameter space

dimension, and slowly increases with experimental design space dimension

• Ongoing/future work:

– Employ Laplace approximation as control variate for expectation over data – For Helmholtz, rank of Hessian grows with increasing frequency; pursue compression schemes that do not depend on global low rank (H-matrix, product-convolution) – Pursue Newton methods for OED problem to achieve sensor-dimension independence

Conclusions & Ongoing Work