Fast Sampling in a Linear-Gaussian Inverse Problem

We solve the inverse problem of deblurring a pixelized image of Jupiter using regularized decon-volution and by sample-based Bayesian inference. By efficiently sampling the marginal posterior distribution for hyperparameters, then the full conditional for the deblurred image, we find that we can evaluate the posterior mean faster than regularized inversion, when selection of the regularizing parameter is considered. To our knowledge, this is the first demonstration of sampling and inference being computed in less time than regularized inversion, in a significant inverse problem. Comparison to random-walk Metropolis Hastings and block Gibbs Markov chain Monte Carlo shows that marginal then conditional sampling also outperforms these more common sampling algorithms. The asymptotic cost of an independent sample is one linear solve, implying that, when problem-specific computations are feasible, sample-based Bayesian inference may be performed directly over function spaces when that limit exists.