如何在scikit-learn的GaussianProcessRegressor中使用预计算核？

问题原因 #

1. GaussianProcessRegressor（以下简称GPR）的API设计和SVM系列不同，不支持直接给kernel参数传入'precomputed'字符串。GPR要求kernel参数必须是sklearn.gaussian_process.kernels下的核类实例，传入字符串后sklearn克隆参数时找不到get_params方法，就会触发你遇到的TypeError。
2. 你测试用的随机核矩阵K_train = np.random.randn(80, 80)不是半正定矩阵，高斯过程的核矩阵必须满足半正定性要求，这也是后续拟合异常的原因，也是你换用SVR后出现无限运行的核心原因：非半正定核会导致SVR的拟合迭代收敛速度极慢，看起来就像卡住。

解决方法 #

我们可以自定义一个适配GPR接口的预计算核类，实现必要的__call__和get_params方法，同时生成半正定的核矩阵用于测试，修正后的可运行代码如下：

from sklearn.model_selection import KFold
from sklearn.utils import check_random_state
from sklearn.gaussian_process import GaussianProcessRegressor as GPR
from sklearn.gaussian_process.kernels import Kernel
from sklearn.metrics import mean_squared_error
import numpy as np

# 适配GPR接口的预计算核类
class PrecomputedKernel(Kernel):
    def __init__(self, K_train, K_test=None):
        self.K_train = K_train
        self.K_test = K_test
    
    def __call__(self, X, Y=None, eval_gradient=False):
        # 训练阶段Y为None，返回训练核矩阵
        if Y is None:
            return self.K_train
        # 预测阶段返回测试核矩阵
        return self.K_test
    
    def get_params(self, deep=True):
        return {"K_train": self.K_train, "K_test": self.K_test}

def Kfold_CV_GPR(Gs, y, n_iter=4, n_splits=5, random_state=None):
    random_state = check_random_state(random_state)
    kf = KFold(n_splits=n_splits, random_state=random_state, shuffle=True)
    errors = []
    for train_idxs, test_idxs in kf.split(y):
        # 替换为实际grakel核计算逻辑即可
        # gk = WeisfeilerLehman(n_iter=n_iter, normalize=True)
        # K_train = gk.fit_transform(Gs[train_idxs])
        # K_test = gk.transform(Gs[test_idxs])
        
        # 生成测试用的半正定随机核矩阵
        rand_mat = np.random.randn(80, 80)
        K_train = rand_mat @ rand_mat.T + 1e-6 * np.eye(80) # 加小常数避免数值奇异
        K_test = np.random.randn(20, 80)

        kernel = PrecomputedKernel(K_train, K_test)
        # 关闭核参数优化，预计算核不需要调整参数
        gpr = GPR(kernel=kernel, optimizer=None)
        # 传入任意形状匹配的占位矩阵即可
        gpr.fit(np.zeros((len(train_idxs), 1)), y[train_idxs])
        y_pred = gpr.predict(np.zeros((len(test_idxs), 1)))

        rmse = mean_squared_error(y[test_idxs], y_pred, squared=False)
        errors.append(rmse)
    return -np.mean(errors)

score = Kfold_CV_GPR(Gs=None, y=np.random.randn(100, ), n_iter=4, n_splits=5)
print(score)

注意事项 #

• 实际使用时把测试核生成的代码替换为你注释的grakel核计算逻辑即可，grakel输出的核矩阵本身是半正定的，不需要额外处理。
• 必须设置optimizer=None，否则GPR会尝试优化核参数，既没有意义也会额外消耗计算资源。

内容的提问来源于stack exchange，提问作者Shaun Han