简介：这份PDF文档是概率图模型领域的经典综述，由Wainwright与Michael Jordan合著，面向机器学习、统计学方向的研究生与科研人员，帮助读者系统理解图模型、指数族与变分推断三者的内在联系。文档共1个PDF文件，压缩包约2.06MB，内容源自《Foundations and Trends in Machine Learning》2008年刊载的长文，篇幅达三百余页，结构完整、推导严谨。文中以指数族表示与累积量函数、熵之间的共轭对偶为核心工具，统一刻画似然、边际概率与最大后验配置的计算问题，并系统梳理了和积算法、簇变分方法、期望传播、平均场、最大积以及线性规划松弛等算法的变分解释。读者可借此掌握从精确推断到近似推断的完整方法论，理解变分方法与MCMC的互补关系，并了解其在生物信息学、通信理论、信号与图像处理等场景中的应用。目前已有235人学习，适合作为推断算法方向的核心参考读物。

1. 从一份 Graphical Models 讲义说起：指数族与变分推断到底卡在哪

如果你手头正对着一份名为 Graphical Models, Exponential Families and Variational Inference 的讲义或笔记，大概率已经翻过前几页：概率图模型讲得挺清楚，指数族公式也推得动，可一旦跳到变分推断那一章，符号突然密集起来，ELBO、平均场、坐标上升全堆在一起，读着读着就断了线。这不是你基础差，而是这三块内容天然咬合得很紧——图模型负责把联合分布拆成因子，指数族负责给每个因子一个统一的可微形式，变分推断则负责在这个形式下把难解的后验近似出来。任何一环没接上，后面就是黑匣子。

这份材料真正想解决的问题是：如何用指数族把图模型里的消息传递和变分更新写成统一公式，从而让推断过程可推导、可编程、可调试。它适合已经学过概率论和基础贝叶斯、但一碰到 LDA、隐马尔可夫或贝叶斯网络推断就卡壳的工程师和研究生。读完并动手复现之后，你应该能自己写出一个平均场变分推断的最小实现，而不是只会调库。下面我按“先立住理论、再动手复现、最后排坑”的顺序，把这条路走一遍。

2. 指数族与图模型的接合点：为什么变分推断非用它不可

2.1 指数族的标准形式与充分统计量

指数族不是某个具体分布，而是一大类分布的通用写法。把概率密度写成

p(x | η) = h(x) exp( ηᵀ T(x) − A(η) )

之后，高斯、伯努利、泊松、狄利克雷、伽马全都能塞进同一个框架。这里 η 叫自然参数，T(x) 是充分统计量，A(η) 是对数配分函数，也叫累积量函数。A(η) 的梯度恰好等于 T(x) 的期望，这个性质在推导变分更新时反复用到，是整份材料里最值得先记牢的一条。

为什么图模型非要接指数族？因为图模型把联合分布写成因子的乘积，如果每个因子都是指数族，那么整个联合分布的对数就变成若干线性项的求和。求导、取期望、做坐标上升时，形式统一，不会每换一个分布就重推一遍公式。常见做法是先把每个条件分布写成指数族，再套用共轭关系，这样后验更新就有闭式解可用。

提示：如果你只记一个公式，记 A(η) 的导数是充分统计量的期望。后面 ELBO 对变分参数的梯度几乎都从这里长出来。

2.2 从因子图到指数族因子：把联合分布拆开

图模型的核心操作是把联合概率拆成因子。有向图对应条件概率的乘积，无向图对应势函数的乘积再归一化。以隐马尔可夫模型为例，联合分布是

p(x, z) = p(z₁) ∏ p(zₜ | zₜ₋₁) ∏ p(xₜ | zₜ)

每个条件概率如果都写成指数族，那么整条链上的推断就可以用前向-后向算法统一处理。变分推断在这里做的事，是用一个更简单的分布 q(z) 去近似 p(z | x)，而 q 通常也选指数族，比如平均场假设下 q(z) = ∏ qᵢ(zᵢ)，每个 qᵢ 都是指数族。

选指数族做 q 的理由很实际：指数族的 KL 散度有解析形式，坐标上升时每个 qᵢ 的更新可以写成“期望的充分统计量”形式，不需要数值积分。我一般会先把模型里所有变量按是否观测分成两组，观测变量固定，隐变量用指数族 q 近似，然后逐变量更新。

2.3 变分推断的目标：ELBO 的指数族写法

变分推断最大化的是证据下界 ELBO：

L(q) = E_q[ log p(x, z) ] − E_q[ log q(z) ]

当 p 和 q 都是指数族时，这两项都能写成自然参数和期望充分统计量的内积。具体地，如果 q(z) 的自然参数是 λ，那么 ELBO 对 λ 的梯度就是

∇_λ L = E_q[ T(z) ] 的某种组合 − ∇_λ A(λ)

这个形式让坐标上升变得非常机械：固定其他变量，对当前变量的自然参数求导置零，解出更新式。下面用一个小例子把这条路径跑通。

3. 用 Python 复现平均场变分推断：从公式到可运行代码

3.1 选一个最小模型：贝叶斯高斯混合的简化版

为了不陷进符号泥潭，我选一个能手动推导、又能写成代码的最小模型：两个一维高斯组成的混合，隐变量是每个数据点的簇分配 z，参数是簇均值 μ₁、μ₂ 和精度 τ。观测数据 x 已知，目标是近似 p(z, μ, τ | x)。

平均场假设 q(z, μ, τ) = q(z) q(μ) q(τ)。每个 q 都取指数族：q(z) 是类别分布，q(μ) 是高斯，q(τ) 是伽马。这样每个因子的更新都能写成闭式。

3.2 初始化与数据生成

先造一批可复现的数据，再初始化变分参数。代码里所有随机种子固定，方便你对照结果。

import numpy as np np.random.seed(42) # 生成两个簇的数据，各 100 个点 N = 200 true_mu = np.array([-2.0, 3.0]) true_tau = 1.0 z_true = np.random.choice([0, 1], size=N, p=[0.5, 0.5]) x = np.random.normal(true_mu[z_true], 1.0 / np.sqrt(true_tau)) # 变分参数初始化 # q(z) 用责任度 r[n, k] 表示 r = np.random.rand(N, 2) r = r / r.sum(axis=1, keepdims=True) # q(mu_k) 是高斯，初始化均值和精度 m_mu = np.array([-1.0, 1.0]) beta_mu = np.array([1.0, 1.0]) # 精度 # q(tau) 是伽马，初始化形状和率 a_tau = np.array([1.0, 1.0]) b_tau = np.array([1.0, 1.0])

这段代码做了三件事：生成两个高斯簇的观测数据；用责任度矩阵 r 表示 q(z)；用均值和精度表示 q(μ)，用伽马参数表示 q(τ)。参数含义很直接：m_mu 是每个簇均值的变分均值，beta_mu 是变分精度，a_tau 和 b_tau 控制精度分布的形状。

3.3 坐标上升更新：逐个因子推公式

平均场坐标上升的更新顺序是：固定 q(μ) 和 q(τ)，更新 q(z)；再固定 q(z) 和 q(τ)，更新 q(μ)；最后更新 q(τ)。每步都用到指数族的期望充分统计量。

def update_r(x, m_mu, beta_mu, a_tau, b_tau): # 计算 log rho[n, k] = E[log pi_k] + 0.5 E[log tau_k] - 0.5 E[tau_k] (x_n - mu_k)^2 E_log_tau = np.log(b_tau) - np.digamma(a_tau) # 注意：这里用 digamma 近似 E_tau = a_tau / b_tau log_rho = np.zeros((len(x), 2)) for k in range(2): log_rho[:, k] = ( 0.5 * E_log_tau[k] - 0.5 * E_tau[k] * (x ** 2 - 2 * x * m_mu[k] + m_mu[k] ** 2 + 1.0 / beta_mu[k]) ) # 减去 logsumexp 做归一化 log_rho -= log_rho.max(axis=1, keepdims=True) r = np.exp(log_rho) r = r / r.sum(axis=1, keepdims=True) return r def update_mu(x, r, a_tau, b_tau): E_tau = a_tau / b_tau Nk = r.sum(axis=0) m_mu = (r * x[:, None]).sum(axis=0) / Nk beta_mu = Nk * E_tau return m_mu, beta_mu def update_tau(x, r, m_mu, beta_mu): Nk = r.sum(axis=0) a_tau = Nk / 2.0 + 1.0 b_tau = 0.5 * (r * (x[:, None] ** 2 - 2 * x[:, None] * m_mu + m_mu ** 2 + 1.0 / beta_mu)).sum(axis=0) + 1.0 return a_tau, b_tau

这三个函数对应三个因子的更新。update_r 里用到了 E[log τ] 和 E[τ]，前者用 digamma 函数，后者是 a/b。update_mu 里 m_mu 是加权均值，beta_mu 是责任度之和乘以 E[τ]。update_tau 里 a_tau 和 b_tau 的更新来自伽马分布的共轭更新。注意 b_tau 里多加了 1.0，这是先验的贡献，实际使用时按你的先验参数调整。

3.4 迭代与收敛判断

把三个更新串起来，跑若干轮，观察 ELBO 或参数变化。

def elbo(x, r, m_mu, beta_mu, a_tau, b_tau): E_log_tau = np.log(b_tau) - np.digamma(a_tau) E_tau = a_tau / b_tau Nk = r.sum(axis=0) # 近似 ELBO，忽略常数项 term1 = 0.5 * (E_log_tau * Nk).sum() term2 = -0.5 * (E_tau * (r * (x[:, None] ** 2 - 2 * x[:, None] * m_mu + m_mu ** 2 + 1.0 / beta_mu)).sum(axis=0)).sum() term3 = - (r * np.log(r + 1e-12)).sum() return term1 + term2 + term3 for it in range(50): r = update_r(x, m_mu, beta_mu, a_tau, b_tau) m_mu, beta_mu = update_mu(x, r, a_tau, b_tau) a_tau, b_tau = update_tau(x, r, m_mu, beta_mu) if it % 10 == 0: print(f"iter {it}, ELBO={elbo(x, r, m_mu, beta_mu, a_tau, b_tau):.2f}, m_mu={m_mu}")

跑完 50 轮，m_mu 应该收敛到接近 -2 和 3。如果没收敛，先检查 r 的归一化，再看 b_tau 的更新是否漏了先验项。这个最小实现虽然简单，但把指数族、平均场、坐标上升三件事串起来了，改模型时只需要替换对应的更新函数。

4. 避坑与排查：变分推断落地时最容易翻车的五个地方

4.1 现象：ELBO 震荡不收敛，责任度来回跳

原因通常是更新顺序或初始化不当。平均场坐标上升对初始化敏感，如果两个簇的初始均值太近，责任度会在两个簇之间反复横跳。解决方法是多组初始化取 ELBO 最高的那组，或者先用 k-means 给一个粗略初始。

4.2 现象：更新公式推出来和代码对不上，结果差一个常数

指数族里 A(η) 的常数项很容易漏。比如伽马分布的 E[log τ] 是 digamma(a) − log(b)，不是 log(b) − digamma(a)。符号反了会导致 ELBO 一直下降。建议每推一个期望，先用数值积分验证一遍。

4.3 现象：q(τ) 的更新里 b_tau 越来越小，精度爆炸

这是先验太弱或数据点太少导致的。伽马分布的形状参数 a_tau 如果小于 1，分布会集中在 0 附近，精度估计不稳定。解决方法是给 a_tau 和 b_tau 设一个合理的先验，比如 a_tau=1, b_tau=1，并在更新时保留先验项。

4.4 现象：ELBO 计算出来是正的，但理论上应该是负的

ELBO 的绝对值没有意义，只有相对变化有意义。如果你把常数项随便扔了，ELBO 可能变成正数。只要它在迭代中单调上升，就说明更新方向是对的。不要纠结符号。

4.5 现象：换一个数据集，收敛速度差很多

变分推断的收敛速度和数据的尺度有关。如果 x 的方差很大，E[τ] 的初始值需要相应调整。我一般会先把数据标准化到零均值单位方差，再跑变分更新，最后把参数变换回去。这一步能省掉很多调参时间。

5. 进阶技巧：用指数族形式统一检查你的变分更新

5.1 把更新式写成自然参数形式

当你推导出多个模型的变分更新后，会发现它们长得越来越像。把每个 q 的自然参数写成 λ，更新式往往可以统一成

λ_new = 先验自然参数 + 数据充分统计量的期望

这个形式在指数族共轭模型里几乎总是成立。我习惯在写完更新函数后，回头检查它是否符合这个结构。如果不符合，大概率是某个期望算错了。

5.2 用有限差分验证梯度

如果你不想手推 ELBO 的梯度，可以用有限差分做数值验证。对每个变分参数加一个小扰动，看 ELBO 的变化是否和解析梯度一致。这个方法在调试复杂模型时特别管用，能快速定位是公式错还是代码错。

def numerical_grad(param, func, eps=1e-5): grad = np.zeros_like(param) for i in range(len(param)): param[i] += eps f_plus = func(param) param[i] -= 2 * eps f_minus = func(param) param[i] += eps grad[i] = (f_plus - f_minus) / (2 * eps) return grad

这个函数对一维参数有效，多维参数需要展平后处理。验证时先固定其他变量，只对当前变量做差分，和解析更新式对比。

5.3 一个我常犯的错误

早期我总想把所有隐变量的更新写成矩阵运算，结果维度对不上时排查很久。后来我养成一个习惯：先用循环写一版能跑通的，确认 ELBO 上升后，再改写成向量化形式，每改一步都对比数值。这样虽然慢一点，但不会在维度变换里迷路。变分推断的代码一旦维度错了，报错信息往往指向别处，血泪经验就是先慢后快。

希望帮到你。

本文还有配套的精品资源，点击获取