I am trying to optimize the hyperparameters for a Gaussian process. As a starter I choose the squared exponential function for covariance where iI have to optimize 3 parameters $\sigma_f$, $\sigma_n$ and the length parameter $l$.

$$k_y(x_p,x_q) = \sigma^2_f \exp\left(-\frac{1}{2l^2}(x_p-x_q)^2\right) + \sigma^2_n\delta_{pq}$$

As described in [1] I try to maximize the marginal likelihood using the parameters. To achieve this I use Commons Math optimizers and compute the gradients of the marginal likelihood function.

Everything works fine except for a tiny example if I only use $ \sigma_f $ and $ l$ for the optimization. As soon as I try to optimize $ \sigma_n$ as well, I run into numerical problems and the the covariance matrix of the process is NOT positive-definite anymore after a few iterations of the optimizer and thus produces INFINITY or NAN errors within the optimizer.

Can anyone explain that behavior and the connection to the hyperparameter $ \sigma_n$?

[1] Gaussian Processes for Machine Learning Carl Edward Rasmussen and Christopher K. I. Williams The MIT Press, 2006. ISBN 0-262-18253-X.

  • $\begingroup$ What do you mean by "tiny example"? Small sample size? $\endgroup$
    – jbowman
    Commented Jun 8, 2012 at 13:26
  • $\begingroup$ Yes i just evaluate the example: y = {0.55*-3.0, 0.55*-2.0, 0.55*-0.6, 0.55*0.4, 0.55*1.0, 0.55*1.6} X = {-1.5, -1.0, -.75, -0.4, -0.25 ,0}; $\endgroup$
    – Andreas
    Commented Jun 8, 2012 at 13:41
  • $\begingroup$ Could you try this again with a sample of size 12? And add a little fuzz to the x-values, so they aren't quite so close to being on a grid... I'm thinking your small sample size might, but only might, be the problem. $\endgroup$
    – jbowman
    Commented Jun 8, 2012 at 17:41
  • $\begingroup$ I had the Ininity/NAN problem as well for a long time until I switched to using SVD instead of the native cholesky that GPML use. The problem is that with a smooth kernel like this the singular values die down pretty quickly. SVD can alleviate problem to some level if you ignore eigen values below a certain threshold. $\endgroup$
    – sachinruk
    Commented Nov 15, 2013 at 4:25
  • $\begingroup$ It is very common to get non semidefinite covariance matrix due to numerical issue. you then need to increase the lower bound of noise level for your optimization. $\endgroup$
    – user140868
    Commented Dec 2, 2016 at 12:34

1 Answer 1


It would be a good idea to get the optimisation code to print out the hyper-parameters each time it performs a function evaluation. Usually you can work out what is going wrong once you know what the model decides about the hyper-parameter values.

I suspect what is happening is that the model has decided that an essentially linear classifier would be best, and has set $\sigma_n$ to zero, as the noise component is not seen as necessary (which is a shame as it adds a ridge to the covariance matrix, which helps ensure that it is p.d.). In that case, only the SEiso bit is used, so $\sigma_f$ will probably be much larger than $\sigma_n$, however to make a linear classifier it will try and make $l$ as large as possible, which seems to end up resulting in numerical problems when evaluating the bit inside the exponential. I'm a pretty heavy user of GPML and have seen this a fair bit. One solution is to limit the magnitudes of the logarithms of the hyper-parameters during the search (which is equivalent to having a hyper-prior on the hyper-parameters), which tends to prevent this from happening. If you print out the values of the hyper-parameters, then the last ones that get printed before it goes "bang" will give you a good idea where to place the limits. This tends not to affect performance very much, the generalisation error at such points in hyper-parameter tends to be fairly flat, which causes gradient descent methods to take large steps that put you far enough from the origin that you get numerical accuracy problems.

In short print out the hyper-parameter values at each step, whenever you run into numerical issues in model selection.

  • $\begingroup$ Why should i limit the magnitudes if the logarithms instead of the magnitudes themselfes? As far as i know i never use any log values of the hyperparameters at all and just use them in the squared Exponential formula. $\endgroup$
    – Andreas
    Commented Jun 11, 2012 at 13:40
  • 1
    $\begingroup$ It is convenient to optimise parameters that are strictly positive by taking logs so that the search is unconstrained and then convert back again. If you are using GPML it will already be doing this for you automatically and the hyper-parameters you see are actually the logs IIRC. $\endgroup$ Commented Jun 11, 2012 at 17:47
  • $\begingroup$ So i input them as logs to the algorithm -> within i convert them calculate the gradients and i return the log of the gradients? $\endgroup$
    – Andreas
    Commented Jun 14, 2012 at 15:11
  • $\begingroup$ yes, see this paper jmlr.csail.mit.edu/papers/volume8/cawley07a/cawley07a.pdf (bottom of page 848) for an example. $\endgroup$ Commented Jun 14, 2012 at 18:39
  • 1
    $\begingroup$ If the gradient in log-hyper-parameter space is negative, that is no problem as it doesn't matter if the log-hyper-parameters are negative, once the exponential is taken the hyper-paremeters themselves are guaranteed to be positive as exp(x) > 0 for all x. $\endgroup$ Commented Jun 15, 2012 at 12:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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