I am trying to implement the Karhunen Loeve expansion for a 1-D Gaussian random field with a square-exponential kernel.

Specifically, I know that a Gaussian process has a KL expansion $\hat{U}=\sum_{i-1}^{N_{KL}} \sqrt{\lambda_i}u_i(\mathbf{x})\xi_i(\theta) $, where $\xi_i$ are independent, center-normalized Gaussian variables, and $\lambda_i$ and $u_i$ are the eigenvalue/eigenfunction solutions to the Fredholm equation $\int_\Theta K(\mathbf{x},\mathbf{x}')u_i(\mathbf{x}')d\mathbf{x}'=\lambda_iu_i(\mathbf{x})$.

For the case of a square-exponential kernel of the form $K(\mathbf{x},\mathbf{x}')=\sigma^2\exp(\frac{-(x-x')^2}{2l^2})$, This source gives the solution $\lambda_k=\sqrt{\frac{2a}{A}}B^k$ and $u_k(x)=\exp(-(c-a)x^2)H_k(\sqrt{2c}x)$, where $H_k$ are $k$th order Hermite polynomials, $a^{-1}=4\sigma^2$, $b^{-1}=2l^2$, $c=\sqrt{a^2+2ab}$, $A=a+b+c$, and $B=b/A$.

I am now trying to apply this to generate a Gaussian process for a domain length of 1000 units, with a grid spacing of 1 unit, and correlation length of 20 units. Because random variables in a KL expansion are drawn from a unit-normal distribution, I interpreted that as $sigma^2=1$. My MATLAB code for this is below:

% define domain:
x = linspace(0,999,1000);

% define statistical properties:
correlation_length = 20;
variance = 1;
mean = 0;

% helper variables for eigenvalues/functions:
a = (4*variance.^2).^-1;
b = (2*correlation_length.^2).^-1;
c = sqrt(a.^2 + 2*a*b);
A = a + b + c;
B = b/A;

% set arbitrary truncation point to 20 eigenfunctions:
Nkl = 20;
% plot 4 realizations: 
U = zeros(4,1000);
for i = 1:4
   for k = 1:Nkl 
       U(i,:) = U(i,:) + sqrt(kth_eval(k, a, A, B)) .* ...
           kth_efunc(k, a, c, x) * normrnd(0,1);
plot(x, U);

% eigen value/vector calculations:
function e_val = kth_eval(k, a, A, B)
    % return scalar value of kth eigenvalue
    e_val = sqrt(2*a/A) * B.^k;

function e_func = kth_efunc(k, a, c, x)
    % return 1-D vector of values taken by the kth eigenfunction 
    % over the domain x
    e_func = exp(-(c-a)*x.^2).*hermiteH(k,sqrt(2*c)*x);

This gives me an output with the following appearance:


Which seems incorrect to me, as I want to get a process that is statistically similar to the following:

evenson et al

With this code, I can also plot the $k=1,2,3$ eigenfunctions, giving the plot:

bad eigen fns

1 has a plot of the first $k=0,1,2$ eigenfunctions, as shown below for $a=1, b=3$:

first 3 eigenfunctions

When I plot the first 3 eigenfunctions with my code with the same a and b, I do get a visually similar results:

similar to book

However, even this still generates problems for me, as when I attempt to realize a Gaussian process as with the code above and normalized random variables, I seem to just get the same symmetric shape re-scaled by a random factor, rather than random shapes that are statistically similar:

too symmetric

If anyone knows what's causing the discrepancy between the results of my code and what I am expecting, I would be grateful to know.


I've had issue with this -- please see the errata of the reference. The eigenfunctions are not scaled in the original version of the book. This errata should solve it.

| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.