18.3. 高斯过程推断¶ 在 SageMaker Studio Lab 中打开 Notebook
在本节中，我们将展示如何使用上一节介绍的GP先验执行后验推断和进行预测。我们将从回归开始，在回归中我们可以以*闭合形式*执行推断。这是一个“GP简介”部分，可以让你在实践中快速上手高斯过程。我们将从头开始编写所有基本操作的代码，然后介绍GPyTorch，它将使处理最先进的高斯过程以及与深度神经网络的集成变得更加方便。我们将在下一节中深入探讨这些更高级的主题。在该节中，我们还将考虑需要近似推断的设置——分类、点过程或任何非高斯似然。
18.3.1. 回归的后验推断¶
一个*观测*模型将我们想要学习的函数\(f(x)\)与我们的观测\(y(x)\)联系起来，两者都由某个输入\(x\)索引。在分类中，\(x\)可以是图像的像素，而\(y\)可以是相关的类别标签。在回归中，\(y\)通常表示一个连续的输出，例如地表温度、海平面、\(CO_2\)浓度等。
在回归中，我们通常假设输出由一个潜在的无噪声函数\(f(x)\)加上独立同分布的高斯噪声\(\epsilon(x)\)给出
其中\(\epsilon(x) \sim \mathcal{N}(0,\sigma^2)\)。令\(\mathbf{y} = y(X) = (y(x_1),\dots,y(x_n))^{\top}\)为我们训练观测的向量，\(\textbf{f} = (f(x_1),\dots,f(x_n))^{\top}\)为在训练输入\(X = {x_1, \dots, x_n}\)处查询的潜在无噪声函数值的向量。
我们将假设\(f(x) \sim \mathcal{GP}(m,k)\)，这意味着任何函数值集合\(\textbf{f}\)都有一个联合多元高斯分布，其均值向量为\(\mu_i = m(x_i)\)，协方差矩阵为\(K_{ij} = k(x_i,x_j)\)。RBF核\(k(x_i,x_j) = a^2 \exp\left(-\frac{1}{2\ell^2}||x_i-x_j||^2\right)\)是协方差函数的标准选择。为了符号简洁，我们将假设均值函数\(m(x)=0\)；我们的推导可以稍后轻松推广。
假设我们想要在一组输入上进行预测
然后我们想要找到\(x^2\)和\(p(\mathbf{f}_* | \mathbf{y}, X)\)。在回归设置中，在找到\(\mathbf{f}_* = f(X_*)\)和\(\mathbf{y}\)的联合分布后，我们可以方便地使用高斯恒等式来找到这个分布。
如果我们对训练输入\(X\)评估方程(18.3.1)，我们有\(\mathbf{y} = \mathbf{f} + \mathbf{\epsilon}\)。根据高斯过程的定义（见上一节），\(\mathbf{f} \sim \mathcal{N}(0,K(X,X))\)，其中\(K(X,X)\)是一个\(n \times n\)的矩阵，通过在所有可能的输入对\(x_i, x_j \in X\)上评估我们的协方差函数（也称*核*）形成。\(\mathbf{\epsilon}\)只是一个由来自\(\mathcal{N}(0,\sigma^2)\)的独立同分布样本组成的向量，因此其分布为\(\mathcal{N}(0,\sigma^2I)\)。\(\mathbf{y}\)因此是两个独立的多元高斯变量的和，因此其分布为\(\mathcal{N}(0, K(X,X) + \sigma^2I)\)。还可以证明\(\textrm{cov}(\mathbf{f}_*, \mathbf{y}) = \textrm{cov}(\mathbf{y},\mathbf{f}_*)^{\top} = K(X_*,X)\)，其中\(K(X_*,X)\)是一个\(m \times n\)的矩阵，通过在所有测试输入和训练输入对上评估核函数形成。
然后我们可以使用标准的高斯恒等式从联合分布中找到条件分布（例如，参见Bishop第2章），\(\mathbf{f}_* | \mathbf{y}, X, X_* \sim \mathcal{N}(m_*,S_*)\)，其中\(m_* = K(X_*,X)[K(X,X)+\sigma^2I]^{-1}\textbf{y}\)，并且\(S = K(X_*,X_*) - K(X_*,X)[K(X,X)+\sigma^2I]^{-1}K(X,X_*)\)。
通常，我们不需要使用完整的预测协方差矩阵\(S\)，而是使用\(S\)的对角线来表示每个预测的不确定性。因此，我们经常为单个测试点\(x_*\)写出预测分布，而不是为一组测试点。
核矩阵具有我们希望估计的参数\(\theta\)，例如上面RBF核的幅度\(a\)和长度尺度\(\ell\)。为此，我们使用*边际似然*\(p(\textbf{y} | \theta, X)\)，我们已经在推导联合分布\(\textbf{y},\textbf{f}_*\)时推导出了它。我们将看到，边际似然被分解为模型拟合项和模型复杂性项，并自动编码了奥卡姆剃刀原理来学习超参数。有关全面讨论，请参见MacKay第28章 (MacKay, 2003), 和 Rasmussen and Williams 第5章 (Rasmussen and Williams, 2006)。
import math
import os
import gpytorch
import matplotlib.pyplot as plt
import numpy as np
import torch
from scipy import optimize
from scipy.spatial import distance_matrix
from d2l import torch as d2l
d2l.set_figsize()
18.3.2. 在GP回归中进行预测和学习核超参数的方程¶
我们在此列出了您将在高斯过程回归中用于学习超参数和进行预测的方程。同样，我们假设一个回归目标向量\(\textbf{y}\)，由输入\(X = \{x_1,\dots,x_n\}\)索引，我们希望在一个测试输入\(x_*\)处进行预测。我们假设存在独立同分布的零均值高斯噪声，方差为\(\sigma^2\)。我们使用高斯过程先验\(f(x) \sim \mathcal{GP}(m,k)\)来描述潜在的无噪声函数，其中均值函数为\(m\)，核函数为\(k\)。核函数本身有我们想要学习的参数\(\theta\)。例如，如果我们使用RBF核，\(k(x_i,x_j) = a^2\exp\left(-\frac{1}{2\ell^2}||x-x'||^2\right)\)，我们想要学习\(\theta = \{a^2, \ell^2\}\)。令\(K(X,X)\)表示一个\(n \times n\)的矩阵，对应于对所有可能的\(n\)个训练输入对评估核函数。令\(K(x_*,X)\)表示一个\(1 \times n\)的向量，通过评估\(k(x_*, x_i)\)，\(i=1,\dots,n\)形成。令\(\mu\)为一个均值向量，通过在每个训练点\(x\)处评估均值函数\(m(x)\)形成。
通常在使用高斯过程时，我们遵循一个两步程序。1. 通过最大化边际似然来学习核超参数\(\hat{\theta}\)。2. 使用预测均值作为点预测器，并使用2倍的预测标准差来形成一个95%的置信区间，条件是这些学习到的超参数\(\hat{\theta}\)。
对数边际似然只是一个对数高斯密度，其形式为
预测分布的形式为
18.3.3. 解释学习和预测的方程¶
关于高斯过程的预测分布，有几个关键点需要注意
尽管模型类别非常灵活，但GP回归可以在*闭合形式*中进行*精确*的贝叶斯推断。除了学习核超参数，没有*训练*过程。我们可以精确地写下我们想要用来做预测的方程。高斯过程在这方面相对特殊，这极大地促进了它们的便利性、多功能性和持续的受欢迎程度。
预测均值\(a_*\)是训练目标\(\textbf{y}\)的线性组合，由核函数\(k_{\theta}(x_*,X)[K_{\theta}(X,X)+\sigma^2I]^{-1}\)加权。我们将看到，核函数（及其超参数）因此在模型的泛化性能中扮演着至关重要的角色。
预测均值明确地依赖于目标值\(\textbf{y}\)，但预测方差不依赖。预测不确定性随着测试输入\(x_*\)远离目标位置\(X\)而增加，这由核函数决定。然而，不确定性会通过从数据中学习到的核超参数\(\theta\)间接依赖于目标值\(\textbf{y}\)。
边际似然分解为模型拟合项和模型复杂度（对数行列式）项。边际似然倾向于选择那些能够提供与数据一致的最简单拟合的超参数。
主要的计算瓶颈来自于对\(n\)个训练点的\(n \times n\)对称正定矩阵\(K(X,X)\)求解线性系统和计算对数行列式。朴素地讲，这些操作每个都需要\(\mathcal{O}(n^3)\)的计算量，以及\(\mathcal{O}(n^2)\)的存储空间来存放核（协方差）矩阵的每个元素，通常从Cholesky分解开始。历史上，这些瓶颈将GP限制在训练点少于约10,000个的问题上，并给GP带来了“慢”的声誉，而这个声誉近十年来已经不再准确。在高级主题中，我们将讨论如何将GP扩展到拥有数百万个点的问题上。
对于常见的核函数选择，\(K(X,X)\)通常接近奇异，这在执行Cholesky分解或其他旨在解决线性系统的操作时可能导致数值问题。幸运的是，在回归中我们通常处理\(K_{\theta}(X,X)+\sigma^2I\)，这样噪声方差\(\sigma^2\)就被加到\(K(X,X)\)的对角线上，显著改善了其条件数。如果噪声方差很小，或者我们正在进行无噪声回归，通常的做法是在对角线上添加少量“抖动”，数量级在\(10^{-6}\)左右，以改善条件数。
18.3.4. 从零开始的实例¶
让我们创建一些回归数据，然后用GP拟合数据，从头开始实现每一步。我们将从以下分布中采样数据
其中 \(\epsilon \sim \mathcal{N}(0,\sigma^2)\)。我们希望找到的无噪声函数是 \(f(x) = \sin(x) + \frac{1}{2}\sin(4x)\)。我们将从噪声标准差 \(\sigma = 0.25\) 开始。
def data_maker1(x, sig):
return np.sin(x) + 0.5 * np.sin(4 * x) + np.random.randn(x.shape[0]) * sig
sig = 0.25
train_x, test_x = np.linspace(0, 5, 50), np.linspace(0, 5, 500)
train_y, test_y = data_maker1(train_x, sig=sig), data_maker1(test_x, sig=0.)
d2l.plt.scatter(train_x, train_y)
d2l.plt.plot(test_x, test_y)
d2l.plt.xlabel("x", fontsize=20)
d2l.plt.ylabel("Observations y", fontsize=20)
d2l.plt.show()
这里我们看到带噪声的观测点是圆圈，而我们希望找到的无噪声函数是蓝色的。
现在，让我们在潜在的无噪声函数上指定一个GP先验，\(f(x)\sim \mathcal{GP}(m,k)\)。我们将使用均值函数\(m(x) = 0\)，以及一个RBF协方差函数（核）
mean = np.zeros(test_x.shape[0])
cov = d2l.rbfkernel(test_x, test_x, ls=0.2)
我们从长度尺度0.2开始。在拟合数据之前，重要的是要考虑我们是否指定了一个合理的先验。让我们可视化一些来自这个先验的样本函数，以及95%的置信区间（我们相信真实函数有95%的几率在这个区域内）。
prior_samples = np.random.multivariate_normal(mean=mean, cov=cov, size=5)
d2l.plt.plot(test_x, prior_samples.T, color='black', alpha=0.5)
d2l.plt.plot(test_x, mean, linewidth=2.)
d2l.plt.fill_between(test_x, mean - 2 * np.diag(cov), mean + 2 * np.diag(cov),
alpha=0.25)
d2l.plt.show()
这些样本看起来合理吗？这些函数的高级属性是否与我们试图建模的数据类型一致？
现在让我们在任意测试点\(x_*\)处形成后验预测分布的均值和方差。
在进行预测之前，我们应该学习我们的核超参数\(\theta\)和噪声方差\(\sigma^2\)。让我们将长度尺度初始化为0.75，因为我们的先验函数看起来比我们正在拟合的数据变化得太快。我们还会猜测一个噪声标准差\(\sigma\)为0.75。
为了学习这些参数，我们将最大化关于这些参数的边际似然。
也许我们的先验函数变化太快了。让我们猜测一个长度尺度为0.4。我们还会猜测一个噪声标准差为0.75。这些仅仅是超参数的初始化——我们将从边际似然中学习这些参数。
ell_est = 0.4
post_sig_est = 0.5
def neg_MLL(pars):
K = d2l.rbfkernel(train_x, train_x, ls=pars[0])
kernel_term = -0.5 * train_y @ \
np.linalg.inv(K + pars[1] ** 2 * np.eye(train_x.shape[0])) @ train_y
logdet = -0.5 * np.log(np.linalg.det(K + pars[1] ** 2 * \
np.eye(train_x.shape[0])))
const = -train_x.shape[0] / 2. * np.log(2 * np.pi)
return -(kernel_term + logdet + const)
learned_hypers = optimize.minimize(neg_MLL, x0=np.array([ell_est,post_sig_est]),
bounds=((0.01, 10.), (0.01, 10.)))
ell = learned_hypers.x[0]
post_sig_est = learned_hypers.x[1]
在这种情况下，我们学习到一个长度尺度为0.299，噪声标准差为0.24。请注意，学习到的噪声非常接近真实的噪声，这有助于表明我们的GP对这个问题拟合得非常好。
总的来说，仔细考虑选择核函数和初始化超参数至关重要。虽然边际似然优化对初始化相对稳健，但它并非对差的初始化免疫。尝试用各种初始化运行上述脚本，看看你得到了什么结果。
现在，让我们用这些学习到的超参数进行预测。
K_x_xstar = d2l.rbfkernel(train_x, test_x, ls=ell)
K_x_x = d2l.rbfkernel(train_x, train_x, ls=ell)
K_xstar_xstar = d2l.rbfkernel(test_x, test_x, ls=ell)
post_mean = K_x_xstar.T @ np.linalg.inv((K_x_x + \
post_sig_est ** 2 * np.eye(train_x.shape[0]))) @ train_y
post_cov = K_xstar_xstar - K_x_xstar.T @ np.linalg.inv((K_x_x + \
post_sig_est ** 2 * np.eye(train_x.shape[0]))) @ K_x_xstar
lw_bd = post_mean - 2 * np.sqrt(np.diag(post_cov))
up_bd = post_mean + 2 * np.sqrt(np.diag(post_cov))
d2l.plt.scatter(train_x, train_y)
d2l.plt.plot(test_x, test_y, linewidth=2.)
d2l.plt.plot(test_x, post_mean, linewidth=2.)
d2l.plt.fill_between(test_x, lw_bd, up_bd, alpha=0.25)
d2l.plt.legend(['Observed Data', 'True Function', 'Predictive Mean', '95% Set on True Func'])
d2l.plt.show()
我们看到橙色的后验均值几乎完美地匹配了真实的无噪声函数！请注意，我们显示的95%置信区间是针对潜在的*无噪声*（真实）函数的，而不是数据点。我们看到这个置信区间完全包含了真实函数，并且看起来既不过宽也不过窄。我们既不希望也不期望它包含数据点。如果我们希望得到观测值的置信区间，我们应该计算
lw_bd_observed = post_mean - 2 * np.sqrt(np.diag(post_cov) + post_sig_est ** 2)
up_bd_observed = post_mean + 2 * np.sqrt(np.diag(post_cov) + post_sig_est ** 2)
不确定性有两个来源，*认知*不确定性，代表*可减少*的不确定性，以及*偶然*或*不可减少*的不确定性。这里的*认知*不确定性代表了关于无噪声函数真实值的不确定性。当我们远离数据点时，这种不确定性应该会增加，因为远离数据点时，与我们的数据一致的函数值种类更多。随着我们观察到越来越多的数据，我们对真实函数的信念变得更加自信，认知不确定性就会消失。在这种情况下，*偶然*不确定性是观测噪声，因为数据是以这种噪声给我们的，并且它无法被减少。
数据中的*认知*不确定性由潜在无噪声函数的方差 np.diag(post_cov) 捕获。*偶然*不确定性由噪声方差 post_sig_est**2 捕获。
不幸的是，人们常常不注意他们如何表示不确定性，许多论文展示的误差棒是完全未定义的，没有明确说明我们是在可视化认知不确定性、偶然不确定性还是两者兼而有之，并且混淆了噪声方差与噪声标准差、标准差与标准误、置信区间与可信区间等等。如果不精确地说明不确定性代表什么，它基本上是无意义的。
本着密切关注我们的不确定性代表什么的原则，至关重要的是要注意到我们正在取我们对无噪声函数的方差估计的*平方根*的*两倍*。由于我们的预测分布是高斯的，这个量使我们能够形成一个95%的可信区间，代表了我们对该区间有95%的可能包含真实函数的信念。噪声*方差*处于完全不同的尺度，并且解释性要差得多。
最后，让我们看20个后验样本。这些样本告诉我们，我们认为哪些类型的函数可能拟合我们的数据，即后验地。
post_samples = np.random.multivariate_normal(post_mean, post_cov, size=20)
d2l.plt.scatter(train_x, train_y)
d2l.plt.plot(test_x, test_y, linewidth=2.)
d2l.plt.plot(test_x, post_mean, linewidth=2.)
d2l.plt.plot(test_x, post_samples.T, color='gray', alpha=0.25)
d2l.plt.fill_between(test_x, lw_bd, up_bd, alpha=0.25)
plt.legend(['Observed Data', 'True Function', 'Predictive Mean', 'Posterior Samples'])
d2l.plt.show()
在基本的回归应用中，最常见的是使用后验预测均值和标准差分别作为点预测器和不确定性度量。在更高级的应用中，例如使用蒙特卡洛采集函数的贝叶斯优化，或者用于基于模型强化学习的高斯过程，通常需要采集后验样本。然而，即使在基本应用中并非严格要求，这些样本也能让我们对数据的拟合有更多的直觉，并且通常在可视化中很有用。
18.3.5. 用GPyTorch让事情变简单¶
正如我们所见，从头开始实现基本的高斯过程回归实际上非常简单。然而，一旦我们想要探索各种核函数的选择，考虑近似推断（即使是分类也需要），将GP与神经网络结合，或者甚至有一个超过约10,000个点的数据集，那么从头开始的实现就会变得笨拙和繁琐。一些最有效的可扩展GP推断方法，如SKI（也称为KISS-GP），可能需要数百行代码来实现高级的数值线性代数例程。
在这些情况下，*GPyTorch*库将使我们的生活轻松得多。我们将在未来关于高斯过程数值计算和高级方法的笔记中更多地讨论GPyTorch。GPyTorch库包含许多示例。为了对这个包有一个感觉，我们将通过简单的回归示例，展示如何调整它来使用GPyTorch重现我们上面的结果。这可能看起来是很多代码来简单地重现上面的基本回归，在某种意义上，确实是这样。但是我们可以立即使用各种核函数、可扩展的推断技术和近似推断，只需从下面更改几行代码，而不需要编写可能数千行的新代码。
# First let's convert our data into tensors for use with PyTorch
train_x = torch.tensor(train_x)
train_y = torch.tensor(train_y)
test_y = torch.tensor(test_y)
# We are using exact GP inference with a zero mean and RBF kernel
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.ZeroMean()
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)
此代码块将数据放入 GPyTorch 的正确格式，并指定我们正在使用精确推断，以及我们想要使用的均值函数（零）和核函数（RBF）。我们可以非常容易地使用任何其他核函数，例如，通过调用 gpytorch.kernels.matern_kernel() 或 gpytorch.kernels.spectral_mixture_kernel()。到目前为止，我们只讨论了精确推断，在这种情况下，可以无需任何近似就能推断出预测分布。对于高斯过程，只有当似然函数是高斯时，我们才能执行精确推断；更具体地说，当我们假设我们的观测是由高斯过程表示的无噪声函数加上高斯噪声生成时。在未来的笔记中，我们将考虑其他情况，例如分类，在这些情况下我们不能做出这些假设。
# Initialize Gaussian likelihood
likelihood = gpytorch.likelihoods.GaussianLikelihood()
model = ExactGPModel(train_x, train_y, likelihood)
training_iter = 50
# Find optimal model hyperparameters
model.train()
likelihood.train()
# Use the adam optimizer, includes GaussianLikelihood parameters
optimizer = torch.optim.Adam(model.parameters(), lr=0.1)
# Set our loss as the negative log GP marginal likelihood
mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)
在这里，我们明确指定了我们想要使用的似然（高斯），用于训练核超参数的目标（这里是边际似然），以及我们想要用于优化该目标的过程（在这种情况下是Adam）。我们注意到，虽然我们正在使用Adam，这是一个“随机”优化器，但在这种情况下，它是全批量Adam。因为边际似然不能在数据实例上分解，我们不能在数据的“小批量”上使用优化器并保证收敛。GPyTorch也支持其他优化器，例如L-BFGS。与标准深度学习不同，优化边际似然与良好的泛化性能密切相关，这通常使我们倾向于使用强大的优化器，如L-BFGS，前提是它们不是过于昂贵。
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 % 10 == 0:
print(f'Iter {i+1:d}/{training_iter:d} - Loss: {loss.item():.3f} '
f'squared lengthscale: '
f'{model.covar_module.base_kernel.lengthscale.item():.3f} '
f'noise variance: {model.likelihood.noise.item():.3f}')
optimizer.step()
Iter 1/50 - Loss: 1.000 squared lengthscale: 0.693 noise variance: 0.693
Iter 11/50 - Loss: 0.711 squared lengthscale: 0.490 noise variance: 0.312
Iter 21/50 - Loss: 0.451 squared lengthscale: 0.506 noise variance: 0.127
Iter 31/50 - Loss: 0.330 squared lengthscale: 0.485 noise variance: 0.055
Iter 41/50 - Loss: 0.344 squared lengthscale: 0.472 noise variance: 0.038
在这里，我们实际运行优化过程，每10次迭代输出一次损失值。
# Get into evaluation (predictive posterior) mode
test_x = torch.tensor(test_x)
model.eval()
likelihood.eval()
observed_pred = likelihood(model(test_x))
上面的代码块使我们能够在测试输入上进行预测。
with torch.no_grad():
# Initialize plot
f, ax = d2l.plt.subplots(1, 1, figsize=(4, 3))
# Get upper and lower bounds for 95\% credible set (in this case, in
# observation space)
lower, upper = observed_pred.confidence_region()
ax.scatter(train_x.numpy(), train_y.numpy())
ax.plot(test_x.numpy(), test_y.numpy(), linewidth=2.)
ax.plot(test_x.numpy(), observed_pred.mean.numpy(), linewidth=2.)
ax.fill_between(test_x.numpy(), lower.numpy(), upper.numpy(), alpha=0.25)
ax.set_ylim([-1.5, 1.5])
ax.legend(['True Function', 'Predictive Mean', 'Observed Data',
'95% Credible Set'])
最后，我们绘制拟合图。
我们看到拟合结果几乎完全相同。有几点需要注意：GPyTorch处理的是*平方*长度尺度和观测噪声。例如，我们在从零开始的代码中学到的噪声标准差约为0.283。GPyTorch找到的噪声方差是\(0.81 \approx 0.283^2\)。在GPyTorch的图中，我们还显示了*观测空间*中的置信区间，而不是潜在函数空间，以证明它们确实覆盖了观测到的数据点。
18.3.6. 总结¶
我们可以将高斯过程先验与数据结合形成后验，我们用它来进行预测。我们还可以形成边际似然，这对于自动学习核超参数很有用，这些超参数控制着高斯过程的变化率等属性。为回归形成后验和学习核超参数的机制很简单，大约只需要十几行代码。这个笔记对于任何希望快速“上手”高斯过程的读者来说是一个很好的参考。我们还介绍了 GPyTorch 库。虽然用于基本回归的 GPyTorch 代码相对较长，但它可以很容易地修改以用于其他核函数，或我们将在未来笔记中讨论的更高级功能，例如可扩展的推断或用于分类的非高斯似然。
18.3.7. 练习¶
我们强调了*学习*核超参数的重要性，以及超参数和核对高斯过程泛化性能的影响。尝试跳过我们学习超参数的步骤，而是猜测各种长度尺度和噪声方差，并检查它们对预测的影响。当你使用大的长度尺度时会发生什么？小的长度尺度呢？大的噪声方差？小的噪声方差？
我们已经说过，边际似然不是一个凸目标函数，但在GP回归中，像长度尺度和噪声方差这样的超参数可以被可靠地估计出来。这通常是正确的——事实上，边际似然在学习长度尺度超参数方面比空间统计中的传统方法要好得多，后者涉及拟合经验自相关函数（“协变图”）。可以说，机器学习对高斯过程研究的最大贡献，至少在最近关于可扩展推断的工作之前，是引入了边际似然用于超参数学习。
然而，即使是这些参数的不同配对也会为许多数据集提供可解释的不同合理解释，导致我们的目标函数出现局部最优解。如果我们使用一个大的长度尺度，那么我们假设真实的底层函数是缓慢变化的。如果观察到的数据确实变化显著，那么我们能够合理解释大长度尺度的唯一方法就是有一个大的噪声方差。另一方面，如果我们使用一个小的长度尺度，我们的拟合将对数据的变化非常敏感，几乎没有空间用噪声（偶然不确定性）来解释变化。
试着看看你是否能找到这些局部最优解：用非常大的长度尺度和大的噪声进行初始化，以及用小的长度尺度和小的噪声进行初始化。你会收敛到不同的解吗？
我们已经说过，贝叶斯方法的一个基本优势是能够自然地表示*认知*不确定性。在上面的例子中，我们不能完全看到认知不确定性的影响。尝试用
test_x = np.linspace(0, 10, 1000)来进行预测。当你的预测超出数据范围时，95%的置信区间会发生什么变化？它是否覆盖了那个区间的真实函数？如果你在该区域只可视化偶然不确定性会发生什么？尝试运行上述示例，但分别使用10,000、20,000和40,000个训练点，并测量运行时间。训练时间如何扩展？或者，运行时间如何随着测试点的数量扩展？对于预测均值和预测方差，情况是否不同？通过理论上计算训练和测试的时间复杂度，以及通过用不同数量的点运行上述代码来回答这个问题。
尝试使用不同的协方差函数运行GPyTorch示例，例如Matern核。结果如何变化？GPyTorch库中的谱混合核怎么样？有些核比其他核更容易训练边际似然吗？有些核对长程预测比短程预测更有价值吗？
在我们的GPyTorch示例中，我们绘制了包含观测噪声的预测分布，而在我们的“从零开始”的示例中，我们只包含了认知不确定性。重新做一遍GPyTorch示例，但这次只绘制认知不确定性，并与从零开始的结果进行比较。预测分布现在看起来是否相同？（它们应该相同。）