贝叶斯优化循环

第 1.3 节指出了贝叶斯优化的两个要素：一是能够刻画自身不确定性的目标函数模型，二是把这种认识转化为下一次实验的规则。第二部分构建了第一个要素。以已有的评估为条件，高斯过程给出后验均值和后验标准差：前者描述函数大致的形状，后者描述模型在每个输入处的不确定程度（第 8 章）。本章把这一模型放进循环之中。

循环本身十分简短，一张卡片即可写下：拟合模型；选出再评估一次看来最有价值的输入；评估该输入；把结果加入数据；重复。本书其余部分的大多数思想，都是对这四步中某一步的改进。循环的关键在第二步，因为“最有价值”没有唯一正确的答案。一个输入值得评估，可能是因为模型预期该处取值较高，也可能是因为模型对它所知甚少，高值或许就藏在那里。选择下一个输入的每条规则，都要在这两个理由之间取得某种平衡；本章让读者亲自体会，这种平衡偏向任何一方会带来什么后果。

本章先精确地表述问题，再把循环写成算法，并在贯穿全书的示例目标函数上观察它的运行。本章中间部分讨论探索与利用之间的权衡，最后讨论两个实际问题：模型尚无任何依据时如何选取最初的几个点，以及这一方法从何而来。之后，第 12 章将逐一推导选择下一个输入的各种规则。

11.1 问题 #

设 $f$ 是从输入定义域 $\X$ 到实数的函数，目标是找到使 $f$ 最大的输入：

\vx^\star \in \argmax_{\vx \in \X} f(\vx).

如此表述，这一问题与其他优化问题并无二致。贝叶斯优化的不同之处在于它的适用条件，Frazier（2018）明确列出了这些条件。

每次评估都很昂贵。在一个输入处评估 $f$ 需要几分钟乃至几小时，或者需要花钱，或者需要人付出一定的代价。评估次数有限，通常至多几百次。允许的评估次数称为预算（budget），记作 $N$ 。
函数是黑箱。 $f$ 可以在任意位置评估，但没有解析式，也没有凸性、线性等可供专门方法利用的已知结构。
没有导数。一次评估只返回一个函数值，不提供其他信息，因此无法使用梯度下降及其同类方法。
定义域简单且规模不大。 $\X$ 通常是一个箱形区域，即 $d$ 个输入各有一个取值范围，可缩放为 $[0, 1]^d$ 。大多数成功的应用满足 $d \le 20$ （Frazier，2018）；超出这一范围时情况如何变化，见第 14.6 节。
函数较为光滑。相近的输入往往给出相近的值。高斯过程先验编码的正是这一假设；没有这一假设，有限次评估无论多少，都无法提供评估点之间任何输入的信息。

两个已发表的问题可以说明这些条件在实践中的具体情形。Snoek 等人（2012）为一个图像分类卷积神经网络调节了 9 个超参数，包括学习率、训练轮数和 4 个权重惩罚项。每次评估都是一次完整的训练，找到的最佳设定在 CIFAR-10 基准上的测试误差为 14.98%，而人类专家调出的设定为 18%。Shields 等人（2021）优化化学反应，每次评估都是一次实验室实验。在一项以在线游戏形式进行的基准测试中，与化学和工程领域的专家相比，优化器平均所需的实验更少，结果对初始数据的依赖也更小。第 22 章和第 23 章完整演示了这两类问题。

评估也可能带有噪声。用不同的随机种子把同一个神经网络训练两次，得到的准确率不同；同一个人给同一杯咖啡评两次分，分数也不同。此时观测值为 $y = f(\vx) + \varepsilon$ ，其中 $\varepsilon$ 为噪声，即第 8.3 节中的模型；任务仍是最大化 $f$ ，而不是带噪声的 $y$ 。

预算既然有限，方法就不必精确地找到 $\vx^\star$ 。 $N$ 次评估之后，方法必须返回一个推荐点（recommendation） $\hat\vx_N$ ，评判标准是 $f(\hat\vx_N)$ 与可达到的最优值 $f^\star = f(\vx^\star)$ 的接近程度。差距 $f^\star - f(\hat\vx_N)$ 称为简单遗憾（simple regret），第 13.1 节将仔细定义它及相关概念。此处只需明白两点：这一差距越小越好；有些评估本身得分很低，却提供了大量信息，只要它们最终带来更好的推荐点，就不算浪费。

本部分的图都使用 $[0, 1]$ 上的同一个示例目标函数，以青绿色虚线表示。它在 $x = 0.23$ 附近有一个宽峰，高约 $0.53$ ；在 $x = 0.73$ 附近有一个更窄也更高的峰，高 $0.82$ ；两峰坐落在一条缓缓下降、略有起伏的基线上。这一形状是有意设置的陷阱：先找到宽峰的搜索很容易停留在那里，再也发现不了更高的峰。图中之所以能画出这条虚线，是因为作图时已知其解析式。循环却始终看不到它，只能看到自己选择评估的那些输入处的函数值。

第 11.1 节引用的文献 3

Frazier（2018）A Tutorial on Bayesian Optimization
Snoek 等人（2012）Practical Bayesian Optimization of Machine Learning Algorithms
Shields 等人（2021）Bayesian reaction optimization as a tool for chemical synthesis

11.2 循环 #

贝叶斯优化的核心思想，是把一个困难的问题换成一系列容易的问题。直接最大化 $f$ 很难，因为每次评估都很昂贵。因此，方法在每一步都根据模型构造一个廉价的函数，转而最大化这个函数。廉价函数为每个候选输入打分，分数表示下一次在该处评估 $f$ 的用处有多大。求它的最大值可能需要数千次求值，但每次只需几微秒，而不是几小时。

记 $\D_n = \{(\vx_1, y_1), \dots, (\vx_n, y_n)\}$ 为 $n$ 次评估后的数据， $\mu_n(\vx)$ 和 $\sigma_n(\vx)$ 为给定 $\D_n$ 时 $f(\vx)$ 的后验均值和后验标准差，由式（8.6）计算。这里的下标 $n$ 表示评估次数， $\sigma_n(\vx)$ 总是连同自变量一起书写，以区别于第 8.3 节中的噪声标准差 $\sigma_n$ ：后者不带自变量，下标代表噪声（noise）。

定义 11.1 采集函数

采集函数（acquisition function）是由给定 $\D_n$ 的后验计算得到的函数 $a_n : \X \to \R$ ，它为每个输入打分，衡量下一次在该处评估 $f$ 的用处。循环在采集函数取最大值处评估 $f$ ：

\vx_{n+1} \in \argmax_{\vx \in \X} a_n(\vx).

大多数采集函数只通过 $\mu_n(\vx)$ 和 $\sigma_n(\vx)$ 依赖于 $\vx$ ，有时还用到迄今观测到的最优值，因此在任意位置求值的代价都只是一次后验预测。两个要素齐备之后，整个方法就是一个循环。

算法 11.1 贝叶斯优化

输入：定义域 $\X$ ；预算 $N$ ；高斯过程先验（核函数与均值）；采集函数 $a$ ；初始设计的规模 $n_0 < N$ 。

初始设计。不借助模型选出 $\vx_1, \dots, \vx_{n_0}$ （第 11.4 节），在各点评估 $f$ ，得到 $\D_{n_0}$ 。
对 $n = n_0, n_0 + 1, \dots, N - 1$ $n = n_{0}, n_{0} + 1, \dots, N - 1$ ：
1. 拟合。以 $\D_n$ 为条件更新高斯过程，得到 $\mu_n(\vx)$ 和 $\sigma_n(\vx)$ ；如有需要，重新拟合核函数的超参数（第 9.4 节）。
2. 决策。用普通的数值优化器求出 $\vx_{n+1} \in \argmax_{\vx} a_n(\vx)$ （第 12.9 节）。
3. 查询。在 $\vx_{n+1}$ 处评估目标函数。
4. 回答。得到 $y_{n+1} = f(\vx_{n+1}) + \varepsilon_{n+1}$ 。
5. 更新。令 $\D_{n+1} = \D_n \cup \{(\vx_{n+1}, y_{n+1})\}$ 。
返回推荐点 $\hat\vx_N$ （第 11.2.3 节）。

这一结构与 Frazier（2018）的算法 1、Garnett（2023）的算法 1.1 相同。

查询、回答与更新在这里只是一次函数调用，算法却把它们分为三步，这样区分的好处在后文才会显现。在偏好贝叶斯优化（第 19 章）中，查询变为展示给人的一对输入，回答变为表示其偏好哪一个的单个比特；由于一次比较并不是带噪声的数值，更新也不能再使用高斯过程的闭式公式（第 18 章）。这些步骤外面的循环则保持不变。

循环中三类工作的代价相差悬殊，方法的设计正是由此而来。评估 $f$ 是昂贵的一步，其余一切安排都因它而起。拟合模型需要做 Cholesky 分解，代价为 $O(n^3)$ （第 8.4 节），典型运行只有几百个观测，用时远不到一秒。最大化采集函数需要多次计算 $\mu_n(\vx)$ 和 $\sigma_n(\vx)$ ，分解完成后每次代价为 $O(n^2)$ 。与耗时一小时的评估相比，这些计算都可以忽略不计，因此值得花费大量计算来选好每一次评估。

11.2.1 第一个采集函数 #

运行循环需要一个具体的采集函数。标准的采集函数留待第 12 章推导，这里用一条简单明了的规则即可。以乐观估计为每个输入打分，即后验均值加上后验标准差的若干倍：

a_n(\vx) = \mu_n(\vx) + \beta^{1/2}\, \sigma_n(\vx),

其中 $\beta \ge 0$ 是选定的常数。取 $\beta^{1/2} = 2$ 时，这一分数接近图中 95% 可信带的上边缘（该带在均值之上延伸 $1.96$ 个标准差），因此这条规则选出的是合理范围内最好情形最高的输入。这条规则称为上置信界（upper confidence bound，UCB），第 12.4 节将再次讨论它，并介绍 $\beta$ 应当如何随时间增大的理论。

式（11.3）的两项作用方向不同，这正是第 11.3 节的主题。均值项偏好模型已认为较好的输入，标准差项偏好模型所知甚少的输入。权重 $\beta^{1/2}$ 规定了两者之间的兑换率：一个单位的不确定性折合多少个单位的期望值。

11.2.2 在示例目标函数上运行循环 #

图 11.1 在示例目标函数上运行算法 11.1，采集函数取式（11.3）， $\beta^{1/2} = 2$ 。运行从两个随机输入开始，其中一个恰好落在宽峰顶部 $x \approx 0.23$ 处，因此循环一开始就落入了目标函数有意设下的陷阱。

图 11.1 示例目标函数上的贝叶斯优化。时间轴的每一步对应算法 11.1 的一次迭代。上：隐藏的目标函数（虚线）；根据迄今评估（圆点，最新一次加圈）得到的后验均值与 95% 可信带；下一次查询（橙色线）。下：采集函数及其最大值。图的初始状态为式（11.3）的上置信界，“上置信界权重 √β”滑块设定标准差的乘数

\beta^{1/2}

。其他选项在第 12 章中推导。核函数及其长度尺度固定不变、不经拟合，以便各图之间可以比较。

逐步查看时间轴，可以看出大多数运行遵循的规律。

最初的查询落在可信带最宽处。只有两个观测时，式（11.3）几乎处处由标准差项主导，因此循环的头几步都用在远离这两个观测的输入上，包括定义域的两端。这些评估并未浪费：每一次都使附近的可信带收窄，从而排除一片区域。

乐观找到了高峰。大约在第五或第六步，可信带剩下的宽阔部分正位于高峰上方，在那里的一次评估返回的值远高于此前见过的任何值。均值随之跃升，此后均值项便把查询拉回这一区域。

运行末段做精细调整。一旦可信带除最佳区域附近外处处变窄，查询便聚集在 $x \approx 0.72$ 周围，找到的最优值逼近真实最大值，采集函数变成一个窄尖峰。

换一条规则，换一个起点。将采集函数切换为“期望改进”“改进概率”或“Thompson 采样”（三者均在第 12 章中推导），并按“新一轮运行”更换初始点。细节会变，但先大范围搜索、后精细调整的模式不变。切换回“上置信界”并把“上置信界权重 √β”滑块设为 0，循环便再也不会离开宽峰。

三维中的同一个循环。算法 11.1 并不要求输入是单个数。下图在 Hartmann 函数上运行该算法；这是一个标准的三维测试问题，有一个全局最大值和若干局部极大值。后验已无法用曲线表示，因此图中给出单位立方体内的已评估点（可拖动旋转），以及穿过当前最佳点的三条模型切片，每个输入方向一条。

图 11.2 三维 Hartmann 函数上的循环：先取 5 个初始点，再用期望改进，共 25 次评估。左：单位立方体中的已评估点，值越高，点越大、越不透明，并以落到底面的垂线显示深度；最新的点加圈，最佳点为橙色，星号为真正的最大值点。右：穿过最佳点的模型切片，每个输入一条，显示后验均值与 95% 可信带（蓝色），以及同一直线上的真实函数（虚线）。下：找到的最优值。可以逐步查看这次运行、旋转立方体，或按“新一轮运行”尝试其他初始设计。

有两点值得留意。其一，早期的点散布在整个立方体中，后期的点则聚集起来，与一维中先广后窄的模式相同。其二，默认的运行在大部分预算内停滞不前：从第 9 次到第 23 次评估，最佳点一直是 $(0.83, 0.56, 0.87)$ ，值为 3.59，而最大值为 3.86。退回到第 23 次评估，沿 $x_1$ 的切片显示真实函数向 $x_1 \approx 0.1$ 方向上升，那里恰恰是模型的可信带仍然很宽的地方。第 24 次评估朝这个方向移动，在 $(0.33, 0.51, 0.85)$ 处达到 3.76，但仍未到达位于 $x_1 = 0.11$ 的最大值点。更多的预算、不同的起点或更偏重探索的规则能够找到全局最大值；要看出一次运行尚未找到它，查看切片往往是最快的办法。

预算用完时，循环必须给出一个输入。常见的选择有两种：观测值最好的已评估输入，或后验均值最高的输入（Frazier，2018）。评估精确时，第一种选择是稳妥的：它的值已经测得，确定无疑。

评估带有噪声时，最好的观测值是有偏的依据。在许多带噪声的测量中，最大的那个往往恰好带有正的噪声，因此产生它的输入很可能不如测量值显示的那样好。推荐后验均值最高的输入（无论是在已评估的输入中，还是在整个定义域上），利用的是某个输入附近的全部评估，而不是单独一次侥幸的评估。第 14.2 节将再次讨论这一选择；第 12.6 节则表明，推荐方式不同，有原则的采集函数也随之不同。

代码实现 NumPy

import numpy as np
# rbf() 和 gp_posterior() 来自高斯过程那一章

def f(x):  # 示例目标函数；假装每次调用要花一小时
    return (0.62 * np.exp(-(x - 0.25) ** 2 / (2 * 0.1**2))
            + np.exp(-(x - 0.73) ** 2 / (2 * 0.055**2))
            + 0.1 * np.sin(11 * x + 0.6) - 0.35 * x)

rng = np.random.default_rng(0)
grid = np.linspace(0, 1, 501)      # 内层搜索的候选点
X = rng.uniform(0, 1, size=2)      # 初始设计
Y = f(X)
for n in range(10):
    m = Y.mean()                   # 常数先验均值
    mean, var = gp_posterior(X, Y - m, grid, ell=0.08)
    sd = np.sqrt(np.maximum(var, 0.0))
    acq = mean + m + 2.0 * sd      # 均值 + sqrt(beta) * 标准差
    x_next = grid[np.argmax(acq)]  # 决策
    X = np.append(X, x_next)       # 查询……
    Y = np.append(Y, f(x_next))    # ……回答，并更新

print("recommend x =", X[np.argmax(Y)], "with f =", Y.max())

函数 gp_posterior 即第 8.4 节中的实现。在网格上取最大值只适用于一维，更高维中的替代方法见第 12.9 节。在这一随机种子下，循环在第 9 次查询时找到高峰，并推荐 $x = 0.72$ ，该处 $f \approx 0.81$ 。附录 C 把这段示意代码扩展到期望改进。

第 11.2 节引用的文献 2

Frazier（2018）A Tutorial on Bayesian Optimization
Garnett（2023）Bayesian Optimization

11.3 探索与利用 #

每个采集函数都必须回答本章开头提出的问题：一个输入值得评估，是因为模型预期该处取值较高，还是因为模型对该处并不了解？这两种回答各有名称。利用（exploitation）指在后验均值高的地方评估，以精细调整已经看来不错的结果。探索（exploration）指在后验标准差高的地方评估，以了解模型所知甚少的区域。这两个术语来自赌博机问题的研究（第 13.2 节）：赌徒必须在两台赌博机之间做出选择，一台迄今回报最好，另一台尝试次数太少，尚无法判断。

两种纯粹的策略都行不通，示例目标函数可以说明原因。

纯粹的利用会陷入停滞。在式（11.3）中令 $\beta = 0$ ，循环便总在后验均值最高处评估。假设最初的几次评估落在宽峰上。此时均值在其中最好的点附近最高，于是下一次评估落在附近，确认这一区域不错，并进一步抬高那里的均值。循环沿宽峰爬升，到达约 $0.53$ 的峰顶后便停在那里。它的规则不会把它引向定义域的右半边；那里的均值较低，仅仅是因为尚无任何观测。观测精确时，它甚至会一次又一次地评估同一个输入，每次都一无所获（习题 11.2）。

纯粹的探索始终无法聚焦。把 $\beta$ 设得很大，使标准差项占主导。循环便总在模型最不确定处评估，在区间上，这意味着填补已有评估之间最大的空隙，结果近似于逐点搭建的网格。它最终会落到高峰附近，但对这一区域的关注并不多于右端那片较差的区域，因此找到的最优值受限于预算用完时网格的疏密。

纯粹的探索随维度增长表现也很差。假设一次评估能让模型在距其 $0.2$ 以内（约一个长度尺度）有把握。半径为 $0.2$ 的球覆盖单位区间的 40%、单位正方形的 13%、单位立方体的 3.4%，但只覆盖 6 维单位立方体的 0.033%。以这种方式覆盖 6 维空间至少需要 3,000 次评估，覆盖 Snoek 等人（2012）所调网络的 9 维空间至少需要 590,000 次。由于这些球彼此重叠，还会伸出立方体之外，实际数字还要更大。没有哪种预算承受得起。维度一多，探索就必须有所取舍：只能在仍可能出现高值的地方减少不确定性。

两个极端之间，存在一段表现良好的权重范围。这一范围有多宽？两个极端又失败到什么程度？图 11.3 通过多次运行循环来回答。

图 11.3 式（11.3）中的探索权重

\beta^{1/2}

，由滑块设定。上：一次共 12 次评估的运行（两个随机初始点，空心；其后 10 个点由规则选出）及最终的后验，找到的最优值加圈标出。左下：该次运行中每次评估后找到的最优值。右下：对 0 至 8 之间的每个权重，从 32 个随机初始设计重复同一实验，显示真实最大值与找到的最优值之间的平均差距；橙色标记为当前权重。图的初始状态为贪心规则，即

\beta = 0

。目标函数、核函数与预算仅作示意，重要的是曲线的形状而非具体数值。

借助这张图做以下四个实验，可以具体地看到这种权衡。

从贪心开始。权重为 0 时，上方面板中的运行沿宽峰爬升，停在 $0.53$ 。32 个起点中，只有 12 个到达距最大值 $0.05$ 以内；其余的始终没有离开起步时所在的峰，平均差距为 $0.18$ 。

增加少许乐观。权重在 1 至 2 之间时，全部 32 个起点都到达距最大值 $0.05$ 以内，平均差距至多为 $0.003$ 。即使权重只取 $0.5$ ，帮助也很大：32 个起点中有 24 个成功。

探索过度。权重为 8 时，上方面板显示评估几乎均匀地散布在整个定义域上，这次运行的最优值为 $0.62$ 。在 32 个起点上，平均差距升至约 $0.07$ 。探索者确实找到了高峰的邻域，但预算分散在各处，很少有评估靠近峰顶本身。

更换起点。按几次“换一个起点”。对某些初始设计，连贪心规则也能成功，因为两个随机点之一恰好落在高峰的山坡上。初始设计中的运气能让任何规则在单次运行中显得不错，因此比较优化器时要对多次运行取平均（第 31.4 节）。

右下面板的形状在这一领域中反复出现：两种失效模式之间的一道谷。在这个问题上，谷底很宽，权重无须精细调节，但两个极端的代价都很大。谷的位置取决于预算。只有剩余的评估足以利用探索的发现，探索才有回报，因此预算短时宜取较小的权重，预算长时可以容忍较大的权重（推断）。

要点采集函数为不确定性定价

每个采集函数都在期望值与不确定性之间权衡。上置信界通过显式的兑换率 $\beta^{1/2}$ 实现这种权衡；第 12 章中的规则则从对评估目的的刻画出发，推导出这一兑换率。

式（11.3）中的加权和已经做到了固定安排做不到的事。像“前一半预算随机探索，然后转为利用”这样的安排，会把探索分摊到各处，包括模型已知较差的区域。上置信界只在不确定性与可能的高值同时存在之处探索：如果一个输入的均值加两个标准差仍低于见过的最优值，那么无论它多么不确定，都不会被选中。并非所有不确定性都值得减少，值得减少的只是可能改变最终推荐结果的那部分不确定性。第 12 章中那些更有原则的采集函数正是从这一观察出发，直接追问一次评估预期能把结果改进多少。

第 11.3 节引用的文献 1

Snoek 等人（2012）Practical Bayesian Optimization of Machine Learning Algorithms

11.4 启动循环 #

模型需要先有数据，才能给出有用的预测；而最初的几个点不宜交给采集函数选择，理由有二（Garnett，2023，第 9.3 节）。

第一个理由是，采集函数此时无所依据。在没有任何数据时，若高斯过程采用常数先验均值和平稳核（协方差只取决于输入之间距离的核函数），则它在每个输入处给出相同的均值和相同的标准差。由此构造的任何采集函数都是常数，任何位置都是它的最大值点（习题 11.1）。

第二个理由与模型的超参数有关。长度尺度、信号幅度和噪声水平通常由数据拟合得到（第 9.4 节），两三个点不足以确定它们。长度尺度拟合不当，模型要么在点与点之间过度自信，要么处处缺乏信息量，据此做出的早期决策可能把搜索引向错误的方向。先用几个不借助模型选出的点，拟合才有可用的数据。

因此，循环从包含 $n_0$ 个点的初始设计（initial design）开始，这些点由一条与目标函数无关的规则选出。其目标是覆盖：各点应散布在定义域上，不留下大片未经考察的区域。常见的布点方法有四种。

网格对每个输入取 $k$ 个值，评估所有组合，在 $d$ 维中共有 $k^d$ 个点。网格容易描述，但有一种特定的浪费：每个输入只用到 $k$ 个不同的值。如果目标函数实际上主要依赖于某一个输入（这很常见），网格花费 $k^d$ 次评估，对这个输入却只了解到 $k$ 个值。点数本身也增长很快：在 6 维中，每个输入只取三个值的网格就需要 $3^6 = 729$ 次评估；在 Snoek 等人（2012）所调网络的 9 维中，则需要 $3^9 = 19{,}683$ 次。Bergstra 与 Bengio（2012）发现，在他们研究的大多数数据集上，学习算法的超参数中真正重要的只有少数几个，而且不同数据集上重要的超参数各不相同，因此网格不宜作为默认选择。

均匀随机采样使每个点在每个输入上都取不同的值，而且无须规划。它的弱点是扎堆：有些点挤在一起，别处却留下空隙。把一个输入的取值范围等分为 $n$ 格， $n$ 个随机点平均留下比例为 $(1 - 1/n)^n$ 的空格， $n$ 较大时约为 37%。

拉丁超立方同时消除沿各坐标轴的扎堆（McKay 等，1979）。把每个输入的取值范围等分为 $n$ 格，放置 $n$ 个点，使每个输入的每一格恰好含有一个点：对每个输入独立地随机排列各格，把第 $i$ 个点放入排列中的第 $i$ 格，格内位置随机。这一名称来自拉丁方，即每个符号在每行、每列中都恰好出现一次的方阵。拉丁超立方保证对每个输入分别覆盖，却不保证覆盖整个空间：所有点都落在对角线上，同样满足定义。一种补救办法是抽取许多个拉丁超立方，保留最近两点间距离最大的那一个。

Sobol 序列是确定性的（Sobol'，1967）。它的构造使每个新点落入先前各点留下的最大空隙，这种性质称为低差异（low discrepancy）：定义域中任一箱形子区域所含的点数都接近其应得的份额。在下图的二维版本中，序列的前 $2^k$ 个点恰好在每个输入的 $2^k$ 格中各放一个点，与拉丁超立方相同。与拉丁超立方不同的是，序列可以延长。向拉丁超立方添加点会破坏每格一点的性质，而继续添加 Sobol 序列的后续点，到下一个 2 的幂时又会恢复这一性质。

图 11.4 单位正方形中的初始设计。底部和左侧的条带显示各点在单个输入上的投影，该输入的取值范围被等分为与点数相同的格数；阴影格中没有点。橙色线段连接距离最近的两点。可以切换设计、改变点数，并按“重新抽取”得到另一个随机实例。

主要结论体现在图 11.4 的投影条带上。

网格。16 个点时，网格是 $4 \times 4$ 的阵列，每个输入的 16 格中有 12 格为空。如果只有第一个输入重要，这 16 次评估只相当于 4 次。

随机。16 个随机点通常在每个输入上留下五六个空格，接近上文预测的 37%；距离最近的一对点往往比其他设计中任何一对点都近得多。按“重新抽取”，可以看到图形变化有多大。

拉丁超立方。由构造方式可知，不会出现空格。最近的一对点通常比随机设计中相距更远，但并非总是如此：反复按“重新抽取”，直到出现两点几乎相接的情形。

Sobol。16 个或 32 个点时没有空格。把滑块移到 20，第二个输入就有四格变空；在 16 至 32 之间的点数中，只有 24 能填满每一格。因此 Sobol 设计通常按 2 的幂取点。

初始点的数量本身也需要权衡，因为设计中的每个点都不是由模型选出的。计算机实验设计中有一条常用规则：每个输入维度取 10 个点，Jones 等人（1998）采用了这条规则，Loeppky 等人（2009）为它提供了理由和证据。他们的准则是高斯过程在整个定义域上预测函数的准确程度，这超出了优化的需要，因为优化器只需在顶部附近预测准确。在 5 维中预算为 30 次评估时，这条规则要求 50 个点，因此预算较小时初始设计也应更小，把其余的探索交给采集函数（推断）。本章的一维图都从两个随机点开始。

第 11.4 节引用的文献 7

Garnett（2023）Bayesian Optimization
Snoek 等人（2012）Practical Bayesian Optimization of Machine Learning Algorithms
Bergstra 与 Bengio（2012）Random Search for Hyper-Parameter Optimization
McKay 等人（1979）A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code
Sobol'（1967）On the Distribution of Points in a Cube and the Approximate Evaluation of Integrals
Jones 等人（1998）Efficient Global Optimization of Expensive Black-Box Functions
Loeppky 等人（2009）Choosing the Sample Size of a Computer Experiment: A Practical Guide

11.5 简史 #

贝叶斯优化的历史比它的名称所暗示的更长。它的思想来自统计学、运筹学和工程设计，其中有几个思想曾被不止一次地提出。

模型与决策规则（20 世纪 60 年代）。自 20 世纪 40 年代起，统计学家就在研究如何序贯地设计实验，即根据上一次实验的结果选择下一次实验（Garnett，2023，第 12.2 节）。Kushner（1964）把这一思想用于寻找带噪声观测的一维函数的最大值。他用 Wiener 过程为函数建模；这种高斯过程的样本路径形如随机游走的轨迹，处处连续而处处不光滑。早期工作偏爱这类过程，因为其更新代价低，当时的计算机足以承受（Garnett，2023，第 12.3 节）。Kushner 认为最优序贯策略在计算上不可行，因而将其搁置，转而提出更简单的规则，其中包括最大化超过迄今最优值的概率（第 12.2 节）。他的论文还讨论了人类专家如何在搜索过程中调整这条规则的改进阈值（Garnett，2023，第 12.3 节），这是人在回路的早期例子，这一主题将在第四部分中再次出现。

一步前瞻（20 世纪 70 年代）。苏联的一个研究方向发展出了恰好向前看一次评估的采集函数。期望改进（第 12.3 节）通常归功于 Močkus 及其同事（Močkus，1975；Jones 等，1998；Frazier，2018；Brochu 等，2010）。Garnett 的史述把期望改进的显式公式追溯到 Šaltenis 1971 年的工作，并把 Močkus 的一步准则（以一次评估能把任意位置的最优期望值提高多少来衡量其价值）解读为如今所说的知识梯度（第 12.6 节）；在 Wiener 过程下，两个准则是一致的（Garnett，2023，第 12.3 节）。无论如何，到 20 世纪 70 年代初，这两个一步准则都已发表。

克里金法与计算机实验（20 世纪 50 年代至 90 年代）。在另一条独立的线索上，矿山中的矿石品位估计以克里金法（kriging）之名催生了高斯过程回归（Krige，1951；Matheron，1963）。Sacks 等人（1989）把它引入计算机实验的设计与分析，在这类实验中，昂贵的模拟取代了物理实验。Jones 等人（1998）把这种模型与闭式的期望改进结合起来，同时细致地处理了模型验证，并用分支定界法最大化采集函数，称之为 Efficient Global Optimization，简称 EGO。EGO 使这一方法受到广泛关注，首先是在工程设计领域（Frazier，2018）。

来自赌博机的理论保证（2010 年）。式（11.3）中的上置信界由 Kushner 提出，此后又多次被重新发现（Garnett，2023，第 12.5 节）。Srinivas 等人（2010）把它与多臂赌博机文献联系起来，并证明了它在高斯过程模型下的收敛速度，第 13.4 节将解释这一分析。

机器学习（2012 年以后）。Snoek 等人（2012）表明，只要仔细选择先验及其超参数的处理方式，贝叶斯优化为机器学习算法（包括卷积神经网络）调参的效果就能与人类专家相当，甚至更好。这篇论文在机器学习界引发了研究热潮：Garnett 2023 年的教科书所引文献中，一半以上发表于 2012 年之后（Garnett，2023，第 12.4 节），随后又出现了 BoTorch 等软件框架（Balandat 等，2020）。同一时期，另一个研究方向发展出基于信息论的采集函数：最初由 Villemonteix 及其同事提出，由 Hennig 与 Schuler（2012）命名为熵搜索（entropy search）（Garnett，2023，第 12.4 节），此后又经 Hernández-Lobato 等人（2014）和 Wang 与 Jegelka（2017）改进（第 12.7 节）。

偏好（2007 年以后）。Brochu 等人（2007）使用同样的机制设计渲染材质的外观，只是由人在选项之间做出选择，而不是报告数值。这一研究方向就是偏好贝叶斯优化，即第四部分的主题。

表 11.1 贝叶斯优化发展中的里程碑。

年份	工作	贡献
1933	Thompson（1933）	按各治疗方案更优的后验概率分配治疗；Thompson 采样的起源（第 12.5 节）
1964	Kushner（1964）	采用 Wiener 过程模型的一维优化；改进概率
1975	Močkus（1975）	贝叶斯方法中的一步前瞻；期望改进通常归功于此
1989	Sacks 等人（1989）	用于昂贵计算机模拟的高斯过程模型
1998	Jones 等人（1998）	EGO：基于拟合高斯过程的闭式期望改进
2009	Frazier 等人（2009）	相关信念下的知识梯度
2010	Srinivas 等人（2010）	GP-UCB 及其遗憾界
2012	Hennig 与 Schuler（2012）	熵搜索：按关于最大值点的信息选择评估
2012	Snoek 等人（2012）	机器学习算法的超参数优化
2020	Balandat 等人（2020）	BoTorch：结合自动微分的蒙特卡洛采集函数

第 11.5 节引用的文献 18

Garnett（2023）Bayesian Optimization
Kushner（1964）A New Method of Locating the Maximum Point of an Arbitrary Multipeak Curve in the Presence of Noise
Močkus（1975）On Bayesian Methods for Seeking the Extremum
Jones 等人（1998）Efficient Global Optimization of Expensive Black-Box Functions
Frazier（2018）A Tutorial on Bayesian Optimization
Brochu 等人（2010）A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning
Krige（1951）A Statistical Approach to Some Basic Mine Valuation Problems on the Witwatersrand
Matheron（1963）Principles of Geostatistics
Sacks 等人（1989）Design and Analysis of Computer Experiments
Srinivas 等人（2010）Gaussian Process Optimization in the Bandit Setting: No Regret and Experimental Design
Snoek 等人（2012）Practical Bayesian Optimization of Machine Learning Algorithms
Balandat 等人（2020）BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization
Hennig 与 Schuler（2012）Entropy Search for Information-Efficient Global Optimization
Hernández-Lobato 等人（2014）Predictive Entropy Search for Efficient Global Optimization of Black-box Functions
Wang 与 Jegelka（2017）Max-value Entropy Search for Efficient Bayesian Optimization
Brochu 等人（2007）Active Preference Learning with Discrete Choice Data
Thompson（1933）On the Likelihood that One Unknown Probability Exceeds Another in View of the Evidence of Two Samples
Frazier 等人（2009）The Knowledge-Gradient Policy for Correlated Normal Beliefs

11.6 习题 #

习题 11.1

设高斯过程先验具有常数均值 $m$ 和平稳核 $k(\vx, \vx') = \kappa(\vx - \vx')$ ，其中 $\kappa$ 为某个函数。证明在没有任何数据时， $\mu_0(\vx)$ 和 $\sigma_0(\vx)$ 不依赖于 $\vx$ 。若采集函数只通过 $\mu_0(\vx)$ 和 $\sigma_0(\vx)$ 依赖于 $\vx$ ，这一结论意味着什么？

解答

没有数据时，后验即先验，因此 $\mu_0(\vx) = m$ ， $\sigma_0^2(\vx) = k(\vx, \vx) = \kappa(\mathbf{0})$ ，在每个输入处都相同。于是形如 $a_0(\vx) = h(\mu_0(\vx), \sigma_0(\vx))$ 的采集函数为常数，每个输入都是它的最大值点。第一次查询因而是任意的，这正是用设计规则选取最初几个点的理由之一。

习题 11.2

设观测精确，贪心规则 $a_n(\vx) = \mu_n(\vx)$ 选中了已评估过的输入 $\vx_i$ 。证明再次评估它不会改变后验，因此循环将一直选择同一个输入。

解答

观测精确时，已评估输入处的后验方差为零，后验均值等于观测值： $\sigma_n(\vx_i) = 0$ ， $\mu_n(\vx_i) = y_i$ （第 8.2 节）。在 $\vx_i$ 处重新评估仍返回 $y_i$ ，而模型早已确定地预测了这个值。以概率为一的事件为条件不会改变分布，因此在每个 $\vx$ 处都有 $\mu_{n+1}(\vx) = \mu_n(\vx)$ 和 $\sigma_{n+1}(\vx) = \sigma_n(\vx)$ 。于是贪心规则再次选择 $\vx_i$ ，如此反复，直到预算用完。（在浮点运算中，重复的输入会使核矩阵奇异；第 8.4 节中的小抖动项能使分解继续进行，对后验的影响可以忽略不计。）

习题 11.3

随机搜索从定义域中均匀抽取 $n$ 个输入。设 $p$ 为定义域中 $f$ 与其最大值之差在某一容差以内的区域所占的体积比例。证明 $n$ 个输入中至少有一个落入该区域的概率为 $1 - (1 - p)^n$ ，并求出 $p = 0.05$ 时使这一概率至少为 $0.95$ 的最小 $n$ 。为什么这个数不依赖于维度 $d$ ？为什么这一点并不像听上去那样令人放心？

解答

每个输入独立地以概率 $1 - p$ 错过该区域，因此 $n$ 个输入全部错过的概率为 $(1 - p)^n$ ，至少一个命中的概率为 $1 - (1 - p)^n$ 。令 $1 - 0.95^n \ge 0.95$ ，得 $n \ge \log 0.05 / \log 0.95 \approx 58.4$ ，故 $n = 59$ 。这一计算只用到体积比例 $p$ ，与 $d$ 无关。但在 $d$ 维中，若一个区域在每个输入上都占其取值范围的比例 $q$ ，则其体积比例为 $p = q^d$ 。取 $q = 0.5$ 、 $d = 10$ ，则 $p \approx 0.001$ ，所需的 $n$ 增长到约 3,000。只有少数几个输入重要时，随机搜索是很强的基线（Bergstra 与 Bengio，2012），因为此时 $p$ 只由这几个输入决定。

第 11.6 节引用的文献 1

Bergstra 与 Bengio（2012）Random Search for Hyper-Parameter Optimization

参考文献

Balandat, M., Karrer, B., Jiang, D. R., Daulton, S., Letham, B., Wilson, A. G., and Bakshy, E. (2020). BoTorch: A Framework for Efficient Monte-Carlo Bayesian Optimization. Advances in Neural Information Processing Systems 33 (NeurIPS 2020). 引用于 §11.5
Bergstra, J., and Bengio, Y. (2012). Random Search for Hyper-Parameter Optimization. Journal of Machine Learning Research. 引用于 §11.4 §11.6
Brochu, E., de Freitas, N., and Ghosh, A. (2007). Active Preference Learning with Discrete Choice Data. Advances in Neural Information Processing Systems. 引用于 §11.5
Brochu, E., Cora, V. M., and de Freitas, N. (2010). A Tutorial on Bayesian Optimization of Expensive Cost Functions, with Application to Active User Modeling and Hierarchical Reinforcement Learning. arXiv preprint. 预印本引用于 §11.5
Frazier, P. I. (2018). A Tutorial on Bayesian Optimization. arXiv. 预印本引用于 §11.1 §11.2 §11.5
Frazier, P., Powell, W., and Dayanik, S. (2009). The Knowledge-Gradient Policy for Correlated Normal Beliefs. INFORMS Journal on Computing. 引用于 §11.5
Garnett, R. (2023). Bayesian Optimization. Cambridge University Press. 引用于 §11.2 §11.4 §11.5
Hennig, P., and Schuler, C. J. (2012). Entropy Search for Information-Efficient Global Optimization. Journal of Machine Learning Research. 引用于 §11.5
Hernández-Lobato, J. M., Hoffman, M. W., and Ghahramani, Z. (2014). Predictive Entropy Search for Efficient Global Optimization of Black-box Functions. Advances in Neural Information Processing Systems 27 (NeurIPS 2014). 引用于 §11.5
Jones, D. R., Schonlau, M., and Welch, W. J. (1998). Efficient Global Optimization of Expensive Black-Box Functions. Journal of Global Optimization. 引用于 §11.4 §11.5
Krige, D. G. (1951). A Statistical Approach to Some Basic Mine Valuation Problems on the Witwatersrand. Journal of the Southern African Institute of Mining and Metallurgy. 引用于 §11.5
Kushner, H. J. (1964). A New Method of Locating the Maximum Point of an Arbitrary Multipeak Curve in the Presence of Noise. Journal of Basic Engineering. 引用于 §11.5
Loeppky, J. L., Sacks, J., and Welch, W. J. (2009). Choosing the Sample Size of a Computer Experiment: A Practical Guide. Technometrics. 引用于 §11.4
Matheron, G. (1963). Principles of Geostatistics. Economic Geology. 引用于 §11.5
McKay, M. D., Beckman, R. J., and Conover, W. J. (1979). A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code. Technometrics. 引用于 §11.4
Močkus, J. (1975). On Bayesian Methods for Seeking the Extremum. Optimization Techniques IFIP Technical Conference. 引用于 §11.5
Sacks, J., Welch, W. J., Mitchell, T. J., and Wynn, H. P. (1989). Design and Analysis of Computer Experiments. Statistical Science. 引用于 §11.5
Shahriari, B., Swersky, K., Wang, Z., Adams, R. P., and de Freitas, N. (2016). Taking the Human Out of the Loop: A Review of Bayesian Optimization. Proceedings of the IEEE.
Shields, B. J., Stevens, J., Li, J., Parasram, M., Damani, F., Alvarado, J. I. M., … Doyle, A. G. (2021). Bayesian reaction optimization as a tool for chemical synthesis. Nature. 引用于 §11.1
Snoek, J., Larochelle, H., and Adams, R. P. (2012). Practical Bayesian Optimization of Machine Learning Algorithms. Advances in Neural Information Processing Systems 25 (NeurIPS 2012). 引用于 §11.1 §11.3 §11.4 §11.5
Sobol', I. M. (1967). On the Distribution of Points in a Cube and the Approximate Evaluation of Integrals. USSR Computational Mathematics and Mathematical Physics. 引用于 §11.4
Srinivas, N., Krause, A., Kakade, S. M., and Seeger, M. (2010). Gaussian Process Optimization in the Bandit Setting: No Regret and Experimental Design. ICML 2010. 引用于 §11.5
Thompson, W. R. (1933). On the Likelihood that One Unknown Probability Exceeds Another in View of the Evidence of Two Samples. Biometrika. 引用于 §11.5
Wang, Z., and Jegelka, S. (2017). Max-value Entropy Search for Efficient Bayesian Optimization. Proceedings of the 34th International Conference on Machine Learning (ICML 2017). 引用于 §11.5