Vtyjkyuk

Assume $X(omega, t) sim mathcalN(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 sqrtlambda_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.

%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);