目录
  1. 1. 一、引言:当数据不完整时
  2. 2. 二、问题形式化
    1. 2.1. 2.1 观测变量与隐变量
    2. 2.2. 2.2 为什么不直接做梯度上升
  3. 3. 三、EM 算法的核心思想
    1. 3.1. 3.1 Jensen 不等式:EM 的数学基石
    2. 3.2. 3.2 构造证据下界(ELBO)
    3. 3.3. 3.3 E 步与 M 步
  4. 4. 四、Q 函数形式
  5. 5. 五、收敛性证明
    1. 5.1. 5.1 对数似然单调不减
    2. 5.2. 5.2 收敛到局部极值
  6. 6. 六、高斯混合模型(GMM)实例
    1. 6.1. 6.1 模型定义
    2. 6.2. 6.2 完全数据似然
    3. 6.3. 6.3 E 步
    4. 6.4. 6.4 M 步
    5. 6.5. 6.5 GMM EM 算法伪代码
  7. 7. 七、抛硬币数值例子
    1. 7.1. 7.1 数据生成
    2. 7.2. 7.2 迭代过程
    3. 7.3. 7.3 数值稳定性:Log-Sum-Exp 技巧
  8. 8. 八、EM 算法的变体与扩展
    1. 8.1. 8.1 广义 EM 算法(GEM)
    2. 8.2. 8.2 随机 EM(Stochastic EM)
    3. 8.3. 8.3 变分 EM(Variational EM)
    4. 8.4. 8.4 在线 EM(Online EM)
    5. 8.5. 8.5 EM 与 k-means 的关系
  9. 9. 九、EM 的局限性与实践建议
    1. 9.1. 9.1 主要局限
    2. 9.2. 9.2 实践建议
  10. 10. 十、其他应用场景
  11. 11. 十一、面试高频问题
  12. 12. 十二、总结
【统计学习方法死磕系列】EM算法

一、引言:当数据不完整时

在许多机器学习问题中,我们面对的并非完整的观测数据。想象一下:你有一个包含男性和女性身高的数据集,但性别这一列是缺失的。你知道身高服从两个不同的高斯分布,但你不知道每个样本来自哪个分布。这类含有隐变量(latent variable)或缺失数据的问题,正是 EM 算法(Expectation-Maximization Algorithm)大显身手的舞台。

EM 算法并非单一算法,而是一类算法的统称。它是一种通过迭代寻找含有隐变量的概率模型参数极大似然估计(MLE)或极大后验概率估计(MAP)的通用框架。自 1977 年 Dempster、Laird 和 Rubin 将其形式化以来,EM 算法已成为统计学习中最重要、应用最广泛的算法之一。

本文将系统地推导 EM 算法的原理、证明其收敛性,并通过高斯混合模型和抛硬币等具体例子带读者深入理解这一经典算法。

二、问题形式化

2.1 观测变量与隐变量

设我们有一个概率模型,其中:

  • 观测变量(observed variable)$X$:我们能直接观测到的数据
  • 隐变量(latent variable)$Z$:我们无法直接观测,但影响观测数据的变量
  • 参数 $\theta$:模型需要估计的参数

我们的目标是:给定观测数据 $X = {x_1, x_2, \dots, x_N}$,找到参数 $\theta$ 使得观测数据的对数似然最大:

$$ \theta^* = \arg\max_\theta \log P(X|\theta) = \arg\max_\theta \log \sum_Z P(X, Z|\theta) $$

这里的关键难点在于:似然函数中包含了对隐变量 $Z$ 的求和(或积分)操作,再取对数,通常无法直接解析求解。

2.2 为什么不直接做梯度上升

你可能会问:为什么不对 $\theta$ 直接做梯度上升来优化 $\log P(X|\theta)$ ?理论上可行,但在实践中:

  1. 求和在对数内部使得梯度计算复杂:$\nabla_\theta \log \sum_Z P(X,Z|\theta) = \frac{\sum_Z \nabla_\theta P(X,Z|\theta)}{\sum_Z P(X,Z|\theta)}$
  2. 许多模型(如 GMM)的对数似然是非凸的,直接优化容易陷入局部最优
  3. 当 $Z$ 是离散变量且可能取值很多时,$\sum_Z$ 的计算量极大

EM 算法巧妙地绕开这些困难,通过引入一个辅助函数(Q 函数)来迭代优化。

三、EM 算法的核心思想

3.1 Jensen 不等式:EM 的数学基石

EM 算法的推导依赖于一个关键的数学工具——Jensen 不等式

定理(Jensen 不等式):对于凸函数 $f$,有:

$$ f(\mathbb{E}[Y]) \leq \mathbb{E}[f(Y)] $$

特别地,由于 $\log$ 是凹函数(concave),对于任意分布 $Q(Z)$,有:

$$ \log \mathbb{E}_{Q}[g(Z)] \geq \mathbb{E}_{Q}[\log g(Z)] $$

这个不等式告诉我们:先求期望再取对数,大于等于先取对数再求期望。这正是 EM 算法构造下界的关键。

3.2 构造证据下界(ELBO)

对于任意关于隐变量 $Z$ 的分布 $Q(Z)$(满足 $\sum_Z Q(Z) = 1, Q(Z) \geq 0$),我们可以将对数似然重写:

$$ \begin{aligned} \log P(X|\theta) &= \log \sum_Z P(X, Z|\theta) \\ &= \log \sum_Z Q(Z) \frac{P(X, Z|\theta)}{Q(Z)} \\ &= \log \mathbb{E}_{Q(Z)}\left[\frac{P(X, Z|\theta)}{Q(Z)}\right] \\ &\geq \mathbb{E}_{Q(Z)}\left[\log \frac{P(X, Z|\theta)}{Q(Z)}\right] \quad \text{(Jensen不等式)} \end{aligned} $$

这个下界被称为证据下界(Evidence Lower BOund, ELBO),记为 $\mathcal{L}(Q, \theta)$:

$$ \mathcal{L}(Q, \theta) = \sum_Z Q(Z) \log \frac{P(X, Z|\theta)}{Q(Z)} $$

3.3 E 步与 M 步

EM 算法的精妙之处在于:它交替进行两步操作,每一步都使目标函数不降:

E 步(Expectation Step):固定当前参数 $\theta^{(t)}$,选择分布 $Q(Z)$ 使下界 $\mathcal{L}(Q, \theta^{(t)})$ 尽可能大。可以证明,当 $Q(Z) = P(Z|X, \theta^{(t)})$ 时,下界紧贴,即:

$$ \mathcal{L}(P(Z|X, \theta^{(t)}), \theta^{(t)}) = \log P(X|\theta^{(t)}) $$

M 步(Maximization Step):固定 $Q(Z) = P(Z|X, \theta^{(t)})$,优化 $\theta$ 使下界最大化:

$$ \theta^{(t+1)} = \arg\max_\theta \mathcal{L}(P(Z|X, \theta^{(t)}), \theta) $$

由于下界在 $\theta^{(t)}$ 处紧贴,提升下界必然提升真正的对数似然。这一巧妙设计保证了每次迭代对数似然单调不减。

四、Q 函数形式

在实际计算中,我们通常使用 Q 函数 来表达 M 步。Q 函数定义为完全数据对数似然在隐变量后验分布下的期望:

$$ Q(\theta | \theta^{(t)}) = \mathbb{E}_{Z|X,\theta^{(t)}}[\log P(X, Z|\theta)] $$

展开为:

$$ Q(\theta | \theta^{(t)}) = \sum_{i=1}^{N} \mathbb{E}_{z_i|x_i,\theta^{(t)}}[\log P(x_i, z_i|\theta)] $$

M 步等价于最大化 Q 函数:

$$ \theta^{(t+1)} = \arg\max_\theta Q(\theta | \theta^{(t)}) $$

为什么 ELBO 最大化等价于 Q 函数最大化?因为:

$$ \mathcal{L}(Q, \theta) = \sum_Z P(Z|X,\theta^{(t)}) \log P(X,Z|\theta) - \sum_Z P(Z|X,\theta^{(t)}) \log P(Z|X,\theta^{(t)}) $$

第二项(熵项)在 M 步中与 $\theta$ 无关,因此最大化 ELBO 等价于最大化 Q 函数。

五、收敛性证明

5.1 对数似然单调不减

EM 算法最重要的性质是:每次迭代后,观测数据的对数似然不会下降。我们来证明这一点。

定理:对于 EM 算法生成的参数序列 ${\theta^{(t)}}$,有:

$$ \log P(X|\theta^{(t+1)}) \geq \log P(X|\theta^{(t)}) $$

证明

首先,将边缘似然分解:

$$ \log P(X|\theta) = \mathcal{L}(Q, \theta) + \text{KL}(Q(Z) \| P(Z|X,\theta)) $$

其中:

$$ \mathcal{L}(Q, \theta) = \sum_Z Q(Z) \log \frac{P(X,Z|\theta)}{Q(Z)} $$

$$ \text{KL}(Q(Z) \| P(Z|X,\theta)) = -\sum_Z Q(Z) \log \frac{P(Z|X,\theta)}{Q(Z)} \geq 0 $$

这一分解可以验证:

$$ \begin{aligned} \mathcal{L}(Q, \theta) + \text{KL}(Q \| P) &= \sum_Z Q(Z) \log \frac{P(X,Z|\theta)}{Q(Z)} - \sum_Z Q(Z) \log \frac{P(Z|X,\theta)}{Q(Z)} \\ &= \sum_Z Q(Z) \log P(X,Z|\theta) - \sum_Z Q(Z) \log P(Z|X,\theta) \\ &= \sum_Z Q(Z) \log P(X|\theta) = \log P(X|\theta) \end{aligned} $$

因为 $P(X,Z|\theta) = P(Z|X,\theta) P(X|\theta)$。

现在,考虑第 $t$ 次迭代:

  • E 步设 $Q(Z) = P(Z|X, \theta^{(t)})$,则 $\text{KL}(Q | P) = 0$,所以:
    $$ \log P(X|\theta^{(t)}) = \mathcal{L}(P(Z|X, \theta^{(t)}), \theta^{(t)}) $$

  • M 步最大化 $\mathcal{L}$,所以:
    $$ \mathcal{L}(P(Z|X, \theta^{(t)}), \theta^{(t+1)}) \geq \mathcal{L}(P(Z|X, \theta^{(t)}), \theta^{(t)}) $$

  • 一般地,对任意 $Q$ 和 $\theta$:
    $$ \log P(X|\theta) \geq \mathcal{L}(Q, \theta) $$

于是:

$$ \log P(X|\theta^{(t+1)}) \geq \mathcal{L}(P(Z|X, \theta^{(t)}), \theta^{(t+1)}) \geq \mathcal{L}(P(Z|X, \theta^{(t)}), \theta^{(t)}) = \log P(X|\theta^{(t)}) $$

证毕。

5.2 收敛到局部极值

虽然 EM 保证似然单调不减,但一般只保证收敛到局部极大值鞍点,而不一定是全局最大值。收敛速度和最终结果高度依赖于初始值的选择。实践中通常采用:

  • 多次随机初始化,选择似然最高的结果
  • 使用 k-means 等方法提供较好的初始值

六、高斯混合模型(GMM)实例

高斯混合模型是 EM 算法的经典应用。让我们完整推导 GMM 的 EM 算法。

6.1 模型定义

GMM 假设数据由 $K$ 个高斯分布的混合生成:

$$ P(x|\theta) = \sum_{k=1}^{K} \pi_k \mathcal{N}(x | \mu_k, \Sigma_k) $$

其中:

  • $\pi_k$ 是第 $k$ 个组分的混合系数,满足 $\pi_k \geq 0, \sum_k \pi_k = 1$
  • $\mu_k$ 和 $\Sigma_k$ 是第 $k$ 个高斯分布的均值和协方差矩阵
  • 隐变量 $z_i \in {1, 2, \dots, K}$ 表示样本 $i$ 来自哪个组分

完整模型参数:$\theta = {\pi_k, \mu_k, \Sigma_k}_{k=1}^{K}$

6.2 完全数据似然

观测数据 $x_i$ 和隐变量 $z_i$ 的联合分布为:

$$ P(x_i, z_i = k|\theta) = \pi_k \mathcal{N}(x_i | \mu_k, \Sigma_k) $$

完全数据对数似然为:

$$ \log P(X, Z|\theta) = \sum_{i=1}^{N} \sum_{k=1}^{K} \mathbb{I}(z_i = k) [\log \pi_k + \log \mathcal{N}(x_i | \mu_k, \Sigma_k)] $$

6.3 E 步

计算每个样本属于各组分的后验概率(称为 responsibility):

$$ \begin{aligned} \gamma_{ik} &= P(z_i = k | x_i, \theta^{(t)}) \\ &= \frac{P(x_i | z_i = k, \theta^{(t)}) P(z_i = k | \theta^{(t)})}{\sum_{j=1}^{K} P(x_i | z_i = j, \theta^{(t)}) P(z_i = j | \theta^{(t)})} \\ &= \frac{\pi_k^{(t)} \mathcal{N}(x_i | \mu_k^{(t)}, \Sigma_k^{(t)})}{\sum_{j=1}^{K} \pi_j^{(t)} \mathcal{N}(x_i | \mu_j^{(t)}, \Sigma_j^{(t)})} \end{aligned} $$

$\gamma_{ik}$ 可以理解为样本 $i$ 属于组分 $k$ 的”软分配”概率。

6.4 M 步

最大化 Q 函数得到参数更新公式:

混合系数

$$ \pi_k^{(t+1)} = \frac{1}{N} \sum_{i=1}^{N} \gamma_{ik} $$

即组分 $k$ 的平均 responsibility。

均值

$$ \mu_k^{(t+1)} = \frac{\sum_{i=1}^{N} \gamma_{ik} x_i}{\sum_{i=1}^{N} \gamma_{ik}} $$

即所有样本的加权平均,权重为 responsbility。

协方差

$$ \Sigma_k^{(t+1)} = \frac{\sum_{i=1}^{N} \gamma_{ik} (x_i - \mu_k^{(t+1)})(x_i - \mu_k^{(t+1)})^T}{\sum_{i=1}^{N} \gamma_{ik}} $$

6.5 GMM EM 算法伪代码

Algorithm: EM for Gaussian Mixture Model
------------------------------------------------------
Input: data X = {x_1, ..., x_N}, number of components K
Output: parameters {π_k, μ_k, Σ_k} for k = 1..K

1. Initialize π_k, μ_k, Σ_k (e.g., use k-means result)

2. Repeat until convergence:
// E-step: compute responsibilities
For each i = 1..N, k = 1..K:
γ_ik = π_k * N(x_i | μ_k, Σ_k) / Σ_j π_j * N(x_i | μ_j, Σ_j)
// use log-sum-exp trick for numerical stability

// M-step: update parameters
For each k = 1..K:
N_k = Σ_i γ_ik
π_k = N_k / N
μ_k = (1/N_k) * Σ_i γ_ik * x_i
Σ_k = (1/N_k) * Σ_i γ_ik * (x_i - μ_k)(x_i - μ_k)^T

3. Check convergence: |log L(t+1) - log L(t)| < ε

七、抛硬币数值例子

一个更直观的例子:假设两种硬币 A 和 B,各自正面概率为 $\theta_A$ 和 $\theta_B$。每次实验先随机选一枚硬币(等概率),然后抛 10 次。但我们只记录每次实验的结果(正反面序列),不记录用了哪枚硬币。

7.1 数据生成

假设真实参数 $\theta_A = 0.6, \theta_B = 0.5$,进行了 5 次实验。观测数据为正面次数:

实验 正面次数 总次数 实际硬币
1 5 10 B
2 9 10 A
3 8 10 A
4 4 10 B
5 7 10 A

EM 算法不知道实际硬币,仅凭正面次数来推断 $\theta_A$ 和 $\theta_B$。

7.2 迭代过程

初始化:设 $\theta_A^{(0)} = 0.6, \theta_B^{(0)} = 0.5$

第 1 轮 E 步:计算每组数据来自 A 或 B 的概率(给定当前参数估计)

对于实验 1(5 次正面):
$$ \begin{aligned} P(5H,5T|A) &= \binom{10}{5} (0.6)^5 (0.4)^5 \propto (0.6)^5(0.4)^5 = 7.96 \times 10^{-4} \\ P(5H,5T|B) &= \binom{10}{5} (0.5)^5 (0.5)^5 \propto (0.5)^5(0.5)^5 = 9.76 \times 10^{-4} \end{aligned} $$

$$ P(A|5H) \propto 7.96 \times 10^{-4} \times 0.5 = 3.98 \times 10^{-4} $$
$$ P(B|5H) \propto 9.76 \times 10^{-4} \times 0.5 = 4.88 \times 10^{-4} $$

归一化后:
$$ P(A|5H) = \frac{3.98}{3.98 + 4.88} = 0.45 $$
$$ P(B|5H) = \frac{4.88}{3.98 + 4.88} = 0.55 $$

对于实验 2(9 次正面):
$$ \begin{aligned} P(9H,1T|A) &\propto (0.6)^9(0.4)^1 = 1.03 \times 10^{-3} \\ P(9H,1T|B) &\propto (0.5)^9(0.5)^1 = 9.76 \times 10^{-4} \end{aligned} $$

$$ P(A|9H) = \frac{1.03 \times 0.5}{1.03 \times 0.5 + 0.976 \times 0.5} = 0.51 $$
$$ P(B|9H) = 0.49 $$

类似地计算所有实验的后验概率,得到完整的 responsibility。

第 1 轮 M 步:用加权平均更新参数

设计算得到的 $\gamma_{iA}$ 和 $\gamma_{iB}$ 为:

实验 正面 $\gamma_{iA}$ $\gamma_{iB}$
1 5 0.45 0.55
2 9 0.51 0.49
3 8 0.52 0.48
4 4 0.46 0.54
5 7 0.50 0.50

更新 A:
$$ \theta_A^{(1)} = \frac{0.45 \times 5 + 0.51 \times 9 + 0.52 \times 8 + 0.46 \times 4 + 0.50 \times 7}{0.45 \times 10 + 0.51 \times 10 + 0.52 \times 10 + 0.46 \times 10 + 0.50 \times 10} $$
$$ = \frac{16.5}{24.4} \approx 0.68 $$

更新 B:
$$ \theta_B^{(1)} = \frac{0.55 \times 5 + 0.49 \times 9 + 0.48 \times 8 + 0.54 \times 4 + 0.50 \times 7}{0.55 \times 10 + 0.49 \times 10 + 0.48 \times 10 + 0.54 \times 10 + 0.50 \times 10} $$
$$ = \frac{16.1}{25.6} \approx 0.63 $$

后续迭代

重复 E 步和 M 步。通常 5-10 轮后收敛到接近真实值 $\theta_A \approx 0.8, \theta_B \approx 0.45$。

(注意:由于样本量小,EM 估计可能稍偏离真实值。)

7.3 数值稳定性:Log-Sum-Exp 技巧

在实际计算中,如上例的 $P(X|A)$ 使用 $(0.6)^5(0.4)^5$,当正面次数很大时,连乘会产生极小的浮点数导致下溢。因此实际实现中必须使用对数域计算:

# 计算 log responsibility
log_gamma_ik = log(π_k) + sum_d [x_id * log(p_kd) + (1-x_id) * log(1-p_kd)]

# log-sum-exp 归一化
log_sum = logsumexp(log_gamma_i1, log_gamma_i2, ..., log_gamma_iK)
log_gamma_ik = log_gamma_ik - log_sum

# 转回原始域
gamma_ik = exp(log_gamma_ik)

其中 logsumexp(a_1, ..., a_K) = a_max + log(Σ_k exp(a_k - a_max)) 是数值稳定的实现。

八、EM 算法的变体与扩展

8.1 广义 EM 算法(GEM)

在某些模型中,M 步的精确最大化难以实现。广义 EM(Generalized EM)放宽要求:M 步只需要增加 Q 函数而不是最大化它:

$$ Q(\theta^{(t+1)} | \theta^{(t)}) \geq Q(\theta^{(t)} | \theta^{(t)}) $$

GEM 仍能保证似然单调不减,但收敛速度可能较慢。

8.2 随机 EM(Stochastic EM)

当隐变量后验分布比较复杂时,精确计算期望可能困难。随机 EM 在 E 步使用 Monte Carlo 采样来近似期望:

  1. 从 $P(Z|X, \theta^{(t)})$ 中采样 $Z^{(s)}$
  2. 用样本均值近似期望:$\hat{Q}(\theta|\theta^{(t)}) = \frac{1}{S} \sum_s \log P(X, Z^{(s)} | \theta)$

Monte Carlo EM(MCEM)是这类方法的总称。

8.3 变分 EM(Variational EM)

当 $P(Z|X, \theta)$ 难以精确计算时,使用一个参数化的变分分布 $Q(Z|\phi)$ 来近似。E 步变为优化 $\phi$ 使 KL 散度最小:

$$ \phi^{(t+1)} = \arg\min_\phi \text{KL}(Q(Z|\phi) \| P(Z|X, \theta^{(t)})) $$

这在 LDA(潜在狄利克雷分配)等复杂模型中被广泛使用。

8.4 在线 EM(Online EM)

在大规模数据中,每次遍历全部数据成本太高。在线 EM 每看到一个(或一小批)样本就更新参数:

  • 步进式 EM(Incremental EM):对新增样本计算 responsibility,逐步更新参数
  • 随机梯度 EM:对参数做随机梯度上升

8.5 EM 与 k-means 的关系

k-means 聚类可以看作 GMM EM 算法的一个特例(极限情况)。当 GMM 中各组分协方差趋近于 $\sigma^2 I$ 且 $\sigma^2 \to 0$ 时:

  • E 步的软分配 $\gamma_{ik}$ 退化为硬分配(每个样本只属于最近的簇中心)
  • M 步的加权平均退化为簇内样本的算术平均

因此 k-means 等价于协方差极小的 GMM 的硬 EM(hard EM)。

九、EM 的局限性与实践建议

9.1 主要局限

问题 说明
局部最优 EM 只保证收敛到局部极值,结果依赖初始化
收敛速度 在高维或接近收敛时可能很慢(线性收敛)
模型选择 需要预先指定组分数量 K,不能自动确定
奇异解 GMM 中,当某组分只有一个样本时,协方差矩阵可能奇异
隐变量推断 对复杂模型,精确的 E 步计算可能不可行

9.2 实践建议

  1. 多次初始化:运行多次,从不同起点出发,选择似然最高的解
  2. 正则化:对协方差矩阵添加小的对角项 $\epsilon I$ 防止奇异
  3. 模型选择:使用 BIC(Bayesian Information Criterion)或 AIC 选择 K:
    • $\text{BIC} = -2\log L + d \log N$
    • $\text{AIC} = -2\log L + 2d$
      其中 $d$ 是参数个数,$N$ 是样本数
  4. 收敛判据:监控对数似然的相对变化:$|\log L_{t+1} - \log L_t| / |\log L_t| < 10^{-6}$
  5. k-means 初始化:用 k-means 结果作为 GMM 的初始均值和组分分配

十、其他应用场景

EM 算法的应用极其广泛,包括但不限于:

应用领域 具体示例
聚类分析 GMM 聚类
主题模型 LDA、pLSA
隐马尔可夫模型 HMM 参数学习(Baum-Welch算法)
缺失值填补 多重插补
因子分析 FA 和概率 PCA 的参数估计
图像分割 基于 GMM 的前景/背景分割
遗传学 基因型推断(haplotype phasing)
推荐系统 矩阵分解中的缺失值处理
医学影像 PET/CT 图像重建
信号处理 噪声环境下的信号参数估计

十一、面试高频问题

Q1:EM 算法和梯度上升有什么区别?什么时候该用 EM?

EM 利用问题的概率结构(隐变量)来构造一个每次都易优化的 Q 函数,而不是直接对复杂的边缘似然做优化。当存在隐变量且完全数据似然易于最大化时,EM 通常比梯度方法更稳定且无需调学习率。但当模型参数多且完全数据似然不易最大化(如深度生成模型),VAE 等基于梯度的方法更合适。

Q2:E 步和 M 步的直观理解是什么?

E 步是”猜测”隐变量——在给定当前参数估计下,推断每个隐变量取各个值的概率。可以理解为”在有噪声的数据中猜测哪些数据点属于哪个组分”。

M 步是”更新参数”——假设 E 步的”猜测”是正确的,用这些加权数据(软标签)去最大化似然来更新参数。可以理解为”根据猜测重新拟合各组分的参数”。

两者交替:E 步用参数去猜隐变量,M 步用猜出的隐变量去更新参数,如此循环直到收敛。

Q3:EM 算法一定收敛到全局最优吗?

不。EM 算法只能保证收敛到局部极值或鞍点。对于非凸目标函数(如 GMM),EM 的结果高度依赖于初始化。实践中通过多次随机初始化和选择最佳似然来缓解。此外,EM 的收敛速度在接近极值时是线性的(即 $|\theta^{(t+1)} - \theta^| \approx \rho |\theta^{(t)} - \theta^|$),可能比牛顿法等二阶方法慢。

Q4:如果 GMM 中某组分退化为单个样本(协方差矩阵奇异),怎么办?

这种情况称为”退化解”或”collapsing”。解决方案包括:

  1. 正则化:约束协方差矩阵,如 $\Sigma_k \succeq \epsilon I$,或对特征值设下限
  2. 先验:使用共轭先验(Normal-Inverse-Wishart)进行 MAP 估计而非 MLE
  3. 重新初始化:检测到奇异时,将该组分重新随机初始化
  4. 增大样本量:确保组分样本数不会太小

Q5:请推导 EM 算法收敛性的关键不等式。

核心是使用 KL 散度和 Jensen 不等式。考虑分解:

$$ \log P(X|\theta) = \underbrace{\mathbb{E}_Q[\log P(X,Z|\theta)] + H(Q)}_{\mathcal{L}(Q,\theta)} + \underbrace{\text{KL}(Q \| P(Z|X,\theta))}_{\geq 0} $$

E 步设置 $Q = P(Z|X,\theta^{(t)})$ 使 KL 项为零,M 步最大化 $\mathcal{L}$。由于 $\mathcal{L}$ 是 $\log P(X|\theta)$ 的下界(在 $\theta^{(t)}$ 处紧贴),所以提升下界必然提升似然。

或者使用链式不等式:

$$ \log P(X|\theta^{(t+1)}) \geq Q(\theta^{(t+1)}|\theta^{(t)}) + H(\theta^{(t)}) \geq Q(\theta^{(t)}|\theta^{(t)}) + H(\theta^{(t)}) = \log P(X|\theta^{(t)}) $$

其中第二步是因为 M 步保证了 $Q(\theta^{(t+1)}|\theta^{(t)}) \geq Q(\theta^{(t)}|\theta^{(t)})$。

十二、总结

EM 算法是统计学习中最优美的算法之一。它将一个包含隐变量的难以直接优化的问题,分解为两个更简单的子问题交替求解。其核心数学工具是 Jensen 不等式和 KL 散度,保证了似然函数的单调收敛。

EM 算法已经被广泛应用于聚类(GMM)、主题模型(LDA)、隐马尔可夫模型(Baum-Welch)、缺失数据填补等多个领域。深入理解 EM 算法不仅是掌握统计学习的关键一步,也是后续学习变分推断、MCMC 等高级方法的基础。

从算法设计的角度看,EM 教给我们的不仅是数学推导,更是一种哲学:面对一个困难的问题,不要硬攻,而是构建一个更简单的辅助问题,通过迭代地解决辅助问题来逼近原问题的解。这也许是 EM 算法留给我们的最深刻的启示。

参考文献

  1. Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society: Series B, 39(1), 1-38.
  2. Bishop, C. M. (2006). Pattern Recognition and Machine Learning. Chapter 9.
  3. 李航. (2019). 《统计学习方法》(第2版). 第9章.
  4. Murphy, K. P. (2012). Machine Learning: A Probabilistic Perspective. Chapter 11.
文章作者: Leo·Cheung
文章链接: http://tufusi.com/2022/01/10/%E3%80%90%E7%BB%9F%E8%AE%A1%E5%AD%A6%E4%B9%A0%E6%96%B9%E6%B3%95%E6%AD%BB%E7%A3%95%E7%B3%BB%E5%88%97%E3%80%91EM%E7%AE%97%E6%B3%95/
版权声明: 本博客所有文章除特别声明外,均采用 CC BY-NC-SA 4.0 许可协议。转载请注明来自 ONE·PIECE
打赏
  • 微信
  • 支付宝

评论