Johnathan M. Bardsley

Arvind K. Saibaba, Johnathan Bardsley, D. Andrew Brown et al.

Hierarchical models in Bayesian inverse problems are characterized by an assumed prior probability distribution for the unknown state and measurement error precision, and hyper-priors for the prior parameters. Combining these probability models using Bayes' law often yields a posterior distribution that cannot be sampled from directly, even for a linear model with Gaussian measurement error and Gaussian prior. Gibbs sampling can be used to sample from the posterior, but problems arise when the dimension of the state is large. This is because the Gaussian sample required for each iteration can be prohibitively expensive to compute, and because the statistical efficiency of the Markov chain degrades as the dimension of the state increases. The latter problem can be mitigated using marginalization-based techniques, but these can be computationally prohibitive as well. In this paper, we combine the low-rank techniques of Brown, Saibaba, and Vallelian (2018) with the marginalization approach of Rue and Held (2005). We consider two variants of this approach: delayed acceptance and pseudo-marginalization. We provide a detailed analysis of the acceptance rates and computational costs associated with our proposed algorithms, and compare their performances on two numerical test cases---image deblurring and inverse heat equation.

Kevin T. Joyce, Johnathan M. Bardsley, Aaron Luttman

The point spread function (PSF) of a translation invariant imaging system is its impulse response, which cannot always be measured directly. This is the case in high energy X-ray radiography, and it must be estimated from images of calibration objects indirectly related to the impulse response. When the PSF is assumed to have radial symmetry, it can be estimated from an image of an opaque straight edge. We use a non-parametric Bayesian approach, where the prior probability density for the PSF is modeled as a Gaussian Markov random field and radial symmetry is incorporated in a novel way. Markov Chain Monte Carlo posterior estimation is carried out by adapting a recently developed improvement to the Gibbs sampling algorithm, referred to as partially collapsed Gibbs sampling. Moreover, the algorithm we present is proven to satisfy invariance with respect to the target density. Finally, we demonstrate the efficacy of these methods on radiographic data obtained from a high-energy X-ray diagnostic system at the U.S. Department of Energy's Nevada National Security Site.

Johnathan Bardsley, Tiangang Cui

In many hierarchical inverse problems, not only do we want to estimate high- or infinite-dimensional model parameters in the parameter-to-observable maps, but we also have to estimate hyperparameters that represent critical assumptions in the statistical and mathematical modeling processes. As a joint effect of high-dimensionality, nonlinear dependence, and non-concave structures in the joint posterior posterior distribution over model parameters and hyperparameters, solving inverse problems in the hierarchical Bayesian setting poses a significant computational challenge. In this work, we aim to develop scalable optimization-based Markov chain Monte Carlo (MCMC) methods for solving hierarchical Bayesian inverse problems with nonlinear parameter-to-observable maps and a broader class of hyperparameters. Our algorithmic development is based on the recently developed scalable randomize-then-optimize (RTO) method [4] for exploring the high- or infinite-dimensional model parameter space. By using RTO either as a proposal distribution in a Metropolis-within-Gibbs update or as a biasing distribution in the pseudo-marginal MCMC [2], we are able to design efficient sampling tools for hierarchical Bayesian inversion. In particular, the integration of RTO and the pseudo-marginal MCMC has sampling performance robust to model parameter dimensions. We also extend our methods to nonlinear inverse problems with Poisson-distributed measurements. Numerical examples in PDE-constrained inverse problems and positron emission tomography (PET) are used to demonstrate the performance of our methods.