从EM算法的典型应用GMM说起,需要知悉的几个点
因为毕设一半的工作量是基于GMM的,代码方面也是用c++复现了一遍GMM模型的参数更新公式。大概知道GMM模型的参数更新公式是基于EM算法得出的。但之前一直没有把GMM模型和EM算法的联系真正建立起来,一来也是懒,二来是以前对贝叶斯条件概率那一套没有很好理解。
其实将GMM模型用条件概率的思想进行理解,这就是GMM模型的贝叶斯理解,而利用贝叶斯条件概率得出的生成模型,如果含有隐变量,就自然可以用EM算法来解,(是不是联想到了LDA)(当然有凸函数的限制),而EM算法能保证训练参数的收敛是靠数学上凸函数的Jensen不等式来证明的。
下面大致从6个方面来复习EM算法和复习GMM
- 高斯混合模型(GMM)的表现形式
- GMM的应用
- GMM的贝叶斯理解
- EM算法估计给出的GMM参数更新闭式表达式。
- EM算法收敛性:Jensen不等式
- kmeans和GMM的潜在联系,附带kmeans的优化
高斯混合模型(GMM)的表现形式
设有随机变量
X,则混合高斯模型可以用下式表示:
p(x)=k=1∑KπkN(x∣μk,Σk)
其中
N(x∣μk,Σk)称为混合模型中的第k个分量(component)。
πk是混合系数(mixture coefficient)
满足
∑k=1Kπk=1 ,
0<=πk<=1 ,
实际上,可以认为
πk就是每个分量
N(x∣μk,Σk)的权重。
GMM的应用
GMM常用于聚类。如果要从 GMM 的分布中随机地取一个点的话,实际上可以分为两步:首先随机地在这 K 个 Component 之中选一个,每个 Component 被选中的概率实际上就是它的系数
πk,选中 Component 之后,再单独地考虑从这个 Component 的分布中选取一个点就可以了──这里已经回到了普通的 Gaussian 分布,转化为已知的问题。
将GMM用于聚类时,假设数据服从混合高斯分布(Mixture Gaussian Distribution),那么只要根据数据推出 GMM 的概率分布来就可以了;然后 GMM 的 K 个 Component 实际上对应K个 cluster 。根据数据来推算概率密度通常被称作 density estimation 。特别地,当我已知(或假定,这里就是假设高斯分布,LDA假设狄利克雷分布)概率密度函数的形式,而要估计其中的参数的过程被称作『参数估计』。
假设有两个聚类,可以定义K=2. 那么对应的GMM形式如下:
p(x)=π1N(x∣μ1,Σ1)+π2N(x∣μ2,Σ2)
上式中未知的参数有六个:
(π1,μ1,Σ1;π2,μ2,Σ2)之前提到GMM聚类时分为两步,第一步是随机地在这K=2分量中选一个,每个分量被选中的概率即为
π1=π2=0.5混合系数,表示每个分量被选中的概率是0.5,即从中抽出一个点,这个点属于第一类的概率和第二类的概率各占一半。但实际应用中事先指定
πk的值是很笨的做法,只是第一轮得不得这样初始化而已。当问题一般化后,会出现一个问题:当从图2中的集合随机选取一个点,怎么知道这个点是来
自N(x∣μ1,Σ1)还是
N(x∣μ2,Σ2)呢?换言之怎么根据数据自动确定
π1,π2的值?这就是GMM参数估计的问题。要解决这个问题,可以使用EM算法。通过EM算法,我们可以迭代计算出GMM中的参数:
(πk,μk,Σk)
GMM的贝叶斯理解
上面提到的,“先随机地在这K分量中选一个”,“每个分量被选中的概率” ,这就已经明显感觉到了条件概率的味道了,对,这里我们着重改写GMM的形式,方便我们引出GMM背后的条件概率,后验概率这些东西,从而得到GMM的贝叶斯理解。
GMM的原始形式如下:
p(x)=k=1∑KπkN(x∣μk,Σk)(1)
前面提到
πk可以看成是第k类被选中的概率。我们引入一个新的K维随机变量
z.(贝叶斯公式的味道)。但这个随机变量很特殊,是个二值分布
zk(1≤k≤K)只能取0或1两个值。
zk=1表示第k类被选中的概率,即:
p(zk)=πk.
如果
zk=0表示第k类没有被选中的概率。更数学化一点,
zk要满足以下两个条件:
zk∈{0,1},K∑zk=1
对k=2有两类的例子,则
z的维数是2.
如果从第一类中取出一个点,则
z=(1,0),
如果从第二类中取出一个点,则
z=(0,1),
Zk是用来表示第k类被选中的随机变量,只能取0和1两个值(随机变量是要取值的)。
所以变量
Zk=1应该是表示“第k类被选中”,而非“第k类被选中的概率”。
但
P(Zk=1)表示第k类被选中的概率。
zk=1的概率
P(Zk=1)是
πk,我们记为
p(zk)=πk=πkz1=1
假设
zk之间是独立同分布的(iid),我们可以写
z的联合概率分布形式,就是连乘:
p(z)=p(z1)p(z2)...p(zK)=k=1∏Kπkzk(2).
因为
zk只能取0或1,且
z中只能有一个
zk为1而其它
(j=k)全为0,所以上式是成立的。(?,这是在解释独立iid吗?解释为什么可以连乘是吗?)
以上其实是我认为 GMM的贝叶斯理解中最难的部分。,因为假设看不懂,后面越糊涂
类确定了(条件),我们可以在类下看数据但高斯分布分布,这个叙述其实可以用条件概率来表示了:
p(x∣zk=1)=N(x∣μk,Σk)
即每一类下数据服从高斯分布,进一步还可以写成:
p(x∣z)=N(x∣μk,Σk)zk(3)
上面式2和式3分别给出了
p(z) 和
p(x∣z)的公式,根据条件概率公式,我们得出观测数据的联合分布:
p(x)=z∑p(z)∗p(x∣z)=z∑(k=1∏KπkzkN(x∣μk,Σk)zk)=K=1∑KπkN(x∣μk,Σk)(4)
上式从第二个等号到第三个等号其实并不难,第二个等号对z求和,其实就是
∑K=1K,而里面,只要
i=k,则
zi=0,所以z_k=0$的项为1,可忽略,最终得到了第三个等式
可以看到GMM模型的(1)式与(4)式有一样的形式,且(4)式中引入了一个新的变量
z,通常称为隐含变量(latent variable)。对于一个要聚类的数据来说,『隐含』的意义是:我们知道数据可以分成两类,但是随机抽取一个数据点,我们不知道这个数据点属于第一类还是第二类,它的归属我们观察不到,因此引入一个隐含变量
z来描述这个现象。
再把视角拉回贝叶斯,注意到在贝叶斯的思想下,
p(z)是先验概率,
p(x∣z)是似然概率,很自然我们会想到求出后验概率
p(z∣x):
γ(zk)=p(zk=1∣x)==p(x,zk=1)p(zk=1)p(x∣zk=1)==∑j=1Kp(zj=1)p(x∣zj=1)p(zk=1)p(x∣zk=1)==∑j=1KπjN(x∣μj,Σj)πkN(x∣μk,Σk)(5)
(上面的第二行,贝叶斯定理,关于这一行的分母,很多人会疑问,应该是
p(x,zk=1)还是
p(x),按照正常的写法,肯定是
p(x),但为了强调
zk的取值,有的书会写成
p(x,zk=1),比如李航的小蓝书,这里就约定
p(x)与
p(x,zk=1)是等同的。
第三行全概率公式,最后一个等号,结合三四式.
上式中,我们定义了符号
γ(zk)来表示第k个分量的后验概率。在贝叶斯的观点下,
πk可视为
zk=1的先验概率.
上述内容改写了GMM的形式,并引入了隐含变量
z和已知
x后的的后验概率
γ(zk),这样做是为了方便我们使用EM算法来估计GMM的参数。其实EM算法 过程中,我们能看到除了GMM中高斯分布的参数在更新,
隐含变量
z的先验概率
πk和已知
x后的后验概率
γ(zk)也在不断更新。
EM算法估计给出的GMM参数更新闭式表达式。
注意,式8和式9中,N代表点的数量,
γ(znk)表示点
n(xn)属于聚类k的后验概率,则,
Nk代表属于第
k个聚类的点的数量。
uk代表属于第
k个聚类的所有点的加权平均,每个点的权重是
γ(znk),即这个点属于聚类k的后验概率大小,跟第k个聚类有关。
EM算法估计GMM参数即最大化(8),(10)和(12)。需要用到(5),(8),(10)和(12)四个公式。我们先指定
π,
μ和
Σ的初始值,带入(5)中计算出
γ(znk),然后再将
γ(znk)带入(8),(10)和(12),求得
πk,
μk和
Σk;接着用求得的
πk,
μk和
Σk;再带入(5)得到新的
γ(znk),再将更新后的
γ(znk)带入(8),(10)和(12),如此往复,直到算法收敛。
好了,看到这里了,GMM和贝叶斯的关系说清楚了,GMM的参数求解用EM框架给出了流程和闭式表达式了,还有什么呢?有没有看到流程的第五个的“收敛”两个字?
我么难道就不怕这样迭代下去不收敛吗?
所以我们接下来分析下EM 算法的收敛性,准确来说或者具体来说,是含隐变量问题求解过程中的参数更新的收敛性问题。
EM算法收敛性:Jensen不等式
上面谈收敛,首先要弄清楚 ,收敛是对谁说的,即是什么东西收敛,其实这里的收敛是针对参数更新过程中的目标函数而言的。
其实我们在上面讨论时,忽略了谈目标函数的问题,
但其实也提到了目标函数的影子,那就是最大化似然函数。
我们的参数更新过程中,似然函数是会不断变大直至收敛,还是一会儿变大,一会儿变小呢?
从数学的角度分析一下,哎,凸优化知识还有五秒钟抵达战场。。。
还是回忆一下凸函数的图吧,
上面如果t=1/2,就很 直观了
f(21x1+21x2)<=21f(x1)+21f(x2)(1)
一般的凸函数定义是
f(λx1+(1−λ)x2)<=λf(x1)+(1−λ)f(x2)(2)
jensen不等式,是对凸函数,把上式扩展到了,无穷多点,即有无穷多个
λ
f(i=1∑Mλixi)<=i=1∑Mλif(xi)(3)
i=1∑Mλi=1,λi>0
公式(3)被称为 Jensen 不等式,它是式(2)的泛化形式。
因为
∑i=1Mλi=1,λi>0,从概率论的角度 来说,如果把
λi看成是取值为
xi的离散变量x的概率分布的话,那么就可以写成期望的形式,如下:
f(E[x])<=E[f(x)](4)
E[.]表示期望。
同样当函数变为凹函数时,上式的不等号只是变了个方向而已,
f(i=1∑Mλixi)>=i=1∑Mλif(xi)(5)
f(E[x])>=E[f(x)](6)
为什么要这个,因为在EM算法或者GMM算法的利用Jensen不等式证明收敛性用到的是式子5和式子6这两个不等式。因为EM算法或者GMM算法的利用Jensen证明不等式证明收敛性用到的是Log函数的凹凸性。Log函数是凹函数,二阶导恒小于0.Jensen不等式在Log函数的直观表现如下所示,
下面开始证明之旅。
给定训练样本
{x(1),x(2),,...,,x(m)},样本间独立,我们想找到每个样本的隐含类别,能使得
p(x,z)最大。于是
p(x,z)的最大似然估计是:
l(θ)=i=1∑mlogp(x;θ)i=1∑mlogz∑p(x,z;θ)
第一步是对极大似然取对数,第二步是对每个样例的每个可能类别z求联合分布概率和。但是直接求
l(θ)一般比较困难,因为有隐藏变量z存在,但是一般确定了z后,求解就容易了。
EM是一种解决存在隐含变量优化问题的有效方法。既然不能直接最大化
l(θ),我们可以不断地建立
l(θ)的下界(E步),然后优化下界(M步)。这句话比较抽象,看下面的
l(θ)=i=1∑logp(xi;θ)=i=1∑logz(i)∑p(x(i),z(i);θ)=i=1∑logz(i)∑Qi(z(i))Qi(z(i))p(x(i),z(i);θ)≥i=1∑z(i)∑Qi(z(i))logQi(z(i))p(x(i),z(i);θ)(1)
(1)到(2)比较直接,就是分子分母同乘以一个相等的函数。(2)到(3)利用了Jensen不等式,考虑到Log函数是凹函数(二阶导数小于0),而且
这个过程可以看作是对
l(θ)求了下界。对于
Qi的选择,有多种可能,那种更好的?假设
θ已经给定,那么
l(θ)的值就决定于
Qi(z(i))和
p(x(i),z(i))了。我们可以通过调整这两个概率使下界J不断上升,
Qi(z(i))是可变的,不同的
Qi(z(i))和
p(x(i),z(i))对应了不同的下界J),以逼近
l(θ)的真实值,那么什么时候算是调整好了呢?当不等式变成等式时,说明我们调整后的概率能够等价于
l(θ)了。按照这个思路,我们要找到等式成立的条件。根据Jensen不等式,要想让等式成立,需要让随机变量变成常数值,这里得到:
Qi(z(i))p(x(i),z(i);θ)=常数c,Qi(z(i))=cp(x(i),z(i);θ)
c为常数,不依赖
z(i),又因为,
∑zQi(z(i))=1,所以,
p(x(i),z(i);θ)=c∗Qi(z(i))
∑zp(x(i),z(i);θ)=∑zc∗Qi(z(i))=c∗∑zQi(z(i))=c
带入
Qi(z(i))=cp(x(i),z(i);θ)
得到:
Qi(z(i))=cp(x(i),z(i);θ)=∑zp(x(i),z(i);θ)p(x(i),z(i);θ)=p(x(i),θ)p(x(i),z(i);θ)=p(z(i)∣x(i),θ)
至此,我们推出了在固定其他参数
θ后,最优的
Qi(z(i))(也就是让Jensen不等式等号成立的
Qi(z(i)))的计算公式就是后验概率,
解决了
Qi(z(i))如何选择的问题。这一步就是E步,建立
l(θ)的一个比较好比较紧的下界J。
接下来的M步,就是在给定
Qi(z(i)),也就是给定了下界函数J后,调整
θ,去极大化下界函数J(从而达到间接极大化
l(θ)的目的),(在固定
Qi(z(i))后,下界还可以调整的更大,靠调整
θ)。
那么究竟怎么确保EM收敛?假定
θ(t)和
θ(t+1)是EM第t次和t+1次迭代后的结果。如果我们证明了
l(θ(t))≤l(θ(t+1)),也就是说极大似然估计单调增加,那么最终我们会到达最大似然估计的最大值。
其实综合上面的知识就比较好理解了,
从上面的推导,我一步我们用了最优的
Qi(z(i)),也就是后验概率
p(z(i)∣x(i),θ),能让Jensen不等式等号成立。
l(θ)=i=1∑z(i)∑Qi(z(i))logQi(z(i))p(x(i),z(i);θ)
所以:
l(θ(t+1))≥i=1∑z(i)∑Qi(z(i))logQi(z(i))p(x(i),z(i);θ(t+1))(4)
≥i=1∑z(i)∑Qi(z(i))logQi(z(i))p(x(i),z(i);θ(t))(5)
=l(θ(t))(6)
上式还是可读,不过是不是有疑问,为什么4式子能大于等于,而6式子就只能等于(中间的5式子是普通的参数估计极大化过程,好理解),因为6式子用了最优的
Qi(z(i)),也就是后验概率
p(z(i)∣x(i),θ),能让Jensen不等式等号成立。而6式最优的
Qi(z(i))不一定在下一次迭代中,同样是最优的,所以需要调整,所以不是等号,而是大于等于。
这样就证明了
l(θ(t))会单调增加。一种收敛方法是
l(θ(t))不再变化,还有一种就是变化幅度很小。
再次解释一下(4)、(5)、(6)。首先(4)对所有的参数都满足,而其等式成立条件只是在固定
θ,并调整好Q时成立,而第(4)步只是固定Q,调整
θ,不能保证等式一定成立。(4)到(5)就是M步的定义,(5)到(6)是前面E步所保证等式成立条件。也就是说E步会将下界拉到与
l(θ)一个特定值(这里
θ)一样的高度,而此时发现下界仍然可以上升,因此经过M步后,下界又被拉升,但达不到与
l(θ)另外一个特定值一样的高度,之后E步又将下界拉到与这个特定值一样的高度,重复下去,直到最大值。
再啰嗦一句,见上图,在
θt处我们又多种
Qi(z(i))可选,不同的
Qi(z(i))对应不同的下界函数J(图中绿色和蓝色就是两种
Qi(z(i))选择下不同的下界函数J),我们会选出最优的
Qi(z(i)),让这个
Qi(z(i))对应的下界函数正好能在
θt处和原目标函数
l(θt)相等(其他地方不能保证),这样的
Qi(z(i))被我们推导出就是当前的后验概率,然后我们对这个下界函数J,通过参数估计,求极大似然,这里的参数就是
θ,于是我们来到了一个新的地方
θ(t+1).并继续我们的更新之路。
总结
如果将样本看作观察值,潜在类别看作是隐藏变量,那么聚类问题也就是参数估计问题,只不过聚类问题中参数分为隐含类别变量和其他参数,这犹如在x-y坐标系中找一个曲线的极值,然而曲线函数不能直接求导,因此什么梯度下降方法就不适用了。但固定一个变量后,另外一个可以通过求导得到,因此可以使用坐标上升法,一次固定一个变量,对另外的求极值,最后逐步逼近极值。对应到EM上,E步估计隐含变量,M步估计其他参数,交替将极值推向最大。EM中还有“硬”指定和“软”指定的概念,“软”指定看似更为合理,但计算量要大,“硬”指定在某些场合如K-means中更为实用(要是保持一个样本点到其他所有中心的概率,就会很麻烦)。
另外,EM的收敛性证明方法确实很牛,能够利用log的凹函数性质,还能够想到利用创造下界,拉平函数下界,优化下界的方法来逐步逼近极大值。而且每一步迭代都能保证是单调的。最重要的是证明的数学公式非常精妙,硬是分子分母都乘以z的概率变成期望来套上Jensen不等式,前人都是怎么想到的。