# 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.

• Gaussian variograms (at least in one and two dimensions) are problematic because they are too smooth at the origin. At a minimum, include a small nugget term. Generally they should be avoided altogether.
– whuber
Commented Feb 19, 2023 at 16:29

## 2 Answers

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. Adding a larger, positive number to the diagonal can suppress this.

If two feature vectors are identical, or nearly so, this can also cause numerical issues.

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.

• 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. Commented Nov 16, 2015 at 14:51
• @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. Likewise, we know $K$ is symmetric, so you only need to compute one triangle and then fill in the other triangle's values.
– Sycorax
Commented Nov 16, 2015 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.

• Just to chime in and mention what I would do precisely: I would add the number 1e-6 to the diagonal of the correlation matrix. Commented Sep 30, 2021 at 14:55