机器学习进阶（六）EM算法

前提知识

1. Jensen不等式
   ${\theta }_1+{\theta }_2+,\dots ,+{\theta }_n=1,f(x)$ 为凸函数，则有
   $$f({\theta }_1{x}_1+\cdots +{\theta }_n{x}_n)\le {\theta }_{1}f(x_1)+\cdots +{\theta }_{n}f(x_n)$$
   $$\begin{gathered} ifp(x)\ge 0onS\subset domf,{\int }_{S}p(x)=1, \\ thenf({\int }_{S}p(x)xdx)\le {\int }_{S}p(x)f(x)dx \end{gathered}$$
2. EM算法要解决的问题
   模型存在一个不能被观察到的潜变量latent variable,但是该变量会影响其他变量的取值。

GMM高斯混合模型

背景：随机变量$X$是有$k$个高斯分布混合而成，取各个高斯分布的概率为${\pi }_1,{\pi }_2...{\pi }_k$。
第$i$个高斯分布的均值为${\mu }_i$，方差为${\Sigma }_i$。
设观测到随机变量$X$的一系列样本$x_1,x_2,\dots ,x_n$。若得到的观察数据有未观察到的隐含数据${\pi }_1,{\pi }_2...{\pi }_k$ ，即上文中每个样本属于哪个分布是未知的则极大似然为$$L_{\pi ,\mu ,\Sigma }(x)=\prod _{i=1}^{N}p(x_i)=\prod _{i=1}^N\sum _{k=1}^{K}p(x_i,{\pi }_k)=\prod _{i=1}^N\left(\sum _{k=1}^K{\pi }_{i}N(x_i|{\mu }_i,{\Sigma }_i)\right)$$
第二个等号是根据$x_i$的边缘概率为${\sum }_{k=1}^{K}p(x_i,{\pi }_k)$计算得来，第三个等号是条件分布$p(x_i,{\pi }_k)={\pi }_{i}p(x_i,|{\pi }_k)$。

从而对数似然则为$$l_{\pi ,\mu ,\Sigma }(x)=\sum _{i=1}^N\log \sum _{k=1}^K{\pi }_{i}N(x_i|{\mu }_i,{\Sigma }_i)$$

为了解决这个问题，分成两步：

1. 估算数据来自的组份E-step
   估计数据由每个组份生成的概率：对于每个样本$x_i$，它由第$k$个组份生成的概率为
   $$\gamma (i,k)=\frac{{\pi }_{k}N(x_i|{\mu }_k,{\Sigma }_k)}{{\pi }_{j}N(x_i|{\mu }_j,{\Sigma }_j)}$$
   上式中的$\mu$和$\Sigma$同样也是待估计的值。
   使用采样迭代法：在计算$\gamma (i,k)$时假定$\mu$和$\Sigma$已知，即需要先验给定$\mu$和$\Sigma$（对初值选择是敏感的，需要一些其他知识）。$\gamma (i,k)$亦可看成组份$k$在生成数据$x_i$时所做的贡献。
2. 估计每个组份的参数M-step
   对于样本点$x_i$而言，可看成是每个组份$k$生成 { γ ( i , k ) x i } \{\gamma(i,k)x_i\} {γ(i,k)xi​}共同组成了$x_i$。其中，组份$k$是一个标准的高斯分布。
   $$\begin{array}{rl}{N}_k& =\sum _{i=1}^N\gamma (i,k)\\ {\mu }_k& =\frac{1}{N_k}\sum _{i=1}^N\gamma (i,k)x_i\\ {\Sigma }_k& =\frac{1}{N_k}\sum _{i=1}^N\gamma (i,k)\left(x_i-{\mu }_k\right){\left(x_i-{\mu }_k\right)}^T\\ {\pi }_k& =\frac{N_k}{N}=\frac{1}{N}\sum _{i=1}^N\gamma (i,k)\end{array}$$

EM算法

令$X$成为观察到的变量，其密度函数为$p_{\theta }(x)$。
令$Z$成为缺失变量或潜在变量。引入其分布函数为$Q(z)$，密度函数为$q(z)$。
令$p_{\theta }(X,Z)$为$(X,Z)$的真实联合分布，$p_{\theta }(Z|X)$为给定$X,Z$的条件分布。

对数似然函数即为：
$l(\theta )={\sum }_i^n\log p_{\theta }(x)={\sum }_i^n\log {\sum }_z{p}_{\theta }(x,z)$
z是隐变量，不方便直接找到参数估计。
策略： 计算$l(\theta )$下界，求该下界的最大值；重复该过程，直到收敛到局部最大值。这一迭代过程，实质是在不断地提高下界。

计算下界

令$q_i$是$Z$的某一个分布 , $q_i\ge 0$，有:
$$\begin{array}{rl}l(\theta )=& \sum _{i=1}^m\log \sum _{z_i=1}^{k}p(x_i,z_i;\theta )\\ =& \sum _{i=1}^m\log \sum _{z_i=1}^{k}p\left(x_i,z_i;\theta \right)\\ =& \sum _{i=1}^m\log \sum _{z_i=1}^k{q}_i\left(z_i\right)\frac{p\left(x_i,z_i;\theta \right)}{q_i\left(z_i\right)}\\ (Jenseninequation)\ge & \sum _{i=1}^m\sum _{z_i=1}^k{q}_i\left(z_i\right)\log \frac{p\left(x_i,z_i;\theta \right)}{q_i\left(z_i\right)}\end{array}$$
在$\frac{p\left(x_i,z_i;\theta \right)}{q_i\left(z_i\right)}=c,\forall i$时，等号可以取到。

E-step

由于$$\frac{p\left(x_i,z_i;\theta \right)}{q_i\left(z_i\right)}=c,\sum _{z_i=1}^k{q}_i\left(z_i\right)=1$$
可得(更新了$k\times m$个参数)
$$\begin{array}{rl}{q}_i\left(z_i\right)& =\frac{p\left(x_i,z_i;\theta \right)}{{\sum }_{z_i}p\left(x_i,z_i;\theta \right)}\\ & =\frac{p\left(x_i,z_i;\theta \right)}{p\left(x_i;\theta \right)}\\ & =p\left(z_i\mid x_i;\theta \right)\end{array}$$

M-step

将更新的$q_i\left(z_i\right)$带入${\sum }_{i=1}^m{\sum }_{z_i}{q}_i\left(z_i\right)\log \frac{p\left(x_i,z_i;\theta \right)}{q_i\left(z_i\right)}$在最大化该式的过程中更新参数$\theta$（若为GMM模型，$\theta$则包含${\mu }_j,{\Sigma }_j,{\pi }_j$）。
$$\begin{array}{rl}& \sum _{i=1}^m\sum _{z_i=1}^k{q}_i(z_i)\log p\left(x_i,z_i;\theta \right)\\ =& \sum _{i=1}^m\sum _{z_i=1}^k{q}_i(z_i)\log p\left(x_i|z_i;\theta \right)p(z_i;\theta )\\ =& \sum _{i=1}^m\sum _{z_i=1}^k{q}_i(z_i)\log p\left(x_i|z_i;\theta \right)+q_i(z_i)\log p(z_i;\theta )\end{array}$$

M-step 中的GMM

$$\begin{array}{rl}& \sum _{i=1}^m\sum _{z_i}{q}_i(z_i)\log p\left(x_i,z_i;\theta \right)\\ =& \sum _{i=1}^m\sum _{z_i}{q}_i(z_i)\log p\left(x_i|z_i;\theta \right)+\sum _{i=1}^m\sum _{z_i}{q}_i(z_i)\log p(z_i;\theta )\end{array}$$
注意到，在GMM模型中，加号左边的参数只有${\mu }_j,{\Sigma }_j$，加号右边的参数只有${\pi }_j$，可以分别最大化求解：
$$\begin{array}{rl}({\mu }_j,{\Sigma }_j)=& \mathrm{argmax}\sum _{i=1}^m\sum _{z_i=1}^k{q}_i(z_i)\log p\left(x_i|z_i;\theta \right)\\ {\pi }_j=& \mathrm{argmax}\sum _{i=1}^m\sum _{z_i=1}^k{q}_i(z_i)\log p(z_i;\theta )\end{array}$$
第一个式子对${\mu }_j,{\Sigma }_j$分别求偏导，令导函数为0，可得：

$$\begin{array}{rl}{\mu }_j=& \frac{{\sum }_{i=1}^m{q}_i(z_j)x_i}{{\sum }_{i=1}^m{q}_i(z_j)}\\ {\Sigma }_j=& \frac{{\sum }_{i=1}^m{q}_i(z_j)\left(x_i-{\mu }_j\right){\left(x_i-{\mu }_j\right)}^T}{{\sum }_{i=1}^m{q}_i(z_j)}\end{array}$$
第二个式子$p(z_i;\theta )$求和为0的条件（$\ge 0$的条件，由于对数的定义域即大于零，可省略）,用拉格朗日乘子法，再求偏导
$$L=\sum _{i=1}^m\sum _{z_i=1}^k{q}_i(z_i)\log {\pi }_{z_i}+\beta (\sum _{z_i=1}^k{\pi }_{z_i}-1)$$
可得
$$\begin{gathered} \frac{\partial L}{\partial {\pi }_{z_i}}=\sum _{i=1}^m\frac{q_i(z_i)}{{\pi }_{z_i}}+\beta \to \beta {\pi }_{z_i}=-\sum _{i=1}^m{q}_i(z_i)\to \beta =-m \\ {\pi }_{z_i}=\frac{1}{m}\sum _{i=1}^m{q}_i(z_i) \end{gathered}$$

步骤

交替更新，坐标上升的过程。
E-step：$$q_{new}(Z)=p_{{\theta }_{old}}(Z\mid X)$$
M-step:
$${\theta }_{new}=\underset{\theta }{\mathrm{argmax}}\sum q(Z)\log \frac{p_{\theta }(X,Z)}{q(Z)}$$

对于多个样本$(X_1,Z_1),(X_2,Z_2),\dots ,(X_n,Z_n)$,时，
$$\begin{array}{r}\sum _i^n\log p_{\theta }(X)=\sum _i^n\left(\sum _j^{K}q(Z)\log \frac{p_{\theta }(X,Z)}{q(Z)}+D_{\mathrm{k}\mathrm{l}}\left(Q(Z)\Vert P_{\theta }(Z\mid X)\right)\right)\end{array}$$

E-step：用$p_{{\theta }_{old}}(Z_i=j\mid X_i)$更新$q_{new}(Z_i=j)$
$$q_{new}(Z_i=j)=p_{{\theta }_{old}}(Z_i=j\mid X_i)=\frac{p_{{\theta }_{old}}(X_i|Z_i=j){\pi }_j}{{\sum }_j{p}_{{\theta }_{old}}(X_i|Z_i=j){\pi }_j}$$
M-step:更新${\theta }_j,{\pi }_j$
$$\begin{gathered} {\theta }_{new}=\underset{\theta }{\mathrm{argmax}}\sum _i\sum _j{q}_{new}(Z_i=j)\log \frac{p_{\theta }(X,Z_i=j)}{q_{new}(Z_i=j)} \\ =\underset{\theta }{\mathrm{argmax}}\sum _i\sum _j{q}_{new}(Z_i=j)\log \frac{p_{{\theta }_{old}}(X_i|Z_i=j){\pi }_j}{q_{new}(Z_i=j)} \end{gathered}$$

KL散度角度的EM算法

令$X$成为观察到的变量，其密度函数为$p_{\theta }(x)$。
令$Z$成为缺失变量或潜在变量。引入其分布函数为$Q(z)$，密度函数为$q(z)$。
令$p_{\theta }(X,Z)$为$(X,Z)$的真实联合分布，$p_{\theta }(Z|X)$为给定$X,Z$的条件分布。

EM算法的目标是最大化$E_X(\log p_{\theta }(X))$，引入潜变量$Z$可将式子写成
max ⁡ θ E X { log ⁡ p θ ( X ) } = max ⁡ θ { ∬ p θ 0 ( X ) q ( Z ) log ⁡ p θ ( X ) d Z d X } \max _{\theta} E_{X}\left\{\log p_{\theta}(X)\right\}=\max _{\theta}\left\{\iint p_{\theta_{0}}(X) q(Z) \log p_{\theta}(X) d Z d X\right\} θmax​EX​{logpθ​(X)}=θmax​{∬pθ0​​(X)q(Z)logpθ​(X)dZdX}
理想状态下，$\theta$ 的最佳选择${\theta }_0$，从而可得 max ⁡ θ E X { log ⁡ p θ ( X ) } = E X { log ⁡ p θ 0 ( X ) } \max _{\theta} E_{X}\left\{\log p_{\theta}(X)\right\}=E_{X}\left\{\log p_{\theta_{0}}(X)\right\} maxθ​EX​{logpθ​(X)}=EX​{logpθ0​​(X)}，等号右边可以写成
E X { log ⁡ p θ 0 ( X ) } = ∫ p θ 0 ( X ) log ⁡ p θ 0 ( X ) d X = ∬ p θ 0 ( X ) p θ 0 ( Z ∣ X ) d Z log ⁡ p θ 0 ( X ) d X = ∬ p θ 0 ( X ) p θ 0 ( Z ∣ X ) log ⁡ p θ 0 ( X ) d Z d X \begin{aligned} E_{X}\left\{\log p_{\theta_{0}}(X)\right\} &=\int p_{\theta_{0}}(X) \log p_{\theta_{0}}(X) d X \\ &=\iint p_{\theta_{0}}(X) p_{\theta_{0}}(Z \mid X) d Z \log p_{\theta_{0}}(X) d X \\ &=\iint p_{\theta_{0}}(X) p_{\theta_{0}}(Z \mid X) \log p_{\theta_{0}}(X) d Z d X \end{aligned} EX​{logpθ0​​(X)}​=∫pθ0​​(X)logpθ0​​(X)dX=∬pθ0​​(X)pθ0​​(Z∣X)dZlogpθ0​​(X)dX=∬pθ0​​(X)pθ0​​(Z∣X)logpθ0​​(X)dZdX​
可以看出$q(Z)$ 的最佳替换是 $p_{{\theta }_0}(Z\mid X)$。
对于最大化的目标$E_X(\log p_{\theta }(X))$，引入$p_{{\theta }_0}(Z\mid X)$，
E X { log ⁡ p θ ( X ) } = E X { ∫ q ( Z ∣ X ) log ⁡ p θ ( X ) d Z } ( ∫ q ( Z ∣ X ) d Z = 1 ) = E X { ∫ q ( Z ∣ X ) log ⁡ p θ ( X , Z ) p θ ( Z ∣ X ) d Z } ( p θ ( Z ∣ X ) = p θ ( X , Z ) p θ ( X ) ) = E X { ∫ q ( Z ∣ X ) log ⁡ p θ ( X , Z ) / q ( Z ∣ X ) p θ ( Z ∣ X ) / q ( Z ∣ X ) d Z } = E X { ∫ q ( Z ∣ X ) log ⁡ p θ ( X , Z ) q ( Z ∣ X ) d Z − ∫ q ( Z ∣ X ) log ⁡ p θ ( Z ∣ X ) q ( Z ∣ X ) d Z } = E X { ∫ q ( Z ∣ X ) log ⁡ p θ ( X , Z ) q ( Z ∣ X ) d Z + D k l ( Q ( Z ∣ X ) ∥ P θ ( Z ∣ X ) } \begin{aligned} & E_{X}\left\{\log p_{\theta}(X)\right\} \\ =& E_{X}\left\{\int q(Z \mid X) \log p_{\theta}(X) d Z\right\} \quad\left(\int q(Z \mid X) d Z=1\right) \\ =& E_{X}\left\{\int q(Z \mid X) \log \frac{p_{\theta}(X, Z)}{p_{\theta}(Z \mid X)} d Z\right\}\left(p_{\theta}(Z \mid X)=\frac{p_{\theta}(X, Z)}{p_{\theta}(X)}\right) \\ =& E_{X}\left\{\int q(Z \mid X) \log \frac{p_{\theta}(X, Z) / q(Z \mid X)}{p_{\theta}(Z \mid X) / q(Z \mid X)} d Z\right\} \\ =& E_{X}\left\{\int q(Z \mid X) \log \frac{p_{\theta}(X, Z)}{q(Z \mid X)} d Z-\int q(Z \mid X) \log \frac{p_{\theta}(Z \mid X)}{q(Z \mid X)} d Z\right\} \\ =& E_{X}\left\{\int q(Z \mid X) \log \frac{p_{\theta}(X, Z)}{q(Z \mid X)} d Z+D_{\mathrm{kl}}\left(Q(Z \mid X) \| P_{\theta}(Z \mid X)\right\}\right. \end{aligned} =====​EX​{logpθ​(X)}EX​{∫q(Z∣X)logpθ​(X)dZ}(∫q(Z∣X)dZ=1)EX​{∫q(Z∣X)logpθ​(Z∣X)pθ​(X,Z)​dZ}(pθ​(Z∣X)=pθ​(X)pθ​(X,Z)​)EX​{∫q(Z∣X)logpθ​(Z∣X)/q(Z∣X)pθ​(X,Z)/q(Z∣X)​dZ}EX​{∫q(Z∣X)logq(Z∣X)pθ​(X,Z)​dZ−∫q(Z∣X)logq(Z∣X)pθ​(Z∣X)​dZ}EX​{∫q(Z∣X)logq(Z∣X)pθ​(X,Z)​dZ+Dkl​(Q(Z∣X)∥Pθ​(Z∣X)}​
给定 ${\theta }_{\text{old }}$，E-Step实际做的事情是用$q(Z\mid X)$替代$P_{\theta }(Z\mid X)$，即在KL散度意义下，用$q(Z\mid X)$逼近$P_{\theta }(Z\mid X)$。
在理想状态下

EM算法是一个交替更新的过程，先更新红色部分，再更新蓝色的部分。在更新蓝色部分时，对数中的分母不含需要最大化的参数，可忽略，从而M-step从KL散度看，如下，M-step也是一个KL散度逼近的问题。
argmax ⁡ θ E X { ∫ p θ old  ( Z ∣ X ) log ⁡ p θ ( X ∣ Z ) p θ ( Z ) d Z } = argmax ⁡ θ E X { ∫ q new  ( Z ∣ X ) log ⁡ p θ ( X ∣ Z ) p θ ( Z ) p θ 0 ( X ) q new  ( Z ∣ X ) d Z } = argmax ⁡ θ { ∫ p θ 0 ( X ) q new  ( Z ∣ X ) log ⁡ p θ ( X , Z ) p θ 0 ( X ) q new  ( Z ∣ X ) d Z d X } = argmin ⁡ θ D k l { Q new  ( Z ∣ X ) P θ 0 ( X ) ∥ P θ ( X , Z ) } \begin{aligned} & \underset{\theta}{\operatorname{argmax}} E_{X}\left\{\int p_{\theta_{\text {old }}}(Z \mid X) \log p_{\theta}(X \mid Z) p_{\theta}(Z) d Z\right\} \\ =& \underset{\theta}{\operatorname{argmax}} E_{X}\left\{\int q_{\text {new }}(Z \mid X) \log \frac{p_{\theta}(X \mid Z) p_{\theta}(Z)}{p_{\theta_{0}}(X) q_{\text {new }}(Z \mid X)} d Z\right\} \\ =& \underset{\theta}{\operatorname{argmax}}\left\{\int p_{\theta_{0}}(X) q_{\text {new }}(Z \mid X) \log \frac{p_{\theta}(X, Z)}{p_{\theta_{0}}(X) q_{\text {new }}(Z \mid X)} d Z d X\right\} \\ =& \underset{\theta}{\operatorname{argmin}} D_{\mathrm{kl}}\left\{Q_{\text {new }}(Z \mid X) P_{\theta_{0}}(X) \| P_{\theta}(X, Z)\right\} \end{aligned} ===​θargmax​EX​{∫pθold ​​(Z∣X)logpθ​(X∣Z)pθ​(Z)dZ}θargmax​EX​{∫qnew ​(Z∣X)logpθ0​​(X)qnew ​(Z∣X)pθ​(X∣Z)pθ​(Z)​dZ}θargmax​{∫pθ0​​(X)qnew ​(Z∣X)logpθ0​​(X)qnew ​(Z∣X)pθ​(X,Z)​dZdX}θargmin​Dkl​{Qnew ​(Z∣X)Pθ0​​(X)∥Pθ​(X,Z)}​

总结

1. E-step是一个无约束的优化问题，
   $$Q_{\text{new }}(Z\mid X)=\underset{Q(Z\mid X)}{\mathrm{argmin}}{E}_X\left[D_{\mathrm{k}\mathrm{l}}\left(Q(Z\mid X)\Vert P_{\theta \text{ old }}(Z\mid X)\right)\right]$$
   给定${\theta }^{\text{old }}$，直接令$q(Z\mid X)=P_{\theta }(Z\mid X)$即可。

2. M-step是一个有约束的优化问题。
   $P_{{\theta }^{\text{new }}}(X\mid Z)$和$P_{{\theta }^{\text{new }}}(Z)$均为事先确定的分布簇，如$GMM$中前者为正态分布，后者为多项分布。所以这里是在优化分布，使得这两个分布的成绩和$Z\mid X)P_{{\theta }_0}(X)$的分布在KL散度意义下很接近。
   P θ new  ( X ∣ Z ) P θ new  ( Z ) = argmin ⁡ P θ ( X ∣ Z ) P θ ( Z ) ∈ F θ D k l { Q new  ( Z ∣ X ) P θ 0 ( X ) ∥ P θ ( X ∣ Z ) P θ ( Z ) } \begin{array}{l} P_{\theta^{\text {new }}}(X \mid Z) P_{\theta^{\text {new }}}(Z) \\ =\underset{P_{\theta}(X \mid Z) P_{\theta}(Z) \in \mathcal{F}_{\theta}}{\operatorname{argmin}} D_{\mathrm{kl}}\left\{Q_{\text {new }}(Z \mid X) P_{\theta_{0}}(X) \| P_{\theta}(X \mid Z) P_{\theta}(Z)\right\} \end{array} Pθnew ​(X∣Z)Pθnew ​(Z)=Pθ​(X∣Z)Pθ​(Z)∈Fθ​argmin​Dkl​{Qnew ​(Z∣X)Pθ0​​(X)∥Pθ​(X∣Z)Pθ​(Z)}​