机器学习数学推导（十四）：变分推断与变分 EM

后验 $p(\mathbf{z}\mid\mathbf{x})$ 无法直接计算时，我们面临两条路径。采样方法（MCMC）通过构造一条马尔可夫链，使其平稳分布恰好等于目标后验——理论上最终能精确逼近，但收敛缓慢且诊断困难。变分推断（VI）则另辟蹊径：先选定一个结构简单的分布族 $\mathcal{Q}$ ，再从中找出最接近真实后验的那个成员 $q^\star$ 。如此一来，推断问题就转化为优化问题——训练神经网络的那一套工具，现在也能用来拟合贝叶斯模型了。

本文将从一个核心恒等式出发，推导出平均场假设与坐标上升变分推断（CAVI）算法，阐明经典 EM 与变分 EM 如何作为特例自然浮现，并最终揭示重参数化技巧如何将 ELBO 转化为兼容自动微分的随机目标函数——这正是所有变分自编码器（VAE）背后的核心引擎。

你将学到什么#

• ELBO 恒等式如何将推断问题转化为优化问题
• 平均场假设及其导出的闭式坐标上升更新规则
• 变分 EM 如何通过优化 $q$ 来替代经典 EM 中的精确 E 步
• 重参数化技巧：一种低方差梯度估计器，支持对采样操作进行反向传播
• 何时应选用变分推断——以及它在哪些情况下会悄然误导你（模式寻求、方差低估）

前置知识#

• EM 算法与 ELBO（来自第十三篇 ）
• KL 散度与 Jensen 不等式
• 多元微积分与指数族分布
• 随机梯度估计的基本概念

后验瓶颈#

分子很简单，就是模型定义的联合密度；而分母——即证据 $p(\mathbf{x})$ ——需要在整个隐变量空间上对联合分布积分。除非模型属于共轭指数族，否则该积分通常无法解析求解。

学术界主要采用两种策略应对这一挑战：

变分推断以偏差为代价换取速度。本文后续将量化这一权衡，并展示如何在实践中加以管理。

ELBO 恒等式#

由于 $\log p(\mathbf{x})$ 与 $q$ 无关，最大化 ELBO 等价于最小化 $q$ 到真实后验的 KL 散度。

图 1： 一旦数据和模型确定，$\log p(\mathbf{x})$ 就固定不变。ELBO 每提升一点，KL 间隙就相应缩小同等幅度；完美推断时，间隙完全闭合。又因 KL 散度非负，ELBO 自然构成边缘似然的一个严格下界——这使其成为模型比较的有力工具。

值得驻足思考的是：我们用一个优化问题（在函数空间中寻找最优 $q$ ）替换了原本无法计算的积分（证据）。变分推断的所有工程实践，本质上都是围绕这一关键转换展开的。

平均场近似#

每个因子 $q_j$ 可属于不同的分布族，无需额外参数化约束。

图 2： 中间图虽匹配了真实后验（左图）的边缘分布，却丢失了所有变量间的相关性。右图展示了 $\mathrm{KL}(q\|p)$ 的全局最优解：方差按 $1-\rho^2$ 的比例收缩，导致近似分布在两个维度上都显著低估了不确定性。这种系统性的方差压缩，正是平均场 VI 最常见的失效模式。

最优因子的形式#

这是平均场 VI 的核心公式：第 $j$ 个坐标的最优变分因子，等于联合分布在其他所有因子下的几何平均（再归一化）。

坐标上升变分推断（CAVI）#

由于每个因子的最优解依赖于其他因子，我们采用循环更新策略：

1
2
3
4
5
6

初始化 q_1, ..., q_M
重复
    for j = 1, ..., M:
        log q_j(z_j) <- E_{q_{-j}}[log p(x, z)] + 常数
        归一化 q_j
直到 ELBO 变化 < 阈值

每次更新都是对一个关于单个坐标的凹函数进行坐标上升，因此 ELBO 单调不减。算法必定收敛至某个局部最优解，但无法保证达到全局最优。

图 5： 对一个相关二维高斯分布执行 8 轮 CAVI。对角化的 $q$ 不仅均值被迅速拉向原点，方差也随之塌缩至 $1/\text{precision}_{ii}$ 。右图展示了 ELBO 的单调上升轨迹——这是任何 VI 实现中极为实用的收敛诊断指标。

共轭指数族#

当模型属于共轭指数族时——即每个条件分布 $p(z_j\mid \mathbf{z}_{-j},\mathbf{x})$ 均属于某指数族——CAVI 更新存在闭式解。此时最优 $q_j$ 与条件分布同属一个指数族，更新操作等价于在 $q_{-j}$ 下对自然参数取平均。这类模型包括贝叶斯高斯混合模型（GMM）、LDA、贝叶斯线性回归、带 Dirichlet 先验的隐马尔可夫模型（HMM）等。对于非共轭模型，则需借助第 6 节 所述的黑盒方法。

把变分族当作逼近器#

若抛开平均场假设的束缚，我们会对 VI 有更深刻的理解。不妨任选一个参数化分布族 $q_\phi(\mathbf{z})$ ——例如带可学习均值与协方差的高斯分布、归一化流，或摊还推断网络——然后直接对参数 $\phi$ 最小化 $\mathrm{KL}(q_\phi\,\|\,p(\cdot\mid\mathbf{x}))$ 。

图 3： 面对双峰目标分布，反向 KL（紫色）的最优解会锁定单一峰值——因其具有“模式寻求”特性：只要 $q>0$ 而 $p\approx 0$ ，$\mathrm{KL}(q\|p)$ 就会施加无穷大的惩罚。相比之下，正向 KL（橙色，期望传播 EP 所用）则倾向于匹配矩并“覆盖”所有模式，即便这意味着在低概率区域分配密度也在所不惜。VI 之所以采用反向 KL，是因为仅需从 $q$ 中采样即可计算；但这种不对称性确实会对下游的不确定性估计产生实质性影响。

图 6： 在对称混合分布上，反向 KL 存在多个局部最优，各自对应一个峰值；而正向 KL 则给出唯一、宽泛且居中的解。设计 VI 系统时，明确下游任务能容忍哪种行为，往往能省去大量调试功夫。

平均场的隐患#

图 2 中的方差低估并非实现缺陷，而是反向 KL 散度与因子分解假设共同作用下的必然结果。若下游任务依赖后验方差——如贝叶斯模型平均、校准预测或决策理论——平均场 VI 必然会系统性地低估不确定性。以下是三种常见应对策略：

1. 结构化变分推断：保留部分变量间的依赖关系（例如采用树结构的 $q$ ）；
2. 归一化流：通过一系列可逆变换赋予 $q$ 足够灵活性以建模相关性；
3. 摊还推断：使用灵活的神经网络作为编码器（如 VAE 中的做法）。

变分 EM#

第十三篇 指出，经典 EM 算法本身也建立在 ELBO 恒等式之上。EM 交替执行以下两步：

• E 步：令 $q(\mathbf{z}) = p(\mathbf{z}\mid\mathbf{x};\boldsymbol{\theta}^{(t)})$ ，即精确后验——此时 KL 间隙为零，ELBO 恰好等于 $\log p(\mathbf{x};\boldsymbol{\theta}^{(t)})$ ；
• M 步：固定 $q$ ，对 $\boldsymbol{\theta}$ 最大化 ELBO，这等价于最大化 $\mathbb{E}_q[\log p(\mathbf{x},\mathbf{z};\boldsymbol{\theta})]$ 。

当精确后验不可得时，E 步便失效了。变分 EM 的解决方案很简单：用一个 VI 子程序替代 E 步：

在变分 EM 中，E 步完成后 ELBO 不再紧贴对数似然（KL 间隙非零），因此算法实际优化的是 $\log p(\mathbf{x};\boldsymbol{\theta})$ 的一个下界，而非似然本身。不过，该下界仍会单调递增——这在实践中是一个非常有用的保证。

变分贝叶斯 EM（VBEM）更进一步，将 $\boldsymbol{\theta}$ 也视为随机变量，并为其引入独立的变分因子，从而得到 $\mathbf{z}$ 与 $\boldsymbol{\theta}$ 的联合后验近似。变分贝叶斯 GMM 与变分 LDA 正是基于此框架求解的。

黑盒变分推断与重参数化技巧#

此处的期望依赖于 $q_\phi$ ，而 $q_\phi$ 本身又由 $\phi$ 决定，因此梯度不能直接移入期望内部。

该估计器无偏，但方差极高，通常需配合控制变量及大样本量才能实用。

以对角高斯分布 $q_\phi(\mathbf{z}) = \mathcal{N}(\boldsymbol{\mu}_\phi,\,\mathrm{diag}(\boldsymbol{\sigma}_\phi^2))$ 为例，其变换为 $g_\phi(\boldsymbol{\epsilon}) = \boldsymbol{\mu}_\phi + \boldsymbol{\sigma}_\phi\odot\boldsymbol{\epsilon}$ ，其中 $\boldsymbol{\epsilon}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$ 。整个计算图完全可微，自动微分可处理其余部分，所得梯度方差远低于 REINFORCE。这正是变分自编码器（VAE）及几乎所有现代连续隐变量变分模型的核心驱动力。

对于离散隐变量，该技巧失效（无法将 $z\in\{0,1\}$ 表示为 $\epsilon$ 的光滑函数）。常用替代方案包括：带控制变量的 REINFORCE、Gumbel-Softmax / Concrete 松弛，或先做连续松弛再采用直通估计（straight-through estimation）。

如何选择#

图 4： 典型的速度-精度权衡图。VI 能快速达到较小误差，随后便停滞在其偏差下限（即 $\mathcal{Q}$ 与真实后验之间的差距）；MCMC 则需经历 burn-in 阶段，之后以经典的 $1/\sqrt{T}$ 速率逐步逼近零误差。决策准则十分简洁：

• 选用 VI：适用于大数据集、在线学习、探索性建模，或作为神经网络的内部组件；
• 选用 MCMC：当后验分布本身至关重要时——如药物剂量响应分析、科学推断，或任何将公开发表的结果；
• 两者结合：用 VI 的均值为 HMC 链提供暖启动，兼顾快速初始化与渐近无偏性。

应用：变分 LDA#

其中 $\theta_d$ 与 $\beta_k$ 采用狄利克雷变分因子，而每个词的主题分配 $z_{d,n}$ 则使用分类（categorical）因子。

图 7： 变分 LDA 的两大输出。左：八篇文档的主题比例后验均值 $\mathbb{E}_q[\theta_d]$ 。每篇文档均为四个主题的软混合——文档 1 和 5 主要属于“机器学习”，文档 3 则主导“金融”，依此类推。右：各主题的词分布后验均值 $\mathbb{E}_q[\beta_k]$ 。每个主题的概率质量集中于少数特征词汇——正是这种可解释性，使 LDA 成为 2000 至 2010 年代文本分析的主力工具。

随机变分推断（SVI，Hoffman 等 2013）通过小批量采样与自然梯度更新 $\lambda_k$ ，将相同更新机制扩展至十亿级文档规模，完美体现了图 4 所示的速度优势。

实现：变分贝叶斯 GMM#

这是一个简洁的共轭贝叶斯 GMM 的 CAVI 实现。更新公式的数学细节可参考 Bishop PRML §10.2；此处我们更关注循环结构的清晰表达。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71

import numpy as np
from scipy.special import digamma

class VariationalGMM:
    """平均场变分贝叶斯高斯混合（CAVI 版）。"""

    def __init__(self, K=3, max_iter=100, tol=1e-3):
        self.K = K
        self.max_iter = max_iter
        self.tol = tol

    def fit(self, X):
        N, d = X.shape
        # 弱共轭先验
        alpha0, beta0, nu0 = 1.0, 1.0, float(d)
        m0, W0 = X.mean(0), np.eye(d)

        # 初始化变分参数为先验加上均匀分布的质量
        self.alpha = np.full(self.K, alpha0 + N / self.K)
        self.beta = np.full(self.K, beta0 + N / self.K)
        self.nu = np.full(self.K, nu0 + N / self.K)
        self.m = np.array([m0 + 0.1 * np.random.randn(d) for _ in range(self.K)])
        self.W = np.array([W0.copy() for _ in range(self.K)])

        r = np.random.dirichlet([1] * self.K, N)
        for _ in range(self.max_iter):
            r_old = r.copy()

            # ---- E 步：计算变分责任度 q(z_n) ----
            r = self._update_r(X, N, d)

            # ---- M 步：更新 Dirichlet 和 Normal-Wishart 因子 ----
            N_k = r.sum(0)
            x_bar = (r.T @ X) / N_k[:, None]
            self.alpha = alpha0 + N_k
            self.beta = beta0 + N_k
            self.nu = nu0 + N_k
            self.m = (beta0 * m0 + N_k[:, None] * x_bar) / self.beta[:, None]
            for k in range(self.K):
                diff = X - x_bar[k]
                S = (r[:, k, None] * diff).T @ diff / N_k[k]
                dm = x_bar[k] - m0
                self.W[k] = np.linalg.inv(
                    np.linalg.inv(W0)
                    + N_k[k] * S
                    + beta0 * N_k[k] / (beta0 + N_k[k]) * np.outer(dm, dm)
                )

            if np.max(np.abs(r - r_old)) < self.tol:
                break
        return self

    def _update_r(self, X, N, d):
        """计算 rho_{nk} = exp( E_q[log pi_k] + 0.5 E_q[log|Lambda_k|]
                                 - 0.5 E_q[(x-mu_k)^T Lambda_k (x-mu_k)] )"""
        r = np.zeros((N, self.K))
        for k in range(self.K):
            E_lp = digamma(self.alpha[k]) - digamma(self.alpha.sum())
            E_ld = sum(digamma((self.nu[k] + 1 - i) / 2) for i in range(1, d + 1))
            E_ld += d * np.log(2) + np.log(max(np.linalg.det(self.W[k]), 1e-10))
            for n in range(N):
                diff = X[n] - self.m[k]
                r[n, k] = (E_lp
                           + 0.5 * E_ld
                           - 0.5 * (self.nu[k] * diff @ self.W[k] @ diff
                                    + d / self.beta[k]))
        r = np.exp(r - r.max(1, keepdims=True))
        return r / r.sum(1, keepdims=True)

    def predict(self, X):
        return np.argmax(self._update_r(X, len(X), X.shape[1]), axis=1)

该算法的整体流程与第十三篇的 EM-GMM 相似——交替更新责任度与成分统计量——但所有量均为 $q$ 下的后验矩，而非点估计。一个额外优势是自动成分剪枝：空成分的 $\alpha_k$ 会逐渐收缩至先验值 $\alpha_0$ ，使得无用主题悄然淡出，避免拟合噪声。

常见问题#

Q1：为何使用反向 KL 而非正向 KL？
反向 $\mathrm{KL}(q\|p)$ 仅需在可控的 $q$ 下计算期望；正向 $\mathrm{KL}(p\|q)$ 则需在难以处理的真实后验 $p$ 下求期望。代价是反向 KL 会表现出模式寻求行为——参见图 6。

Q2：我的 ELBO 在下降，是否出错了？
是的。CAVI 的 ELBO 按构造应单调非减。若出现下降，几乎可以肯定是代码存在 bug，常见原因包括：(i) 责任度归一化错误；(ii) 熵项符号弄反；(iii) 使用了尚未更新的变分参数。

Q3：ELBO 下界有多紧？
$\log p(\mathbf{x}) - \mathcal{L}(q^\star) = \mathrm{KL}(q^\star\|p(\cdot\mid\mathbf{x})) \geq 0$ 。该间隙精确反映了变分族的建模误差。因此，仅当不同模型的间隙量级相近时，基于 ELBO 的模型比较才是公平的。

Q4：平均场方法何时会彻底失效？
当后验存在强相关性（图 2）或多峰结构（图 6）时。典型症状包括：后验方差被严重低估、预测过度自信、对初始化高度敏感。解决方法包括结构化 VI、归一化流或全协方差摊还推断。

Q5：重参数化与 REINFORCE，如何选择？
对连续可微的 $q$ （如高斯、Logistic，或经 Gumbel-Softmax 处理的混合分布），优先使用重参数化；对离散隐变量，则使用 REINFORCE。当两者皆适用时，重参数化的梯度方差通常比 REINFORCE 低 100–1000 倍。

Q6：变分 EM 与完整变分贝叶斯（VBEM），何时选用后者？
若需获取 $\boldsymbol{\theta}$ 的不确定性（如小数据场景、模型选择或决策理论），应选用 VBEM；若仅需良好点估计（如大数据、追求预测精度），变分 EM 更快且更简单。

练习题#

E1. 从零开始证明 $\mathcal{L}(q) \leq \log p(\mathbf{x})$ 。
提示： 由 $\log p(\mathbf{x}) = \mathcal{L}(q) + \mathrm{KL}(q\|p(\cdot\mid\mathbf{x}))$ 及 Jensen 不等式可知 KL 散度非负。

E2. 对二变量平均场 $q(z_1,z_2) = q_1(z_1)q_2(z_2)$ ，写出最优 $q_1^\star$ 。
答案： $\log q_1^\star(z_1) = \mathbb{E}_{q_2}[\log p(\mathbf{x},z_1,z_2)] + \text{常数}$ 。

E3. 推导零均值二维高斯分布（精度矩阵为 $\Lambda$ ）的 CAVI 闭式更新，并证明收敛时 $\mathrm{Var}_q[z_j] = 1/\Lambda_{jj}$ ——即精度矩阵对角元的倒数，而非协方差矩阵对角元。
提示： 参考图 5 的轨迹进行推导。

E4. 证明变分 EM 的 ELBO 在迭代中单调不减，即使 E 步不再精确。
提示： E 步通过优化 $q$ 提升 ELBO；M 步通过优化 $\boldsymbol{\theta}$ 提升 ELBO；交替结构保证了整体单调性。

E5. 为何直接用蒙特卡洛估计 $\nabla_\phi \mathbb{E}_{q_\phi}[f(\mathbf{z})] \approx \frac{1}{S}\sum_s \nabla_\phi f(\mathbf{z}^{(s)})$ （$\mathbf{z}^{(s)}\sim q_\phi$ ）会出错？
答案： 样本 $\mathbf{z}^{(s)}$ 本身依赖于 $\phi$ ，而朴素梯度忽略了这一依赖。REINFORCE 与重参数化是两种正确的修正方法。

下一步#

EM 和变分 EM 处理的是"独立同分布数据 + 隐变量"。但现实里有一类数据天然带有时间顺序——语音、文本、动作序列、生物序列——同一时刻的观测依赖于前一时刻的隐状态。把变分 EM 直接套上去会丢掉时间结构。下一章是这类问题的经典模型：隐马尔可夫模型（HMM）。

HMM 把隐变量串成一条马尔可夫链：当前隐状态只依赖前一隐状态，当前观测只依赖当前隐状态。这个看似很硬的假设，给了我三个超级好用的精确推断算法——前向算法（边缘似然）、Viterbi（最优路径）、Baum-Welch（即 EM 在 HMM 上的特化）。我会把它们当作"动态规划在概率图模型上的具体实现"来推，而不是当作三个独立公式来背。理解了 HMM 的精确推断，后面 CRF 的判别式版本、RNN/Transformer 的连续版本、状态空间模型的现代复兴，都会变得容易看穿。