高斯过程--GPyTorch回归教程

0、前言

本文是基于官方文档的教程的

知识储备：初步了解高斯过程

1. 高斯过程很好的教程
2. 这个知乎专栏也很好，讲的挺透的

1、导入

我们将用RBF（Radial Basis Function）核的高斯过程对下面这个函数建模：
$$y=\sin (2\pi x)+ϵ,ϵ\sim \mathcal{N}(0,0.04)$$
用100个样本训练，51个样本测试。

import os
import math

from matplotlib import pyplot as plt
import torch
import gpytorch
from gpytorch import kernels, means, models, mlls, settings
from gpytorch import distributions as distr

2、生成训练数据

选取100个在$[0,1]$上等距分布的点作为 { x 1 , x 2 , … , x 100 } \{x_1,x_2,\ldots,x_{100}\} {x1​,x2​,…,x100​}。

# Training data is 100 points in [0,1] inclusive regularly spaced
train_x = torch.linspace(0, 1, 100)
# True function is sin(2*pi*x) with Gaussian noise
train_y = torch.sin(train_x * (2 * math.pi)) + torch.randn(train_x.size()) * math.sqrt(0.04)

3、模型建立

与大多数GP的package不同，GPyTorch不向用户提供完整的GP模型。GPyTorch提供的是建立一个GP模型所需的必要工具。这么做给用户留存了很大的灵活性，方便定制个性化的模型。

对大多数的GP回归，我们需要下面这些GPyTorch对象：

1. GP模型GP models，gpytorch.models.ExactGP–这负责大部分的推理
2. 似然Likelihoodgpytorch.likelihoods.GaussianLikelihood–GP回归最常用的似然
3. 均值Mean，定义GP的先验均值。如果你不知道用什么均值，可以先从gpytorch.means.ConstantMean()开始。
4. 核kernel，定义GP的先验协方差。如果你不知道用什么核，可以先从gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel())开始。
5. 多元正态分布Multivariate Normal Distributiongpytorch.distributions.MultivariateNormal。

3.1、GP模型

总的来说，用户在GPyTorch建立的GP模型包括以下组件：

• __init__方法 ，需输入训练数据和似然，同时定义modelforward方法需要的必要对象，主要包括均值mean模组与核kernel模组
• forward方法，需输入$n\times d$的数据x，返回带有在x处计算的先验均值与协方差的MultivariateNormal对象。也就是说，返回的是GP的先验均值向量$\mu (x)$与协方差矩阵$K_{xx}$

定义模型时，用户有很大的弹性空间。比如，如需通过叠加来组合两个kernels，既可以把两个模组modules直接相加：

self.covar_module = kernels.ScaleKernel(RBFKernel() + kernels.LinearKernel())

也可以在forward方法里把kernels的输出相加：

covar_x = self.rbf_kernel_module(x) + self.white_noise_module(x)

3.2、似然

用于回归最简单的似然是gpytorch.likelihoods.GaussianLikelihood，这基于噪声模型同方差homoskedastic的假设（即所有输入具有相同的观测噪声）。

此外还有别的针对exact GP regression的似然，如FixedNoiseGaussianLikelihood，可对不同的训练输入赋不同的观测噪声。

# We will use the simplest form of GP model, exact inference
class ExactGPModel(models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super(ExactGPModel, self).__init__(train_x, train_y, likelihood)
        self.mean_module = means.ConstantMean()
        self.covar_module = kernels.ScaleKernel(kernels.RBFKernel())

    def forward(self, x):
        mean_x = self.mean_module(x)
        covar_x = self.covar_module(x)
        return distr.MultivariateNormal(mean_x, covar_x)


# initialize likelihood and model
likelihood = gpytorch.likelihoods.GaussianLikelihood()
model = ExactGPModel(train_x, train_y, likelihood)

3.3、模型的模式

正如大多数PyTorch模组，ExactGP也有train()和eval()两种模式

• train()用于优化模型超参数
• eval()用于通过模型后验做预测

4、训练模型

下面将用Type-Ⅱ MLE训练模型。

这里与许多别的GP实现最显著的差异在于，核心训练代码是要用户自己写的，正如PyTorch一样。在GPyTorch中，我们可以调用来自torch.optim的PyTorch优化器，所有可学习的参数都应该属于torch.nn.Parameter的类型。由于GP模型直接从torch.nn.Modules衍生而来，我们可以像在PyTorch里面一样调用model.parameters()或model.named_parameters()。

大多数情况下，下面的套路效果还不错。它跟PyTorch中的模型训练循环差不多。

1. 所有参数梯度置零
2. 调用模型，计算损失loss
3. loss反向传播填充梯度
4. 优化器优化

考虑到GPyTorch较大的弹性空间，我们也可以定义个性化的训练loops。比如，可以很容易地保留每一步训练得到的参数，抑或针对不同参数使用不同的学习率（在deep kernel learning中挺有用）。

smoke_test = ('CI' in os.environ)
training_iter = 2 if smoke_test else 50

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

# Use the adam optimizer
optimizer = optim.Adam(model.parameters(), lr=0.1)  # Includes GaussianLikelihood parameters

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

for i in range(training_iter):
    # Zero gradients from previous iteration
    optimizer.zero_grad()
    # Output from model
    output = model(train_x)
    # Calc loss and backprop gradients
    loss = -mll(output, train_y)
    loss.backward()
    if i % 5 == 4:
        print('Iter %d/%d - Loss: %.3f   lengthscale: %.3f   noise: %.3f' % (
            i + 1, training_iter, loss.item(),
            model.covar_module.base_kernel.lengthscale.item(),
            model.likelihood.noise.item()
        ))
    optimizer.step()

前两行smoke_test部分我没看懂要干啥。。。期待大佬指点！

Iter 5/50 - Loss: 0.811   lengthscale: 0.514   noise: 0.513
Iter 10/50 - Loss: 0.549   lengthscale: 0.342   noise: 0.339
Iter 15/50 - Loss: 0.360   lengthscale: 0.240   noise: 0.215
Iter 20/50 - Loss: 0.194   lengthscale: 0.205   noise: 0.133
Iter 25/50 - Loss: 0.047   lengthscale: 0.208   noise: 0.082
Iter 30/50 - Loss: -0.053   lengthscale: 0.234   noise: 0.052
Iter 35/50 - Loss: -0.081   lengthscale: 0.265   noise: 0.036
Iter 40/50 - Loss: -0.064   lengthscale: 0.279   noise: 0.029
Iter 45/50 - Loss: -0.061   lengthscale: 0.265   noise: 0.028
Iter 50/50 - Loss: -0.075   lengthscale: 0.239   noise: 0.031

5、用模型做预测

为了用训练好的模型做预测，我们只需把model和likelihood都设成eval()模式，然后在测试数据上调用它们。

正如用户定义的GP模型从forward方法返回的是包含先验统计量的MultivariateNormal，训练完成的GP模型在eval()模式下则返回包含后验统计量的MultivariateNormal。因此，获取预测均值与方差，然后从给定的测试点采样函数可以通过下面的代码来实现：

f_preds = model(test_x)
y_preds = likelihood(model(test_x))

f_mean = f_preds.mean
f_var = f_preds.variance
f_covar = f_preds.covariance_matrix
f_samples = f_preds.sample(sample_shape=torch.Size((1000,)))

gpytorch.settings.fast_pred_var的部分可以不要，但这儿我们还是给出一个例子，用LOVE来做快速的分布预测。

# Get into evaluation (predictive posterior) mode
model.eval()
likelihood.eval()

# Test points are regularly spaced along [0,1]
# Make predictions by feeding model through likelihood
with torch.no_grad(), settings.fast_pred_var():
    test_x = torch.linspace(0, 1, 51)
    observed_pred = likelihood(model(test_x))

6、拟合模型可视化

下面我们要画出GP模型的均值与置信区域。confidence_region方法可以返回均值之上和之下的两个标准偏差。

with torch.no_grad():
    # Initialize plot
    f, ax = plt.subplots(1, 1, figsize=(4, 3))

    # Get upper and lower confidence bounds
    lower, upper = observed_pred.confidence_region()
    # Plot training data as black stars
    ax.plot(train_x.numpy(), train_y.numpy(), 'k*')
    # Plot predictive means as blue line
    ax.plot(test_x.numpy(), observed_pred.mean.numpy(), 'b')
    # Shade between the lower and upper confidence bounds
    ax.fill_between(test_x.numpy(), lower.numpy(), upper.numpy(), alpha=0.5)
    ax.set_ylim([-3, 3])
    ax.legend(['Observed Data', 'Mean', 'Confidence'])
plt.show()