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

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;

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?

    Alternatively, this is potentially a numerical issue. I'm not sure how you've concluded that the result is not p.s.d., but, for example, spectral decomposition algorithms will sometimes commit errors due to finite-precision arithmetic, so some of the smallest eigenvalues will be slightly negative (on the order of whatever the error is in your numerical software). These minute errors can be safely ignored in most practical applications.

  2. Testing if an arbitrary covariance function is a valid kernel relies on Mercer's theorem, which is a bit involved for this forum. I'd recommend referring to his research.
  • $\begingroup$ Thanks a lot! It seems to be a numerical issue. I got an error when calculating the cholesky decomposition of K, saying that input is not positive definite. When I check the eigenvalues they're slightly below zero. $\endgroup$ – stedec Nov 16 '15 at 14:51
  • 1
    $\begingroup$ @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. $\endgroup$ – Sycorax says Reinstate Monica Nov 16 '15 at 15:57

I believe I know that example. I worked it in R and got the same problem with singularity. You might try adding some white noise to the model. The squared exponential kernel is very smooth, so if your data points are close together, the covariance matrix goes singular. Recall that at short distances, differentiable functions are approximately linear and the columns of the matrix will be collinear.

To add noise, add a constant to $K(i,i)$, the diagonal element.


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.