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

by

Colin Fox

## Abstract

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