图模型学习算法解析：迭代缩放、伪似然与评分匹配的工程实践

1. 图模型学习：从理论到实践的深度解析

在概率图模型的世界里，我们常常面对一个核心挑战：如何从一个复杂的、高维的联合概率分布中，学习到其内在的结构和参数？这不仅仅是理论上的优雅问题，更是许多实际应用——从自然语言处理中的依存句法分析，到计算机视觉中的图像分割，再到生物信息学里的基因调控网络推断——所必须解决的工程难题。图模型，特别是马尔可夫随机场，为我们提供了一种强大的框架来描述变量之间的依赖关系，但随之而来的学习过程却充满了计算上的“陷阱”。

传统的最大似然估计在面对图模型时，往往会撞上“配分函数”这堵高墙。这个用于归一化概率分布的和或积分，在模型复杂时几乎无法精确计算。这就好比你想知道一个巨大迷宫里所有路径的总长度，却只能一条条去走完。因此，从业者们发展出了一系列巧妙的“迂回战术”。迭代缩放算法像一位耐心的雕塑家，通过逐步、局部的调整，将一块均匀的“概率黏土”雕刻成符合数据统计特征（矩匹配）的形状，其核心思想是最大熵原理——在满足已知约束的条件下，选择最“平坦”、最不确定的分布。伪似然方法则更“务实”一些，它承认直接搞定整个联合分布太难，于是转而优化一系列局部的条件概率，这些局部问题往往容易得多。而评分匹配则为连续变量领域打开了一扇窗，它不直接比较概率密度本身（因为密度值可能难以计算），而是去匹配概率对数密度的梯度场（即“评分函数”），这在处理基于能量的模型时尤为有效。

这些方法不仅仅是数学公式的堆砌，它们代表了我们在计算可行性与统计效率之间寻找平衡的持续努力。接下来，我们将深入这些算法的内部，拆解其设计思路、实现细节，并分享在实际应用中积累下来的那些“教科书上不会写”的经验与教训。

2. 核心原理与算法思想拆解

2.1 最大熵与迭代缩放：约束下的最均匀分布

迭代缩放算法的根源在于最大熵原理。它的哲学很直观：当我们从数据中观测到一些统计特征（例如，某些特征的期望值）后，在所有满足这些特征的候选概率分布中，我们应该选择那个“最均匀”、最没有额外假设的分布。这个分布就是最大熵分布。

为什么选择最大熵？因为它是最保守的推断。在没有更多信息的情况下，引入任何超出约束条件的结构都是没有根据的。在指数族分布的框架下，给定一组充分统计量 $U(x) = (U_1(x), ..., U_d(x))$ 及其在数据上的经验期望 $u = (u_1, ..., u_d)$，最大熵分布恰好具有指数族的形式：$\pi_\theta(x) \propto \exp(-\theta^T U(x))$。这里的参数 $\theta$ 就是拉格朗日乘子，其物理意义是让分布的期望 $E_{\pi_\theta}[U(X)]$ 等于观测值 $u$。

注意：最大熵分布的存在性是有条件的，即观测到的矩 $u$ 必须在指数族模型所能实现的期望值集合的内部。如果 $u$ 落在边界上，最大似然估计可能趋于无穷，这在算法实现中需要特别检查。

迭代缩放算法就是为了求解这个最大熵分布（即最大似然估计）而设计的。它的巧妙之处在于处理了 $U(x)$ 的一个特殊性质：可以将其视为一个随机概率向量。也就是说，对于每一个样本 $x$，$U_j(x)$ 可以被解释为一个概率值，满足 $\sum_{j=1}^d U_j(x) = 1$ 且 $U_j(x) \ge 0$。许多常见的图模型特征（如边上的势函数、节点的状态指示函数）都可以通过归一化满足这个条件。

算法的核心迭代步骤（公式18.13的精简表述）是： $$\theta^{(n+1)}j = \theta^{(n)}j + \log \left( \frac{E{\pi{\theta^{(n)}}}[U_j]}{u_j} \right) + \text{KL}(u | E_{\pi_{\theta^{(n)}}}[U])$$

这个公式的直观理解是什么？

1. 比例项：$\log(E[U_j] / u_j)$ 是驱动更新的主要力量。如果模型当前对第 $j$ 个特征的期望 $E[U_j]$ 小于观测值 $u_j$，那么这一项为负，$\theta_j$ 会减小。在指数族中，减小 $\theta_j$ 会使分布更倾向于产生较大的 $U_j(x)$，从而在下次迭代中提升 $E[U_j]$，使其向 $u_j$ 靠近。
2. KL散度项：$\text{KL}(u | E[U])$ 是一个与 $j$ 无关的常数。它的作用是作为一个整体的偏移量，确保在迭代过程中参数向量 $\theta$ 满足一个额外的线性约束（例如 $\sum_j u_j \theta_j = 0$），从而解决模型因 $U(x)$ 归一化而带来的参数冗余问题。没有这一项，算法可能无法收敛到唯一解。

每一次迭代，算法都保证会减少当前分布 $\pi_{\theta^{(n)}}$ 与真实最大熵分布 $\pi^*$ 之间的KL散度。从几何上看，它是在沿着一个保证下降的方向，在概率分布的流形上行走，最终收敛到目标。

2.2 伪似然：化整为零的局部策略

当图模型变得复杂（比如带有环的马尔可夫随机场），计算配分函数 $Z(\theta)$ 和全局的似然函数变得不可行。伪似然方法提供了一个计算上可行的替代方案。

它的思想是用一系列局部条件概率的乘积来近似联合概率。对于一个定义在顶点集 $V$ 上的随机场 $X$，其伪似然函数定义为： $$PL(\theta) = \prod_{s \in V} \prod_{k=1}^N \pi_\theta (x_k^{(s)} | x_k^{(t)}, t \neq s)$$ 其中，$\pi_\theta (x^{(s)} | x^{(t)}, t \neq s)$ 是给定其他所有顶点取值时，顶点 $s$ 取值为 $x^{(s)}$ 的条件概率。

为什么这样做是合理的？

1. 计算可行性：对于大多数图模型，尤其是成对交互的马尔可夫随机场，这个条件概率的形式非常简单。它只依赖于顶点 $s$ 的邻居变量，其归一化分母只涉及顶点 $s$ 所有可能取值的求和，计算量远小于全局的配分函数。例如，对于一个二值随机场，这个求和只有两项。
2. 统计一致性：在样本量 $N \to \infty$ 时，最大化伪似然得到的估计量 $\hat{\theta}{PL}$ 通常是真实参数 $\theta^*$ 的一致估计量。也就是说，只要有足够的数据，伪似然估计会收敛到正确的参数。其理论依据在于，伪似然实际上最小化了真实分布 $\pi{\text{true}}$ 与模型 $\pi_\theta$ 在所有局部条件分布上的KL散度之和。
3. 保持凹性：对于指数族分布，对数伪似然函数 $\log PL(\theta)$ 关于参数 $\theta$ 仍然是凹函数。这意味着优化问题没有讨厌的局部最优解，我们可以使用梯度上升等凸优化方法高效地找到全局最优解。

实操心得：伪似然估计的方差通常比最大似然估计要大。这意味着在有限样本下，它的估计可能不够精确。因此，它更适合作为大规模、复杂模型的一个初始估计，或者在其他方法（如MCMC-MLE）中作为重要采样提议分布的参数。在模型选择（比如结构学习）的初期，伪似然因其速度优势也是一个很好的工具。

2.3 评分匹配：绕过配分函数的梯度匹配

对于连续变量模型，特别是形式为 $\pi_\theta(x) = \frac{1}{Z(\theta)} \exp(-F(x; \theta))$ 的基于能量的模型，评分匹配提供了一种独特的学习思路。它完全避开了对配分函数 $Z(\theta)$ 的计算。

评分函数定义为概率对数密度的负梯度：$s(x; \theta) = -\nabla_x \log \pi_\theta(x) = \nabla_x F(x; \theta)$。它刻画了在数据点 $x$ 处，概率密度的局部变化方向。评分匹配的目标是让模型评分函数 $s(x; \theta)$ 尽可能接近真实数据分布评分函数 $s_{\text{true}}(x)$。

其目标函数是最小化平方误差的期望： $$J(\theta) = \frac{1}{2} \int p_{\text{true}}(x) | s(x; \theta) - s_{\text{true}}(x) |^2 dx$$ 其中 $p_{\text{true}}(x)$ 是真实数据的未知密度。

关键技巧在于，经过分部积分（或散度定理），这个依赖于未知 $s_{\text{true}}(x)$ 的目标函数，可以转化为一个只依赖于模型和数据的表达式： $$J(\theta) = \int p_{\text{true}}(x) \left[ \text{trace}(\nabla_x s(x; \theta)) + \frac{1}{2} | s(x; \theta) |^2 \right] dx + \text{const}$$ 这里 $\nabla_x s(x; \theta)$ 是评分函数 $s$ 的雅可比矩阵，$\text{trace}$ 是求迹运算。常数项与 $\theta$ 无关，可以忽略。

因此，我们可以用样本平均来近似这个期望，得到最终的优化目标： $$\hat{J}(\theta) = \frac{1}{N} \sum_{k=1}^N \left[ \text{trace}(\nabla_x s(x_k; \theta)) + \frac{1}{2} | s(x_k; \theta) |^2 \right]$$

评分匹配的优势与挑战：

• 优势：完全不需要计算或估计配分函数 $Z(\theta)$，也无需从模型中进行采样（与对比散度等算法不同）。它只需要计算评分函数 $s(x; \theta)$ 及其雅可比矩阵的迹，对于许多参数化模型（如神经网络），这些可以通过自动微分高效完成。
• 挑战：计算雅可比矩阵的迹 $\text{trace}(\nabla_x s(x; \theta))$ 在数据维度 $d$ 很高时，计算成本是 $O(d^2)$，这可能非常昂贵。为此，后续发展出了切片评分匹配、去噪评分匹配等变体来缓解这个问题。

3. 算法实现细节与实操要点

3.1 迭代缩放算法的工程实现

理论很优美，但将迭代缩放算法投入实际运行，需要处理一系列工程细节。

1. 特征函数 $U(x)$ 的设计与归一化： 算法的前提是 $U_j(x) \ge 0$ 且 $\sum_j U_j(x) = 1$。对于一般的特征，我们需要进行变换。假设原始特征向量为 $\tilde{U}(x)$，其每个分量有下界 $\tilde{u}^-$ 和上界 $\tilde{u}^+$（或估计一个上界）。我们可以进行如下线性变换：

1. 平移：$U'(x) = \tilde{U}(x) - \tilde{u}^- \cdot \mathbf{1}$，确保非负。
2. 缩放：$U''(x) = U'(x) / C$，其中 $C$ 是一个足够大的常数（如 $\max_x \sum_j U'_j(x)$），确保 $\sum_j U''_j(x) \le 1$。
3. 补全：定义第 $d+1$ 个特征 $U_{d+1}(x) = 1 - \sum_{j=1}^d U''_j(x)$。 这样，新的特征向量 $(U''_1(x), ..., U''d(x), U{d+1}(x))$ 就满足了算法要求。对应的经验矩 $u$ 也需要进行同样的变换。

2. 期望 $E_{\pi_\theta}[U]$ 的计算——MCMC采样： 迭代缩放每一步都需要计算当前模型下特征 $U$ 的期望 $E_{\pi_\theta}[U]$。对于复杂的图模型，这个期望没有解析解，必须通过蒙特卡洛方法估计。最常用的是马尔可夫链蒙特卡洛方法。

• 吉布斯采样：适用于条件分布容易采样的模型，如受限玻尔兹曼机、条件随机场。它逐个变量地依据其条件分布进行更新。
• Metropolis-Hastings采样：更通用，需要设计提议分布。对于离散图模型，常采用“单个变量翻转”或“块翻转”作为提议。
• 实操要点：
  • 预热：每次参数 $\theta$ 更新后，马尔可夫链需要重新达到平稳分布。可以从上一轮采样的最终状态继续运行，但需要足够的预热步数以消除初始状态的影响。
  • 采样数：需要在估计精度和计算时间之间权衡。通常，在迭代初期，参数变化大，可以少用一些采样样本；接近收敛时，应增加采样数以获得更精确的梯度方向。
  • 诊断：必须监控MCMC的收敛性，例如通过计算不同链之间的方差、自相关函数等，确保采样结果可靠。

3. 迭代停止准则： 迭代缩放是一个迭代过程，需要明确的停止条件。

1. 参数变化：$|\theta^{(n+1)} - \theta^{(n)}| < \epsilon$。但参数本身可能缩放不同，更推荐使用相对变化。
2. 矩匹配误差：$| E_{\pi_{\theta^{(n)}}}[U] - u | < \epsilon$。这是最直接的准则，反映了算法的核心目标。
3. 似然增长：对数似然（或伪似然）在连续几次迭代中的增长量小于某个阈值。但计算真实似然可能困难。
4. 最大迭代次数：设置一个安全上限，防止无限循环。

一个稳健的实现通常会结合2和4，并辅以详细的日志记录，以便追踪算法进程。

3.2 伪似然估计的优化技巧

最大化伪似然是一个无约束的凸优化问题，可以使用标准的梯度上升或拟牛顿法（如L-BFGS）。

对数伪似然的梯度计算： 对于指数族分布 $\pi_\theta(x) \propto \exp(-\theta^T U(x))$，顶点 $s$ 的条件概率为： $$\pi_\theta(x^{(s)} | x^{(\setminus s)}) = \frac{\exp(-\theta^T U(x))}{\sum_{y^{(s)}} \exp(-\theta^T U(y^{(s)}, x^{(\setminus s)}))}$$ 其对参数 $\theta$ 的梯度为： $$\nabla_\theta \log \pi_\theta(x^{(s)} | x^{(\setminus s)}) = -E_{\pi_\theta(\cdot | x^{(\setminus s)})}[U(X)] + U(x)$$ 其中，期望 $E_{\pi_\theta(\cdot | x^{(\setminus s)})}[\cdot]$ 是在给定其他所有顶点取值时，顶点 $s$ 的条件分布下的期望。这个期望的计算是局部的，只涉及顶点 $s$ 的所有可能状态，计算代价很小。

因此，整个对数伪似然函数的梯度是数据中所有顶点-样本对上的梯度之和。对于大规模数据，我们通常使用随机梯度上升：每次随机选取一个样本（或一小批样本），计算其梯度并更新参数。

正则化与模型选择： 为了避免过拟合，尤其是在变量多、样本少的情况下，需要在目标函数中加入正则化项。

• L2正则化（岭回归）：在目标函数中增加 $-\frac{\lambda}{2}|\theta|^2$，倾向于产生较小、较分散的参数，提高模型稳定性。
• L1正则化（Lasso）：增加 $-\lambda |\theta|1$。这倾向于产生稀疏解，即许多参数被精确地推向零。在图模型学习中，这对应着结构学习——如果一个交互项的系数 $\theta{st}$ 为零，就意味着图中节点 $s$ 和 $t$ 之间没有边。L1正则化的伪似然最大化是学习图结构的一种常用且高效的方法。

避坑指南：伪似然估计的梯度计算涉及大量条件期望。在实现时，务必利用模型的局部结构进行向量化操作。例如，对于成对交互的随机场，每个条件期望的计算可以并行化。同时，注意数值稳定性，在计算 softmax（即指数归一化）时，使用log-sum-exp技巧来避免数值上溢或下溢。

3.3 评分匹配的实现与变体

实现评分匹配，核心是高效计算目标函数 $\hat{J}(\theta)$ 及其梯度。

目标函数计算： 对于模型 $\pi_\theta(x) = \frac{1}{Z(\theta)} \exp(-F(x; \theta))$，评分函数 $s(x; \theta) = \nabla_x F(x; \theta)$。我们需要计算：

1. $s(x_k; \theta)$：一次前向传播和梯度计算（相对于输入 $x$）。
2. $| s(x_k; \theta) |^2$：向量点积。
3. $\text{trace}(\nabla_x s(x_k; \theta))$：这是主要计算瓶颈。$\nabla_x s(x_k; \theta)$ 是 $d \times d$ 的雅可比矩阵。

高效计算迹的技巧： 直接计算雅可比矩阵再求迹的代价是 $O(d^2)$。有两种主流方法可以将其降至 $O(d)$：

1. Hutchinson 迹估计器：利用恒等式 $\text{trace}(A) = E_{v}[v^T A v]$，其中 $v$ 是一个随机向量，满足 $E[vv^T] = I$（例如，$v$ 的元素为独立同分布的 Rademacher 随机变量，即等概率取 $+1$ 或 $-1$）。这样，我们只需要计算矩阵 $A$ 作用于向量 $v$ 的结果 $Av$，然后计算 $v^T (Av)$ 即可。这只需要一次额外的向量-雅可比乘积计算。

   # 伪代码示例 (使用 PyTorch 风格) def score_matching_loss(model, x_batch): x_batch.requires_grad_(True) energy = model(x_batch) # F(x; theta) score = torch.autograd.grad(energy.sum(), x_batch, create_graph=True)[0] # s(x; theta) # 计算 ||s||^2 loss1 = 0.5 * (score ** 2).sum(dim=1).mean() # 使用 Hutchinson 估计器计算 trace(Jacobian) v = torch.randn_like(x_batch) # 随机向量 vJ = torch.autograd.grad(score, x_batch, grad_outputs=v, create_graph=True)[0] # 计算 Jacobian-vector product: J * v loss2 = (v * vJ).sum(dim=1).mean() # 估计 trace(J) = E[v^T J v] total_loss = loss1 + loss2 return total_loss

2. 自动微分库的迹计算：现代深度学习框架（如 PyTorch、JAX）的自动微分系统在计算标量函数对向量的梯度时，其反向传播过程本质上就计算了雅可比矩阵的迹。具体来说，$\text{trace}(\nabla_x s(x; \theta))$ 等价于对 $s(x; \theta)$ 的每一个分量 $s_i$ 求其对 $x_i$ 的偏导数，然后求和。这可以通过构造一个特殊的“求和”操作来实现。

评分匹配的变体：

• 切片评分匹配：为了进一步降低计算和统计方差，不匹配整个评分向量，而是随机投影到某个方向 $v$ 上，匹配投影后的评分：$E_{p_{\text{true}}}[(v^T s(x; \theta) - v^T s_{\text{true}}(x))^2]$。通过多次随机投影取平均来近似原目标。
• 去噪评分匹配：这是当前生成模型（如扩散模型）的理论基石之一。其核心思想是，对一个干净数据点 $x$ 添加一个小的高斯噪声 $\epsilon \sim N(0, \sigma^2 I)$ 得到噪声数据 $\tilde{x} = x + \epsilon$。然后，模型学习去预测这个噪声，即学习一个函数 $s_\theta(\tilde{x}, \sigma)$ 来近似 $\nabla_{\tilde{x}} \log p_\sigma(\tilde{x})$，其中 $p_\sigma$ 是加噪后的数据分布。通过在不同噪声水平 $\sigma$ 上训练，模型可以学习到数据分布在整个噪声尺度上的评分函数，进而通过逆向随机微分方程从噪声中生成数据。

4. 算法对比、选择与混合策略

在实际项目中，没有一种算法是万能的。选择哪种学习算法，取决于模型的具体形式、数据特性以及计算资源。

4.1 方法对比速查表

4.2 混合与进阶策略

资深从业者不会死守一种方法，而是根据情况灵活组合。

1. 伪似然预热 + 迭代缩放精调： 伪似然估计速度快，但可能不够精确。迭代缩放精度高，但需要好的初始点且计算慢。一个常见的策略是：

• 阶段一：使用伪似然估计快速得到一个接近最优解的参数 $\theta_{PL}$。
• 阶段二：以 $\theta_{PL}$ 作为初始值，运行迭代缩放算法进行精调。由于起点好，迭代缩放所需的MCMC采样链可以更短，收敛更快。

2. 对比散度作为迭代缩放的快速期望估计： 在迭代缩放中，最耗时的步骤是用MCMC估计 $E_{\pi_\theta}[U]$。对比散度提供了一种快速的近似方法。CD-$k$ 算法不是运行MCMC链至平稳，而是只运行固定的 $k$ 步（通常 $k=1$ 就有效），从训练数据本身开始（称为“正相”），然后从这 $k$ 步后的状态计算特征期望（“负相”）。虽然是有偏估计，但在实践中对于获取参数更新方向非常有效，广泛应用于受限玻尔兹曼机的训练。

3. 随机梯度上升处理大规模数据与隐变量： 当数据量巨大或模型包含隐变量时，标准的批量算法内存可能不足。如原文第18.3.2节所述，可以采用随机梯度上升。其更新规则为： $$\theta^{(n+1)} = \theta^{(n)} + \gamma_n \left[ U(\tilde{x}^{(n+1)}) - \frac{1}{N} \sum_{k=1}^N E_{\theta^{(n)}}[U(X) | X^{(S)} = x_k^{(S)}] \right]$$ 其中 $\tilde{x}^{(n+1)}$ 是从当前模型 $\pi_{\theta^{(n)}}$ 中采样的一个完整配置（可通过MCMC获得），而条件期望 $E[U | X^{(S)}]$ 则需要针对每个观测数据 $x_k^{(S)}$ 进行推断（例如，运行条件MCMC）。这本质上是一种在线EM算法。学习率 $\gamma_n$ 需要满足Robbins-Monro条件（$\sum \gamma_n = \infty, \sum \gamma_n^2 < \infty$），例如 $\gamma_n = a / (b + n)$。

实操心得：学习率调参：SGA的学习率调度至关重要。开始时可设置较大的学习率以快速接近解，后期必须减小以保证稳定收敛。监控目标函数（如伪似然）在验证集上的变化是调整学习率计划的主要依据。当目标函数剧烈震荡时，说明学习率太大；当长期停滞时，可能太小或已收敛。

5. 实战中的常见问题与排查技巧

5.1 数值不稳定与溢出/下溢

指数族模型的计算大量涉及指数运算 $exp(-\theta^T U(x))$，极易导致数值问题。

• 问题：exp(1000)导致上溢（Inf），exp(-1000)导致下溢（0），进而导致对数运算log(0)产生-Inf。
• 解决方案：
  1. Log-Sum-Exp 技巧：计算 $\log \sum_i \exp(z_i)$ 时，使用公式：$\log \sum_i \exp(z_i) = m + \log \sum_i \exp(z_i - m)$，其中 $m = \max_i(z_i)$。这确保了指数运算的参数范围可控。
  2. 参数标准化：在迭代过程中，定期对参数 $\theta$ 进行平移。因为对于指数族 $\exp(-\theta^T U(x))$，同时给所有 $\theta$ 加上一个常数 $c$，相当于给配分函数乘以 $\exp(c)$，不影响概率分布。可以强制让 $\theta$ 的均值或某一部分为0，以控制其幅度。
  3. 使用高精度浮点数：在关键计算部分（如计算配分函数的对数）使用float64（双精度）而非float32（单精度）。

5.2 模型不收敛或收敛缓慢

• 可能原因1：学习率不当（对于梯度类方法）。
  • 排查：绘制目标函数（似然、伪似然）随迭代次数的变化图。如果曲线剧烈震荡，学习率太大；如果几乎是一条水平线，学习率太���或算法已卡住。
  • 解决：使用自适应学习率算法（如Adam, Adagrad），或实现学习率衰减计划。对于迭代缩放，检查MCMC采样是否充分，不准确的期望估计会导致更新方向错误，类似学习率过大。
• 可能原因2：MCMC混合时间过慢（影响迭代缩放和SGA）。
  • 现象：特征期望 $E[U]$ 的估计值波动很大，即使增加采样数也不稳定。
  • 排查：计算MCMC采样序列的自相关性。如果自相关衰减很慢，说明链混合不好。
  • 解决：
    • 改进提议分布：在Metropolis-Hastings中，尝试不同的提议机制（如块更新、Swendsen-Wang算法用于Potts模型）。
    • 提高温度：暂时在采样时使用一个“退火”版本，即从分布 $\pi_\theta(x)^{1/T}$ 中采样，$T>1$ 会使分布更平坦，混合更快，但会引入偏差。可以逐渐将 $T$ 降至1。
    • 并行多条链：从不同的随机初始状态运行多条MCMC链，比较它们的统计量，以确保采样来自平稳分布。
• 可能原因3：模型不可识别或数据不足。
  • 现象：参数变化很大，但似然值几乎不变；或者不同运行得到完全不同的参数估计。
  • 排查：检查特征 $U(x)$ 是否线性相关，或者某些特征的观测方差是否为0。这会导致最大似然解不唯一（似然函数存在平坦区域）。
  • 解决：加入正则化项（L1/L2）；如果可能，增加数据量；重新设计或筛选特征。

5.3 隐变量模型与EM算法的困境

当图模型包含隐变量时，标准的EM算法（如原文18.3.1节）需要在E步计算隐变量的后验期望，这通常涉及难以处理的推断。

• 近似E步：不要追求精确的后验计算，转而使用近似方法。
  • 变分推断：用另一个简单的分布 $q(H)$ 来近似真实后验 $p(H|S)$，并通过优化KL散度 $KL(q|p)$ 来更新 $q$。这将E步转化为一个优化问题。
  • MCMC-EM：在E步，用MCMC从后验 $p(H|S, \theta^{old})$ 中采样，然后用样本平均来近似期望。这就是随机EM的一种形式。
  • 对比散度：对于像RBM这样的特定模型，CD算法可以看作是对对数似然梯度中难处理项的一个快速近似，本质上是EM算法中E步的一种近似。
• 绕过EM：直接优化边缘似然：对于连续隐变量模型，如果隐变量与观测变量的联合分布是易处理的，有时可以通过重参数化技巧和自动微分，直接对边缘似然 $p(S|\theta)$ 进行随机梯度上升。这即是变分自编码器中使用的策略。

5.4 评估与诊断：如何知道模型学得好不好？

学习完成后，必须评估模型质量。

1. 似然与伪似然：在测试集上计算对数似然（如果可算）或对数伪似然。更高的值通常意味着更好的拟合。但要注意过拟合：训练集似然高，测试集似然低。
2. 特征矩匹配：计算模型生成的数据（通过采样获得）的特征期望 $E_{\text{model}}[U]$，与训练集的经验期望 $E_{\text{data}}[U]$ 进行比较。对于迭代缩放和最大熵模型，这是直接的评估指标。可以绘制散点图或计算均方误差。
3. 生成样本可视化：对于图像、文本等数据，最直观的评估是查看模型生成的样本。检查样本是否清晰、多样，且与训练数据类似。如果样本模糊或模式单一，说明模型可能欠拟合或崩溃。
4. 下游任务性能：如果模型用于某个具体任务（如分类、分割），最终应以其在任务上的性能（准确率、F1分数等）为评判标准。

图模型的学习是一个融合了统计理论、优化算法和工程实践的领域。从最大熵的哲学出发，到迭代缩放、伪似然、评分匹配这些精巧的算法，再到处理数值问题、加速收敛、应对隐变量的各种技巧，每一步都需要对模型本质和计算代价有深刻的理解。没有最好的算法，只有最适合当前问题和资源的策略。掌握这些方法的原理、实现细节以及它们之间的权衡，是构建高效、鲁棒的概率图模型应用的关键。在实际操作中，我习惯从简单的伪似然估计开始，快速验证模型结构和特征的有效性，然后再根据需求决定是否投入更多计算资源，使用迭代缩放或MCMC-based方法进行精调。对于连续数据生成任务，评分匹配及其现代变体（扩散模型）已成为我的首选工具。记住，耐心调试和系统性评估与算法设计本身同等重要。