统计学习方法 第 09 章 EM 算法及其推广

EM 算法可以用于求解存在隐变量的非监督模型生成问题：已知模型 $P(y \mid \theta)$ 和数据集 ${y_i}$，通过学习求模型参数 $\theta$。或者已知模型 $P(y, z \mid \theta)$ 和数据 ${y_i}$，${z_i}$ 未知，称为隐变量，求解模型参数 $\theta$。EM 算法有时还能估计隐变量 ${z_i}$，如隐马尔可夫模型。

EM 算法

$y_i$ 代表一个数据集中一个随机变量取值，$Y$ 代表数据集。$z_i, Z$ 与之类似。

【输入】观测变量数据 $Y={y_i}$，联合分布 $P(y_i,z_i \mid \theta)$，条件分布 $P(z_i \mid y_i, \theta)$，循环终止条件。 【输出】模型参数估计 $\theta$.

1. 初始化 $\theta = \theta^{(1)}$
2. 循环 $\text{for } i = 1 \cdots:$
   1. E 步骤，求 $Q$ 函数： \(Q(\theta, \theta^{(i)})=\sum_Z P(Z\mid Y,\theta^{(i)})\log P(Y,Z\mid \theta)\)
   2. M 步骤，求使 $Q$ 函数最大化的 $\theta$： \(\theta^{(i+1)} = \arg\max_\theta Q(\theta, \theta^{(i)})\)
   3. 检查循环终止条件。终止条件一般是 $Q$ 函数趋于稳定或者 $\theta$ 趋于稳定： \(\|\theta^{(i+1)}-\theta^{(i)}\| < \epsilon_1 \text{ or } \| Q(\theta^{(i+1)}, \theta^{(i)}) - Q(\theta^{(i)}, \theta^{(i)}) \| < \epsilon_2\)

对于具体的模型，循环中的 E 步骤和 M 步骤的最优化问题可以预先进行符号运算求解，编程实现时合并为一步迭代步骤。

EM 算法对初值敏感，可能收敛到局部最优点。

---

课本本节开头使用三硬币问题作为 EM 算法的引入，这里从 EM 算法一般形式导出三硬币问题的迭代公式。

推导承接 P156 公式 9.5 继续，问题在于从公式 9.5 如何得到公式 9.6. 以参数 $\pi$ 为例，按照 EM 算法：

\[\begin{align*} \pi &= \arg\max_\pi \sum_Z P(Z \mid Y, \theta') \log P(Y, Z \mid \theta) \end{align*}\]

式中使用 $\theta’$ 代表上一轮迭代的到的参数。主要困难在于 $\sum_Z$ 是关于隐变量序列所有可能求和（共有指数级别的可能），而非关于序列下标求和。

简单起见记 $P_Z = P(Z \mid Y, \theta’)$. 继续推导：

\[\begin{align*} \pi &= \arg\max_\pi \sum_Z P_Z \log P(Y, Z \mid \theta) \\ &= \arg\max_\pi \sum_Z P_Z \log \left[ P(Z \mid \theta) P(Y \mid Z, \theta) \right]\\ &= \arg\max_\pi \sum_Z P_Z \log \left[ P(Z \mid \pi) P(Y \mid Z, \{p, q\}) \right]\\ &= \arg\max_\pi \sum_Z P_Z \log \left[ P(Z \mid \pi) \right]\\ &= \arg\max_\pi \sum_Z P_Z \log \left[ \pi^{\sum z_i} (1-\pi)^{\sum (1-z_i)} \right]\\ &= \arg\max_\pi \sum_Z P_Z \left[ \sum z_i \log \pi + \sum (1 - z_i) \log (1 - \pi) \right]\\ &= \arg\max_\pi \sum_Z P_Z \left[ \sum_i z_i \log \pi + \sum_i (1 - z_i) \log (1 - \pi) \right] \\ &= \arg\max_\pi\left[ \log \pi \sum_Z P_Z \sum_i z_i + \log (1 - \pi) \sum_Z P_Z \sum_i (1 - z_i) \right] \\ &= \arg\max_\pi\left[ \log \pi \sum_i \sum_Z P_Z z_i + \log (1 - \pi) \sum_i \sum_Z P_Z (1 - z_i) \right] \\ \end{align*}\]

注意 $\sum_Z P_Z z_i = \sum_Z P(Z \mid Y, \theta’) z_i$ 就是 $z_i$ 在 $Y, \theta’$ 给定时的期望，也即课本中公式 9.5，因为 $z_i$ 即代表硬币 A 的到正面，也即 $y_i$ 来自硬币 $B$. 于是：

\[\begin{align*} \sum_Z P_Z z_j &= \mu_j \\ \sum_Z P_Z (1 - z_j) &= 1 - \mu_j \end{align*}\]

继续推导：

\[\begin{align*} \pi &= \arg\max_\pi\left[ \log \pi \sum_i \mu_i + \log (1 - \pi) \sum_i (1 - \mu_i) \right] \\ \end{align*}\]

这个最大化问题可以借助导数求解。最终得到：

\[\begin{align*} \pi = \frac 1 N \sum_i \mu_i \end{align*}\]

EM 算法的导出

以下内容来自自己的理解。首先解释怎么做，然后解释为什么这么做。

EM 算法可以由对数似然函数极大化问题导出：

\[\begin{align*} L(\theta) =\log P(Y \mid \theta) =\log \sum_Z P(Y,Z\mid \theta) \end{align*}\]

尝试求解 $\theta = \arg\max_\theta L(\theta)$, 但 $\log\sum$ 的形式使得求解变得困难。为此，使用 Jensen 不等式进行松弛：

\[\begin{align*} L(\theta) &= \log\left[ \sum_Z P(Y,Z \mid \theta) \right] \\ &= \log\left[ \sum_Z P(Z\mid Y, \theta^{(i)}\frac{P(Y, Z \mid \theta)}{P(Z\mid Y,\theta^{(i)})} \right] \\ &\geq \sum_Z \left[ P(Z\mid Y,\theta^{(i)}) \log \frac{P(Y, Z \mid \theta)}{P(Z\mid Y,\theta^{(i)})} \right] \\ &:= Q(\theta, \theta^{(i)})\\ \end{align*}\]

$Q$ 是 $L$ 的某个下界，因此可以尝试通过最大化 $Q$ 来逼近 $L(\theta)$ 的最大值（实际上，分号上下引入的辅助项的取值使得不等号取到等号，下文有推导）：

\[\begin{align*} \theta^{(i+1)} &= \arg\max_\theta Q(\theta, \theta^{(i)}) \\ &= \arg\max_\theta \sum_Z \left[ P(Z\mid Y,\theta^{(i)}) \log \frac{P(Y, Z \mid \theta)}{P(Z \mid Y, \theta^{(i)})} \right] \\ &= \arg\max_\theta \sum_Z P(Z\mid Y,\theta^{(i)}) \log P(Y, Z \mid \theta) \\ &= \arg\max_\theta Q(\theta, \theta^{(i+1)}) \end{align*}\]

---

接下来讨论使用 Jensen 不等式时引入的辅助函数为什么是 $P(Z\mid Y \theta^{(i)})$. Jensen 不等式：

\[\begin{align*} \phi\left[\sum_Z P(Z) f(Z)\right] &\geq \sum_Z P(Z) \phi\left[f(Z)\right] \\ \sum_Z P(Z) &= 1 \end{align*}\]

离散情境下推导最大化不等式右侧项的条件。选取恰当的 $P$ 数列最大化以下求和式：

\[\begin{align*} L &= \sum_z P_z \log \frac {F_z} {P_z} \\ &= P_1 \log \frac {F_1} {P_1} + P_2 \log \frac {F_2} {P_2} \cdots + \left(1 - \sum_{z=1}^{n-1} P_z \right) \log \frac {F_n} {1 - \sum_{z=1}^{n-1} P_z} \\ \frac{\partial L}{\partial P_1} &= \log\frac{F_1}{P_1} - \log\frac{F_n}{1 - \sum_{z=1}^{n-1} P_z} \end{align*}\]

利用 $\frac{\partial L}{\partial P_1} = 0$ 得到：

\[\begin{align*} \frac{F_1}{F_n} = \frac{P_1}{1 - \sum_{z=1}^{n-1} P_z} = \frac{P_1}{P_n} \end{align*}\]

一般的，可以选取恰当的 $P(Z)$ 最大化不等式右侧：

\[\begin{align*} P(Z) = c \cdot f(Z) = \frac{f(Z)}{\sum_Z f(Z)} \end{align*}\]

将该表达式带入原式可以发现，该条件能够使得等号成立。

于是，EM 算法的推导中，为了的到 $L$ 尽量紧的下界，引入的辅助函数应该满足：

\[\begin{align*} P &= c \cdot P(Y,Z\mid \theta) \\ \sum_Z P &= 1 \end{align*}\]

容易证明这样的函数是：$P = P(Z \mid Y, \theta)$.

于是，原始最优化目标函数被转换为

\[\begin{align*} \theta &= \arg\max_\theta \sum_Z P(Z\mid Y,\theta) \log\frac{P(Y, Z \mid \theta)}{P(Z\mid Y,\theta)} \\ &= \arg\max_\theta\sum_Z L'(\theta) \end{align*}\]

该问题仍然不易求解。EM 算法实际上是解该问题的数值解法，只对 $P(Y,Z\mid \theta)$。 项进行优化，$P(Y,Z\mid \theta)$ 中的 $\theta$ 取上一轮迭代值 $\theta^{(i)}$.

一般的，Jensen 不等式的另一个形式：

\[\begin{align*} E(f(X)) \geq f(E(X)) \end{align*}\]

等号成立的条件是 $X = E(X)$ 以概率一成立，对于性质足够好的函数，这限制了 $X$ 是常数。

EM 算法的推广：GEM 算法

GEM 不是 EM 的「强化版」，相反，其实是个「弱化版」的算法，因此比 EM 更简单。

有时 $Q$ 函数的最大化问题不易求解，可以退而求其次，分别对 $\theta$ 每个维度进行优化：

\[\theta^{(i+1)}_k = \arg\max_{\theta_k} Q(\theta, \theta^{(i)})\]

习题

9.1

Output:

9.2

【引理 9.2】若 $\tilde P(Z) = P(Z\mid Y,\theta)$，则：$F(\tilde P, \theta)=\log P(Y\mid\theta)$.

【证明】

\[\begin{align*} F(\tilde P,\theta) &= E_{\tilde P}\log P(Y,Z\mid\theta) + H({\tilde P}) \\ &= E_{\tilde P}\log\frac{P(Y,Z\mid\theta)}{ {\tilde P}_\theta(Z)} \\ &= \sum_Z P(Z\mid Y,\theta)\log\frac{P(Y,Z\mid\theta)}{P(Z\mid Y,\theta)} \\ &= \sum_Z P(Z\mid Y,\theta)\log P(Y\mid\theta) \\ &= \log P(Y\mid\theta) \end{align*}\]

9.3

以下代码尽管可以收敛，但是容易 NAN.

9.4

朴素贝叶斯模型的参数是 $\theta = {P(y=\alpha), P(x^d=\beta^d \mid y = \alpha)}$, 上标 $d$ 代表维度。下文中，下标 $n$ 代表第 $n$ 个样本。为了简便起见，记 $P(\cdot \mid \theta^{(i)}) = P’(\cdot), P(\cdot \mid \theta) = P(\cdot)$.

\[\begin{align*} Q(\theta, \theta') &= \sum_Y P'(Y \mid X) \log P(Y, X) \\ P(y=\alpha) &= \arg\max_{P(y=\alpha)}\sum_Y P'(Y \mid X) \log P(Y) P(X \mid Y) \\ &= \arg\max_{P(y=\alpha)}\sum_Y P'(Y \mid X) \log P(Y) \\ &= \arg\max_{P(y=\alpha)}\sum_Y P'(Y \mid X) \sum_n \log P(y_n) \\ &= \arg\max_{P(y=\alpha)} A \\ \frac{\partial A}{\partial P(y=\alpha_1)} &= \sum_Y P'(Y \mid X) \sum_n \left[ \frac {I(y_n = \alpha_1)} {P(y=\alpha_1)} - \frac {I(y_n = \alpha_K)} {P(y=\alpha_K)} \right] \end{align*}\]

注意，这里利用了 $P(y=\alpha_K) = 1 - \sum_{k=1}^{K-1}P(y=\alpha_k)$.

令导函数为零的到：

\[\begin{align*} \sum_n \sum_Y P'(Y \mid X) \frac{I(y_n = \alpha_1)} {P(y=\alpha_1)} &= \sum_n \sum_Y P'(Y \mid X) \frac{I(y_n = \alpha_K)} {P(y=\alpha_K)} \\ \sum_n \sum_Y P'(Y \mid X) \frac{I(y_n = \alpha_k)} {P(y=\alpha_k)} &= C_1 \\ P(y = \alpha_k) &= \frac 1 {C_1} \sum_n \sum_Y I(y_n = \alpha_k) P'(Y \mid X) \\ &= \frac 1 {C_1} \sum_n P'(y_n = \alpha_k \mid X) \\ &= \frac 1 {C_1} \sum_n \frac{P'(X \mid y_n = \alpha_k)P'(y_n = \alpha_k)}{P'(X)} \\ &= \frac 1 {C_2} \sum_n P'(x_n \mid y_n = \alpha_k)P'(y_n = \alpha_k) \end{align*}\]

$C_2$ 是一个归一化系数。类似可以得到：

\[\begin{align*} P(x^d = \beta_l \mid y = \alpha_k) = \frac 1 C \sum_{n, x_n^d = \beta_l } P'(x_n^d \mid y_n = \alpha_k) P'(y_n = \alpha_k) \end{align*}\]