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