机器学习（十一）——高斯混合模型和EM算法

高斯混合模型和EM算法

GMM与概率分布（续）

Andrew讲义在介绍GMM的时候，由于要和K-means做对比，因此主要介绍了GMM在聚类中的应用，但GMM的应用并不仅于此。

我们知道函数有泰勒级数和傅立叶级数等展开方式，同样的任意随机变量的概率分布，也可以展开为若干高斯分布的和：

p(x)=∑m=1McmN(x;μm,σm2)

不难看出上式和上节的公式5是等价的。

GMM的这种用法在语音识别领域用的较多，是GMM-HMM算法的核心思想之一。

http://cseweb.ucsd.edu/~atsmith/project1_253.pdf

Clustering With EM and K-Means

https://mp.weixin.qq.com/s/wk3_wG1xSuMX1HjJeEgErQ

kmeans聚类理论篇K的选择（轮廓系数）

https://mp.weixin.qq.com/s/tLcF7_jjl3_FmjM7L1kcfw

一文详解高斯混合模型原理

https://mp.weixin.qq.com/s/RpEgTMRUM4GWhIdIPvkGTg

优于VAE，为万能近似器高斯混合模型加入Wasserstein距离

https://zhuanlan.zhihu.com/p/85338773

高斯混合模型（GMM）推导及实现

https://mp.weixin.qq.com/s/Uayd-KiXeoL5vBUyC1-oDQ

小孩都看得懂的GMM

本节将进一步讨论EM算法的性质，并将之应用到使用latent random variables的一大类估计问题中。

Jensen不等式

首先我们给出凸函数的定义：

设f是定义域为实数的函数，如果对于所有的实数x，f″(x)≥0，那么f是凸函数。当x是向量时，如果其Hessian矩阵H是半正定的（H≥0），那么f是凸函数。如果f″(x)>0或H>0，那么称f是严格凸函数。

Jensen不等式表述如下:

如果f是凸函数，X是随机变量，那么E[f(X)]≥f(EX)。

特别地，如果f是严格凸函数，那么E[f(X)]=f(EX)，当且仅当p(X=EX)=1。也就是说X是常量。

在不引起误会的情况下，E[X]简写做EX。

这个不等式的含义如下图所示：

Jensen不等式应用于凹函数时，不等号方向反向，也就是E[f(X)]≤f(EX)

注：Johan Ludwig William Valdemar Jensen，1859～1925，丹麦人。主业工程师，副业数学家。

EM算法的一般形式

EM算法一般形式的似然函数为：

(1)ℓ(θ)=∑i=1mlog⁡p(x(i);θ)=∑i=1mlog⁡∑zp(x(i),z(i);θ)

根据这个公式直接求θ一般比较困难，但是如果确定了z之后，求解就容易了。

EM算法是一种解决存在隐含变量优化问题的有效方法。既然不能直接最大化ℓ(θ)，我们可以不断地建立ℓ(θ)的下界(E-Step),然后优化下界(M-Step)。

这里首先假设z(i)的分布为Qi。显然：

∑zQi(z)=1,Qi(z)≥0ℓ(θ)=∑i=1mlog⁡∑z(i)p(x(i),z(i);θ)=∑i=1mlog⁡∑z(i)Qi(z(i))p(x(i),z(i);θ)Qi(z(i))(2)≥∑i=1m∑z(i)Qi(z(i))log⁡p(x(i),z(i);θ)Qi(z(i))

这里解释一下最后一步的不等式变换。

首先，根据数学期望的定义公式：

E[X]=∑i=1nxipi∑z(i)Qi(z(i))p(x(i),z(i);θ)Qi(z(i))=E[p(x(i),z(i);θ)Qi(z(i))]

又因为f(x)=log⁡x是凹函数，根据Jensen不等式可得：

f(Ez(i)∼Qi[p(x(i),z(i);θ)Qi(z(i))])≥Ez(i)∼Qi[f(p(x(i),z(i);θ)Qi(z(i)))]

综上，公式2给出了ℓ(θ)的下界。对于Qi的选择，有多种可能，那种更好的？

假设θ已经给定，那么ℓ(θ)的值就取决于Qi(z(i))和p(x(i),z(i))了。我们可以通过调整这两个概率，使下界不断上升，以逼近ℓ(θ)的真实值。当不等式变成等式时，下界达到最大值。

由Jensen不等式相等的条件可知，ℓ(θ)的下界达到最大值的条件为：

p(x(i),z(i);θ)Qi(z(i))=c

其中c为常数。

这实际上表明：

Qi(z(i))∝p(x(i),z(i);θ)

其中的∝符号是两者成正比例的意思。

从中还可以推导出：

p(x(i),z(i);θ)Qi(z(i))=c=∑zp(x(i),z(i);θ)∑zQi(z(i))

因为∑zQi(z(i))=1，所以上式可变形为：

Qi(z(i))=p(x(i),z(i);θ)∑zp(x(i),z(i);θ)=p(x(i),z(i);θ)p(x(i);θ)=p(z(i)∣x(i);θ)

可见，当Qi(z(i))为z(i)的后验分布时，ℓ(θ)的下界达到最大值。

因此，EM算法的过程为：

Repeat until convergence {

(E-step) For each i：

Qi(z(i)):=p(z(i)∣x(i);θ)

(M-step) Update the parameters：

θ:=arg⁡maxθ∑i∑zQi(z(i))log⁡p(x(i),z(i);θ)Qi(z(i))

如何保证算法的收敛性呢？

如果我们用θ(t)和θ(t+1)表示EM算法第t次和第t+1次迭代后的结果，那么我们的任务就是证明ℓ(θ(t))≤ℓ(θ(t+1))。

由公式1和2可得：

(3)ℓ(θ)≥∑i∑z(i)Qi(z(i))log⁡p(x(i),z(i);θ)Qi(z(i))

由之前的讨论可以看出，E-Step中的步骤是使上式的等号成立的条件，即：

ℓ(θ(t))=∑i∑z(i)Qi(t)(z(i))log⁡p(x(i),z(i);θ(t))Qi(t)(z(i))

因为公式3对于任意Qi和θ都成立，因此：

ℓ(θ(t+1))≥∑i∑z(i)Qi(t)(z(i))log⁡p(x(i),z(i);θ(t+1))Qi(t)(z(i))

因为M-Step的最大化过程，可得：

∑i∑z(i)Qi(t)(z(i))log⁡p(x(i),z(i);θ(t+1))Qi(t)(z(i))≥∑i∑z(i)Qi(t)(z(i))log⁡p(x(i),z(i);θ(t))Qi(t)(z(i))

综上可得：ℓ(θ(t))≤ℓ(θ(t+1))

事实上，如果我们定义：

J(Q,θ)=∑i∑z(i)Qi(z(i))log⁡p(x(i),z(i);θ)Qi(z(i))

则EM算法可以看作是J函数的坐标上升法。E-Step固定θ，优化Q；M-Step固定Q，优化θ。

https://mp.weixin.qq.com/s/NbM4sY93kaG5qshzgZzZIQ

EM算法的九层境界：Hinton和Jordan理解的EM算法

https://mp.weixin.qq.com/s/FSID_FLMVxrKEj_padg_bQ

EM算法是炼金术吗？

https://mp.weixin.qq.com/s/JSnPTVRgXO3aih5srU1eZQ

EM算法原理总结

https://mp.weixin.qq.com/s/192sLXAvLKzwsTKCZs-AuA

一文详尽系列之EM算法

https://mp.weixin.qq.com/s/D2Am1tyCt3DQLvF0ylLJrg

深入浅出聚类算法！如何对王者英雄聚类分析，探索英雄之间的秘密

https://mp.weixin.qq.com/s/PKfsSMMtwks_KqRkpKJbKA

理解高斯混合模型中期望最大化的M-Step

重新审视混合高斯模型

下面给出混合高斯模型各参数的推导过程。

E-Step很简单：

wj(i)=Qi(z(i)=j)=p(z(i)=j∣x(i);ϕ,μ,Σ)

在M-Step中，我们将各个变量和分布的概率密度函数代入，可得：

∑i=1m∑z(i)Qi(z(i))log⁡p(x(i),z(i);θ)Qi(z(i))=∑i=1m∑j=1kQi(z(i)=j)log⁡p(x(i)∣z(i)=j;μ,Σ)p(z(i)=j;ϕ)Qi(z(i)=j)=∑i=1m∑j=1kwj(i)log⁡1(2π)n/2|Σj|1/2exp⁡(−12(x(i)−μj)TΣj−1(x(i)−μj))⋅ϕjwj(i)

对μl求导，可得：

∇μl∑i=1m∑j=1kwj(i)log⁡1(2π)n/2|Σj|1/2exp⁡(−12(x(i)−μj)TΣj−1(x(i)−μj))⋅ϕjwj(i)=−∇μl∑i=1m∑j=1kwj(i)12(x(i)−μj)TΣj−1(x(i)−μj)(4)=12∑i=1mwl(i)∇μl(2μlTΣl−1x(i)−μlTΣl−1μl)=∑i=1mwl(i)(Σl−1x(i)−Σl−1μl)

这里对最后一步的推导，做一个说明。为了便于以下的讨论，我们引入符号“tr”，该符号表示矩阵的主对角线元素之和，也被叫做矩阵的“迹”（Trace）。按照通常的写法，在不至于误会的情况下，tr后面的括号会被省略。

tr的其他性质还包括：

(5.1)tra=a,a∈R(5.2)tr⁡(A+B)=tr⁡(A)+tr⁡(B)(5.3)tr⁡(cA)=ctr⁡(A)(5.4)tr⁡(A)=tr⁡(AT)(5.5)tr⁡(AB)=tr⁡(BA)(5.6)tr⁡(ABCD)=tr⁡(BCDA)=tr⁡(CDAB)=tr⁡(DABC)(5.7)∇Atr⁡(AB)=BT(5.8)∇ATf(A)=(∇Af(A))T(5.9)∇Atr⁡(ABATC)=CAB+CTABT(5.10)∇ATtr⁡(ABATC)=BTATCT+BATC(5.11)∇A∣A∣=∣A∣(A−1)T

因为μlTΣl−1μl是实数，由公式5.1可得：

∇μlμlTΣl−1μl=∇μltr⁡(μlTΣl−1μl)

由公式5.10可得：

∇μltr⁡(μlTΣl−1μl)=(Σl−1)T(μlT)TIT+Σl−1(μlT)TI=(Σl−1)Tμl+Σl−1μl

因为Σl−1是对称矩阵，因此，综上可得：

(5.12)∇μlμlTΣl−1μl=2Σl−1μl

回到正题，令公式4等于0，可得：

μj:=∑i=1mwj(i)x(i)∑i=1mwj(i)

同样的，对ϕj求导，可得：

∑i=1m∑j=1kwj(i)log⁡ϕj

因为存在∑j=1kϕj=1这样的约束，因此需要使用拉格朗日函数：

L(ϕ)=∑i=1m∑j=1kwj(i)log⁡ϕj+β(∑j=1kϕj−1)

这里的β是拉格朗日乘子，ϕj≥0的约束条件不用考虑，因为对log⁡ϕj求导已经隐含了这个条件。

∂L(ϕ)∂ϕj=∑i=1mwj(i)ϕj+β

令上式等于0，可得：

∑i=1mwj(i)ϕj=−β=∑i=1m∑j=1kwj(i)∑j=1kϕj

因为∑j=1kϕj=1,∑j=1kwj(i)=1，所以：

−β=∑i=1m11=mϕj:=1m∑i=1mwj(i)