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

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.

-

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

-

Hello There is a video for mathmonk in youtube which explain how to do that: https://www.youtube.com/watch?v=e7as_wY0hUI&list=HL1361329139&feature=mh_lolz Regards

-
Please pay your attention to the part starting from "Also, ..." in my question. This is exactly what Mathematical Monk is showing in the video. I would like to do the same, but using the decomposition above. –  Ivan Feb 20 '13 at 6:47
Please try to describe as much here as possible in order to make the answer self-contained. Links are fine as support, but they can go stale and then an answer which is nothing more than a link loses its value. –  robjohn Feb 20 '13 at 9:21
I thought that the Mathematical Monk used this composition, as I mentioned above. He obtained the eigenvalues and vectors and did the sampling for finite value of n. –  ibrahim Feb 20 '13 at 10:25
%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); 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); –  ibrahim Feb 22 '13 at 0:20