# Kernel matrix is not positive definite

I tried to implement a Gaussian process in octave.

As a starting point I used the algorithm described on page 19 of Rasmussens GP book (http://www.gaussianprocess.org/gpml/).

As a covariance matrix I used the squared exponential function (as it is used in the book as well):

function re = k(x1, x2)
re = exp(-(1/2.0) * abs(x1.-x2).^2);
endfunction

And calculate the covariance matrix (of the training inputs with):

# Calculate covariance matrix
s = size(X);
K = [];
for i = 1:s
for j = 1:i
re = k(X(i), X(j));
K(i,j) = re;
K(j,i) = re;
endfor
endfor

But for some reason the resulting covariance matrix K sometimes is not positive definite (depending on inputs X).

So can anyone tell me what I'm doing wrong here, please? And is there a way to test whether a covariance function results in a positive definite covariance matrix? Since the squared exponential function seems to be a covariance function, I assumed it should create a positive definite matrix.

1. Most research in kernel methods focuses on Mercer kernels, which have two properties: (1) the function is symmetric: $K(x,y)=K(y,x)$ and (2) the function is positive semi-definite (p.s.d.). The Gaussian covariance function is certainly p.s.d., but I can't recall if it is also p.d. -- perhaps you've mistakenly omitted the "semi-" from your mental definition?
• @stedec Point of interest: it looks like you're computing the diagonal of the kernel matrix manually (i.e. visiting $K_{ii}$ for each $i$) -- but in fact, we know that $x_i-x_i=0,$ so you could save $s$ computations by just setting the diagonal to the appropriate constant and not looping over it. – Sycorax says Reinstate Monica Nov 16 '15 at 15:57
To add noise, add a constant to $K(i,i)$, the diagonal element.