Randomize-Then-Optimize for Sampling and Uncertainty Quantification in Electrical Impedance Tomography

Abstract

Abstract. In a typical inverse problem, a spatially distributed parameter in a physical model is estimated from measurements of model output. Since measurements are stochastic in nature, so is any parameter estimate. Moreover, in the Bayesian setting, the choice of regularization corresponds to the definition of the prior probability density function, which in turn is an uncertainty model for the unknown parameters. For both of these reasons, significant uncertainties exist in the solution of an inverse problem. Thus to fully understand the solution, quantifying these uncertainties is important. When the physical model is linear and the error model and prior are Gaussian, the posterior density function is Gaussian with a known mean and covariance matrix. However, the electrical impedance tomography inverse problem is nonlinear, and hence no closed form expression exists for the posterior density. The typical approach for such problems is to sample from the posterior and then use the samples to compute statistics (such as the mean and variance) of the unknown parameters. Sampling methods for electrical impedance tomography have been studied by various authors in the inverse problems community. However, up to this point the focus has been on the development of increasingly sophisticated implementations of the Gibbs sampler, whose samples are known to converge very slowly to the correct density for large-scale problems. In this paper, we implement a recently developed sampling method called randomize-then-optimize (RTO), which provides nearly independent samples for each application of an appropriate numerical optimization algorithm. The sample density for RTO is not the posterior density, but RTO can be used as a very effective proposal within a Metropolis-Hastings algorithm to obtain samples from the posterior. Here our focus is on implementing the method on synthetic examples from electrical impedance tomography, and we show that it is both computationally efficient and provides good results. We also compare RTO performance with the Metropolis adjusted Langevin algorithm and find RTO to be much more efficient. 1. Introduction. Solving an inverse problem requires estimating unknown parameters in a physical model from measurements of model output. The unknown parameter vector typically corresponds to a numerical discretization of a continuously defined, spatially distributed parameter in the physical model and hence is high-dimensional. This leads to an overparameterized statistical model, and hence estimates have large variance or, in other words, are unstable with respect to perturbations in the observations. Regularization methods [36], which have been studied for decades in the inverse problems literature, are the typical fix for such instability.Over the past ten years or so, an increasing number of researchers in inverse problems have drawn the connection between regularized solutions of inverse problems and Bayesian statistics [20]. Since the measured data in an inverse problem will always contain noise, a ...

A nonlinear approach to difference imaging in EIT; assessment of the robustness in the presence of modelling errors

“…In the case of the selected regularization, the feasible set of parameters is in the same range for both absolute and difference imaging, and all the estimates are relatively tolerant with respect to the choices of parameters, given that the parameters are in the feasible range-these observations were made based on simulations, but they can be justified based on the Bayesian interpretations of the prior models. For the Bayesian interpretation of the parameters of the model (22), see [59,60]. For a Bayesian intepretation of the TV model ( 23) and an approximation for the selection of the parameter α based on the expected magnitude of the changes in f, see [61].…”

Efficient Krylov subspace methods for uncertainty quantification in large Bayesian linear inverse problems

“…The use of a low‐rank surrogate of H Q has also been explored in References 9,10 and is similar to Method 1 (cf, Section 4.1) that we propose. Other approaches to sampling from the posterior distribution include randomize‐then‐optimize 30,31 and randomized MAP approach 32 . However, none of these methods can handle the case where Q −1 or Q −1/2 is not available.…”

Uncertainty quantification in large Bayesian linear inverse problems using Krylov subspace methods

“…The use of a low-rank surrogate of H Q has also been explored in [10,11] and is similar to Method 1 (c.f., subsection 4.1) that we propose. Other approaches to sampling from the posterior distribution include randomize-then-optimize (RTO) [5,4] and randomized MAP approach [45]. However, none of these methods can handle the case where Q −1 or Q −1/2 are not available.…”