Simulating Gaussian processes on a lattice

How can I simulate a Gaussian Processes on a lattice with a given covariance?

The general (non-lattice) case is given in historical overview in Liu et al. (2019), but in this notebook we are interested in specialising a little. Following the introduction in Dietrich and Newsam (1993), let’s say we wish to generate a stationary Gaussian process \(Y(x)\) on a points \(\Omega\). \(\Omega=(x_0, x_1,\dots, x_m)\).

Stationary in this context means that the covariance function \(r\) is translation-invariance and depend only on distance, so that it may be given \(r(|x|)\). Without loss of generality, we assume that \(\mathbb E[Y(x)]=0\) and \(\var[Y(x)]=1\).

The problem then reduces to generating a vector \(\vv y=(Y(x_0), Y(x_1), \dots, Y(x_m) )\sim \mathcal{N}(0, R)\) where \(R\) has entries \(R[p,q]=r(|x_p-x_q|).\)

Note that if \(\mathbb \varepsilon\sim\mathcal{N}(0, I)\) is an \(m+1\)-dimensional normal random variable, and \(AA^T=R\), then \(\vv y=\mm A \vv \varepsilon\) has the required distribution.