混合模型与期望最大化算法

本部分内容是《模式识别与机器学习》第九章，主要回答了2个问题：
1 什么是EM期望最大化算法
2 EM算法在混合模型中的应用，如高斯混合、伯努利混合以及贝叶斯混合

目录

看目录的话，本章节分为四个小节：
1 K-means聚类：作为理解EM算法的简单例子
2 高斯混合：主要是梯度固定和后验概率的概念
3 EM算法的另一观点：EM算法能够使用使用潜在变量补全数据
4 EM算法的一般形式：将EM理解为坐标上升法

介绍

• 若我们定义了关于观测值和潜在变量的联合分布，引入潜在变量和扩展空间能够从简单分布中表达出相对复杂的边缘分布。
• 最大似然估计在该空间中就理解为期望最大化算法，但最大似然法也有很大的限制。
• 第10章将会介绍更有效的变分推断的贝叶斯处理方式

K-means聚类

失真度

给定数据集{$x_1,...,x_N$}, 我们想要将它分割为K个簇。
根据K-means的原理，我们将使用每个簇的中心点${\mu }_k$替代整个簇的(即数据集中的真实点，使簇中所有其他数据点到${\mu }_k$的距离平方和为最小的那个点)
然后对每个数据点$x_n$，使用二元指示变量(即1-K的编码策略，one-hot)$r_{nk}\in (0,1)$，即$r_{nk}$=1 且$r_{nj}=1,j̸=0$，k=1,…,K,(它只有两种类型的值)用来描述该数据点所属簇。
然后介绍失真度公式如下：
$$J=\sum _{n=1}^N\sum _{k=1}^K{r}_{nk}||x_n-{\mu }_k|{|}^2$$ (式9.1)
我们计算每个簇中的数据点和作为簇中心点距离的平方和，也称为欧氏距离。

期望最大化

• 找出{$r_{nk}$}和{${\mu }_k$}的值用于最小化失真度(式9.1)
• 迭代过程：
  • 1 求期望(确定属于哪个簇)：固定${\mu }_k$，根据$r_{nk}$最小化J
    $$r_{nk}=\left(\genfrac{}{}{0ex}{}{1，若k=argmi{n}_j||x_n-{\mu }_k|{|}^2}{0，其他}\right)$$ (式9.2)
    即确定每个数据点所属簇。
  • 2 最大化(重新找出簇中心点)：固定$r_{nk}$，根据${\mu }_k$最小化J
    $$2\sum _{n=1}^N{r}_{nk}(x_n-{\mu }_k)=0$$ (式9.3)
    要使J为最小值，需要对J的函数求导，使其导数为0的点就是最小值点。
    $${\mu }_k=\sum _n{r}_{nk}{x}_n\sum _n{r}_{nk}$$（式9.4）
    即，根据第1步分好的簇来重新选择作为中心的点。

示例

• 每个E或M步骤都会降低目标函数J的值。
• 则会收敛到局部或全局最大值。

结束

直接实现K-Means算法可能会运行缓慢，可以使用在线随机算法替代：
${\mu }_k^{new}={\mu }_k^{old}+{\eta }_n(x_n-{\mu }_k^{old})$ (式9.5)
此处的${\eta }_n为学习速率，考虑数据点增加时仍能够单调减的情况$
K-mediods(中心点)，一般的失真度量为：
$\tilde{J}={\sum }_{n=1}^N{\sum }_{k=1}^K{r}_{nk}v(x_n,{\mu }_k)$ (式9.6)
此处的$v$(., .)表示任意相异度测量，如欧氏距离、曼哈顿距离、余弦距离等。
K-Means算法可以用于图像分割和压缩的例子：

高斯混合

潜在变量

• 潜在变量与可观测变量相反，不能直接观察，但可以通过观察到的其他变量通过数学模型进行推断得到。

• 高斯混合分布可写为一般高斯分布(正态分布)的线性叠加(第二章提过)：
  $p(x)={\sum }_{k=1}^K{\pi }_k\mathcal{N}(x|{\mu }_k,{\sum }_k)$

• 引入潜在变量z

  • z是二元1-of-K编码变量
  • 联合分布：p(x, z) = p(z) p(x|z)
    

• 混合系数${\pi }_k$有：
  当使潜在变量的为1时的概率值，p($z_k$ = 1) = ${\pi }_k$，显然可得概率值的限制条件是：$0\le {\pi }_k\le 1$ 且 ${\sum }_k{\pi }_k=1$。
  因为z为one-hot型变量，该分布也可以写为：$p(z)={\prod }_k{\pi }_k^{z_k}$ (式9.10)

• 对特定潜在变量$z_k$，x的条件分布为：
  $p(x|z_k=1)=\mathcal{N}(x|{\mu }_k,{\sum }_k)$同理也可以写为：$p(x|z)={\prod }_k\mathcal{N}(x|{\mu }_k,{\sum }_k{)}^z{}_k$(式9.11)

• x的边缘分布可以通过对联合分布对z的所有可能状态求和得出(结合式9.10和9.11)：
  $p(x)={\sum }_{z}p(x,z)={\sum }_{z}p(z)p(x|z)={\sum }_k{\pi }_k\mathcal{N}(x|{\mu }_k,{\sum }_k)$
  同时，x的边缘分布符合式9.7的高斯混合形式。

• 通过利用了联合概率p(x, z)，使得原式9.7得到了大大的简化。

• 利用关于k的后验概率来生成观测x(贝叶斯原理)：
  $$p(z_k=1|x)=\gamma (z_k)\equiv p(z_k=1|x)=\frac{p(z_k=1)p(x|z_k=1)}{{\sum }_{k}p(z_k=1)p(x|z_k=1)}$$
  $$=\frac{{\pi }_k\mathcal{N}(x|{\mu }_k,{\sum }_k)}{{\sum }_k{\pi }_k\mathcal{N}(x|{\mu }_k,{\sum }_k}$$ (式9.13)
  此处${\pi }_k$表示$z_k=1$的先验概率。

• 使用祖先取样生成随机样本：
  首先由p(z)生成$\hat{z}$)，随后从p(x | $\hat{z}$)(即高斯混合模型中生成x的值，具体见第11章。
  
  如图a，高斯部分负责生成x。
  图b，边缘分布p(x)的相似样本通过从联合分布中采样并忽略z值得到。
  图c，可以表示与数据点$x_n$联系的先验概率值$\gamma (z_{nk})$，并绘出对应点使用红绿蓝对应$\gamma (z_{nk})$的k=1, 2, 3时。如$\gamma (z_{n1})=1$显示为红色，$\gamma (z_{n2})=\gamma (z_{n3})=0.5$分别显示为绿色和蓝色。与图a相比，数据点的标记使用生成的真实component密度。

最大似然法

• 最大似然法也称为极大似然估计，是一种点估计法。基本思想是当从模型总体随机抽取n组样本观测值后，最合理的参数估计量应该是使从模型中抽取的这n组样本的观测值的概率最大。

• log似然：
  假定观测数据集{$x_1,...,x_N$}，然后使用高斯混合对该数据建模。将该数据表示为一个N * D矩阵X，即第n行表示为$x_n^T$。相似的，对应潜在变量为N * K的矩阵Z和行$z_n^T$。
  若假定数据点与分布相互独立，则对该独立同分布数据的高斯混合模型可以表示为下图。

$lnp(X|\pi ,\mu ,\sum )={\sum }_{n=1}^{N}ln\left\{\begin{array}{c}{\sum }_{k=1}^K{\pi }_k\mathcal{N}(x|{\mu }_k,{\sum }_k\end{array}\right\}$(式9.14)
其中一组分为协方差矩阵${\sum }_k={\sigma }_k^{2}I，I$为单位矩阵。假定混合模型的第j个组分有均值${\mu }_j$，它实际上等于其中一个数据点，即对一些n值有${\mu }_j=x_n$，使得该组分坍塌于一个点。$\pi$为$z_n$组分的混合系数。

• 奇异性
  
  如上图所示，当混合中有至少2个组分时，其中一个含有限方差，因此可以分配有限概率给所有数据点，但当其他组分快速上升到一个特定数据点时会给log似然造成一个不断增加的相加值，则该点就称为高斯混合的奇点。
  此外，对高斯混合使用最大似然法可能会导致较差的解，若使用贝叶斯方法就不会出现该问题。(第11章)

• 唯一性
  在k-组分混合中对一个机器学习的解，有$K!$个等价解。对应于有K!种方式去分配K个数据集的参数给K的组分。即，对任意参数空间中的点(非退化点)，将至多有K!-1个其他点也能准确达到相同分布，称之为识别能力，当解释从模型中发现的参数时是一个重要的问题。将在第十二章详细讨论。

EM用于高斯混合

• Dempster et al., 1977正式引入期望最大化算法

• 最大化log似然：
  使$lnp(X|\pi ,\mu ,\sum )$的导数相关参数为0：
  根据式9.14:$lnp(X|\pi ,\mu ,\sum )={\sum }_{n=1}^{N}ln\left\{\begin{array}{c}{\sum }_{k=1}^K{\pi }_k\mathcal{N}(x|{\mu }_k,{\sum }_k)\end{array}\right\}$

• 对${\mu }_k$求导有：
  $0=-{\sum }_{n=1}^N\frac{{\pi }_k\mathcal{N}(x|{\mu }_k,{\sum }_k}{\begin{array}{c}\sum _k{\pi }_k\mathcal{N}(x|{\mu }_k,\sum _k\\ \gamma (z_k)\end{array}}{\sum }_k^{-1}(x_n-{\mu }_k)$
  通过乘${\sum }_k^{-1}$(非奇异矩阵)得，${\mu }_k=\frac{1}{{\sum }_n\gamma (z_k)}{\sum }_n\gamma (z_k)x_n$ (式9.17)
  第k个高斯组分的均值${\mu }_k$通过数据集中所有数据点的权重均值得到，而数据点$x_n$的权重通过先验概率$\gamma (z_{nk})$给定，即组分k负责生成$x_n$。

• 对${\sum }_k$求导有：
  ${\sum }_k=\frac{1}{{\sum }_n\gamma (z_k)}{\sum }_n\gamma (z_k)(x_n-{\mu }_k)(x_n-{\mu }_k{)}^T$ (式9.19)

• 对混合系数${\pi }_k$有：

  • 加入限制条件${\sum }_k{\pi }_k=1$
  • 使用拉格朗日乘数法并最大化如下式子：
    $lnp(X|\pi ,\mu ,\sum )+\lambda ({\sum }_k{\pi }_k-1)$
    可以得到，$0={\sum }_n\frac{\mathcal{N}(x|{\mu }_k,{\sum }_k)}{{\sum }_k{\pi }_k\mathcal{N}(x|{\mu }_k,{\sum }_k)}+\lambda$
    ${\pi }_k=\frac{{\sum }_n\gamma (z_k)}{\mathcal{N}}$
    此处，令$\lambda =-N$。即第k个组分的混合系数有每个组分用于解释数据点的平均后验概率给定。

EM用于高斯混合的例子

• 没有封闭的解：因为先验概率$\gamma (z_k)依赖于参数$（见式9.13）
• 但上述等式也暗示了找出最大似然的简化迭代策略：
  首先选择均值、协方差和混合系数的初始值，然后交替进行计算当前$\gamma (z_k)$以及更新参数{${\mu }_k,{\sum }_k,{\pi }_k$}
  
  如上图，
  • a中以绿色展示了数据点，和混合模型的初始配置并使用蓝色和红色的圈表示标准差轮廓。
  • b展示了初始化E步骤的结果，蓝色部分表示等于由蓝色组分生成的后验概率，红色同理。此外，拥有显著概率的点属于2个簇显示为紫色。
  • c展示了第一次M步骤后的结果，蓝色高斯分布的均值移动到数据集的均值，概率的权重来自属于蓝色簇的每个数据点，即移动到蓝色区域的质量中心。相似的，蓝色高斯分布的协方差设置为蓝色区域的协方差。红色同理。
  • d, e, f分别展示了完成2，5，20个EM完整循环后的效果，其中f更接近于收敛。
• 但相比于K-Means，高斯混合需要进行更多的迭代才能收敛，以及每次循环都得进行更多的计算。
• 因此一般先使用K-Means初始化参数，协方差矩阵可以通过从K-Means算法发现的集群中的样本协方差得到，然后混合系数可以由对应集群的数据点分数得到，均值同理。

EM用于高斯混合总结

1. 初始化{${\mu }_k,{\sum }_k,{\pi }_k$}并计算log似然
2. E步骤计算先验概率$\gamma (z_k)$
3. M步骤使用当前先验概率重新计算参数：
   ${\mu }_k^{new}=\frac{1}{{\sum }_n\gamma (z_k)}{\sum }_n\gamma (z_k)x_n$
   ${\sum }_k^{new}=\frac{1}{{\sum }_n\gamma (z_k)}{\sum }_n\gamma (z_k)(x_n-{\mu }_k)(x_n-{\mu }_k{)}^T$
   ${\pi }_k^{new}=\frac{{\sum }_n\gamma (z_k)}{N}$
4. 计算Log似然$lnp(X|\pi ,\mu ,\sum )并检查收敛(返回第二步)$

EM的另一观点

潜在变量

• 令X为观测数据，Z为潜在变量， $\theta$为所有参数
• 目标：最大化观测值的边缘log似然
  $lnp(X|\theta )=ln\left\{\begin{array}{c}{\sum }_{z}p(X,Z|\theta )\end{array}\right\}$ (式9.29)
• 但由于log求和导致优化问题，因为联合分布$p(X,Z|\theta )$属于指数族，则边缘分布$p(X|\theta )$并不作为求和的结果，而因为和的存在是算法不能直接作用于联合分布，导致了最大似然解的表达更加复杂。
• 假设对完整数据的直接最大化：
  $lnp(X,Z|\theta )$，但一般通过观测只能知道X的值，而潜在的Z只能通过后验分布$p(Z|X,\theta )$来得到，考虑潜在变量的先验分布期望值也对应于EM算法中的E步骤。接下来将进行M步骤，最大化期望。

算法

• 初始化：选择初始参数集${\theta }^{old}$
• E步骤：使用当前参数${\theta }^{lod}$来计算$p(Z|X,{\theta }^{old})$来找出对于一般化参数$\theta$的完整数据的log似然的期望：
  $\mathcal{Q}(\theta ,{\theta }^{old})={\sum }_{Z}p(Z|X,{\theta }^{old})lnp(X,Z|\theta )$(式9.30)
• M步骤：通过最大化式9.30来确定${\theta }^{new}$
  ${\theta }^{new}=argma{x}_{\theta }\mathcal{Q}(\theta ,{\theta }^{old}$(式9.31)
• 然后检查收敛：收敛则停止，或KaTeX parse error: Expected 'EOF', got '\the' at position 17: …theta^{old} <- \̲t̲h̲e̲^{new}并进入E步骤。

高斯混合回顾

• 对混合的每个x分配潜在变量$z_k$
  
• 完整数据(log-)似然(9.36)和期望(9.40)
  似然函数：
  $p(x,z|\theta )={\prod }_{k=1}^K{\pi }_k^{z_k}\mathcal{N}(x;{\mu }_k,{\sum }_k{)}^{z_k}$
  取对数得：
  $lnp(x,z|\theta )={\sum }_{k=1}^K{z}_{k}ln{\pi }_k+ln\mathcal{N}(x;{\mu }_k,{\sum }_k)$
  对比于非完整数据的似然函数对k的求和，log型函数是可以互换的，即可以直接作用于高斯分布（它本身就是指数族的成员），这使得对最大似然法有了更容易的解。
  求期望：
  $\mathcal{Q}(\theta )=I{E}_z[lnp(x,z|\theta )]={\sum }_{k=1}^K\gamma (z_k)ln{\pi }_k+ln\mathcal{N}(x;{\mu }_k,{\sum }_k)$

上星期我们介绍了K-means算法、EM一般过程及其在高斯混合中的应用，接下来继续介绍几个EM的应用

EM算法示例

伯努利混合

• 至今我们介绍高斯分布描述的连续变量，接下来介绍不同背景的EM算法：伯努利混合描述的离散二元变量，其分布如下(该模型也称为潜在类别分析，它也为考虑离散变量的隐马可夫模型奠定了基础)：
  KaTeX parse error: Expected 'EOF', got '}' at position 50: …{x_i}(1 - \mu_i}̲^{(1 - x_i)}(式9.44)
  即为D个二元变量(矢量)，i= 1到D，每个都由参数mu(矢量)控制，独立变量xi独立于mu，显然期望为mu(矢量)， 协方差为对角mu i(1 - mu i)。
• 伯努利混合还可以可以对变量相关性进行建模，也称为潜在类别分析。
• 与高斯混合相反，在本模型中没有奇点(即会使似然函数变得无穷大的点)，因为$p ( x_n | \mu_k) $大于0小于1使得似然函数有上界，存在的奇点变为0。而且模型是线性且数据完整的log似然。
• 用简便的EM算法去寻找ML参数
  • E步骤：计算混合组分k对于数据点xn的先验概率$\gamma (z_{n}k)\propto {\pi }_{k}p(x_n|{\mu }_k)$
  • M步骤：更新参数${\pi }_k=N^{-1}{\sum }_n\gamma (z_{nk})$和${\mu }_k=(N{\pi }_k{)}^{-1}{\sum }_n\gamma (z_{nk})x_n$
  • 示例：如下图，将该模型用于手写数字。数字图片已经被转换为二元矢量(将超过0.5的所有值设置为1，其他设置为0)。然后去拟合个数N=600的数据集，使用混合组分个数K=3的伯努利混合进行10次EM算法迭代来比较数字2， 3， 4。混合系数被初始化为${\pi }_k=1/K$，参数${\mu }_{k_{kj}}$设置为从(0.25, 0.75)中随机选取的值，然后通过标准化来满足限制条件${\sum }_j{\mu }_{kj}=1$。我们可以看到3个伯努利分布的混合已经能够找出数据集中的3个簇，这3个簇也对应于3个不同的数字。
    中间4个展示的数字是使用阈值0.5将灰度图转换为二元的效果，最右侧3个数字展示了混合模型中对于3个组分不同mu k(作为显示像素值)的效果，作为对比左侧的一个数字为使用单个多元伯努利分布(同样使用最大似然法)的去拟合相同数据集的效果，可见它只是简单的重复每个像素值。
    

贝叶斯线性回归

• 作为介绍EM应用的第三个例子，我们回顾第三章贝叶斯线性回顾的证据近似，我们是先通过计算证据，然后再设置证据表达式的导数为0，才能获得超参$\alpha$和$\beta$的值，因此它也是一个潜在变量(潜在变量就是$\alpha$和$\beta$)模型，我们的目标就是最大化证据函数p(t | alpha, beta)，又因为参数矢量w是边缘化的，我们将其也作为潜在变量：
  $p(t|w,\beta ,X)={\prod }_{n=1}^N\mathcal{N}(t_n;w^T\varphi (x_n),{\beta }^{-1})$
  $p(w|\alpha )=\mathcal{N}(w;0,{\alpha }^{-1}I)$
  $p(t|\alpha ,\beta ,X)=\int p(t|w,\beta )p(w|\alpha )dw$

• 简便的EM算法来找出机器学习参数($\alpha$和$\beta$)

  • E步骤：计算关于潜在变量w的先验概率
    $p(w|t,\alpha ,\beta )=\mathcal{N}(w;m,S),m=\beta S{\Phi }^{T}t,S^{-1}=\alpha I+\beta {\Phi }^T\Phi$，此处
    -M步骤：使用完整数据的似然估计更新参数
    ${\alpha }^{-1}=(1/M)(m^{T}m+TrS)$
    ${\beta }^{-1}=(1/N){\sum }_{n=1}^N{t_n-m^T\varphi (x_n)}^2$

一般化的EM算法

总的来说，EM算法即是是找出拥有潜在变量概率模型的最大似然解的手段，现在提供一个一般化的处理过程：

• 令X为观测数据，Z为潜在变量， $\theta$为参数， 联合分布为p(X, Z | theta)。
• 我们的目标：最大化观测数据的边缘log似然
  $lnp(X|\theta )=ln\left\{\begin{array}{c}{\sum }_{Z}p(X,Z|\theta )\end{array}\right\}$ (式9.29)
• 但问题在于最大化完整数据(即若潜在变量Z事先知道)的似然函数$p(X,Z|\theta )$简单，但最大化$p(X|\theta )$却很难。
• 引入关于潜在变量的分布q(Z)，对任意q(Z)都可以将将似然分解为：
  $lnp(X|\theta )=\mathcal{L}(q,\theta )+KL(q(Z)||p(Z|X,\theta )),$
  此处KaTeX parse error: Expected group after '\frac' at end of input: …X, Z | \theta)},
  $KL(q||p)是q(Z)和后验分布p(Z|X,\theta )的KL散度，包含Z对给定X的条件分布，式子-{\sum }_{Z}q(Z)ln\frac{p(Z|X,\theta )}{q(Z)}.$
  $\mathcal{L}(q,\theta )={\sum }_{Z}q(Z)ln\frac{p(X,Z|\theta )}{q(Z)},$
  $KL(q(Z)||p(Z|X,\theta ))=-{\sum }_{Z}q(Z)ln\frac{p(Z|X,\theta )}{q(Z)}\ge 0.$

EM边界

• $\mathcal{L}(q,\theta )$是该数据log似然的下界，式子如下
  • $-\mathcal{L}(q,\theta )$称为变分自由能，
    $\mathcal{L}(q,\theta )=lnp(X|\theta )-KL(q(Z)||p(Z|X,\theta )\le lnp(X|\theta )$
• EM算法对$\mathcal{L}$进行了坐标上升法
  • E步骤：固定$\theta$，使用q去最大化L
  • M步骤：固定q，使用$\theta$去最大化L
    如下图所示为分解结果，此处q(Z)可以为任意分布，因为KL散度满足条件$KL(q||p)\ge 0$，则$\mathcal{L}(q,\theta )$为似然函数$lnp(X|\theta )$的下界：
    

E步骤

• E步骤：固定$\theta$，使用q去最大化$\mathcal{L}(q,\theta )$
  $\mathcal{L}(q,\theta )=lnp(X|\theta )-KL(q(Z)||p(Z|X,\theta ))$
• 使用$q(Z)<-p(Z|X,\theta )$最大化$\mathcal{L}$
  如下图为EM算法的E步骤，q分布设置等于对当前参数${\theta }^{old}$的后验分布，这使得KL散度消失，进而下界移动到log似然函数的相同值。
  
  如下图，展示了EM算法的M步骤，分布q(Z)保持固定，然后下界L(q, theta)对应修正值theta new的theta参数矢量的最大值。因为KL散度为非负值，这将造成log似然ln p(X | theta)至少增加到下界。
  

参数空间的图片

• E步骤重置边界$\mathcal{L}(q,\theta )为lnp(X|\theta )(\theta ={\theta }^{(}old))$，它即是：
  • 固定在$\theta ={\theta }^{old}$
  • 正切于点$\theta ={\theta }^{old}$
  • 对指数族混合成分在$\theta$是凸的
    如下图所示，多次计算log似然函数对当前参数值的下界然后最大化该边界去获得新的参数值(为似然函数和下界的切点)：
    

最后的思考

• $\mathcal{L}(q,\theta )$的(局部)最大值取决于 $lnp(X|\theta )$
• EM收敛于似然(局部最大值)
  • 对$\mathcal{L}(q,\theta )$的局部上升法，并在E步骤之后$\mathcal{L}=lnp(X|\theta )$
• 多种策略来优化边界：
  • 归纳的EM：从最大化放松M步骤以提升$\mathcal{L}$
  • 期望条件最大化：M步骤轮流使用不同参数组合最大化
  • 增加EM：对每个数据点进行E步骤，进而增加了M步骤
  • 变分EM：从最大化放松E步骤来提升$\mathcal{L}$
    • E步骤之后不再使$\mathcal{L}=lnp(X|\theta )$
• 同样应用于MAP估计 $p(\theta |X)=p(\theta )p(X|\theta )/p(X)$
  • 边界第二项：$lnp(\theta |X)\ge lnp(\theta )+\mathcal{L}(q,\theta )-lnp(X)$
    MAP是Maximum posterior即最大后验。