超参数边际似然优化

高斯过程（GP）提供了一种灵活的非参数方法来对函数进行建模，直接在函数空间上设定先验。协方差函数（或核函数）编码了对函数属性（如平滑性或周期性）的设想。GP回归允许在新的点 $x_*$ 处预测 $f_*$，并获得不确定性估计 $\mathbb{E}[f_* | \mathbf{y}]$ 和 $\text{Var}[f_* | \mathbf{y}]$.

然而，核函数本身（例如流行的径向基函数（RBF）核函数：$k(x, x') = \sigma_f^2 \exp\left(-\frac{\|x - x'\|^2}{2l^2}\right)$）包含参数：信号方差 $\sigma_f^2$ 和长度尺度 $l$。类似地，GP回归中常用的高斯似然 $p(y|f) = \mathcal{N}(y|f, \sigma_n^2)$ 引入了一个噪声方差参数 $\sigma_n^2$。这些参数统称为 $\boldsymbol{\theta} = \{ \sigma_f^2, l, \sigma_n^2 \}$，称为超参数。我们如何为它们选择合适的值呢？

不同于通常的有监督学习中我们可能使用交叉验证来调整超参数，GP的贝叶斯特性提供了一种更统一的方法：优化边际似然（也称为证据）。

边际似然：对函数进行积分

其核心思路是找到超参数值 $\boldsymbol{\theta}$，使得在给定输入 $\mathbf{X}$ 和超参数 $\boldsymbol{\theta}$ 的情况下，观测到实际数据 $\mathbf{y}$ 的概率达到最大。这个概率是通过对训练点上未知的函数值 $\mathbf{f}$ 进行积分（或边际化）来获得的：

$$p(\mathbf{y} | \mathbf{X}, \boldsymbol{\theta}) = \int p(\mathbf{y} | \mathbf{f}, \sigma_n^2) p(\mathbf{f} | \mathbf{X}, \boldsymbol{\theta}_{\text{kernel}}) d\mathbf{f}$$

这里，$\boldsymbol{\theta}_{\text{kernel}}$ 表示核超参数（例如 $l$ 和 $\sigma_f^2$），并且我们明确显示了似然项中的噪声方差 $\sigma_n^2$。

由于先验 $p(\mathbf{f} | \mathbf{X}, \boldsymbol{\theta}_{\text{kernel}})$ 和似然 $p(\mathbf{y} | \mathbf{f}, \sigma_n^2)$ 都是高斯分布（假设回归中使用高斯似然），这个积分可以解析求解。其结果是 $\mathbf{y}$ 本身的高斯分布：

$$\mathbf{y} | \mathbf{X}, \boldsymbol{\theta} \sim \mathcal{N}(\mathbf{0}, K_y) = \mathcal{N}(\mathbf{0}, K_f + \sigma_n^2 I)$$

其中 $K_f$ 是使用核参数在训练输入 $\mathbf{X}$ 处计算得到的 $N \times N$ 协方差矩阵，$K_{ij} = k(x_i, x_j | \boldsymbol{\theta}_{\text{kernel}}$，而 $I$ 是单位矩阵。$K_y = K_f + \sigma_n^2 I$ 是观测数据点 $\mathbf{y}$ 的协方差矩阵，它结合了函数变异性（来自 $K_f$）和观测噪声（来自 $\sigma_n^2 I$）。

用于优化的对数边际似然

为了数值稳定性和计算便利，我们通常使用对数边际似然：

$$\log p(\mathbf{y} | \mathbf{X}, \boldsymbol{\theta}) = -\frac{1}{2} \mathbf{y}^T K_y^{-1} \mathbf{y} - \frac{1}{2} \log |K_y| - \frac{N}{2} \log(2\pi)$$

让我们分解一下这些重要项：

1. 数据拟合项： $-\frac{1}{2} \mathbf{y}^T K_y^{-1} \mathbf{y}$ 此项衡量模型解释观测数据 $\mathbf{y}$ 的效果。它鼓励超参数设置，使得在由协方差 $K_y$ 定义的高斯分布下，观测输出 $\mathbf{y}$ 具有较高概率。如果基于 $K_y$ 的预测与实际的 $\mathbf{y}$ 相差较大，则此项会变成一个较大的负数（表示惩罚）。

2. 复杂性惩罚项： $-\frac{1}{2} \log |K_y|$ 此项涉及协方差矩阵 $K_y$ 的行列式的对数。行列式 $|K_y|$ 衡量分布所覆盖的“体积”。导致点之间高度相关（例如，长度尺度非常大）或方差非常大（$\sigma_f^2$ 或 $\sigma_n^2$ 很大）的核函数会产生更大的行列式，从而导致此项更大的惩罚。此项起到了自动“奥卡姆剃刀”的作用，惩罚过于复杂的模型（例如，那些过于灵活或假设过多噪声的模型）。

3. 归一化常数： $-\frac{N}{2} \log(2\pi)$ 此项对于超参数 $\boldsymbol{\theta}$ 是常数，在优化过程中可以忽略。

最大化对数边际似然会促成一种权衡：找到能够使模型很好地拟合数据（第一项），同时避免不必要的复杂性（第二项）的超参数。

优化步骤

目标是找到使对数边际似然最大化的超参数 $\boldsymbol{\theta}^*$：

$$\boldsymbol{\theta}^* = \arg \max_{\boldsymbol{\theta}} \log p(\mathbf{y} | \mathbf{X}, \boldsymbol{\theta})$$

由于对数边际似然函数通常对超参数 $\boldsymbol{\theta}$（包括核参数如 $l, \sigma_f^2$ 和噪声方差 $\sigma_n^2$）可导，我们可以使用标准的基于梯度的优化算法。常用选项包括：

• 共轭梯度法
• L-BFGS-B（带边界的有限内存Broyden–Fletcher–Goldfarb–Shanno，对于像方差和长度尺度这样必须为正的受限超参数很有用）。

要使用这些方法，我们需要对数边际似然对每个超参数 $\theta_j \in \boldsymbol{\theta}$ 的偏导数。该导数具有相对简洁的解析形式：

$$\frac{\partial}{\partial \theta_j} \log p(\mathbf{y} | \mathbf{X}, \boldsymbol{\theta}) = \frac{1}{2} \mathbf{y}^T K_y^{-1} \frac{\partial K_y}{\partial \theta_j} K_y^{-1} \mathbf{y} - \frac{1}{2} \text{Tr}\left( K_y^{-1} \frac{\partial K_y}{\partial \theta_j} \right)$$

这里，$\frac{\partial K_y}{\partial \theta_j}$ 是包含核函数（如果 $\theta_j = \sigma_n^2$ 则加上噪声项）对超参数 $\theta_j$ 的逐元素偏导数的矩阵。对于标准核函数，这些导数通常易于计算。例如，如果 $K_y = K_f + \sigma_n^2 I$，并且 $\theta_j$ 是像长度尺度 $l$ 这样的核参数，那么 $\frac{\partial K_y}{\partial \theta_j} = \frac{\partial K_f}{\partial l}$。如果 $\theta_j = \sigma_n^2$，那么 $\frac{\partial K_y}{\partial \theta_j} = I$。

现代的GP库（如GPy、GPflow、scikit-learn的 GaussianProcessRegressor、GPytorch）会自动计算这些梯度并为您执行优化。通常您只需定义核结构并提供数据。

# 使用类似 scikit-learn 接口的示例
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import RBF, ConstantKernel as C, WhiteKernel

# 定义带有初始超参数的核函数（这些参数将被优化）
# 核函数: k(x, x') = C(常数值) * RBF(长度尺度) + WhiteKernel(噪声水平)
# 我们通常优化对数变换后的参数以获得稳定性
kernel = C(1.0) * RBF(length_scale=1.0) + WhiteKernel(noise_level=1.0)

# 创建 GP 模型
# n_restarts_optimizer 通过从不同点重新开始有助于找到更好的最优解
gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10, random_state=42)

# 拟合模型：此步骤执行对数边际似然最大化
# X_train, y_train 是您的训练数据
# gp.fit(X_train, y_train)

# 拟合后，优化的超参数会被保存
# print(f"优化后的核函数: {gp.kernel_}")
# print(f"对数边际似然: {gp.log_marginal_likelihood(gp.kernel_.theta)}")

# 示例输出（示意性）：
# 优化后的核函数: 1.5**2 * RBF(length_scale=0.8) + WhiteKernel(noise_level=0.1)
# 对数边际似然: -15.23

挑战与考量

尽管功能强大，但优化边际似然并非没有挑战：

1. 计算成本： 主要的瓶颈在于计算 $K_y^{-1}$ 和 $\log |K_y|$，这两项的复杂度通常都为 $O(N^3)$，其中 $N$ 是训练点的数量。这使得大型数据集的直接优化变得困难（从而促使我们稍后讨论可扩展的近似方法）。梯度计算也涉及矩阵乘法，这增加了成本，但每次迭代仍主要由 $O(N^3)$ 项决定。
2. 非凸性： 对数边际似然曲面可能不是凸的，这意味着优化算法可能会收敛到局部最优而不是全局最大值。这就是为什么实际实现中通常会从不同的初始超参数值进行多次重启（如上例中的 n_restarts_optimizer）。
3. 初始化： 初始超参数值的选择可能会显著影响优化结果，尤其是在复杂的曲面上。基于领域知识或数据属性（如输入和输出的范围）的合理初始猜测会有帮助。

让我们看看在一个简单的1D回归问题中，对于RBF核函数，对数边际似然（LML）如何随长度尺度（$l$）和信号方差（$\sigma_f^2$）这两个超参数（为简化起见，噪声 $\sigma_n^2$ 固定）的变化而变化。

等高线图显示了RBF核函数长度尺度和信号方差不同组合下的对数边际似然（LML）值。优化算法在此曲面上寻找峰值（最高的LML值）。这种可能复杂的形状表明为什么可能需要多次重启才能找到全局最优解。

通过找到使观测数据概率最大化的超参数，边际似然优化提供了一种有原则且通常有效的方法来调整GP模型，将模型选择直接融入贝叶斯推理框架中。这与交叉验证等方法不同，它提供了一种更节省数据的方法，尤其是在数据稀缺时。一旦找到最优超参数 $\boldsymbol{\theta}^*$，它们将被代入GP预测方程，用于在新数据上进行预测。