Accelerated Gibbs sampling of normal distributions using matrix splittings and polynomials

TL;DR

This paper links Gibbs sampling for multivariate normals to iterative linear solvers, showing how matrix splittings and polynomial acceleration can significantly improve convergence rates in high-dimensional Bayesian problems.

Contribution

It establishes a theoretical connection between Gibbs sampling and matrix iterative methods, enabling the use of polynomial acceleration techniques to enhance sampling efficiency.

Findings

Polynomial acceleration improves convergence rates

Matrix splitting methods can be applied to Gibbs sampling

Numerical examples demonstrate practical benefits

Abstract

Standard Gibbs sampling applied to a multivariate normal distribution with a specified precision matrix is equivalent in fundamental ways to the Gauss-Seidel iterative solution of linear equations in the precision matrix. Specifically, the iteration operators, the conditions under which convergence occurs, and geometric convergence factors (and rates) are identical. These results hold for arbitrary matrix splittings from classical iterative methods in numerical linear algebra giving easy access to mature results in that field, including existing convergence results for antithetic-variable Gibbs sampling, REGS sampling, and generalizations. Hence, efficient deterministic stationary relaxation schemes lead to efficient generalizations of Gibbs sampling. The technique of polynomial acceleration that significantly improves the convergence rate of an iterative solver derived from a \emph{symmetric} matrix splitting may be applied to accelerate the equivalent generalized Gibbs sampler. Identicality of error polynomials guarantees convergence of the inhomogeneous Markov chain, while equality of convergence factors ensures that the optimal solver leads to the optimal sampler. Numerical examples are presented, including a Chebyshev accelerated SSOR Gibbs sampler applied to a stylized demonstration of low-level Bayesian image reconstruction in a large 3-dimensional linear inverse problem.