Multigrid Monte Carlo Revisited: Theory and Bayesian Inference

TL;DR

This paper revisits the Multigrid Monte Carlo algorithm, demonstrating its potential to efficiently sample high-dimensional Gaussian random fields, including non-stationary cases, with proven grid-size-independent convergence and optimality in Bayesian inverse problems.

Contribution

It establishes a new convergence theory for MGMC applicable to non-stationary fields and extends multigrid methods to Bayesian inverse problems with low-rank updates.

Findings

MGMC can significantly accelerate sampling in spatial statistics.

Theoretical proof of grid-size-independent convergence for MGMC.

Numerical results show MGMC outperforms alternative methods in Bayesian settings.

Abstract

Gaussian random fields play an important role in many areas of science and engineering. In practice, they are often simulated by sampling from a high-dimensional multivariate normal distribution, which arises from the discretisation of a suitable precision operator. Existing methods such as Cholesky factorization and Gibbs sampling become prohibitively expensive on fine meshes due to their high computational cost. In this work, we revisit the Multigrid Monte Carlo (MGMC) algorithm developed by Goodman & Sokal (Physical Review D 40.6, 1989) in the quantum physics context. While the authors of this paper conclude that MGMC does not overcome critical slowing down in simulations of field theories near phase transitions, we demonstrate here that it has the potential to significantly accelerate sampling in spatial statistics. The class of Gaussian Random Fields we consider includes those with Mat\'{e}rn covariance, but is more general in that it also allows for non-stationary covariance functions. To show that MGMC can overcome the limitation of existing methods, we establish a grid-size-independent convergence theory based on the link between linear solvers and samplers for multivariate normal distributions, drawing on standard multigrid convergence arguments. We then apply this theory to linear Bayesian inverse problems. This application is achieved by extending the standard multigrid theory to operators with a low-rank perturbation. Moreover, we develop a novel bespoke random smoother which takes care of the low-rank updates that arise in constructing posterior moments. In particular, we prove that Multigrid Monte Carlo is algorithmically optimal in the limit of the grid-size going to zero. Numerical results support our theory, demonstrating that Multigrid Monte Carlo can be significantly more efficient than alternative methods when applied in a Bayesian setting.