贝叶斯优化实战（三）

原文：zh.annas-archive.org/md5/5ea9246d4fedff957a1c887275e38f01

译者：飞龙

协议：CC BY-NC-SA 4.0

第十一章：通过偏好优化进行成对比较学习

本章涵盖

仅使用成对比较数据学习和优化偏好的问题
在成对比较上训练 GP
成对比较的优化策略

你是否曾经发现难以为某物（食物、产品或体验）打分？在 A/B 测试和产品推荐工作流程中，询问客户对产品的数值评分是一个常见任务。

定义术语 A/B 测试 指的是通过随机实验在两个环境（称为 A 和 B）中测量用户体验，并确定哪个环境更理想的方法。 A/B 测试通常由技术公司进行。

A/B 测试人员和产品推荐工程师经常需要处理从客户收集的反馈中的高水平噪音。通过噪音，我们指的是客户反馈所受到的任何类型的损坏。产品评分中的噪音示例包括在线流媒体服务上提供的广告数量，包裹的送货服务质量，或客户在消费产品时的一般心情。这些因素影响客户对产品的评分方式，可能会破坏客户对产品的真实评价。

不可控的外部因素使客户难以报告他们对产品的真实评价。因此，当在评分时难以选择数值评分作为对产品的评价时，客户通常会发现这很困难。在 A/B 测试和产品推荐中的反馈噪音的普遍存在意味着服务平台不能仅依靠从用户那里收集到的少量数据来了解其偏好。相反，平台需要从客户那里收集更多数据，以更加确定客户真正想要的是什么。

然而，正如在其他黑盒优化设置中一样，比如超参数调整和药物发现，查询目标函数是昂贵的。在产品推荐中，每当我们询问客户对产品的评分时，我们都面临着侵入客户体验和阻止他们继续使用平台的风险。因此，需要大量数据来更好地了解客户的偏好和具有侵入性之间存在自然紧张关系，可能导致客户流失。

幸运的是，有一种方法可以解决这个问题。心理学领域的研究（mng.bz/0KOl）发现了一个直观的结果，即我们人类在进行成对比较的偏好反应方面要比在评分产品时更擅长（例如，“产品 A 比产品 B 更好”）。

定义成对比较是一种收集偏好数据的方法。每次我们想要获取有关客户偏好的信息时，我们都会要求客户从两个项目中选择他们更喜欢的项目。成对比较不同于数值评分，其中我们要求客户在一定比例上对项目进行评分。

成对比较和评分之间难度差异的原因在于，比较两个项目是一项认知要求较低的任务，因此，我们可以在比较两个对象时更好地与我们的真实偏好保持一致，而不是提供数值评分。在图 10.1 中，考虑一个在线购物网站的两个示例界面，该网站正在尝试了解您对夏威夷衬衫的偏好：

第一个界面要求您按照从 1 到 10 的比例为衬衫评分。这可能很难做到，特别是如果您没有一个参考框架的话。
第二个界面要求您选择您更喜欢的衬衫。这个任务更容易完成。

图 10.1 生产推荐中用户偏好引导的示例。左侧，用户被要求对推荐产品进行评分。右侧，用户被要求选择他们更喜欢的产品。后者有助于更好地引导用户的偏好。

鉴于我们可以使用成对比较收集高质量数据的潜力，我们希望将这种偏好引导技术应用于用户偏好的 BayesOpt。问题是，“我们如何在成对比较数据上训练 ML 模型，然后，如何向用户呈现新的比较以最好地学习和优化他们的偏好？”我们在本章中回答了这些问题，首先使用一个能够有效地从成对比较中学习的 GP 模型。然后，我们开发了策略，将迄今为止我们找到的最佳数据点（代表一个产品）与一个有希望的竞争对手相比较，从而使我们能够尽快优化用户的偏好。换句话说，我们假设用户的偏好是定义在搜索空间上的目标函数，我们希望优化这个目标函数。

从成对比较中学习和优化用户偏好的这种设置是一个独特的任务，位于黑盒优化和产品推荐的交集处，两个社区都对此产生了兴趣。通过本章末尾，我们了解了如何从 BayesOpt 的角度解决这个问题，通过从用户那里收集数据来权衡开发和探索。

10.1 使用成对比较进行黑盒优化

在这一部分中，我们进一步讨论了成对比较在引导偏好任务中的有用性。然后，我们研究了为这种基于偏好的优化设置修改过的 BayesOpt 循环。

除了将精确的数字评估作为评级外，成对比较还提供了一种在生产推荐应用中收集有关客户信息的方法。与数字评级相比，成对比较对用户的认知负担较小，因此可能会产生更高质量的数据（即与用户真实偏好一致的反馈）。

多目标优化中的成对比较

成对比较特别有用的一个场景是在需要考虑多个标准的决策中。例如，假设你想买一辆车，正在选择 A 车和 B 车。为了做出决定，你列出了你在一辆车上关心的不同特征：外观、实用性、能效、成本等。然后你为两辆车在每个标准上评分，希望找到一个明显的赢家。不幸的是，A 车在某些标准上得分比 B 车高，但并不是所有标准都是如此，而 B 车在其余标准上的得分高于 A 车。

所以，在这两辆车之间没有明显的赢家，将不同标准的分数合并成一个单一分数可能会很困难。你关心某些标准胜过其他标准，所以在将这些标准与其他标准相结合以产生单一数字时，需要更加重视这些标准的权重。然而，确定这些权重的确切值可能比选择这两辆车本身更具挑战性！忽视具体细节，将每辆车作为一个整体，并将两辆车进行“头对头”比较可能更容易一些。

因此，在需要考虑多个标准的优化情况下，利用成对比较的便利性已经被利用起来。例如，Edward Abel、Ludmil Mikhailov 和 John Keane 的一个研究项目 (mng.bz/KenZ) 使用成对比较来解决群体决策问题。

当然，成对比较并不一定比数字评估更好。虽然前者更容易从用户那里获取，但它们所包含的信息明显比后者少得多。比如，在图 10.1 中，你喜欢橙色衬衫胜过红色衬衫的回答正好包含一位信息（比较的结果是二元的；要么橙色比红色好，要么红色比橙色好，所以观察结果信息理论上构成了一位信息）。而如果你报告说你给橙色衬衫评了 8 分，红色衬衫评了 6 分，那么我们获得的信息就比仅仅知道橙色被高估要多得多。

换句话说，在选择从用户那里引出反馈的方法时总是存在权衡。数值评估包含更多信息，但容易受到噪声影响，并且可能给用户带来更大的认知负担。另一方面，两两比较提供的信息较少，但用户报告起来更容易。这些优缺点在图 10.2 中总结。

图 10.2 数值评分和两两比较在信息量和报告难度方面的差异。每种偏好引出方法都有其优缺点。

在考虑信息和报告难度之间的权衡时，如果我们愿意让用户完成更具认知要求的任务以获取更多信息，并且可以考虑到噪声，那么我们应该坚持使用数值评估。另一方面，如果我们更注重客户准确表达他们真实的偏好，并且愿意获取更少的信息，那么两两比较应该是我们引出客户反馈的首选方法。

其他引出客户偏好的方法

两两比较并不是减轻数值评估认知负担的唯一形式。例如，在线流媒体服务 Netflix 通过要求观众在三个选项中进行选择来收集观众的评分：“向下拇指”表示他们不喜欢某物，“向上拇指”表示他们喜欢某物，“双向上拇指”表示他们喜欢某物 (mng.bz/XNgl)。这种设置构成了一种有序分类问题，其中项目被分类到不同的类别中，并且类别之间存在固有的顺序。在这种情况下，产品推荐问题同样值得考虑，但在本章中我们将重点放在两两比较上。

在本章中，我们学习如何利用 BayesOpt 来促进使用两两比较来学习和优化客户偏好的任务。首先，我们研究了一个修改过的 BayesOpt 循环版本，如图 1.6 所示，如图 10.3 所示：

在第一步中，GP 是根据两两比较数据而不是数值评估进行训练的。关键挑战在于确保 GP 对于目标函数（用户真实偏好函数）的信念反映了观察到的比较中的信息。
在第二步中，BayesOpt 策略计算获取分数，以量化对用户每个潜在新查询的有用程度。用户的查询需要以一对产品的形式提供给用户进行比较。就像在其他情况下一样，策略需要平衡利用我们知道用户偏好高的区域和探索我们对用户偏好了解不多的其他区域。
在第 3 步中，用户比较了由贝叶斯优化策略呈现给他们的两种产品，并报告他们更喜欢的产品。然后，将此新信息添加到我们的训练集中。

图 10.3 使用成对比较进行偏好优化的贝叶斯优化循环。高斯过程根据成对比较数据进行训练，而贝叶斯优化策略决定应该要求用户比较哪一对数据点。

我们在本章的剩余部分试图解决两个主要问题：

我们如何仅根据成对比较训练高斯过程？高斯过程在数值响应上进行训练时，会产生具有量化不确定性的概率预测，这在决策中至关重要。我们能否在这里使用相同的模型来处理成对比较响应？
我们应该如何生成新的产品对供用户比较，以便尽快确定用户偏好的最大化者？也就是说，我们如何通过成对比较最好地引出用户的反馈以优化他们的偏好？

10.2 制定偏好优化问题和格式化成对比数据

在我们开始解决这些问题之前，本节介绍了我们将在整章中解决的产品推荐问题以及我们如何在 Python 中模拟这个问题。正确设置问题将帮助我们更轻松地整合我们将在随后章节学习到的贝叶斯优化工具。我们在此处使用的代码包含在 CH10/01 - 从成对比较中学习.ipynb Jupyter 笔记本的第一部分中。

正如图 10.1 和 10.3 所示的，我们现在面临的情景是夏威夷衬衫的产品推荐问题。也就是说，想象我们经营一家夏威夷衬衫的在线购物网站，我们试图确定一款特定客户在购物时最大化偏好的产品。

为了简单起见，让我们假设在简要调查之后，我们得知对客户最重要的因素是衬衫上印花的数量。其他因素，如款式和颜色，也很重要，但对于这个客户来说，夏威夷衬衫最重要的是衬衫上的花朵有多少。此外，假设我们库存中有许多夏威夷衬衫，花朵数量各异，因此我们大致可以找到任何指定“花朵程度”的衬衫。因此，我们的目标是找到符合客户偏好的衬衫，这对我们来说是未知的。我们在一维搜索空间中进行这一搜索，其中空间的下限对应于没有花纹的衬衫，空间的上限包含覆盖着花朵的衬衫。

图 10.4 更详细地展示了我们的设置。在图的顶部部分，显示了客户的真实偏好以及随着衬衫花朵程度的变化而变化的情况：

x轴表示衬衫的花朵数量。在光谱的一端，我们有没有花朵的衬衫；另一端是满是花朵的衬衫。
y轴是每个衬衫的客户偏好度。客户对衬衫的偏好度越高，表示客户越喜欢这件衬衫。

图 10.4 在一个产品推荐问题中搜索具有最佳花朵数量的衬衫。我们的搜索空间是一维的，因为我们只搜索衬衫上花朵的数量。一件半面覆盖花朵的衬衫是一个局部最优点，而几乎完全覆盖的衬衫最大化了用户的偏好。

我们可以看到这个客户喜欢花纹衬衫：在衬衫的中间点过后有一个局部最优点，而偏好函数的全局最优点位于搜索空间的上界附近。这意味着一件有很多花朵但不完全覆盖的衬衫最大化了客户的偏好。

由于我们处理的是一个黑盒优化问题，在实际世界中我们实际上无法获得图 10.4 中客户的偏好曲线，我们需要使用成对比较来学习这个偏好函数，并尽快对其进行优化。现在，让我们看看如何在 Python 中设置这个优化问题。

你可能已经注意到，在图 10.4 中我们使用了前几章中使用的 Forrester 函数来模拟客户的目标函数，也就是客户的真实偏好。因此，这个函数的代码与前几章没有任何区别，其定义如下（定义范围是我们搜索空间的下界-5 和上界 5 之间）：

def objective(x):                                       ❶
    y = -((x + 1) ** 2) * torch.sin(2 * x + 2) /
    ➥5 + 1 + x / 3                                     ❶
    return y                                            ❶

lb = -5                                                 ❷
ub = 5                                                  ❷
bounds = torch.tensor([[lb], [ub]], dtype=torch.float)  ❷

❶ 目标函数

❷ 搜索空间的边界

从前几章可以记得，当我们的数据标签具有数值值时，在变量train_x中的每个数据点都有对应的train_y标签。我们当前的设置有点不同。由于我们的数据以成对比较的形式存在，每个观察结果都来自于对train_x中的两个数据点进行比较，并且观察的标签指示了客户更看重哪个数据点。

注意：我们遵循 BoTorch 的规定，用一个包含两个元素的 PyTorch 张量来编码在 train_x 中每对数据点之间的每个配对比较的结果：第一个元素是在 train_x 中被偏好的数据点的索引，第二个元素是未被偏好的数据点的索引。

举个例子，假设根据两次用户查询，我们知道用户更喜欢x = 0 而不是x = 3（也就是说，f（0）>f（3），其中f（x）是目标函数），用户也更喜欢x = 0 而不是x = -4（所以f（0）>f（-4））。我们可以用train_x来表示这两个信息作为训练数据集，train_x的取值如下:

tensor([[ 0.],   ❶
        [ 3.],   ❷
        [-4.]])  ❸

❶ 表示 x = 0

❷ 表示 x = 3

❸ 表示 x = −4

这些值是我们用来查询用户的三个x值。而训练标签train_comp则应该

tensor([[0, 1],   ❶
        [0, 2]])  ❷

❶ 表示 f(0) > f(3)

❷ 表示 f(0) > f(−4)

train_comp中的每一行都是表示成对比较结果的两个元素张量。在第一行中，[0,``1]表示train_x中索引为0的数据点（即x = 0）优先于索引为1的点（即x = 3）。同样，第二行[0,``2]编码了比较f(0) > f(-4)。

为了简化在我们的搜索空间内比较任意一对数据点的过程，我们编写了一个辅助函数，该函数接受任意两个数据点的目标值，并在第一个目标值大于第二个值时返回[0,``1]，否则返回[1,``0]：

def compare(y):
    assert y.numel() == 2                     ❶

    if y.flatten()[0] > y.flatten()[1]:       ❷
        return torch.tensor([[0, 1]]).long()  ❷
    else:                                     ❷
        return torch.tensor([[1, 0]]).long()  ❷

❶ 确保我们只有两个目标值进行比较

❷ 如果第一个值较大

❸ 如果第二个值较大

让我们使用这个函数来生成一个样本训练集。我们首先在我们的搜索空间内随机绘制两个数据点：

torch.manual_seed(0)                                              ❶
train_x = bounds[0] + (bounds[1] - bounds[0]) * torch.rand(2, 1)  ❷

❶ 为了可重现性，修正随机种子

❷ 在 0 和 1 之间绘制两个数字，并将它们缩放到我们的搜索空间

此处的变量train_x包含以下两个点：

tensor([[-0.0374],
        [ 2.6822]])

现在，我们通过评估用户的真实偏好函数并调用compare()来获得这两个点之间的比较结果：

train_y = objective(train_x)    ❶
train_comp = compare(train_y)   ❷

❶ 计算实际的目标值，这些值对我们是隐藏的

❷ 获取比较结果

将train_x中数据点的目标值的比较结果存储在train_comp中，即

tensor([[0, 1]])

这个结果意味着train_x中的第一个数据点比第二个点受到客户的更高评价。

我们还编写了另一个名为observe_and_append_data()的辅助函数，其作用是接受一对数据点，比较它们，并将比较结果添加到运行的训练集中：

该函数首先调用辅助函数compare()来获得[0,``1]或[1,``0]，然后调整存储在两个元素张量中的索引值，以便这些索引指向训练集中数据点的正确位置：

def observe_and_append_data(x_next, f, x_train, comp_train, tol=1e-3):
    x_next = x_next.to(x_train)              ❶
    y_next = f(x_next)                       ❶
    comp_next = compare(y_next)              ❶

    n = x_train.shape[-2]                    ❷
    new_x_train = x_train.clone()            ❷
    new_comp_next = comp_next.clone() + n    ❷

❶ 根据用户的偏好评估比较

❷ 跟踪索引

该函数还检查训练集中彼此接近到可以视为相同点的数据点（例如，x = 1 和 x = 1.001）。这些非常相似的数据点可能会导致我们在下一节学习的基于偏好的高斯过程的训练变得数值不稳定。我们的解决方案是标记这些相似的数据点，将它们视为重复项，并删除其中一个：

n_dups = 0

  dup_ind = torch.where(                                   ❶
      torch.all(torch.isclose(x_train, x_next[0],
      ➥atol=tol), axis=1)                                 ❶
  )[0]                                                     ❶
  if dup_ind.nelement() == 0:                              ❷
      new_x_train = torch.cat([x_train, x_next[0]
      ➥.unsqueeze(-2)])                                   ❷
  else:                                                    ❸
      new_comp_next = torch.where(                         ❸
          new_comp_next == n, dup_ind, new_comp_next - 1   ❸
      )                                                    ❸
      n_dups += 1

  dup_ind = torch.where(                                   ❹
      torch.all(torch.isclose(new_x_train, x_next[1],
      ➥atol=tol), axis=1)                                 ❹
  )[0]                                                     ❹
  if dup_ind.nelement() == 0:                              ❷
      new_x_train = torch.cat([new_x_train, x_next[1]
      ➥.unsqueeze(-2)])                                   ❷
  else:                                                    ❺
      new_comp_next = torch.where(                         ❺
          new_comp_next == n + 1 - n_dups, dup_ind,
          ➥new_comp_next                                  ❺
      )                                                    ❺

  new_comp_train = torch.cat([comp_train,
  ➥new_comp_next])                                        ❻
  return new_x_train, new_comp_train                       ❻

❶ 检查新对中第一个数据点的重复项

❷ 如果没有重复，则将数据点添加到 train_x 中

❸ 如果至少有一个重复项，则跟踪重复项的索引

❹ 检查新对中第二个数据点的重复项

❺ 如果至少有一个重复项，请跟踪重复项的索引

❻ 返回更新后的训练集

我们在训练 GP 和优化用户偏好函数的下游任务中利用这两个辅助函数，我们将在下一节中探讨第一个辅助函数。

10.3 训练基于偏好的高斯过程

我们将继续使用 CH10/01 - 从成对比较中学习.ipynb 笔记本中的代码，在本节中实现我们的 GP 模型。

我们在第 2.2.2 节中学到，在贝叶斯更新规则下（这使我们能够根据数据更新我们对数据的信念），我们可以在观察到一些变量的值的情况下获得 MVN 分布的精确后验形式。准确计算后验 MVN 分布的能力是根据新观测更新 GP 的基础。不幸的是，这种精确更新仅适用于数值观测。也就是说，我们只能使用形式为 y = f(x) 的观测精确更新 GP，其中 x 和 y 是实数。

在我们当前的设置下，观测结果以成对比较的形式出现，当以这种类型的基于偏好的数据为条件时，GP 的后验形式再也不是 GP，这排除了我们在本书中开发的大部分依赖于我们的预测模型是 GP 这一事实的方法。然而，这并不意味着我们必须放弃整个项目。

在成对比较下近似后验 GP

机器学习（以及计算机科学一般）中的一个共同主题是在无法准确完成任务时尝试近似解决任务。在我们的上下文中，这种近似等同于为我们的 GP 找到一个后验形式，该后验形式为我们观察到的成对比较提供了最高的可能性。对此方法感兴趣的读者可以在 Wei Chu 和 Zoubin Ghahramani 提出的研究论文中找到更多细节：mng.bz/9Dmo。

当然，真正最大化数据可能性的分布是非 GP 后验分布。但是由于我们希望将 GP 作为我们的预测模型，从而实现我们已学到的贝叶斯优化策略，我们的目标是找到具有最高数据可能性的 GP。请注意，找到最大化数据可能性的 GP 也是我们训练 GP 时所做的事情：我们找到最佳的 GP 超参数（例如，长度尺度和输出尺度），以最大化数据可能性。（请参见第 3.3.2 节，我们首先讨论了这种方法。）

在实现方面，我们可以使用以下代码对成对比较进行初始化和训练 GP：

BoTorch 为这个 GP 模型提供了一个特殊的类实现，命名为 PairwiseGP，可以从 botorch.models.pairwise_gp 模块中导入。
两两比较数据的可能性需要与实值数据的可能性不同的计算。对于这种计算，我们使用从同一模块导入的 PairwiseLaplaceMarginalLogLikelihood。
为了能够可视化和检查 GP 进行的预测，我们固定其输出比例，使其在训练期间保持其默认值 1。我们通过使用 model.covar_module.raw_outputscale.requires_grad_(False) 来禁用其梯度来实现这一点。这一步仅用于可视化目的，因此是可选的；在本章后面运行优化策略时我们不会这样做。
最后，我们使用 botorch.fit 中的辅助函数 fit_gpytorch_mll 来获得最大化我们训练数据可能性的后验 GP：

from botorch.models.pairwise_gp import PairwiseGP,        ❶
➥ PairwiseLaplaceMarginalLogLikelihood                   ❶
from botorch.fit import fit_gpytorch_mll                  ❶

model = PairwiseGP(train_x, train_comp)                   ❷
model.covar_module.raw_outputscale.requires_grad_(False)  ❸
mll = PairwiseLaplaceMarginalLogLikelihood(model)         ❹
fit_gpytorch_mll(mll);                                    ❺

❶ 导入必要的类和辅助函数

❷ 初始化 GP 模型

❸ 固定输出比例以获得更易读的输出（可选）

❹ 初始化（对数）可能性对象

❺ 通过最大化可能性训练模型

使用这个训练好的 GP 模型，我们现在可以在图 10.5 中跨我们的搜索空间进行预测并可视化。关于这些预测，我们注意到一些有趣的点：

均值预测遵循训练数据中表示 f(-0.0374) > f(2.6822) 的关系，在这个关系中，在 x = –0.0374 处的均值预测大于 0，而在 x = 2.6822 处小于 0。
在–0.0374 和 2.6822 处的预测不确定性也比其余预测低。这种不确定性的差异反映了观察到 f(-0.0374) > f(2.6822) 后，我们对 f(-0.0374) 和 f(2.6822) 有了一些信息，我们对这两个目标值的了解应该增加。

然而，在这些点上的不确定性并没有显著减少到零，正如我们在训练数值观测时看到的情况（例如，图 2.14 中）。这是因为，正如我们在第 10.1 节中所述，两两比较没有提供与数值评估相同数量的信息，因此仍然存在显著水平的不确定性。图 10.5 显示了我们训练的 GP 可以有效地从两两比较中学习，其中均值函数遵循观察到的比较，并且不确定性是良好校准的。

在进行预测时的 BoTorch 警告

当使用我们刚刚训练的 GP 进行预测时，您可能会遇到 BoTorch 类似以下的警告：

NumericalWarning: A not p.d., added jitter of 1.0e-06 to the diagonal
  warnings.warn(

这个警告表示 GP 生成的协方差矩阵不是正定的，导致数值稳定性相关的问题，BoTorch 已经自动向矩阵的对角线添加了“抖动”作为修复措施，所以我们用户不需要再做进一步的操作。有关我们遇到此警告的示例，请参阅第 5.3.2 节。

图 10.5 展示了由 GP 训练出的对比 f(–0.0374) > f(2.6822) 的预测。后验均值反映了这个比较的结果，而围绕两个数据点的后验标准偏差从先验中略微减小。

要进一步玩弄这个模型，并看看它如何从更复杂的数据中学习，让我们创建一个稍大一点的训练集。具体来说，假设我们想要训练 GP 在三个单独的比较上：f(0) > f(3)，f(0) > f(–4)，和 f(4) > f(–0)，所有这些对于我们在图 10.5 中有的目标函数都是正确的。为此，我们将我们的训练数据点存储在 train_x 中

train_x = torch.tensor([[0.], [3.], [-4.], [4.]])

这个集合包含了前面观察到的所有比较中涉及的所有数据点。至于 train_comp，我们使用我们在第 10.2 节中讨论过的方式，使用两元张量来编码这三个比较：

train_comp = torch.tensor(
    [
        [0, 1],    ❶
        [0, 2],    ❷
        [3, 0],    ❸
    ]
)

❶ [0, 1] 表示 f(train_x[0]) > f(train_x[1])，或者 f(0) > f(3)。

❷ [0, 2] 表示 f(train_x[0]) > f(train_x[2])，或者 f(0) > f(−4)。

❸ [3, 0] 表示 f(train_x[3]) > f(train_x[0])，或者 f(4) > f(0)。

现在，我们简单地重新声明 GP 并在这个新的训练数据上重新拟合它：

model = PairwiseGP(train_x, train_comp)              ❶
mll = PairwiseLaplaceMarginalLogLikelihood(model)    ❷
fit_gpytorch_mll(mll)                                ❸

❶ 初始化 GP 模型

❷ 初始化（对数）似然对象

❸ 通过最大化似然来训练模型

GP 模型产生了图 10.6 中显示的预测，在这里我们看到训练数据中的所有三个比较结果都反映在平均预测中，并且不确定性再次在训练数据点周围减小。

图 10.6 展示了由对比训练的 GP 进行的预测，右侧显示了后验均值反映了这个比较的结果，而在训练集中数据点周围的后验标准偏差从先验中减小了。

图 10.6 显示了我们的 GP 模型可以有效地在对比数据上进行训练。我们现在有了一种方法来从基于偏好的数据中学习，并对用户的偏好函数进行概率预测。这引导我们来到本章的最后一个话题：偏好优化中的决策制定。也就是说，我们应该如何选择数据对让用户将它们进行比较以尽快找到最受欢迎的数据点？

优先学习中目标函数的范围

在将 GP 训练成对比较的方式与使用数值评估进行比较时，一个有趣的优势是，在训练过程中不需要考虑目标函数的范围。这是因为我们只关心目标值之间的相对比较。换句话说，了解 f(x) 等同于了解 f(x) + 5，或者 2 f(x)，或者 f(x) / 10。

与此同时，当训练传统的 GP 时，使用数值评估是至关重要的，因为只有这样我们才能拥有一个具有良好校准的不确定性量化的模型。（例如，要对范围从-1 到 1 的目标函数建模，适当的输出尺度为 1，而对于范围从-10 到 10 的目标函数，我们需要更大的输出尺度。）

10.4 通过玩“山顶之王”进行偏好优化

在本节中，我们学习将 BayesOpt 应用于偏好学习。我们使用的代码包含在 CH10/02 - 优化偏好.ipynb 笔记本中。

我们需要解决的问题是如何选择最佳的一对数据点，呈现给用户，并询问他们的偏好，以找到用户最喜欢的数据点。与任何 BayesOpt 优化策略一样，我们的策略需要在利用（将搜索空间中用户价值高的区域归零）和探索（检查我们不太了解的区域）之间取得平衡。

我们在第四章到第六章学到的 BayesOpt 策略有效地使用各种启发式方法来解决利用-探索的权衡。因此，我们将开发一种策略来重新利用这些策略，以适应我们基于偏好的设置。请记住，在前几章中，BayesOpt 策略为搜索空间中的每个数据点计算一个收获分数，量化帮助我们优化目标函数的数据点的价值。通过找到最大化此收获分数的数据点，我们获得下一个要评估目标函数的点。

使用 BayesOpt 策略建议成对比较

在我们当前的基于偏好的设置中，我们需要向用户展示一对数据点以供他们比较。在优化循环的每次迭代中，我们首先组装这一对数据点，第一个是最大化给定 BayesOpt 策略的收获分数的数据点，第二个是我们迄今为止看到的最佳点。

我们使用的策略类似于流行的儿童游戏“山顶之王”，在每次迭代中，我们试图“击败”迄今为止收集到的最佳数据点（当前的“山顶之王”），使用一个由 BayesOpt 策略选择的挑战者，如图 10.7 所示。

图 10.7 在贝叶斯偏好优化中“山顶之王”策略的示意图。我们将迄今为止看到的最佳点与由 BayesOpt 策略确定的一个有希望的候选点进行比较。

通过使用这种“山顶之王”策略，我们将构造一对数据点的任务外包给了一个常规的 BayesOpt 策略，该策略能够很好地平衡利用-探索的权衡，并且我们已经知道如何使用它了。

从代码的角度来看，这个策略实现起来很简单。我们只需声明一个 BayesOpt 策略对象，并使用辅助函数optimize_acqf()优化其收获分数。例如，以下代码使用了我们在 5.2 节中学到的上置信度界（UCB）策略。UCB 策略使用由 GP 生成的预测正态分布的上界作为收获分数，以量化检查数据点的价值：

policy = UpperConfidenceBound(model, beta=2)   ❶

challenger, acq_val = optimize_acqf(           ❷
    policy,                                    ❷
    bounds=bounds,                             ❷
    q=1,                                       ❷
    num_restarts=50,                           ❷
    raw_samples=100,                           ❷
)                                              ❷

❶ 初始化 BayesOpt 策略

❷ 找到最大化收获分数的数据点

另一个我们使用的策略是期望改善（EI），我们在 4.3 节中学到了这个策略。EI 的一个特点使得它适用于我们的环境，就是该策略的动机与我们采用的“山顶之王”策略完全匹配。也就是说，EI 旨在搜索数据点，这些数据点平均而言可以从迄今为止看到的最佳点导致最大的改进（就目标函数的值而言，我们的优化目标）。超过迄今为止找到的最佳值恰恰是“山顶之王”策略的全部内容。为了在我们的环境中实现 EI，我们使用了一种不同的类实现，它可以处理嘈杂的观测值，命名为qNoisyExpectedImprovement。

BayesOpt 中的嘈杂观测

BayesOpt 中的嘈杂观测一词指的是我们怀疑观察到的标签可能受到与本章开头描述的方式相同的噪声的污染。

正如图 10.5 和 10.6 所示，在包括在我们的训练数据train_x中的位置上，我们的 GP 预测仍然存在相当大的不确定性。在这里应该使用嘈杂的 EI 版本，因为这个策略处理这种类型的不确定预测比常规 EI 策略更好。我们实现嘈杂的 EI 如下：

policy = qNoisyExpectedImprovement(model, train_x)  ❶

challenger, acq_val = optimize_acqf(                ❷
    policy,                                         ❷
    bounds=bounds,                                  ❷
    q=1,                                            ❷
    num_restarts=50,                                ❷
    raw_samples=100,                                ❷
)                                                   ❷

❶ 初始化 BayesOpt 策略

❷ 找到最大化收获分数的数据点

作为比较的一点，让我们还包括一个简单的策略，即在搜索空间内均匀随机选择挑战者来挑选到目前为止看到的最佳点：

challenger = bounds[0] + (bounds[1] - bounds[0]) * torch.rand(1, 1)   ❶

❶ 在 0 和 1 之间随机选择一个点，并将该点缩放到我们的搜索空间

这个随机策略作为一个基准，用来确定我们手头的 BayesOpt 策略是否比随机选择更好。有了这些策略，我们现在准备好运行我们的 BayesOpt 循环，以优化我们示例问题中用户的偏好。这个循环的代码类似于我们在之前章节中使用的，除了将数据点对呈现给用户以获取他们的反馈并将结果附加到我们的训练集的步骤。这是通过我们在 10.2 节中编写的observe_and_ append_data()辅助函数完成的：

incumbent_ind = train_y.argmax()                   ❶

next_x = torch.vstack([train_x[incumbent_ind,
➥:], challenger])                                 ❷

train_x, train_comp = observe_and_append_data(     ❸
    next_x, objective, train_x, train_comp         ❸
)                                                  ❸
train_y = objective(train_x)                       ❸

❶ 找到迄今为止看到的最佳点

❷ 组装最佳点和策略建议的点的批次

❸ 更新我们的训练数据

在 CH10 / 02-Optimizing preferences.ipynb 笔记本中的代码中，每个 BayesOpt 运行都始于随机生成的一对数据点，以及比较这两个点的目标函数的反馈。然后，每个运行都进行 20 个成对比较（即，向用户查询 20 个问题）。我们还为每个策略重复实验 10 次，以便观察每种策略的汇总表现。

图 10.8 显示了我们使用的优化策略找到的平均最佳值（和误差条）。EI 性能最佳，不断发现全局最优解。也许 EI 的成功很大程度上归功于我们的“国王山”方法与 EI 背后的算法动机之间的一致性。更令人惊讶的是，UCB 未能优于随机策略；也许对于权衡参数β的不同值可以改善 UCB 的性能。

图 10.8 汇总了 10 个实验的各种 BayesOpt 策略的优化性能。EI 性能最佳，不断发现全局最优解。令人惊讶的是，UCB 未能优于随机策略。

注意，UCB 的权衡参数β直接控制策略在探索和利用之间的平衡。有关此参数的更多讨论，请参见第 5.2.2 节。

在本章中，我们介绍了使用成对比较进行偏好学习和优化的问题。我们了解了数据收集的这种特定方法背后的动机以及它优于要求用户报告数字评估的优点。然后，我们使用 BayesOpt 解决了优化问题，首先使用近似方法在成对比较上训练 GP。该 GP 模型可以有效地了解在训练集中表达的数据点之间的关系，同时仍然提供良好校准的不确定性量化。最后，我们学习将 BayesOpt 策略应用于该问题，进行最佳数据点与给定 BayesOpt 策略推荐的点之间的竞争。在下一章中，我们将了解一个黑盒优化问题的多目标变体，其中我们在优化过程中需要平衡多个竞争目标函数。

摘要

在生产推荐应用中，比较两个物品可以帮助我们获得比数字评级更符合用户真实偏好的反馈。这是因为前者提出的任务量较小。
成对比较包含较少的信息，因此在选择两种引出偏好的方法时，存在减轻用户认知负担和获得信息之间的权衡。
可以训练 GP，使其最大化成对比较数据集的似然。当在成对比较数据上进行条件化时，此模型近似为真实的后验非 GP 模型。
在成对比较上训练的高斯过程产生的均值预测与训练集中的比较结果一致。特别是在首选位置的均值预测大于非首选位置的均值预测。
对于在成对比较上训练的高斯过程，其不确定性略微减小于先验高斯过程，但并未降至零，这恰如其分地反映了我们对用户偏好函数的不确定性，因为成对比较比数值评估提供的信息较少。
使用贝叶斯优化优化用户偏好的策略涉及将找到的最佳数据点与由贝叶斯优化策略推荐的候选数据点进行比较。这一策略的动机是不断尝试从迄今为止找到的最佳点进行改进。
成对比较的结果在 BoTorch 中表示为一个两元张量，其中第一个元素是首选的数据点在训练集中的索引，第二个元素是不被偏好的数据点的索引。
在使用成对比较的优化设置中使用 EI 策略时，我们使用可以更好处理训练的高斯过程中高不确定性的噪声版本的策略，而不是常规的 EI。

第十二章：同时优化多个目标

本章涵盖内容

同时优化多个目标的问题
训练多个 GP 同时学习多个目标
共同优化多个目标

每天，我们都面临着优化的权衡：

“这杯咖啡尝起来不错，但糖太多了。”
“那件衬衫看起来很棒，但超出了我的价格范围。”
“我刚训练的神经网络准确率很高，但太大了，训练时间太长。”

为了在某个目标上取得良好的性能，我们牺牲了另一个同样重要的标准：一个喜欢甜食的咖啡饮用者可能会优化咖啡的味道，同时使用不健康数量的糖；购物者在外观上打分高而在价格上打分低的服装；ML 工程师开发出具有良好预测性能但太大以至于无法在实时应用中使用的神经网络。通过专注于一个优化目标，我们可能在需要考虑的另一个目标上表现不佳。相反，我们应该将所有要优化的目标函数建模到我们的优化过程中，并尝试联合优化它们所有。例如，我们应该寻找既美味又低糖的咖啡配方，时尚又实惠的服装，或者性能良好且实用的 ML 模型。这种类型的优化问题称为多目标优化。

定义多目标优化问题，正如其名称所示，涉及到多个要同时优化的目标函数。其目标是找到在所有目标上都达到高值的数据点。

当然，在任何非平凡的多目标优化问题中，我们可能会有竞争的目标，为了在一个目标函数上获得良好的性能，唯一的方法是牺牲另一个目标上的性能。这种优化目标之间的固有冲突导致了需要平衡这些目标的需求（非常类似于需要平衡探索和开发，讨论在第 4.1.2 节中，这是 BayesOpt 循环内需要优化的两个“目标”）。

在本章中，我们学习多目标优化，如何通过找到在一个目标上表现无法改善而不牺牲另一个目标的数据点成功地解决它，以及如何在目标函数是昂贵的黑箱的情况下将 BayesOpt 应用于这个问题。多目标优化是许多领域都面临的共同问题，到本章结束时，我们将使用贝叶斯方法来解决这个问题的能力添加到我们的工具包中。

11.1 使用 BayesOpt 平衡多个优化目标

多目标优化的应用无处不在：

在工程和制造领域，工程师经常面临多个目标之间的权衡，比如产品质量与制造成本之间的权衡。例如，汽车制造商不断优化生产线以最大化质量，同时最小化成本。
在资源分配问题中，如在贫困社区或自然灾害受灾人群之间分配货币和医疗援助，决策者需要在对这些社区产生最大影响和各种物流分配困难之间取得平衡。
与我们在第 8.1.1 节中讨论的成本受限优化问题类似，开发治疗某种疾病的药物的科学家需要在最大化疗效和最小化对患者的副作用之间取得平衡。
对于机器学习工程师更相关的是，一个可以在现实世界中部署的实用机器学习模型需要在保持低训练成本的同时实现良好的性能。

与之前章节讨论的优化设置不同，我们不再有单一的优化目标可以专注。在许多这些问题中，我们需要优化的目标之间存在冲突：只有牺牲一个指标的性能我们才能在另一个指标上取得提高。思考这些优化目标之间固有的冲突的一种方式是，我们必须同时“权衡”各种目标：我们不能简单地专注于某些目标而忽略其他目标。这种同时权衡多个目标的需求在图 11.1 中可视化。幸运的是，现在有多个目标函数不会影响我们在本书中开发的大部分贝叶斯优化工作流程。

图 11.1：漫画说明我们在多目标优化中需要取得的平衡，我们需要同时权衡不同的目标函数。

使用多个高斯过程模型建模多个目标函数

在之前的章节中，我们根据观察到的数据训练了一个高斯过程模型以建模我们对单一优化目标函数的信念。在本章中，我们需要建模多个目标函数，但每个这些目标仍然可以被建模为高斯过程。通过维护这些多个高斯过程，我们有一种以概率方式推理所有目标函数的方法。

图 11.2 展示了贝叶斯优化循环，其中有两个需要优化的目标函数。与图 1.6 相比，第 1 步现在针对每个目标都有一个高斯过程，而贝叶斯优化策略识别的每个数据点在第 3 步都要对所有目标函数进行评估。

图 11.2：多目标贝叶斯优化循环，具有两个目标函数。每个目标函数的数据进行高斯过程训练，而贝叶斯优化策略决定下一步评估目标函数的数据点。

对每个目标的数据进行 GP 训练非常容易实现；事实上，我们已经在第八章对约束优化进行了这样的训练，在那里我们对目标函数进行了一次 GP 训练，对约束函数进行了另一次 GP 训练。换句话说，我们只需要专注于图 11.2 中第 2 步中贝叶斯优化策略的设计，以帮助我们在优化过程中做出有效的决策。我们重点研究了贝叶斯优化策略如何应对多个目标之间的平衡，以便在本章的其余部分尽快找到性能优异的数据点。

11.2 寻找最优数据点的边界

在本节中，我们学习了在多目标优化中常用于量化我们在优化过程中取得多大进展的数学概念。这些概念帮助我们建立本章后面我们开发的优化策略的目标。为了使我们的讨论具体化，我们使用了 CH11/01 - 计算超体积.ipynb 笔记本中的代码。

我们从需要同时优化的两个目标函数开始：

第一个目标是在以前章节中使用的熟悉的 Forrester 函数。该目标函数的全局最优解位于搜索空间的右侧。该函数在以下代码中实现为objective1()。
我们还有另一个目标函数，实现为objective2()，其具有与 Forrester 不同的功能形式和行为。关键是，该目标的全局最优解位于搜索空间的左侧——两个目标函数的全局最优解位置的不匹配模拟了多目标优化问题中常见的权衡。
我们编写一个辅助函数joint_objective()，它返回给定输入数据点x的两个目标函数的值的 PyTorch 张量。这个函数有助于保持我们的代码简洁。
最后，我们将我们的优化问题的搜索空间定义为-5 到 5 之间。

def objective1(x):                                     ❶
    return -((x + 1) ** 2) * torch.sin(2 * x + 2)
    ➥/ 5 + 1 + x / 20                                 ❶

def objective2(x):                                     ❷
    return (0.1 * objective1(x) + objective1(x - 4))
    ➥/ 3 - x / 3 + 0.5                                ❷

def joint_objective(x):                                ❸
    y1 = objective1(x)                                 ❸
    y2 = objective2(x)                                 ❸
    return torch.vstack([y1.flatten(), y2.flatten()])  ❸
    ➥.transpose(-1, -2)                               ❸

lb = -5                                                ❹
ub = 5                                                 ❹
bounds = torch.tensor([[lb], [ub]], dtype=torch.float) ❹

❶ 第一个目标函数

❷ 第二个目标函数

❸ 调用两个目标函数的辅助函数

❹ 搜索空间的边界

图 11.3 显示了我们搜索空间中的这两个目标函数。我们看到最大化这两个目标的数据点彼此不同：实线曲线在x=4.5 附近最大化，而虚线曲线在x=-4.5 附近最大化。这种差异意味着我们有两个相互冲突的目标，并且这两个函数的联合优化需要权衡它们的目标值。

图 11.3 我们当前的多目标优化问题的两个目标函数。最大化这两个目标的数据点彼此不同，因此在优化这两个目标时存在权衡。

通过“权衡”，我们指的是存在搜索空间中的点x，其第一个目标的值（表示为f1）无法提高，除非第二个目标的值（表示为f2）降低。换句话说，存在一些数据点在优化一个目标函数方面，其值无法超过，除非我们牺牲另一个目标函数。

例如，考虑图 11.3 中标示的x = –5。这是搜索空间的最左边的数据点。这个点的目标值f1 大约为 4，目标值f2 大约为 1.5。现在，x = –5 是一个数据点，如果我们想在第一个目标f1 上做得比 4 更好，我们将不得不在第二个目标f2 上做得不如 1.5。事实上，我们要实现高于 4 的f1 值的唯一方式是查询空间的最右侧，其中x > 4。在这里，f2 的值下降到 0 以下。

相反，右侧的区域（x > 4）也是f1 和f2 之间的紧张存在的地方：要增加f2 的值，我们必须向空间左侧移动，这样f1 的值就会受到影响。

定义：数据点其一个目标的值无法超过，除非另一个目标的值降低，称为非支配。相反，支配点x[1]是指存在另一个点x[2]，其所有目标值都超过x[1]的点。非支配点也可以称为帕累托最优、帕累托有效或非劣。

因此，点x = –5 是一个非支配点，一些点当x > 4 也是非支配点。图 11.3 中的一个支配点的例子是x = –1.9，它给出了f1 ≈ f2 ≈ 1。这个点被x = –5 支配，因为前者的目标值低于后者：f1 < f1 和f2 < f2。

在许多情况下，我们有无限多个非支配点。图 11.4 展示了我们当前问题中的非支配点作为虚线阴影区域（我们将在本节稍后讨论如何找到这些非支配点；现在，让我们关注这些非支配点的行为）：

我们发现x = –5 确实是一个非支配点，以及该区域周围许多点给出了较高的第二个目标f2 的值。超出此区域的点不会产生更高的f2 值，因此该区域内的点是非支配的。我们将这些点称为group 1。
右侧的小区域为第一个目标f1 提供了较高的值，同样是非支配的。这些点称为group 2。
在x=4 周围还有一个第三个最小区域，它也是非支配的，在搜索空间的左侧没有被非支配点的f1 值超过。尽管这个区域不包含任何目标函数的全局最优解，但该区域在两个目标的值之间进行权衡，因此是非支配的。我们称这些点为“群组 3”。

图 11.4 两个目标函数和非支配点。在这个多目标优化问题中有无穷多个非支配点。

非支配点在多目标优化中非常有价值，因为它们本身就是优化问题的解，除非牺牲至少一个目标，否则我们无法改进它们。通过研究非支配点，它们之间的关系以及它们在搜索空间中的分布，我们可以更加了解优化问题中多个目标之间的权衡。因此，多目标优化的一个合理目标是找到尽可能多的非支配点。

然而，我们并不立即清楚如何具体量化找到非支配点的目标。我们不应该简单地寻求揭示尽可能多的非支配点，因为它们可能是无限多的。相反，我们使用一个更容易思考的量，如果我们通过将数据点可视化到一个不同的空间来观察它。

在图 11.3 和 11.4 中，x-轴对应于数据点本身，y-轴对应于这些数据点的目标值。为了研究两个冲突目标之间的权衡，我们还可以使用散点图，其中给定数据点x的x-坐标是第一个目标函数f1 的值，y-坐标是第二个目标函数f2 的值。

图 11.5 展示了在-5 和 5 之间的 201 个等距点的密集网格中每个点的散点图，被支配的点用点表示，非支配点用星号表示。我们可以看出，在这个空间中，一个点是否被支配更容易确定：对于每个数据点x[1]，如果存在另一个数据点x[2]，它位于x[1]的上方且右侧，则x[1]是一个被支配的点；相反，如果不存在同时位于x[1]的上方且右侧的点x[2]，则x[1]是非支配的。在图 11.5 中我们还可以看到三个非支配点的组，与图 11.4 相关的讨论相符。

图 11.5 基于两个目标函数值的数据点的散点图。被支配的点用点表示，非支配点用星号表示。非支配点的三个组对应于图 11.4 中的讨论。

从可视化的目标值空间中的非支配点集合中，我们现在引入另一个概念：Pareto 前沿。图 11.6 可视化了当前优化问题的 Pareto 前沿。

图片 11.6

图 11.6 通过非支配点绘制的 Pareto 前沿。没有数据点位于此 Pareto 前沿的右侧（上方和右侧）。

定义跟踪非支配点的曲线称为 Pareto 前沿。它被称为前沿，因为当我们将所有数据点视为一个集合时，非支配点的这条曲线构成了集合的边界或前沿，在该边界或前沿之外没有数据点。

Pareto 前沿的概念在多目标优化中至关重要，因为前沿直接导致可以量化多目标优化问题中的进展的度量标准。特别是，我们关注 Pareto 前沿覆盖的空间——由多个目标的收集目标值定义的空间——的大小；也就是说，Pareto 前沿内部（下方和左侧）的区域。该区域显示为图 11.7 的左侧面板中的阴影区域。

定义我们使用术语 支配超体积（有时简称为 超体积）来表示 Pareto 前沿覆盖了多少空间。在我们的示例中有两个目标函数，因此空间是二维的，支配超体积是支配区域的面积。当有两个以上的目标时，支配超体积以更高维度度量相同的数量。

左侧图 11.7 中显示的分散点是使用密集的网格在搜索空间中生成的，以便我们可以详细研究 Pareto 前沿及其超体积的行为。换句话说，这个密集的网格代表了对空间的穷尽搜索，以完全绘制出 Pareto 前沿。

作为比较的一点，图 11.7 的右侧面板显示了在 -5 到 5 之间均匀选择的 20 个点的结果。从所选点中，我们再次找到在 20 个点集内没有被任何其他点支配的点，并绘制此第二个数据集的 Pareto 前沿。

图片 11.7

图 11.7 密集网格的支配超体积（左侧），等效于穷举搜索，以及随机选择的 20 个数据点（右侧）。第一个数据集具有更大的支配体积，因此在多目标优化方面的表现比第二个数据集好。

与左侧情况不同，左侧完全覆盖搜索空间，而右侧的小数据集只有四个非支配点。在多目标优化的背景下，我们使用第一个数据集（来自穷举搜索）比使用第二个数据集（来自随机搜索）取得更多进展。

与穷尽搜索得到的数据集相比，这四个非支配点构成了一个更加锯齿状的帕累托前沿，进而拥有一个更小的被支配超体积。换句话说，这个被支配超体积的度量可以用来量化多目标优化问题中的优化进展。

注意在多目标优化问题中，我们通过当前收集的数据产生的被支配区域的超体积来衡量优化进展。我们收集的数据的被支配超体积越大，我们同时优化目标函数的进展就越大。

根据超体积度量，图 11.7 显示，穷尽搜索比随机搜索（查询更少）在优化方面做得更好，这是一个预期结果。但是要量化前者搜索策略比后者好多少，我们需要一种计算这个超体积度量的方法。对于这个计算，需要一个参考点；这个参考点充当被支配区域的终点，为该区域设置了一个左下边界。我们可以将这个参考点视为在多目标优化设置下我们能观察到的最差结果，因此该区域与帕累托前沿之间的超体积量化了我们从这个最差结果改进了多少。（如果我们不知道每个目标函数的最差结果，BayesOpt 用户可以将每个查询可以达到的最低值设定为我们认为的最低值。）

注意在多目标优化中，一个常见的参考点是一个数组，其中的每个元素对应于要最大化的目标函数的最低值。

例如，我们当前优化问题的参考点为[–2.0292, –0.4444]，因为第一个元素–2.0292 是第一个目标函数的最小值（图 11.3 中的实线曲线），–0.4444 是第二个目标的最小值（图 11.3 中的虚线曲线）。这个参考点在图 11.8 中被可视化为星号，再次为被支配空间设定了一个下限。

图 11.8 在多目标优化问题中的参考点，它为被支配空间设定了一个下限。超体积计算为参考点与帕累托前沿之间区域的体积。

有了这个参考点，我们可以计算由多目标优化策略收集的数据集占优区域的超体积。完成这个计算的算法涉及将占优区域分成多个不相交的超矩形，这些超矩形共同构成了占优区域。从那里，我们可以容易地计算每个超矩形的超体积，并将它们相加以获得整个区域的超体积。有兴趣的读者可以参考 Renaud Lacour、Kathrin Klamroth 和 Carlos M. Fonseca 等人提出的此算法的研究论文（mng.bz/jPdp）。

使用 BoTorch，我们可以导入和运行此算法，而无需实现底层细节。具体而言，假设我们已将优化期间找到的收集标签存储在变量train_y中。因为我们的示例中有两个目标函数，所以train_y的形状应为n-by-2，其中n是收集集合中数据点的数量。然后，我们可以使用下面的代码来计算超体积度量，其中

DominatedPartitioning类实现了占优区域的分区。为了初始化此对象，我们传入参考点和收集标签train_y。
然后我们调用占优区域对象的compute_hypervolume()方法来计算其超体积度量：

from botorch.utils.multi_objective
➥.box_decompositions.dominated import                ❶
➥DominatedPartitioning                               ❶

dominated_part = DominatedPartitioning
➥(ref_point, train_y)                                ❷
volume = dominated_part.compute_hypervolume().item()  ❷

❶ 导入占优区域类的实现

❷ 计算相对于参考点的占优区域的超体积度量

使用此方法，我们可以计算完全搜索和随机搜索的超体积度量，如图 11.9 左侧和中间面板所示。我们看到，与随机搜索的超体积度量（25.72）相比，完全搜索确实实现了更高的超体积度量（31.49）。

图 11.9 各种搜索策略的多目标优化结果及相应的超体积度量。BayesOpt 几乎达到了完全搜索的超体积度量，但查询少了很多。

在图 11.9 的右侧面板中，我们还可以看到我们在下一节学习的 BayesOpt 策略仅使用了 20 个数据点就实现了相应的结果。只有预算的十分之一（20 与 201），BayesOpt 就几乎达到了完全搜索的超体积度量。与具有相同预算的随机搜索相比，BayesOpt 能够更全面地映射出真正的 Pareto 前沿，并实现更高的超体积度量。

11.3 寻求改进最佳数据边界

贝叶斯优化策略应如何最大化其收集数据中受支配区域的超体积？一个简单的策略是在迭代方式下交替优化每个目标：在贝叶斯优化循环的这一次迭代中，我们试图最大化第一个目标 f1；在下一次迭代中，我们则试图最大化第二个目标 f2；依此类推。在一个迭代中，我们有一个特定的目标要优化，我们可以通过使用我们在第 4 至 6 章学到的各种贝叶斯优化策略来实现这一点。在本章的其余部分中，我们使用了期望改进（EI），这是我们在第 4.3 节中学到的。EI 是一种在实践中常用的策略，因为它的算法简单且性能稳定。

假设在我们的多目标优化问题中，我们观察到图 11.10 顶部面板中 X 标记的数据点。通过对属于每个目标函数的数据集进行 GP 训练，我们得到了第一个目标的 GP 预测（左上角面板）和第二个目标的 GP 预测（右上角面板）。

在图 11.10 的底部面板中，我们展示了各个 EI 策略在相应目标函数上的收购分数。左下角的 EI 试图最大化第一个目标 f1，而右下角的 EI 则寻找第二个目标 f2 的最优解。我们可以看到，当第一个 EI 关注搜索空间的右侧区域，即 f1 最大化的区域时，而第二个 EI 则关注左侧区域，即 f2 最大化的区域时，两个目标之间的冲突在这里是明显的。

图 11.10：关于两个目标函数的当前 GP 信念（顶部）和相应的 EI 收购分数（底部）。每个 EI 策略都试图优化自己的目标函数，并关注不同的区域。

注意：由贝叶斯优化策略计算的数据点的收购分数 quantifies 了数据点对我们搜索目标函数最优解的价值。收购分数越高，数据点越有价值，而给出最高收购分数的点是策略建议查询的点。

在我们之前提出的交替策略中，我们要么遵循第一个 EI 策略并查询 x = 4.5 附近的点，要么遵循第二个 EI 并查询 x = –4.5 附近的点，这取决于是 f1 还是 f2 被优化的轮次。我们将这种交替策略作为基准，用来与我们最终的解决方案进行比较。

为了比简单交替不同目标函数的策略做得更好，这个解决方案应该是什么？我们注意到，通过让每个要最大化的目标都有一个 GP 模型，我们可以概率地推理出每个潜在新查询在每个目标上的值同时。具体来说，我们知道每个潜在新查询在每个目标上的价值都遵循一个已知的正态分布；这个正态分布是我们对查询价值的预测。

这个预测让我们能够推理出每个潜在新查询是否是一个无支配点，如果是的话，它将如何增加被支配区域的超体积。每个新观察到的无支配点都会延伸被支配区域的边界（即 Pareto 边界），因此增加了支配超体积。因此，我们可以使用每个新查询导致的超体积增加的期望值作为收购分数，来量化查询的价值。我们能够期望从查询获得的超体积增加越大，它对我们进行优化的帮助就越大。

当然，我们无法确定通过查询会获得多少超体积的增加，直到我们真正对目标函数进行查询。然而，我们可以通过一种概率方式来推理这个超体积增加。也就是说，我们可以计算可能查询产生的超体积增加的期望值。

同样，类似于确定被支配区域的超体积的算法，这个对超体积增加的期望计算涉及将被支配区域分成超矩形，非常复杂。再一次地，我们这里不会详细介绍数学细节，但是你可以参考杨凯锋、米歇尔·埃默里奇、安德烈·德茨和托马斯·贝克提出的相应 BayesOpt 策略的研究论文，该论文称为预期超体积增加（EHVI）以获取更多详情（mng.bz/WzYw）。

定义预期超体积增加策略使用新数据点导致的被支配区域超体积增加的期望值作为该数据点的收购分数。这个策略是将 EI 推广到了多目标设置的结果，其中我们的目标是最大化被支配超体积。

图 11.11 显示了与图 11.10 中相同数据集的 EHVI 的收购分数在底部右侧面板中。我们可以看到，与单个 EI 策略相比，EHVI 通过将高收购分数分配给可能扩展 Pareto 边界的多个区域来很好地平衡了两个目标：搜索空间的最左侧区域的得分最高，但最右侧区域以及中间的其他区域也有非常可观的收购分数。

图 11.11 当前 GP 对每个目标函数的信念（顶部）、对应的 EI 获取分数（左下角）和 EHVI 获取分数（右下角）的看法。EHVI 平衡了这两个目标，将高获取分数分配给可能延伸 Pareto 边界的多个区域。

为了验证这个 EHVI 策略确实在多目标优化中给我们带来了优势，我们实现了这个策略并在当前问题上运行它。我们使用的代码包含在 CH11/02 - Multi-objective BayesOpt loop.ipynb 笔记本中。

首先，我们需要 GP 模型的类实现和一个帮助函数 fit_gp_model()，它有助于在观察到的数据上训练每个 GP。由于我们在前几章中已经实现了这些组件，所以我们不会再次在这里展示它们的代码；您可以参考第 4.1.1 节来复习这些代码。在 BayesOpt 循环的每一步中，我们调用帮助函数在每个目标函数的数据上初始化和训练一个 GP。在我们的情况下，我们有两个目标函数，所以我们分别调用帮助函数两次，每次分别使用 train_y[:, 0]（从第一个目标 f1 观察到的标签）或 train_y[:, 1]（从第二个目标 f2 观察到的标签）。

model1, likelihood1 = fit_gp_model(train_x, train_y[:, 0])
model2, likelihood2 = fit_gp_model(train_x, train_y[:, 1])

然后我们使用 botorch.acquisition.multi_objective.analytic 模块中的 ExpectedHypervolumeImprovement 类来实现 EHVI 策略。为了初始化策略对象，我们设置以下参数：

参数 model 接受一系列的 GPs，每个 GP 建模一个目标函数。这个 GP 列表被实现为 ModelListGP 类的一个实例，接受单独的 GP 对象 (model1, model2)。
参数 ref_point 接受参考点，这对于计算 HV 和潜在 HV 增加量是必要的。
最后，参数 partitioning 接受 FastNondominatedPartitioning 类的一个实例，它有助于计算 HV 增加量。这个对象的初始化与我们之前看到的 DominatedPartitioning 对象类似，接受一个参考点和观察标签 train_y：

from botorch.acquisition.multi_objective
➥.analytic import                                    ❶
➥ExpectedHypervolumeImprovement                      ❶
from botorch.utils.multi_objective.box_decompositions
➥.non_dominated import                               ❶
➥FastNondominatedPartitioning                        ❶
from botorch.models.model_list_gp_regression import ModelListGP

policy = ExpectedHypervolumeImprovement(
    model=ModelListGP(model1, model2),                ❷
    ref_point=ref_point,                              ❸
    partitioning=FastNondominatedPartitioning
    ➥(ref_point, train_y)                            ❹
)

❶ 导入必要的类

❷ GP 模型列表，每个模型对应一个目标函数

❸ 参考点

❹ 无支配分区对象用于计算 HV 增加量

使用 EHVI 策略的 policy 对象，我们可以计算获取分数，表示由潜在新观测引起的预期 HV 增加量。然后我们可以使用辅助函数 optimize_acqf() 找到给出最高分数的数据点：

next_x, acq_val = optimize_acqf(
    policy,
    bounds=bounds,
    q=1,
    num_restarts=20,
    raw_samples=50
)

变量 next_x 存储我们将在下一步中使用的查询位置：next_y = joint_objective(next_x)。

这就是我们在当前优化问题上运行 EHVI 所需的一切。作为参考，我们还测试了之前讨论过的交替优化策略，在这种策略中，我们使用常规的 EI 来优化选择的目标函数。由于我们有两个目标，我们只需在两个目标之间来回切换（这里的num_queries是贝叶斯优化运行中可以进行的总评估次数）：

for i in range(num_queries):
    if i % 2 == 0:                        ❶
        model = model1                    ❶
        best_f = train_y[:, 0].max()      ❶
    else:                                 ❷
        model = model2                    ❷
        best_f = train_y[:, 1].max()      ❷

    policy = ExpectedImprovement(model=model,
    ➥best_f=best_f)                      ❸

❶ 如果当前迭代次数为偶数，则优化第一个目标

❷ 如果当前迭代次数为奇数，则优化第二个目标

❸ 相应地创建 EI 策略

最后，为了量化我们的优化进展，我们记录了由当前数据集在搜索过程中收集的支配区域的超体积。这个记录是用一个名为hypervolumes的张量完成的，它在实验过程中在每一步存储当前的支配超体积，跨多个实验。总的来说，我们的贝叶斯优化循环如下，对于每个策略，我们运行实验多次，每次都使用均匀随机选择的初始数据集：

hypervolumes = torch.zeros((num_repeats, num_queries))   ❶

for trial in range(num_repeats):
  torch.manual_seed(trial)                               ❷
  train_x = bounds[0] + (bounds[1] - bounds[0]) * torch  ❷
  ➥.rand(1, 1)                                          ❷
  train_y = joint_objective(train_x)                     ❷
  for i in range(num_queries):
    dominated_part = DominatedPartitioning(ref_point,
    ➥train_y)                                           ❸
    hypervolumes[trial, i] = dominated_part
    ➥.compute_hypervolume().item()                      ❸

    ...                                                  ❹

❶ 优化过程中发现的超体积历史

❷ 初始化一个随机初始训练集

❸ 记录当前的超体积

❹ 重新训练模型，初始化一个策略，并找到下一个查询

CH11/02 - 多目标贝叶斯优化循环.ipynb 笔记本对我们有两个 20 个查询的实验，每个实验都有 10 次贝叶斯优化策略进行运行。图 11.12 显示了两个策略所进行的查询次数与平均超体积和误差棒的关系。我们看到 EHVI 一直优于交替 EI 策略，这说明了基于超体积的方法的好处。

图 11.12 两个贝叶斯优化策略所进行的查询次数与平均超体积和误差棒的关系。EHVI 一直优于交替 EI 策略。

在本章中，我们学习了多目标优化问题以及如何使用贝叶斯优化方法来解决它。我们讨论了超体积的概念作为优化性能的衡量标准，量化了我们在优化目标函数方面取得的进展。通过使用 EI 策略的变体来优化超体积的增加，我们得到了一个表现强劲的 EHVI 策略。

不幸的是，本章无法涵盖多目标贝叶斯优化的其他方面。具体来说，除了 EHVI 之外，我们还可以考虑其他优化策略。一种常见的技术是标量化，它通过取加权和将多个竞争目标合并为一个目标。该策略是交替 EI 策略的一般化，我们可以认为在每次迭代中将一个目标的权重设置为 1，另一个目标的权重设置为 0。感兴趣的读者可以参考 BoTorch 文档（请参阅botorch.org/docs/multi_objective和botorch.org/tutorials/multi_objective_bo），该文档提供了 BoTorch 提供的不同多目标优化策略的简要摘要。

11.4 练习：飞机设计的多目标优化

在这个练习中，我们将所学的多目标优化技术应用于优化飞机的航空结构设计问题。这个问题首次在第七章的练习 2 中介绍，并在第八章的练习 2 中修改为成本约束问题。我们在这里重用第八章的代码。这个练习使我们能够观察到期望超体积改进（EHVI）策略在多维问题中的性能。解决方案包含在 CH11/03 - Exercise 1.ipynb 笔记本中。

执行以下步骤：

从第八章的练习 2 中复制目标函数flight_utility()和flight_cost()的代码。取反第二个函数flight_cost()返回值的符号。我们将这两个函数用作多目标优化问题的目标。
编写一个辅助函数，该函数接受一个输入X（可能包含多个数据点），并返回在两个目标函数上评估的X的值。返回的值应该是一个大小为n-by-2 的张量，其中n是X中数据点的数量。
声明搜索空间为四维单位正方形。即，四个下限为 0，四个上限为 1。
要计算由优化算法收集的数据集的超体积，我们需要一个参考点。声明此参考点为[–1.5, –2]，这是两个目标函数的对应最低值。
实现 GP 模型的类，该类应具有常数均值和一个四维 Matérn 2.5 核，并具有自动相关性确定（ARD；参见 3.4.2 节），以及一个辅助函数fit_gp_model()，该函数在训练集上初始化和训练 GP 模型。有关实现这些组件的详细信息，请参见 4.1.1 节。
将要运行的实验次数设置为 10，每次实验的预算（要进行的查询数量）设置为 50。
运行 EHVI 策略来优化我们拥有的两个目标函数，以及第 11.3 节中讨论的交替 EI 策略。绘制这两个策略所实现的平均超体积和误差条形图（类似于图 11.2），并比较它们的性能。

总结

当存在多个潜在的冲突目标需要同时优化时，多目标优化问题就会出现。这个问题在现实世界中很常见，因为我们经常在许多真实任务中与多个竞争目标相争。
在使用 BayesOpt 进行多目标优化时，我们使用多个高斯过程来模拟我们对目标函数的信念（每个目标函数使用一个模型）。我们可以使用这些高斯过程以概率方式同时推理目标函数。
非支配点实现的目标值不能得到改善，除非我们牺牲至少一个目标函数的性能。发现非支配数据点是多目标优化的目标，因为它们允许我们研究目标函数之间的权衡。
无支配数据点构成帕累托前沿，它设置了代表多目标优化中最优的边界。没有数据点位于所有非支配点的帕累托前沿之外。
支配空间的超体积（即由帕累托前沿覆盖的区域）测量算法收集的数据集的优化性能。超体积越大，算法的性能就越好。可以通过在 BoTorch 的"DominatedPartitioning"类的实例上调用compute_hypervolume()方法来计算数据集的超体积。
要计算数据集的超体积，我们需要一个作为支配空间结束点的参考点。我们通常将参考点设置为要优化的目标函数的最低值。
由于高斯过程允许我们对目标函数进行预测，因此我们可以寻求改进当前数据集的超体积。这种策略对应于 EHVI 策略，是多目标优化中 EI 的一种变体。这种策略成功平衡了竞争性目标。

第四部分：特殊高斯过程模型

高斯过程（GPs），在 BayesOpt 的背景之外，本身就是一类强大的 ML 模型。虽然本书的主要主题是 BayesOpt，但如果不多关注 GPs，那就是一个错失的机会。这一部分向我们展示了如何扩展 GPs，并使它们在各种 ML 任务中更实用，同时保留了它们最有价值的特性：对预测中不确定性的量化。

第十二章中，我们学习如何加速训练高斯过程（GPs）并将其扩展到大数据集。这一章帮助我们解决了 GPs 的一个最大劣势：它们的训练成本。

第十三章展示了如何通过将 GPs 与神经网络结合，将 GP 的灵活性提升到另一个水平。这种组合提供了两全其美的好处：神经网络近似任何函数的能力和 GPs 对不确定性的量化。这一章还让我们真正欣赏到了在 PyTorch、GPyTorch 和 BoTorch 中拥有一个简化的软件生态系统的好处，这使得同时使用神经网络和 GPs 变得无缝。

第十三章：将高斯过程扩展到大型数据集

本章介绍：

训练大型数据集上的 GP
在训练 GP 时使用小批量梯度下降。
采用高级梯度下降技术来更快地训练 GP。

到目前为止，我们已经看到 GP 提供了极高的建模灵活性。在第三章中，我们学习了如何使用 GP 的均值函数来模拟高级别趋势，并使用协方差函数来模拟变异性。GP 还提供了校准的不确定性量化。也就是说，训练数据集中接近观测值的数据点的预测比远离观测值的点的预测具有更低的不确定性。这种灵活性使 GP 与其他只产生点估计（如神经网络）的 ML 模型区别开来。然而，它也导致了速度问题。

训练和预测 GP（具体来说，计算协方差矩阵的逆）与训练数据的规模呈立方级扩展关系。也就是说，如果我们的数据集大小翻倍，GP 将需要花费八倍的时间进行训练和预测。如果数据集增加十倍，GP 将需要花费 1,000 倍的时间。这给将 GP 扩展到大型数据集的应用带来了挑战。

如果我们的目标是对整个国家（例如美国）的房价进行建模，其中每个数据点表示给定时间的单个住宅的价格，则我们的数据集大小将包含数亿个点。例如，在线数据库 Statista 记录了自 1975 年至 2021 年美国住房单位数量的变化；该报告可在 www.statista.com/statistics/240267/number-of-housing-units-in-the-united-states/ 上访问。我们可以看到，自 1975 年以来，这个数字一直在稳步增长，1990 年超过了 1 亿，现在已经超过 1.4 亿。
在我们在第 1.1.3 节中讨论的药物发现应用程序中，一个可能被合成为药物的可能分子的数据库可能会拥有数十亿个条目。
在天气预报中，低成本的监测设备使得大规模收集天气数据变得容易。数据集可以包含跨多年的每分钟测量结果。

考虑到正常 GP 模型的立方运行时间，将其应用于这种规模的数据集是不可行的。在本章中，我们将学习如何利用一类称为“变分高斯过程”（VGPs）的 GP 模型来解决从大型数据中学习的问题。

定义变分高斯过程选择一小部分数据，很好地表示整个集合。它通过寻求最小化它本身和完整数据集上训练的普通 GP 之间的差异来实现这一点。术语“变分”是指研究函数式优化的数学子领域。

选择仅对这些代表性点的小型子集进行训练的想法是相当自然和直观的。图 12.1 展示了 VGP 的运行情况，通过从少数精选数据点中学习，该模型产生了几乎与普通 GP 产生的预测相同的预测。

图 12.1 显示了普通 GP 和 VGP 的预测。VGP 产生了几乎与 GP 相同的预测，同时训练时间显著缩短。

我们在本章介绍如何实现这个模型，并观察其在计算上的优势。此外，当使用 VGP 时，我们可以使用更高级的梯度下降版本，正如我们在第 3.3.2 节中看到的，它用于优化 GP 的超参数。我们学会使用这个算法的版本来更快、更有效地训练，并最终将我们的 GPs 扩展到大型数据集。本章附带的代码可以在 CH11/01 - 近似高斯过程推理.ipynb 笔记本中找到。

12.1 在大型数据集上训练高斯过程模型

在本节中，为了直观地看到在大型数据集上训练高斯过程模型面临的困难挑战，我们试图将我们在第二章和第三章中使用的 GP 模型应用于一个包含 1,000 个点的中型数据集。这个任务将清楚地表明使用普通 GP 是不可行的，并激发我们在下一节学习的内容：变分 GPs。

12.1.1 设置学习任务

在这个小节中，我们首先创建我们的数据集。我们重新使用了在第二章和第三章中看到的一维目标函数，即 Forrester 函数。我们再次按照以下方式实现它：

def forrester_1d(x):
    y = -((x + 1) ** 2) * torch.sin(2 * x + 2) / 5 + 1
    return y.squeeze(-1)

类似于我们在第 3.3 节中所做的，我们还将拥有一个辅助函数，该函数接收一个 GP 模型，并在整个域上可视化其预测。该函数具有以下头部，并接受三个参数——GP 模型、相应的似然函数和一个布尔标志，表示模型是否为 VGP：

此辅助函数的逻辑概述在图 12.2 中草绘出来，它包括四个主要步骤：计算预测、绘制真实函数和训练数据、绘制预测，最后，如果模型是 VGP，则绘制诱导点。

图 12.2 是可视化 GP 预测的辅助函数的流程图。该函数还显示了 VGP 的诱导点（如果传入的是该模型）。

定义：诱导点是 VGP 模型选择的表示整个数据集的小型子集，用于训练。顾名思义，这些点旨在“诱导”关于所有数据的知识。

现在我们更详细地介绍这些步骤。在第一步中，我们使用 GP 计算均值和 CI 预测：

with torch.no_grad():
    predictive_distribution = likelihood(model(xs))
    predictive_mean = predictive_distribution.mean
    predictive_upper, predictive_lower =
    ➥predictive_distribution.confidence_region()

在第二步中，我们制作 Matplotlib 图并显示存储在 xs 和 ys 中的真实函数（我们稍后生成）以及我们的训练数据 train_x 和 train_y：

plt.figure(figsize=(8, 6))

plt.plot(xs, ys, label="objective", c="r")  ❶
plt.scatter(                                ❷
    train_x,                                ❷
    train_y,                                ❷
    marker="x",                             ❷
    c="k",                                  ❷
    alpha=0.1 if variational else 1,        ❷
    label="observations",                   ❷
  )                                         ❷

❶ 绘制真实目标函数

❷ 为训练数据制作散点图

在这里，如果模型是一个 VGP（如果 variational 设置为 True），那么我们会用较低的不透明度绘制训练数据（通过设置 alpha = 0.1），使它们看起来更透明。这样我们可以更清楚地绘制后面学习到的 VGP 的代表性点。

GP 所做的预测随后在第三步中以实线均值线和阴影 95% CI 区域的形式显示：

plt.plot(xs, predictive_mean, label="mean")
plt.fill_between(
    xs.flatten(),
    predictive_upper,
    predictive_lower,
    alpha=0.3,
    label="95% CI"
)

最后，我们通过提取 model.variational_strategy.inducing_points 来绘制 VGP 选择的代表性点：

if variational:
  inducing_points =
  ➥model.variational_strategy.inducing_points.detach().clone()
  with torch.no_grad():
      inducing_mean = model(inducing_points).mean

  plt.scatter(
      inducing_points.squeeze(-1),
      inducing_mean,
      marker="D",
      c="orange",
      s=100,
      label="inducing pts"
  )                         ❶

❶ 散布感应点

现在，为了生成我们的训练和数据集，我们在-5 和 5 之间随机选择了 1,000 个点，并计算了这些点的函数值：

torch.manual_seed(0)
train_x = torch.rand(size=(1000, 1)) * 10 - 5
train_y = forrester_1d(train_x)

为了生成我们的测试集，我们使用 torch.linspace() 函数在-7.5 和 7.5 之间计算了一个密集的网格。该测试集包括-7.5、7.4、-7.3 等，直到 7.5：

xs = torch.linspace(-7.5, 7.5, 151).unsqueeze(1)
ys = forrester_1d(xs)

要可视化我们的训练集的样子，我们可以再次制作一个散点图如下：

plt.figure(figsize=(8, 6))
plt.scatter(
    train_x,
    train_y,
    c="k",
    marker="x",
    s=10,
    label="observations"
)
plt.legend();

这段代码产生了图 12.3，其中黑点表示我们训练集中的个别数据点。

图 12.3 我们学习任务的训练数据集，包含 1,000 个数据点。在这个集合上训练一个常规的 GP 需要相当长的时间。

12.1.2 训练常规 GP

我们现在准备在这个数据集上实现并训练一个 GP 模型。首先，我们实现 GP 模型类，其具有常数函数（gpytorch.means .ConstantMean 的一个实例）作为其均值函数，以及具有输出尺度的 RBF 核函数（使用 gpytorch.kernels.ScaleKernel(gpytorch.kernels.RBFKernel()) 实现）作为其协方差函数：

class GPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super().__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)                               ❸

❶ 常数均值函数

❷ 具有输出尺度的 RBF 核函数

❸ 创建 MVN 分布作为预测

现在，我们使用我们的训练数据和一个 GaussianLikelihood 对象初始化了这个 GP 模型：

likelihood = gpytorch.likelihoods.GaussianLikelihood()
model = GPModel(train_x, train_y, likelihood)

最后，我们通过运行梯度下降来训练我们的 GP，以最小化由数据的可能性定义的损失函数。在训练结束时，我们得到了模型的超参数（例如，均值常数、长度尺度和输出尺度），这些超参数给出了一个较低的损失值。梯度下降是使用 Adam 优化器（torch .optim.Adam）实现的，这是最常用的梯度下降算法之一：

optimizer = torch.optim.Adam(model.parameters(), lr=0.01)          ❶
mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)  ❷

model.train()                                                      ❸
likelihood.train()                                                 ❸

for i in tqdm(range(500)):                                         ❹
    optimizer.zero_grad()                                          ❹

    output = model(train_x)                                        ❹
    loss = -mll(output, train_y)                                   ❹

    loss.backward()                                                ❹
    optimizer.step()                                               ❹

model.eval()                                                       ❺
likelihood.eval()                                                  ❺

❶ 梯度下降算法 Adam

❷ 损失函数，计算由超参数确定的数据的可能性

❸ 启用训练模式

❹ 运行 500 次梯度下降迭代

❺ 启用预测模式

注意：作为提醒，当训练 GP 时，我们需要为模型和可能性都启用训练模式（使用 model.train() 和 likelihood .train()）。在训练后和进行预测之前，我们需要启用预测模式（使用 model.eval() 和 likelihood.eval()）。

使用 GPU 训练 GPs

一个将 GP 扩展到大型数据集的方法，本章不重点讨论，是使用图形处理单元（GPU）。GPU 通常用于并行化矩阵乘法并加速训练神经网络。

同样的原则也适用于此处，而 GPyTorch 通过遵循 PyTorch 的语法将 GP 训练简化到了 GPU 上（通过在对象上调用 cuda() 方法）。具体来说，我们调用 train_x = train_x.cuda() 和 train_y = train_y.cuda() 将我们的数据放到 GPU 上，然后调用 model = model.cuda() 和 likelihood = likelihood.cuda() 将 GP 模型和其似然放到 GPU 上。

您可以在 GPyTorch 的文档中找到有关此主题的更多详细信息，链接地址为 mng.bz/lW8B。

我们还对梯度下降进行了 500 次迭代，但由于我们当前的数据集明显更大，这个循环可能需要一段时间才能完成（所以在等待时来杯咖啡吧！）。训练完成后，我们调用了之前编写的 visualize_gp_belief() 辅助函数来显示我们训练好的 GP 所做的预测，生成了图 12.4：

visualize_gp_belief(model, likelihood)

图 12.4 由普通 GP 进行的预测。预测很好地匹配了训练数据，但训练时间很长。

我们看到我们的 GP 的预测很好地匹配了训练数据点——这是一个鼓舞人心的迹象，表明我们的模型成功地从数据中学习了。然而，这个过程中存在一些问题。

12.1.3 普通 GP 训练中的问题

在这一小节中，我们讨论了在大型数据集上训练 GP 面临的一些挑战。首先，正如我们之前提到的，训练需要相当长的时间。在我的 MacBook 上，500 次梯度下降可能需要长达 45 秒的时间，这明显比我们在第二章和第三章观察到的情况要长。这直接是我们之前提到的 GP 的立方运行时间的结果，随着数据集的不断增大，这种长时间的训练只会变得更加禁锢，正如表 12.1 所示。

表 12.1 给定训练数据集大小的 GP 预计训练时间。训练很快变得困难。

训练集大小	训练时间
500 个点	45 秒
2,000 个点	48 分钟
3,000 个点	2.7 小时
5,000 个点	12.5 小时
10,000 个点	4 天

第二个，也许更令人担忧的问题源于这样一个事实，即随着训练数据的规模增加，计算损失函数（用于梯度下降的训练数据的边际对数似然）变得越来越困难。这可以通过 GPyTorch 在训练过程中打印的警告信息来看出：

NumericalWarning: CG terminated in 1000 iterations with average
  residual norm...

这些信息告诉我们，在计算损失时遇到了数值不稳定性。

注意，计算许多数据点之间的损失是一个计算上不稳定的操作。

数值不稳定性阻止我们正确计算损失，因此无法有效地最小化该损失。这在梯度下降的 500 次迭代中损失的变化中得到了说明，如图 12.5 所示。

图 12.5 在梯度下降过程中普通高斯过程的逐渐损失。由于数值不稳定性，损失曲线崎岖不平，无法有效地最小化。

与第二章和第三章中所见不同，我们这里的损失上下波动，表明梯度下降未能很好地最小化该损失。事实上，随着我们进行更多的迭代，我们的损失实际上增加了，这意味着我们得到了一个次优模型！这种现象是可以理解的：如果我们误计算了模型的损失，那么通过使用该误计算项来指导梯度下降中的学习，我们可能得到一个次优解。

你可能对梯度下降类比为下山的比喻很熟悉。假设你在山顶，想下山。沿途的每一步，你找到一个朝向能让你到达更低处的方向（即下降）。最终，经过足够的步骤，你抵达山脚。类似地，在梯度下降中，我们从一个相对较高的损失开始，并在每次迭代中调整我们模型的超参数，逐步降低损失。经过足够的迭代，我们到达最佳模型。

注涵盖了对梯度下降以及它如何类比为下山的出色讨论，可参考路易斯·塞拉诺的《深入理解机器学习》。

只有在我们能够准确计算损失时，这个过程才能正常运行，也就是说，我们可以清楚地看到哪个方向会让我们到达山上的更低处。然而，如果这个计算容易出错，我们自然无法有效地最小化模型的损失。这就好比戴着眼罩下山一样！正如我们在图 12.5 中看到的，我们实际上停留在山上的更高处（我们的损失高于梯度下降之前的值）。

图 12.6 使用数值不稳定的损失计算运行梯度下降类似于戴着眼罩下山。

总的来说，在大型数据集上训练常规高斯过程并不是一个好的方法。训练不仅随着训练数据规模的立方级增长，而且损失值的计算也不稳定。在本章的其余部分，我们将了解到变分高斯过程或 VGPs 是解决这个问题的方案。

12.2 从大型数据集中自动选择代表性点

VGPs 的思想是选择一组代表整个数据集的点，并在这个较小的子集上训练 GP。我们已经学会了如何在小数据集上训练 GP。希望这个较小的子集能捕捉到整个数据集的一般趋势，这样当在子集上训练 GP 时，仅有最少的信息会丢失。

这个方法非常自然。大数据集通常包含冗余信息，如果我们只从最信息丰富的数据点中学习，就可以避免处理这些冗余性。我们在 2.2 节中指出，像任何 ML 模型一样，GP 工作的假设是相似的数据点会产生相似的标签。当大数据集包含许多相似的数据点时，GP 只需要关注其中一个数据点来学习其趋势。例如，即使有按分钟的天气数据可用，天气预报模型也可以从仅有的小时测量中有效地进行学习。在这个小节中，我们将学习如何通过确保从小子集中学习相对于从大集合中学习时信息损失最小的方式来自动实现这一点，以及如何使用 GPyTorch 实现这个模型。

12.2.1 最小化两个 GP 之间的差异

我们如何选择这个小子集，使得最终的 GP 模型能够从原始数据集中获得最多的信息。在这个小节中，我们讨论了如何通过 VGP 来实现这一高级想法。这个过程等同于找到诱导点的子集，当在这个子集上训练 GP 时，诱导出的后验 GP 应该尽可能接近在整个数据集上训练的后验 GP。

当训练 VGP 时，深入一些数学细节，我们的目标是最小化在诱导点上条件化的后验 GP 和在整个数据集上条件化的后验 GP 之间的差异。这需要一种方法来衡量两个分布（两个 GP）之间的差异，而为此选择的衡量标准是 Kullback-Leibler（KL）散度。

定义Kullback-Leibler（KL）散度是一种统计距离，用于衡量两个分布之间的距离，也就是 KL 散度计算概率分布与另一个分布之间的不同程度。

KL 散度的补充材料

Will Kurt 的博客文章“Kullback-Leibler Divergence Explained”中提供了 KL 散度的直观解释（www.countbayesie.com/blog/2017/5/9/kullback-leibler-divergence-explained）。有数学背景的读者可以参考 David MacKay 的《信息论、推理和学习算法》第二章（剑桥大学出版社，2003 年）。

正如欧几里得空间中点 A 和点 B 之间的欧几里得距离（即连接两点的线段的长度）衡量了这两点在欧几里得空间中的距离一样，KL 散度衡量了概率分布空间中两个给定分布之间的距离，即它们彼此之间的差异有多大。这在图 12.7 中有所说明。

图 12.7 欧几里得距离衡量了平面上两点之间的距离，而 KL 散度衡量了两个概率分布之间的距离。

注作为一个数学上有效的距离度量，KL 散度是非负的。换句话说，任意两个分布之间的距离至少为零，当距离等于零时，两个分布完全匹配。

因此，如果我们能够轻松计算在诱导点上训练的后验 GP 与整个数据集上训练的后验 GP 之间的 KL 散度，我们应该选择使得 KL 散度为零的诱导点。不幸的是，类似于计算边际对数似然的计算不稳定性，计算 KL 散度也不容易。然而，由于其数学特性，我们可以将 KL 散度重写为两个量之间的差异，如图 12.8 所示。

图 12.8 KL 散度被分解为边际对数似然和证据下界（ELBO）之间的差异。ELBO 易于计算，因此被选择为要优化的度量。

这个方程中的第三项，也就是证据下界（ELBO），是边际对数似然和 KL 散度之间的精确差异。尽管边际对数似然和 KL 散度这两项很难计算，但 ELBO 具有简单的形式并且可以轻松计算。因此，我们可以最大化 ELBO 来间接最大化边际对数似然，而不是最小化 KL 散度，使得在诱导点上训练的后验 GP 尽可能接近在完整数据集上训练的后验 GP。

综上所述，为了找到一组诱导点，使得后验 GP 尽可能与我们能够在大数据集上训练时获得的 GP 相似，我们的目标是最小化两个 GP 之间的 KL 散度。然而，这个 KL 散度很难计算，所以我们选择优化 KL 散度的代理，即模型的 ELBO，这样更容易计算。正如我们在下一小节中看到的，GPyTorch 提供了一个方便的损失函数，用于计算这个 ELBO 项。在我们讨论实现之前，还有一件事情需要讨论：在最大化 ELBO 项时如何考虑大型训练集中的所有数据点。

12.2.2 在小批量中训练模型

由于我们的目标是找到最能代表整个训练数据集的一组感兴趣的点，我们仍然需要在计算 ELBO 时包含训练集中的所有点。但是我们之前说过，跨多个数据点计算边际对数似然是数值不稳定的，因此梯度下降变得无效。我们在这里面对相同的问题吗？在本小节中，我们看到当通过优化 ELBO 项来训练 VGP 时，我们可以通过使用更适合大型数据集的梯度下降的修改版本来避免这个数值不稳定性问题。

在许多数据点上计算 ML 模型的损失函数的任务并不是 GP 独有的。例如，神经网络通常在数千和数百万的数据点上进行训练，计算网络的损失函数对于所有数据点也是不可行的。对于这个问题的解决方法，对于神经网络和 VGP 都是近似通过对随机点的损失值计算跨所有数据点的真实损失值。例如，以下代码片段来自官方 PyTorch 文档，并显示了如何在图像数据集上训练神经网络（mng.bz/8rBB）。在这里，内部循环迭代训练数据的小子集，并在这些子集上计算的损失值上运行梯度下降：

for epoch in range(2):                        ❶

    running_loss = 0.0
    for i, data in enumerate(trainloader, 0):
        inputs, labels = data                 ❷

        optimizer.zero_grad()                 ❸

        outputs = net(inputs)                 ❹
        loss = criterion(outputs, labels)     ❹
        loss.backward()                       ❹
        optimizer.step()                      ❹

❶ 对数据集进行多次循环

❷ 获取输入；数据是一个 [输入，标签] 的列表

❸ 归零参数梯度

❹ 前向 + 反向 + 优化

当我们在少量点上计算模型的损失时，计算可以稳定有效地进行。此外，通过多次重复这个近似，我们可以很好地近似真实损失。最后，我们在这个近似损失上运行梯度下降，希望也最小化所有数据点上的真实损失。

定义在使用数据的随机子集计算的损失上运行梯度下降的技术有时被称为小批量梯度下降。在实践中，我们通常不会在梯度下降的每次迭代中随机选择一个子集，而是将训练集分成小的子集，并通过每个这些小的子集迭代计算近似损失。

例如，如果我们的训练集包含 1,000 个点，我们可以将其分成 10 个小子集，每个子集包含 100 个点。然后，我们对每个包含 100 个点的子集计算梯度下降的损失，并迭代重复所有 10 个子集。（这恰好是我们后面代码示例中所做的。）同样，虽然从数据子集计算的这种近似损失并不完全等于真实损失，但在梯度下降中，我们重复这个近似很多次，这在聚合中指引我们朝着正确的下降方向。

图 12.9 中说明了梯度下降最小化真实损失和小批量梯度下降最小化近似损失之间的区别。与梯度下降相比（再次强调，无法在大数据上运行），小批量版本可能不会指向最有效的下降方向，但通过多次重复近似，我们仍然能够达到目标。

图 12.9 梯度下降和小批量梯度下降在损失“谷底”中的示意图，其中谷底的中心给出最低的损失。梯度下降，如果可以计算，直接导向目标。小批量梯度下降朝着不是最优的方向前进，但最终仍然到达目标。

如果我们用盲目下山的比喻来思考，小批量梯度下降类似于戴着一块可以部分看穿的薄布的盲目。并不能保证我们每迈出一步就到达一个更低的位置，但是给定足够的时间，我们将能够成功下降。

注意并非所有的损失函数都可以通过对数据子集的损失来近似。换句话说，并非所有的损失函数都可以通过小批量梯度下降来最小化。GP 的负边际对数似然就是一个例子；否则，我们可以在这个函数上运行小批量梯度下降。幸运的是，小批量梯度下降适用于 VGP 的 ELBO。

总之，训练一个 VGP 遵循大致相似的程序，就像训练一个常规的 GP 一样，我们使用梯度下降的一个版本来最小化模型的适当损失。表 12.2 总结了两个模型类之间的关键差异：常规 GP 应该通过运行梯度下降来最小化精确的负边际对数似然在小数据集上训练，而 VGP 可以通过运行小批量梯度下降来优化 ELBO，在大数据集上训练，这是真实对数似然的一个近似。

表 12.2 训练一个 GP 与训练一个 VGP。高层次的过程是相似的；只有具体的组件和设置被替换。

训练过程	GP	VGP
训练数据大小	小	中等到大
训练类型	精确训练	近似训练
损失函数	负边际对数似然	ELBO
优化	梯度下降	小批量梯度下降

12.2.3 实现近似模型

我们现在准备在 GPyTorch 中实现一个 VGP。我们的计划是编写一个 VGP 模型类，这类似于我们已经使用的 GP 模型类，并使用小批量梯度下降最小化其 ELBO。我们在表 12.2 中描述的工作流程的不同之处反映在了本小节中的代码中。表 12.3 显示了在 GPyTorch 中实现 GP 与 VGP 时所需的组件。除了均值和协方差模块外，VGP 还需要另外两个组件：

变分分布—定义了 VGP 中诱导点的分布。正如我们在上一节中学到的，此分布要进行优化，以使 VGP 类似于在完整数据集上训练的 GP。
变分策略—定义了如何从诱导点产生预测。在第 2.2 节中，我们看到多元正态分布可以根据观测结果进行更新。这种变分策略促使对变分分布进行相同的更新。

表 12.3 在 GPyTorch 中实现 GP 与 VGP 时所需的组件。VGP 需要像 GP 一样的均值和协方差模块，但还需要变分分布和变分策略。

组件	GP	VGP
均值模块	是	是
协方差模块	是	是
变分分布	否	是
变分策略	否	是

考虑到这些组件，我们现在实现了 VGP 模型类，我们将其命名为 ApproximateGPModel。我们不再在 __init__() 方法中接受训练数据和似然函数。相反，我们接受一组诱导点，这些点将用于表示整个数据集。 __init__() 方法的其余部分包括声明将用于学习哪组诱导点最佳的学习流程：

变分分布 variational_distribution 变量是 CholeskyVariationalDistribution 类的一个实例，在初始化期间接受诱导点的数量。 变分分布是 VGP 的核心。
变分策略 variational_strategy 变量是 VariationalStrategy 类的一个实例。它接受一组诱导点以及变分分布。我们将 learn_inducing_locations = True 以便在训练过程中学习这些诱导点的最佳位置。如果将此变量设置为 False，则传递给 __init__() 的点（存储在 inducing 中）将用作诱导点：

class ApproximateGPModel(gpytorch.models.ApproximateGP):             ❶
  def __init__(self, inducing_points):                               ❷
    variational_distribution =                                       ❸
    ➥gpytorch.variational.CholeskyVariationalDistribution(          ❸
        inducing_points.size(0)                                      ❸
    )                                                                ❸
    variational_strategy = gpytorch.variational.VariationalStrategy( ❸
        self,                                                        ❸
        inducing_points,                                             ❸
        variational_distribution,                                    ❸
        learn_inducing_locations=True,                               ❸
    )                                                                ❸
    super().__init__(variational_strategy)                           ❸

    ...                                                              ❹

❶ 我们的 VGP 不是一个 ExactGP 对象，而是一个近似 GP。

❷ 接受一组初始诱导点

❸ 设置训练所需的变分参数

❹ 待续

在 __init__() 方法的最后一步中，我们声明了 VGP 的均值和协方差函数。它们应该与在数据上训练的普通 GP 中要使用的函数相同。在我们的情况下，我们使用常数均值和具有输出比例的 RBF 核：

class ApproximateGPModel(gpytorch.models.ApproximateGP):
    def __init__(self, inducing_points):
        ...

        self.mean_module = gpytorch.means.ConstantMean()
        self.covar_module = gpytorch.kernels.ScaleKernel(
            gpytorch.kernels.RBFKernel()
        )

我们还以与常规 GP 相同的方式声明 forward() 方法，这里不再展示。现在，让我们用训练集中的前 50 个数据点来初始化这个模型作为引入点：

model = ApproximateGPMod1el(train_x[:50, :])            ❶
likelihood = gpytorch.likelihoods.GaussianLikelihood()

❶ 切片张量 train_x[:50, :] 给出了 train_x 中的前 50 个数据点。

前 50 个数据点没有什么特别之处，它们的值，在 VGP 模型内部存储，将在训练过程中被修改。这种初始化最重要的部分是，我们指定模型应该使用 50 个引入点。如果我们想要使用 100 个，我们可以将 train_x[:100, :] 传递给初始化。

很难准确地说，多少个引入点足以用于 VGP。我们使用的点越少，模型训练就越快，但是这些引入点在表示整个集合时的效果就越不好。随着点数的增加，VGP 有更多的自由度来展开引入点以覆盖整个集合，但训练会变慢。

一般规则是不超过 1,000 个引入点。正如我们马上要讨论的，50 个点足以让我们以高保真度近似前一小节中训练的 GP。

要设置小批量梯度下降，我们首先需要一个优化器。我们再次使用 Adam 优化器：

optimizer = torch.optim.Adam(
    [
        {"params": model.parameters()},
        {"params": likelihood.parameters()}    ❶
    ],
    lr=0.01
)

❶ 优化似然函数的参数以及 GP 的参数

要优化的参数

之前，我们只需要将 model.parameters() 传递给 Adam。在这里，似然函数没有与 VGP 模型耦合——常规 GP 使用似然函数初始化，而 VGP 则没有。因此，在这种情况下，有必要将 likelihood.parameters() 传递给 Adam。

对于损失函数，我们使用 gpytorch.mlls.VariationalELBO 类，该类实现了我们希望通过 VGP 优化的 ELBO 量。在初始化期间，此类的一个实例接受似然函数、VGP 模型和完整训练集的大小（我们可以通过 train_y.size(0) 访问）。有了这些，我们声明这个对象如下：

mll = gpytorch.mlls.VariationalELBO(
    likelihood,
    model,
    num_data=train_y.size(0)    ❶
)

❶ 训练数据的大小

有了模型、优化器和损失函数设置好了，我们现在需要运行小批量梯度下降。为此，我们将训练数据集分成批次，每个批次包含 100 个点，使用 PyTorch 的 TensorDataset 和 DataLoader 类：

train_dataset = torch.utils.data.TensorDataset(train_x, train_y)
train_loader = torch.utils.data.DataLoader(train_dataset, batch_size=100)

这个 train_loader 对象允许我们以干净的方式迭代大小为 100 的数据集的小批量，当运行梯度下降时。损失—即 ELBO—通过以下语法计算：

output = model(x_batch)
loss = -mll(output, y_batch)

这里，x_batch 和 y_batch 是完整训练集的一个给定批次（小子集）。总的来说，梯度下降的实现如下：

model.train()                                 ❶
likelihood.train()                            ❶

for i in tqdm(range(50)):                     ❷
    for x_batch, y_batch in train_loader:     ❸
        optimizer.zero_grad()                 ❹

        output = model(x_batch)               ❹
        loss = -mll(output, y_batch)          ❹

        loss.backward()                       ❹
        optimizer.step()                      ❹

model.eval()                                  ❺
likelihood.eval()                             ❺

❶ 启用训练模式

❷ 迭代整个训练数据集 50 次

❸ 在每次迭代中，迭代 train_loader 中的小批次

❹ 小批量梯度下降，在批次上运行梯度下降

❺ 启用预测模式

在运行这个小批量梯度下降循环时，你会注意到它比使用普通 GP 的循环要快得多。（在同一台 MacBook 上，这个过程只需不到一秒钟，速度大幅提升！）

VGP 的速度

你可能认为我们在普通 GP 中进行的 500 次迭代与 VGP 中进行的 50 次小批量梯度下降的比较不公平。但是请记住，在小批量梯度下降的外部for循环的每次迭代中，我们还要在train_loader中迭代 10 个小批量，所以最终总共进行了 500 次梯度步骤。此外，即使我们运行了 500 次小批量梯度下降的迭代，也只需要不到 1 秒乘以 10，仍然比 45 秒快 4 倍。

因此，通过小批量梯度下降，我们的 VGP 模型可以更高效地训练。但是对于训练质量如何呢？图 12.10 的左侧面板可视化了我们小批量梯度下降运行期间的逐渐 ELBO 损失。与图 12.5 相比，尽管损失并没有在每一步时持续减小（存在锯齿趋势），但损失在整个过程中有效地最小化了。

图 12-10.png

图 12.10 VGP 在小批量梯度下降期间逐步损失服从的长度尺度和输出尺度的对应关系

这表明，在优化过程中的每一步可能不是最小化损失的最佳方向，但小批量梯度下降确实是有效地减小了损失。这在图 12.11 中更清楚地显示出来。

图 12-11.png

图 12.11 过程中 VGP 的逐步损失在小批量梯度下降中逐渐减小。尽管存在一定的变异性，但损失有效地最小化了。

现在，让我们可视化这个 VGP 模型产生的预测结果，看看它是否产生了合理的结果。使用visualize_gp_belief()助手函数，我们得到图 12.12，显示我们以较小的时间成本获得了对真实损失进行训练的 GP 的高质量近似。

图 12-12.png

图 12.12 由 GP 和 VGP 进行的预测。VGP 进行的预测与 GP 的预测大致相同。

结束我们对 VGP 的讨论，让我们可视化 VGP 模型学到的诱导点的位置。我们已经说过这些诱导点应该代表整个数据集并很好地捕捉其趋势。要绘制诱导点，我们可以使用model.variational_strategy.inducing_points.detach()访问它们的位置，并将它们作为散点沿着均值预测进行绘制。当调用这个函数时，我们只需要将variational设置为True：

visualize_gp_belief(model, likelihood, variational=True)

这产生了图 12.13，在其中我们看到这些感应点的非常有趣的行为。它们并不均匀分布在我们的训练数据之间；相反，它们聚集在数据的不同部分。这些部分是目标函数上升或下降或呈现某些非平凡行为的地方。通过将感应点分配到这些位置，VGP 能够捕捉嵌入大型训练数据集中最重要的趋势。

图 12.13 VGP 的感应点。这些点被放置在整个数据中，捕捉最重要的趋势。

我们已经学习了如何使用小批量梯度下降来训练 VGP，并且已经看到这有助于以更低的成本近似于不可训练的常规 GP。在接下来的部分，我们将学习另一种梯度下降算法，可以更有效地训练 VGP。

12.3 通过考虑损失表面的几何特性来进行更好的优化

在本节中，我们将学习一种名为自然梯度下降的算法，这是梯度下降的另一种版本，在计算下降步骤时更仔细地推理损失函数的几何结构。正如我们很快会看到的那样，这种谨慎的推理使我们能够快速降低损失函数，最终导致更有效的优化，迭代次数更少（也就是更快的收敛）。

要理解自然梯度下降的动机以及为什么它比我们已经拥有的更好，我们首先区分 VGP 的两种参数类型：

第一种类型是 GP 的常规参数，例如均值常数和协方差函数的长度和输出比例。这些参数取得常规的数值，存在于欧几里得空间中。
第二种类型包括只有 VGP 才有的变分参数。这些与感应点和促使变分分布近似所需的各种组件有关。换句话说，这些参数与概率分布相关，并且具有无法在欧几里得空间内很好表示的值。

注意：这两种参数之间的差异有些相似，尽管不完全类似于欧几里得距离可以衡量该空间中两点之间的距离，但无法衡量两个概率分布之间的差异。

尽管我们在前一节中使用的小批量梯度下降效果已经足够好，但该算法假设所有参数存在于欧几里得空间中。例如，从算法的角度来看，长度尺度从 1 变到 2 的差异与诱导点的均值从 1 变到 2 的差异是相同的。然而，事实并非如此：从长度尺度从 1 变到 2 会对 VGP 模型产生非常不同的影响，而从诱导点的均值从 1 变到 2 的影响也不同。这在图 12.14 的示例中有所体现，其中损失对长度尺度的行为与对诱导均值的行为非常不同。

图 12.14 一个示例，说明了要最小化的损失可能在正常参数和变分参数方面的行为非常不同。这就需要考虑损失的几何形状。

正常参数和 VGP 的变分参数的损失函数的几何形状之间存在这种行为差异是因为。如果小批量梯度下降能够在计算损失的下降方向时考虑到这种几何差异，那么该算法在最小化损失时将更有效。这就是自然梯度下降的作用所在。

自然梯度下降 的定义利用了关于损失函数对于变分参数的几何特性的信息，以计算这些参数更好的下降方向。

通过采用更好的下降方向，自然梯度下降可以更有效地帮助我们优化 VGP 模型，而且更快。最终的结果是我们能够在更少的步骤中收敛到我们的最终模型。继续我们对不同梯度下降算法的二维图示，图 12.15 展示了这种几何推理如何帮助自然梯度下降比小批量梯度下降更快地达到目标。也就是说，自然梯度下降在训练过程中往往需要更少的步骤来达到与小批量梯度下降相同的损失。在我们下山的类比中，使用自然梯度下降时，我们在试图下山时仍然被一块薄布蒙住了眼睛，但现在我们穿着特制的登山鞋，可以更有效地穿越地形。

图 12.15 展示了在损失“谷底”中进行梯度下降、小批量梯度下降和自然梯度下降的插图，在谷底的中心给出了最低的损失。通过考虑损失函数的几何特性，自然梯度下降比小批量梯度下降更快地达到了损失最小值。

自然梯度下降的补充材料

对于更加数学化的自然梯度下降解释，我推荐阅读 Agustinus Kristiadi 的优秀博文“自然梯度下降”：mng.bz/EQAj。

注意需要注意的是，自然梯度下降算法仅优化 VGP 的变分参数。常规参数，例如长度和输出尺度，仍然可以通过常规小批量梯度下降算法优化。在接下来实现新的训练过程时，我们将看到这一点。

因此，让我们使用自然梯度下降来训练我们的 VGP 模型。与之前部分的相同一维目标函数一样，我们实现了一个能够与自然梯度下降一起工作的 VGP 模型。这种情况下的模型类似于我们在上一节为小批量梯度下降实现的ApproximateGPModel，它

仍然扩展gpytorch.models.ApproximateGP
需要一个变分策略来管理学习过程
具有类似常规 GP 模型的均值函数、协方差函数和forward()方法

这里唯一的区别是，当训练模型时，变分分布需要是gpytorch.variational.NaturalVariationalDistribution的一个实例，以便我们在使用自然梯度下降时进行训练。整个模型类实现如下：

class NaturalGradientGPModel(gpytorch.models.ApproximateGP):
  def __init__(self, inducing_points):
    variational_distribution =                         ❶
      gpytorch.variational.                            ❶
      ➥NaturalVariationalDistribution(                ❶
        inducing_points.size(0)                        ❶
    )                                                  ❶

    variational_strategy = gpytorch.variational.
    ➥VariationalStrategy(                             ❷
        self,                                          ❷
        inducing_points,                               ❷
        variational_distribution,                      ❷
        learn_inducing_locations=True,                 ❷
    )                                                  ❷
    super().__init__(variational_strategy)             ❷
    self.mean_module = gpytorch.means.ConstantMean()   ❷
    self.covar_module = gpytorch.kernels.ScaleKernel(  ❷
        gpytorch.kernels.RBFKernel()                   ❷
    )                                                  ❷
  def forward(self, x):
        ...                                            ❸

❶ 变分分布需要是自然的才能与自然梯度下降一起工作。

❷ 声明剩余的变分策略与以前相同。

❸ forward()方法与以前相同。

我们再次使用 50 个引导点初始化此 VGP 模型：

model = NaturalGradientGPModel(train_x[:50, :])         ❶
likelihood = gpytorch.likelihoods.GaussianLikelihood()

❶ 50 个引导点

现在是重要的部分，我们为训练声明优化器。请记住，我们使用自然梯度下降算法来优化模型的变分参数。然而，其他参数，例如长度和输出尺度，仍然必须由 Adam 优化器优化。因此，我们使用以下代码：

ngd_optimizer = gpytorch.optim.NGD(                  ❶
  model.variational_parameters(), num_data=train_y.  ❶
  ➥size(0), lr=0.1                                  ❶
)                                                    ❶

hyperparam_optimizer = torch.optim.Adam(             ❷
  [{"params": model.parameters()}, {"params":        ❷
  ➥likelihood.parameters()}],                       ❷
  lr=0.01                                            ❷
)                                                    ❷
mll = gpytorch.mlls.VariationalELBO(
  likelihood, model, num_data=train_y.size(0)
)

❶ 自然梯度下降接受 VGP 的变分参数，model.variational_parameters()。

❷ Adam 接受 VGP 的其他参数，model.parameters()和likelihood.parameters()。

现在，在训练期间，我们仍然使用以下方式计算损失

output = model(x_batch)
loss = -mll(output, y_batch)

在计算损失时，我们循环遍历我们训练数据的小批量（x_batch和y_batch）。然而，现在我们有两个优化器同时运行，我们需要通过在每次训练迭代时调用zero_grad()（清除前一步的梯度）和step()（执行一步下降）来管理它们：

model.train()                              ❶
likelihood.train()                         

for i in tqdm(range(50)):
    for x_batch, y_batch in train_loader:
        ngd_optimizer.zero_grad()          ❷
        hyperparam_optimizer.zero_grad()   ❷

        output = model(x_batch)
        loss = -mll(output, y_batch)

        loss.backward()

        ngd_optimizer.step()               ❸
        hyperparam_optimizer.step()        ❸

model.eval()                               ❹
likelihood.eval()                          ❹

❶ 启用训练模式

❷ 清除前一步的梯度

❸ 使用每个优化器执行下降步骤

❹ 启用预测模式

注意像往常一样，在梯度下降之前，我们需要调用model.train()和likelihood.train()，在训练完成后需要调用model.eval()和likelihood.eval()。

请注意，我们在自然梯度下降优化器和 Adam 上都调用了 zero_grad() 和 step()，以优化 VGP 模型的相应参数。训练循环再次只需很短的时间即可完成，并且训练得到的 VGP 产生的预测结果如图 12.16 所示。我们看到的预测结果与图 12.4 中的普通 GP 和使用小批量梯度下降训练的 VGP 在图 12.13 中的预测结果非常相似。

图 12.16 由 VGP 和通过自然梯度下降训练的感应点所做出的预测。这些预测的质量很高。

我们可以进一步检查训练过程中 ELBO 损失的逐步进展。其进展在图 12.17 的左侧面板中可视化。

图 12.17 VGP 在自然梯度下降期间的逐步损失。经过少数迭代后，损失得到了有效最小化。

令人惊讶的是，在训练过程中，我们的 ELBO 损失几乎立即降至一个较低的值，这表明自然梯度下降能够帮助我们快速收敛到一个良好的模型。这说明了当训练 VGP 时，这种梯度下降算法的变体的好处。

我们现在已经到达第十二章的结尾。在本章中，我们学习了如何通过使用感应点将 GP 模型扩展到大型数据集，感应点是一组代表性点，旨在捕获大型训练集所展现的趋势。由此产生的模型称为 VGP，它适用于小批量梯度下降，因此可以在不计算所有数据点的模型损失的情况下进行训练。我们还研究了自然梯度下降作为小批量算法的更有效版本，使我们能够更有效地优化。在第十三章中，我们将涵盖 GP 的另一种高级用法，将其与神经网络结合以建模复杂的结构化数据。

12.4 练习

此练习演示了在加利福尼亚房价的真实数据集上，从普通 GP 模型转换为 VGP 模型时效率的提高。我们的目标是观察 VGP 在实际环境中的计算优势。

完成以下步骤：

使用 Pandas 库中的 read_csv() 函数读取存储在名为 data/housing.csv 的电子表格中的数据集，该数据集来自 Kaggle 的加利福尼亚房价数据集（mng.bz/N2Q7），使用的是创意共用公共领域许可证。一旦读取，Pandas dataframe 应该看起来类似于图 12.18 中的输出。

图 12.18 作为 Pandas dataframe 显示的房价数据集。这是本练习的训练集。
在散点图中可视化median_house_value列，这是我们的预测目标，其x-和y-轴对应于longitude和latitude列。一个点的位置对应于一个房子的位置，点的颜色对应于价格。可视化效果应该类似于图 12.19。

图 12.19 房价数据集的散点图表示
提取除最后一列（median_house_value）以外的所有列，并将它们存储为 PyTorch 张量。这将用作我们的训练特征，train_x。
提取median_house_value列，并将其对数变换存储为另一个 PyTorch 张量。这是我们的训练目标，train_y。
通过减去均值并除以标准差来归一化训练标签train_y。这将使训练更加稳定。
使用具有恒定均值函数和具有自动相关性确定（ARD）的 Matérn 5/2 核实现常规 GP 模型。关于 Matérn 核和 ARD 的详细内容，请参阅第 3.4.2 和第 3.4.3 节。
使用以下代码创建一个噪声至少为 0.1 的似然性：
```
likelihood = gpytorch.likelihoods.GaussianLikelihood(
    noise_constraint=gpytorch.constraints.GreaterThan(1e-1)  ❶
)
```
❶该约束强制噪声至少为 0.1。

此约束有助于通过提高噪声容限来平滑训练标签。
初始化之前实现的 GP 模型，并使用梯度下降对其进行 10 次迭代训练。观察总训练时间。
实现具有与 GP 相同的均值和协方差函数的变分 GP 模型。这个模型看起来类似于我们在本章中实现的ApproximateGPModel类，只是现在我们需要 Matérn 5/2 核和 ARD。
使用类似初始化的似然性和 100 个诱导点对此 VGP 进行训练，使用自然梯度下降进行 10 次迭代。对于小批量梯度下降，您可以将训练集分成大小为 100 的批次。
验证训练 VGP 所需的时间是否少于训练 GP 的时间。对于计时功能，您可以使用time.time()来记录每个模型训练的开始和结束时间，或者您可以使用tqdm库来跟踪训练持续时间，就像我们一直在使用的代码一样。

解决方案包含在CH11/02 - Exercise.ipynb笔记本中。

摘要

GP 的计算成本与训练数据集的大小呈立方比例。因此，随着数据集的大小增长，训练模型变得不可行。
在大量数据点上计算 ML 模型的损失在数值上是不稳定的。以不稳定的方式计算损失可能会误导梯度下降过程中的优化，导致预测性能不佳。
VGP 通过仅对一小组诱导点进行训练来扩展到大型数据集。这些诱导点需要代表数据集，以便训练的模型尽可能地与在整个数据集上训练的 GP 相似。
为了产生一个尽可能接近在所有训练数据上训练的模型的近似模型，Kullback-Leibler 散度，它度量两个概率分布之间的差异，被用于 VGP 的制定中。
当训练一个 VGP 时，证据下界（ELBO）充当真实损失的代理。更具体地说，ELBO 限制了模型的边际对数似然，这是我们的优化目标。通过优化 ELBO，我们间接优化了边际对数似然。
训练一个 VGP 可以通过小批量进行，从而实现更稳定的损失计算。在该过程中使用的梯度下降算法是小批量梯度下降。
尽管小批量梯度下降的每一步都不能保证完全最小化损失，但是当运行大量迭代时，该算法可以有效降低损失。这是因为许多小批量梯度下降的步骤在聚合时可以指向最小化损失的正确方向。
自然梯度下降考虑了相对于 VGP 的变分参数的损失函数的几何性质。这种几何推理使得算法能够更新训练模型的变分参数并更有效地最小化损失，从而实现更快的收敛。
自然梯度下降优化了 VGP 的变分参数。诸如长度和输出比例等常规参数由小批量梯度下降进行优化。

第十四章：将高斯过程与神经网络结合

本章涵盖

使用常见协方差函数处理复杂结构化数据的困难
使用神经网络处理复杂结构化数据
将神经网络与 GP 结合

在第二章中，我们学到了高斯过程（GP）的均值和协方差函数作为我们希望在模型中融入的先验信息，当进行预测时。因此，这些函数的选择极大地影响了训练后的 GP 的行为。因此，如果均值和协方差函数被错误地指定或不适用于手头的任务，那么得到的预测就不会有用。

例如，记住covariance function，或kernel，表示两个点之间的相关性——即相似性。两个点越相似，它们的标签值可能越相似，我们试图预测的标签值。在我们的房价预测示例中，相似的房子可能会有类似的价格。

核到底如何计算两个给定房子之间的相似性？我们考虑两种情况。在第一种情况中，核函数仅考虑前门的颜色，并对于任何具有相同门颜色的两个房子输出 1，否则输出 0。换句话说，如果两个房子的前门颜色相同，这个核函数认为它们相似。

正如图 13.1 所示，这个核函数对于房价预测模型来说是一个糟糕的选择。核函数认为左边的房子和中间的房子应该有相似的价格，而右边的房子和中间的房子应该有不同的价格。这是不合适的，因为左边的房子比其他两个大得多，而其他两个房子的大小相似。误差发生的原因是核函数错误地判断了房子的哪个特征是房子成本的良好预测特征。

图 13.1 由不合适的核函数计算的房屋之间的协方差。因为它只关注前门的颜色，所以这个核函数无法产生合适的协方差。

另一个核函数更加复杂，并考虑了相关因素，比如位置和居住面积。这个核函数更加合适，因为它可以更合理地描述两个房子之间的价格相似性。拥有合适的核函数——即正确的相似度度量——对于 GP 来说至关重要。如果核函数能够正确地描述给定数据点之间的相似性或差异性，那么使用协方差的 GP 将能够产生良好校准的预测。否则，预测将具有较低的质量。

你可能会认为一个只考虑门颜色的房屋核函数是不合适的，并且在 ML 中没有合理的核函数会表现出这种行为。然而，正如我们在本章中所展示的，到目前为止我们使用的一些常见核函数（例如，RBF 和 Matérn）在处理结构化输入数据（例如图像）时会出现相同的问题。具体来说，它们未能充分描述两个图像之间的相似性，这给在这些结构化数据类型上训练 GPs 带来了挑战。我们采取的方法是使用神经网络。神经网络是灵活的模型，可以在有足够数据的情况下很好地逼近任何函数。我们学会使用神经网络来转换 GP 的核函数无法很好地处理的输入数据。通过这样做，我们既得到了神经网络的灵活建模，又从 GP 中获得了不确定性校准的预测。

在本章中，我们展示了我们通常的 RBF 核函数不能很好地捕捉常见数据集的结构，从而导致 GP 的预测不佳。然后我们将一个神经网络模型与此 GP 结合起来，看到新的核函数可以成功地推理相似性。到本章结束时，我们获得了一个框架，帮助 GP 处理结构化数据类型并提高预测性能。

13.1 包含结构的数据

在本节中，我们解释了结构化数据的确切含义。与我们在之前章节中用来训练 GPs 的数据类型不同，在那些数据类型中，数据集中的每个特征（列）可以取得连续范围内的值，而在许多应用中，数据具有更复杂性。比如说：

房子的楼层数只能是正整数。
在计算机视觉任务中，图像中的像素值是 0 到 255 之间的整数。
在分子 ML 中，分子通常被表示为图形。

那就是，这些应用中的数据点中嵌入了结构，或者需要数据点遵循的要求：没有房子可以有负数的楼层；像素不能以分数作为其值；表示分子的图形将具有表示化学物质和结合物的节点和边缘。我们称这些类型的数据为结构化数据。在本章中，我们将使用流行的 MNIST 手写数字数据集（见huggingface.co/datasets/mnist）作为我们讨论的案例研究。

定义修改后的美国国家标准与技术研究所（MNIST）数据集包含手写数字的图像。每个图像是一个 28×28 的整数矩阵，取值范围在 0 到 255 之间。

这个数据集中的一个示例数据点如图 13.2 所示，其中像素的阴影对应于其值；0 对应于白色像素，255 对应于黑色像素。我们看到这个数据点是数字五的图像。

图 13.2 来自 MNIST 数据集的数据点，这是一个由 28 行和 28 列像素组成的图像，表示为一个 PyTorch 张量

注意虽然这个手写数字识别任务在技术上是一个分类问题，但我们使用它来模拟一个回归问题（这是我们在 BayesOpt 中要解决的问题类型）。由于每个标签都是一个数字（一个数字），我们假装这些标签存在于一个连续的范围内，并直接将它们用作我们的预测目标。

我们的任务是在一个图像标签数据集上训练一个 GP（每个标签都是对应图像中写的数字的值），然后在一个测试集上进行预测。这个区别在图 13.3 中有所体现，它显示与分类不同，在分类中，我们选择一个类作为每个数据点的预测，而在回归任务中，这里的每个预测是一个连续范围内的数字。

图 13.3 在 MNIST 数据的上下文中的分类与回归。每个预测是一个分类任务，对应于其中的一个类别；在回归中的每个预测是一个连续范围内的数字。

有许多现实世界的应用遵循这种结构化数据的回归问题的形式：

在产品推荐中，我们希望预测某人点击自定义广告的概率。广告是可以自定义的图片，是结构化数据，点击概率是预测目标。这个概率可以是 0 到 1 之间的任何数字。
在材料科学中，科学家可能希望在实验室中合成某种分子组合时预测其能量水平。每种分子组合都可以表示为具有节点和边的特定结构的图，并且能量水平可以是理论最小和最大能量水平之间的任何数字，一个组合可能表现出的。
在药物发现中，我们希望预测可能产生的药物的有效性。如图 13.4 所示，每种药物对应于一种化合物，它也可以表示为一个图。其有效性可以是一个实数，介于某个范围内（比如从 0 到 10）。

图 13.4 药物发现作为一个结构化回归问题的例子。每个化合物都表示为一个结构化图，并且我们的目标是预测这种化合物在治疗某种疾病方面的有效性，其范围从 0 到 10。

在所有这些应用中，我们想要进行预测的输入数据是结构化的，我们的预测目标是一个实数值。简而言之，它们是结构化数据的回归问题。使用 MNIST 数据集，我们模拟了这样一个问题。

13.2 捕捉结构化数据内的相似性

在本节中，我们将探讨常见内核（如径向基函数内核）如何无法描述结构化数据中的相似性。量化两个输入的协方差的内核输出，x[1] 和 x[2] 的输出定义如下：

这个输出是两个变量之间的协方差，是两个输入之间差异的负指数除以一个标度。输出始终在 0 和 1 之间，而且越大的差异，输出越小。

这在许多情况下是有道理的，因为如果两个输入具有类似的值，因此差异很小，则它们的协方差将很高；而如果它们具有不同的值，则协方差将很低。两栋面积大致相等的房屋可能会有类似的价格，也就是说，它们的价格具有高的协方差；另一方面，非常大和非常小的房子的价格可能会具有较低的协方差。

13.2.1 使用 GPyTorch 的内核

让我们用代码验证一下。当我们创建 GP 模型时，通常会初始化一个 RBFKernel 对象。这里，我们直接使用这个内核对象进行工作。为此，我们首先使用 GPyTorch 创建一个 RBF 内核对象：

import gpytorch

rbf_kernel = gpytorch.kernels.RBFKernel()

请注意，作为 Python 中实现 GP 相关对象的首选库，我们一如既往地使用 GPyTorch。有关如何在 GPyTorch 中使用内核对象的详细信息，请参见第 2.4 节。

要计算两个输入之间的协方差，我们只需将它们传递给该内核对象即可。例如，让我们计算 0 和 0.1 之间的协方差：

>>> rbf_kernel(torch.tensor([0.]), torch.tensor([0.1])).evaluate().item()
0.9896470904350281

这两个数字在实数线上非常接近（也就是说，它们是相似的），因此它们的协方差非常高，几乎为 1。现在让我们计算 0 和 10 之间的协方差，这是两个不同的数字：

>>> rbf_kernel(torch.tensor([0.]), torch.tensor([10.])).evaluate().item()
0.0

这次，由于两个数字之间的差异要大得多，它们的协方差降为 0。这种对比是合理的行为，并且通过图 13.5 进行说明。

图 13.5 各种数字之间的协方差。两个数字之间的差异较小时，协方差增加；差异较大时，协方差降低。

当两个输入之间的值差异不能捕捉到数据结构差异时，问题就出现了。这通常是对结构化数据（如图像）的情况。接下来，我们将看到像径向基函数（RBF）这样的常见内核如何无法描述结构化数据中的相似性。

13.2.2 在 PyTorch 中处理图像

在这一小节中，我们将看到如何将图像导入和存储为 PyTorch 张量，以及如何处理此类数据时，基于值的相似度度量（如 RBF 内核）如何失效。首先，我们重新定义我们的 RBF 内核，使其具有较大的长度尺度，因此更有可能出现较高的协方差：

rbf_kernel = gpytorch.kernels.RBFKernel()
rbf_kernel.lengthscale = 100               ❶

❶ 长度尺度越大，协方差越高。

现在，我们需要将 MNIST 数据集中的图像导入到我们的 Python 代码中。我们可以使用 PyTorch 及其流行的附加库 torchvision 来实现：

import torch
from torchvision import datasets, transforms

transform = transforms.Compose([                               ❶
    transforms.ToTensor(),                                     ❶
    transforms.Normalize((0.1307,), (0.3081,))                 ❶
])                                                             ❶

dataset = datasets.MNIST(                                      ❷
    "../data", train=True, download=True, transform=transform  ❷
)                                                              ❷

train_x = dataset.data.view(-1, 28 * 28)                       ❸

❶ 定义规范化像素值的转换

❷ 下载并导入数据集

❸ 提取像素值作为一个扁平化的张量

我们不会深入研究这段代码，因为它不是我们讨论的重点。我们只需要知道 train_x 包含 MNIST 数据集中的图像，每个图像都存储为一个 PyTorch 张量，其中包含表示手写数字图像的像素值。

由于数据点是图像，我们可以将它们可视化为热图，使用 Matplotlib 中熟悉的 imshow() 函数。例如，以下代码可视化了 train_x 中的第一个数据点：

plt.figure(figsize=(8, 8))

plt.imshow(train_x[0, :].view(28, 28));    ❶

❶ 每个图像有 28 行和 28 列的像素，因此我们需要将其重塑为一个 28×28 的方形张量。

这段代码生成了图 13.2，我们看到它是数字 5 的图像。当我们打印出这个第一个数据点的实际值时，我们看到它是一个 28 × 28 = 784 元素的 PyTorch 张量：

>>> train_x[0, :]
tensor([ 0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   3,  18,
        18,  18, 126, 136, 175,  26, 166, 255, 247, 127,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,  30,  36,  94, 154, 170, 253,
       253, 253, 253, 253, 225, 172, 253, 242, 195,  64,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,  49, 238, 253, 253, 253, 253, 253,
       253, 253, 253, 251,  93,  82,  82,  56,  39,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,  18, 219, 253, 253, 253, 253, 253,
       198, 182, 247, 241,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,  80, 156, 107, 253, 253, 205,
        11,   0,  43, 154,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,  14,   1, 154, 253,  90,
[output truncated]

此张量中的每个元素范围在 0 到 255 之间，表示我们在图 13.2 中看到的像素。值为 0 对应于最低信号，即图中的背景，而较高的值对应于亮点。

13.2.3 计算两个图像的协方差

这就是我们探索普通 GP 核处理结构化数据时所需要的所有背景信息。为了突出问题，我们单独提出了三个特定的数据点，分别称为点 A、点 B 和点 C，它们的索引如下：

ind1 = 304    ❶
ind2 = 786    ❷
ind3 = 4      ❸

❶ 点 A

❷ 点 B

❸ 点 C

在检查这些图像实际显示的数字之前，让我们使用我们的 RBF 核来计算它们的协方差矩阵：

>>> rbf_kernel(train_x[[ind1, ind2, ind3], :]).evaluate()
tensor([[1.0000e+00, 4.9937e-25, 0.0000e+00],
        [4.9937e-25, 1.0000e+00, 0.0000e+00],
        [0.0000e+00, 0.0000e+00, 1.0000e+00]], ...)

这是一个 3×3 协方差矩阵，具有熟悉的结构：对角线元素取值为 1，表示各个变量的方差，而非对角线元素表示不同的协方差。我们看到点 A 和 C 完全不相关，协方差为零，而点 A 和 B 稍微相关。根据 RBF 核，点 A 和 B 相似，并且与点 C 完全不同。

我们应该期望点 A 和 B 具有相同的标签。然而，事实并非如此！再次将这些数据点可视化为热图，我们得到图 13.6。

图 13.6 MNIST 数据集中的三个特定数据点。第一个和第二个点具有非零协方差，尽管标签不同。第一个和第三个点具有零协方差，尽管标签相同。

在这里，点 A 和 C 共享相同的标签（数字 9）。那么为什么 RBF 核会认为点 A 和 B 有相关性呢？看图 13.6，我们可以猜测，虽然点 A 和 B 有不同的标签，但图像本身在很多像素上是相似的。事实上，构成数字尾巴的笔画在这两幅图像中几乎完全相同。因此，在某种程度上，RBF 核正在做它的工作，根据这种差异计算图像之间的差异并输出代表它们协方差的数字。然而，这种差异是通过比较像素本身来计算的，这不是我们试图学习的指标：数字的值。

通过仅仅查看像素值，RBF 核高估了点 A 和 B 之间的协方差，这两个点具有不同的标签，并低估了点 A 和 C 之间的协方差，它们具有相同的标签，正如图 13.7 所示。这里可以使用类比来演示我们在本章开头提到的不恰当的 house 核：这个核只看前门的颜色来决定两个房屋是否相关，导致对它们价格的不准确预测。类似地（但不如此极端），RBF 核在比较两幅图像时只考虑像素的值，而不考虑更高级别的模式，这导致了较差的预测性能。

图 13.7 由 RBF 核计算的手写数字之间的协方差。因为它只看像素值，所以 RBF 核无法产生适当的协方差。

13.2.4 在图像数据上训练 GP

通过使用错误的相似度度量标准，RBF 混淆了点 B 和 C 中哪个与点 A 相关联，这导致在训练 GP 时产生了不良结果。我们再次使用 MNIST 数据集，这次提取 1,000 个数据点作为训练集，另外 500 个数据点作为测试集。我们的数据准备和学习工作流程总结在图 13.8 中，我们将详细介绍其中的不同步骤。

图 13.8 在 MNIST 上进行 GP 学习案例的流程图。我们提取 1,000 个数据点作为训练集，另外 500 个数据点作为测试集。

首先，我们导入 PyTorch 和 torchvision——后者是 PyTorch 的一个扩展，管理与计算机视觉相关的功能和数据集，如 MNIST。从 torchvision 中，我们导入模块 datasets 和 transforms，它们帮助我们下载和操作 MNIST 数据，分别是：

import torch
from torchvision import datasets, transforms

在第二个数据准备步骤中，我们再次使用将图像转换为 PyTorch 张量的对象（这是 GPyTorch 中实现的 GP 可以处理的数据结构）并规范化像素值。此规范化通过将像素值减去 0.1307（数据的平均值）并将值除以 0.3081（数据的标准差）来完成。这种规范化被认为是 MNIST 数据集的常见做法，有关此步骤的更多详细信息可以在 PyTorch 的官方论坛讨论中找到（mng.bz/BmBr）：

transform = transforms.Compose([
    transforms.ToTensor(),                       ❶
    transforms.Normalize((0.1307,), (0.3081,))   ❷
])

❶ 将数据转换为 PyTorch 张量

❷ 规范化张量

存储在 transform 中的此转换对象现在可以传递给对任何 torchvision 数据集初始化的调用，并且将应用转换（转换为 PyTorch 张量和规范化）到我们的数据上。我们使用此转换对象初始化 MNIST 数据集如下。请注意，我们创建数据集两次，一次将 train 设置为 True 以创建训练集，另一次将 train 设置为 False 以创建测试集：

dataset1 = datasets.MNIST(                                       ❶
    "../data", train=True, download=True, transform=transform    ❶
)                                                                ❶

dataset2 = datasets.MNIST(                                       ❷
    "../data", train=False, download=True, transform=transform   ❷
)                                                                ❷

❶ 下载并导入训练集

❷ 下载并导入测试集

作为数据准备的最后一步，我们从训练集中提取前 1,000 个数据点和测试集中的 500 个点。我们通过从数据集对象 dataset1 和 dataset2 中访问来实现这一点：

使用 data 属性获取特征，即构成每个数据点图像的像素值
使用 targets 属性获取标签，即手写数字的值：

train_x = dataset1.data[:1000, ...].view(1000, -1)
➥.to(torch.float)                   ❶
train_y = dataset1.targets[:1000]    ❶

test_x = dataset2.data[:500, ...].view(500, -1)
➥.to(torch.float)                   ❷
test_y = dataset2.targets[:500]      ❷

❶ 获取训练集中的前 1,000 个点

❷ 获取测试集中的前 500 个点

我们还实现了一个简单的 GP 模型，具有恒定均值和带有输出比例的 RBF 核：

class GPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super().__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) ❸

❶ 一个恒定均值函数

❷ 具有输出比例的 RBF 协方差函数

❸ 以输入 x 的预测制作 MVN 分布

注意 GPyTorch GP 模型的 forward() 方法首次讨论于第 2.4 节。

然后，我们初始化我们的 GP 并在 1,000 个点的训练集上进行训练，使用 Adam 优化器的梯度下降。此代码将优化 GP 的超参数值（例如，均值常量和长度和输出比例），以便我们获得观察到的数据的高边际似然度：

likelihood = gpytorch.likelihoods.GaussianLikelihood()     ❶
model = GPModel(train_x, train_y, likelihood)              ❶

optimizer = torch.optim.Adam(model.parameters(), lr=0.01)  ❷
mll = gpytorch.mlls.ExactMarginalLogLikelihood             ❷
➥(likelihood, model)                                      ❷

model.train()                                              ❸
likelihood.train()                                         ❸

for i in tqdm(range(500)):                                 ❹
    optimizer.zero_grad()                                  ❹

    output = model(train_x)                                ❹
    loss = -mll(output, train_y)                           ❹

    loss.backward()                                        ❹
    optimizer.step()                                       ❹

model.eval()                                               ❺
likelihood.eval()                                          ❺

❶ 声明似然函数和 GP 模型

❷ 声明梯度下降算法和损失函数

❸ 启用训练模式

❹ 运行五百次梯度下降迭代

❺ 启用预测模式

注意参考第 2.3.2 节，了解梯度下降如何优化我们观察到的数据的似然度，即梯度下降如何训练 GP。

最后，为了查看我们的模型在测试集上的表现如何，我们计算 GP 预测值与地面实况（每个数据点的标签值）之间的平均绝对差异。这个指标通常被称为 平均绝对误差。

注意 MNIST 数据集的典型指标是该模型正确预测测试集的百分比（即准确度），这是分类问题的规范。由于我们使用这个数据集来模拟一个回归问题，因此均方误差是合适的。

这是通过将均值预测与存储在test_y中的真实标签进行比较来完成的：

with torch.no_grad():
    mean_preds = model(test_x).mean

print(torch.mean(torch.abs(mean_preds - test_y)))

Output: 2.7021167278289795

这个输出意味着，平均而言，高斯过程对图像中描绘的数字的值的预测误差几乎达到了 3。考虑到这项任务只有 10 个值需要学习，这个表现相当差。这个结果强调了常规高斯过程模型处理结构化数据（如图像）的无能。

13.3 使用神经网络处理复杂的结构化数据

我们所看到的高斯过程表现较差的根本原因是核不具备处理输入数据的复杂结构的装备，从而导致协方差计算不良。特别是，径向基核具有一个简单的形式，只考虑两个输入之间的数字值的差异。在本节中，我们学习如何通过使用神经网络处理结构化数据，然后将处理后的数据馈给高斯过程的均值函数和核来解决这个问题。

13.3.1 为什么使用神经网络进行建模？

我们在本书开始时指出，神经网络在进行昂贵的数据获取时，特别是在进行不确定性校准的预测方面表现不佳。（这是为什么 BayesOpt 中使用高斯过程的全部原因。）然而，神经网络擅长学习复杂结构。这种灵活性是由于神经网络中有多个计算层（具体来说，是矩阵乘法），如图 13.9 所示。

图 13.9 神经网络是一组层计算的集合。通过将多个计算层链接在一起，神经网络可以很好地模拟复杂函数。

在神经网络中，每个层都对应于一个矩阵乘法，其输出然后经过非线性激活函数处理。通过在一次前向传递中将多个这样的层链在一起，可以以灵活的方式处理和操作网络的输入。最终结果是神经网络可以很好地模拟复杂函数。有关神经网络及其用法的详细解释，请参阅 François Chollet 的优秀著作Deep Learning with Python, Second Edition（Manning, 2021）。

神经网络具有的灵活性可以帮助我们解决上一节中描述的问题。如果高斯过程的核，如径向基核，不能很好地处理复杂的数据结构，我们可以让神经网络来处理这项工作，并将处理后的输入仅馈送给高斯过程的核。这个过程在图 13.10 中进行了可视化，其中输入的x首先通过神经网络层，然后再传递给高斯过程的均值函数和核。

图 13.10 结合了神经网络和 GP。神经网络首先处理结构化数据输入x，然后将输出馈送到 GP 的均值函数和核。

尽管最终结果仍然是一个 MVN 分布，但均值函数和核函数的输入现在是由神经网络产生的处理过的输入。这种方法很有前途，因为它具有灵活的建模能力，神经网络将能够从结构化输入数据中提取重要特征（在提供相似性计算信息方面很重要）并将其化简为适合 GP 核的数值。

定义神经网络通常被称为组合模型的特征提取器，因为网络从结构化数据中提取有利于 GP 建模的特征。

通过这种方式，我们可以利用神经网络的灵活学习能力，同时保持使用 GP 进行不确定性校准预测的能力。这是两全其美！此外，训练这个组合模型的过程与训练常规 GP 的过程相同：我们定义我们的损失函数，即负对数似然，然后使用梯度下降来找到最能解释我们的数据的超参数值（通过最小化损失）。现在，我们不仅优化均值常数、长度和输出比例，还要额外优化神经网络的权重。在下一小节中，我们将看到，使用 GPyTorch 实现这个学习过程几乎不需要做任何改动。

注意：这种组合框架是一种动态学习如何处理结构化数据的方法，纯粹来自我们的训练数据集。以前，我们只使用固定的核来处理数据，在多个应用程序中以相同的方式进行处理。在这里，我们“动态地”学习处理我们的输入数据的最佳方式，这对于手头的任务是独一无二的。这是因为神经网络的权重是相对于训练数据进行优化的。

13.3.2 在 GPyTorch 中实现组合模型

最后，我们现在实现这个框架并将其应用于我们的 MNIST 数据集。在这里，定义我们的模型类更加复杂，因为我们需要实现神经网络并将其连接到 GP 模型类中。让我们先解决第一部分——先实现神经网络。我们设计一个简单的神经网络，其架构如图 13.11 所示。此网络具有四个层，节点数分别为 1,000、5,000、50 和 2，如图中所示。这是一个常见的 MNIST 数据集架构。

图 13.11 要实现的神经网络架构。它有四层，并为每个输入数据点产生一个大小为两个的数组。

注意，我们需要关注最后一层（2）的大小，它表示要馈入高斯过程的均值函数和核函数的处理输出的维度。将该层的大小设置为 2，目的是学习存在于二维空间中的图像表示。其他值也可以使用，但为了可视化的目的，我们选择了 2。

我们使用 PyTorch 中的Linear()和ReLU()类实现该体系结构。在这里，我们的网络的每一层都被实现为一个带有相应大小的torch.nn.Linear模块，如图 13.11 所定义的。每个模块还与一个torch.nn.ReLU激活函数模块相耦合，该模块实现了前面提到的非线性变换。这在图 13.12 中得到了说明，其中注释了网络体系结构的每个组件对应于实现它的相应代码。

图 13.12 中实现的神经网络体系结构及其相应的 PyTorch 代码。每个层都使用torch.nn.Linear实现，每个激活函数都使用torch.nn.ReLU实现。

通过使用方便的add_module()方法，我们隐含定义了神经网络模型的forward()方法的逻辑。接下来，我们使用LargeFeatureExtractor类实现模型。该类将其输入x依次通过我们在__init__()方法中实现的层中，该方法接收data_dim，即输入数据的维数。在我们的情况下，该数字为 28×28=784，我们使用train_x.size(-1)进行计算：

data_dim = train_x.size(-1)                            ❶

class LargeFeatureExtractor(torch.nn.Sequential):
    def __init__(self, data_dim):
        super(LargeFeatureExtractor, self).__init__()

        self.add_module('linear1', torch.nn.Linear
        ➥(data_dim, 1000))                            ❷
        self.add_module('relu1', torch.nn.ReLU())      ❷
        self.add_module('linear2', torch.nn.Linear
        ➥(1000, 500))                                 ❸
        self.add_module('relu2', torch.nn.ReLU())      ❸

        self.add_module('linear3', torch.nn.Linear
        ➥(500, 50))                                   ❹
        self.add_module('relu3', torch.nn.ReLU())      ❹

        self.add_module('linear4', torch.nn.Linear
        ➥(50, 2))                                     ❺

feature_extractor = LargeFeatureExtractor(data_dim)    ❻

❶ 数据的维度

❷ 网络的第一层

❸ 第二层

❹ 第三层

❺ 第四层

❻ 初始化网络

接下来，我们讨论组合模型——一种利用我们刚刚初始化的神经网络特征提取器feature_extractor的高斯过程模型类。我们首先实现它的__init__()方法，该方法由几个组件组成：

协方差模块被包装在gpytorch.kernels.GridInterpolationKernel对象中，为我们的中等大小训练集（1,000 个点）提供计算速度加速。我们声明输入数据的维数为二，因为这是特征提取器生成的输出的维度。
特征提取器本身就是我们之前声明的feature_extractor变量。
如果神经网络的权重初始化得不好，特征提取器的输出值可能会取极端值（负无穷或正无穷）。为解决这个问题，我们使用gpytorch.utils.grid.ScaleToBounds模块将这些输出值缩放到-1 和 1 之间的范围内。

__init__()方法的实现如下：

class GPRegressionModel(gpytorch.models.ExactGP):
  def __init__(self, train_x, train_y, likelihood):
      super(GPRegressionModel, self).__init__(train_x, train_y, likelihood)

      self.mean_module = gpytorch.means.ConstantMean()

      self.covar_module = gpytorch.kernels
      ➥.GridInterpolationKernel(                    ❶
          gpytorch.kernels.ScaleKernel(              ❶
              gpytorch.kernels.RBFKernel             ❶
              ➥(ard_num_dims=2)                     ❶
          ),                                         ❶
          num_dims=2,                                ❶
          grid_size=100                              ❶
      )                                              ❶

      self.feature_extractor = feature_extractor     ❷

      self.scale_to_bounds = gpytorch.utils.grid
      ➥.ScaleToBounds(-1., 1.)                      ❸

❶ 具有两个维度的 RBF 核函数，具有计算速度加速

❷ 神经网络特征提取器

❸ 一个用于将神经网络的输出缩放至合理值的模块

在我们的forward()方法中，我们将所有这些组件结合在一起。首先，我们使用我们的神经网络特征提取器处理输入。然后，我们将处理后的输入馈送到我们的 GP 模型的平均值和协方差模块中。最后，我们仍然得到一个 MVN 分布，就是forward()方法返回的结果：

class GPRegressionModel(gpytorch.models.ExactGP):
  def forward(self, x):
    projected_x = self.feature_extractor(x)            ❶
    projected_x = self.scale_to_bounds(projected_x)    ❶

    mean_x = self.mean_module(projected_x)             ❷
    covar_x = self.covar_module(projected_x)           ❷
    return gpytorch.distributions.MultivariateNormal   ❷
    ➥(mean_x, covar_x)                                ❷

❶ 缩放后的神经网络特征提取器的输出

❷ 从处理后的输入创建一个 MVN 分布对象

最后，为了使用梯度下降训练这个组合模型，我们声明了以下对象。在这里，除了常规的 GP 超参数，如平均常数和长度和输出尺度，我们还想优化神经网络特征提取器的权重，这些权重存储在model.feature_extractor.parameters()中：

likelihood = gpytorch.likelihoods.GaussianLikelihood()             ❶
model = GPRegressionModel(train_x, train_y, likelihood)            ❶
mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)  ❶

optimizer = torch.optim.Adam([
    {'params': model.feature_extractor.parameters()},              ❷
    {'params': model.covar_module.parameters()},                   ❷
    {'params': model.mean_module.parameters()},                    ❷
    {'params': model.likelihood.parameters()},                     ❷
], lr=0.01)

❶ 同之前一样，似然函数、GP 模型和损失函数保持不变。

❷ 现在，梯度下降优化器 Adam 需要优化特征提取器的权重和 GP 的超参数。

现在我们可以像之前一样运行梯度下降：

model.train()                      ❶
likelihood.train()                 ❶

for i in tqdm(range(500)):         ❷
    optimizer.zero_grad()          ❷

    output = model(train_x)        ❷
    loss = -mll(output, train_y)   ❷

    loss.backward()                ❷
    optimizer.step()               ❷

model.eval()                       ❸
likelihood.eval()                  ❸

❶ 启用训练模式

❷ 运行 500 次梯度下降迭代

❸ 启用预测模式

注意：提醒一下，当训练 GP 模型时，我们需要同时为模型和似然函数启用训练模式（使用model.train()和likelihood.train()）。在训练之后并在进行预测之前，我们需要启用预测模式（使用model.eval()和likelihood.eval()）。

现在，我们已经训练了与神经网络特征提取器结合的 GP 模型。在使用该模型对测试集进行预测之前，我们可以查看模型的内部，看看神经网络特征提取器是否学会了很好地处理我们的数据。请记住，每个图像都被特征提取器转化为一个二元数组。因此，我们可以将训练数据通过该特征提取器，并使用散点图可视化输出。

注意：训练这个组合模型比训练普通的 GP 模型需要更多时间。这是因为我们现在要优化的参数更多。然而，正如我们马上会看到的那样，这个成本是非常值得的，因为我们获得了更高的性能提升。

在这个散点图中，如果我们看到相同标签的点（即，描绘相同数字的图像）聚在一起，这将表明特征提取器能够有效地从数据中学习。同样，我们通过将训练数据通过特征提取器的方式进行处理来做到这一点，这与模型类的forward()方法中的数据处理方式相同：

with torch.no_grad():
    extracted_features = model.feature_extractor(train_x)
    extracted_features = model.scale_to_bounds(extracted_features)

在这里，extracted_features是一个大小为 1,000x2 的 PyTorch 张量，存储了我们训练集中 1,000 个数据点的二维提取特征。为了在散点图中可视化这个张量，我们使用 Matplotlib 库的plt.scatter()方法，确保每个标签对应一个颜色：

for label in range(10):
    mask = train_y == label           ❶

    plt.scatter(                      ❷
        extracted_features[mask, 0],  ❷
        extracted_features[mask, 1],  ❷
        c=train_y[mask],              ❷
        vmin=0,                       ❷
        vmax=9,                       ❷
        label=label,                  ❷
    )                                 ❷

❶ 过滤具有特定标签的数据点

❷ 为当前数据点创建一个散点图，它们具有相同的颜色

此代码生成图 13.13，尽管你的结果可能会有所不同，这取决于库版本和代码运行的系统。正如我们所预期的，相同标签的数据点围绕在一起。这意味着我们的神经网络特征提取器成功地将具有相同标签的点分组在一起。经过网络处理后，具有相同标签的两个图像变成了二维平面上彼此靠近的两个点，如果由 RBF 核计算，则它们将具有高的协方差。这正是我们希望我们的特征提取器帮助我们做的事情！

图 13.13 由神经网络从 MNIST 数据集中提取的特征。不仅相同标签的数据点聚集在一起，而且在图中还存在一个标签梯度：从底部到顶部，标签值逐渐增加。

图 13.13 另一个有趣的方面是，与标签值相关的梯度明显：从底部到顶部的聚类，相应标签的值从 0 逐渐增加到 9。这是特征提取器中很好的特性，因为它表明我们的模型已经找到了一种平滑的 MNIST 图像表示，符合标签值。

例如，考虑图 13.14 中的比较，左侧面板显示图 13.13，右侧面板显示相同散点图标签的随机交换，使特征变得“粗糙”。所谓“粗糙”，是指标签值在不规律地跳动：底部聚类包含 0，中间某些聚类对应于 7 和 9，顶部聚类包含 5。换句话说，具有粗糙特征的标签趋势不是单调的，这使得训练 GP 更加困难。

图 13.14 图 13.13 中提取的平滑特征与标签随机交换的比较，使特征变得不那么平滑。平滑的特征比粗糙的特征更容易通过 GP 学习。

看起来神经网络在从图像中提取有用特征方面表现不错。为了确定这是否确实导致更好的预测性能，我们再次计算平均绝对误差（MAE）：

with torch.no_grad():
    mean_preds = model(test_x).mean

print(torch.mean(torch.abs(mean_preds - test_y)))

Output: 0.8524129986763

这个结果告诉我们，平均而言，我们的预测偏差为 0.85；这是对前一节中我们拥有的普通 GP 的显着改进，其 MAE 大约为 2.7。这种改进说明了联合模型的卓越性能，这源于神经网络灵活的建模能力。

正如我们在开始时所说，这个框架不仅适用于手写数字，还适用于各种类型的结构化数据，神经网络可以从中学习，包括其他类型的图像和图形结构，如分子和蛋白质。我们所要做的就是定义一个合适的 DL 架构，从这些结构化数据中提取特征，然后将这些特征传递给 GP 的均值函数和核函数。

这就结束了第十二章。在本章中，我们了解了从结构化数据中学习的困难，例如图像，在这些数据中，常见的核无法有效地计算数据点之间的协方差。通过在 GP 前面附加一个神经网络特征提取器，我们学会将这些结构化数据转换成 GP 的核函数可以处理的形式。最终结果是一个结合模型，可以灵活地从结构化数据中学习，但仍然产生具有不确定性量化的概率预测。

总结

结构化数据是指其特征需要满足约束条件的数据，例如必须是整数或者非负数，并且不能被视为连续的实值数据。例如，常见应用程序中的数据，如计算机视觉中的图像和药物发现中的蛋白质结构。
结构化数据对于 GP 的常见核构成挑战。这是因为这些核只考虑输入数据的数值，这可能是不良的预测特征。
使用错误特征计算协方差的核可能会导致生成的 GP 的预测质量低下。使用错误特征在结构化数据中特别常见。
对于图像数据特别是，像素的原始值不是一个信息量丰富的特征。使用原始像素值计算协方差的核可能导致低质量的 GP。
由于具有多层非线性计算，神经网络能够有效地学习复杂函数，并且可以从结构化数据中提取特征。通过使用神经网络从结构化数据中提取连续的实值特征，GP 仍然可以有效地学习。
在将神经网络与 GP 结合时，我们动态学习一种处理问题的数据方式。这种灵活性使得该模型可以推广到许多种结构化数据。
在将输出缩放到小范围之前，将神经网络的输出传递给 GP 是很重要的。通过这样做，我们避免了由于神经网络特征提取器初始化不良而导致的极端值。
从神经网络特征提取器学习到的表示对标签具有平滑的梯度。这种平滑的梯度使得提取的特征更容易通过 GP 学习。