Randomize-then-optimize for sampling and uncertainty quantification in electrical impedance tomography

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.