Multigrid Monte Carlo Revisited: Theory and Bayesian Inference

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é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.

💡 Research Summary

This paper revisits the Multigrid Monte Carlo (MGMC) algorithm originally introduced by Goodman and Sokal in 1989 and demonstrates its relevance for modern spatial statistics and Bayesian inference. The authors focus on the problem of sampling high‑dimensional Gaussian random fields (GRFs) that arise from discretising a precision operator A on a fine mesh. Traditional approaches such as Cholesky factorisation scale poorly (O(n¹⁺ζ) memory and time) and suffer from numerical bias for very large systems, while stationary iterative samplers (Gibbs, random‑walk Metropolis, preconditioned Crank–Nicolson, Langevin, HMC) exhibit critical slowing down as the grid spacing h→0, leading to exploding autocorrelation times.

The key insight exploited in MGMC is the equivalence between a Gibbs sampler for a multivariate normal N(μ, A⁻¹) and a Gauss–Seidel iteration for solving Ax = f, differing only by an additive zero‑mean noise term. By embedding this stochastic iteration within a multigrid hierarchy—random smoothing on fine levels and coarse‑grid corrections—the algorithm targets low‑frequency components of the error distribution, accelerating convergence of both the mean and the covariance. The authors extend the classical multigrid convergence theory (which normally guarantees only mean convergence) to the stochastic setting, proving that the first two moments converge exponentially at a rate independent of the mesh size. This yields grid‑size‑independent autocorrelation decay and root‑mean‑square convergence of Monte‑Carlo estimators.

A major contribution is the treatment of linear Bayesian inverse problems, where the precision matrix takes the form A + U Uᵀ, i.e., a low‑rank perturbation of a differential operator. The paper introduces a bespoke random smoother that exactly accounts for the low‑rank update while preserving the optimal O(n) computational complexity. The authors prove that MGMC remains algorithmically optimal in the continuum limit (h→0), meaning that the cost per independent sample grows linearly with the number of degrees of freedom, matching the best possible complexity for such problems.

Numerical experiments cover a range of test cases: stationary Matérn fields in two and three dimensions, non‑stationary covariance structures, and linear Bayesian inverse problems (e.g., parameter estimation from noisy observations). MGMC is benchmarked against direct Cholesky sampling, FFT‑based SPDE methods, and low‑rank Cholesky approximations. Results consistently show that MGMC achieves the same statistical accuracy with 5–10× less wall‑clock time, and its autocorrelation times remain essentially constant as the mesh is refined. Parallel scalability is demonstrated using domain‑decomposition techniques, confirming that MGMC can efficiently exploit distributed‑memory architectures.

In summary, the paper provides a rigorous theoretical foundation for MGMC, establishes grid‑independent convergence of both mean and covariance, extends the method to low‑rank Bayesian settings, and validates its practical superiority through extensive experiments. It positions MGMC as a versatile, scalable alternative for sampling Gaussian random fields in high‑resolution spatial statistics and Bayesian inference, and suggests future work on non‑Gaussian targets, algebraic multigrid extensions, and hybrid MCMC schemes.