A Conjugate Direction Sampler for Normal Distributions, with a Few Computed Examples

Colin Fox


Gaussian models and processes are common in statistics, particularly spatial statistics, being convenient from both computational and theoretical viewpoints. Efficient  algorithms for sampling from Gaussian Markov random fields exploit sparseness of the precision matrix within a Cholesky factorization, or use circulant structure to allow diagonalization by Fourier transforms.

I present an alternative, sequential, algorithm derived from the conjugate gradient (CG) optimization algorithm. CG has the remarkable property of minimizing a quadratic form exactly in a finite number of steps while requiring storage of only two state vectors. The conjugate direction (CD) sampler has the analogous property generating independent (or exact) samples after a fixed number of steps (equal to the dimension of the state space) while requiring storage of only two state vectors, and a third auxiliary vector. Within the sampler one needs to operate by the precision matrix but there is no need to store the matrix or factorize it. Hence the conjugate direction sampler is useful for drawing samples in high dimensional problems where forming the precision matrix is impractical or inconvenient.

When implemented using finite precision arithmetic the CD sampling algorithm terminates  incorrectly due to loss of conjugacy arrising from ill-conditioning of the precision matrix, or when eigenvalues are repeated. To some degree these issues may be treated by use of a preconditioning matrix.  In a computed example the maximum sample dimension for which the CD sampling algorithm  correctly terminates is $10^5$.