高斯过程分类（Gaussian Process Classification, GPC）的拉普拉斯近似推断过程

我将为您详细讲解高斯过程分类中的拉普拉斯近似推断过程。这是一个非常经典的概率机器学习算法，用于解决分类问题。

一、问题描述与背景

高斯过程分类是高斯过程回归在分类任务上的推广。但分类问题与回归问题有一个根本区别：分类输出是离散的（如0/1），而高斯过程模型本身是定义在连续值函数上的。

核心挑战：在回归中，我们假设观测数据y来自一个带有高斯噪声的潜在函数f(x)：y = f(x) + ε。这使得y的分布也是高斯的，我们可以直接计算边缘似然。但在二分类中，观测标签y ∈ {0, 1}，我们通常通过一个非线性函数（如逻辑函数）将潜在函数f(x)映射到[0,1]区间：

p(y=1 | f(x)) = σ(f(x))

其中σ(·)是逻辑函数（或其他单调函数，如标准正态的累积分布函数Φ）。这里，f(x)仍然被建模为一个高斯过程。

关键问题：这个非线性连接使得模型的后验分布 p(f | X, y) 不再是高斯的（f是所有潜在函数值的向量）。这使得精确推断变得不可能。我们需要一个高效的近似方法来处理这个非高斯的后验。

拉普拉斯近似就是为了解决这个问题而引入的。其核心思想是：用一个高斯分布来近似这个非高斯后验。这个高斯分布的均值被设在真实后验分布的众数（最大值点）处，协方差由该点处对数后验的负Hessian矩阵的逆给出。

二、模型定义与贝叶斯框架

1. 模型层级

对于一个二分类问题，有数据集 D = {(x_i, y_i)}_{i=1}^n，其中 y_i ∈ {0, 1}。

• 先验：f ~ GP(0, k(x, x’))。在有限的n个数据点上，f = [f(x_1), …, f(x_n)]^T 服从多元高斯分布：

  p(f | X) = N(f | 0, K)
  

  这里K是n×n的协方差矩阵，K_{ij} = k(x_i, x_j)，k是核函数（如RBF核）。
• 似然：观测y_i是条件独立的，给定f_i = f(x_i)：

  p(y_i | f_i) = σ(y_i * f_i)  (注意：为简化，有时定义 t_i = 2y_i - 1 ∈ {-1, 1})
  

  更常见的是使用伯努利分布：

  p(y_i | f_i) = σ(f_i)^{y_i} * (1 - σ(f_i))^{1-y_i}
  

• 后验：根据贝叶斯定理：

  p(f | X, y) = p(y | f) * p(f | X) / p(y | X)
  

  分母p(y|X)是边缘似然，是我们进行模型比较（超参数优化）的目标，但计算它需要积分∫ p(y|f) p(f|X) df，这在高维空间中难以处理。

2. 拉普拉斯近似的目标

我们想用一个高斯分布q(f | X, y) = N(f | \hat{f}, A^{-1})来近似真实后验p(f | X, y)，其中：

• \hat{f}是后验分布的众数（mode），即最大化p(f|X, y)的点，等价于最大化对数联合分布：

  \hat{f} = argmax_f log p(y|f) + log p(f|X)
  

• A是在\hat{f}处对数后验的负Hessian矩阵：A = -∇∇ log p(f | X, y) |_{f=\hat{f}}。由于我们使用高斯近似，A^{-1}就是近似后验的协方差矩阵。

三、拉普拉斯近似的具体步骤

我们可以将整个过程分解为以下几步：

步骤1：定义目标函数（对数后验）

我们最大化对数联合分布（忽略常数）：

Ψ(f) = log p(f|X, y) = log p(y|f) + log p(f|X) + const.
          = ∑_{i=1}^n log p(y_i | f_i) - (1/2) f^T K^{-1} f - (1/2) log|K| + const.

注意：log p(f|X)来自于多元高斯分布的对数概率密度。

步骤2：寻找后验众数 \hat{f}

我们需要找到使Ψ(f)最大化的\hat{f}。由于Ψ(f)关于f是凹的（对于逻辑似然和GP先验），我们可以使用牛顿法进行迭代优化。

牛顿法更新公式：f^{new} = f^{old} - (∇∇Ψ)^{-1} ∇Ψ

让我们计算梯度和Hessian：

• 梯度：∇Ψ(f) = ∇ log p(y|f) - K^{-1} f
  其中∇ log p(y|f)是一个n维向量，第i个分量为 ∂ log p(y_i | f_i)/∂f_i。对于逻辑似然p(y_i|f_i) = 1/(1+exp(-y_i f_i))（这里y_i ∈ {-1,1}），有：

  ∂ log p(y_i | f_i) / ∂f_i = t_i * π_i
  其中 t_i = (y_i+1)/2 ∈ {0,1} 的常见伯努利形式下，梯度为：∂/∂f_i = y_i - σ(f_i)
  更常用：设 π_i = p(y_i=1|f_i) = σ(f_i)，则 ∂ log p(y_i | f_i)/∂f_i = y_i - π_i
  

  所以∇ log p(y|f) = y - π，其中π = [π_1, …, π_n]^T。

• Hessian矩阵：∇∇Ψ(f) = ∇∇ log p(y|f) - K^{-1}
∇∇ log p(y|f)是一个对角矩阵，因为似然对数据点是独立的。第i个对角元素为：

  ∂² log p(y_i | f_i)/∂f_i² = - π_i (1 - π_i)
  

  因为对于逻辑函数，有σ’(z) = σ(z)(1-σ(z))。记W = diag(π_i (1-π_i))，则∇∇ log p(y|f) = -W。
  因此，∇∇Ψ(f) = -W - K^{-1}。

牛顿法迭代更新为：

f^{new} = f^{old} + (W + K^{-1})^{-1} (y - π - K^{-1} f^{old})

但通常为了数值稳定性，我们使用等价形式：

f^{new} = (K^{-1} + W)^{-1} W (f^{old} + W^{-1}(y - π))

迭代直到收敛，得到众数\hat{f}。

步骤3：计算近似后验的协方差

在收敛点\hat{f}处，我们计算负Hessian的逆作为近似后验的协方差矩阵：

A = -∇∇Ψ(f) |_{f=\hat{f}} = W + K^{-1}

其中W是在\hat{f}处计算的。那么，近似后验为：

q(f | X, y) = N(f | \hat{f}, A^{-1}) = N(f | \hat{f}, (W + K^{-1})^{-1})

更常见的写法是利用矩阵恒等式，写成：

q(f | X, y) = N(f | \hat{f}, (K^{-1} + W)^{-1})

或者等价地，A^{-1} = K - K (K + W^{-1})^{-1} K。

四、预测新样本

当我们有一个新测试点x_*时，我们想要预测其标签y_*的概率p(y_*=1 | x_*, X, y)。这需要两个步骤：

步骤1：计算潜在函数f_*的后验分布

在GP中，训练潜在变量f和测试点潜在变量f_* = f(x_*)的联合先验是高斯分布：

[ f; f_* ] ~ N(0, [ K   k_*; k_*^T  k_** ])

其中k_* = [k(x_1, x_*), …, k(x_n, x_*)]^T，k_** = k(x_*, x_*)。

由于我们已经用高斯分布q(f) = N(\hat{f}, Σ)近似了p(f|X,y)（其中Σ = (W + K^{-1})^{-1}），那么f_*的条件后验p(f_* | x_*, X, y) ≈ ∫ p(f_* | f, x_*) q(f) df也是高斯的。计算高斯条件分布得到：

均值： μ_* = k_*^T K^{-1} \hat{f}   (更常见的是 μ_* = k_*^T (K + W^{-1})^{-1} W^{-1} (y - π)? 我们推导一下)

实际上，更直接的公式是：利用q(f)的均值和协方差，f_*的后验均值为：

E[f_*] = k_*^T ∇ log p(y | \hat{f})   (一种形式)

经过推导，常用的公式为：

μ_* = k_*^T K^{-1} \hat{f} = k_*^T (K + W^{-1})^{-1} W^{-1} (y - π)   （不一定精确，但思想类似）

为了避免混淆，我们使用标准结果：μ_* = k_*^T (K + W^{-1})^{-1} W^{-1} \tilde{y}，其中\tilde{y} = W \hat{f} + (y - π)。

方差为：

σ^2_* = k_** - k_*^T (K + W^{-1})^{-1} k_*

所以q(f_* | x_*, X, y) = N(f_* | μ_*, σ^2_*)。

步骤2：通过积分得到预测概率

预测概率为：

p(y_*=1 | x_*, X, y) ≈ ∫ σ(f_*) q(f_* | x_*, X, y) df_*

这个积分没有解析解（因为σ是非线性的）。常用的近似方法是：

• 近似方法1：蒙特卡洛采样：从q(f_*)中采样f_*，计算σ(f_*)再平均。
• 近似方法2：解析近似：使用probit似然Φ（标准正态CDF）代替逻辑函数σ，这样积分有闭合形式。或者对σ进行局部近似。
• 最常用的方法：直接使用μ_*进行预测，即p(y_*=1) ≈ σ(μ_*)。虽然不够精确，但计算简单，实践中常被采用。

五、边缘似然的近似与超参数学习

模型的超参数（如核函数的长度尺度、方差等）通常通过最大化对数边缘似然 log p(y|X, θ) 来选择。在拉普拉斯近似下，我们可以得到这个量的近似：

log p(y|X) ≈ log q(y|X) = Ψ(\hat{f}) - (1/2) log |A| + (n/2) log(2π)

其中Ψ(\hat{f})是在众数处的对数联合分布值。这个近似允许我们使用梯度下降法来优化超参数。

六、总结与要点

1. 核心动机：高斯过程分类中，由于非高斯似然，后验分布p(f|X,y)没有解析形式。拉普拉斯近似用一个高斯分布来近似它。
2. 关键步骤：
   • 寻找众数：通过牛顿法优化，找到使对数联合分布Ψ(f)最大的\hat{f}。
   • 计算协方差：在\hat{f}处计算对数后验的负Hessian矩阵的逆，作为近似后验的协方差。
3. 预测：测试点的潜在函数值f_*的后验也是高斯的（均值和方差由训练数据决定），预测概率需要对σ(f_*)在这个高斯分布上积分，通常用近似方法计算。
4. 优势与局限：
   • 优势：相对简单，计算效率比MCMC采样高，是GPC中最常用的近似推断方法之一。
   • 局限：假设后验是单峰高斯的，当真实后验是多峰或严重非对称时，近似可能不准确。对于高度不确定性的分类问题，期望传播（EP）通常能提供更精确的近似。

拉普拉斯近似为高斯过程分类提供了一个坚实的计算框架，使得我们可以将强大的高斯过程非参数建模能力应用于分类问题。