Bayesian Inverse Problems Research Articles

On the local Lipschitz stability of Bayesian inverse problems

In this note we consider the stability of posterior measures occurring in Bayesian inference w.r.t. perturbations of the prior measure and the log-likelihood function. This extends the well-posedness analysis of Bayesian inverse problems. In particular, we prove a general local Lipschitz continuous dependence of the posterior on the prior and the log-likelihood w.r.t. various common distances of probability measures. These include the total variation, Hellinger, and Wasserstein distance and the Kullback–Leibler divergence. We only assume the boundedness of the likelihoods and measure their perturbations in an Lp-norm w.r.t. the prior. The obtained stability yields under mild assumptions the well-posedness of Bayesian inverse problems, in particular, a well-posedness w.r.t. the Wasserstein distance. Moreover, our results indicate an increasing sensitivity of Bayesian inference as the posterior becomes more concentrated, for example due to more or more accurate data. This confirms and extends previous observations made in the sensitivity analysis of Bayesian inference.

A Bayesian filtering approach to layer stripping for electrical impedance tomography

Layer stripping is a method for solving inverse boundary value problems for elliptic PDEs, originally proposed in the literature for solving the Calderón problem of electrical impedance tomography (EIT), where the data consist of the Neumann-to-Dirichlet operator on the boundary. Defining a tangent–normal coordinate system near the boundary, the data are extended to a family of boundary operators on tangential surfaces inside the body, and it is shown that the operators satisfy a non-linear Riccati type differential equation with respect to the normal coordinate. The layer stripping process consists of a sequence of two alternating steps: the conductivity near the current boundary is estimated from the spatial high-frequency limit of the boundary data, and the boundary operator is propagated through a thin layer further into the domain via the Riccati equation. This way, the unknown conductivity in the interior of the domain is estimated layer by layer starting from the boundary and moving inward. The ill-posedness of the EIT problem manifests itself in such high sensitivity of the backwards Riccati equation to errors in the boundary data to cause the solutions to blow up in finite time, thus requiring regularization. In this article, we formulate the layer stripping process in the framework of Bayesian inverse problems, and we revisit the implementation in the light of Bayesian filtering. More specifically, we recast the related inverse boundary value problem as a state estimation problem, and propose an algorithm for its numerical solution based on ensemble Kalman filtering (EnKF). The new Bayesian layer stripping approach that we propose is quite robust, derivative-free and intrinsically suited for the quantification of uncertainties in the estimate. Furthermore, we show that the algorithm can be extended to realistic data collected by using a finite number of contact electrodes.

Semivariogram methods for modeling Whittle–Matérn priors in Bayesian inverse problems

We present a new technique, based on semivariogram methodology, for obtaining point estimates for use in prior modeling for solving Bayesian inverse problems. This method requires a connection between Gaussian processes with covariance operators defined by the Matérn covariance function and Gaussian processes with precision (inverse-covariance) operators defined by the Green’s functions of a class of elliptic stochastic partial differential equations (SPDEs). We present a detailed mathematical description of this connection. We will show that there is an equivalence between these two Gaussian processes when the domain is infinite—for us, —which breaks down when the domain is finite due to the effect of boundary conditions on Green’s functions of PDEs. We show how this connection can be re-established using extended domains. We then introduce the semivariogram method for estimating the Matérn covariance hyperparameters, which specify the Gaussian prior needed for stabilizing the inverse problem. Results are extended from the isotropic case to the anisotropic case where the correlation length in one direction is larger than another. Finally, we consider the situation where the correlation length is spatially dependent rather than constant. We implement each method in two-dimensional image inpainting test cases to show that it works on practical examples.

Perturbation Monte Carlo Method for Quantitative Photoacoustic Tomography.

Quantitative photoacoustic tomography aims at estimating optical parameters from photoacoustic images that are formed utilizing the photoacoustic effect caused by the absorption of an externally introduced light pulse. This optical parameter estimation is an ill-posed inverse problem, and thus it is sensitive to measurement and modeling errors. In this work, we propose a novel way to solve the inverse problem of quantitative photoacoustic tomography based on the perturbation Monte Carlo method. Monte Carlo method for light propagation is a stochastic approach for simulating photon trajectories in a medium with scattering particles. It is widely accepted as an accurate method to simulate light propagation in tissues. Furthermore, it is numerically robust and easy to implement. Perturbation Monte Carlo maintains this robustness and enables forming gradients for the solution of the inverse problem. We validate the method and apply it in the framework of Bayesian inverse problems. The simulations show that the perturbation Monte Carlo method can be used to estimate spatial distributions of both absorption and scattering parameters simultaneously. These estimates are qualitatively good and quantitatively accurate also in parameter scales that are realistic for biological tissues.

Combining CSEM or gravity inversion with seismic AVO inversion, with application to monitoring of large-scale CO2 injection

A sequential inversion methodology for combining geophysical data types of different resolutions is developed and applied to monitoring of large-scale CO2 injection. The methodology is a two-step approach within the Bayesian framework where lower resolution data are inverted first, and subsequently used in the generation of the prior model for inversion of the higher resolution data. For the application of CO2 monitoring, the first step is done with either controlled source electromagnetic (CSEM) or gravimetric data, while the second step is done with seismic amplitude-versus-offset (AVO) data. The Bayesian inverse problems are solved by sampling the posterior probability distributions using either the ensemble Kalman filter or ensemble smoother with multiple data assimilation. A model-based parameterization is used to represent the unknown geophysical parameters: electric conductivity, density, and seismic velocity. The parameterization is well suited for identification of CO2 plume location and variation of geophysical parameters within the regions corresponding to inside and outside of the plume. The inversion methodology is applied to a synthetic monitoring test case where geophysical data are made from fluid-flow simulation of large-scale CO2 sequestration in the Skade formation. The numerical experiments show that seismic AVO inversion results are improved with the sequential inversion methodology using prior information from either CSEM or gravimetric inversion.

Implicit sampling for hierarchical Bayesian inversion and applications in fractional multiscale diffusion models

This paper presents an implicit sampling method for hierarchical Bayesian inverse problems. A widely used approach for sampling posterior distribution is Markov chain Monte Carlo (MCMC). However, the samples generated by MCMC are usually strongly correlated, which may lead to a small size of effective samples. The implicit sampling method proposed in Chorin and Tu (2009) can generate independent samples and capture inherent non-Gaussian features of the posterior. In implicit sampling method, the posterior samples are generated by a map and distribution around the maximum a posterior (MAP) point. For practical applications, some parameters in prior density are often unknown and the hierarchical Bayesian formulation is used to estimate the MAP point effectively. It is applied to the Bayesian inverse problems of multi-term time fractional diffusion models in heterogeneous media. To effectively capture the heterogeneity effect, we present a mixed generalized multiscale finite element method (mixed GMsFEM) to substantially speed up the Bayesian inversion. Finally, we carry out a few numerical examples to effectively identify various unknown inputs.

Analysis of a multilevel Markov chain Monte Carlo finite element method for Bayesian inversion of log-normal diffusions

We develop the multilevel Markov chain Monte Carlo finite element method (MLMCMC-FEM) to sample from the posterior density of the Bayesian inverse problems. The unknown is the diffusion coefficient of a linear, second-order divergence form, elliptic equation in a bounded, polytopal subdomain of . We provide a convergence analysis with absolute mean convergence rate estimates for the proposed modified MLMCMC-FEM showing in particular error versus work bounds, which are explicit in the discretization parameters. This work generalizes the MLMCMC-FEM algorithm and the error versus work analysis for the uniform prior measure from Hoang et al (2013 Inverse Problems 29), which we also review here, to linear, elliptic, divergence-form PDEs with a log-Gaussian uncertain coefficient and Gaussian prior measure. In comparison to Hoang et al (2013 Inverse Problems 29), we show by mathematical proofs and numerical examples that the unboundedness of the parameter range under Gaussian prior and the non-uniform ellipticity of the forward model require essential modifications in the MCMC sampling algorithm and in the error analysis. The proposed novel multilevel MCMC sampler applies to general Bayesian inverse problems for linear, second order elliptic divergence-form PDEs with log-Gaussian coefficients. It only requires a numerical forward solver with essentially optimal complexity for producing an approximation of the posterior expectation of a quantity of interest within a prescribed accuracy. Numerical examples using independence and pCN samplers are in agreement with our error versus work analysis.

Convergence of Gaussian Process Regression with Estimated Hyper-Parameters and Applications in Bayesian Inverse Problems

This work is concerned with the convergence of Gaussian process regression. A particular focus is on hierarchical Gaussian process regression, where hyper-parameters appearing in the mean and covariance structure of the Gaussian process emulator are a-priori unknown and are learned from the data, along with the posterior mean and covariance. We work in the framework of empirical Bayes, where a point estimate of the hyper-parameters is computed, using the data, and then used within the standard Gaussian process prior to posterior update. We provide a convergence analysis that (i) holds for a given, deterministic function $f$ to be emulated, and (ii) shows that convergence of Gaussian process regression is unaffected by the additional learning of hyper-parameters from data and is guaranteed in a wide range of scenarios. As the primary motivation for the work is the use of Gaussian process regression to approximate the data likelihood in Bayesian inverse problems, we provide a bound on the error introduced in the Bayesian posterior distribution in this context.

Transport Map Accelerated Adaptive Importance Sampling, and Application to Inverse Problems Arising from Multiscale Stochastic Reaction Networks

In many applications, Bayesian inverse problems can give rise to probability distributions which contain complexities due to the Hessian varying greatly across parameter space. This complexity ofte...

An Adaptive Surrogate Modeling Based on Deep Neural Networks for Large-Scale Bayesian Inverse Problems

In Bayesian inverse problems, surrogate models are often constructed to speed up the computational procedure, as the parameter-to-data map can be very expensive to evaluate. However, due to the curse of dimensionality and the nonlinear concentration of the posterior, traditional surrogate approaches (such us the polynomial-based surrogates) are still not feasible for large scale problems. To this end, we present in this work an adaptive multi-fidelity surrogate modeling framework based on deep neural networks (DNNs), motivated by the facts that the DNNs can potentially handle functions with limited regularity and are powerful tools for high dimensional approximations. More precisely, we first construct offline a DNNs-based surrogate according to the prior distribution, and then, this prior-based DNN-surrogate will be adaptively \& locally refined online using only a few high-fidelity simulations. In particular, in the refine procedure, we construct a new shallow neural network that view the previous constructed surrogate as an input variable -- yielding a composite multi-fidelity neural network approach. This makes the online computational procedure rather efficient. Numerical examples are presented to confirm that the proposed approach can obtain accurate posterior information with a limited number of forward simulations.

Kernel Methods for Bayesian Elliptic Inverse Problems on Manifolds

This paper investigates the formulation and implementation of Bayesian inverse problems to learn input parameters of partial differential equations (PDEs) defined on manifolds. Specifically, we stu...

Randomization and Reweighted $\ell_1$-Minimization for A-Optimal Design of Linear Inverse Problems

We consider optimal design of PDE-based Bayesian linear inverse problems with infinite-dimensional parameters. We focus on the A-optimal design criterion, defined as the average posterior variance ...

Sparsity Promoting Hybrid Solvers for Hierarchical Bayesian Inverse Problems

The recovery of sparse generative models from few noisy measurements is an important and challenging problem. Many deterministic algorithms rely on some form of $\ell_1$-$\ell_2$ minimization to combine the computational convenience of the $\ell_2$ penalty and the sparsity promotion of the $\ell_1$. It was recently shown within the Bayesian framework that sparsity promotion and computational efficiency can be attained with hierarchical models with conditionally Gaussian priors and gamma hyperpriors. The related Gibbs energy function is a convex functional, and its minimizer, which is the maximum a posteriori (MAP) estimate of the posterior, can be computed efficiently with the globally convergent Iterated Alternating Sequential (IAS) algorithm [D. Calvetti, E. Somersalo, and A. Strang, Inverse Problems, 35 (2019), 035003]. Generalization of the hyperpriors for these sparsity promoting hierarchical models to a generalized gamma family either yield globally convex Gibbs energy functionals or can exhibit local convexity for some choices for the hyperparameters [D. Calvetti et al., Inverse Problems, 36 (2020), 025010]. The main problem in computing the MAP solution for greedy hyperpriors that strongly promote sparsity is the presence of local minima. To overcome the premature stopping at a spurious local minimizer, we propose two hybrid algorithms that first exploit the global convergence associated with gamma hyperpriors to arrive in a neighborhood of the unique minimizer and then adopt a generalized gamma hyperprior that promotes sparsity more strongly. The performance of the two algorithms is illustrated with computed examples.

Affine Invariant Interacting Langevin Dynamics for Bayesian Inference

We propose a computational method (with acronym ALDI) for sampling from a given target distribution based on first-order (overdamped) Langevin dynamics which satisfies the property of affine invariance. The central idea of ALDI is to run an ensemble of particles with their empirical covariance serving as a preconditioner for their underlying Langevin dynamics. ALDI does not require taking the inverse or square root of the empirical covariance matrix, which enables application to high-dimensional sampling problems. The theoretical properties of ALDI are studied in terms of nondegeneracy and ergodicity. Furthermore, we study its connections to diffusion on Riemannian manifolds and Wasserstein gradient flows. Bayesian inference serves as a main application area for ALDI. In case of a forward problem with additive Gaussian measurement errors, ALDI allows for a gradient-free approximation in the spirit of the ensemble Kalman filter. A computational comparison between gradient-free and gradient-based ALDI is provided for a PDE constrained Bayesian inverse problem.

On the Well-posedness of Bayesian Inverse Problems

The subject of this article is the introduction of a new concept of well-posedness of Bayesian inverse problems. The conventional concept of (Lipschitz, Hellinger) well-posedness in [Stuart 2010, Acta Numerica 19, pp. 451-559] is difficult to verify in practice and may be inappropriate in some contexts. Our concept simply replaces the Lipschitz continuity of the posterior measure in the Hellinger distance by continuity in an appropriate distance between probability measures. Aside from the Hellinger distance, we investigate well-posedness with respect to weak convergence, the total variation distance, the Wasserstein distance, and also the Kullback--Leibler divergence. We demonstrate that the weakening to continuity is tolerable and that the generalisation to other distances is important. The main results of this article are proofs of well-posedness with respect to some of the aforementioned distances for large classes of Bayesian inverse problems. Here, little or no information about the underlying model is necessary; making these results particularly interesting for practitioners using black-box models. We illustrate our findings with numerical examples motivated from machine learning and image processing.

Multilevel Adaptive Sparse Leja Approximations for Bayesian Inverse Problems

Deterministic interpolation and quadrature methods are often unsuitable to address Bayesian inverse problems depending on computationally expensive forward mathematical models. While interpolation may give precise posterior approximations, deterministic quadrature is usually unable to efficiently investigate an informative and thus concentrated likelihood. This leads to a large number of required expensive evaluations of the mathematical model. To overcome these challenges, we formulate and test a multilevel adaptive sparse Leja algorithm. At each level, adaptive sparse grid interpolation and quadrature are used to approximate the posterior and perform all quadrature operations, respectively. Specifically, our algorithm uses coarse discretizations of the underlying mathematical model to investigate the parameter space and to identify areas of high posterior probability. Adaptive sparse grid algorithms are then used to place points in these areas, and ignore other areas of small posterior probability. The points are weighted Leja points. As the model discretization is coarse, the construction of the sparse grid is computationally efficient. On this sparse grid, the posterior measure can be approximated accurately with few expensive, fine model discretizations. The efficiency of the algorithm can be enhanced further by exploiting more than two discretization levels. We apply the proposed multilevel adaptive sparse Leja algorithm in numerical experiments involving elliptic inverse problems in 2D and 3D space, in which we compare it with Markov chain Monte Carlo sampling and a standard multilevel approximation.

Posterior contraction rates for non-parametric state and drift estimation

<p style='text-indent:20px;'>We consider a combined state and drift estimation problem for the linear stochastic heat equation. The infinite-dimensional Bayesian inference problem is formulated in terms of the Kalman–Bucy filter over an extended state space, and its long-time asymptotic properties are studied. Asymptotic posterior contraction rates in the unknown drift function are the main contribution of this paper. Such rates have been studied before for stationary non-parametric Bayesian inverse problems, and here we demonstrate the consistency of our time-dependent formulation with these previous results building upon scale separation and a slow manifold approximation.

Enhancing industrial X-ray tomography by data-centric statistical methods

AbstractX-ray tomography has applications in various industrial fields such as sawmill industry, oil and gas industry, as well as chemical, biomedical, and geotechnical engineering. In this article, we study Bayesian methods for the X-ray tomography reconstruction. In Bayesian methods, the inverse problem of tomographic reconstruction is solved with the help of a statistical prior distribution which encodes the possible internal structures by assigning probabilities for smoothness and edge distribution of the object. We compare Gaussian random field priors, that favor smoothness, to non-Gaussian total variation (TV), Besov, and Cauchy priors which promote sharp edges and high- and low-contrast areas in the object. We also present computational schemes for solving the resulting high-dimensional Bayesian inverse problem with 100,000–1,000,000 unknowns. We study the applicability of a no-U-turn variant of Hamiltonian Monte Carlo (HMC) methods and of a more classical adaptive Metropolis-within-Gibbs (MwG) algorithm to enable full uncertainty quantification of the reconstructions. We use maximum a posteriori (MAP) estimates with limited-memory BFGS (Broyden–Fletcher–Goldfarb–Shanno) optimization algorithm. As the first industrial application, we consider sawmill industry X-ray log tomography. The logs have knots, rotten parts, and even possibly metallic pieces, making them good examples for non-Gaussian priors. Secondly, we study drill-core rock sample tomography, an example from oil and gas industry. In that case, we compare the priors without uncertainty quantification. We show that Cauchy priors produce smaller number of artefacts than other choices, especially with sparse high-noise measurements, and choosing HMC enables systematic uncertainty quantification, provided that the posterior is not pathologically multimodal or heavy-tailed.

Scalable Optimization-Based Sampling on Function Space

Optimization-based samplers such as randomize-then-optimize (RTO) [2] provide an efficient and parallellizable approach to solving large-scale Bayesian inverse problems. These methods solve randomly perturbed optimization problems to draw samples from an approximate posterior distribution. "Correcting" these samples, either by Metropolization or importance sampling, enables characterization of the original posterior distribution. This paper focuses on the scalability of RTO to problems with high- or infinite-dimensional parameters. We introduce a new subspace acceleration strategy that makes the computational complexity of RTO scale linearly with the parameter dimension. This subspace perspective suggests a natural extension of RTO to a function space setting. We thus formalize a function space version of RTO and establish sufficient conditions for it to produce a valid Metropolis-Hastings proposal, yielding dimension-independent sampling performance. Numerical examples corroborate the dimension-independence of RTO and demonstrate sampling performance that is also robust to small observational noise.

Nonlocal TV-Gaussian prior for Bayesian inverse problems with applications to limited CT reconstruction

Bayesian inference methods have been widely applied in inverse problems due to the ability of uncertainty characterization of the estimation. The prior distribution of the unknown plays an essential role in the Bayesian inference, and a good prior distribution can significantly improve the inference results. In this paper, we propose a hybrid prior distribution on combining the nonlocal total variation regularization (NLTV) and the Gaussian distribution, namely NLTG prior. The advantage of this hybrid prior is two-fold. The proposed prior models both texture and geometric structures present in images through the NLTV. The Gaussian reference measure also provides a flexibility of incorporating structure information from a reference image. Some theoretical properties are established for the hybrid prior. We apply the proposed prior to limited tomography reconstruction problem that is difficult due to severe data missing. Both maximum a posteriori and conditional mean estimates are computed through two efficient methods and the numerical experiments validate the advantages and feasibility of the proposed NLTG prior.