Simulation of a Gaussian process on $R^2$ with a stationary kernel using the Karhunen-Loève expansion

Question

Assume $X(\omega, t) \sim \mathcal{N}(0, K(\cdot, \cdot))$ is a real-valued, centered Gaussian process on $R^2$, i.e., $X: \Omega \times R^2 \to R$. Let the covariance function of the process be stationary, i.e., $K(t, s) = K(t - s)$, $\forall t, s \in R^2$. More specifically, let it be exponential: $K(t, s) = \exp(-\|t - s\|/L)$ where $\|\cdot\|$ stands for the Euclidean distance, and $L$ denotes the correlation length. Using the Karhunen-Loève expansion, $X$ can be represented as the following summation:

$X(t, \omega) = \sum_{i=1}^\infty \sqrt{\lambda_i} \phi_i(t) \xi_i(\omega)$

where $\xi_i$ are independent, standard Gaussians, and $\phi_i$ and $\lambda_i$ are the eigenfunctions and eigenvalues of the covariance function $K$, which, in our case, can be computed analytically. For practical computations, the expansion above is truncated to keep, say, $n$ terms.

The question is: How to draw samples of the process on $R^2$, i.e., on a two-dimensional grid of fixed points, using the first $n$ eigenfunctions and eigenvalues of $K$? The problem is that the eigenfunctions are somewhat one-dimensional, and I cannot figure out how to properly use them to get realizations of $X$ on a plane. To be honest, the analytical solution that I claimed previously was initially obtained for processes on $R$, so I am not even sure that it applies directly to $R^2$. When I plug in the norms of points' coordinates into $\phi_i$, I get radial patterns, which is presumably not correct.

Also, it is easy to simulate $X$ by constructing the corresponding covariance matrix for a set of points, which boils down to drawing samples from a multivariate Gaussian distribution. I am not interested in this approach.

I would be very grateful for any ideas. Thank you.

Answer 1

%choose a kernel (covariance function)

kernel = 2;
switch kernel
    case 1; k = @(x,y) 1*x'*y; %linear
    case 2; k = @(x,y) exp(-100*(x-y)'*(x-y));
    case 3; k = @(x,y) exp(-1*sqrt((x-y)'*(x-y)));
end

%choose points at which to sample
points = (0:0.02:1)';
[U,V] = meshgrid(points, points);
x=[U(:) V(:)]';
n=size(x,2);

%construct the covariance matrix
C = zeros(n,n);
for i = 1:n
    for j = 1:n
        C(i,j) = k(x(:,i), x(:,j));
    end
end
%sample from the Gaussian process at these points
u = randn(n,1);
[A,S,B]= svd(C);
z = A*sqrt(S)*u;

%plot
figure(2); clf;
Z = reshape(z, sqrt(n), sqrt(n));
surf(U,V,Z);