Bug in scikit-learn GaussianProcessRegressor

[dreamer2368]

GaussianProcessRegressor exhibits an erroneous behavior depending on the data. Specifically, optimal length scale for kernel shrinks down to minimum value when fitted with the data that has a sign change. Following is the minimum failing example.

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

from sklearn.gaussian_process.kernels import ConstantKernel, Matern, RBF
from sklearn.gaussian_process import GaussianProcessRegressor

x = np.array([[0.], [1.0]])
xpred = np.linspace(x[0][0], x[-1][0], 1000)
xpred = xpred.reshape(-1, 1)

# this is a working dataset, which does not have a sign change.
y1 = np.array([-4, -1e-3])
y1 = y1.reshape(-1, 1)

kernel = ConstantKernel() * RBF(length_scale_bounds = (1e-5, 1e5))
gp = GaussianProcessRegressor(kernel = kernel, n_restarts_optimizer = 100, random_state = 1, alpha=1e-10)
gp = gp.fit(x, y1)

print(gp.kernel_.get_params()['k2__length_scale'])
yavg1, ystd1 = gp.predict(xpred, return_std = True)

plt.figure(1)
plt.plot(x, y1, 'ok')
plt.plot(xpred, yavg1, '-r')
plt.plot(xpred, yavg1 + 2 * ystd1, '--r')
plt.plot(xpred, yavg1 - 2 * ystd1, '--r')
plt.title('With (0, -4) and (1, -1e-3)')

# this is a failing dataset, which is close to the dataset above but has a sign change.
y2 = np.array([-4, 1e-3])
y2 = y2.reshape(-1, 1)

kernel = ConstantKernel() * RBF(length_scale_bounds = (1e-5, 1e5))
gp = GaussianProcessRegressor(kernel = kernel, n_restarts_optimizer = 100, random_state = 1, alpha=1e-10)
gp = gp.fit(x, y2)

print(gp.kernel_.get_params()['k2__length_scale'])
yavg2, ystd2 = gp.predict(xpred, return_std = True)

plt.figure(2)
plt.plot(x, y2, 'ok')
plt.plot(xpred, yavg2, '-r')
plt.plot(xpred, yavg2 + 2 * ystd2, '--r')
plt.plot(xpred, yavg2 - 2 * ystd2, '--r')
plt.title('With (0, -4) and (1, +1e-3)')

The reason for this behavior is not known yet. Running the same example with other packages, such as GPy, does not show this behavior, so this is a bug from scikit-learn package.

So far, the only so-called 'fix' is to set the lower bound to a 'reasonable' value. However, this does not really fix the root cause and forfeits the reason of using GP for properly tuning length scale hyperparameter under statistical principle. At this point, scikit-learn is better to be replaced with other GP packages.

[CBonneville45]

This is actually not a bug and isn't due to the sign of the data (in general). It is an unfortunate but quite common instability inherent to Gaussian processes. In GP inference and training, you need at some point to invert a kernel matrix (this is all done under the hood in sklearn and most common GP packages), but this matrix can be poorly conditioned depending on your data. To make the GP stable, a common way around is to add jitter on the kernel matrix diagonal. "Physically" this jitter is equivalent to the inherent noise in the data (so for virtually noiseless data you should be able to get away with adding any jitter). In sklearn you can manually tune the jitter with the alpha parameter. With a larger alpha value you shouldn't observe the instability anymore (but with a value too large, the GP is essentially gonna assume that the data is all noise and underfit)

…

________________________________ De : dreamer2368 ***@***.***> Envoyé : Saturday, December 21, 2024 7:15:21 AM À : LLNL/GPLaSDI ***@***.***> Cc : Subscribed ***@***.***> Objet : [LLNL/GPLaSDI] Bug in scikit-learn GaussianProcessRegressor (Issue #22) GaussianProcessRegressor exhibits an erroneous behavior depending on the data. Specifically, optimal length scale for kernel shrinks down to minimum value when fitted with the data that has a sign change. Following is the minimum failing example. import numpy as np import matplotlib.pyplot as plt %matplotlib inline from sklearn.gaussian_process.kernels import ConstantKernel, Matern, RBF from sklearn.gaussian_process import GaussianProcessRegressor x = np.array([[0.], [1.0]]) xpred = np.linspace(x[0][0], x[-1][0], 1000) xpred = xpred.reshape(-1, 1) y1 = np.array([-4, -1e-3]) y1 = y1.reshape(-1, 1) kernel = ConstantKernel() * RBF(length_scale_bounds = (1e-5, 1e5)) gp = GaussianProcessRegressor(kernel = kernel, n_restarts_optimizer = 100, random_state = 1, alpha=1e-10) gp = gp.fit(x, y1) print(gp.kernel_.get_params()['k2__length_scale']) yavg1, ystd1 = gp.predict(xpred, return_std = True) plt.figure(1) plt.plot(x, y1, 'ok') plt.plot(xpred, yavg1, '-r') plt.plot(xpred, yavg1 + 2 * ystd1, '--r') plt.plot(xpred, yavg1 - 2 * ystd1, '--r') plt.title('With (0, -4) and (1, -1e-3)') y2 = np.array([-4, 1e-3]) y2 = y2.reshape(-1, 1) kernel = ConstantKernel() * RBF(length_scale_bounds = (1e-5, 1e5)) gp = GaussianProcessRegressor(kernel = kernel, n_restarts_optimizer = 100, random_state = 1, alpha=1e-10) gp = gp.fit(x, y2) print(gp.kernel_.get_params()['k2__length_scale']) yavg2, ystd2 = gp.predict(xpred, return_std = True) plt.figure(2) plt.plot(x, y2, 'ok') plt.plot(xpred, yavg2, '-r') plt.plot(xpred, yavg2 + 2 * ystd2, '--r') plt.plot(xpred, yavg2 - 2 * ystd2, '--r') plt.title('With (0, -4) and (1, +1e-3)') The reason for this behavior is not known yet. Running the same example with other packages, such as GPy, does not show this behavior, so this is a bug from scikit-learn package. So far, the only so-called 'fix' is to set the lower bound to a 'reasonable' value. However, this does not really fix the root cause and forfeits the reason of using GP for properly tuning length scale hyperparameter under statistical principle. At this point, scikit-learn is better to be replaced with other GP packages. — Reply to this email directly, view it on GitHub<#22>, or unsubscribe<https://github.com/notifications/unsubscribe-auth/ALEYFFYRZGEMYHLVLRHXPCD2GWAYTAVCNFSM6AAAAABUAWUEWWVHI2DSMVQWIX3LMV43ASLTON2WKOZSG42TIMJTGUZTGMY>. You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[dreamer2368]

I've been looking into this deeper, and I agree that this is not a bug. But this is not simply a matter of the common instability. You'll see that it cannot be solved by simply increasing the measurement noise (the alpha value). Rather, it is an inherent issue of maximum likelihood estimator (MLE) with very few data.

Below is the calibrated GP with the second dataset, where the correlation scale shrank down to 0 and the function is essentially independent Gaussian noise:

This is the log-likelihood depending on the correlation scale hyperparameter, which monotonically decreases:

For the first dataset, the correlation scale looks "reasonable":

The log-likelihood looks almost similar to the first one, but luckily it has a very slight bump in the middle:

The main takeaway of these figures is that GP chooses the best parameters given the data and its likelihood. There is nothing wrong with it. Rather, the problem is that there is too few data provided. MLE calibrates hyperparameters solely from data, so if the data size is too small then it leads to these problematic values. This is actually a common situation in the beginning of the greedy sampling.

If we want to handle this more stably, we'll probably need a more formal Bayesian calibration, where the likelihood is compensated with prior knowledge (AKA maximum a posteriori (MAP) estimation). When the data size is small, the prior dominates over the likelihood. As the data increases, the likelihood dominates over the prior (and at that point we do not have this issue any more). A caveat is that scikit-learn may not be sufficient for a custom estimation. I'm looking around to see how I may implement this.