Bayesian Solution Of Inverse Problems Research Articles

HIPPYlib-MUQ: A Bayesian Inference Software Framework for Integration of Data with Complex Predictive Models under Uncertainty

Bayesian inference provides a systematic framework for integration of data with mathematical models to quantify the uncertainty in the solution of the inverse problem. However, the solution of Bayesian inverse problems governed by complex forward models described bypartial differential equations (PDEs)remains prohibitive with black-boxMarkov chain Monte Carlo (MCMC)methods. We present hIPPYlib-MUQ, an extensible and scalable software framework that contains implementations of state-of-the art algorithms aimed to overcome the challenges of high-dimensional, PDE-constrained Bayesian inverse problems. These algorithms accelerate MCMC sampling by exploiting the geometry and intrinsic low-dimensionality of parameter space via derivative information and low rank approximation. The software integrates two complementary open-source software packages, hIPPYlib and MUQ. hIPPYlib solves PDE-constrained inverse problems using automatically-generated adjoint-based derivatives, but it lacks full Bayesian capabilities. MUQ provides a spectrum of powerful Bayesian inversion models and algorithms, but expects forward models to come equipped with gradients and Hessians to permit large-scale solution. By combining these two complementary libraries, we created a robust, scalable, and efficient software framework that realizes the benefits of each and allows us to tackle complex large-scale Bayesian inverse problems across a broad spectrum of scientific and engineering disciplines. To illustrate the capabilities of hIPPYlib-MUQ, we present a comparison of a number of MCMC methods available in the integrated software on several high-dimensional Bayesian inverse problems. These include problems characterized by both linear and nonlinear PDEs, various noise models, and different parameter dimensions. The results demonstrate that large (∼ 50×) speedups over conventional black box and gradient-based MCMC algorithms can be obtained by exploiting Hessian information (from the log-posterior), underscoring the power of the integrated hIPPYlib-MUQ framework.

A probabilistic finite element method based on random meshes: A posteriori error estimators and Bayesian inverse problems

We present a novel probabilistic finite element method (FEM) for the solution and uncertainty quantification of elliptic partial differential equations based on random meshes, which we call random mesh FEM (RM-FEM). Our methodology allows to introduce a probability measure on classical FEMs to quantify the uncertainty due to numerical errors either in the context of a-posteriori error quantification or for FE based Bayesian inverse problems. The new approach involves only a perturbation of the mesh and an interpolation that are very simple to implement We present a posteriori error estimators and a rigorous a posteriori error analysis based uniquely on probabilistic information for standard piecewise linear FEM. A series of numerical experiments illustrates the potential of the RM-FEM for error estimation and validates our analysis. We furthermore demonstrate how employing the RM-FEM enhances the quality of the solution of Bayesian inverse problems, thus allowing a better quantification of numerical errors in pipelines of computations.

Bayesian inference of random fields represented with the Karhunen–Loève expansion

The integration of data into engineering models involving uncertain and spatially varying parameters is oftentimes key to obtaining accurate predictions. Bayesian inference is effective in achieving such an integration. Uncertainties related to spatially varying parameters are typically represented through random fields discretized into a finite number of random variables. The prior correlation length and variance of the field, as well as the number of terms in the random field discretization, have a considerable impact on the outcome of the Bayesian inference, which has received little attention in the literature. Here, we investigate the implications of different choices in the prior random field model on the solution of Bayesian inverse problems. We employ the Karhunen–Loève expansion for the representation of random fields. We show that a higher-order Karhunen–Loève discretization is required in Bayesian inverse problems as compared to standard prior random field representations, since the updated fields are non-homogeneous. Furthermore, the smoothing effect of the forward operator has a large influence on the posterior solution, particularly if we are analyzing quantities of interest that are sensitive to local random fluctuations of the inverse quantity. This is also reflected in posterior predictions, such as the estimation of rare event probabilities. We illustrate these effects analytically through a 1D cantilever beam with spatially varying flexibility, and numerically using a 2D linear elasticity example where the Young’s modulus is spatially variable.

Sequential Design of Computer Experiments for the Solution of Bayesian Inverse Problems

We present a sequential design strategy for efficient sampling of model functions during the solution of Bayesian inverse problems. The model function is assumed to be computationally expensive and therefore is described by a random field (such as a Gaussian process emulator). The sequential design strategy is a greedy one-step look ahead method, minimizing the Bayes risk with respect to a loss function measuring the quadratic $L^2$-error in the likelihood estimate. Four numerical examples demonstrate that the proposed sampling method is more efficient than space-filling, prior-based designs.

Stochastic Collocation Algorithms Using $l_1$-Minimization for Bayesian Solution of Inverse Problems

The Bayesian methodology has proven to be a convenient framework to solve inverse problems from available data with limited information. The main computational challenges in the Bayesian solution of inverse problems arise from the need for repeated evaluations of the forward model, as required by Markov chain Monte Carlo (MCMC) methods for posterior sampling. In this paper, we present an efficient technique for constructing stochastic surrogate models to accelerate the Bayesian inference approach for statistical inverse problems. The stochastic collocation algorithms using $l_1$-minimization (L1-SCM), based on generalized polynomial chaos, are used to construct a polynomial approximation of the forward solution over the support of the prior distribution. This approximation then defines a surrogate posterior probability density that is inexpensive to evaluate. A rigorous error analysis shows that the convergence rate of the posterior density is at least twice as fast as the convergence rate of the gPC expa...

Iterative regularization for ensemble data assimilation in reservoir models

We propose the application of iterative regularization for the development of ensemble methods for solving Bayesian inverse problems. In concrete, we construct (i) a variational iterative regularizing ensemble Levenberg-Marquardt method (IR-enLM) and (ii) a derivative-free iterative ensemble Kalman smoother (IR-ES). The aim of these methods is to provide a robust ensemble approximation of the Bayesian posterior. The proposed methods are based on fundamental ideas from iterative regularization methods that have been widely used for the solution of deterministic inverse problems (Katltenbacher et al. de Gruyter, Berlin 2008). In this work, we are interested in the application of the proposed ensemble methods for the solution of Bayesian inverse problems that arise in reservoir modeling applications. The proposed ensemble methods use key aspects of the regularizing Levenberg-Marquardt scheme developed by Hanke (Inverse Problems 13,79–95 1997) and that we recently applied for history matching in Iglesias (Comput. Geosci. 1–21 2013). Unlike most existing methods where the stopping criteria and regularization parameters are typically selected heuristically, in the proposed ensemble methods, the discrepancy principle is applied for (i) the selection of the regularization parameters and (ii) the early termination of the scheme. The discrepancy principle is key for the theory of iterative regularization, and the purpose of the present work is to apply this principle for the development of ensemble methods defined as iterative updates of solutions to linear ill-posed inverse problems. The regularizing and convergence properties of iterative regularization methods for deterministic inverse problems have long been established. However, the approximation properties of the proposed ensemble methods in the context of Bayesian inverse problems is an open problem. In the case where the forward operator is linear and the prior is Gaussian, we show that the tunable parameters of the proposed IR-enLM and IR-ES can be chosen so that the resulting schemes coincide with the standard randomized maximum likelihood (RML) and the ensemble smoother (ES), respectively. Therefore, the proposed methods sample from the posterior in the linear-Gaussian case. Similar to RML and ES methods, in the nonlinear case, one may not conclude that the proposed methods produce samples from the posterior. The present work provides a numerical investigation of the performance of the proposed ensemble methods at capturing the posterior. In particular, we aim at understanding the role of the tunable parameters that arise from the application of iterative regularization techniques. The numerical framework for our investigations consists of using a state-of-the art Markov chain Monte Carlo (MCMC) method for resolving the Bayesian posterior from synthetic experiments. The resolved posterior via MCMC then provides a gold standard against to which compare the proposed IR-enLM and IR-ES. Our numerical experiments show clear indication that the regularizing properties of the regularization methods applied for the computation of each ensemble have significant impact of the approximation properties of the proposed ensemble methods at capturing the Bayesian posterior. Furthermore, we provide a comparison of the proposed regularizing methods with respect to some unregularized methods that have been typically used in the literature. Our numerical experiments showcase the advantage of using iterative regularization for obtaining more robust and stable approximation of the posterior than unregularized methods.

Limitations of polynomial chaos expansions in the Bayesian solution of inverse problems

Polynomial chaos expansions are used to reduce the computational cost in the Bayesian solutions of inverse problems by creating a surrogate posterior that can be evaluated inexpensively. We show, by analysis and example, that when the data contain significant information beyond what is assumed in the prior, the surrogate posterior can be very different from the posterior, and the resulting estimates become inaccurate. One can improve the accuracy by adaptively increasing the order of the polynomial chaos, but the cost may increase too fast for this to be cost effective compared to Monte Carlo sampling without a surrogate posterior.

Data‐driven model reduction for the Bayesian solution of inverse problems

SummaryOne of the major challenges in the Bayesian solution of inverse problems governed by partial differential equations (PDEs) is the computational cost of repeatedly evaluating numerical PDE models, as required by Markov chain Monte Carlo (MCMC) methods for posterior sampling. This paper proposes a data‐driven projection‐based model reduction technique to reduce this computational cost. The proposed technique has two distinctive features. First, the model reduction strategy is tailored to inverse problems: the snapshots used to construct the reduced‐order model are computed adaptively from the posterior distribution. Posterior exploration and model reduction are thus pursued simultaneously. Second, to avoid repeated evaluations of the full‐scale numerical model as in a standard MCMC method, we couple the full‐scale model and the reduced‐order model together in the MCMC algorithm. This maintains accurate inference while reducing its overall computational cost. In numerical experiments considering steady‐state flow in a porous medium, the data‐driven reduced‐order model achieves better accuracy than a reduced‐order model constructed using the classical approach. It also improves posterior sampling efficiency by several orders of magnitude compared with a standard MCMC method. Copyright © 2014 John Wiley & Sons, Ltd.

Adaptive Bayesian Inference for Discontinuous Inverse Problems, Application to Hyperbolic Conservation Laws

AbstractVarious works from the literature aimed at accelerating Bayesian inference in inverse problems. Stochastic spectral methods have been recently proposed as surrogate approximations of the forward uncertainty propagation model over the support of the prior distribution. These representations are efficient because they allow affordable simulation of a large number of samples from the posterior distribution. Unfortunately, they do not perform well when the forward model exhibits strong nonlinear behavior with respect to its input.In this work, we first relate the fast (exponential) L2-convergence of the forward approximation to the fast (exponential) convergence (in terms of Kullback-Leibler divergence) of the approximate posterior. In particular, we prove that in case the prior distribution is uniform, the posterior is at least twice as fast as the convergence rate of the forwardmodel in those norms. The Bayesian inference strategy is developed in the framework of a stochastic spectral projectionmethod. The predicted convergence rates are then demonstrated for simple nonlinear inverse problems of varying smoothness.We then propose an efficient numerical approach for the Bayesian solution of inverse problems presenting strongly nonlinear or discontinuous systemresponses. This comes with the improvement of the forward model that is adaptively approximated by an iterative generalized Polynomial Chaos-based representation. The numerical approximations and predicted convergence rates of the former approach are compared to the new iterative numericalmethod for nonlinear time-dependent test cases of varying dimension and complexity, which are relevant regarding our hydrodynamics motivations and therefore regarding hyperbolic conservation laws and the apparition of discontinuities in finite time.

Adaptive Construction of Surrogates for the Bayesian Solution of Inverse Problems

The Bayesian approach to inverse problems typically relies on posterior sampling approaches, such as Markov chain Monte Carlo, for which the generation of each sample requires one or more evaluations of the parameter-to-observable map or forward model. When these evaluations are computationally intensive, approximations of the forward model are essential to accelerating sample-based inference. Yet the construction of globally accurate approximations for nonlinear forward models can be computationally prohibitive and in fact unnecessary, as the posterior distribution typically concentrates on a small fraction of the support of the prior distribution. We present a new approach that uses stochastic optimization to construct polynomial approximations over a sequence of measures adaptively determined from the data, eventually concentrating on the posterior distribution. The approach yields substantial gains in efficiency and accuracy over prior-based surrogates, as demonstrated via application to inverse problems in partial differential equations.

A Stochastic Collocation Approach to Bayesian Inference in Inverse Problems

We present an efficient numerical strategy for the Bayesian solution of inverse problems. Stochastic collocation methods, based on generalized polynomial chaos (gPC), are used to construct a polynomial approximation of the forward solution over the support of the prior distribution. This approximation then defines a surrogate posterior probability density that can be evaluated repeatedly at minimal computational cost. The ability to simulate a large number of samples from the posterior distribution results in very accurate estimates of the inverse solution and its associated uncertainty. Combined with high accuracy of the gPC-based forward solver, the new algorithm can provide great efficiency in practical applications. A rigorous error analysis of the algorithm is conducted, where we establish convergence of the approximate posterior to the true posterior and obtain an estimate of the convergence rate. It is proved that fast (exponential) convergence of the gPC forward solution yields similarly fast (exponential) convergence of the posterior. The numerical strategy and the predicted convergence rates are then demonstrated on nonlinear inverse problems of varying smoothness and dimension. AMS subject classifications: 41A10, 60H35, 65C30, 65C50

Stochastic spectral methods for efficient Bayesian solution of inverse problems

We present a reformulation of the Bayesian approach to inverse problems, that seeks to accelerate Bayesian inference by using polynomial chaos (PC) expansions to represent random variables. Evaluation of integrals over the unknown parameter space is recast, more efficiently, as Monte Carlo sampling of the random variables underlying the PC expansion. We evaluate the utility of this technique on a transient diffusion problem arising in contaminant source inversion. The accuracy of posterior estimates is examined with respect to the order of the PC representation, the choice of PC basis, and the decomposition of the support of the prior. The computational cost of the new scheme shows significant gains over direct sampling.

Maximum likelihood estimation of covariance parameters for Bayesian atmospheric trace gas surface flux inversions

This paper introduces a Maximum Likelihood (ML) approach for estimating the statistical parameters required for the covariance matrices used in the solution of Bayesian inverse problems aimed at estimating surface fluxes of atmospheric trace gases. The method offers an objective methodology for populating the covariance matrices required in Bayesian inversions, thereby resulting in better estimates of the uncertainty associated with derived fluxes and minimizing the risk of inversions being biased by unrealistic covariance parameters. In addition, a method is presented for estimating the uncertainty associated with these covariance parameters. The ML method is demonstrated using a typical inversion setup with 22 flux regions and 75 observation stations from the National Oceanic and Atmospheric Administration‐Climate Monitoring and Diagnostics Laboratory (NOAA‐CMDL) Cooperative Air Sampling Network with available monthly averaged carbon dioxide data. Flux regions and observation locations are binned according to various characteristics, and the variances of the model‐data mismatch and of the errors associated with the a priori flux distribution are estimated from the available data.