高斯过程分类的期望传播（EP）算法原理与计算过程
我将为你讲解高斯过程分类（GPC）中一种重要的近似推断方法——期望传播（Expectation Propagation, EP）算法。这个过程涉及将高斯过程先验与分类似然结合，并通过消息传递来近似后验分布。
题目描述
高斯过程分类的核心问题：对于二分类任务，我们有一个高斯过程先验 \(f \sim \mathcal{GP}(0, k(\cdot, \cdot))\) 作用于隐函数 \(f\)，观测数据 \(D = \{(x_i, y_i)\}_{i=1}^n\)，其中 \(y_i \in \{-1, +1\}\)。分类似然为 \(p(y_i | f_i) = \Phi(y_i f_i)\)（对于probit模型）或 \(\sigma(y_i f_i)\)（对于logit模型），其中 \(f_i = f(x_i)\)，\(\Phi\) 是标准正态CDF，\(\sigma\) 是sigmoid函数。目标是计算隐函数值的后验分布 \(p(f | D)\) 和新输入的预测分布 \(p(y_* | x_*, D)\)。由于似然非高斯，后验无法精确计算，EP算法通过局部近似和迭代优化来逼近后验。
解题过程
步骤1：问题形式化与后验分解
后验分布为：
\[p(f | D) \propto p(f) \prod_{i=1}^n p(y_i | f_i) \]
其中 \(p(f) = \mathcal{N}(f | 0, K)\) 是先验，\(K\) 是核矩阵 \(K_{ij} = k(x_i, x_j)\)。每个数据点贡献一个非高斯因子 \(p(y_i | f_i)\)。EP的核心思想是用高斯分布 \(\tilde{t}_i(f_i) = \tilde{Z}_i \mathcal{N}(f_i | \tilde{\mu}_i, \tilde{\sigma}_i^2)\) 来近似每个似然项，其中 \(\tilde{Z}_i, \tilde{\mu}_i, \tilde{\sigma}_i\) 是待优化的局部参数。近似后验 \(q(f)\) 取为：
\[q(f) \propto p(f) \prod_{i=1}^n \tilde{t}_i(f_i) \]
由于所有因子都是高斯，\(q(f)\) 也是高斯分布。
步骤2：期望传播的迭代框架
EP通过迭代地优化每个近似因子 \(\tilde{t}_i\) 来工作，基本迭代步骤（对第 \(i\) 个数据点）如下：
- 移除近似因子：从当前近似后验 \(q(f)\) 中移除 \(\tilde{t}_i\) 的影响，得到“去除了第 \(i\) 个数据点”的分布，称为消去后验（cavity distribution）：
\[ q^{\setminus i}(f) \propto \frac{q(f)}{\tilde{t}_i(f_i)} \]
由于高斯分布的性质，这对应于从 \(q(f)\) 的精度矩阵中减去 \(\tilde{t}_i\) 的贡献。
- 计算真后验投影：结合真实的似然 \(p(y_i | f_i)\) 与消去后验在 \(f_i\) 上的边缘分布 \(q^{\setminus i}(f_i) = \mathcal{N}(f_i | \mu^{\setminus i}, \nu^{\setminus i})\)，形成“ tilted distribution”：
\[ \hat{p}(f_i) \propto q^{\setminus i}(f_i) p(y_i | f_i) \]
这是一个一维分布（非高斯），计算其均值和方差：
\[ \hat{\mu}_i = \mathbb{E}_{\hat{p}}[f_i], \quad \hat{\nu}_i = \mathbb{V}_{\hat{p}}[f_i] \]
- 更新近似因子：选择新的近似因子 \(\tilde{t}_i^{\text{new}}\) 使得将 \(q^{\setminus i}(f_i) \tilde{t}_i^{\text{new}}(f_i)\) 归一化后，其均值和方差等于 \(\hat{\mu}_i\) 和 \(\hat{\nu}_i\)。这通过匹配矩来实现：设 \(\tilde{t}_i^{\text{new}}(f_i) = \tilde{Z}_i^{\text{new}} \mathcal{N}(f_i | \tilde{\mu}_i^{\text{new}}, \tilde{\sigma}_i^{\text{new}2})\)，则：
\[ \frac{1}{Z} q^{\setminus i}(f_i) \tilde{t}_i^{\text{new}}(f_i) = \mathcal{N}(f_i | \hat{\mu}_i, \hat{\nu}_i) \]
其中 \(Z\) 是归一化常数。通过高斯分布的乘法公式，可解出：
\[ \tilde{\sigma}_i^{\text{new}2} = \left( \frac{1}{\hat{\nu}_i} - \frac{1}{\nu^{\setminus i}} \right)^{-1}, \quad \tilde{\mu}_i^{\text{new}} = \tilde{\sigma}_i^{\text{new}2} \left( \frac{\hat{\mu}_i}{\hat{\nu}_i} - \frac{\mu^{\setminus i}}{\nu^{\setminus i}} \right), \quad \tilde{Z}_i^{\text{new}} = Z \cdot \sqrt{2\pi \tilde{\sigma}_i^{\text{new}2}} \exp\left( \frac{(\tilde{\mu}_i^{\text{new}})^2}{2\tilde{\sigma}_i^{\text{new}2}} \right) \]
- 更新全局后验：用 \(\tilde{t}_i^{\text{new}}\) 替换旧的 \(\tilde{t}_i\)，重新计算全局近似后验 \(q(f)\)。
步骤3：一维矩计算（核心数值步骤）
对于probit似然 \(p(y_i | f_i) = \Phi(y_i f_i)\)，消去后验边缘 \(q^{\setminus i}(f_i) = \mathcal{N}(f_i | \mu^{\setminus i}, \nu^{\setminus i})\)， tilted distribution 为：
\[\hat{p}(f_i) = \frac{1}{\hat{Z}_i} \mathcal{N}(f_i | \mu^{\setminus i}, \nu^{\setminus i}) \Phi(y_i f_i) \]
其中 \(\hat{Z}_i = \int \mathcal{N}(f_i | \mu^{\setminus i}, \nu^{\setminus i}) \Phi(y_i f_i) df_i = \Phi\left( \frac{y_i \mu^{\setminus i}}{\sqrt{1 + \nu^{\setminus i}}} \right)\)。均值和方差可通过以下公式计算：
\[\hat{\mu}_i = \mu^{\setminus i} + \nu^{\setminus i} \cdot \frac{\mathcal{N}(z_i)}{ \Phi(z_i) \sqrt{1 + \nu^{\setminus i}} } \cdot y_i, \quad \hat{\nu}_i = \nu^{\setminus i} - (\nu^{\setminus i})^2 \cdot \frac{\mathcal{N}(z_i)}{ \Phi(z_i) (1 + \nu^{\setminus i}) } \cdot \left( z_i + \frac{\mathcal{N}(z_i)}{ \Phi(z_i) } \right) \]
其中 \(z_i = \frac{y_i \mu^{\setminus i}}{\sqrt{1 + \nu^{\setminus i}}}\)，\(\mathcal{N}\) 是标准正态PDF。对于logit似然，矩计算需用数值积分或近似（如基于二次近似的解析近似）。
步骤4：全局后验的参数化与更新
近似后验 \(q(f) = \mathcal{N}(f | \mu, \Sigma)\)。在EP中，通常参数化为：
\[\Sigma = (K^{-1} + \tilde{\Sigma}^{-1})^{-1}, \quad \mu = \Sigma \tilde{\Sigma}^{-1} \tilde{\mu} \]
其中 \(\tilde{\Sigma} = \text{diag}(\tilde{\sigma}_1^2, \dots, \tilde{\sigma}_n^2)\)，\(\tilde{\mu} = (\tilde{\mu}_1, \dots, \tilde{\mu}_n)^\top\)。迭代中更新 \(\tilde{\mu}_i, \tilde{\sigma}_i^2\) 后，需同步更新 \(\mu\) 和 \(\Sigma\)。为效率，常维护和更新 \(\Sigma\) 的Cholesky分解。
步骤5：预测分布计算
对于新输入 \(x_*\)，隐函数值的后验预测分布为：
\[q(f_*) = \int p(f_* | f) q(f) df = \mathcal{N}(f_* | k_*^\top K^{-1} \mu, k_{**} - k_*^\top (K^{-1} - K^{-1} \Sigma K^{-1}) k_*) \]
其中 \(k_* = [k(x_*, x_1), \dots, k(x_*, x_n)]^\top\)，\(k_{**} = k(x_*, x_*)\)。分类预测概率为：
\[p(y_* = +1 | x_*, D) \approx \int \Phi(f_*) q(f_*) df_* = \Phi\left( \frac{\mathbb{E}[f_*]}{\sqrt{1 + \mathbb{V}[f_*]}} \right) \]
步骤6：算法初始化与收敛
初始化所有近似因子 \(\tilde{t}_i\) 为 \(\tilde{\mu}_i = 0, \tilde{\sigma}_i^2 = \infty\)（即无信息），此时 \(q(f) = p(f)\)。然后迭代更新所有数据点（顺序或随机），直到 \(\tilde{\mu}_i, \tilde{\sigma}_i^2\) 变化小于阈值。EP通常收敛快，但可能不保证全局收敛，实践中常采用阻尼更新（如 \(\theta^{\text{new}} = \theta^{\text{old}} + \epsilon (\theta^{\text{target}} - \theta^{\text{old}})\)）来稳定。
总结
期望传播为高斯过程分类提供了一个高效且通常准确的确定性近似。它通过局部矩匹配迭代地优化近似后验，平衡了计算复杂度和精度，尤其适合中小规模数据集，是贝叶斯分类中广泛使用的方法之一。