Efficient sampling approaches based on generalized Golub-Kahan methods for large-scale hierarchical Bayesian inverse problems

Uncertainty quantification for large-scale inverse problems remains a challenging task. For linear inverse problems with additive Gaussian noise and Gaussian priors, the posterior is Gaussian but sampling can be challenging, especially for problems with a very large number of unknown parameters (e.g., dynamic inverse problems) and for problems where computation of the square root and inverse of the prior covariance matrix are not feasible. Moreover, for hierarchical problems where several hyperparameters that define the prior and the noise model must be estimated from the data, the posterior distribution may no longer be Gaussian, even if the forward operator is linear. Performing large-scale uncertainty quantification for these hierarchical settings requires new computational techniques. In this work, we consider a hierarchical Bayesian framework where both the noise and prior variance are modeled as hyperparameters. Our approach uses Metropolis-Hastings independence sampling within Gibbs where the proposal distribution is based on generalized Golub-Kahan based methods. We consider two proposal samplers, one that uses a low rank approximation to the conditional covariance matrix and another that uses a preconditioned Lanczos method. Numerical examples from seismic imaging, dynamic photoacoustic tomography, and atmospheric inverse modeling demonstrate the effectiveness of the described approaches.