Richard A. Norton
Langevin and Hamiltonian proposals are matrix splittings

Department of Physics
University of Otago
PO Box 56
Dunedin 9054
New Zealand
richard.norton@otago.ac.nz
Colin Fox

Stationary iterative methods based on matrix splittings, such as Gauss-Seidel and symmetric successive over-relaxation (SSOR), are well understood for solving linear systems

$$A x = b$$

where $A \in \mathbb{R}^{d\times d}$ and $x, b \in \mathbb{R}^d$ . The matrix $A$ is split $A = M-N$ and the solution is computed by solving

$$M x^{k+1} = N x^{k} + b.$$

The convergence rate is determined by the spectral radius of $M^{-1} N$ , but can be improved by using a polynomial acceleration technique, such as Chebyshev acceleration.

We investigate how this technology can be transferred to the problem of constructing efficient Markov chain Monte Carlo (MCMC) methods to sample from high dimensional Gaussian distributions,

$$\pi(x) \propto \exp \left( -\frac{1}{2} x^T A x + b^T x \right)$$

where the matrix $A \in \mathbb{R}^{d\times d}$ is the precision matrix for the distribution. We can use exactly the same splitting for a MCMC method as for a linear solver $A = M-N$ to obtain an iterative procedure

$$M x^{k+1} = N x^{k} + b + \nu^k,$$

where $\nu^k$ is an independent sample from $N(0,M^T+N)$ . This iterative procedure computes a Markov chain that will converge in distribution to a sample from $\pi$ if the corresponding iterative solver also converges. The rate of convergence is again dependent on the spectral radius of $M^{-1} N$ .

The Gibbs sampling algorithm corresponds to the Gauss-Seidel and SSOR iterative solvers, as studied by Fox and Parker 2013.

Other sampling methods also correspond to matrix splittings, and we use this correspondence to help us understand the convergence properties of these sampling methods. For example Crank-Nicolson algorithms, described in Cotter et al. 2013, can be expressed as matrix splittings.

The Metropolis-Hastings algorithm is a popular method for constructing MCMC methods, where the accept/reject step in the algorithm is used because a proposal must be corrected so that the Markov chain converges to the desired target distribution. For Metropolis-adjusted Langevin (MALA) and Hamiltonian/hybrid Monte Carlo (HMC) proposals we can express the proposal as a matrix splitting of an ``effective'' precision matrix $\tilde{A} = M-N$ . The difference between the target precision matrix $A$ and the proposal precision matrix $\tilde{A}$ is what the Metropolis-Hastings algorithm corrects. We explore the consequences of the matrix splitting correspondence for these proposals, in particular the opportunity for using polynomial acceleration to construct more efficient sampling algorithms.

• C. Fox and A. Parker, ``Convergence in variance of Chebyshev accelerated Gibbs samplers", SIAM J. Sci. Comput., 2013.

• S. L. Cotter, G. O. Roberts, A. M. Stuart and D. White, ``MCMC methods for functions: modifying old algorithms to make them faster'', Statistical Science, 2013, Vol. 28, No. 3, pp. 424-446.