本内容主要介绍 EM 算法 以及其公式的详细推导过程。
EM 算法 全称 Expectation Maximization Algorithm,译作 期望最大算法 或最大期望化算法,它是一种迭代算法,用于含有 隐变量(hidden variable)的概率模型参数的极大似然估计,或极大后验概率估计。1977 年由 Dempster 等人总结提出。
1.1 Jensen 不等式
在开始介绍 EM 算法前,我们先介绍和了解一下凹凸函数和 Jensen 不等式。如果您已经了解它们,可以直接跳过。
1.1.1 凸函数和凹函数
凸函数 是一个定义在某个向量空间的子集 C C C(区间)上的实值函数 f f f,如果在其定义域 C C C 上的任意两点 x 1 x_1 x1, x 2 x_2 x2, 0 < t < 1 0<t<1 0<t<1,有
t f ( x 1 ) + ( 1 − t ) f ( x 2 ) ≥ f ( t x 1 + ( 1 − t ) x 2 ) (1.1) t f(x_1) + (1-t)f(x_2) \geq f\Big(t x_1 + (1-t)x_2\Big) \tag{1.1} tf(x1)+(1−t)f(x2)≥f(tx1+(1−t)x2)(1.1)
也就是说凸函数任意两点的割线位于函数图形上面,这也是 Jensen 不等式的两点形式。如果总有 t f ( x 1 ) + ( 1 − t ) f ( x 2 ) > f ( t x 1 + ( 1 − t ) x 2 ) t f(x_1) + (1-t)f(x_2) > f\Big(t x_1 + (1-t)x_2\Big) tf(x1)+(1−t)f(x2)>f(tx1+(1−t)x2),则称函数 f f f 为严格凸的。
凹函数 满足
t
f
(
x
1
)
+
(
1
−
t
)
f
(
x
2
)
≤
f
(
t
x
1
+
(
1
−
t
)
x
2
)
(1.2)
t f(x_1) + (1-t)f(x_2) \leq f\Big(t x_1 + (1-t)x_2\Big) \tag{1.2}
tf(x1)+(1−t)f(x2)≤f(tx1+(1−t)x2)(1.2)
同理,如果总有
t
f
(
x
1
)
+
(
1
−
t
)
f
(
x
2
)
<
f
(
t
x
1
+
(
1
−
t
)
x
2
)
t f(x_1) + (1-t)f(x_2) < f\Big(t x_1 + (1-t)x_2\Big)
tf(x1)+(1−t)f(x2)<f(tx1+(1−t)x2),则称函数
f
f
f 为严格凹的。
从导数的角度理解,若在其定义域 C C C 上,函数 f f f 二次可微,如果 f ′ ′ ( x ) ≥ 0 f^{''}(x) \geq 0 f′′(x)≥0,那么 f f f 为凸函数;反之,如果 f ′ ′ ( x ) ≤ 0 f^{''}(x) \leq 0 f′′(x)≤0,那么 f f f 为凹函数。
![](https://i-blog.csdnimg.cn/blog_migrate/e54d22933ff90a23f8f5b44835e41726.png)
1.1.2 Jensen 不等式
1.1.2.1 Jensen 不等式
对于任意点集 { x i } \{x_i\} {xi},若 λ i ≥ 0 \lambda_i \geq 0 λi≥0 且 ∑ i λ i = 1 \sum_{i}\lambda_i = 1 ∑iλi=1,函数 f ( x ) f(x) f(x) 满足:
f ( ∑ i = 1 m λ i x i ) ≤ ∑ i = 1 m λ i f ( x i ) , 如 果 f ( x ) 为 凸 函 数 f ( ∑ i = 1 m λ i x i ) ≥ ∑ i = 1 m λ i f ( x i ) , 如 果 f ( x ) 为 凹 函 数 (1.3) \begin{aligned} f(\sum_{i=1}^{m} \lambda_i x_i) \leq \sum_{i=1}^{m} \lambda_i f(x_i), \quad 如果 f(x) 为凸函数 \\ f(\sum_{i=1}^{m} \lambda_i x_i) \geq \sum_{i=1}^{m} \lambda_i f(x_i), \quad 如果 f(x) 为凹函数 \end{aligned} \tag{1.3} f(i=1∑mλixi)≤i=1∑mλif(xi),如果f(x)为凸函数f(i=1∑mλixi)≥i=1∑mλif(xi),如果f(x)为凹函数(1.3)
上式被称为 Jensen 不等式。
在概率论中,如果把
λ
i
\lambda_i
λi 看成取值为
x
i
x_i
xi 的离散变量
X
X
X 的概率分布,那么 Jensen 不等式可以写成
f
(
E
[
X
]
)
≤
E
[
f
(
X
)
]
,
如
果
f
(
x
)
为
凸
函
数
f
(
E
[
X
]
)
≥
E
[
f
(
X
)
]
,
如
果
f
(
x
)
为
凹
函
数
(1.4)
\begin{aligned} f(E[X]) \leq E[f(X)], \quad 如果 f(x) 为凸函数 \\ f(E[X]) \geq E[f(X)], \quad 如果 f(x) 为凹函数 \end{aligned} \tag{1.4}
f(E[X])≤E[f(X)],如果f(x)为凸函数f(E[X])≥E[f(X)],如果f(x)为凹函数(1.4)
其中,
E
[
⋅
]
E[\cdot]
E[⋅] 表示期望。此外,如果
f
(
x
)
f(x)
f(x) 为严格凸函数或者凹函数,当且仅当
X
=
E
[
X
]
X=E[X]
X=E[X]时,
E
[
f
(
X
)
]
=
f
(
E
[
X
]
)
E[f(X)] = f(E[X])
E[f(X)]=f(E[X]) 成立(即当
X
X
X 为一个常数时)。
对于连续变量,Jensen 不等式给出了积分的函数值和函数的积分值间的关系:
f ( ∫ x p ( x ) d x ) ≤ ∫ f ( x ) p ( x ) d x , 如 果 f ( x ) 为 凸 函 数 f ( ∫ x p ( x ) d x ) ≥ ∫ f ( x ) p ( x ) d x , 如 果 f ( x ) 为 凹 函 数 (1.5) \begin{aligned} f(\int x p(x)dx) \leq \int f(x)p(x)dx, \quad 如果 f(x) 为凸函数 \\ f(\int x p(x)dx) \geq \int f(x)p(x)dx, \quad 如果 f(x) 为凹函数 \end{aligned} \tag{1.5} f(∫xp(x)dx)≤∫f(x)p(x)dx,如果f(x)为凸函数f(∫xp(x)dx)≥∫f(x)p(x)dx,如果f(x)为凹函数(1.5)
1.1.2.2 Jensen 不等式的证明过程
可以使用数学归纳法证明 Jensen 不等式成立,下面我们来证明当 f ( x ) f(x) f(x) 为凸函数的情况,为凹函数时可以采用同样的方法进行证明。
当 i = 1 , 2 i=1,2 i=1,2 时,由凸函数的定义式(1.1),可知其成立。
假设当 i = m i=m i=m 时,其成立;当 i = m + 1 i = m+1 i=m+1 时,得
f ( ∑ i = 1 m + 1 λ i x i ) = f ( λ m + 1 x m + 1 + ∑ i = 1 m λ i x i ) (1.6) f(\sum_{i=1}^{m+1} \lambda_i x_i) =f(\lambda_{m+1}x_{m+1} + \sum_{i=1}^{m}\lambda_i x_i) \tag{1.6} f(i=1∑m+1λixi)=f(λm+1xm+1+i=1∑mλixi)(1.6)
我们令 η i = λ i 1 − λ m + 1 \eta_i = \frac{\lambda_i}{1-\lambda_{m+1}} ηi=1−λm+1λi,得
f ( ∑ i = 1 m + 1 λ i x i ) = f ( λ m + 1 x m + 1 + ( 1 − λ m + 1 ) ∑ i = 1 m η i x i ) (1.7) f(\sum_{i=1}^{m+1} \lambda_i x_i) =f(\lambda_{m+1}x_{m+1} + (1-\lambda_{m+1})\sum_{i=1}^{m}\eta_i x_i) \tag{1.7} f(i=1∑m+1λixi)=f(λm+1xm+1+(1−λm+1)i=1∑mηixi)(1.7)
由凸函数的定义得
f ( ∑ i = 1 m + 1 λ i x i ) ≤ λ m + 1 f ( x m + 1 ) + ( 1 − λ m + 1 ) f ( ∑ i = 1 m η i x i ) (1.8) f(\sum_{i=1}^{m+1} \lambda_i x_i) \leq \lambda_{m+1} f(x_{m+1}) +(1-\lambda_{m+1})f(\sum_{i=1}^{m}\eta_i x_i) \tag{1.8} f(i=1∑m+1λixi)≤λm+1f(xm+1)+(1−λm+1)f(i=1∑mηixi)(1.8)
因为 ∑ i = 1 m + 1 λ i = 1 \sum_{i=1}^{m+1} \lambda_i = 1 ∑i=1m+1λi=1,可得 ∑ i = 1 m λ i = 1 − λ m + 1 \sum_{i=1}^{m} \lambda_i = 1-\lambda_{m+1} ∑i=1mλi=1−λm+1,所以
∑ i = 1 m η i = ∑ i = 1 m λ i 1 − λ m + 1 = 1 (1.9) \sum_{i=1}^{m}\eta_i = \frac{\sum_{i=1}^{m} \lambda_i}{1-\lambda_{m+1}} = 1 \tag{1.9} i=1∑mηi=1−λm+1∑i=1mλi=1(1.9)
由式(1.3)和(1.9)得
f ( ∑ i = 1 m η i x i ) ≤ ∑ i = 1 m η i f ( x i ) (1.10) f(\sum_{i=1}^{m}\eta_i x_i) \leq \sum_{i=1}^{m}\eta_i f(x_i) \tag{1.10} f(i=1∑mηixi)≤i=1∑mηif(xi)(1.10)
由式(1.8)和(1.10)得
f ( ∑ i = 1 m + 1 λ i x i ) ≤ λ m + 1 f ( x m + 1 ) + ( 1 − λ m + 1 ) ∑ i = 1 m η i f ( x i ) = λ m + 1 f ( x m + 1 ) + ∑ i = 1 m λ i f ( x i ) = ∑ i = 1 m + 1 λ i f ( x i ) (1.11) \begin{aligned} f(\sum_{i=1}^{m+1} \lambda_i x_i) &\leq \lambda_{m+1} f(x_{m+1}) +(1-\lambda_{m+1})\sum_{i=1}^{m}\eta_i f(x_i) \\ &=\lambda_{m+1} f(x_{m+1}) +\sum_{i=1}^{m}\lambda_i f(x_i) \\ &=\sum_{i=1}^{m+1} \lambda_i f(x_i) \end{aligned} \tag{1.11} f(i=1∑m+1λixi)≤λm+1f(xm+1)+(1−λm+1)i=1∑mηif(xi)=λm+1f(xm+1)+i=1∑mλif(xi)=i=1∑m+1λif(xi)(1.11)
因此当 i = m + 1 i=m+1 i=m+1 时,Jensen 不等式也成立。
综上,Jensen 不等式成立。
1.2 EM 算法
阅读下面的内容,需要您对似然函数和极大似然估计有一定了解,如果您还不了解或者想温习一下,可以参考 这里。
概率模型有时既含有观测变量(observable variable),又含有隐变量或潜在变量(latent variable)。如果概率模型的变量都是观测变量,那么给定数据,可以直接用极大似然估计法,或贝叶斯估计法估计模型参数。但是,当模型含有隐变量时,就不能简单地使用这些估计方法。EM 算法就是含有隐变量的概率模型参数的极大似然估计法,或极大后验概率估计法。
下面我们从 Chuong B Do & Serafim Batzoglou 在论文《What is the expectation maximization algorithm?》中的抛硬币的例子开始。论文中的截图如下:
![](https://img-blog.csdnimg.cn/20200308140604776.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L2JlbnpodWppZTEyNDVjb20=,size_16,color_FFFFFF,t_70)
我们将论文中的实验简单整理后如下:
场景:假设有两枚硬币(硬币 A A A 和硬币 B B B),硬币 A A A 正面朝上的概率为 θ A \theta_A θA,硬币 B B B 正面朝上的概率为 θ B \theta_B θB。
实验:从两枚硬币中随机选择 1 枚硬币(选中的概率是相等的),使用选中的硬币抛十次。重复 5 次这样的操作,最终将得到 50 次的抛硬币结果。
目的:根据实验结果,推断出 θ A \theta_A θA 和 θ B \theta_B θB 的值。
如果我们知道每次抛的硬币是硬币 A A A 或者硬币 B B B,那我们可以直接使用极大似然估计求解 θ A \theta_A θA 和 θ B \theta_B θB 的值。但是当前我们不知道每次抛的到底是哪个硬币,即存在隐变量,这个时候我们需要使用 EM 算法求解 θ A \theta_A θA 和 θ B \theta_B θB 的值。
EM 算法流程:
(1)参数初始化:假设 θ ^ A ( 0 ) = 0.60 \hat{\theta}_A^{(0)} = 0.60 θ^A(0)=0.60, θ ^ B ( 0 ) = 0.50 \hat{\theta}_B^{(0)} = 0.50 θ^B(0)=0.50。
(2)E 步:计算后验概率,即每一次抛的硬币分别为 A A A 和 B B B 的概率。
比如第一次抛硬币的结果为 HTTTHHTHTH,即 5 正 5 反,该硬币分别为
A
A
A 和
B
B
B 的概率为
P
(
A
)
=
θ
A
5
(
1
−
θ
A
)
5
θ
A
5
(
1
−
θ
A
)
5
+
θ
B
5
(
1
−
θ
B
)
5
=
0.
6
5
∗
0.
4
5
0.
6
5
∗
0.
4
5
+
0.
5
5
∗
0.
5
5
≈
0.45
P(A) = \frac{\theta_A^5 (1-\theta_A)^5} {\theta_A^5(1-\theta_A)^5 + \theta_B^5 (1-\theta_B)^5} =\frac{0.6^5 * 0.4^5}{0.6^5 * 0.4^5 + 0.5^5*0.5^5} \approx 0.45
P(A)=θA5(1−θA)5+θB5(1−θB)5θA5(1−θA)5=0.65∗0.45+0.55∗0.550.65∗0.45≈0.45
P ( B ) = 1 − P ( A ) ≈ 0.55 P(B) = 1 - P(A) \approx 0.55 P(B)=1−P(A)≈0.55
使用上面的计算方法,可以求出 5 次实验中,硬币分别为 A A A 和 B B B 的概率如表 1.1 所示。
次数 | 为硬币 A A A 的概率 | 为硬币 B B B 的概率 |
---|---|---|
1 | 0.45 | 0.55 |
2 | 0.80 | 0.20 |
3 | 0.73 | 0.27 |
4 | 0.35 | 0.65 |
5 | 0.65 | 0.35 |
然后根据上面得出的概率,计算硬币 A A A 和硬币 B B B 出现正反面的期望。比如第一次抛硬币为 5 正 5 反,则硬币 A A A 分别为正面和反面的期望为
E ( A H ) = P ( A ) ∗ 5 = 0.45 ∗ 5 ≈ 2.2 E(A_{H}) = P(A)*5 = 0.45 * 5 \approx 2.2 E(AH)=P(A)∗5=0.45∗5≈2.2
E ( A T ) = P ( A ) ∗ 5 = 0.45 ∗ 5 ≈ 2.2 E(A_{T}) = P(A)*5 = 0.45 * 5 \approx 2.2 E(AT)=P(A)∗5=0.45∗5≈2.2
使用上面的计算方法,可以求出 5 次实验中,硬币 A A A 和硬币 B B B 出现正反面的期望如表 1.2 所示。
次数 | 硬币 A A A | 硬币 B B B |
---|---|---|
1 | ≈ \approx ≈ 2.2 H, 2.2 T | ≈ \approx ≈ 2.8 H, 2.8 T |
2 | ≈ \approx ≈ 7.2 H, 0.8 T | ≈ \approx ≈ 1.8 H, 0.2 T |
3 | ≈ \approx ≈ 5.9 H, 1.5 T | ≈ \approx ≈ 2.1 H, 0.5 T |
4 | ≈ \approx ≈ 1.4 H, 2.1 T | ≈ \approx ≈ 2.6 H, 3.9 T |
5 | ≈ \approx ≈ 4.5 H, 1.9 T | ≈ \approx ≈ 2.5 H, 1.1 T |
合计 | ≈ \approx ≈ 21.3 H, 8.6 T | ≈ \approx ≈ 11.7 H, 8.4 T |
(3)M 步:计算新的参数 θ ^ A \hat{\theta}_A θ^A 和 θ ^ B \hat{\theta}_B θ^B
θ ^ A ( 1 ) ≈ 21.3 21.3 + 8.6 ≈ 0.71 \hat{\theta}_A^{(1)} \approx \frac{21.3}{21.3 + 8.6} \approx 0.71 θ^A(1)≈21.3+8.621.3≈0.71
θ ^ B ( 1 ) ≈ 11.7 11.7 + 8.4 ≈ 0.58 \hat{\theta}_B^{(1)} \approx \frac{11.7}{11.7 + 8.4} \approx 0.58 θ^B(1)≈11.7+8.411.7≈0.58
(4)进行迭代:重复 E 步和 M 步,直到收敛。
经过十次迭代后,得到
θ
^
A
(
10
)
≈
0.80
\hat{\theta}_A^{(10)} \approx 0.80
θ^A(10)≈0.80,
θ
^
B
(
10
)
≈
0.52
\hat{\theta}_B^{(10)} \approx 0.52
θ^B(10)≈0.52。
1.3 EM 算法的推导
通过上面抛硬币的例子,我们已经大致了解了 EM 算法,下面我们开始详细介绍 EM 算法的推导过程。
给定 m m m 个训练样本 { x ( 1 ) , ⋯ , x ( m ) } \{x^{(1)}, \cdots, x^{(m)}\} {x(1),⋯,x(m)},假设样本间相互独立,我们希望将模型 p ( x , z ) p(x,z) p(x,z) 的参数与数据进行拟合,其似然函数为:
l ( θ ) = ∑ i = 1 m log p ( x ; θ ) = ∑ i = 1 m log ∑ z p ( x , z ; θ ) (1.12) \begin{aligned} l(\theta) &= \sum_{i=1}^{m} \log p(x;\theta) \\ &= \sum_{i=1}^{m} \log \sum_{z} p(x,z;\theta) \end{aligned} \tag{1.12} l(θ)=i=1∑mlogp(x;θ)=i=1∑mlogz∑p(x,z;θ)(1.12)
但是,直接求解参数 θ \theta θ 的极大似然估计一般会比较困难,因为上式存在一个隐变量 z z z。通常情况下,如果确定 z z z 后,求解 θ \theta θ 就很容易了。
针对存在含有隐变量的情况下,EM 算法提供了一种有效的极大似然估计方法。因为无法直接最大化 l ( θ ) l(\theta) l(θ),所以采用此方法:不断地建立 l ( θ ) l(\theta) l(θ) 的下界(E步),然后优化下界(M步)。这句话比较抽象,我们继续往下看。
对每一个样例 i i i,让 Q i Q_i Qi 表示表示该样例隐变量 z z z 的某种分布(存在 ∑ z Q i ( z ) = 1 \sum_{z} Q_i(z) =1 ∑zQi(z)=1, Q i ( z ) ≥ 0 Q_i(z) \geq 0 Qi(z)≥0)需要注意:如果 Q i Q_i Qi 是连续性的,则 Q i Q_i Qi 表示概率密度函数,需要将求和符号换成积分符号。
对式(1.12)进行变换得:
∑ i log p ( x ( i ) ; θ ) = ∑ i log ∑ z ( i ) p ( x ( i ) , z ( i ) ; θ ) = ∑ i log ∑ z ( i ) Q i ( z ( i ) ) p ( x ( i ) , z ( i ) ; θ ) Q i ( z ( i ) ) ≥ ∑ i ∑ z ( i ) Q i ( z ( i ) ) log p ( x ( i ) , z ( i ) ; θ ) Q i ( z ( i ) ) (1.13) \begin{aligned} \sum_{i} \log p(x^{(i)};\theta) &= \sum_{i} \log \sum_{z^{(i)}} p(x^{(i)},z^{(i)};\theta) \\ &= \sum_{i} \log \sum_{z^{(i)}} Q_i(z^{(i)}) \frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \\ &\geq \sum_{i} \sum_{z^{(i)}} Q_i(z^{(i)}) \log \frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \end{aligned} \tag{1.13} i∑logp(x(i);θ)=i∑logz(i)∑p(x(i),z(i);θ)=i∑logz(i)∑Qi(z(i))Qi(z(i))p(x(i),z(i);θ)≥i∑z(i)∑Qi(z(i))logQi(z(i))p(x(i),z(i);θ)(1.13)
对分子和分母同时乘以 Q i ( z ( i ) ) Q_{i}(z^{(i)}) Qi(z(i)),所以式(1.13)中的第二个等号成立。根据 Jensen 不等式,式(1.13)中的不等式成立。在这里, f ( x ) = log ( x ) f(x) = \log(x) f(x)=log(x),由于 log ( x ) \log(x) log(x) 的二阶导数为 − 1 x 2 < 0 -\frac{1}{x^2} < 0 −x21<0,所以其为凹函数。我们可以把 Q i ( z ( i ) ) Q_i(z^{(i)}) Qi(z(i)) 看做概率分布 p p p,那么 ∑ z ( i ) Q i ( z ( i ) ) p ( x ( i ) , z ( i ) ; θ ) Q i ( z ( i ) ) \sum_{z^{(i)}} Q_i(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} ∑z(i)Qi(z(i))Qi(z(i))p(x(i),z(i);θ) 可以看做是 p ( x ( i ) , z ( i ) ; θ ) Q i ( z ( i ) ) \frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} Qi(z(i))p(x(i),z(i);θ) 的期望。根据 Jensen 不等式可得
f ( E z ( i ) ∼ Q i [ p ( x ( i ) , z ( i ) ; θ ) Q i ( z ( i ) ) ] ) ≥ E z ( i ) ∼ Q i [ f ( p ( x ( i ) , z ( i ) ; θ ) Q i ( z ( i ) ) ) ] (1.14) f\left(E_{z^{(i)}\thicksim Q_i} \left[\frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})}\right]\right) \geq E_{z^{(i)}\thicksim Q_i}\left[f\left( \frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \right)\right] \tag{1.14} f(Ez(i)∼Qi[Qi(z(i))p(x(i),z(i);θ)])≥Ez(i)∼Qi[f(Qi(z(i))p(x(i),z(i);θ))](1.14)
这样就得到式(1.13)中的不等式成立。
我们可以把上式写成: L ( θ ) ≥ J ( z , Q ) L(\theta) \geq J(z, Q) L(θ)≥J(z,Q)( z z z 为隐变量),那么我们可以通过不断的最大化 J ( z , Q ) J(z,Q) J(z,Q),来使得 L ( θ ) L(\theta) L(θ) 不断提高,最终达到它的最大值。图 1.2 更形象地描述这个过程:
![](https://i-blog.csdnimg.cn/blog_migrate/a59daf24c440024d011fbecaabe50cd4.png)
这里来说说上图的内在含义:首先我们固定 θ \theta θ,调整 Q ( z ) Q(z) Q(z) 使下界 J ( z , Q ) J(z,Q) J(z,Q) 与 L ( θ ) L(\theta) L(θ) 在此点 θ \theta θ 处相等(即绿色曲线到蓝色曲线),然后固定 Q ( z ) Q(z) Q(z),调整 θ \theta θ 使下界 J ( z , Q ) J(z,Q) J(z,Q) 达到最大值(即 θ ( t ) \theta^{(t)} θ(t) 到 θ ( t + 1 ) \theta^{(t+1)} θ(t+1));然后再固定 θ \theta θ,调整 Q ( z ) Q(z) Q(z),……,直到收敛到 L ( θ ) L(\theta) L(θ) 的最大值处 θ ∗ \theta^{*} θ∗。
在上面的迭代过程中,存在以下两个问题:
- 什么时候下界 J ( z , Q ) J(z, Q) J(z,Q) 与 L ( θ ) L(\theta) L(θ) 在此点 θ \theta θ 处相等?
- 为什么一定会收敛?
1.3.1 什么时候下界 J ( z , Q ) J(z, Q) J(z,Q) 与 L ( θ ) L(\theta) L(θ) 在此点 θ \theta θ 处相等?
在前面介绍 Jensen 不等式时提到,当自变量 X = E ( x ) X=E(x) X=E(x) 时,即为常数的时候,等式成立。换言之,为了使式(1.13)中的不等式取等号,需要满足
p ( x ( i ) , z ( i ) ; θ ) Q i ( z ( i ) ) = c (1.15) \frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} = c \tag{1.15} Qi(z(i))p(x(i),z(i);θ)=c(1.15)
其中, c c c 为常数,不依赖于 z ( i ) z^{(i)} z(i)。对上面的等式做一下变换,得
p ( x ( i ) , z ( i ) ; θ ) = c Q i ( z ( i ) ) (1.16) p(x^{(i)},z^{(i)};\theta) = cQ_i(z^{(i)}) \tag{1.16} p(x(i),z(i);θ)=cQi(z(i))(1.16)
对上面的等式两边对 z z z 求和,得
∑ z p ( x ( i ) , z ( i ) ; θ ) = ∑ z c Q i ( z ( i ) ) = c ∑ Q i ( z ( i ) ) (1.17) \begin{aligned} \sum_z p(x^{(i)},z^{(i)};\theta) &= \sum_z cQ_i(z^{(i)}) \\ &= c \sum Q_i(z^{(i)}) \end{aligned} \tag{1.17} z∑p(x(i),z(i);θ)=z∑cQi(z(i))=c∑Qi(z(i))(1.17)
因为 ∑ Q i ( z ( i ) ) = 1 \sum Q_i(z^{(i)}) = 1 ∑Qi(z(i))=1(概率之和为 1),得
∑ z p ( x ( i ) , z ( i ) ; θ ) = c (1.18) \sum_z p(x^{(i)},z^{(i)};\theta) = c \tag{1.18} z∑p(x(i),z(i);θ)=c(1.18)
由式(1.15)和式(1.18)得
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 ) ; θ ) (1.19) \begin{aligned} Q_i(z^{(i)}) &= \frac{p(x^{(i)},z^{(i)};\theta)}{c} \\ &= \frac{p(x^{(i)},z^{(i)};\theta)}{\sum_z p(x^{(i)},z^{(i)};\theta)} \\ &= \frac{p(x^{(i)},z^{(i)};\theta)}{p(x^{(i)};\theta)} \\ &= p(z^{(i)} | x^{(i)}; \theta) \end{aligned} \tag{1.19} 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);θ)(1.19)
至此,我们推出了在固定参数 θ \theta θ 后, J ( z , Q ) J(z,Q) J(z,Q) 与 L ( θ ) L(\theta) L(θ) 相等时, Q i ( z ( i ) ) Q_{i}(z^{(i)}) Qi(z(i)) 的取值就是其后验概率(在给定 x ( i ) x^{(i)} x(i) 和 θ \theta θ 后),这样我们同时解决了 Q i ( z ( i ) ) Q_{i}(z^{(i)}) Qi(z(i)) 的选值问题。此步就是 E 步,即建立 L ( θ ) L(\theta) L(θ) 的下界。接下来是 M 步,即在给定 Q i ( z ( i ) ) Q_i(z^{(i)}) Qi(z(i)) 后,调整 θ \theta θ,从而极大化 L ( θ ) L(\theta) L(θ) 的下界 J ( z , Q ) J(z,Q) J(z,Q)。不断地重复 E 步和 M 步,直至收敛,这就是 EM 算法。
EM 算法的完整步骤如下:
- 参数初始化:随机初始化参数 θ ( 0 ) \theta^{(0)} θ(0)。
- E 步:根据当前参数 θ ( t ) \theta^{(t)} θ(t)(初始值 θ ( 0 ) \theta^{(0)} θ(0) 或上一次迭代中 M 步求得的 θ \theta θ 值)求隐变量的后验概率 Q i ( z ( i ) ) = p ( z ( i ) ∣ x ( i ) ; θ ( t ) ) Q_i(z^{(i)})=p(z^{(i)} | x^{(i)}; \theta^{(t)}) Qi(z(i))=p(z(i)∣x(i);θ(t)),即式(1.19)。
- M 步:固定 Q i ( z ( i ) ) Q_i(z^{(i)}) Qi(z(i)),通过极大化 J ( z , Q ) J(z,Q) J(z,Q) 求得新的参数 θ ( t ) \theta^{(t)} θ(t)。
- 进行迭代:重复 E 步和 M 步,直到收敛。
1.3.2 EM 算法的收敛性
我们怎么确保 EM 算法一定会收敛呢?首先,假设 θ ( t ) \theta^{(t)} θ(t) 和 θ ( t + 1 ) \theta^{(t+1)} θ(t+1) 是 EM 算法第 t t t 次和 t + 1 t+1 t+1 次迭代后的结果。如果我们证明了 l ( θ ( t ) ) ≤ l ( θ ( t + 1 ) ) l(\theta^{(t)}) \leq l(\theta^{(t+1)}) l(θ(t))≤l(θ(t+1)),也就是说对数似然函数单调递增,那么最终就会得到最大值。
证明过程:
在选定
θ
(
t
)
\theta^{(t)}
θ(t) 后,由 E 步得
Q
i
(
t
)
(
z
(
i
)
)
:
=
p
(
z
(
i
)
∣
x
(
i
)
;
θ
(
t
)
)
(1.20)
Q_i^{(t)}(z^{(i)}) := p(z^{(i)}|x^{(i)};\theta^{(t)}) \tag{1.20}
Qi(t)(z(i)):=p(z(i)∣x(i);θ(t))(1.20)
当
Q
i
Q_i
Qi 为后验概率时保证了 Jensen 不等式中的等号成立,即得
l ( θ ( t ) ) = ∑ i ∑ z ( i ) Q i ( t ) ( z ( i ) ) log p ( x ( i ) , z ( i ) ; θ ( t ) ) Q i ( t ) ( z ( i ) ) (1.21) l(\theta^{(t)}) = \sum_i \sum_{z^{(i)}} Q_i^{(t)}(z^{(i)}) \log \frac{p(x^{(i)},z^{(i)};\theta^{(t)})}{Q_i^{(t)}(z^{(i)})} \tag{1.21} l(θ(t))=i∑z(i)∑Qi(t)(z(i))logQi(t)(z(i))p(x(i),z(i);θ(t))(1.21)
参数 θ ( t + 1 ) \theta^{(t+1)} θ(t+1) 是通过极大化式(1.21)得到。经过一些推导会有一下式子成立
l ( θ ( t + 1 ) ) ≥ ∑ i ∑ z ( i ) Q i ( t ) ( z ( i ) ) log p ( x ( i ) , z ( i ) ; θ ( t + 1 ) ) Q i ( t ) ( z ( i ) ) ≥ ∑ i ∑ z ( i ) Q i ( t ) ( z ( i ) ) log p ( x ( i ) , z ( i ) ; θ ( t ) ) Q i ( t ) ( z ( i ) ) = l ( θ ( t ) ) (1.22) \begin{aligned} l(\theta^{(t+1)}) &\geq \sum_i \sum_{z^{(i)}} Q_i^{(t)}(z^{(i)}) \log \frac{p(x^{(i)},z^{(i)};\theta^{(t+1)})}{Q_i^{(t)}(z^{(i)})} \\ &\geq \sum_i \sum_{z^{(i)}} Q_i^{(t)}(z^{(i)}) \log \frac{p(x^{(i)},z^{(i)};\theta^{(t)})}{Q_i^{(t)}(z^{(i)})} \\ &= l(\theta^{(t)}) \end{aligned} \tag{1.22} l(θ(t+1))≥i∑z(i)∑Qi(t)(z(i))logQi(t)(z(i))p(x(i),z(i);θ(t+1))≥i∑z(i)∑Qi(t)(z(i))logQi(t)(z(i))p(x(i),z(i);θ(t))=l(θ(t))(1.22)
下面来具体解释一下式(1.22)。对下面的式子(即式(1.13))
l ( θ ) ≥ ∑ i ∑ z ( i ) Q i ( z ( i ) ) log p ( x ( i ) , z ( i ) ; θ ) Q i ( z ( i ) ) (1.23) l(\theta) \geq \sum_{i} \sum_{z^{(i)}} Q_i(z^{(i)}) \log \frac{p(x^{(i)},z^{(i)};\theta)}{Q_i(z^{(i)})} \tag{1.23} l(θ)≥i∑z(i)∑Qi(z(i))logQi(z(i))p(x(i),z(i);θ)(1.23)
中的参数分别取值 Q i = Q i ( t ) Q_i = Q_i^{(t)} Qi=Qi(t) 和 θ = θ ( t + 1 ) \theta = \theta^{(t+1)} θ=θ(t+1) 得
l ( θ ( t + 1 ) ) ≥ ∑ i ∑ z ( i ) Q i ( t ) ( z ( i ) ) log p ( x ( i ) , z ( i ) ; θ ( t + 1 ) ) Q i ( t ) ( z ( i ) ) (1.24) l(\theta^{(t+1)}) \geq \sum_{i} \sum_{z^{(i)}} Q_i^{(t)}(z^{(i)}) \log \frac{p(x^{(i)},z^{(i)};\theta^{(t+1)})}{Q_i^{(t)}(z^{(i)})} \tag{1.24} l(θ(t+1))≥i∑z(i)∑Qi(t)(z(i))logQi(t)(z(i))p(x(i),z(i);θ(t+1))(1.24)
从而式(1.22)中的第一个不等号成立。
因为参数 θ ( t + 1 ) \theta^{(t+1)} θ(t+1) 是通过极大化式(1.21)得到,所以可得
∑ i ∑ z ( i ) Q i ( t ) ( z ( i ) ) log p ( x ( i ) , z ( i ) ; θ ( t + 1 ) ) Q i ( t ) ( z ( i ) ) ≥ ∑ i ∑ z ( i ) Q i ( t ) ( z ( i ) ) log p ( x ( i ) , z ( i ) ; θ ( t ) ) Q i ( t ) ( z ( i ) ) (1.25) \sum_{i} \sum_{z^{(i)}} Q_i^{(t)}(z^{(i)}) \log \frac{p(x^{(i)},z^{(i)};\theta^{(t+1)})}{Q_i^{(t)}(z^{(i)})} \geq\sum_i \sum_{z^{(i)}} Q_i^{(t)}(z^{(i)}) \log \frac{p(x^{(i)},z^{(i)};\theta^{(t)})}{Q_i^{(t)}(z^{(i)})} \tag{1.25} i∑z(i)∑Qi(t)(z(i))logQi(t)(z(i))p(x(i),z(i);θ(t+1))≥i∑z(i)∑Qi(t)(z(i))logQi(t)(z(i))p(x(i),z(i);θ(t))(1.25)
从而式(1.22)中的第二个不等号成立。
由式(1.21)可得式(1.22)中的等号成立。
至此,我们证明了 l ( θ ( t ) ) ≤ l ( θ ( t + 1 ) ) l(\theta^{(t)}) \leq l(\theta^{(t+1)}) l(θ(t))≤l(θ(t+1))。这样也就证明了 EM 算使得似然函数 l ( θ ) l(\theta) l(θ) 单调递增,并收敛。在实际应用中,可以使用如下方法判断是否收敛: l ( θ ) l(\theta) l(θ) 不再变化或者变化很小。
从上面的推导可以看出,EM 算法可以保证收敛到一个稳定点,但是却不能保证收敛到全局的极大值点,因此它是局部最优的算法。所以在应用中,初值的选择变得非常重要,常用的办法是选取几个不同的初值进行迭代,然后对得到的各个估计值加以比较,从中选择最好的。
从前面的推导中我们知道
l
(
θ
)
≥
J
(
Q
,
θ
)
l(\theta) \geq J(Q,\theta)
l(θ)≥J(Q,θ),EM 算法可以看做是
J
(
Q
,
θ
)
J(Q,\theta)
J(Q,θ) 的坐标上升法。E 步固定
θ
\theta
θ,优化
Q
Q
Q;M 步固定
Q
Q
Q,优化
θ
\theta
θ。
1.4 EM 算法的应用
K-means 算法是 EM 算法思想的体现,E 步骤为聚类过程,M 步骤为更新类簇中心。高斯混合模型(GMM)也是 EM 算法的一个应用。
缺点:对初始值敏感,EM 算法需要初始化参数
θ
\theta
θ,而参数
θ
\theta
θ 的选择直接影响收敛效率以及能否得到全局最优解。
参考
[1] 李航《统计学习方法》
[2] 周志华《机器学习》
[3] http://cs229.stanford.edu/notes/cs229-notes8.pdf
[4] Chuong B Do & Serafim Batzoglou《What is the expectation maximization algorithm?》
[5] 维基百科:凸函数
[6] 维基百科:凹函数
[7] Jensen不等式初步理解及证明
[8] 如何通俗理解EM算法
[9] 从最大似然到EM算法浅解
[10] EM-最大期望算法
[11] (EM算法)The EM Algorithm