★ FEATURED ARTICLE 简介这份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 的指数族写法变分推断最大化的是证据下界 ELBOL(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], sizeN, 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(axis1, keepdimsTrue) # 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(axis1, keepdimsTrue) r np.exp(log_rho) r r / r.sum(axis1, keepdimsTrue) return r def update_mu(x, r, a_tau, b_tau): E_tau a_tau / b_tau Nk r.sum(axis0) m_mu (r * x[:, None]).sum(axis0) / Nk beta_mu Nk * E_tau return m_mu, beta_mu def update_tau(x, r, m_mu, beta_mu): Nk r.sum(axis0) 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(axis0) 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(axis0) # 近似 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(axis0)).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(fiter {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_tau1, b_tau1并在更新时保留先验项。4.4 现象ELBO 计算出来是正的但理论上应该是负的ELBO 的绝对值没有意义只有相对变化有意义。如果你把常数项随便扔了ELBO 可能变成正数。只要它在迭代中单调上升就说明更新方向是对的。不要纠结符号。4.5 现象换一个数据集收敛速度差很多变分推断的收敛速度和数据的尺度有关。如果 x 的方差很大E[τ] 的初始值需要相应调整。我一般会先把数据标准化到零均值单位方差再跑变分更新最后把参数变换回去。这一步能省掉很多调参时间。5. 进阶技巧用指数族形式统一检查你的变分更新5.1 把更新式写成自然参数形式当你推导出多个模型的变分更新后会发现它们长得越来越像。把每个 q 的自然参数写成 λ更新式往往可以统一成λ_new 先验自然参数 数据充分统计量的期望这个形式在指数族共轭模型里几乎总是成立。我习惯在写完更新函数后回头检查它是否符合这个结构。如果不符合大概率是某个期望算错了。5.2 用有限差分验证梯度如果你不想手推 ELBO 的梯度可以用有限差分做数值验证。对每个变分参数加一个小扰动看 ELBO 的变化是否和解析梯度一致。这个方法在调试复杂模型时特别管用能快速定位是公式错还是代码错。def numerical_grad(param, func, eps1e-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 上升后再改写成向量化形式每改一步都对比数值。这样虽然慢一点但不会在维度变换里迷路。变分推断的代码一旦维度错了报错信息往往指向别处血泪经验就是先慢后快。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
