Fast and Exact Simulation of Stationary Gaussian Processes through Circulant Embedding of the Covariance Matrix

Abstract

Geostatistical simulations often require the generation of numerous realizations of a stationary Gaussian process over a regularly meshed sample grid $\Omega$. This paper shows that for many important correlation functions in geostatistics, realizations of the associated process over $m+1$ equispaced points on a line can be produced at the cost of an initial FFT of length $2m$ with each new realization requiring an additional FFT of the same length. In particular, the paper first notes that if an $(m+1)\times(m+1) $ Toeplitz correlation matrix R can be embedded in a nonnegative definite $2M\times2M$ circulant matrix S, exact realizations of the normal multivariate $y \sim {\cal N}(0,R)$ can be generated via FFTs of length $2M$. Theoretical results are then presented to demonstrate that for many commonly used correlation structures the minimal embedding in which $M = m$ is nonnegative definite. Extensions to simulations of stationary fields in higher dimensions are also provided and illustrated.