Optimization-Based Markov Chain Monte Carlo Methods for Nonlinear Hierarchical Statistical Inverse Problems

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 nonconcave structures in the joint posterior distribution over model parameters and hyperparameters, solving inverse problems in the hierarchical Bayesian setting poses a significant computational challenge. In this work, we 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 [J. M. Bardsley et al., SIAM J. Sci. Comput., 42 (2016), pp. A1317-A1347] for exploring the high- or infinite-dimensional parameter space. We first extend the RTO machinery to the Poisson likelihood and discuss the implementation of RTO in the hierarchical setting. Then, by using RTO either as a proposal distribution in a Metropolis-within-Gibbs update or as a biasing distribution in the pseudomarginal MCMC [C. Andrieu and G. O. Roberts, Ann. Statist., 37 (2009), pp. 697-725], we present efficient sampling tools for hierarchical Bayesian inversion. In particular, the integration of RTO and the pseudomarginal MCMC has sampling performance robust to model parameter dimensions. Numerical examples in PDE-constrained inverse problems and positron emission tomography are used to demonstrate the performance of our methods.