Noemi Petra

A computational framework for infinite-dimensional Bayesian inverse problems, part II: Stochastic Newton mcmc with application to ice sheet flow inverse problems

We address the numerical solution of infinite-dimensional inverse problems in the framework of Bayesian inference. In the Part I companion to this paper (arXiv.org:1308.1313), we considered the linearized infinite-dimensional inverse problem. Here in Part II, we relax the linearization assumption and consider the fully nonlinear infinite-dimensional inverse problem using a Markov chain Monte Carlo (MCMC) sampling method. To address the challenges of sampling high-dimensional pdfs arising from Bayesian inverse problems governed by PDEs, we build on the stochastic Newton MCMC method. This method exploits problem structure by taking as a proposal density a local Gaussian approximation of the posterior pdf, whose construction is made tractable by invoking a low-rank approximation of its data misfit component of the Hessian. Here we introduce an approximation of the stochastic Newton proposal in which we compute the low-rank-based Hessian at just the MAP point, and then reuse this Hessian at each MCMC step. We compare the performance of the proposed method to the original stochastic Newton MCMC method and to an independence sampler. The comparison of the three methods is conducted on a synthetic ice sheet inverse problem. For this problem, the stochastic Newton MCMC method with a MAP-based Hessian converges at least as rapidly as the original stochastic Newton MCMC method, but is far cheaper since it avoids recomputing the Hessian at each step. On the other hand, it is more expensive per sample than the independence sampler; however, its convergence is significantly more rapid, and thus overall it is much cheaper. Finally, we present extensive analysis and interpretation of the posterior distribution, and classify directions in parameter space based on the extent to which they are informed by the prior or the observations.

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

The majority of research on efficient and scalable algorithms in computational science and engineering has focused on the forward problem: given parameter inputs, solve the governing equations to determine output quantities of interest. In contrast, here we consider the broader question: given a (large-scale) model containing uncertain parameters, (possibly) noisy observational data, and a prediction quantity of interest, how do we construct efficient and scalable algorithms to (1) infer the model parameters from the data (the deterministic inverse problem), (2) quantify the uncertainty in the inferred parameters (the Bayesian inference problem), and (3) propagate the resulting uncertain parameters through the model to issue predictions with quantified uncertainties (the forward uncertainty propagation problem)?We present efficient and scalable algorithms for this end-to-end, data-to-prediction process under the Gaussian approximation and in the context of modeling the flow of the Antarctic ice sheet and its effect on loss of grounded ice to the ocean. The ice is modeled as a viscous, incompressible, creeping, shear-thinning fluid. The observational data come from satellite measurements of surface ice flow velocity, and the uncertain parameter field to be inferred is the basal sliding parameter, represented by a heterogeneous coefficient in a Robin boundary condition at the base of the ice sheet. The prediction quantity of interest is the present-day ice mass flux from the Antarctic continent to the ocean.We show that the work required for executing this data-to-prediction process-measured in number of forward (and adjoint) ice sheet model solves-is independent of the state dimension, parameter dimension, data dimension, and the number of processor cores. The key to achieving this dimension independence is to exploit the fact that, despite their large size, the observational data typically provide only sparse information on model parameters. This property can be exploited to construct a low rank approximation of the linearized parameter-to-observable map via randomized SVD methods and adjoint-based actions of Hessians of the data misfit functional.

MODELING AND DESIGN OPTIMIZATION OF A RESONANT OPTOTHERMOACOUSTIC TRACE GAS SENSOR

Trace gas sensors that are compact and portable are being deployed for use in a variety of applications including disease diagnosis via breath analysis, monitoring of atmospheric pollutants and greenhouse gas emissions, control of industrial processes, and for early warning of terrorist threats. One such sensor is based on optothermal detection and uses a modulated laser source and a quartz tuning fork resonator to detect trace gases. In this paper we introduce the first mathematical model of such a resonant optothermoacoustic sensor. The model is solved via the finite element method and couples heat transfer and thermoelastic deformation to determine the strength of the generated signal. Numerical simulations validate the experimental observation that the source location that produces the maximum signal is near the junction of the tines of the tuning fork. Determining an optimally designed sensor requires maximizing the signal as a function of the geometry of the quartz tuning fork (length and width of the tines, etc). To avoid difficulties from numerical differentiation we chose to solve the optimization problem using the derivative-free mesh adaptive direct search algorithm. An optimal tuning fork constrained to resonate at a frequency close to the 32.8 kHz resonance frequency of many commercially available tuning forks produces a signal that is three times larger than the one obtained with the current experimental design. Moreover, the optimal tuning fork found without imposing any constraint on the resonance frequency produces a signal that is 24 times greater than that obtained with the current sensor.

Numerical and experimental investigation for a resonant optothermoacoustic sensor

A theoretical study of a resonant optothermoacoustic sensor employing a laser source and a quartz tuning fork receiver validates experimental results showing that the source should be positioned near the base of the receiver.

Mean-Variance Risk-Averse Optimal Control of Systems Governed by PDEs with Random Parameter Fields Using Quadratic Approximations

We present a method for optimal control of systems governed by partial differential equations (PDEs) with uncertain parameter fields. We consider an objective function that involves the mean and variance of the control objective, leading to a risk-averse optimal control problem. To make the problem tractable, we invoke a quadratic Taylor series approximation of the control objective with respect to the uncertain parameter. This enables deriving explicit expressions for the mean and variance of the control objective in terms of its gradients and Hessians with respect to the uncertain parameter. The risk-averse optimal control problem is then formulated as a PDE-constrained optimization problem with constraints given by the forward and adjoint PDEs defining these gradients and Hessians. The expressions for the mean and variance of the control objective under the quadratic approximation involve the trace of the (preconditioned) Hessian and are thus prohibitive to evaluate. To address this, we employ trace estimators that only require a modest number of Hessian-vector products. We illustrate our approach with two problems: the control of a semilinear elliptic PDE with an uncertain boundary source term, and the control of a linear elliptic PDE with an uncertain coefficient field. For the latter problem, we derive adjoint-based expressions for efficient computation of the gradient of the risk-averse objective with respect to the controls. Our method ensures that the cost of computing the risk-averse objective and its gradient with respect to the control, measured in the number of PDE solves, is independent of the (discretized) parameter and control dimensions, and depends only on the number of random vectors employed in the trace estimation. Finally, we present a comprehensive numerical study of an optimal control problem for fluid flow in a porous medium with uncertain permeability field.

A Bayesian Approach for Parameter Estimation With Uncertainty for Dynamic Power Systems

We address the problem of estimating the uncertainty in the solution of power grid inverse problems within the framework of Bayesian inference. We investigate two approaches, an adjoint-based method and a stochastic spectral method. These methods are used to estimate the maximum a posteriori point of the parameters and their variance, which quantifies their uncertainty. Within this framework, we estimate several parameters of the dynamic power system, such as generator inertias, which are not quantifiable in steady-state models. We illustrate the performance of these approaches on a 9-bus power grid example and analyze the dependence on measurement frequency, estimation horizon, perturbation size, and measurement noise. We assess the computational efficiency, and discuss the expected performance when these methods are applied to large systems.

hIPPYlib: An Extensible Software Framework for Large-Scale Inverse Problems

We present an extensible software framework, hIPPYlib, for solution of large-scale deterministic and Bayesian inverse problems governed by partial differential equations (PDEs) with infinite-dimensional parameter fields (which are high-dimensional after discretization). hIPPYlib overcomes the prohibitive nature of Bayesian inversion for this class of problems by implementing state-of-the-art scalable algorithms for PDE-based inverse problems that exploit the structure of the underlying operators, notably the Hessian of the log-posterior. The key property of the algorithms implemented in hIPPYlib is that the solution of the deterministic and linearized Bayesian inverse problem is computed at a cost, measured in linearized forward PDE solves, that is independent of the parameter dimension. The mean of the posterior is approximated by the MAP point, which is found by minimizing the negative log-posterior. This deterministic nonlinear least-squares optimization problem is solved with an inexact matrix-free Newton-CG method. The posterior covariance is approximated by the inverse of the Hessian of the negative log posterior evaluated at the MAP point. This Gaussian approximation is exact when the parameter-to-observable map is linear; otherwise, its logarithm agrees to two derivatives with the log-posterior at the MAP point, and thus it can serve as a proposal for Hessian-based MCMC methods. The construction of the posterior covariance is made tractable by invoking a low-rank approximation of the Hessian of the log-likelihood. Scalable tools for sample generation are also implemented. hIPPYlib makes all of these advanced algorithms easily accessible to domain scientists and provides an environment that expedites the development of new algorithms. hIPPYlib is also a teaching tool to educate researchers and practitioners who are new to inverse problems and the Bayesian inference framework.

Scalable Algorithms for Bayesian Inference of Large-Scale Models from Large-Scale Data

One of the greatest challenges in computational science and engineering today is how to combine complex data with complex models to create better predictions. This challenge cuts across every application area within CS&E, from geosciences, materials, chemical systems, biological systems, and astrophysics to engineered systems in aerospace, transportation, structures, electronics, biomedicine, and beyond. Many of these systems are characterized by complex nonlinear behavior coupling multiple physical processes over a wide range of length and time scales. Mathematical and computational models of these systems often contain numerous uncertain parameters, making high-reliability predictive modeling a challenge. Rapidly expanding volumes of observational data—along with tremendous increases in HPC capability—present opportunities to reduce these uncertainties via solution of large-scale inverse problems.

Towards Adjoint-Based Inversion of the Lamé Parameter Field for Slender Structures With Cantilever Loading

Continuum models of slender structures are effective in simulating the mechanics of nano-scale filaments. However, the accuracy of these simulations strictly depends on the knowledge of the constitutive laws that may in general be non-homogeneous. It necessitates an inverse problem framework that can leverage the data provided by physical experiments and molecular dynamics simulations to estimate the unknown parameters in the constitutive law. In this paper, we formulate a simple but representative inverse problem as a nonlinear least-squares optimization problem whose cost functional is the misfit between synthetic observations of a cantilever displacement field and model predictions. A Tikhonov regularization term is added to the cost functional to render the problem well-posed and account for observational error. We solve this optimization problem with an adjoint-based inexact Newton-conjugate gradient method. We show that the reconstruction of the Lame parameter field converges to the exact coefficient as the observation error decreases.Copyright