为 GPR 创建自定义内核

Question

我想编写一个仅在 X 轴特定范围内工作的 RBF 内核。我尝试编写一个包含 RBF 内核的类来测试代码

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()

实现和适配是这样的

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)

如果我运行代码，我会得到以下输出

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

如果我尝试使用通常的 RBF 内核，代码可以正常工作。我还尝试禁用优化器

optimizer=None

并且代码可以工作，但我得到了非常大的错误。

Answer 1

似乎添加一个属性装饰器

hyperparameter_length_scale

（如RBF的源代码）可以解决问题：

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()

使用来自 scikit-learn 的 GaussianProcessRegressor 示例的小型合成数据进行测试：

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

def train_pred_plot(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")

使用 RBF 内核

kernel = 1.0 * RBF(length_scale=1.0, length_scale_bounds=(1e-2, 1e2))
train_pred_plot(X_train, y_train, 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_train, y_train, kernel)

我们得到以下预测置信区间：