高斯过程回归（Gaussian Process Regression）的协方差矩阵计算与预测过程

好的，我将为你讲解高斯过程回归中协方差矩阵的计算与基于该矩阵的预测过程。这不仅仅是之前GPR主题的简单重复，而是深入其计算核心，详细拆解如何从数据和协方差函数（核函数）出发，构建预测所需的全部矩阵运算，最终得到预测分布。

题目描述

假设我们有一个回归问题：给定一组训练数据点 (X, y)，其中 X 是 (n×d) 的特征矩阵（n 个样本，每个样本 d 维），y 是 (n×1) 的观测值向量。高斯过程回归模型不对函数形式做参数化假设，而是将函数 f(x) 视为一个从高斯过程中抽取的随机函数。我们的目标是在给定训练数据后，对任意新的测试输入 X*（(m×d) 的矩阵，m 个测试点），预测其对应的函数值 f*（或含噪声的观测值 y*）的联合分布。

核心问题：如何具体地计算训练样本之间、训练与测试样本之间、测试样本之间的协方差矩阵？如何利用这些矩阵和贝叶斯公式，推导出预测均值向量和预测协方差矩阵的闭式解？

解题过程

步骤 1：理解高斯过程（GP）先验

一个高斯过程完全由其均值函数 m(x) 和协方差函数（核函数）k(x, x’) 定义。为简化，通常假设先验均值函数为0（数据可先中心化处理）。因此，对于任意一组输入点（例如训练点 X 和测试点 X* 的集合），其对应的函数值服从联合高斯分布：

[ f; f* ] ~ N( 0, [ K(X, X)   K(X, X*);
                    K(X*, X)  K(X*, X*) ] )

这里的 K(A, B) 表示一个矩阵，其第 i, j 个元素是 k(a_i, b_j)，其中 a_i 和 b_j 分别是矩阵 A 和 B 的行向量。

步骤 2：引入观测噪声模型

在实际中，我们观测到的是含噪声的 y，通常假设为：y = f(x) + ε，其中 ε ~ N(0, σ_n^2)，且各观测噪声独立。因此，训练观测值 y 的协方差 不仅来自函数值的协方差 K(X, X)，还要加上噪声方差项。
定义：

• K = K(X, X)：一个 (n×n) 的矩阵，表示训练点自身之间的先验协方差。
• 那么 y 的协方差矩阵为 Cov(y) = K + σ_n^2 * I，其中 I 是 (n×n) 的单位矩阵。记 Ky = K + σ_n^2 * I。

步骤 3：构建训练-测试联合分布

现在，将训练观测值 y 和测试点的函数值 f*（我们关心的是无噪声的函数值）组成一个联合分布。根据GP的定义和噪声假设：

[ y;  f* ] ~ N( 0, [ Ky                K(X, X*);
                     K(X*, X)          K(X*, X*) ] )

这里：

• K(X, X*)：一个 (n×m) 的矩阵，我通常记作 K_*，即 K_*[i, j] = k(x_i, x*_j)。
• K(X*, X)：是 K_* 的转置，一个 (m×n) 的矩阵，记作 K_*^T。
• K(X*, X*)：一个 (m×m) 的矩阵，记作 K_**，表示测试点自身之间的先验协方差。

所以联合分布简写为：

[ y;  f* ] ~ N( 0, [ Ky     K_*;
                     K_*^T  K_** ] )

步骤 4：应用条件分布公式进行预测

我们想要的是在已知训练数据 y 的条件下，f* 的分布，即 p(f* | X, y, X*)。对于联合高斯分布，条件分布有标准的闭式解。

设联合分布为：

[ u; v ] ~ N( [ μ_u; μ_v ], [ Σ_uu  Σ_uv;
                               Σ_vu  Σ_vv ] )

则条件分布 p(v | u) 也是高斯的，其均值和协方差为：

μ_{v|u} = μ_v + Σ_vu * Σ_uu^{-1} * (u - μ_u)
Σ_{v|u} = Σ_vv - Σ_vu * Σ_uu^{-1} * Σ_uv

在我们的设定中：

• u = y, v = f*
• μ_u = 0, μ_v = 0
• Σ_uu = Ky
• Σ_uv = Σ_vu^T = K_*
• Σ_vv = K_**

代入公式，得到高斯过程回归的预测方程：

预测均值向量（长度为 m）：

μ_f* = K_*^T * Ky^{-1} * y

这可以理解为：预测是训练观测值 y 的线性组合，权重由核函数度量的相似性（K_*^T）和训练点自身的协方差结构（Ky^{-1}）共同决定。

预测协方差矩阵（m×m）：

Σ_f* = K_** - K_*^T * Ky^{-1} * K_*

第一项 K_** 是测试点自身的先验不确定性。减去第二项 K_*^T * Ky^{-1} * K_* 代表了由于观察到训练数据 y 而减少的不确定性（信息增益）。这个矩阵的对角线元素就是每个测试点的预测方差。

步骤 5：含噪声的测试观测值预测

如果我们想预测的是带噪声的测试观测值 y*（其中 y* = f* + ε*, ε* ~ N(0, σ_n^2)），那么只需在预测函数值 f* 的协方差上加上噪声方差项：

μ_y* = μ_f*  （均值不变）
Σ_y* = Σ_f* + σ_n^2 * I

步骤 6：计算流程总结与复杂度分析

1. 计算核矩阵：根据选择的核函数（如径向基函数RBF），计算：
   • K = k(X, X)：O(n^2 * d) 计算，O(n^2) 存储。
   • K_* = k(X, X*)：O(n * m * d)。
   • K_** = k(X*, X*)：O(m^2 * d)。
2. 构造 Ky：Ky = K + σ_n^2 * I。
3. 求解线性系统：这是计算的核心和瓶颈。需要计算 Ky^{-1} * y 和 Ky^{-1} * K_*。
   • 不直接求逆：通常通过求解线性方程组来完成。例如，对 Ky 进行Cholesky分解（因为 Ky 是对称正定的）：Ky = L * L^T，其中 L 是下三角矩阵。
   • 然后：
     • 解 L * z = y 得 z，再解 L^T * α = z，得到的 α 就是 Ky^{-1} * y。
     • 解 L * A = K_* 得 A，再解 L^T * B = A，得到的 B 就是 Ky^{-1} * K_*。
   • Cholesky分解的复杂度是 O(n^3)，后续回代求解是 O(n^2) 和 O(n^2 * m)。这是GPR应用于大规模数据的主要瓶颈。
4. 计算预测均值和方差：
   • μ_f* = K_*^T * α （O(n * m)）
   • v = L^{-1} * K_* （实际上已在上一步得到 B）
   • Σ_f* = K_** - v^T * v （利用 v^T * v = (K_*^T * L^{-T}) * (L^{-1} * K_*) = K_*^T * Ky^{-1} * K_*，计算复杂度 O(m^2 * n)）
   • 预测方差是 Σ_f* 的对角线元素。如果只需要预测方差（而不是完整的协方差矩阵），可以更高效地计算，避免 O(m^2) 的存储和计算。

关键点与意义

• 协方差矩阵 K, K_*, K_** 是连接所有点的“关系网”，核函数定义了如何衡量两个输入点之间的相似性。
• 预测均值是训练值的加权平均，权重由测试点与所有训练点的相似性，并通过训练点之间的相互关系进行校正（Ky^{-1}）后决定。
• 预测协方差不仅给出了预测的不确定性（方差），还给出了不同测试点预测值之间的相关性（协方差）。这是GPR相比其他回归方法（如仅给出点估计的模型）的一大优势。
• 计算核心在于稳定高效地求解以 Ky 为系数矩阵的线性系统，通常通过矩阵分解实现。

这个过程清晰地展示了高斯过程回归如何从贝叶斯框架出发，通过严谨的矩阵运算，将先验分布更新为后验预测分布，从而提供完整的概率预测。