A Computational Framework for Infinite-Dimensional Bayesian Inverse Problems, Part II: Stochastic Newton MCMC with Application to Ice Sheet Flow Inverse Problems

Abstract

Abstract. We address the numerical solution of infinite-dimensional inverse problems in the framework of Bayesian inference. In the Part I [11] companion to this paper, 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 probability density functions (pdfs) arising upon discretization of 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 covariance operator is given by the inverse of the local Hessian of the negative log posterior pdf. The construction of the covariance is made tractable by invoking a low-rank approximation of the 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.

The Bayesian Approach to Inverse Problems

“…For example, on the analysis side, the idea of MAP estimators, which links the Bayesian approach with classical regularization, developed for Gaussian priors in [30], has recently been extended to other prior models in [47]; the study of contraction of the posterior distribution to a Dirac measure on the truth underlying the data is undertaken in [3,4,99]. On the algorithmic side, algorithms for Bayesian inversion in geophysical applications are formulated in [16,81], and on the computational statistics side, methods for optimal experimental design are formulated in [5,6]. All of these cited papers build on the framework developed in detail here and first outlined in [92].…”

<i>Albany/FELIX</i>: a parallel, scalable and robust, finite element, first-order Stokes approximation ice sheet solver built for advanced analysis

“…Recently, these approaches have been extended to simultaneously optimize both model parameter fields and uncertain initial condition fields, while also accounting for forcing from climate models in order to minimize transient shocks when coupling to climate forcing (Perego et al, 2014). Other recent and noteworthy optimization improvements include the assimilation of time-dependent observations (e.g., Goldberg and Heimbach, 2013) and the estimation of formal uncertainties for optimized parameter fields (Petra et al, 2015).…”

Inversion of geothermal heat flux in a thermomechanically coupled nonlinear Stokes ice sheet model

“…An approximation of the inverse Hessian can be computed even for large-scale inverse problems by exploiting low-rank properties that are typical for many illposed inverse problems (Flath et al, 2011;Petra et al, 2014;Kalmikov and Heimbach, 2014).…”