现代高斯过程回归速览 | 西山晴雪的知识笔记

1 无限模型表达 + 现代计算

有没有想过如何创建具有 无限表达能力 的非参数监督学习模型？看看 **高斯过程回归 (GPR)**，这是一种几乎完全根据数据本身学习做出预测的算法（在超参数的帮助下）。将此算法与自动微分等最新的计算进展相结合，可以应用高斯过程回归近乎实时地解决各种受监督的机器学习问题。

在本文中，我们将讨论：

高斯过程回归理论的简要概述/回顾
我们可以使用高斯过程回归解决的问题类型，以及一些例子
高斯过程回归与其他监督学习算法的比较
可以用来实现高斯过程回归的现代编程包和工具

这是我的高斯过程回归系列中的第二篇文章。如需从头开始对高斯过程回归进行严格的介绍，请查看我之前的文章此处。

2 高斯过程回归的概念

在深入研究如何实现和使用高斯过程回归之前，先快速回顾一下这个监督机器学习算法背后的机制和理论。关于以下概念的详细推导/讨论，请查看我之前的文章《高斯过程回归》的文章。

（i）以 观测到的 训练点为条件，预测测试点的 条件后验分布：

Fig

（ii）将测试点目标的均值被预测为 已观测到的目标值的线性组合，这些线性组合的权重，则由从训练输入到测试点的核距离确定：

fig

（iii）使用 协方差函数 来测量输入之间的 核距离：

Fig

（iv）通过将每个新点视为高斯过程的一部分，从现有点中插值新点，即将新点也参数化为一个高斯分布：

Fig

使用含噪声的正弦时间序列数据集进行一维插值的示例。

3 高斯过程回归能解决的问题

高斯过程回归可以应用于各种 有监督 机器学习问题（在某些情况下，也可以用作 无监督 机器学习中的子程序）。以下是可以使用这种机器学习技术解决的几类问题：

（A）插值/克里金法

插值是信号处理、空间统计和控制等多个领域的关键任务。此应用在利用空间统计的领域中尤为常见，例如 地统计学。作为一个具体示例，请考虑生成与下面的山坡相对应的表面问题，只给定山坡上有限数量的定义点。如果您有兴趣查看此的具体实现，请查看我的文章此处。

克里金法和插值在地统计学中经常用到，可以用于高维空间的曲面插值！ Markos Mant 在 Unsplash 上拍摄的照片

（B）时间序列预测

这类问题着眼于使用历史数据将时间序列投射到未来。与克里金法一样，时间序列预测允许预测看不见的值。然而，这个问题不是预测不同位置的看不见的值，而是应用高斯过程回归来预测未来看不见的点的均值和方差。这非常适用于预测电力需求、股票价格或线性动力系统的状态空间演化等任务。

此外，高斯过程回归不仅预测未来点的均值，而且还输出预测方差，使决策系统能够将不确定性纳入决策。

使用含噪声的正弦曲线进行时间序列预测的示例。深蓝色线代表预测均值，而浅蓝色区间代表模型的置信区间。

（C）预测不确定性

更一般地说，因为高斯过程回归允许预测测试点的方差，高斯过程回归可用于各种不确定性量化任务——即任何与估计期望值和与此 预期值 相关的不确定性或方差。

您可能想知道：_为什么不确定性很重要_？为了激发这个答案，考虑为自主导航安全系统预测行人的轨迹。如果行人的预测轨迹具有高预测不确定性，则自动驾驶车辆应更加谨慎，以应对对行人意图的信心不足。另一方面，如果自动驾驶汽车对行人的轨迹具有低预测方差，那么自动驾驶汽车将能够更好地预测行人的意图，并且可以更轻松地按照当前的驾驶计划进行。

从某种意义上说，通过预测不确定性，决策系统可以根据他们预测这些预期值的 不确定性 来 “加权” 他们估计的 预期值。

4 为什么高斯过程回归优于其他监督学习模型？

你可能想知道：为什么我应该考虑使用高斯过程回归而不是不同的监督学习模型？下面，我列举几个对比原因。

高斯过程回归是 非参数模型。这意味着它主要从数据本身学习，而不是通过学习大量参数。这是特别有利的，因为这导致高斯过程回归模型不像神经网络等高参数模型那样 **数据饥渴**，即它们不需要那么多的样本来实现强大的通用性。
对于插值和预测任务，高斯过程回归估计 预期值 和 不确定性。这对于在决策时考虑到这种不确定性的决策系统尤其有益。

高斯过程回归是一个 线性平滑器 ^[5]，从 **监督学习** 角度来看，这可以概念化为一种 **正则化** 技术。从 **贝叶斯** 的角度来看，这相当于在您的模型上强加一个 **先验**，即测试点上的所有目标必须是现有训练目标的 **_线性组合_**。这个性质有助于高斯过程回归泛化到未观测到的数据，只要真正为观测到的目标可以表示为训练目标的线性组合。 4. 通过集成 `torch` 和 `tensorflow` 等 **自动微分** 后端框架，`gpytorch` 和 `gpflow` 等高斯过程回归的软件包正在变得 **更快** 和 **根据可扩展性**。对于 **批次** 模型尤其如此。有关这方面的示例案例研究，请参阅我之前关于小批量、多维高斯过程回归的文章[**此处**](https://towardsdatascience.com/batched-multi-dimensional-gaussian-process-regression-with-gpytorch-3a6425185109)！ ## 5 如何实现高斯过程回归？下面，我们介绍几个 Python 机器学习包，用于 **可扩展的**、**高效的**、**模块化的** 高斯过程回归实现。 ### **Scikit-Learn** [Scikit-Learn](https://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.GaussianProcessRegressor.html) ^[1] 是一个很好的 GPR 入门包。它允许一定的模型灵活性，并且能够在引擎盖下进行超参数优化和定义可能性。要将 `sklearn` 与您的数据集一起使用，请确保您的数据集可以用 `np.array` 对象以数字表示。使用`sklearn`进行高斯过程回归的主要步骤： 1. 预处理您的数据。训练数据 (`np.array`) 可以表示为具有 `x_train` 形状 **(N, D)** 和 `y_train` 形状 **(N, 1)** 的 `(x_train, y_train)` 元组，其中 **N** 是样本数，**D** 是特征的维度。测试点 (`np.array`) 可以表示为形状为 **(N, D)** 的 `x_test`。 2. 定义您的 **协方差函数**。在下面的代码段中，我们使用径向基函数 (RBF) 核 “RBF” 以及使用 “WhiteKernel” 的加性噪声。 3. 使用您的协方差函数定义 “GaussianProcessRegressor” 对象，以及为您的 GPR 设定种子的随机状态。这种 “random_state” 对于确保再现性很重要。 4. 使用方法 `gpr.fit(x_train, y_train)` 拟合您的 `gpr` 对象。这会 “训练您的模型”，并使用梯度方法（例如基于 Hessian 的二阶优化例程 lbfgs）优化 gpr 对象的超参数。 5. 使用方法 `gpr.predict(x_test, return_std=True)` 预测测试点 `x_test` 上目标的 **均值** 和 **协方差**。这为您提供了预测值以及该预测点的不确定性度量。要使用 `pip` 为下面的示例安装依赖项：

1	pip install scikit-learn numpy matplotlib

这是一个使用 `sklearn` **拟合** 和 **预测** 一维正弦曲线的示例：

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import WhiteKernel, DotProduct, RBF
import matplotlib.pyplot as plt
import numpy as np

# Random seeds
np.random.seed(seed=0)  # Set seed for NumPy
random_state = 0

# Generate features, and take norm for use with target
x = np.random.normal(loc=0, scale=1, size=(50, 1))
y = np.sin(x)

# Create kernel and define GPR
kernel = RBF() + WhiteKernel()
gpr = GaussianProcessRegressor(kernel=kernel, random_state=random_state)

# Fit GPR model
gpr.fit(x, y)

# Create test data
x_test = np.random.normal(loc=0, scale=1, size=(50, 1))
y_test = np.sin(x_test)
 
# Predict mean
y_hat, y_sigma = gpr.predict(x_test, return_std=True)

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

# Squeeze data
x = np.squeeze(x)
y = np.squeeze(y)
x_test = np.squeeze(x_test)
y_test = np.squeeze(y_test)

# Plot the training data
ax.scatter(x, y)

# Plot predictive means as blue line
ax.plot(x_test, y_hat, 'b')

# Shade between the lower and upper confidence bounds
lower = x_test - y_sigma
upper = x_test + y_sigma
ax.fill_between(x_test, lower, upper, alpha=0.5)
ax.set_ylim([-3, 3])
plt.title("GPR Model Predictions")
plt.show()

### GPyTorch [GPyTorch](https://gpytorch.ai/)^[undefined] 非常适用于创建完全可定制的、先进的、可加速的高斯过程回归模型。这个软件包支持从高斯过程回归模型优化（通过自动分化）到硬件加速（通过 CUDA 和 [PyKeOps](https://pypi.org/project/pykeops/) 的一切。建议您在使用 GPyTorch 之前熟悉 PyTorch 和/或 python 中的自动微分包，但是 [**教程**](https://docs.gpytorch.ai/en/v1.1.1/examples/ 01_Exact_GPs/Simple_GP_Regression.html) 使该框架易于学习和使用。 GPyTorch 中 GPR 的数据表示为 “torch.tensor” 对象。以下是在 GPyTorch 中拟合高斯过程回归模型的步骤： 1. 预处理您的数据。训练数据可以表示为具有 `x_train` 形状 **(B, N, D)** 和 `y_train` 形状 **(B, N, 1)** 的 `(x_train, y_train)` 元组，其中 **B** 是批量大小，**N** 是样本数，**D** 是特征的维度。您的测试点可以表示为形状为 **(B, N, D)** 的 “x_test”。 2. 通过子类化 `gpytorch.models.ExactGP` 类来定义您的 `ExactGPModel`。要对该模型进行子类化，您需要定义：(i) 构造函数方法，它指定模型的均值和协方差函数，(ii) `forward` 方法，它描述了高斯过程回归模型如何进行预测。要使用 **批处理**，请查看教程 [**此处**](https://docs.gpytorch.ai/en/v1.1.1/examples/08_Advanced_Usage/index.html?highlight=batch#batch-gpr)。要在超参数上使用 **先验** 分布，请查看教程 [**此处**](https://docs.gpytorch.ai/en/v1.1.1/examples/00_Basic_Usage/Hyperparameters.html?highlight=priors)。 3. 指定您的“似然”函数，您的模型使用该函数将潜在变量 **f** 与观察到的目标 **y** 相关联。 4. 使用您的“可能性”和训练数据“(x_train, y_train)”实例化您的“模型”。 5. 使用 `pytorch` 自动微分对您的 `model` 进行超参数优化（“训练”）。完成后，确保您的 `model` 和 `likelihood` 使用 `model.eval()` 和 `likelihood.eval()` 置于后验模式。 6. 通过调用 `likelihood(model(x_test))` 使用您的 `model` 计算测试点的 **mean** 和 **variance** 预测。内部函数从测试输入“**x**” 预测潜在测试值 f，外部函数从潜在测试值预测均值和方差。要使用 `pip` 为下面的示例安装依赖项：

1	pip install gpytorch torch matplotlib numpy# (Optional) - Installs pykeopspip install pykeops

下面是一个使用 `gpytorch` 拟合含噪声一维正弦曲线的示例：

import gpytorch
import numpy as np
import matplotlib.pyplot as plt
import torch

# Create features
X = torch.tensor(np.arange(100)).float()

# Create targets as noisy function of features
Y = torch.tensor(np.add(np.sin(X / 10), np.random.normal(loc=0, scale=0.5, size=100))).float()

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

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

# Initialize likelihood and model
likelihood = gpytorch.likelihoods.GaussianLikelihood()
model = ExactGPModel(X, Y, likelihood)

# this is for running the notebook in our testing framework
import os
training_iter = 50

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

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

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

# Performs "training" via hyperparameter optimization
for i in range(training_iter):
    # Zero gradients from previous iteration
    optimizer.zero_grad()
    # Output from model
    output = model(X)
    # Calc loss and backprop gradients
    loss = -mll(output, Y)
    loss.backward()
    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()

# Places model in "posterior mode"
model = model.eval()

# Create test data
X_test = torch.tensor(np.arange(0, 100, 0.01)).float()

# Makes predictions without noise
f_preds = model(X_test)

# Attributes of predictive function
f_mean = f_preds.mean
f_var = f_preds.variance
f_covar = f_preds.covariance_matrix

# Run with the no_grad and LOVE (fast variance prediction) context manager
with torch.no_grad(), gpytorch.settings.fast_pred_var():
    
    # Initialize plot
    f, ax = plt.subplots(1, 1, figsize=(4, 3))
    
    # Make predictions with noise model
    observed_pred = likelihood(model(X_test))
    
    # Get upper and lower confidence bounds
    lower, upper = observed_pred.confidence_region()
    
    # Plot training data as black stars
    ax.scatter(X.numpy(), Y.numpy())
    
    # Plot predictive means as blue line
    ax.plot(X_test.numpy(), observed_pred.mean.numpy(), 'b')
    
    # Shade between the lower and upper confidence bounds
    ax.fill_between(X_test.numpy(), lower.numpy(), upper.numpy(), alpha=0.5)
    ax.set_ylim([-3, 3])

### GPFlow [GPFlow](https://gpflow.readthedocs.io/en/master/index.html) ^[2] 是另一个支持自动微分的高斯过程回归包，只是其后端采用了 `tensorflow`。GPFlow 具有广泛的内置功能，可用于创建完全可定制的模型、似然函数、内核以及优化和推断例程。除了 GPR，GPFlow 还内置了针对贝叶斯优化中其他各种最先进问题的功能，例如 [**变分傅里叶特征**](https://gpflow.readthedocs.io/en/master/notebooks/advanced/variational_fourier_features.html) 和 [**卷积高斯过程**](https://gpflow.readthedocs.io/en/master/notebooks/advanced/convolutional.html)。在使用 GPFlow 之前，建议您对 Python 中的 TensorFlow 和/或自动微分包有一定的了解。 GPFlow 中 GPR 的数据表示为 “tf.tensor” 对象。要开始使用 GPFlow，请查看 [**此示例链接**](https://gpflow.readthedocs.io/en/master/notebooks_file.html)。 ### GPy [GPy](https://gpy.readthedocs.io/en/deploy/GPy.models.html) ^[4] 含有大量高斯过程回归模型、似然函数和推断过程的 Python 实现。尽管这个包不支持自动微分的后端，但其多功能性、模块化和可定制性使其成为实现 GPR 的宝贵资源。 ### Pyro [Pyro](http://pyro.ai/) ^[6] 是一个概率编程包，可以与 Python 集成，还支持高斯过程回归，以及高级应用程序，例如 [深度核学习](https://pyro.ai/examples/dkl.html)。 ### Gen [Gen](https://www.gen.dev/) ^[7] 是一个建立在 Julia 之上的概率编程包。 Gen 通过高斯过程回归提供了几个优势：（i）它建立在提议分布中，这可以通过有效地对可能的解集合施加先验来帮助缩小搜索空间，（ii）它有一个简单的 API 用于采样跟踪来自适合高斯过程回归模型，(iii) 作为许多概率编程语言的目标，它能够轻松创建分层模型以调整高斯过程回归超参数的先验。 ### Stan [Stan](https://mc-stan.org/) ^[8] 是另一个可以与 Python 集成的概率编程包，但也支持其他语言，如 R、MATLAB、Julia 和 Stata。除了具有内置的高斯过程回归功能外，Stan 还支持各种其他贝叶斯推断和采样功能。 ### BoTorch [BoTorch](https://botorch.org/) ^[9] 由 GPyTorch 的创建者构建，是一个贝叶斯优化库，它支持许多与 GPyTorch 相同的高斯过程回归技术，以及高级贝叶斯优化技术和分析测试套件。 ## 6 总结在本文中，我们回顾了高斯过程回归 (GPR) 背后的理论，介绍并讨论了高斯过程回归可用于解决的问题类型，讨论了高斯过程回归与其他监督学习算法的比较，并介绍了我们如何能够使用 `sklearn`、`gpytorch` 或 `gpflow` 实现高斯过程回归。 ## 致谢 Thank you to Carl Edward Rasmussen for open-sourcing the textbook _Gaussian Processes for Machine Learning_ ^[5], and for [Scikit-Learn](https://scikit-learn.org/stable/), [GPyTorch](https://gpytorch.ai/), [GPFlow](https://gpflow.readthedocs.io/en/master/intro.html), and [GPy](https://gpy.readthedocs.io/en/deploy/GPy.models.html) for open-sourcing their Gaussian Process Regression Python libraries. ## 参考文献

[1] Pedregosa, Fabian, et al. “Scikit-learn: Machine learning in Python.” _the Journal of machine Learning research_ 12 (2011): 2825–2830.
[2] Gardner, Jacob R., et al. “Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration.” _arXiv preprint arXiv:1809.11165_ (2018).
[3] Matthews, Alexander G. de G., et al. “GPflow: A Gaussian Process Library using TensorFlow.” _J. Mach. Learn. Res._ 18.40 (2017): 1–6.
[4] GPy, “GPy.” [_http://github.com/SheffieldML/GPy_](http://github.com/SheffieldML/GPy)_._
[5] Carl Edward Rasmussen and Christopher K. I. Williams. 2005. Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning). The MIT Press.
[6] Eli Bingham, Jonathan P. Chen, Martin Jankowiak, Fritz Obermeyer, Neeraj Pradhan, Theofanis Karaletsos, Rohit Singh, Paul Szerlip, Paul Horsfall, and Noah D. Goodman. 2019. _Pyro: deep universal probabilistic programming._ J. Mach. Learn. Res. 20, 1 (January 2019), 973–978.
[7] Gen: A General-Purpose Probabilistic Programming System with Programmable Inference. Cusumano-Towner, M. F.; Saad, F. A.; Lew, A.; and Mansinghka, V. K. In Proceedings of the 40th ACM SIGPLAN Conference on Programming Language Design and Implementation (PLDI ‘19).
[8] Stan Development Team. 2021. Stan Modeling Language Users Guide and Reference Manual, VERSION. [https://mc-stan.org](https://mc-stan.org).
[9] Balandat, Maximilian, et al. “BoTorch: A framework for efficient Monte-Carlo Bayesian optimization.” _Advances in Neural Information Processing Systems_ 33 (2020).