高斯过程分类（Gaussian Process Classification, GPC）的期望传播（Expectation Propagation, EP）算法原理与计算过程
题目描述
在高斯过程分类（GPC）中，我们通常使用高斯过程（GP）先验对潜在函数进行建模，并通过一个非线性激活函数（如逻辑函数）将潜在函数的值映射到类别概率。然而，由于这个非线性映射，后验分布不再是高斯的，使得精确推断变得难以处理。期望传播（Expectation Propagation, EP）算法是一种近似推断方法，它通过迭代地将每个非高斯似然项用一个非归一化的高斯分布（即一个高斯“因子”）来近似，从而得到整个后验的一个高斯近似。本题目要求详细解释在二分类GPC中，EP算法的原理、推导步骤和迭代计算过程。
解题过程
1. 高斯过程分类模型设定
我们考虑一个二分类问题，训练数据集为 \(D = \{(\mathbf{x}_i, y_i)\}_{i=1}^n\)，其中 \(y_i \in \{-1, +1\}\)。我们引入一个隐函数 \(f: \mathcal{X} \to \mathbb{R}\)，并假设其服从一个零均值的高斯过程先验：
\[f \sim \mathcal{GP}(0, k(\cdot, \cdot)) \]
其中 \(k\) 是一个协方差函数（核函数），例如径向基函数（RBF）核。
对于分类，我们通过一个sigmoid函数（如逻辑函数）将隐函数值 \(f_i = f(\mathbf{x}_i)\) 映射到类别概率。常用的选择是累积高斯函数（probit函数）或逻辑函数。这里以累积高斯函数为例：
\[p(y_i | f_i) = \Phi(y_i f_i) = \int_{-\infty}^{y_i f_i} \mathcal{N}(z | 0, 1) dz \]
其中 \(\Phi\) 是标准正态分布的累积分布函数。这个似然函数是非高斯的。
整个模型的联合分布为：
\[p(\mathbf{y}, \mathbf{f} | X) = p(\mathbf{f} | X) \prod_{i=1}^n p(y_i | f_i) \]
其中 \(\mathbf{f} = [f_1, \dots, f_n]^T\)，先验 \(p(\mathbf{f} | X) = \mathcal{N}(\mathbf{f} | \mathbf{0}, K)\)，\(K\) 是核矩阵，\(K_{ij} = k(\mathbf{x}_i, \mathbf{x}_j)\)。
我们的目标是计算后验分布 \(p(\mathbf{f} | D) \propto p(\mathbf{f} | X) \prod_i p(y_i | f_i)\)，以及在新点 \(\mathbf{x}_*\) 上的预测概率 \(p(y_* = +1 | D, \mathbf{x}_*)\)。
2. 精确后验的困难与非高斯性
由于每个似然项 \(p(y_i | f_i)\) 是非高斯的（它不是一个高斯分布），而先验是高斯的，整个后验 \(p(\mathbf{f} | D)\) 不是高斯分布。这使得计算后验的均值和协方差、以及预测分布变得解析不可行（无法写成闭式解）。我们需要一种近似推断方法。
3. 期望传播（EP）的核心思想
EP算法的核心思想是：用一个更简单的分布 \(q(\mathbf{f})\) 来近似真实后验 \(p(\mathbf{f} | D)\)。具体地，我们用一个高斯分布来近似，因为高斯分布具有良好的解析性质。
我们注意到，后验可以写成：
\[p(\mathbf{f} | D) \propto p(\mathbf{f} | X) \prod_{i=1}^n p(y_i | f_i) = \mathcal{N}(\mathbf{f} | \mathbf{0}, K) \prod_{i=1}^n t_i(f_i) \]
其中 \(t_i(f_i) = p(y_i | f_i)\) 是每个数据点的似然项，它是一个关于 \(f_i\) 的（非高斯）因子。
EP算法用一个非归一化的高斯分布（称为“局部近似因子”） \(\tilde{t}_i(f_i) = \tilde{Z}_i \cdot \mathcal{N}(f_i | \tilde{\mu}_i, \tilde{\sigma}_i^2)\) 来近似每个真实似然因子 \(t_i(f_i)\)。注意，这里的 \(\tilde{Z}_i\) 是一个标量，用于调整归一化常数；而 \(\tilde{\mu}_i, \tilde{\sigma}_i^2\) 是这个局部高斯因子的均值和方差参数。
那么，近似后验 \(q(\mathbf{f})\) 就是先验与所有这些局部高斯因子的乘积：
\[q(\mathbf{f}) \propto \mathcal{N}(\mathbf{f} | \mathbf{0}, K) \prod_{i=1}^n \tilde{t}_i(f_i) \]
由于每个 \(\tilde{t}_i\) 是 \(f_i\) 的高斯函数，且先验是高斯的，所以整个 \(q(\mathbf{f})\) 也是一个高斯分布（多元高斯）。我们可以解析地计算其均值和协方差。
4. EP算法的迭代更新步骤
EP算法通过迭代地优化每个局部近似因子 \(\tilde{t}_i\) 来工作。算法的每一步针对一个数据点 \(i\) 执行以下操作：
a. 构造“去除近似”的分布（cavity distribution）
首先，从当前的近似后验 \(q(\mathbf{f})\) 中“移除”第 \(i\) 个局部近似因子 \(\tilde{t}_i\) 的影响，得到一个称为 cavity distribution 的分布：
\[q^{\setminus i}(\mathbf{f}) \propto \frac{q(\mathbf{f})}{\tilde{t}_i(f_i)} \]
由于 \(q(\mathbf{f})\) 是高斯分布，且 \(\tilde{t}_i\) 是只关于 \(f_i\) 的高斯函数，所以 \(q^{\setminus i}(\mathbf{f})\) 也是一个高斯分布。具体地，我们可以通过计算高斯条件分布的性质得到 \(q^{\setminus i}\) 的边际分布 \(q^{\setminus i}(f_i)\)，它是一个一维高斯分布，记其均值和方差为 \(\mu_{-i}\) 和 \(\sigma_{-i}^2\)。
b. 构造“倾斜分布”（tilted distribution）
将 cavity distribution 与真实的似然因子 \(t_i\) 结合，得到一个非高斯的“倾斜分布”：
\[\hat{p}_i(f_i) \propto q^{\setminus i}(f_i) \cdot t_i(f_i) \]
这个分布是精确的（在给定其他因子固定时），但它不再是高斯的。我们的目标是让近似后验 \(q(f_i)\) 尽可能接近这个倾斜分布。
c. 投影（矩匹配）
我们希望找到一个高斯分布 \(q^{\text{new}}(f_i)\) 来近似 \(\hat{p}_i(f_i)\)，使得两者的前两阶矩（均值和方差）匹配。这称为矩匹配（moment matching），是一种信息投影（将复杂分布投影到指数族分布上）。
计算 \(\hat{p}_i(f_i)\) 的均值和方差：
\[\mathbb{E}_{\hat{p}_i}[f_i] = \frac{\int f_i \cdot q^{\setminus i}(f_i) t_i(f_i) df_i}{\int q^{\setminus i}(f_i) t_i(f_i) df_i} \]
\[ \mathbb{V}_{\hat{p}_i}[f_i] = \frac{\int (f_i - \mathbb{E}_{\hat{p}_i}[f_i])^2 \cdot q^{\setminus i}(f_i) t_i(f_i) df_i}{\int q^{\setminus i}(f_i) t_i(f_i) df_i} \]
对于 probit 似然 \(t_i(f_i) = \Phi(y_i f_i)\)，这两个积分可以解析计算（它们可以表示为标准正态分布的累积分布函数和概率密度函数的组合）。具体公式如下：
设 \(v = \sigma_{-i}^2\)，\(m = \mu_{-i}\)，定义：
\[z_i = y_i \cdot \frac{m}{\sqrt{1 + v}} \]
（这里除以 \(\sqrt{1+v}\) 是因为 probit 函数中的参数是标准正态的累积分布函数，而 cavity 分布是 \(\mathcal{N}(m, v)\)，结合 probit 后需要归一化。）
然后，计算：
\[\alpha = \frac{\phi(z_i)}{\Phi(z_i) \sqrt{1+v}} \]
其中 \(\phi\) 是标准正态概率密度函数。那么，倾斜分布的均值和方差为：
\[\mathbb{E}_{\hat{p}_i}[f_i] = m + v \cdot y_i \cdot \alpha \]
\[ \mathbb{V}_{\hat{p}_i}[f_i] = v - v^2 \cdot \alpha \cdot (\alpha + z_i / \sqrt{1+v}) \]
（注：这是一个标准结果，推导涉及对 probit 模型的高斯积分。）
d. 更新局部近似因子
现在，我们有了新的高斯边际分布 \(q^{\text{new}}(f_i) = \mathcal{N}(f_i | \mu_i^{\text{new}}, \sigma_i^{2, \text{new}})\)，其中 \(\mu_i^{\text{new}} = \mathbb{E}_{\hat{p}_i}[f_i]\)，\(\sigma_i^{2, \text{new}} = \mathbb{V}_{\hat{p}_i}[f_i]\)。
我们需要反过来更新局部近似因子 \(\tilde{t}_i\)，使得当我们将 cavity distribution \(q^{\setminus i}\) 与新的 \(\tilde{t}_i^{\text{new}}\) 相乘时，能得到 \(q^{\text{new}}\)。即：
\[q^{\text{new}}(f_i) \propto q^{\setminus i}(f_i) \cdot \tilde{t}_i^{\text{new}}(f_i) \]
由于 \(q^{\setminus i}(f_i) = \mathcal{N}(f_i | \mu_{-i}, \sigma_{-i}^2)\)，而 \(\tilde{t}_i^{\text{new}}(f_i) = \tilde{Z}_i^{\text{new}} \cdot \mathcal{N}(f_i | \tilde{\mu}_i^{\text{new}}, \tilde{\sigma}_i^{2, \text{new}})\)，根据高斯分布的乘积公式，我们可以解出新的局部近似因子参数：
首先，由高斯乘积公式，有：
\[\frac{1}{\sigma_i^{2, \text{new}}} = \frac{1}{\sigma_{-i}^2} + \frac{1}{\tilde{\sigma}_i^{2, \text{new}}} \]
\[ \frac{\mu_i^{\text{new}}}{\sigma_i^{2, \text{new}}} = \frac{\mu_{-i}}{\sigma_{-i}^2} + \frac{\tilde{\mu}_i^{\text{new}}}{\tilde{\sigma}_i^{2, \text{new}}} \]
从中可以解出：
\[\tilde{\sigma}_i^{2, \text{new}} = \left( \frac{1}{\sigma_i^{2, \text{new}}} - \frac{1}{\sigma_{-i}^2} \right)^{-1} \]
\[ \tilde{\mu}_i^{\text{new}} = \tilde{\sigma}_i^{2, \text{new}} \cdot \left( \frac{\mu_i^{\text{new}}}{\sigma_i^{2, \text{new}}} - \frac{\mu_{-i}}{\sigma_{-i}^2} \right) \]
而新的归一化常数 \(\tilde{Z}_i^{\text{new}}\) 可以通过匹配 \(q^{\text{new}}(f_i)\) 与 \(q^{\setminus i}(f_i) \tilde{t}_i^{\text{new}}(f_i)\) 的归一化常数得到，但实际计算中通常我们只需要更新 \(q(\mathbf{f})\) 的全局参数，而不需要显式计算 \(\tilde{Z}_i\)。
e. 更新全局近似后验 \(q(\mathbf{f})\)
当我们更新了第 \(i\) 个局部近似因子后，需要更新整个高斯后验 \(q(\mathbf{f})\)。回忆 \(q(\mathbf{f}) \propto \mathcal{N}(\mathbf{f} | \mathbf{0}, K) \prod_j \tilde{t}_j(f_j)\)。我们可以将所有的局部近似因子写成一个统一的表达：设 \(\tilde{\boldsymbol{\mu}} = [\tilde{\mu}_1, \dots, \tilde{\mu}_n]^T\)，\(\tilde{\Sigma} = \operatorname{diag}(\tilde{\sigma}_1^2, \dots, \tilde{\sigma}_n^2)\)。那么，这些因子的乘积正比于 \(\mathcal{N}(\mathbf{f} | \tilde{\boldsymbol{\mu}}, \tilde{\Sigma})\)（注意这里的“正比”意味着可能忽略归一化常数，但形状是高斯）。
那么，\(q(\mathbf{f})\) 是一个高斯分布，其均值和协方差可以通过高斯分布相乘的公式得到：
\[q(\mathbf{f}) = \mathcal{N}(\mathbf{f} | \boldsymbol{\mu}, \Sigma) \]
其中
\[\Sigma = (K^{-1} + \tilde{\Sigma}^{-1})^{-1} \]
\[ \boldsymbol{\mu} = \Sigma \cdot \tilde{\Sigma}^{-1} \tilde{\boldsymbol{\mu}} \]
或者等价地，利用 Woodbury 矩阵恒等式，可以写成：
\[\Sigma = K - K (K + \tilde{\Sigma})^{-1} K \]
\[ \boldsymbol{\mu} = \Sigma \tilde{\Sigma}^{-1} \tilde{\boldsymbol{\mu}} = K (K + \tilde{\Sigma})^{-1} \tilde{\boldsymbol{\mu}} \]
这样，在每次更新一个 \(\tilde{t}_i\) 后，我们可以增量更新 \(\Sigma\) 和 \(\boldsymbol{\mu}\)，而不必重新计算整个矩阵求逆（有高效的 rank-1 更新方法）。
5. EP算法的整体流程
- 初始化：初始化所有局部近似因子 \(\tilde{t}_i\)，通常设为 \(\tilde{\mu}_i = 0\)，\(\tilde{\sigma}_i^2 = \infty\)（或一个很大的数），这等价于初始近似后验等于先验 \(q(\mathbf{f}) = \mathcal{N}(\mathbf{0}, K)\)。
- 迭代：对 \(i = 1, \dots, n\) 进行循环（通常多次扫描整个数据集），执行步骤 (a)-(e)：
a. 计算 cavity 参数 \(\mu_{-i}, \sigma_{-i}^2\) 从当前的 \(q(\mathbf{f})\)。
b. 计算倾斜分布 \(\hat{p}_i(f_i)\) 的矩 \(\mathbb{E}_{\hat{p}_i}[f_i]\)，\(\mathbb{V}_{\hat{p}_i}[f_i]\)。
c. 更新局部近似因子参数 \(\tilde{\mu}_i, \tilde{\sigma}_i^2\)。
d. 更新全局后验 \(q(\mathbf{f})\) 的 \(\boldsymbol{\mu}\) 和 \(\Sigma\)。 - 收敛：当局部近似因子的变化小于某个阈值，或达到最大迭代次数时停止。
- 预测：对于新输入 \(\mathbf{x}_*\)，计算隐函数 \(f_*\) 的后验预测分布 \(q(f_*) = \mathcal{N}(f_* | \mu_*, \sigma_*^2)\)，其中
\[ \mu_* = \mathbf{k}_*^T (K + \tilde{\Sigma})^{-1} \tilde{\boldsymbol{\mu}} \]
\[ \sigma_*^2 = k(\mathbf{x}_*, \mathbf{x}_*) - \mathbf{k}_*^T (K + \tilde{\Sigma})^{-1} \mathbf{k}_* \]
其中 \(\mathbf{k}_* = [k(\mathbf{x}_*, \mathbf{x}_1), \dots, k(\mathbf{x}_*, \mathbf{x}_n)]^T\)。然后，类别概率为：
\[ p(y_* = +1 | D, \mathbf{x}_*) = \int \Phi(f_*) q(f_*) df_* = \Phi\left( \frac{\mu_*}{\sqrt{1 + \sigma_*^2}} \right) \]
这里积分再次利用了 probit 函数与高斯卷积的性质。
总结
期望传播（EP）算法通过用局部高斯因子近似每个非高斯似然项，将后验近似为一个高斯分布。它通过迭代的矩匹配（投影）来优化这些局部近似，使得整体近似后验尽可能接近真实后验。在 Gaussian Process Classification 中，EP 提供了一个高效且通常很准确的后验高斯近似，使得预测可以解析计算。与拉普拉斯近似（Laplace Approximation）相比，EP 通常更精确，但计算量也稍大。