# Gaussian Process: why is my data only explained by noise?

I'm trying to explain a very simple 10-pt dataset by a Gaussian Process (part of a larger bayesian optimization framework) and I don't understand why it is only being explained by noise.

Here is the data, on a 2x5 grid.

I am able to perfectly fit the data if only using an RBF kernel, which gets optimized to

6.54**2 * RBF(length_scale=[0.01, 0.01])

But if I add a noise WhiteKernel, then the data gets explained by the noise

4.85**2 * RBF(length_scale=[146, 6.6]) + WhiteKernel(noise_level=17.2)

and the mean prediction looks like:

I need to add noise because the real data I am working with is noisy (EMG responses to brain stimulation are not deterministic). However I noticed it was mostly explained by noise, so tried to see what the GP would learn with the mean (deterministic) function of the data, and it is again only explained by noise.

I do know about the 5th chapter of Rasmussen and Williams' Gaussian Process Textbook, and the fact that the marginal likelihood can be multimodel, since some functions can equally well be explained by noise or quickly-varying functions, but in this case only the noise version seems to happen. Here is the code for reference.

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

x = np.arange(0,2)
y = np.arange(0,5)
xx,yy = np.meshgrid(x,y)
xr,yr = xx.ravel(), yy.ravel()
xtrain = np.c_[xr,yr]

ytrain = np.array([[ 0.87611],
[ 3.58377],
[ 1.32999],
[ 7.63373],
[ 1.09371],
[ 9.18058],
[13.91528],
[ 3.21335],
[ 3.26006],
[ 7.30752]])

fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.scatter(xx, yy, ytrain.reshape((5,2)), antialiased=False)

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, WhiteKernel

kernel = 1.0 * RBF(length_scale=[1.,1.],
length_scale_bounds=[(1e-2, 1e3), (1e-2, 1e3)]) \
+ WhiteKernel(noise_level=1, noise_level_bounds=(1e-10, 1e+2))
gp = GaussianProcessRegressor(kernel=kernel,
alpha=1e-10,
n_restarts_optimizer=5).fit(xtrain, ytrain)
print(gp.kernel_)

fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.scatter(xx, yy, gp.predict(xtrain).reshape((5,2)), antialiased=False)
plt.show()


I've played around with this a little bit, and it looks like you have two local minima: one explains almost everything with noise, and the other explains almost everything by high-frequency oscillations. This is easier to see if you do a couple of things:

1. Normalize your training data to have mean 0 and variance 1. Because your GP has a zero-mean prior, shifting your data to 0 mean helps. Normalizing the variance to 1 makes the noise variance more interpretable. You can do this by setting
gp = GaussianProcessRegressor(kernel=kernel,
alpha=1e-10,
n_restarts_optimizer=5,
normalize_y=True).fit(xtrain, ytrain)

1. Change the bounds on your hyperparameters to reflect this. You'll need to set the lower bounds on the lengthscale much lower (I used 1e-10), and I'd upper bound the noise at 1 (it doesn't make sense for the noise variance to be greater than the data variance). I'd also do 100 restarts instead of 5.

Even after doing this, the hyperparameters are pretty unstable: I get different results running it multiple times, with a cluster where the noise is 0.3 - 0.6 and the lengthscales are 1e-10, and another where the length scales are more moderate but the noise is 1e-6 - 1e-9. If you want to a more complete picture of what's going on, you probably need to set proper priors over the hyperparameters and use Hamiltonian Monte Carlo to estimate the posterior. This is a pretty good guide.

• Hi Kevin. Thanks a lot for taking the time to do this! A few questions: 1) Never understood why normalizing the input should matter. According to Williams/Rasmussen (ch 2.7) it shouldn't matter. 2) normalize_y=True doesn't actually normalize y, but only centers it. so the variance here would still be large. 3) I don't understand how really small lengthscale (1e-10!) make sense since the domain is a 2x5 integer grid (so shouldn't only lengthscales of 1 make sense...?) 4) Is the prior+HMC necessary or is using bounds to force it to converge to where I want enough? – samlaf Feb 21 '19 at 18:03
• Williams/Ramussen: "It is common but by no means necessary to consider GPs with a zero mean function. Note that this is not necessarily a drastic limitation, since the mean of the posterior process is not confined to be zero." – samlaf Feb 21 '19 at 18:07
• In the limit of large data, the prior mean has no effect on the posterior mean. However, we're not quite there yet, so using a zero-mean prior to model data with a significantly non-zero mean can lead to misbehaving posteriors. Rasmussen and Williams are talking about the mean function of the GP there, not whether we should normalize the data. – Kevin Yang Feb 22 '19 at 19:05
• Is your predictive domain limited to that grid as well? – Kevin Yang Feb 22 '19 at 19:05
• Bounding the hyperparameters is roughly equivalent to using a uniform prior. The benefit of HMC is that it allows you to sample from the posterior over the hyperparameters instead of from a maximum likelihood estimate of the hyperparameters. – Kevin Yang Feb 22 '19 at 19:07