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

变分推断(Variational Inference)将贝叶斯推断转化为优化问题——当后验分布难以精确计算时，变分推断通过优化一个易处理的分布族来近似真实后验，将积分问题转化为优化问题。从变分 EM 到变分自编码器，从主题模型到深度生成模型，变分推断已成为现代机器学习的核心技术。ELBO（证据下界）提供了优化目标，平均场近似简化了分布族，坐标上升算法实现了高效求解，而黑盒变分推断则通过蒙特卡洛梯度估计扩展到任意模型。

Figure 4

Figure 5

贝叶斯推断与后验难题

贝叶斯推断框架

观测数据：

隐变量：

参数：

目标：计算后验分布

困难：边缘似然(证据)

通常无法解析计算，也难以数值积分（高维）。

精确推断 vs 近似推断

精确推断： - 共轭先验：某些模型后验有闭式解 - 图模型：变量消除、信念传播（树结构）

近似推断（大多数情况需要）： 1. 采样方法： MCMC（马尔可夫链蒙特卡洛） - 优点：渐近精确 - 缺点：收敛慢，难以诊断 2. 变分方法：将推断转化为优化 - 优点：快速，确定性 - 缺点：有偏近似

变分推断的基本原理

Figure 2

ELBO 推导

Figure 1

想法：用简单分布近似复杂后验

优化目标：最小化 KL 散度

问题：包含未知的

转换：

其中证据下界(ELBO)：

关键关系：

变分推断目标：

平均场近似

Figure 3

假设：变分分布完全分解

或更简洁地，假设隐变量和参数划分为组：

优化：对每个因子，固定其他因子，最大化 ELBO

坐标上升变分推断

ELBO 展开：

其中是熵。

对优化：固定

最优：

算法：循环更新每个因子直至收敛

变分 EM 算法

EM 与变分推断的联系

标准 EM： - E 步：（精确后验） - M 步：（点估计）

变分 EM： - E 步：（变分近似） - M 步：

变分贝叶斯 EM： - VE 步：变分更新 - VM 步：变分更新（贝叶斯推断参数而非点估计）

变分贝叶斯 GMM

模型：

• 先验：，，
• 似然：，

变分分布：

更新公式（共轭性质）：

：

：，其中，

：正态-Wishart 分布，参数由充分统计量更新（详见 Bishop PRML 10.2 节）

黑盒变分推断(BBVI)

梯度估计问题

ELBO：

梯度：

困难：梯度与期望不可交换（依赖）

REINFORCE 梯度估计器

对数导数技巧：

ELBO 梯度：

蒙特卡洛估计：

其中

问题：高方差

重参数化技巧

想法：将随机性从分离出来

重参数化：，其中是固定分布

示例（高斯）：，

ELBO 梯度：

蒙特卡洛估计：

优势：低方差，可自动微分

实现示例

import numpy as np
from scipy.stats import multivariate_normal, dirichlet
from scipy.special import digamma

class VariationalGMM:
    """变分贝叶斯高斯混合模型（简化版）"""
    def __init__(self, n_components=3, max_iter=100, tol=1e-3):
        self.K = n_components
        self.max_iter = max_iter
        self.tol = tol
        
    def fit(self, X):
        N, d = X.shape
        K = self.K
        
        # 初始化超参数
        alpha0 = 1.0
        m0 = np.mean(X, axis=0)
        beta0 = 1.0
        nu0 = d
        W0 = np.eye(d)
        
        # 初始化变分参数
        self.alpha = np.ones(K) * alpha0 + N / K
        self.beta = np.ones(K) * beta0 + N / K
        self.nu = np.ones(K) * nu0 + N / K
        self.m = np.array([np.mean(X, axis=0) + 0.1 * np.random.randn(d) for _ in range(K)])
        self.W = np.array([W0 for _ in range(K)])
        
        # 初始化责任度
        r = np.random.dirichlet([1] * K, N)
        
        for iteration in range(self.max_iter):
            r_old = r.copy()
            
            # 更新责任度
            r = self._update_r(X, N, d, K)
            
            # 更新参数
            N_k = np.sum(r, axis=0)
            x_bar_k = (r.T @ X) / N_k[:, np.newaxis]
            
            self.alpha = alpha0 + N_k
            self.beta = beta0 + N_k
            self.m = (beta0 * m0 + N_k[:, np.newaxis] * x_bar_k) / self.beta[:, np.newaxis]
            self.nu = nu0 + N_k
            
            for k in range(K):
                S_k = np.zeros((d, d))
                for i in range(N):
                    diff = X[i] - x_bar_k[k]
                    S_k += r[i, k] * np.outer(diff, diff)
                
                diff_m = x_bar_k[k] - m0
                self.W[k] = np.linalg.inv(
                    np.linalg.inv(W0) + N_k[k] * S_k / N_k[k] +
                    (beta0 * N_k[k]) / (beta0 + N_k[k]) * np.outer(diff_m, diff_m)
                )
            
            # 检查收敛
            if np.max(np.abs(r - r_old)) < self.tol:
                break
        
        return self
    
    def _update_r(self, X, N, d, K):
        """更新责任度"""
        r = np.zeros((N, K))
        
        for k in range(K):
            # E[log pi_k]
            E_log_pi = digamma(self.alpha[k]) - digamma(np.sum(self.alpha))
            
            # E[log |Lambda_k|]
            E_log_det = np.sum([digamma((self.nu[k] + 1 - i) / 2) for i in range(1, d + 1)])
            E_log_det += d * np.log(2) + np.log(np.linalg.det(self.W[k]))
            
            # 马氏距离期望
            for i in range(N):
                diff = X[i] - self.m[k]
                E_dist = self.nu[k] * diff @ self.W[k] @ diff + d / self.beta[k]
                r[i, k] = E_log_pi + 0.5 * E_log_det - 0.5 * E_dist
        
        # 归一化
        r = np.exp(r - np.max(r, axis=1, keepdims=True))
        r /= np.sum(r, axis=1, keepdims=True)
        
        return r
    
    def predict(self, X):
        N, d = X.shape
        r = self._update_r(X, N, d, self.K)
        return np.argmax(r, axis=1)

# 黑盒变分推断示例(重参数化)
class BBVI_Gaussian:
    """黑盒变分推断(高斯近似)"""
    def __init__(self, dim, lr=0.01):
        self.mu = np.zeros(dim)
        self.log_sigma = np.zeros(dim)
        self.lr = lr
        
    def sample(self, n_samples=1):
        """重参数化采样"""
        epsilon = np.random.randn(n_samples, len(self.mu))
        return self.mu + np.exp(self.log_sigma) * epsilon
    
    def elbo(self, log_p_func, n_samples=10):
        """估计 ELBO"""
        z_samples = self.sample(n_samples)
        log_p = np.array([log_p_func(z) for z in z_samples])
        log_q = -0.5 * np.sum((z_samples - self.mu) ** 2 / np.exp(2 * self.log_sigma), axis=1)
        log_q -= 0.5 * len(self.mu) * np.log(2 * np.pi) + np.sum(self.log_sigma)
        return np.mean(log_p - log_q)
    
    def step(self, log_p_func, n_samples=10):
        """单步优化(数值梯度)"""
        elbo_current = self.elbo(log_p_func, n_samples)
        
        # 数值梯度(简化)
        eps = 1e-4
        grad_mu = np.zeros_like(self.mu)
        grad_log_sigma = np.zeros_like(self.log_sigma)
        
        for i in range(len(self.mu)):
            self.mu[i] += eps
            grad_mu[i] = (self.elbo(log_p_func, n_samples) - elbo_current) / eps
            self.mu[i] -= eps
            
            self.log_sigma[i] += eps
            grad_log_sigma[i] = (self.elbo(log_p_func, n_samples) - elbo_current) / eps
            self.log_sigma[i] -= eps
        
        # 梯度上升
        self.mu += self.lr * grad_mu
        self.log_sigma += self.lr * grad_log_sigma

if __name__ == '__main__':
    # 变分 GMM 示例
    from sklearn.datasets import make_blobs
    X, _ = make_blobs(n_samples=300, centers=3, n_features=2, random_state=42)
    
    vgmm = VariationalGMM(n_components=3, max_iter=50)
    vgmm.fit(X)
    labels = vgmm.predict(X)
    
    print(f"聚类完成，权重估计: {vgmm.alpha / np.sum(vgmm.alpha)}")

Q&A 精选

Q1: 变分推断 vs MCMC?

A: - 变分: 快速、确定性、有偏(KL 散度非零) - MCMC: 慢、随机、渐近无偏

变分适合大规模数据和在线学习,MCMC 适合精确推断。

---

Q2: 为什么用 KL(q||p)而非 KL(p||q)?

A: KL(q||p)是"反向 KL",使在小的地方也小(零逼近)。 KL(p||q)是"正向 KL",使覆盖的所有模式(矩匹配)。反向 KL 计算上只需可采样,不需要归一化。

---

Q3: 平均场假设何时失效?

A: 当变量强相关时。解决方法: - 结构化变分(保留部分依赖) - 更 rich 的变分族(normalizing flows)

---

Q4: 变分贝叶斯 vs 点估计(MAP/MLE)?

A: 变分贝叶斯保留不确定性,防止过拟合。代价:计算复杂度高。小数据/正则化需求高→变分贝叶斯;大数据/速度需求→点估计。

---

Q5: 重参数化技巧的适用范围?

A: 需要连续可微的分布。适用:高斯、 Logistic 、 Laplace 。不适用:离散分布(需 REINFORCE 或 Gumbel-Softmax)。

---

✏️ 练习题与解答

练习 1：ELBO推导

题目：证明。 解答：，KL非负，得ELBO。

练习 2：平均场近似

题目：，推导。 解答：（固定分布，最大化ELBO）

练习 3：变分EM

题目：变分EM的E步和M步分别优化什么？ 解答：E步：固定优化（最大化ELBO）。M步：固定优化（最大化期望完整对数似然）。

练习 4：VAE重参数化

题目：为什么VAE用而非直接采样？ 解答：直接采样不可微，无法反向传播。重参数化将随机性转移到，可求梯度。

练习 5：变分推断 vs MCMC

题目：何时用变分推断，何时用MCMC？ 解答：变分推断快但有偏（近似假设），适合大数据。MCMC渐近无偏但慢，适合精确推断小规模问题。

参考文献

1. Jordan, M. I., et al. (1999). An introduction to variational methods for graphical models. Machine Learning, 37(2), 183-233.
2. Blei, D. M., Kucukelbir, A., & McAuliffe, J. D. (2017). Variational inference: A review for statisticians. JASA, 112(518), 859-877.
3. Kingma, D. P., & Welling, M. (2014). Auto-encoding variational Bayes. ICLR.
4. Ranganath, R., Gerrish, S., & Blei, D. (2014). Black box variational inference. AISTATS.

---

变分推断将贝叶斯推断的积分难题转化为优化问题,以确定性算法换取计算效率。从经典的平均场近似到现代的黑盒变分推断,从 VAE 到深度生成模型,变分方法已成为机器学习的基础工具。理解变分推断,是通往概率编程、贝叶斯深度学习的必经之路。