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

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

在本文中，我们将讨论：

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

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

2 高斯过程回归的概念

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

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

Fig

Fig

Fig

Fig

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

fig

fig

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

Fig

Fig

Fig

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

Fig

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

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

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

（A）插值/克里金法

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

图

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

（B）时间序列预测

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

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

图

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

（C） 预测不确定性

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

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

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

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

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

  1. 高斯过程回归是 非参数模型。这意味着它主要从数据本身学习，而不是通过学习大量参数。这是特别有利的，因为这导致高斯过程回归模型不像神经网络等高参数模型那样 数据饥渴，即它们不需要那么多的样本来实现强大的通用性。

  2. 对于插值和预测任务，高斯过程回归估计 预期值 和 不确定性。这对于在决策时考虑到这种不确定性的决策系统尤其有益。

  3. 高斯过程回归是一个 线性平滑器 [5]， 从 监督学习 角度来看，这可以概念化为一种 正则化 技术。从 贝叶斯 的角度来看，这相当于在您的模型上强加一个 先验，即测试点上的所有目标必须是现有训练目标的 线性组合。这个性质有助于高斯过程回归泛化到未观测到的数据，只要真正为观测到的目标可以表示为训练目标的线性组合。

  4. 通过集成 torch 和 tensorflow 等 自动微分 后端框架，gpytorch 和 gpflow 等高斯过程回归的软件包正在变得 更快 和 根据可扩展性。对于 批次 模型尤其如此。有关这方面的示例案例研究，请参阅我之前关于小批量、多维高斯过程回归的文章此处！

5 如何实现高斯过程回归？

下面，我们介绍几个 Python 机器学习包，用于 可扩展的、高效的、模块化的 高斯过程回归实现。

Scikit-Learn

Scikit-Learn [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 拟合 和 预测 一维正弦曲线的示例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
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[undefined] 非常适用于创建完全可定制的、先进的、可加速的高斯过程回归模型。这个软件包支持从高斯过程回归模型优化（通过自动分化）到硬件加速（通过 CUDA 和 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 方法，它描述了高斯过程回归模型如何进行预测。要使用 批处理，请查看教程 此处。要在超参数上使用 先验 分布，请查看教程 此处。

  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 拟合含噪声一维正弦曲线的示例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
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 [2] 是另一个支持自动微分的高斯过程回归包，只是其后端采用了 tensorflow。GPFlow 具有广泛的内置功能，可用于创建完全可定制的模型、似然函数、内核以及优化和推断例程。除了 GPR，GPFlow 还内置了针对贝叶斯优化中其他各种最先进问题的功能，例如 变分傅里叶特征 和 卷积高斯过程。

在使用 GPFlow 之前，建议您对 Python 中的 TensorFlow 和/或自动微分包有一定的了解。 GPFlow 中 GPR 的数据表示为 “tf.tensor” 对象。要开始使用 GPFlow，请查看 此示例链接。

GPy

GPy [4] 含有大量高斯过程回归模型、似然函数和推断过程的 Python 实现。尽管这个包不支持自动微分的后端，但其多功能性、模块化和可定制性使其成为实现 GPR 的宝贵资源。

Pyro

Pyro [6] 是一个概率编程包，可以与 Python 集成，还支持高斯过程回归，以及高级应用程序，例如 深度核学习。

Gen

Gen [7] 是一个建立在 Julia 之上的概率编程包。 Gen 通过高斯过程回归提供了几个优势：（i）它建立在提议分布中，这可以通过有效地对可能的解集合施加先验来帮助缩小搜索空间，（ii）它有一个简单的 API 用于采样跟踪来自适合高斯过程回归模型，(iii) 作为许多概率编程语言的目标，它能够轻松创建分层模型以调整高斯过程回归超参数的先验。

Stan

Stan [8] 是另一个可以与 Python 集成的概率编程包，但也支持其他语言，如 R、MATLAB、Julia 和 Stata。除了具有内置的高斯过程回归功能外，Stan 还支持各种其他贝叶斯推断和采样功能。

BoTorch

BoTorch [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, GPyTorch, GPFlow, and GPy 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).