  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. To show that MGMC can
overcome these issues, 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 theory. 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.
