我最近发现了gpytorch (https://github.com/cornellius-gp/gpytorch)。它似乎是一个将GPR集成到pytorch中的优秀工具包。初步测试也非常积极。使用gpytorch，可以利用GPU的强大性能以及智能算法，相较于其他工具包如scikit-learn，性能得到了提升。
然而，我发现估计所需的超参数要困难得多。在scikit-learn中，这部分工作在后台进行，并且非常稳健。我希望从社区获得一些关于原因的反馈，并讨论是否有比gpytorch文档中示例提供的更好的方法来估计这些参数。
为了比较，我使用了gpytorch官方页面上提供的一个示例代码(https://github.com/cornellius-gp/gpytorch/blob/master/examples/03_Multitask_GP_Regression/Multitask_GP_Regression.ipynb)，并在两个部分进行了修改:
- 我使用了不同的核函数（gpytorch.kernels.MaternKernel(nu=2.5) 代替 gpytorch.kernels.RBFKernel()）
- 我使用了不同的输出函数
下面，我首先提供使用gpytorch的代码。随后，我提供scikit-learn的代码。最后，我比较结果
导入（用于gpytorch和scikit-learn）:
import mathimport torchimport numpy as npimport gpytorch生成数据（用于gpytorch和scikit-learn）:
n = 20train_x = torch.zeros(pow(n, 2), 2)for i in range(n): for j in range(n): # 每个坐标从0到1变化，n=100步 train_x[i * n + j][0] = float(i) / (n-1) train_x[i * n + j][1] = float(j) / (n-1)train_y_1 = (torch.sin(train_x[:, 0]) + torch.cos(train_x[:, 1]) * (2 * math.pi) + torch.randn_like(train_x[:, 0]).mul(0.01))/4train_y_2 = torch.sin(train_x[:, 0]) + torch.cos(train_x[:, 1]) * (2 * math.pi) + torch.randn_like(train_x[:, 0]).mul(0.01)train_y = torch.stack([train_y_1, train_y_2], -1)test_x = torch.rand((n, len(train_x.shape)))test_y_1 = (torch.sin(test_x[:, 0]) + torch.cos(test_x[:, 1]) * (2 * math.pi) + torch.randn_like(test_x[:, 0]).mul(0.01))/4test_y_2 = torch.sin(test_x[:, 0]) + torch.cos(test_x[:, 1]) * (2 * math.pi) + torch.randn_like(test_x[:, 0]).mul(0.01)test_y = torch.stack([test_y_1, test_y_2], -1)现在来进行如引用的文档中所述的示例中的估计:
torch.manual_seed(2) # 为了更稳健的比较class MultitaskGPModel(gpytorch.models.ExactGP): def __init__(self, train_x, train_y, likelihood): super(MultitaskGPModel, self).__init__(train_x, train_y, likelihood) self.mean_module = gpytorch.means.MultitaskMean( gpytorch.means.ConstantMean(), num_tasks=2 ) self.covar_module = gpytorch.kernels.MultitaskKernel( gpytorch.kernels.MaternKernel(nu=2.5), num_tasks=2, rank=1 ) def forward(self, x): mean_x = self.mean_module(x) covar_x = self.covar_module(x) return gpytorch.distributions.MultitaskMultivariateNormal(mean_x, covar_x)likelihood = gpytorch.likelihoods.MultitaskGaussianLikelihood(num_tasks=2)model = MultitaskGPModel(train_x, train_y, likelihood)# 寻找最优模型超参数model.train()likelihood.train()# 使用adam优化器optimizer = torch.optim.Adam([ {'params': model.parameters()}, # 包括GaussianLikelihood参数], lr=0.1)# GPs的"损失" - 边缘对数似然mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)n_iter = 50for i in range(n_iter): optimizer.zero_grad() output = model(train_x) loss = -mll(output, train_y) loss.backward() # print('Iter %d/%d - Loss: %.3f' % (i + 1, n_iter, loss.item())) optimizer.step()# 设置为评估模式model.eval()likelihood.eval()# 进行预测with torch.no_grad(), gpytorch.settings.fast_pred_var(): predictions = likelihood(model(test_x)) mean = predictions.mean lower, upper = predictions.confidence_region()test_results_gpytorch = np.median((test_y - mean) / test_y, axis=0)下面，我提供scikit-learn的代码。这稍微方便一些^^:
from sklearn.gaussian_process import GaussianProcessRegressorfrom sklearn.gaussian_process.kernels import WhiteKernel, Maternkernel = 1.0 * Matern(length_scale=0.1, length_scale_bounds=(1e-5, 1e5), nu=2.5) \ + WhiteKernel()gp = GaussianProcessRegressor(kernel=kernel, alpha=0.0).fit(train_x.numpy(), train_y.numpy())# x_interpolation = test_x.detach().numpy()[np.newaxis, :].transpose()y_mean_interpol, y_std_norm = gp.predict(test_x.numpy(), return_std=True)test_results_scitlearn = np.median((test_y.numpy() - y_mean_interpol) / test_y.numpy(), axis=0)最后我比较结果:
comparisson = (test_results_scitlearn - test_results_gpytorch)/test_results_scitlearnprint('变量1: scikit-learn更准确的倍数: ' + str(abs(comparisson[0])))print('变量2: scikit-learn更准确的倍数: ' + str(comparisson[1]))遗憾的是，我没有找到一种简单的方法来为scikit-learn固定种子。最后一次运行代码时，它返回:
变量1: scikit-learn更准确的倍数: 11.362540360431087
变量2: scikit-learn更准确的倍数: 29.64760087022618
对于gpytorch，我认为优化器可能运行在某些局部最优点。但我无法想到任何更稳健的优化算法，同时仍然使用pytorch。
期待大家的建议！
Lazloo
回答：
(我也在你为此创建的GitHub问题中回答了你的问题 这里)
主要原因是你使用了scikit-learn和gpytorch中的不同模型。特别是，scikit-learn在多输出设置中默认学习独立的GP（例如，讨论见这里）。在GPyTorch中，你使用了Bonilla等人2008年引入的多任务GP方法（见Bonilla et al, 2008）。纠正这一差异后得到:
test_results_gpytorch = [5.207913e-04 -8.469360e-05]
test_results_scitlearn = [3.65288816e-04 4.79017145e-05]
