Fast and exact simulation of stationary Gaussian processes through circulant embedding of the covariance matrix

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) x (m + 1) Toeplitz correlation matrix R can be embedded in a nonnegative definite 2M x 2M circulant matrix S, exact realizations of the normal multivariate y similar to N(O, 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.