Simulation of Gaussian Random Fields With One Derivative

Abstract

Circulant embedding provides a fast and exact algorithm for simulating stationary Gaussian random fields on grids. However, the method tends to work poorly for differentiable processes with correlation ranges comparable to the size of the simulation domain. This work describes an algorithm for simulating differentiable stationary Gaussian processes, with the focus on processes possessing exactly one derivative. The algorithm works by first simulating a filtered version of the process using circulant embedding and then recovering the original process from the filtered version. This recovery step is the slowest part of the computation, requiring many solutions of large systems of linear equations. Preconditioned conjugate gradient provides a viable method for obtaining essentially exact solutions to these linear equations. The memory requirements are linear in the number of simulation points and if more than one simulation under a model is needed, subsequent simulations are much faster than the first. The algorithm can be modified to simulate certain differentiable nonstationary intrinsic random functions, to which existing approaches to circulant embedding do not apply. A supplemental file contains code for simulating these intrinsic random functions.