# Fastest way of sampling multivariable Guassian with covariance matrix that is circulant.

A common problem in statistics is to compute sample vectors from a multivariate Gaussian distribution with zero mean and a given covariance matrix $A$. A canonical approach to the problem is to compute vectors of the form $y = Sz$, where $S$ is the Cholesky factor or square root of $A$, and $z$ is a standard normal vector. When $A$ is large, such an approach becomes computationally expensive.

Is there a fast way (below $O(N^3)$) to sample a multivariable Gaussian, given the condition that the covariance matrix is circulant?

Here are two different approaches that can be used to solve this problem:

There's a very large literature in geostatistics on exact and approximate methods for computing Markov random fields according to multivariate normal distributions with covariance matrices which have circulant (in 1-D) or block-circulant (in 2-D or 3-D) structure.

In practice, the covariance matrix is typically Toeplitz, but you embed it in a larger matrix of size $2n$ by $2n$ to produce a circulant matrix. This extra factor of $2$ disappears when we consider an $O(n\log(n))$ time algorithm.

The basic idea is to convert the autocovariance into power spectral density by using the fast Fourier transform and then generate a random field with the desired power spectral density by inverting the PSD and applying a random phase at each spatial frequency. You then invert the Fourier transform to produce your final vector.

In the 1-D case, you can generate a random vector with the desired covariance in $O(n \log(n))$ time. This increases to $O(n^{2}\log(n)^2)$ time in 2-D for an $n$ by $n$ field.

See for example:

C. V. Deutsch and A. G. Journel. GSLIB: Geostatistical Software Library and User's Guide (Applied Geostatistics) 2nd Edition. Oxford University Press, 1997.

This second approach might be of interest if you want or need the Cholesky factorization of your covariance matrix. Many papers have been written on fast algorithms for the Cholesky factorization of positive definite Toeplitz matrices. The Bareiss algorithm is a classic in this area:

Bareiss, Erwin H. "Numerical solution of linear equations with Toeplitz and vector Toeplitz matrices." Numerische Mathematik 13.5 (1969): 404-424.

This algorithm can compute a Cholesky factorization of your positive definite Toeplitz covariance matrix in $O(n^{2})$ time. More recent research has shown that this algorithm is numerically stable (unlike the more well-known Levinson-Durbin algorithm.)

Once you have the Cholesky factorization of your Covariance matrix, you can generate a multivariate normal sample using that Cholesky factorization in $O(n^{2})$ time from a vector of $n$ iid N(0,1) random numbers.