高斯过程（Gaussian Process, GP）的边际似然优化与超参数学习过程

1. 题目描述

高斯过程是一种强大的非参数贝叶斯模型，常用于回归与分类任务。它通过定义函数空间上的先验分布，并利用观测数据更新为后验分布进行预测。高斯过程的核心特性完全由其均值函数和协方差函数（核函数）决定。在实际应用中，协方差函数通常包含超参数（如长度尺度、方差等），这些超参数显著影响模型性能。本题目将详细讲解如何通过最大化边际似然（Marginal Likelihood）来优化高斯过程的超参数，包括边际似然公式的推导、梯度计算以及基于梯度的优化过程。

2. 高斯过程回归回顾

设训练集为 \(\mathcal{D} = \{ (\mathbf{x}_i, y_i) \}_{i=1}^n\)，其中 \(\mathbf{x}_i \in \mathbb{R}^d\)，\(y_i \in \mathbb{R}\)。高斯过程回归模型假设观测值由潜在函数 \(f(\mathbf{x})\) 加上高斯噪声生成：

\[y = f(\mathbf{x}) + \epsilon, \quad \epsilon \sim \mathcal{N}(0, \sigma_n^2) \]

高斯过程将 \(f(\mathbf{x})\) 定义为服从高斯过程的随机函数：

\[f(\mathbf{x}) \sim \mathcal{GP}(m(\mathbf{x}), k(\mathbf{x}, \mathbf{x}')) \]

其中 \(m(\mathbf{x})\) 为均值函数（常设为0），\(k(\mathbf{x}, \mathbf{x}')\) 为协方差函数（如径向基函数核 RBF）：

\[k(\mathbf{x}, \mathbf{x}') = \sigma_f^2 \exp\left(-\frac{\|\mathbf{x} - \mathbf{x}'\|^2}{2\ell^2}\right) \]

超参数集合记为 \(\theta = \{\sigma_f^2, \ell, \sigma_n^2\}\)，分别表示信号方差、长度尺度和噪声方差。

3. 边际似然的定义与推导

3.1 边际似然公式

在贝叶斯框架中，边际似然（或证据） 是在给定超参数下观测数据的概率，通过对潜在函数积分得到：

\[p(\mathbf{y} \mid \mathbf{X}, \theta) = \int p(\mathbf{y} \mid \mathbf{f}, \theta) \, p(\mathbf{f} \mid \mathbf{X}, \theta) \, d\mathbf{f} \]

由于噪声独立同分布，似然项为：

\[p(\mathbf{y} \mid \mathbf{f}, \theta) = \mathcal{N}(\mathbf{y} \mid \mathbf{f}, \sigma_n^2 I) \]

先验项为：

\[p(\mathbf{f} \mid \mathbf{X}, \theta) = \mathcal{N}(\mathbf{f} \mid \mathbf{0}, K) \]

其中 \(K\) 是核矩阵，元素 \(K_{ij} = k(\mathbf{x}_i, \mathbf{x}_j)\)。两个高斯分布的卷积仍为高斯分布，可得：

\[p(\mathbf{y} \mid \mathbf{X}, \theta) = \mathcal{N}(\mathbf{y} \mid \mathbf{0}, K + \sigma_n^2 I) \]

记 \(C = K + \sigma_n^2 I\)，则边际似然的对数形式为：

\[\log p(\mathbf{y} \mid \mathbf{X}, \theta) = -\frac{1}{2} \mathbf{y}^\top C^{-1} \mathbf{y} - \frac{1}{2} \log |C| - \frac{n}{2} \log 2\pi \]

4. 边际似然的三项解释

对数边际似然可分解为三项，每项有直观含义：

1. 数据拟合项 \(-\frac{1}{2} \mathbf{y}^\top C^{-1} \mathbf{y}\)：衡量模型对观测数据的拟合程度（越小表示拟合越好）。
2. 复杂度惩罚项 \(-\frac{1}{2} \log |C|\)：惩罚模型复杂度，防止过拟合（核矩阵行列式随模型复杂度增加而增大）。
3. 常数项 \(-\frac{n}{2} \log 2\pi\)：与超参数无关。

优化超参数时，需在拟合数据与控制模型复杂度之间权衡。

5. 对数边际似然的梯度计算

为了使用梯度下降法等优化算法，需计算对数边际似然对超参数 \(\theta_j\) 的梯度。设 \(C = K + \sigma_n^2 I\)，有：

\[\frac{\partial \log p(\mathbf{y} \mid \mathbf{X}, \theta)}{\partial \theta_j} = \frac{1}{2} \mathbf{y}^\top C^{-1} \frac{\partial C}{\partial \theta_j} C^{-1} \mathbf{y} - \frac{1}{2} \text{tr}\left( C^{-1} \frac{\partial C}{\partial \theta_j} \right) \]

推导过程涉及矩阵微分和迹的性质：

• 第一项来自数据拟合项的微分：\(\frac{\partial}{\partial \theta_j} \left( -\frac{1}{2} \mathbf{y}^\top C^{-1} \mathbf{y} \right) = \frac{1}{2} \mathbf{y}^\top C^{-1} \frac{\partial C}{\partial \theta_j} C^{-1} \mathbf{y}\)。
• 第二项来自复杂度项的微分：\(\frac{\partial}{\partial \theta_j} \left( -\frac{1}{2} \log |C| \right) = -\frac{1}{2} \text{tr}\left( C^{-1} \frac{\partial C}{\partial \theta_j} \right)\)。

5.1 梯度计算示例（RBF核）

对于 RBF 核超参数：

• 信号方差 \(\sigma_f^2\)：\(\frac{\partial C}{\partial \sigma_f^2} = \frac{\partial K}{\partial \sigma_f^2} = \frac{K}{\sigma_f^2}\)（因 \(K = \sigma_f^2 \cdot \text{RBF部分}\)）。
• 长度尺度 \(\ell\)：\(\frac{\partial K_{ij}}{\partial \ell} = K_{ij} \cdot \frac{\|\mathbf{x}_i - \mathbf{x}_j\|^2}{\ell^3}\)。
• 噪声方差 \(\sigma_n^2\)：\(\frac{\partial C}{\partial \sigma_n^2} = I\)。

将这些偏导代入梯度公式即可计算。

6. 超参数优化步骤

步骤1：初始化超参数

• 设置初始值，如 \(\ell = 1.0\), \(\sigma_f^2 = 1.0\), \(\sigma_n^2 = 0.1\)。
• 定义核函数（如 RBF）并计算核矩阵 \(K\)。

步骤2：计算对数边际似然及其梯度

1. 构建矩阵 \(C = K + \sigma_n^2 I\)。
2. 计算 \(C\) 的逆矩阵 \(C^{-1}\) 和行列式 \(|C|\)（常用 Cholesky 分解提高数值稳定性）。
3. 代入公式计算对数边际似然值。
4. 计算 \(C\) 对每个超参数的偏导 \(\frac{\partial C}{\partial \theta_j}\)。
5. 利用梯度公式计算梯度向量。

步骤3：使用优化算法更新超参数

• 采用梯度上升法或更高效的优化器（如 L-BFGS、共轭梯度法）最大化对数边际似然。
• 更新规则（梯度上升）：\(\theta_j \leftarrow \theta_j + \eta \cdot \frac{\partial \log p(\mathbf{y} \mid \mathbf{X}, \theta)}{\partial \theta_j}\)，其中 \(\eta\) 为学习率。
• 重复步骤2-3直至收敛（梯度接近零或似然变化很小）。

步骤4：验证与预测

• 用优化后的超参数重新训练高斯过程模型。
• 进行后验预测（计算测试点的均值和方差）。

7. 关键注意事项

1. 数值稳定性：直接计算 \(C^{-1}\) 和 \(\log |C|\) 易出现数值错误，推荐使用 Cholesky 分解 \(C = LL^\top\)，然后通过回代求解线性系统并利用 \(\log |C| = 2 \sum_i \log L_{ii}\) 计算行列式。
2. 局部最优：边际似然可能非凸，存在多个局部极大值。可尝试多次随机初始化或使用全局优化方法（如贝叶斯优化）。
3. 超参数含义：
   • 长度尺度 \(\ell\) 控制函数平滑度（越大越平滑）。
   • 噪声方差 \(\sigma_n^2\) 反映数据噪声水平。
   • 信号方差 \(\sigma_f^2\) 控制函数幅度。

8. 示例代码框架（Python伪代码）

import numpy as np
from scipy.optimize import minimize

def rbf_kernel(X1, X2, l, sigma_f):
    dist = np.sum(X1**2, axis=1).reshape(-1,1) + np.sum(X2**2, axis=1) - 2 * X1 @ X2.T
    return sigma_f**2 * np.exp(-0.5 * dist / l**2)

def marginal_likelihood(theta, X, y):
    l, sigma_f, sigma_n = theta
    K = rbf_kernel(X, X, l, sigma_f)
    C = K + sigma_n**2 * np.eye(len(X))
    L = np.linalg.cholesky(C)  # Cholesky分解
    alpha = np.linalg.solve(L.T, np.linalg.solve(L, y))
    log_likelihood = -0.5 * y.T @ alpha - np.sum(np.log(np.diag(L))) - len(X)/2 * np.log(2*np.pi)
    return -log_likelihood  # 负值用于最小化

# 优化超参数
theta0 = [1.0, 1.0, 0.1]  # 初始猜测
result = minimize(marginal_likelihood, theta0, args=(X_train, y_train), method='L-BFGS-B')
l_opt, sigma_f_opt, sigma_n_opt = result.x

9. 总结

通过最大化边际似然优化高斯过程的超参数，本质上是寻找最符合数据统计特性的核函数参数。该过程平衡了拟合优度与模型复杂度，并通过梯度计算实现高效优化。掌握这一方法对于实际应用高斯过程至关重要，因为它直接决定了模型的预测性能与泛化能力。