核函数与超参数

第 7 章表明，核函数是对未知函数的陈述；第 8 章说明了数据到来后这一陈述如何变化。两章都手动固定了核函数及其超参数。每张图都有长度尺度滑块，移动滑块，各处的预测和不确定性都随之改变。贝叶斯优化器却没有人替它移动滑块。它要了解目标函数，依靠的是寥寥几次评估；核函数的设定也必须从这几次评估中选出。

本章讨论的就是这一选择，核心是边际似然。边际似然是一个数，表示在给定核函数与超参数下观测数据出现的概率；“选哪个长度尺度？”这一问题，通常的答案就是使边际似然最大的取值。在此之前，先介绍可供选择的核函数；在此之后，讨论这一答案会以哪些方式误导人，以及为什么输入越多，误导越严重。

9.1 核函数族 #

平稳核只通过输入之间的距离 $r = \lVert \vx - \vx' \rVert$ 依赖于输入（第 7.5.1 节）。贝叶斯优化中使用的平稳核彼此只有一处差异：抽出的函数光滑程度不同。

9.1.1 Matérn 阶梯 #

Matérn 族有光滑度参数 $\nu > 0$ （第 7.5.3 节）。 $\nu$ 为半整数时，核函数是指数函数与多项式之积；实践中用这样三个取值就已足够。表 9.1 列出了这三个核函数，以及作为 $\nu \to \infty$ 极限的径向基函数核。

表 9.1 常用的平稳核，取单位幅度，写成距离 r 的函数；乘以幅度的平方即得一般形式。当且仅当 ν 大于 q 时，样本 q 次可微。

核函数	$\nu$	$k(r)$	样本
Matérn 1/2（指数核）	$1/2$	$\exp\!\left(-\dfrac{r}{\ell}\right)$	连续，处处不可微
Matérn 3/2	$3/2$	$\left(1 + \dfrac{\sqrt{3}\,r}{\ell}\right)\exp\!\left(-\dfrac{\sqrt{3}\,r}{\ell}\right)$	一次可微
Matérn 5/2	$5/2$	$\left(1 + \dfrac{\sqrt{5}\,r}{\ell} + \dfrac{5r^2}{3\ell^2}\right)\exp\!\left(-\dfrac{\sqrt{5}\,r}{\ell}\right)$	二次可微
径向基函数核（平方指数核）	$\infty$	$\exp\!\left(-\dfrac{r^2}{2\ell^2}\right)$	无穷次可微

这些公式和可微性规则都是标准结果（Rasmussen 与 Williams，2006，第 4.2.1 节）；可微性指第 7.5.3 节中的均方意义。参数 $\nu$ 将在第 13.4.3 节中再次出现，见于优化器遗憾（与可达到的最优值之间的总差距）的界：贝叶斯优化的理论结论是针对各个核函数分别给出的。

第 7 章中的两个核函数补全了这一组。周期核（式（7.7））抽出严格重复的函数。线性核（linear kernel）是贝叶斯线性回归的核函数（例 7.1），

k(\vx, \vx') = \sigma_0^2 + \sigma_1^2\, (\vx - \mathbf{c})^\T (\vx' - \mathbf{c}),

其样本是截距和斜率随机的直线（多个输入时为平面），绕点 $\mathbf{c}$ 转动。线性核不是平稳的：方差随着与 $\mathbf{c}$ 的距离增大而增长。

9.1.2 和与积 #

核函数可以组合，两条基本规则的证明都很简短。

两个核函数之和是核函数。若 $f_1 \sim \GP(0, k_1)$ 与 $f_2 \sim \GP(0, k_2)$ 相互独立，两者之和是核函数为 $k_1 + k_2$ 的高斯过程（习题 7.3）。核函数之和适合刻画由相互独立的部分组成的函数，例如缓慢的趋势加快速的波动，或信号加相关噪声。

两个核函数之积是核函数。上述两个过程之积 $f_1 f_2$ 的协方差为 $k_1 k_2$ （习题 9.1）。乘积过程不服从高斯分布，但其协方差函数是半正定的，而核函数只需满足这一条件（Rasmussen 与 Williams，2006，第 4.2.4 节）。只有两个因子都大时乘积才大，因此核函数之积刻画的是必须同时在两种意义上成立的结构：周期核乘以径向基函数核表示“会重复，并且相邻的重复比相距遥远的重复更相似”。

乘积也是让核函数覆盖多个输入的方法。若 $k_1$ 作用于第一个输入， $k_2$ 作用于第二个输入，则 $k_1(x_1, x_1')\,k_2(x_2, x_2')$ 就是输入对上的核函数。多输入的径向基函数核恰好就是一维径向基函数核的这种乘积，因为和的指数等于指数的积。跨输入求和 $k_1(x_1, x_1') + k_2(x_2, x_2')$ 则表示函数是若干函数之和，每个函数只依赖一个输入，彼此没有交互，这是强得多的假设。

下图将这一族核函数并排展示。每个面板抽取两个函数，所有面板使用同一组随机数，因此面板之间的差异只来自核函数。

图 9.1 六个核函数，同一组随机数。每个面板显示从高斯过程先验中抽取的两个函数，下方是核函数

k(0.7, x)

。“基本构件”：从 1/2 到径向基函数核的 Matérn 阶梯、周期为 0.25 的周期核，以及线性核。“和与积”：由上述核函数构造的核函数。按“重新抽取”可得到新的随机数；长度尺度滑块作用于径向基函数核和 Matérn 核。组合中的权重和周期仅作示意。

按顺序看前四个面板。由于随机数相同，各 Matérn 样本的大尺度形状一致，粗糙程度则逐级降低。乍看之下，Matérn 5/2 的样本与径向基函数核的样本很接近，差异在于各面板下方核函数曲线的尾部，以及样本的细节。

切换到“和与积”。长短两个径向基函数核相加，得到缓慢的漂移，其上叠加快速的波动。径向基函数核加周期核，得到叠加在趋势上的重复图案；两者之积得到重复出现、形状却缓慢变化的图案。线性核乘以自身抽出抛物线，乘以周期核则抽出幅度随远离中心而增大的振荡。

下面用一个真实的例子说明组合的好处。夏威夷莫纳罗亚（Mauna Loa）测得的二氧化碳月浓度，1958 至 2003 年共 545 个观测，是一个标准的演示例子（Rasmussen 与 Williams，2006，第 5.4.3 节）。其中的核函数由四部分相加而成，每一部分对应记录中的一种特征：径向基函数核刻画长期上升；周期核乘以径向基函数核，刻画可能缓慢变化的季节周期；第三项刻画几年尺度上的不规则变化；最后是噪声项。该核函数有 11 个超参数，全部用第 9.3 节的方法拟合。拟合值读来如同一份数据报告：趋势的长度尺度为 67 年，季节模式在 90 年的尺度上衰减，因而接近严格周期。寻找这类结构的过程本身也可以自动化：以贪心方式对基本核函数做加法和乘法组合，用边际似然为每个候选打分（Duvenaud 等，2013）。

贝叶斯优化很少做到这一步。只有几十次评估时，数据太少，不足以在不同结构之间作出选择，因此通常的模型是单个 Matérn 5/2 核或径向基函数核，每个输入有单独的长度尺度。最后这一点值得用一节专门讨论。

第 9.1 节引用的文献 2

Rasmussen 与 Williams（2006）Gaussian Processes for Machine Learning
Duvenaud 等人（2013）Structure Discovery in Nonparametric Regression through Compositional Kernel Search

9.2 每个输入一个长度尺度 #

目标函数有多个输入时，很少同等地依赖于每一个输入。神经网络的验证误差可能随训练过程中的某个设定（学习率）大幅波动，对另一个设定却几乎没有反应。单一的长度尺度无法表达这一点，因为它假定函数在各个方向上以相同的速率变化。

补救办法是为每个输入单独配一把尺子度量距离，把平方距离 $\lVert \vx - \vx' \rVert^2/\ell^2$ 换成逐输入项之和：

k(\vx, \vx') = \sigma_f^2 \exp\!\left(-\frac12 \sum_{j=1}^{d} \frac{(x_j - x_j')^2}{\ell_j^2}\right),

其中 $d$ 是输入个数；Matérn 核同样如此处理，即在表 9.1 中取 $r^2 = \sum_j (x_j - x_j')^2/\ell_j^2$ 、 $\ell = 1$ 。每个 $\ell_j$ 表示沿输入 $j$ 移动多远，函数才会有明显变化。 $\ell_j$ 较短，说明函数对输入 $j$ 敏感。随着 $\ell_j$ 增大，输入 $j$ 对应的项从和式中消失，核函数不再区分该输入上的差异，先验抽出的每个函数沿这一输入都成为常数。这个输入虽未移除，实际上已经关闭。

长度尺度由数据拟合，因此模型能够发现哪些输入重要。这称为自动相关性确定（automatic relevance determination，ARD），术语出自 Neal（1996）；长度尺度的倒数 $1/\ell_j$ 解读为输入 $j$ 的相关性（Rasmussen 与 Williams，2006，第 5.1 节）。

图 9.2 从高斯过程先验中抽取的一个双输入函数，先验使用式（9.3）的径向基函数核。热图显示函数在单位正方形上的值；两个面板分别沿两个输入方向过中心作切片。调长

\ell_2

，热图变成竖直条纹，沿

x_2

的切片同时变平：函数不再依赖于这一输入。移动滑块时随机数保持不变。这里只使用径向基函数核，其乘积形式使二维采样的代价很低。

把两个长度尺度都设为 0.2。热图呈现由直径约 0.2 的圆形丘陵和洼地构成的地形，两个切片起伏程度相同。

把 $\ell_2$ 调长到 5。丘陵拉伸为沿 $x_2$ 方向延伸的条纹，沿 $x_2$ 的切片几乎是一条水平线。函数仍是随机的，但只依赖于 $x_1$ 。

按“交换长度尺度”。条纹旋转四分之一圈。哪个输入重要，取决于核函数，而不取决于随机数。

真实的调参问题正依赖于此。Snoek 等人（2012）提出用 ARD Matérn 5/2 核调节机器学习模型。他们用论文中的高斯过程模型，在 CIFAR-10 图像基准上调节卷积网络的 9 个超参数，测试误差达到 14.98%，比专家找到的设定低 3 个多百分点。在第 22.4 节的 7 个超参数的问题中，误差的大部分变化由 7 个中的 2 个解释（第 22.4.1 节），这正是 ARD 适用的情形。第 9.5 节中的图 9.4 展示了拟合得到的长度尺度如何把 6 个重要的输入与 14 个不重要的输入区分开。

ARD 也有代价：每个输入多一个超参数，每个超参数都要从同样寥寥几次评估中学习。观测为实数时，这一代价通常可以承受。若每个观测只是两个选项之间的一次比较（第四部分），就未必如此：第 25 章直接固定了六个长度尺度，没有从数据中学习；第 30.3.3 节则指出，尚无研究检验过在这样的预算下能否学到它们。反过来，若模型沿某个属性需要比其他属性短得多的长度尺度，这恰恰是光滑的效用函数在近似一个优先关注该属性的人（第 37.5.2 节）。

第 9.2 节引用的文献 3

Neal（1996）Bayesian Learning for Neural Networks
Rasmussen 与 Williams（2006）Gaussian Processes for Machine Learning
Snoek 等人（2012）Practical Bayesian Optimization of Machine Learning Algorithms

9.3 边际似然 #

记 $\bm{\theta}$ 为全部超参数：各长度尺度、幅度 $\sigma_f$ 和噪声标准差 $\sigma_n$ 。问题是如何根据数据选择 $\bm{\theta}$ 。

最直观的准则行不通。若按后验均值与观测的贴合程度给设定打分，胜出的总是最短的长度尺度和最小的噪声：这样的模型曲折地精确穿过每个点，但在点与点之间什么也预测不了（第 8.2 节）。能拟合任何数据的模型，什么也没有学到。

第 5.6 节在直线的情形中遇到过这一问题，并以模型证据作答：按模型在看到观测数据之前赋予这些数据的概率给模型打分。同样的思路在此也适用，而且对高斯过程，这一概率有闭式解。

推导对数边际似然

第 8.3 节的模型为 $y_i = f(\vx_i) + \varepsilon_i$ 。根据高斯过程的定义（定义 7.1）， $n$ 个已观测输入处的函数值服从高斯分布， $\vf \sim \N(\mathbf{0}, \mK)$ ，其中 $[\mK]_{ij} = k(\vx_i, \vx_j)$ 。
噪声为 $\bm{\varepsilon} \sim \N(\mathbf{0}, \sigma_n^2\mI)$ ，与 $\vf$ 独立。
独立高斯向量之和仍是高斯向量，协方差相加（第 4.6.1 节，经由式（4.9）得到）。因此 $\vy = \vf + \bm{\varepsilon} \sim \N(\mathbf{0}, \mK_y)$ ，其中 $\mK_y = \mK + \sigma_n^2\mI$ 。
数据的概率密度即该高斯密度（式（4.5））在观测值 $\vy$ 处的取值： $p(\vy \given X, \bm{\theta}) = (2\pi)^{-n/2}\,\lvert\mK_y\rvert^{-1/2} \exp\!\left(-\tfrac12 \vy^\T\mK_y^{-1}\vy\right)$ 。
取对数，结果记为 $\mathcal{L}(\bm{\theta}) = \log p(\vy \given X, \bm{\theta})$ 。

\mathcal{L}(\bm{\theta}) = \underbrace{-\tfrac12 \vy^\T \mK_y^{-1} \vy}_{\text{data fit}} \underbrace{-\tfrac12 \log \lvert \mK_y \rvert}_{\text{complexity}} -\tfrac{n}{2}\log 2\pi.

这就是对数边际似然（log marginal likelihood）（Rasmussen 与 Williams，2006，第 5.4.1 节）。称为“边际”，是因为未知的函数值已经通过积分消去：第 3 步用一行完成了积分 $\int p(\vy \given \vf)\, p(\vf \given X, \bm{\theta})\, \dd\vf$ 。它就是第 5.6 节中的公式（式（5.11）），只是以核矩阵代替 $\boldsymbol{\Phi}\mSigma_p\boldsymbol{\Phi}^\T$ ；它只通过 $\mK_y$ 依赖于 $\bm{\theta}$ 。

9.3.1 相互拉扯的两项 #

单个观测即可说明其中的机制。

例 9.1 单个观测

只有一个观测 $y$ 时， $\mK_y$ 是一个数 $s^2 = \sigma_f^2 + \sigma_n^2$ ，即模型预期的总方差，且

\log p(y) = -\frac{y^2}{2s^2} - \frac12\log s^2 - \frac12\log 2\pi.

第一项奖励较大的 $s^2$ ：模型若预期取值很大，远离零的观测就不足为奇。第二项则惩罚它：概率分散在很宽范围上的模型，分给每个特定值的概率就更少。令关于 $s^2$ 的导数为零，得 $s^2 = y^2$ 。最好的模型所预期的取值大小，恰好等于它观测到的大小，不多也不少。至于 $s^2$ 如何分为信号和噪声，数据没有提供任何信息，任何单个观测也都无法提供。

有 $n$ 个观测时，两项的作用不变。数据拟合项（data fit） $-\tfrac12\vy^\T\mK_y^{-1}\vy$ 是唯一含有观测值的项。它是 $\vy$ 到零的平方距离，以模型预期的单位度量（即第 4.2.1 节中的 Mahalanobis 距离）；数据若按核函数预期的方式变化，这一项的负值最小。

复杂度项（complexity） $-\tfrac12\log\lvert\mK_y\rvert$ 完全不含观测值。行列式是模型预期数据所落区域的体积（第 3.6 节），因此这一项会为模型本可以解释的每个数据集计费，无论该数据集是否真的出现。短长度尺度使 $\vy$ 的各元素几乎相互独立， $\mK_y$ 近乎对角，体积在方差允许的范围内达到最大。长长度尺度把各观测联系在一起，使区域变扁、体积缩小。因此在噪声水平固定时，加长 $\ell$ 会减轻复杂度的代价，数据拟合项则变差，因为更僵硬的函数能匹配的数据集更少（Rasmussen 与 Williams，2006，第 5.4.1 节）。两项之和取最大值的长度尺度，就是在数据允许的范围内使模型最简单的长度尺度。这就是第 5.6.1 节中的自动 Occam 剃刀，它无需留出数据。

9.3.2 曲面 #

下图对七个观测，在长度尺度与噪声水平构成的网格上计算式（9.4），幅度固定为 $\sigma_f = 1$ 。

图 9.3 单位幅度径向基函数核在右侧七个观测上的对数边际似然。左：各长度尺度与噪声水平下的取值；每级色阶对应一段对数概率区间（与最大值相差 0.25 以内，其次为 0.5、1、2、4、8、16），十字标出局部极大值。右：所选参数对给出的回归拟合。下：在所选噪声水平下沿长度尺度方向切过曲面的剖面，分解为数据拟合项与复杂度项（式（9.4））。可在曲面上拖动圆环，或使用滑块。点击拟合面板可添加观测，点击某个点可将其删除。数据仅作示意。

从图的初始状态开始。圆环位于一个局部极大值：长度尺度 0.44，噪声 0.74。右侧的拟合是一条平缓的下坡，模型把其余一切都归于噪声。

按“前往最佳参数对”。圆环跳到长度尺度 0.075、噪声 0.23 处。此时拟合曲折地穿过各点，阴影带在点与点之间明显变宽。对数边际似然从 $-9.65$ 升到 $-9.18$ ，因此这种解释的概率是另一种的 $e^{0.46} \approx 1.6$ 倍。仅凭七个点，无法在“精确测量的曲折函数”与“测量粗糙的光滑函数”之间作出裁决。

从最佳参数对出发，把圆环竖直向下拖。数值几乎不变：噪声为 0.01 时为 $-9.23$ 。一旦函数穿过各点，模型就无法区分小噪声与无噪声，曲面成为一道平坦的山脊。

移动长度尺度时，观察底部的剖面。在初始噪声水平下，随着长度尺度从 0.05 增大到 2，数据拟合曲线从 $-2.1$ 降到 $-5.5$ ，复杂度曲线从 $-1.4$ 升到 $0.7$ 。两者之和在 0.1 至 0.44 之间几乎持平。把噪声降到 0.23，数据拟合曲线变为陡崖：长度尺度为 0.5 时降至 $-24$ ，因为噪声很小的僵硬曲线无法靠近这些点。

先按“数据：三个点”，再按“数据：十个点”。只有三个点时，最大值位于曲面的角落，短长度尺度、小噪声的一大片区域与最大值只差零点零几。三个点几乎不提供关于长度尺度的信息。有十个点时，曲面在长度尺度 0.13、噪声 0.27 处形成一个集中的峰。

标准参考书中针对另外七个观测给出了同样的图景，并直接点明了其中的教训：每个局部极大值都是对数据的一种特定解释，点这么少时，模型无法有把握地排除其中任何一种（Rasmussen 与 Williams，2006，第 5.4.1 节）。

要点边际似然依据模型的预测为其打分

边际似然是模型在看到观测数据之前赋予这些数据的概率。仅仅很好地拟合数据还不够；模型还不能预测出许多实际并未出现的数据集。

第 9.3 节引用的文献 1

Rasmussen 与 Williams（2006）Gaussian Processes for Machine Learning

9.4 拟合超参数 #

标准做法是选择使式（9.4）最大的超参数，称为第二类最大似然（type II maximum likelihood）或经验贝叶斯：对先验的设定而非函数本身应用最大似然（第 5.6 节）。

9.4.1 攀登曲面 #

图 9.3 那样的网格适用于两个超参数，却不适用于十个，因此需要沿梯度寻找最大值。

推导对数边际似然的梯度

设 $\theta_j$ 为某个超参数，将逐元素求导得到的矩阵 $\partial\mK_y/\partial\theta_j$ 简记为 $\partial\mK_y$ 。

需要用到两条矩阵求导规则（Petersen 与 Pedersen，2012；Rasmussen 与 Williams，2006，附录 A.3）： $\partial(\mK_y^{-1}) = -\mK_y^{-1}(\partial\mK_y)\mK_y^{-1}$ 和 $\partial \log\lvert\mK_y\rvert = \tr\!\left(\mK_y^{-1}\,\partial\mK_y\right)$ ，其中 $\tr$ 表示迹，即矩阵对角元之和。
将第一条规则用于式（9.4）的数据拟合项： $\partial\!\left(-\tfrac12\vy^\T\mK_y^{-1}\vy\right) = \tfrac12\,\vy^\T\mK_y^{-1}(\partial\mK_y)\mK_y^{-1}\vy = \tfrac12\,\bm{\alpha}^\T(\partial\mK_y)\,\bm{\alpha}$ ，其中 $\bm{\alpha} = \mK_y^{-1}\vy$ 即第 8.2.1 节中的权重。
将第二条规则用于复杂度项： $\partial\!\left(-\tfrac12\log\lvert\mK_y\rvert\right) = -\tfrac12\tr\!\left(\mK_y^{-1}\,\partial\mK_y\right)$ 。
标量等于其自身的迹，且对任意两个使乘积有定义的矩阵有 $\tr(\mA\mathbf{B}) = \tr(\mathbf{B}\mA)$ 。把最后一个因子移到最前面，可知 $\bm{\alpha}^\T(\partial\mK_y)\bm{\alpha}$ 等于 $\tr\!\left(\bm{\alpha}\bm{\alpha}^\T\,\partial\mK_y\right)$ 。
将第 2 步与第 3 步的结果相加。

\frac{\partial\mathcal{L}}{\partial\theta_j} = \frac12 \tr\!\left( \left(\bm{\alpha}\bm{\alpha}^\T - \mK_y^{-1}\right) \frac{\partial\mK_y}{\partial\theta_j} \right).

算法 8.1 中代价为 $O(n^3)$ 的 Cholesky 分解给出 $\bm{\alpha}$ 、对数行列式（第 3.6.1 节）和 $\mK_y^{-1}$ ；此后每个超参数的导数代价为 $O(n^2)$ （Rasmussen 与 Williams，2006，第 5.4.1 节）。实践中，自动微分库可以由计算式（9.4）的代码直接生成梯度。

算法 9.1 用第二类最大似然拟合超参数

输入：缩放到单位立方体的输入 $X$ 、标准化后的观测 $\vy$ 、一个核函数族。

对每个取正值的超参数取对数作为参数，使优化器在无约束的尺度上工作。
选择起点：单位幅度、较小的噪声水平，以及与输入个数相适应的长度尺度（第 9.5 节）。
用 L-BFGS 等拟 Newton 法最大化式（9.4）；这类优化器沿梯度（式（9.5））前进，并根据梯度在各步之间的变化估计曲率。
从其他几个起点重复第 2、3 步，保留边际似然最高的结果。

代码实现 NumPy 与 SciPy

import numpy as np
from scipy.optimize import minimize

# x：已观测的输入，y：标准化后的观测（均为一维数组）
def neg_log_marginal(log_theta, x, y):
    ell, sf, sn = np.exp(log_theta)
    K = sf**2 * np.exp(-0.5 * (x[:, None] - x[None, :]) ** 2 / ell**2)
    L = np.linalg.cholesky(K + (sn**2 + 1e-8) * np.eye(len(x)))
    alpha = np.linalg.solve(L.T, np.linalg.solve(L, y))
    return 0.5 * y @ alpha + np.log(np.diag(L)).sum() + 0.5 * len(x) * np.log(2 * np.pi)

starts = [np.log([0.1, 1.0, 0.1]), np.log([0.5, 1.0, 0.5]), np.log([0.03, 1.0, 0.3])]
fits = [minimize(neg_log_marginal, s, args=(x, y), method="L-BFGS-B") for s in starts]
ell, sf, sn = np.exp(min(fits, key=lambda f: f.fun).x)

Cholesky 因子对角元的对数之和等于对数行列式的一半。不提供梯度函数时，SciPy 采用数值微分；超参数只有三个时，这已经足够。

9.4.2 拟合如何失效 #

图 9.3 已经展示了这一过程误导人的三种方式。

多个极大值。第 4 步要求多次重启，是因为曲面可能有不止一个峰，而梯度方法找到的是离起点最近的峰。数据很少时，各峰可能几乎一样高，优化器报告哪一个，就取决于碰巧从哪里出发。

平坦方向。通向零噪声的山脊，以及只有三个观测时的平台，都是数据无法确定某个超参数的区域，优化器在那里仍会返回一个数值。正因如此，贝叶斯优化循环开始时要先做几次不借助模型选点的评估（第 11.4 节）；否则最初几次拟合都落在平台上。

超参数过多。对 $\bm{\theta}$ 求最大值本身就是一种拟合，同样可能过拟合。每个输入各有一个长度尺度而评估又很少时，关闭大多数输入的设定可能碰巧就能解释数据，模型于是忽略了重要的输入。第 9.5 节展示了 50 维中出现的这种情况。

第四个后果为优化所特有。超参数随新评估的到来而重新拟合，采集函数依据的模型因此不断变化，针对固定核函数的收敛保证也就不再直接适用（第 13.5.2 节）。

9.4.3 先验，以及以平均代替选择 #

有两种改进可以应对这些失效。第一种是在超参数上设定先验 $p(\bm{\theta})$ ，并最大化后验，

\log p(\bm{\theta} \given \vy, X) = \mathcal{L}(\bm{\theta}) + \log p(\bm{\theta}) + \text{const},

这就是最大后验估计（maximum a posteriori estimation，MAP）。先验为平台和山脊加上平缓的坡度，使优化器在数据无从判断之处也有方向可循。BoTorch 的默认模型在长度尺度上就带有这样的先验（Meta Platforms, Inc.，2026k）。

长度尺度的先验需要谨慎设定，因为一个先验是否模糊，取决于它在什么尺度上表述。在 0.01 至 10 之间关于 $\ell$ 平坦的密度，把 90% 的概率放在 $\ell = 1$ 以上；在同一范围内关于 $\log \ell$ 平坦的密度，则在每个十倍区间上各放三分之一的概率。两者都不是中性的（第 2.3.3 节）。因此，长度尺度先验通常在对数尺度上给出，采用对数正态分布，或指定众数的 Gamma 分布。

第二种改进是根本不作选择。全贝叶斯处理在 $\bm{\theta}$ 的后验上对预测取平均：通常用 Markov 链方法抽取 $\bm{\theta}$ 的样本（这种方法能从只确定到一个常数因子的分布中采样，第 17.5 节），再对这些样本上的采集函数取平均。Snoek 等人（2012）采用了这种做法，发现在他们的调参问题上，它优于单一的拟合值。一个合理的解释是，评估很少时，目标函数的不确定性有很大一部分来自超参数的不确定性（推断）。代价是循环的每一步中，每个样本都需要一次 Cholesky 分解。

第 9.4 节引用的文献 4

Petersen 与 Pedersen（2012）The Matrix Cookbook
Rasmussen 与 Williams（2006）Gaussian Processes for Machine Learning
Meta Platforms, Inc.（2026k）botorch/models/utils/gpytorch_modules.py
Snoek 等人（2012）Practical Bayesian Optimization of Machine Learning Algorithms

9.5 长度尺度先验与维度 #

多年来，一直有观点认为贝叶斯优化在输入超过 10 个或 20 个时就不再有效。2024 至 2026 年间，研究者把这一失效追溯到一个平凡的原因：长度尺度的默认先验（Hvarfner 等，2024；Xu 等，2025b；Papenmeier 等，2025b）。

第 3.1.2 节和第 7.5.2 节已经做过这笔算术。单位立方体 $[0, 1]^d$ 中两个随机点的距离约为 $\sqrt{d/6}$ ，即距离随 $\sqrt{d}$ 增长；而典型长度尺度为（例如）0.5 的先验却原地不动。于是在高维中，每一对评估都相距许多个长度尺度，每个核函数值都接近零，模型看到的只是一堆互不相关的点。第 30.1.2 节推导并画出了这一现象。

数据还来不及纠正，损害就已经造成。式（9.5）中要乘以 $\partial\mK_y/\partial\ell$ ，而对径向基函数核，每个元素的导数为 $k \cdot r^2/\ell^3$ 。核函数值接近零时，梯度也接近零，从短长度尺度出发的优化器便无坡可爬（习题 9.4）。Xu 等人（2025b）与 Papenmeier 等人（2025b）认为，初始长度尺度导致的这种梯度消失，是高维中失效的一个主要原因。

解决办法是让长度尺度的先验随维度增长。Hvarfner 等人（2024）让它随 $\sqrt{d}$ 缩放，发现这样一来，在所尝试的 5 个真实任务中，标准贝叶斯优化在 3 个上表现最好，超过了专为高维设计的方法。自 0.12.0 版起，BoTorch 的大多数模型都采用这一先验（Meta Platforms, Inc.，2026e）：位置参数为 $\sqrt{2} + \tfrac12\log d$ 、尺度参数为 $\sqrt{3}$ 的对数正态分布，拟合从先验的众数开始（Meta Platforms, Inc.，2026k）。该众数约为 $0.2\sqrt{d}$ ：10 个输入时为 0.65，20 个输入时为 0.92。

9.5.1 ARD 在 20 维中的表现 #

下图在一个结构已知的测试问题上重放贝叶斯优化。目标函数是 6 维 Hartmann 函数，即有若干局部极大值的标准测试函数，它隐藏在一些不起作用的输入之中：在 20 维中，6 个输入是重要的，14 个是诱饵。模型是采用 ARD 核（式（9.3））的高斯过程，每 5 次评估用最大后验估计重新拟合一次。比较的两种先验分别是：固定的 Gamma(2.4, 2.7) 先验，在每个维度上的众数都为 0.52，这是 BoTorch 偏好模型的默认设定（Meta Platforms, Inc.，2026h）；以及上面的维度缩放先验。

图 9.4 在嵌入 6 至 50 维的 6 维 Hartmann 函数上做贝叶斯优化，由记录的运行重放，每次运行 80 次评估。上：随机搜索、采用固定长度尺度先验的贝叶斯优化和采用维度缩放长度尺度先验的贝叶斯优化的遗憾（与真实最大值的差距），取 6 个随机种子的中位数，阴影表示范围。中：一次运行评估过的所有点，每个输入对应一条竖轴，6 个重要输入所在区域加了阴影。下：该次运行中模型为每个输入学到的长度尺度，对数刻度；短柱表示模型认为相关的输入。可切换显示的运行、维度，以及采集函数的最大化方式。这些是本书自己的运行结果，并非已发表的结果。

先看底部的柱子。在所显示的运行中，缩放先验下的模型最终为 6 个重要输入学到的长度尺度在 0.25 至 0.56 之间，为 14 个诱饵学到的在 2.6 至 30 之间，其中大多数接近 20。模型事先并不知道哪些输入属于哪一类。不用诱饵也能同样好地解释数据，而对于更简单的解释，式（9.4）的复杂度项更大。这正是自动相关性确定名副其实的表现。

切换到固定先验。6 个相关输入的长度尺度相近，在 0.39 至 0.70 之间，但诱饵的长度尺度停留在 1.5 至 2.1 之间。原因在于先验：Gamma(2.4, 2.7) 在长度尺度 2 以上的概率不到 5%，因此在这一先验下，长到足以关闭一个输入的长度尺度不太可能出现；而在 20 维中，缩放先验在该区域放置了约 90% 的概率（推断）。

现在比较遗憾曲线，并改变采集函数的搜索方式。这里的图景不那么整齐，需要如实解读。在 20 维中，若只在均匀随机候选点上最大化采集函数，两种先验最终不相上下：固定先验的中位数遗憾为 0.29，缩放先验为 0.31，随机搜索则为 1.26。在迄今找到的最佳点附近加入候选点后，两者都有所下降，分别降到 0.16 和 0.08。在 50 维中，同样的改变使固定先验从 0.89 降到 0.18，缩放先验从 1.21 降到 0.18。在 10 维中，固定先验领先，0.10 对 0.18。在这些运行中，采集函数的最大化方式比先验的选择更重要（推断）。

转到 50 维，再看一次柱子。在缩放先验下，6 个相关输入中有 2 个的长度尺度此时接近 15：模型关闭了重要的输入。80 次评估太少，不足以确定 50 个长度尺度，这就是第 9.4.2 节中所说的过拟合。

这些运行只是在一个测试函数上、以很小的预算跑了 6 个随机种子，不能据此下定论，但与研究现状一致。上述几篇论文都认同这一诊断：固定的长度尺度先验会随维度增长而失效。这一补救为何有效，则存在争议：Papenmeier 等人（2025b）认为，在极高维度中取得的好结果更多来自局部搜索行为，而非拟合良好的模型。第 30.1 节和第 30.2 节跟进了这一争论。截至 2026 年 9 月，从比较中学习的软件仍未改变默认设置：第 31.2 节考察的每个偏好软件包仍然使用不考虑维度的先验。

第 9.5 节引用的文献 6

Hvarfner 等人（2024）Vanilla Bayesian Optimization Performs Great in High Dimensions
Xu 等人（2025b）Standard Gaussian Process is All You Need for High-Dimensional Bayesian Optimization
Papenmeier 等人（2025b）Understanding High-Dimensional Bayesian Optimization
Meta Platforms, Inc.（2026e）BoTorch CHANGELOG
Meta Platforms, Inc.（2026k）botorch/models/utils/gpytorch_modules.py
Meta Platforms, Inc.（2026h）BoTorch PairwiseGP source code pairwise_gp.py

9.6 检验模型 #

边际似然只给各个模型排出先后，并不说明其中最好的模型是否真的好。即使核函数无法表达目标函数，边际似然仍会在某处取得最大值。在优化器依据模型的不确定性行动之前，有必要检查这种不确定性是否如实；已经收集的数据就能回答这个问题。

所用的检验是留一预测（leave-one-out prediction）：去掉第 $i$ 个观测，用其余 $n - 1$ 个观测预测它，再与实际观测比较。对每个 $i$ 都重复一遍，似乎需要 $n$ 次单独拟合；但对超参数固定的高斯过程，全部 $n$ 个预测都可以由同一个矩阵得到。

推导由一个逆矩阵得到全部留一预测

在该模型下， $\vy \sim \N(\mathbf{0}, \mK_y)$ 。用其余观测预测 $y_i$ ，就是让这一高斯分布以除一个坐标之外的所有坐标为条件。
第 4.5.2 节表明，条件分布可以从精度矩阵 $\bm{\Lambda} = \mK_y^{-1}$ 中读出：坐标 $i$ 的条件方差为 $1/\Lambda_{ii}$ ，条件均值为 $-\Lambda_{ii}^{-1}\sum_{j \ne i}\Lambda_{ij}\,y_j$ 。
对 $j \ne i$ 的求和等于完整求和减去自身对应的一项： $\sum_{j \ne i}\Lambda_{ij}y_j = [\bm{\Lambda}\vy]_i - \Lambda_{ii}\,y_i$ 。
代入，得均值为 $y_i - [\bm{\Lambda}\vy]_i/\Lambda_{ii}$ 。

\mu_{-i} = y_i - \frac{[\mK_y^{-1}\vy]_i}{[\mK_y^{-1}]_{ii}}, \qquad \sigma_{-i}^2 = \frac{1}{[\mK_y^{-1}]_{ii}}.

这就是 $\vx_i$ 处带噪声观测的留一均值与留一方差（Rasmussen 与 Williams，2006，第 5.4.2 节）。尽管从表面看并非如此， $\mu_{-i}$ 并不依赖于 $y_i$ （习题 9.3）。留一预测有三种用途。

校准。标准化残差 $z_i = (y_i - \mu_{-i})/\sigma_{-i}$ 应当近似于标准正态分布的样本：约 95% 落在 $\pm 1.96$ 以内。大残差很多，说明模型过度自信，通常是长度尺度太长或噪声水平太小。残差全都接近零，则说明模型信心不足，优化器会做多于必要的探索。

选择超参数的另一个依据。留一对数密度之和 $\sum_i \log \N(y_i;\, \mu_{-i}, \sigma_{-i}^2)$ 可以代替边际似然，作为最大化的目标。边际似然是假定模型正确时数据的概率；留一得分则无论模型是否正确，都估计预测性能。有观点认为，这使它在核函数设定错误时更稳健（Rasmussen 与 Williams，2006，第 5.4.2 节）。

找出不相符的点。若某个残差远离其余残差，模型就无法让对应的观测与其邻近观测相协调，例如一次失败的运行、一个输错的数值，或函数性质发生改变的区域。

在图 9.3 中打开“留一检验”，即可看到式（9.6）的结果。在最佳参数对处，每个观测都落在（或几乎落在）各自的区间内。拖到长长度尺度、小噪声处，区间缩成一道道短线，偏离了各自的点，这就是过度自信的表现。拖到很短的长度尺度，每个区间都覆盖整个先验范围：模型无法根据邻近观测对某个点作出任何预测，这是如实的，却毫无用处。

在优化中还有一种检验，也是最简单的一种：过迄今最好的点，每次沿一个输入，把后验均值与观测画在同一张图上。若模型在没有评估过的地方阴影带很窄，或在彼此一致的评估之间阴影带很宽，那么在采纳它的下一个建议之前，应当重新检查其超参数。

本章把核函数当作构造协方差矩阵的规则来选择和拟合；第 10 章则把核函数作为独立的数学对象来考察，第 13 章的理论保证正是在这一观点下表述的。

第 9.6 节引用的文献 1

Rasmussen 与 Williams（2006）Gaussian Processes for Machine Learning

9.7 习题 #

习题 9.1

设 $f_1 \sim \GP(0, k_1)$ 与 $f_2 \sim \GP(0, k_2)$ 相互独立，定义 $g(\vx) = f_1(\vx)\,f_2(\vx)$ 。（a）证明 $\Cov[g(\vx), g(\vx')] = k_1(\vx, \vx')\,k_2(\vx, \vx')$ 。（b）尽管 $g$ 不是高斯过程，为什么这仍然证明了两个核函数之积是半正定的？

解答

（a） $g$ 的均值为零，因为由独立性， $\E[f_1(\vx)f_2(\vx)] = \E[f_1(\vx)]\,\E[f_2(\vx)] = 0$ 。因此协方差为 $\E[g(\vx)g(\vx')] = \E[f_1(\vx)f_1(\vx')\,f_2(\vx)f_2(\vx')]$ ，由独立性，这一期望可拆分为 $\E[f_1(\vx)f_1(\vx')]\;\E[f_2(\vx)f_2(\vx')] = k_1(\vx, \vx')\,k_2(\vx, \vx')$ 。（b）第 7.2.2 节中的论证完全没有用到高斯性：对任何方差有限的随机函数， $\sum_{ij} a_i a_j \Cov[g(\vx_i), g(\vx_j)]$ 都是 $\sum_i a_i g(\vx_i)$ 的方差，因而不可能为负。任何随机函数的协方差函数都是有效的核函数；根据高斯过程的定义（定义 7.1），以它为核函数的高斯过程于是存在。

习题 9.2

在单位幅度核函数下，两个无噪声观测的相关系数为 $\rho = k(x_1, x_2)$ 。（a）把对数边际似然写成 $\rho$ 的函数。（b）设两个观测值相同， $y_1 = y_2 = y$ 。证明当 $\rho \to 1$ 时对数边际似然无界增长。拟合出的模型持有什么信念？对于不含噪声项或先验的拟合，这说明了什么？

解答

（a） $\mK = \begin{bmatrix} 1 & \rho \\ \rho & 1 \end{bmatrix}$ 的行列式为 $1 - \rho^2$ ，逆矩阵为 $\frac{1}{1 - \rho^2}\begin{bmatrix} 1 & -\rho \\ -\rho & 1 \end{bmatrix}$ ，所以由式（9.4）

\mathcal{L} = -\frac{y_1^2 - 2\rho\,y_1 y_2 + y_2^2}{2(1 - \rho^2)} - \frac12\log(1 - \rho^2) - \log 2\pi.

（b）当 $y_1 = y_2 = y$ 时，分子为 $2y^2(1 - \rho)$ ，又由于 $1 - \rho^2 = (1 - \rho)(1 + \rho)$ ，数据拟合项为 $-y^2/(1 + \rho)$ ，始终介于 $-y^2$ 与 $-y^2/2$ 之间。当 $\rho \to 1$ 时，复杂度项 $-\tfrac12\log(1 - \rho^2)$ 趋于 $+\infty$ 。因此边际似然在长度尺度为无穷大时取得最大：模型仅凭两个相等的值，就确信函数是常数。而两个相等的值只是这一结论的微弱证据。噪声项（使 $\mK_y$ 远离奇异）或长度尺度上的先验，可以防止优化器跑向这一退化的答案。

习题 9.3

（a）证明式（9.6）中的留一均值 $\mu_{-i}$ 不依赖于 $y_i$ 。（b）两个观测的 $\mK_y = \begin{bmatrix} s^2 & c \\ c & s^2 \end{bmatrix}$ ，由式（9.6）计算 $\mu_{-1}$ 和 $\sigma_{-1}^2$ ；再取 $\sigma_1 = \sigma_2 = s$ 、 $\rho = c/s^2$ 、均值为零，用条件化公式式（4.14）核对结果。

解答

（a）记 $\bm{\Lambda} = \mK_y^{-1}$ 。则 $[\bm{\Lambda}\vy]_i = \Lambda_{ii}y_i + \sum_{j \ne i}\Lambda_{ij}y_j$ ，所以 $\mu_{-i} = y_i - y_i - \Lambda_{ii}^{-1}\sum_{j \ne i}\Lambda_{ij}y_j$ ，两个 $y_i$ 相互抵消。（b）逆矩阵为 $\bm{\Lambda} = \frac{1}{s^4 - c^2}\begin{bmatrix} s^2 & -c \\ -c & s^2 \end{bmatrix}$ 。所以 $\sigma_{-1}^2 = 1/\Lambda_{11} = (s^4 - c^2)/s^2 = s^2 - c^2/s^2$ ， $\mu_{-1} = -\Lambda_{12}\,y_2/\Lambda_{11} = (c/s^2)\,y_2$ 。交换两个坐标的角色后，由式（4.14）得均值 $\rho\,y_2 = (c/s^2)\,y_2$ 和方差 $s^2(1 - \rho^2) = s^2 - c^2/s^2$ ，结果相同。

习题 9.4

对径向基函数核 $k = \exp(-r^2/2\ell^2)$ ，关于长度尺度的导数为 $\partial k/\partial\ell = k\,r^2/\ell^3$ 。（a）把它写成缩放距离 $\rho = r/\ell$ 的函数，并求出使它最大的 $\rho$ 。（b）分别在 $\ell = 0.5$ 和 $\ell = 0.2\sqrt{50}$ 时，计算它在 $[0, 1]^{50}$ 中两个随机点的典型距离处与该最大值之比。

解答

（a） $\partial k/\partial\ell = \rho^2 e^{-\rho^2/2}/\ell$ 。对 $\rho^2 e^{-\rho^2/2}$ 求导得 $(2\rho - \rho^3)e^{-\rho^2/2}$ ，它在 $\rho = \sqrt{2}$ 处为零，此时该因子等于 $2/e \approx 0.74$ 。一对点相距约 1.4 个长度尺度时，为长度尺度提供的信息最多。（b）典型距离为 $\sqrt{50/6} \approx 2.89$ （第 3.1.2 节）。当 $\ell = 0.5$ 时， $\rho = 5.77$ ， $\rho^2 e^{-\rho^2/2} = 33.3 \times e^{-16.7} \approx 2 \times 10^{-6}$ ，约为最大值的百万分之三。当 $\ell = 0.2\sqrt{50} \approx 1.41$ 时， $\rho = 2.04$ ，该因子为 $4.17 \times e^{-2.08} \approx 0.52$ ，约为最大值的 70%。在固定长度尺度下，几乎没有哪对点能提供可用的梯度；在缩放后的长度尺度下，典型的点对就能提供。

参考文献

Duvenaud, D., Lloyd, J., Grosse, R., Tenenbaum, J., and Ghahramani, Z. (2013). Structure Discovery in Nonparametric Regression through Compositional Kernel Search. Proceedings of the 30th International Conference on Machine Learning (ICML 2013). 引用于 §9.1
Garnett, R. (2023). Bayesian Optimization. Cambridge University Press.
Hvarfner, C., Hellsten, E. O., and Nardi, L. (2024). Vanilla Bayesian Optimization Performs Great in High Dimensions. International Conference on Machine Learning. 引用于 §9.5
Meta Platforms, Inc. (2026e). BoTorch CHANGELOG. GitHub. 软件引用于 §9.5
Meta Platforms, Inc. (2026h). BoTorch PairwiseGP source code pairwise_gp.py. GitHub. 软件引用于 §9.5
Meta Platforms, Inc. (2026k). botorch/models/utils/gpytorch_modules.py. GitHub. 软件引用于 §9.4 §9.5
Neal, R. M. (1996). Bayesian Learning for Neural Networks. Springer. 引用于 §9.2
Papenmeier, L., Poloczek, M., and Nardi, L. (2025b). Understanding High-Dimensional Bayesian Optimization. ICML 2025, PMLR 267:47902-47923. 引用于 §9.5
Petersen, K. B., and Pedersen, M. S. (2012). The Matrix Cookbook. Technical University of Denmark. 非同行评审引用于 §9.4
Rasmussen, C. E., and Williams, C. K. I. (2006). Gaussian Processes for Machine Learning. MIT Press. 引用于 §9.1 §9.2 §9.3 §9.4 §9.6
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). 引用于 §9.2 §9.4
Stein, M. L. (1999). Interpolation of Spatial Data: Some Theory for Kriging. Springer.
Xu, Z., Wang, H., Phillips, J. M., and Zhe, S. (2025b). Standard Gaussian Process is All You Need for High-Dimensional Bayesian Optimization. ICLR 2025 (oral). 引用于 §9.5