EM算法可以说是一个非常经典的算法,至今仍在广泛地被使用(如强化学习领域等等)。
网上介绍该算法的文章也很多,比如如何通俗理解EM算法、【机器学习基础】EM算法。但是我认为这些文章讲的太多,反而显得乱。像教科书周志华老师的《机器学习》对于EM算法写得十分简略,李航老师的《统计机器学习》对于 Q ( Z ) Q(Z) Q(Z)的引出没有讲解完整。个人认为在阅读了两位老师对于该算法的讲解后,概率图模型网络参数学习—含隐变量的参数估计(EM算法)这篇文章是最好的补充,如果已经对极大似然估计有了解,可以直接阅读该文。本文是对以上资料的总结和展开,主要参考概率图模型网络参数学习—含隐变量的参数估计(EM算法)这篇文章。
EM 算法(全称为Expectation-Maximization algorithm)是一种对含有隐变量的问题的参数估计算法。
参数估计方法可以参考 抽样分布和参数估计(极大似然估计)。一般而言,在EM算法采用的参数估计方法是极大似然估计。
那什么是含有隐变量的问题?
我们首先来举一个没有隐变量的问题——投硬币问题。我们现在有一个正面出现概率未知的硬币,我们不断地投掷该硬币,得到一系列样本。现在我们需要根据样本估计该硬币正面出现的概率。这是一个没有隐变量的参数估计问题。
含有隐变量的问题: 现在假设有两枚正面出现概率不同的硬币A、B,有个人每次从两个硬币中选择一个硬币进行投掷,得到一系列样本。现在我们需要根据样本估计硬币A、B正面出现的概率。这是一个有隐变量的问题。隐变量就是说一个不确定未知的变量,在这里就是对应每一次投掷对应的硬币 X , X ∈ A , B X, X\in {A,B} X,X∈A,B。因为如果说X是已知的,那么我们可以把已有的属于A和B样本序列分开,然后单独进行估计(即退化为上面所描述的没有隐变量的问题)。
在第三部分我们讲过给出以上两个问题的求解。
这种含有隐变量的问题还有很多,比如NLP方向的同学可能知道的隐马尔科夫模型中的隐状态在估计隐马尔可夫模型参数的时候其实就是隐变量。
参考李航老师的《统计机器学习》。定义X为观测变量集合,Z为隐变量集合, θ \theta θ为参数。
跟不含隐变量的参数估计一样,EM算法其实也采用了极大似然估计。
根据极大似然估计的思想,我们的目标是极大化完全数据X+Z关于参数
θ
\theta
θ的对数似然函数,即极大化
L
(
θ
∣
X
,
Z
)
=
log
P
(
X
,
Z
∣
θ
)
(1)
L(\theta|X,Z) =\log P(X,Z \mid \theta) \tag{1}
L(θ∣X,Z)=logP(X,Z∣θ)(1)
上式L表示Likelihood。由于Z为隐变量,上式无法直接求解,我们只能退而求其次,极大化观测数据X关于参数
θ
\theta
θ的对数似然函数,即极大化
L
(
θ
)
=
log
P
(
X
∣
θ
)
=
log
∑
Z
P
(
X
,
Z
∣
θ
)
(2)
\begin{aligned} L(\theta) &=\log P(X \mid \theta)=\log \sum_{Z} P(X, Z \mid \theta) \end{aligned} \tag{2}
L(θ)=logP(X∣θ)=logZ∑P(X,Z∣θ)(2)
该公式也称为已观测数据的对数“边际似然”(marginal likelihood)。注意,这一极大化的主要困难在于式(2)中有未观测数据Z并包含和的对数。除非
p
(
X
,
Z
∣
θ
)
p(X, Z|\theta)
p(X,Z∣θ)的形式非常简单,否则这个求和难以直接计算。 因此,含有隐变量时,直接进行最大似然估计行不通(如何计算
l
o
g
p
(
X
∣
θ
)
log p(X|\theta)
logp(X∣θ)成为关键)。
为了计算
l
o
g
p
(
X
∣
θ
)
log p(X|\theta)
logp(X∣θ),我们引入一个额外的变分函数Q(Z),Q(Z)为定义在隐变量 Z上的分布。此时对数边际似然函数改写为
log
p
(
X
∣
θ
)
=
log
∑
Z
Q
(
Z
)
P
(
X
,
Z
∣
θ
)
Q
(
Z
)
≥
∑
Z
Q
(
Z
)
log
P
(
X
,
Z
∣
θ
)
Q
(
Z
)
≜
E
L
B
O
(
Q
,
X
∣
θ
)
,
(3)
\begin{aligned} \log p({X} \mid \theta) &=\log \sum_{{Z}} Q({Z}) \frac{P({X}, {Z} \mid \theta)}{Q({Z})} \\ & \geq \sum_{{Z}} Q({Z}) \log \frac{P({X}, {Z} \mid \theta)}{Q({Z})} \\ & \triangleq E L B O(Q, {X} \mid \theta), \end{aligned} \tag{3}
logp(X∣θ)=logZ∑Q(Z)Q(Z)P(X,Z∣θ)≥Z∑Q(Z)logQ(Z)P(X,Z∣θ)≜ELBO(Q,X∣θ),(3)
公式(3)中第二行使用了Jensen不等式。
Jensen 不等式: 即对于凸函数g,有 g ( E [ x ] ) ≤ E [ g ( x ) ] g(E[x])\leq E[g(x)] g(E[x])≤E[g(x)]。最简单的一个例子就是对于凸函数 g ( x ) = x 2 g(x)=x^2 g(x)=x2,对任意 x 1 , x 2 x_1,x_2 x1,x2, 有 g ( x 1 + x 2 2 ) ≤ g ( x 1 ) + g ( x 2 ) 2 g(\frac{x_1+x_2}{2})\leq \frac{g(x_1)+g(x_2)}{2} g(2x1+x2)≤2g(x1)+g(x2)
因此这里引入Q(Z)就是为了能满足Jensen不等式,试想,Q(Z)可视为变量Z的概率分布, p ( X , Z ∣ θ ) Q ( Z ) \frac{p({X}, {Z} \mid \theta)}{Q({Z})} Q(Z)p(X,Z∣θ)可视为一个与Z有关的变量,那么 ∑ Z Q ( Z ) p ( X , Z ∣ θ ) Q ( Z ) = E Z ∼ Q ( Z ) [ p ( X , Z ∣ θ ) Q ( Z ) ] \sum_{{Z}} Q({Z}) \frac{p({X}, {Z} \mid \theta)}{Q({Z})} =E_{Z\sim Q(Z)}[\frac{p({X}, {Z} \mid \theta)}{Q({Z})}] ∑ZQ(Z)Q(Z)p(X,Z∣θ)=EZ∼Q(Z)[Q(Z)p(X,Z∣θ)], 由于log是凹函数,因此公式三成立。
此外,公式3中 E L B O ( Q , X ∣ θ ) E L B O(Q, {X} \mid \theta) ELBO(Q,X∣θ)为对数边际似然函数 l o g p ( X ∣ θ ) log p(X|\theta) logp(X∣θ)的下界,称为证据下界(Evidence Lower Bound, ELBO)。
有Jensen不等式的性质可知,仅当Q(z)=P(Z|X,\theta)时,对数边际似然函数和其下界相等,证明如下
有Jensen不等式可知, ∀ x , g ( x ) = c → g ( E [ x ] ) = E [ g ( x ) ] \forall x ,g(x)=c\rightarrow g(E[x])= E[g(x)] ∀x,g(x)=c→g(E[x])=E[g(x)]
因此,对于公式3,当且仅当 log P ( X , Z ∣ θ ) Q ( Z ) = c \log \frac{P({X}, {Z} \mid \theta)}{Q({Z})}=c logQ(Z)P(X,Z∣θ)=c,即 P ( X , Z ∣ θ ) Q ( Z ) = c \frac{P({X}, {Z} \mid \theta)}{Q({Z})}=c Q(Z)P(X,Z∣θ)=c时等号成立。因此有
∑ Z P ( X , Z ∣ θ ) = ∑ Z Q ( Z ) c ∵ ∑ Z Q ( Z ) = 1 即分布的概率和=1 ∴ ∑ Z P ( X , Z ∣ θ ) = c ∵ P ( X , Z ∣ θ ) Q ( Z ) = c ∴ Q ( Z ) = P ( X , Z ∣ θ ) c = P ( X , Z ∣ θ ) ∑ Z P ( X , Z ∣ θ ) = P ( X , Z ∣ θ ) P ( X ∣ θ ) = P ( Z ∣ X , θ ) \begin{aligned} &\sum_{{Z}}P({X}, {Z} \mid \theta)=\sum_{{Z}}Q({Z})c\\ \because &\sum_{{Z}}Q({Z})=1 \quad \text{即分布的概率和=1} \\ \therefore &\sum_{{Z}}P({X}, {Z} \mid \theta)=c\\ \because &\frac{P({X}, {Z} \mid \theta)}{Q({Z})}=c \quad \\ \therefore &Q(Z)=\frac{P({X}, {Z} \mid \theta)}{c}= \frac{P({X}, {Z} \mid \theta)}{\sum_{{Z}}P({X}, {Z} \mid \theta)}= \frac{P({X}, {Z} \mid \theta)}{P(X|\theta)}=P(Z|X,\theta) \end{aligned} ∵∴∵∴Z∑P(X,Z∣θ)=Z∑Q(Z)cZ∑Q(Z)=1即分布的概率和=1Z∑P(X,Z∣θ)=cQ(Z)P(X,Z∣θ)=cQ(Z)=cP(X,Z∣θ)=∑ZP(X,Z∣θ)P(X,Z∣θ)=P(X∣θ)P(X,Z∣θ)=P(Z∣X,θ)
至此,可将最大化对数边际似然函数
l
o
g
P
(
X
∣
θ
)
log\ P(X|\theta)
log P(X∣θ)的过程分解为两步:
(1) 固定参数
θ
\theta
θ,先找到分布Q(Z)使得
l
o
g
P
(
X
∣
θ
)
=
E
L
B
O
(
Q
,
X
∣
θ
)
log\ P(X|\theta)=E L B O(Q, {X} \mid \theta)
log P(X∣θ)=ELBO(Q,X∣θ)
(2)固定分布Q(Z), 再寻找参数
θ
\theta
θ,极大化
E
L
B
O
(
Q
,
X
∣
θ
)
E L B O(Q, {X} \mid \theta)
ELBO(Q,X∣θ), 由于它是边际似然函数的下界,此时也就是使得
l
o
g
p
(
X
∣
θ
)
log p(X|\theta)
logp(X∣θ)变得更大
以上两个步骤就对应了上面的图像。第一步对应了找到与 L ( θ ) L(\theta) L(θ)相切的ELBO函数,根据以上的证明,使得等号成立的 Q ( Z ) = P ( Z ∣ X , θ ) Q(Z)=P(Z|X,\theta) Q(Z)=P(Z∣X,θ)。第二步对应了找到ELBO函数的极大值( θ 2 \theta_2 θ2处),此时也使得 l o g p ( X ∣ θ ) log\ p(X|\theta) log p(X∣θ)变得更大。(值得注意的是:只有令 Q ( Z ) = P ( Z ∣ X , θ ) Q(Z)=P(Z|X,\theta) Q(Z)=P(Z∣X,θ), 即两函数图像相切,才能保证ELBO函数的极大值处边际似然函数也更大)。
通过迭代以上两个步骤,我们可以不断地寻得更优的 θ \theta θ使得 l o g P ( X ∣ θ ) log\ P(X|\theta) log P(X∣θ)变得更大,并最终寻得最优值(不能保证是全局最优,可能是局部最优)。
最后,我们对以上两个步骤进行形式化
(1)称为E步(expectation)
:该步骤固定参数
θ
t
\theta_t
θt计算分布
Q
(
Z
)
Q(Z)
Q(Z)使得Jensen不等式等号成立, 并计算此时的下界函数。
Q
t
+
1
(
Z
)
=
arg max
Q
E
L
B
O
(
Q
,
X
∣
θ
t
)
=
P
(
Z
∣
X
,
θ
t
)
\begin{aligned}Q_{t+1}(Z)&=\argmax _Q E L B O(Q, {X} \mid \theta_t)\\ &=P(Z|X,\theta_t)\end{aligned}
Qt+1(Z)=QargmaxELBO(Q,X∣θt)=P(Z∣X,θt)
E
L
B
O
(
Q
,
X
∣
θ
)
=
∑
Z
P
(
Z
∣
X
,
θ
t
)
log
P
(
X
,
Z
∣
θ
)
P
(
Z
∣
X
,
θ
t
)
E L B O(Q, {X} \mid \theta) = \sum_{{Z}} P(Z|X,\theta_t) \log \frac{P({X}, {Z} \mid \theta)}{P(Z|X,\theta_t)}
ELBO(Q,X∣θ)=Z∑P(Z∣X,θt)logP(Z∣X,θt)P(X,Z∣θ)
(2) 称为M步(maximization)
:固定
Q
t
+
1
(
Z
)
Q_{t+1}(Z)
Qt+1(Z),找到新的参数
θ
t
+
1
\theta_{t+1}
θt+1极大化下界函数,即
θ
t
+
1
=
arg
max
θ
E
L
B
O
(
Q
t
+
1
,
X
∣
θ
)
\theta_{t+1}=\underset{\theta}{\arg \max } E L B O\left(Q_{t+1}, \mathbf{X} \mid \theta\right)
θt+1=θargmaxELBO(Qt+1,X∣θ)
这里为什么第一步称为expectation(期望)?维基百科称这是一个错误的命名(misnormer),因为这一步实际上是为了简化log求和项的求解难问题而引出的利用Jessen不等式求下界的问题。当然,事实上确实计算了一个期望,因为我们要求的最佳下界其实是 log P ( X , Z ∣ θ ) Q ( Z ) \log \frac{P({X}, {Z} \mid \theta)}{Q({Z})} logQ(Z)P(X,Z∣θ)在 Q ( Z ) = P ( Z ∣ X , θ t ) Q(Z)=P(Z|X,\theta_t) Q(Z)=P(Z∣X,θt)分布下的期望(见公式3)。
周志华老师和李航老师都定义E步为“以当前参数 θ t \theta_t θt推断 P ( Z ∣ X , θ t ) P(Z|X,\theta_t) P(Z∣X,θt),并计算似然对数函数 L ( θ ∣ X , Z ) L(\theta|X,Z) L(θ∣X,Z)关于Z的期望,并将该期望称为Q函数(Q函数 Q ( θ ∣ θ t ) Q(\theta|\theta_t) Q(θ∣θt),它跟前面我们提到的Q(Z)是不同的。Q(Z)是Z的一个概率分布,而Q函数是已知 θ t \theta_t θt、变量为 θ \theta θ的函数)
Q ( θ ∣ θ t ) = E Z ∣ X , θ t L ( θ ∣ X , Z ) = ∑ Z P ( Z ∣ X , θ t ) log P ( X , Z ∣ θ ) Q\left(\theta \mid \theta^{t}\right)=\mathbb{E}_{\mathbf{Z} \mid \mathbf{X}, \theta^{t}} L(\theta \mid \mathbf{X}, \mathbf{Z})= \sum_{{Z}} P(Z|X,\theta_t) \log {P({X}, {Z} \mid \theta)} Q(θ∣θt)=EZ∣X,θtL(θ∣X,Z)=Z∑P(Z∣X,θt)logP(X,Z∣θ)
注意这个Q函数跟我们提到的下界函数 E L B O ( Q , X ∣ θ ) ELBO(Q,X|\theta) ELBO(Q,X∣θ)是等价的。这是因为在第二步中,我们要求得 θ t + 1 \theta_{t+1} θt+1来最大化 E L B O ( Q , X ∣ θ ) ELBO(Q,X|\theta) ELBO(Q,X∣θ), 我们删除对于最大化 θ \theta θ来说无关的项,就得到了Q函数。
θ t + 1 = arg max θ E L B O ( Q t + 1 , X ∣ θ ) = arg max θ ∑ Z P ( Z ∣ X , θ t ) log P ( X , Z ∣ θ ) P ( Z ∣ X , θ t ) = arg max θ ∑ Z [ P ( Z ∣ X , θ t ) log P ( X , Z ∣ θ ) − P ( Z ∣ X , θ t ) log P ( Z ∣ X , θ t ) ] ( 删 除 不 含 θ 的 项 ) = arg max θ ∑ Z P ( Z ∣ X , θ t ) log P ( X , Z ∣ θ ) = arg max θ Q ( θ ∣ θ t ) \begin{aligned} \theta_{t+1}&=\underset{\theta}{\arg \max } E L B O\left(Q_{t+1}, \mathbf{X} \mid \theta\right) \\ &=\underset{\theta}{\arg \max } \sum_{{Z}} P(Z|X,\theta_t) \log \frac{P({X}, {Z} \mid \theta)}{P(Z|X,\theta_t)} \\ &=\underset{\theta}{\arg \max } \sum_{{Z}}[P(Z|X,\theta_t)\log P({X}, {Z} \mid \theta)-P(Z|X,\theta_t)\log P(Z|X,\theta_t)]\quad\quad(删除不含\theta的项)\\ &=\underset{\theta}{\arg \max }\sum_{{Z}} P(Z|X,\theta_t)\log P({X}, {Z} \mid \theta)\\ &=\underset{\theta}{\arg \max }Q(\theta|\theta_t) \end{aligned} θt+1=θargmaxELBO(Qt+1,X∣θ)=θargmaxZ∑P(Z∣X,θt)logP(Z∣X,θt)P(X,Z∣θ)=θargmaxZ∑[P(Z∣X,θt)logP(X,Z∣θ)−P(Z∣X,θt)logP(Z∣X,θt)](删除不含θ的项)=θargmaxZ∑P(Z∣X,θt)logP(X,Z∣θ)=θargmaxQ(θ∣θt)
因此可证下界函数与Q函数等价,为了简化计算,我们可在E步直接求 P ( Z ∣ X , θ t ) P(Z|X,\theta_t) P(Z∣X,θt)和 Q ( θ ∣ θ t ) Q(\theta|\theta_t) Q(θ∣θt),在M步直接极大化 Q ( θ ∣ θ t ) Q(\theta|\theta_t) Q(θ∣θt)(见下).
引自: 如何通俗理解EM算法
举个例子,抛硬币,有两个硬币,但是两个硬币的材质不同导致其出现正反面的概率不一样,目前我们只有一组观测数据,要求出每一种硬币投掷时正面向上的概率。总共投了五轮,每轮投掷五次,现在先考虑一种简单的情况,假设我们知道这每一轮用的是哪一个硬币去投掷的:
那么我们拿着这样的一组数据,就可以很轻松的估计出A硬币和B硬币出现正面的概率,如下:
PA = (3+1+2)/ 15 = 0.4
PB= (2+3)/10 = 0.5
现在把问题变得复杂一点,假设我们不知道每一次投掷用的是哪一种硬币,等于是现在的问题加上了一个隐变量,就是每一次选取的硬币的种类。
那么现在可以想一想,假设我们把每一次硬币的种类设为Z,则这五次实验生成了一个5维的向量 Z = ( z 1 , z 2 , z 3 , z 4 , z 5 ) Z=(z_1,z_2,z_3,z_4,z_5) Z=(z1,z2,z3,z4,z5),现在问题来了,如果我们要根据观测结果去求出PA,PB,那么首先需要知道Z,但是如果用最大似然估计去估计Z,又要先求出PA,PB。这就产生了一个循环。那么这个时候EM算法的作用就体现出来了,EM算法的基本思想是:先初始化一个PA,PB(这里 θ \theta θ就是PA和PB),然后我们拿着这个初始化的PA,PB用最大似然概率估计出Z的分布,接下来有了Z之后就用Z去计算出在当前Z的情况下的PA,PB是多少,然后不断地重复这两个步骤直到收敛。
有了这个思想之后现在用这个思想来做一下这个例子,假设初始状态下PA=0.2, PB=0.7,然后我们根据这个概率去估计出Z的分布(这里
Z
=
(
z
1
,
z
2
,
z
3
,
z
4
,
z
5
)
Z=(z_1,z_2,z_3,z_4,z_5)
Z=(z1,z2,z3,z4,z5)是一个五个变量(下称分量)的集合,每一个分量都是一个二项分布):
按照最大似然估计,Z=(B,A,A,B,A),有了Z之后我们反过来重新估计一下PA,PB:
PA = (2+1+2)/15 = 0.33
PB =(3+3)/10 = 0.6
可以看到PA,PB的值已经更新了,假设PA,PB的真实值0.4和0.5,那么你在不断地重复这两步你就会发现PA,PB在不断地靠近这两个真实值。
注意:以上内容个人认为不太准确
根据公式,这里我们应该考虑Z的所有取值,Z=(B,A,A,B,A)只是其中一种,我们应该依据所有的Z的可能取值求Q函数,然后根据Q函数对PA和PB进行极大似然估计,而不是简单地取Z最优可能的取值对PA、PB进行估计。因此PA、PB应该是五个二项分布的极大似然估计值。