# Estimating the variance of the noise in Gaussian Process prediction

I've been trying to use leave-one-out cross-validation to estimate the $\sigma_n$, the variance of the signal noise when doing prediction according to

$E[f_*] = k_*^T(K+\sigma_n^2I)^{-1}y$ (GPML Equation 2.25)

My questions are:

1. In Chapter 5 of GPML, Rasmussen and Williams suggest using leave one out cross-validation (LOO-CV) to estimate hyperparameters by maximizing the LOO log predictive probability. I think they mean to do this for hyperparameters in the kernel function $k$ -- is this appropriate to use for estimating the $\sigma_n$?

2. If I'm using MSE as my loss function instead of LOO log predictive probability, for some training sets, I get the MSE to be a monotonically increasing function of $\sigma_n$. Does that just mean my training set is just too small? Is there a more appropriate loss function to use?

I presume you mean $\sigma_n$ of the square exponential (now called expodentiated quadratic) kernel:

$k(x,x') = \sigma_n \exp\left(-\frac{(x-x')^2}{2l}\right)$

In which case the LOO negative log marginal likelihood can be used to determine the hyper parameters $\sigma_n$ and $l$. This is done by finding the pair of hyper parameters which maximises the likelihood of the unseen data (the one being left out).

There is a relationship between minimising the MSE and probability of prediction, so I think you should find that you will arrive at the same results either way (although I have not tested this I am making the assumption based on knowledge of fitting linear models and I think the same logic will apply).

What could be happening is that you are being stuck in a local minima. From my experience I have found that there are often three minima in square exponential hyper parameter space. One is when we are overfitting the data, one when we underfit the data and one when it is just right. This relates to balancing the two mane terms of the negative log marginal likelihood:

$\ln P(y|x,θ)=−\frac{1}{2}\ln|K|−\frac{1}{2}y^tK^{−1}y−\frac{N}{2}\ln(2\pi)$

The three components can be seen as balancing the complexity of the GP (to avoid overfit) and the data fit, with a constant on the end. It sounds like your optimisation has found a state where the complexity term is over weighed by the data fit term.

My advice is to perform your optimisation with randomised initialisations of the hyper parameters. Also make sure you don't have any bugs in your code - sometimes they are hard to spot but can cause nightmares (story of my life!)

• Thank you for your detailed response. However, my kernel function actually does not have any hyperparameters. I'm doing learning on amino acid sequences, and the kernel function $k(s,s')$ is the number of amino acids that differ between sequences $s$ and $s'$. – Kevin Yang Apr 18 '15 at 17:47
• @KevinYang Oh Ok in that case what do you mean by $\sigma_n$? Is this meant to be the accuracy of the GP? or is there a single out point variance such as $k(x,x') = \sigma_n |x - x'|$? I might be able to help or amend the answer if you let me know :) – j__ Apr 18 '15 at 18:18
• $\sigma_n$ comes from the Guassian process prediction equations, ie $E[f_*]=k^T_*(K+\sigma_n^2I)^{−1}y$ – Kevin Yang Apr 21 '15 at 0:04
• @KevinYang Ohh that's called jitter and is used to keep K positive definite - it's not really a hyper parameter per say. Jitchol is the function implemented by GPy. It starts by adding a very small constant (the trace times 10^-6) an incrementally multiplies it by 10 until the covariance matrix is positive detonate – j__ Apr 21 '15 at 4:41
• Looking back through my equations, I discovered that I actually do have a hyperparameter in my kernel function. In that case, would I want to optimize that hyperparameter, then go back and increment the jitter until the covariance matrix is positive definite? – Kevin Yang Apr 22 '15 at 18:13