高斯过程回归中核函数选择与超参数边际似然优化

高斯过程回归 (Gaussian Process Regression, GPR) 的核心，在于通过一个核函数（也称为协方差函数）来定义数据点之间的相似性。核函数的选择直接决定了模型的预测能力和泛化性能。而核函数中的超参数，则需要通过优化算法来学习。本指南将手把手教你如何选择和设计核函数，并使用边际似然作为目标函数来优化其超参数。

1. 理解核函数的核心作用

核函数 $k(\mathbf{x}_i, \mathbf{x}_j)$ 衡量两个输入向量 $\mathbf{x}_i$ 和 $\mathbf{x}_j$ 的相似度。它将输入空间映射到一个高维的特征空间，使得在原始空间中线性不可分的问题，在新的特征空间中变得线性可分。

核心结论：核函数决定了模型所能描述的函数空间的“形状”。选择错误的核函数，模型可能完全无法捕捉数据的真实规律。

2. 选择合适的核函数

核函数的选择需要结合你对数据生成过程的先验知识。

1. 识别数据特性。观察你的训练数据。它是否有周期性？是否存在平滑的趋势？数据噪声是大还是小？这些观察是选择核函数的依据。

2. 从标准核函数库中选取。以下是几种常用的基础核函数：

   • 平方指数核 (Squared Exponential Kernel)：最常用的核。它生成非常平滑的函数。适用于绝大多数没有特殊结构的连续数据。其公式中的长度尺度参数 $\ell$ 控制了函数的“摆动”频率。
   • 有理二次核 (Rational Quadratic Kernel)：可以看作是多个不同长度尺度的平方指数核的叠加。它比平方指数核更灵活，能处理具有多尺度变化特性的数据。
   • Matérn 核：通过一个参数 $\nu$ 来控制生成函数的平滑度。$\nu = 1/2$ 时生成粗糙的函数（等价于指数核），$\nu \to \infty$ 时趋于平方指数核。在物理或地理空间数据建模中很常见。
   • 周期核 (Periodic Kernel)：用于捕捉数据中的周期性模式，例如季节性时间序列。
   • 线性核 (Linear Kernel)：生成通过原点的线性函数。

3. 组合核函数以构建复杂结构。将简单的核函数进行组合，可以构建出表达能力更强的核。常用组合方式有：

   • 相加：例如，将平方指数核与周期核相加，可以描述一个有长期平滑趋势并带有周期性波动的数据序列。
   • 相乘：例如，将平方指数核与线性核相乘，可以描述一个幅度随输入线性增长的振荡函数。
   • 缩放：乘以一个常数（信号方差 $\sigma_f^2$）来控制函数的整体幅度。

关键动作：选择一个或组合多个基础核函数，形成一个复合核 $k$。

3. 理解边际似然优化的目标

高斯过程的超参数（例如核函数中的长度尺度 $\ell$、信号方差 $\sigma_f^2$，以及噪声方差 $\sigma_n^2$）需要通过数据来学习。最常用的方法是最大化边际似然（也称为证据）。

边际似然 $p(\mathbf{y}|\mathbf{X}, \theta)$ 表示在给定输入数据 $\mathbf{X}$ 和一组超参数 $\theta$ 的条件下，观测到目标值 $\mathbf{y}$ 的概率。最大化这个概率，意味着寻找一组超参数 $\theta$，使得当前模型（由核函数和超参数定义）能够最“自然”地解释观测到的数据。

通俗理解：边际似然衡量了模型与当前数据的整体匹配程度。优化它就是在自动进行模型选择——匹配度差的复杂模型会被惩罚，从而防止过拟合。

4. 执行边际似然优化（实操步骤）

这里以 Python 的 scikit-learn 库为例，演示如何实现。假设你已经准备好了训练数据 X_train 和 y_train。

1. 定义模型和初始核。导入 GaussianProcessRegressor 和核函数。构建一个复合核。例如，一个带有噪声的平方指数核：

   from sklearn.gaussian_process import GaussianProcessRegressor
   from sklearn.gaussian_process.kernels import RBF, WhiteKernel
   import numpy as np
   
   # 定义核函数：信号方差*平方指数核 + 噪声方差*白噪声核
   kernel = 1.0 * RBF(length_scale=1.0, length_scale_bounds=(1e-2, 1e2)) + WhiteKernel(noise_level=0.1, noise_level_bounds=(1e-5, 1e1))

   length_scale_bounds 和 noise_level_bounds 为超参数指定了优化时的搜索范围。

2. 创建高斯过程回归器实例。设置优化器和迭代次数。

   gp = GaussianProcessRegressor(kernel=kernel, n_restarts_optimizer=10, random_state=42)

   n_restarts_optimizer=10 表示从10个不同的随机初始点开始优化，以降低陷入局部最优的风险。

3. 执行拟合（优化）。调用 fit 方法。此过程会自动执行边际似然优化，寻找最优的超参数。

   gp.fit(X_train, y_train)

   执行后，优化后的超参数存储在 gp.kernel_ 属性中。

4. 检查优化结果。打印优化后的核函数及其超参数。

   print(“优化后的核函数：”, gp.kernel_)

   你会看到类似 1.0**2 * RBF(length_scale=0.5) + WhiteKernel(noise_level=0.01) 的输出，这些数值就是优化得到的最佳超参数。

5. 进行预测与评估。使用优化后的模型进行预测。predict 方法不仅给出预测均值，还可以通过设置 return_std=True 返回预测的不确定性（标准差）。

   y_pred, y_std = gp.predict(X_test, return_std=True)

   预测的标准差 $\sigma(\mathbf{x}_*)$ 能够直观地反映模型在未见过的数据点上的置信度。

5. 手动计算与分析边际似然

为了更深入地理解，你可以手动计算和查看边际似然的值及其梯度。

1. 定义边际似然。在高斯过程回归中，给定超参数 $\theta$，对数边际似然的公式为：
   $$\log p(\mathbf{y}|\mathbf{X}, \theta) = -\frac{1}{2}\mathbf{y}^\top K_\theta^{-1}\mathbf{y} - \frac{1}{2}\log|K_\theta| - \frac{n}{2}\log(2\pi)$$
   其中 $K_\theta = K(\mathbf{X}, \mathbf{X}; \theta) + \sigma_n^2 I$ 是训练数据的核矩阵（包含噪声项）。

2. 理解优化目标。优化器通过梯度上升法最大化对数边际似然。梯度 $\frac{\partial \log p(\mathbf{y}|\mathbf{X}, \theta)}{\partial \theta_j}$ 的计算是自动的。

3. 查看优化过程中的对数边际似然值（高级）。scikit-learn 的优化器内部会计算这个值，但通常不直接暴露给用户。你可以通过设置 optimizer=None 禁用内置优化，然后循环手动计算不同超参数下的对数边际似然来绘制其曲线。

   # 示例：手动计算不同长度尺度下的对数边际似然
   log_likelihoods = []
   length_scales = np.logspace(-1, 1, 50)
   for ls in length_scales:
       temp_kernel = 1.0 * RBF(length_scale=ls) + WhiteKernel(noise_level=0.1)
       temp_gp = GaussianProcessRegressor(kernel=temp_kernel, optimizer=None)
       temp_gp.fit(X_train, y_train)
       # log_marginal_likelihood_value_ 属性存储了当前核参数下的对数边际似然
       log_likelihoods.append(temp_gp.log_marginal_likelihood_value_)

   绘制 length_scales 与 log_likelihoods 的曲线，你能清晰地看到优化器是如何寻找使对数边际似然达到最大值的那个长度尺度的。