为什么gpy手电筒似乎不像科学知识那样准确？

Question

我目前发现了gpy手电筒(https://github.com/cornellius-gp/gpytorch)。它似乎是将探地雷达集成到火把中的一个很好的包。第一次测试也呈阳性。使用gpy手电筒、GPU-Power以及智能算法，可以与其他软件包(如scikit-learn )相比，提高性能。

然而，我发现估计所需的超参数要困难得多。在科学工具包中--学习在后台发生的事情，并且非常健壮。我想从社区获得一些关于原因的信息，并讨论是否有更好的方法来估计这些参数，而不是gpy手电筒文档中的示例所提供的。

为了便于比较，我在gpy手电筒(Regression.ipynb)的正式页面上使用了一个提供的示例代码，并对其进行了两部分的修改：

1. 我使用不同的内核(gpytorch.kernels.MaternKernel(nu=2.5)代替gpytorch.kernels.RBFKernel())
2. 我使用了一个不同的输出函数

在下面的文章中，我首先提供了使用gpy手电筒的代码。随后，我提供了scikit的代码--学习。最后，对实验结果进行了比较。

进口(用于gpy手电和scikit-学习)：

import math
import torch
import numpy as np
import gpytorch

生成数据(用于gpy手电和scikit-学习)：

n = 20
train_x = torch.zeros(pow(n, 2), 2)
for i in range(n):
    for j in range(n):
        # Each coordinate varies from 0 to 1 in n=100 steps
        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))/4
train_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))/4
test_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) # For a more robust comparison
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)

# Find optimal model hyperparameters
model.train()
likelihood.train()

# Use the adam optimizer
optimizer = torch.optim.Adam([
    {'params': model.parameters()},  # Includes GaussianLikelihood parameters
], lr=0.1)

# "Loss" for GPs - the marginal log likelihood
mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)

n_iter = 50
for 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()

# Set into eval mode
model.eval()
likelihood.eval()

# Make predictions
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的代码--学习。比较方便的^^：

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, Matern
kernel = 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_scitlearn
print('Variable 1: scitkit learn is more accurate my factor: ' + str(abs(comparisson[0]))
print('Variable 2: scitkit learn is more accurate my factor: ' + str(comparisson[1]))

不幸的是，我并没有找到一个简单的方法来修复科学种子--学习。上次运行代码时，它返回：

变量1: scitkit学习更准确--我的因素: 11.362540360431087 变量2: scitkit学习更准确--我的因素: 29.64760087022618

在gpy手电筒的情况下，我假设优化器运行在一些局部的optima中。但是，我想不出还有任何更健壮的优化算法仍然使用pytorch。

我期待着你的建议！

拉兹洛

Answer 1

(我还回答了您关于GitHub问题的问题，您为它创建了这里)

这主要是因为您在sklearn和gpy手电中使用了不同的模型。特别是，在默认情况下，sklearn在多输出设置中学习独立的GPs (例如，讨论这里)。在GPyTorch中，您使用了Bonilla等人，2008年中引入的多任务GP方法。纠正这一差额会产生以下结果：

test_results_gpytorch = 5.207913e-04 -8.469360e-05 test_results_scitlearn = 3.65288816e-04 4.79017145e-05