Create custom kernel for GPR

Question

I would like to write a RBF kernel that is working only in a specific range on X axis. I tried to write a class that contains a RBF kernel to test the code

class RangeLimitedRBFTest(Kernel):
    def __init__(self, length_scale=1.0, length_scale_bounds=(1e-5, 1e5), x_min = 0., x_max = 1.):
        self.length_scale = length_scale
        self.length_scale_bounds = length_scale_bounds
        self.rbf_kernel = RBF(length_scale, length_scale_bounds)
        self.x_min = x_min
        self.x_max = x_max

    def __call__(self, X, Y=None, eval_gradient=False):
        if eval_gradient and Y is not None:
            raise ValueError("Gradient can only be evaluated when Y is None.")
        
        X = np.atleast_2d(X)
        if Y is not None:
            Y = np.atleast_2d(Y)

        print(f"X shape: {X.shape}")
        if Y is not None:
            print(f"Y shape: {Y.shape}")
        else:
            print("Y shape: None")

        K_rbf = self.rbf_kernel(X, Y, eval_gradient=eval_gradient)

        if eval_gradient:
            K, K_grad = K_rbf
            print(f"Kernel matrix shape (K): {K.shape}")
            print(f"Kernel gradient matrix shape (K_grad): {K_grad.shape}")
            return K, K_grad
        else:
            K = K_rbf
            return K

    def diag(self, X):
        return self.rbf_kernel.diag(X)

    def is_stationary(self):
        return self.rbf_kernel.is_stationary()

The implementation and the fitting is like this

kernel = 1.0 * RangeLimitedRBFTest(length_scale=0.1, length_scale_bounds=(8e-2, 8e-1), x_min=0., x_max=2.5) + WhiteKernel(noise_level=0.5, noise_level_bounds=(1e-2, 1e1))
gaussian_process = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=1, alpha=1e-5, optimizer='fmin_l_bfgs_b')
gaussian_process.optimizer_kwargs = {"max_iter": 10000} 
gaussian_process.fit(X, T_PMT)

If I run the code I get the following output

X shape: (6248, 1)
Y shape: None
Kernel matrix shape (K): (6248, 6248)
Kernel gradient matrix shape (K_grad): (6248, 6248, 1)
ValueError: 0-th dimension must be fixed to 2 but got 3


The above exception was the direct cause of the following exception:

Traceback (most recent call last):
  File "/home/tdaq/cremonini/pt100_probe/read_temperatures.py", line 97, in <module>
    gaussian_process.fit(X, T_PMT)
  File "/home/tdaq/.local/lib/python3.10/site-packages/sklearn/base.py", line 1389, in wrapper
    return fit_method(estimator, *args, **kwargs)
  File "/home/tdaq/.local/lib/python3.10/site-packages/sklearn/gaussian_process/_gpr.py", line 308, in fit
    self._constrained_optimization(
  File "/home/tdaq/.local/lib/python3.10/site-packages/sklearn/gaussian_process/_gpr.py", line 653, in _constrained_optimization
    opt_res = scipy.optimize.minimize(
  File "/cvmfs/atlas.cern.ch/repo/sw/software/0.3/StatAnalysis/0.3.1/InstallArea/x86_64-el9-gcc13-opt/lib/python3.10/site-packages/scipy/optimize/_minimize.py", line 713, in minimize
    res = _minimize_lbfgsb(fun, x0, args, jac, bounds,
  File "/cvmfs/atlas.cern.ch/repo/sw/software/0.3/StatAnalysis/0.3.1/InstallArea/x86_64-el9-gcc13-opt/lib/python3.10/site-packages/scipy/optimize/_lbfgsb_py.py", line 360, in _minimize_lbfgsb
    _lbfgsb.setulb(m, x, low_bnd, upper_bnd, nbd, f, g, factr,
ValueError: failed in converting 7th argument `g' of _lbfgsb.setulb to C/Fortran array

If i try to use the usual RBF kernel, the code works without any problem. I tried also to disable the optimizer optimizer=None and the code works but I get very large error.

Answer 1

Seems that adding a property decorator hyperparameter_length_scale as in the source code for RBF, resolves the issue (otherwise there is a dimension mismatch in the kernel hyperparameters gaussian_process.kernel.theta ):

class RangeLimitedRBFTest(Kernel):

    @property
    def hyperparameter_length_scale(self):
        return Hyperparameter("length_scale", "numeric", self.length_scale_bounds)
    
    def __init__(self, length_scale=1.0, length_scale_bounds=(1e-5, 1e5), x_min = 0., x_max = 1.):
        self.length_scale = length_scale
        self.length_scale_bounds = length_scale_bounds
        self.rbf_kernel = RBF(length_scale, length_scale_bounds)
        self.x_min, self.x_max = x_min, x_max
    
    def __call__(self, X, Y=None, eval_gradient=False):
        if eval_gradient and Y is not None:
            raise ValueError("Gradient can only be evaluated when Y is None.")
        X = np.atleast_2d(X)
        if Y is not None:
            Y = np.atleast_2d(Y)
        return self.rbf_kernel(X, Y, eval_gradient=eval_gradient) 
        
    def diag(self, X):
         return self.rbf_kernel.diag(X)

    def is_stationary(self):
        return self.rbf_kernel.is_stationary()

Testing with small synthetic data from scikit-learn's GaussianProcessRegressor example:

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

def train_pred_plot(X, y, X_train, y_train, kernel):
    gaussian_process = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=9)
    gaussian_process.fit(X_train, y_train)        
    mean_prediction, std_prediction = gaussian_process.predict(X, return_std=True)        
    plt.plot(X, y, label=r"$f(x) = x \sin(x)$", linestyle="dotted")
    plt.scatter(X_train, y_train, label="Observations")
    plt.plot(X, mean_prediction, label="Mean prediction")
    plt.fill_between(
        X.ravel(),
        mean_prediction - 1.96 * std_prediction,
        mean_prediction + 1.96 * std_prediction,
        alpha=0.5,
        label=r"95% confidence interval",
    )
    plt.legend()
    plt.xlabel("$x$")
    plt.ylabel("$f(x)$")
    plt.title("Gaussian process regression on noise-free dataset")

with RBF kernel

kernel = 1.0 * RBF(length_scale=1.0, length_scale_bounds=(1e-2, 1e2))
train_pred_plot(X, y, X_train, y_train, kernel)

we get the following prediction confidence intervals:

whereas with the custom kernel

kernel = 1.0 * RangeLimitedRBFTest(length_scale=0.1, length_scale_bounds=(8e-2, 8e-1), x_min=0., x_max=2.5) + WhiteKernel(noise_level=0.5, noise_level_bounds=(1e-2, 1e1))
train_pred_plot(X, y, X_train, y_train, kernel)

we get the following prediction confidence intervals: