EM演算法及其應用（一）

07-29

EM演算法及其應用（一）

來自專欄數據科學の雜談

$quad$

EM演算法是期望最大化 (Expectation Maximization) 演算法的簡稱，用於含有隱變數的情況下，概率模型參數的極大似然估計或極大後驗估計。EM演算法是一種迭代演算法，每次迭代由兩步組成：E步，求期望 (expectation)，即利用當前估計的參數值來計算對數似然函數的期望值；M步，求極大 (maximization)，即求參數 $heta$ 來極大化E步中的期望值，而求出的參數 $heta$ 將繼續用於下一個E步中期望值的估計。EM演算法在機器學習中應用廣泛，本篇和下篇文章分別探討EM演算法的原理和其兩大應用 —— K-means和高斯混合模型。

================= 1、先驗知識 ==================

一、凸函數、凹函數和 Jensen不等式

設 $f(x)$ 為定義在區間 $I = [a,b]$ 上的實值函數，對於任意 $forall , x_1, x_2 in I, lambda in [0,1]$ ，有：

$qquadqquadqquad f(lambda ,x_1 + (1-lambda),x_2) leq lambda f(x_1) + (1-lambda)f(x_2)$

則 $f(x)$ 為凸函數 (convex function)，如下圖所示。相應的，若上式中 $leqslant$ 變為 $geqslant$ ，則 $f(x)$ 為凹函數 (concave function)。凸函數的判定條件是二階導 $f^{}(x) geqslant 0$ ，而凹函數為 $f^{}(x) leqslant 0$ 。後文要用到的對數函數$ln(x)$的二階導為 $-frac{1}{x^2} < 0$ ，所以是凹函數。

Jensen不等式就是上式的推廣，設 $f(x)$ 為凸函數， $lambda_i geqslant 0, ;; sum_i lambda_i = 1$ ，則：

$qquad qquad qquad qquad qquad fleft(sumlimits_{i=1}^n lambda_i x_i ight) leq sumlimits_{i=1}^n lambda_i f(x_i)$

如果是凹函數，則將不等號反向，若用對數函數來表示，就是：

$qquadqquadqquadqquadqquad lnleft(sumlimits_{i=1}^n lambda_i x_i ight) geq sumlimits_{i=1}^n lambda_i ln(x_i)$

若將 $lambda_i$ 視為一個概率分布，則可表示為期望值的形式，在後文中同樣會引入概率分布：

$qquadqquadqquadqquadqquadqquad f(mathbb{E}[mathrm{x}]) leq mathbb{E}[f(mathrm{x})]$

$quad$

二、KL散度

KL散度(Kullback-Leibler divergence) 又稱相對熵 (relative entropy)，主要用于衡量兩個概率分布p和q的差異，也可理解為兩個分布對數差的期望。

$mathbb{KL}(p||q) = sum_i p(x_i)log frac{p(x_i)}{q(x_i)}= mathbb{E}_{mathrm{x}sim p}left[log frac{p(x)}{q(x)} ight] = mathbb{E}_{mathrm{x}sim p}left[log,p(x) - log,q(x) ight ]$

KL散度總滿足 $mathbb{KL}(p||q) geqslant 0$ ，而當且僅當 $q=p$ 時， $mathbb{KL}(p||q) = 0$ 。一般來說分布 $p(x)$ 比較複雜，因而希望用比較簡單的 $q(x)$ 去近似 $p(x)$ ，而近似的標準就是KL散度越小越好。

KL散度不具備對稱性，即 $mathbb{KL}(p||q) eq mathbb{KL}(q||p)$ ，因此不能作為一個距離指標。

$quad$

三、極大似然估計和極大後驗估計

極大似然估計 (Maximum likelihood estimation) 是參數估計的常用方法，基本思想是在給定樣本集的情況下，求使得該樣本集出現的「可能性」最大的參數 $heta$ 。將參數 $heta$ 視為未知量，則參數 $heta$ 對於樣本集X的對數似然函數為：

$qquadqquadqquadqquadqquadqquad L( heta) = ln ,P(X| heta)$

這個函數反映了在觀測結果X已知的條件下， $heta$ 的各種值的「似然程度」。這裡是把觀測值X看成結果，把參數 $heta$ 看成是導致這個結果的原因。參數$ heta$雖然未知但是有著固定值 (當然這是頻率學派的觀點)，並非事件或隨機變數，無概率可言，因而改用「似然(likelihood)" 這個詞。

於是通過求導求解使得對數似然函數最大的參數 $heta$ ， $heta = mathop{argmax}limits_{ heta}L( heta)$ ，即為極大似然法。

極大後驗估計 (Maximum a posteriori estimation) 是貝葉斯學派的參數估計方法，相比於頻率學派，貝葉斯學派將參數 $heta$ 視為隨機變數，並將其先驗分布 $P( heta)$ 包含在估計過程中。運用貝葉斯定理，參數 $heta$ 的後驗分布為：

$qquadqquad qquad P( heta|X) = frac{P(X, heta)}{P(X)} = frac{P( heta)P(X| heta)}{P(X)} propto P( heta)P(X| heta)$

上式中 $P(X)$ 不依賴於 $heta$ 因而為常數項可以捨去，則最終結果為

$qquad qquad qquadqquadqquad heta = mathop{argmax}limits_{ heta}P( heta)P(X| heta)$

================ 2、EM演算法初探 =================

概率模型有時既含有觀測變數 (observable variable)，又含有隱變數 (hidden variable)，隱變數顧名思義就是無法被觀測到的變數。如果都是觀測變數，則給定數據，可以直接使用極大似然估計。但如果模型含有隱變數時，直接求導得到參數比較困難。而EM演算法就是解決此類問題的常用方法。

對於一個含有隱變數 $mathbf{Z}$ 的概率模型，一般將 ${mathbf{X}, mathbf{Z}}$ 稱為完全數據，而觀測數據 $mathbf{X}$ 為不完全數據。

我們的目標是極大化觀測數據 $mathbf{X}$ 關於參數 $oldsymbol{ heta}$ 的對數似然函數。由於存在隱變數，因而也可表示為極大化 $mathbf{X}$ 的邊緣分布 (marginal distribution)，即：

$L(oldsymbol{ heta}) = ln,P(mathbf{X}|oldsymbol{ heta}) = ln,sumlimits_{mathbf{Z}}P(mathbf{X},mathbf{Z}|oldsymbol{ heta}) ag{1.1}$

上式中存在「對數的和」—— $lnsum(cdot)$ ，如果直接求導將會非常困難。因而EM演算法採用曲線救國的策略，構建 $(1.1)$ 式的一個下界，然後通過極大化這個下界來間接達到極大化 $(1.1)$ 的效果。

要想構建下界，就需要運用上文中的Jensen不等式。記 $oldsymbol{ heta}^{(t)}$ 為第t步迭代參數的估計值，考慮引入一個分布 $P(mathbf{Z}|mathbf{X},oldsymbol{ heta}^{(t)})$ ，由於

$1. quad P(mathbf{Z}|mathbf{X},oldsymbol{ heta}^{(t)}) geqslant 0$

$2. quad sum_{mathbf{Z}}P(mathbf{Z}|mathbf{X},oldsymbol{ heta}^{(t)}) = 1$

$3. quad ln(cdot)$ 為凹函數

因而可以利用Jensen不等式求出 $L(oldsymbol{ heta})$ 的下界：

$egin{align} L(oldsymbol{ heta}) = ln,sumlimits_{mathbf{Z}}P(mathbf{X},mathbf{Z}|oldsymbol{ heta}) &= ln,sumlimits_{mathbf{Z}}P(mathbf{Z|mathbf{X},oldsymbol{ heta}^{(t)}})frac{P(mathbf{X},mathbf{Z}|oldsymbol{ heta}) }{P(mathbf{Z}|mathbf{X},oldsymbol{ heta}^{(t)})} ag{1.2}\ & geqslant sumlimits_{mathbf{Z}}P(mathbf{Z|mathbf{X},oldsymbol{ heta}^{(t)}}) ,lnfrac{P(mathbf{X},mathbf{Z}|oldsymbol{ heta}) }{P(mathbf{Z}|mathbf{X},oldsymbol{ heta}^{(t)})} ag{1.3} \ & = underbrace{sumlimits_{mathbf{Z}}P(mathbf{Z|mathbf{X},oldsymbol{ heta}^{(t)}}) ,ln{P(mathbf{X},mathbf{Z}|oldsymbol{ heta}) }}_{mathcal{Q}(oldsymbol{ heta},oldsymbol{ heta}^{(t)})} ;;underbrace{- sumlimits_{mathbf{Z}}P(mathbf{Z|mathbf{X},oldsymbol{ heta}^{(t)}}) ,ln{P(mathbf{Z}|mathbf{X},oldsymbol{ heta}^{(t)})}}_{entropy} ag{1.4} end{align}$

$(1.3)$ 式構成了 $L(oldsymbol{ heta})$ 的下界，而 $(1.4)$ 式的右邊為 $P(mathbf{Z}|mathbf{X},oldsymbol{ heta}^{(t)})$ 的熵 $geqslant 0$ ，其獨立於我們想要優化的參數 $oldsymbol{ heta}$ ，因而是一個常數。所以極大化 $L(oldsymbol{ heta})$ 的下界 $(1.3)$ 式就等價於極大化 $mathcal{Q}(oldsymbol{ heta}, oldsymbol{ heta}^{(t)})$ ， $mathcal{Q}(oldsymbol{ heta}, oldsymbol{ heta}^{(t)})$ (Q函數) 亦可表示為 $mathbb{E}_{,mathbf{Z}|mathbf{X},oldsymbol{ heta}^{(t)}},lnP(mathbf{X},mathbf{Z}|oldsymbol{ heta})$ ，其完整定義如下：

基於觀測數據 $mathbf{X}$ 和當前參數 $heta^{(t)}$ 計算未觀測數據 $mathbf{Z}$ 的條件概率分布 $P(mathbf{Z}|mathbf{X}, heta^{(t)})$ ，則Q函數為完全數據的對數似然函數關於 $mathbf{Z}$ 的期望。

此即E步中期望值的來歷。

接下來來看M步。在 $(1.3)$ 式中若令 $oldsymbol{ heta} = oldsymbol{ heta}^{(t)}$ ，則下界 $(1.3)$ 式變為：

$egin{align} & sumlimits_{mathbf{Z}}P(mathbf{Z|mathbf{X},oldsymbol{ heta}^{(t)}}) ,lnfrac{P(mathbf{X},mathbf{Z}|oldsymbol{ heta}^{(t)}) }{P(mathbf{Z}|mathbf{X},oldsymbol{ heta}^{(t)})} \ qquad qquadqquad =;; & sumlimits_{mathbf{Z}}P(mathbf{Z|mathbf{X},oldsymbol{ heta}^{(t)}}) ,lnfrac{P(mathbf{Z|mathbf{X},oldsymbol{ heta}^{(t)}})P(mathbf{X}|oldsymbol{ heta}^{(t)})}{P(mathbf{Z}|mathbf{X},oldsymbol{ heta}^{(t)})} \ = ;; & sumlimits_{mathbf{Z}}P(mathbf{Z|mathbf{X},oldsymbol{ heta}^{(t)}}) ,lnP(mathbf{X}|oldsymbol{ heta}^{(t)}) \ = ;; & lnP(mathbf{X}|oldsymbol{ heta}^{(t)}) ;;=;; L(oldsymbol{ heta}^{(t)}) end{align}$

可以看到在第t步， $L(oldsymbol{ heta}^{(t)})$ 的下界與 $L(oldsymbol{ heta}^{(t)})$ 相等，又由於極大化下界與極大化Q函數等價，因而在M步選擇一個新的 $oldsymbol{ heta}$ 來極大化 $mathcal{Q}(oldsymbol{ heta}, oldsymbol{ heta}^{(t)})$ ，就能使 $L(oldsymbol{ heta}) geqslant mathcal{Q}(oldsymbol{ heta}, oldsymbol{ heta}^{(t)}) geqslant mathcal{Q}(oldsymbol{ heta}^{(t)}, oldsymbol{ heta}^{(t)}) = L(oldsymbol{ heta}^{(t)})$ (這裡為了便於理解就將 $mathcal{Q}(oldsymbol{ heta}, oldsymbol{ heta}^{(t)})$ 與 $(1.3)$ 式等同了)，也就是說 $L(oldsymbol{ heta})$ 是單調遞增的，通過EM演算法的不斷迭代能保證收斂到局部最大值。

$quad$

EM演算法流程：

輸入：觀測數據 $mathbf{X}$ ，隱變數 $mathbf{Z}$ ，聯合概率分布 $P(mathbf{X},mathbf{Z}|oldsymbol{ heta})$

輸出：模型參數 $oldsymbol{ heta}$

1) 初始化參數 $oldsymbol{ heta}^{(0)}$

2) E步：利用當前參數 $oldsymbol{ heta}^{(t)}$ 計算Q函數，

$qquadqquadqquad mathcal{Q}(oldsymbol{ heta}, oldsymbol{ heta}^{(t)}) = sumlimits_{mathbf{Z}}P(mathbf{Z|mathbf{X},oldsymbol{ heta}^{(t)}}) ,ln{P(mathbf{X},mathbf{Z}|oldsymbol{ heta})}$

3) M步：極大化Q函數，求出相應的

$qquadqquadqquadqquadqquad oldsymbol{ heta} = mathop{argmax}limits_{oldsymbol{ heta}}mathcal{Q}(oldsymbol{ heta}, oldsymbol{ heta}^{(t)})$

4) 重複 2. 和3. 步直至收斂。

EM演算法也可用於極大後驗估計，極大後驗估計僅僅是在極大似然估計的基礎上加上參數 $oldsymbol{ heta}$ 的先驗分布，即 $p(oldsymbol{ heta})p(mathbf{X}|oldsymbol{ heta})$ ，則取對數後變為 $ln,p(mathbf{X}|oldsymbol{ heta}) + ln,p(oldsymbol{ heta})$ ，由於後面的 $ln,p(oldsymbol{ heta})$ 不包含隱變數 $mathbf{Z}$ ，所以E步中求Q函數的步驟不變。而在M步中需要求新的參數 $mathbf{ heta}$ ，因此需要包含這一項，所以M步變為

$qquadqquadqquadqquad oldsymbol{ heta} = mathop{argmax}limits_{oldsymbol{ heta}} left[mathcal{Q}(oldsymbol{ heta}, oldsymbol{ heta}^{(t)}) + ln(p(oldsymbol{ heta}) ight]$

================ 3、EM演算法深入 =================

上一節中遺留了一個問題：為什麼式 $(1.2)$ 中引入的分布是 $P(mathbf{Z}|mathbf{X},oldsymbol{ heta}^{(t)})$ 而不是其他分布？下面以另一個角度來闡述。

假設一個關於隱變數 $mathbf{Z}$ 的任意分布 $q(mathbf{Z})$ ，則運用期望值的定義， $(1.1)$ 式變為：

$egin{align*} L(oldsymbol{ heta}) = lnP(mathbf{X}|oldsymbol{ heta}) &= sumlimits_{mathbf{Z}}q(mathbf{Z}),lnP(mathbf{X}|oldsymbol{ heta}) quadqquad \ & = sumlimits_{mathbf{Z}}q(mathbf{Z}) lnfrac{P(mathbf{X}|oldsymbol{ heta})q(mathbf{Z}) P(mathbf{X},mathbf{Z}|oldsymbol{ heta})}{q(mathbf{Z}) P(mathbf{X},mathbf{Z}|oldsymbol{ heta})} \ & = sumlimits_{mathbf{Z}}q(mathbf{Z}) frac{P(mathbf{X},mathbf{Z}|oldsymbol{ heta})}{q(mathbf{Z})} + sumlimits_{mathbf{Z}}q(mathbf{Z}) ln frac{P(mathbf{X}|oldsymbol{ heta})q(mathbf{Z}) }{P(mathbf{X},mathbf{Z}|oldsymbol{ heta})} \ & = sumlimits_{mathbf{Z}}q(mathbf{Z}) frac{P(mathbf{X},mathbf{Z}|oldsymbol{ heta})}{q(mathbf{Z})} + sumlimits_{mathbf{Z}}q(mathbf{Z}) ln frac{q(mathbf{Z}) }{P(mathbf{Z}|mathbf{X},oldsymbol{ heta})} \ & = underbrace{sumlimits_{mathbf{Z}}q(mathbf{Z}) frac{P(mathbf{X},mathbf{Z}|oldsymbol{ heta})}{q(mathbf{Z})}}_{L(q,oldsymbol{ heta})} + mathbb{KL}(q(mathbf{Z})||P(mathbf{Z}|mathbf{X},oldsymbol{ heta}))) ag{2.1} end{align*}$

$(2.1)$ 式的右端為 $q(mathbf{Z})$ 和後驗分布 $P(mathbf{Z}|mathbf{X},oldsymbol{ heta})$ 的KL散度，由此 $lnP(mathbf{X}|oldsymbol{ heta})$ 被分解為 $L(q,oldsymbol{ heta})$ 和 $mathbb{KL}(q||p)$ 。由於KL散度總大於等於0，所以 $L(q,oldsymbol{ heta})$ 是 $lnP(mathbf{X}|oldsymbol{ heta})$ 的下界，如圖：

由此可將EM演算法視為一個坐標提升(coordinate ascent)的方法，分別在E步和M步不斷提升下界 $L(q,oldsymbol{ heta})$ ，進而提升 $lnP(mathbf{X}|oldsymbol{ heta})$ 。

在E步中，固定參數 $oldsymbol{ heta}^{old}$ ，當且僅當 $mathbb{KL}(q||p) = 0$ ，即 $L(q,oldsymbol{ heta}) = lnP(mathbf{X}|oldsymbol{ heta})$ 時， $L(q,oldsymbol{ heta})$ 達到最大，而 $mathbb{KL}(q||p) = 0$ 的條件是 $q(mathbf{Z}) = P(mathbf{Z}|mathbf{X}, oldsymbol{ heta})$ ，因此這就是式 $(1.2)$ 中選擇分布 $P(mathbf{Z}|mathbf{X},oldsymbol{ heta}^{old})$ 的原因，如此一來 $L(q,oldsymbol{ heta})$ 也就與 $(1.3)$ 式一致了。

在M步中，固定分布 $P(mathbf{Z}|mathbf{X},oldsymbol{ heta}^{old})$ ，選擇新的 $oldsymbol{ heta}^{new}$ 來極大化 $L(q,oldsymbol{ heta})$ 。同時由於 $P(mathbf{Z}|mathbf{X},oldsymbol{ heta}^{old}) eq P(mathbf{Z}|mathbf{X},oldsymbol{ heta}^{new})$ ，所以 $mathbb{KL}(P(mathbf{Z}|mathbf{X},oldsymbol{ heta}^{old}) || P(mathbf{Z}|mathbf{X},oldsymbol{ heta}^{new})) > 0$ ，導致 $lnP(mathbf{X}|oldsymbol{ heta})$ 提升的幅度會大於 $L(q,oldsymbol{ heta})$ 提升的幅度，如圖：

因此在EM演算法的迭代過程中，通過交替固定 $oldsymbol{ heta}$ 和 $P(mathbf{Z}|mathbf{X},oldsymbol{ heta}^{old})$ 來提升下界 $L(q,oldsymbol{ heta})$ ，進而提升對數似然函數 $L(oldsymbol{ heta})$ ，從而在隱變數存在的情況下實現了極大似然估計。在下一篇中將探討EM演算法的具體應用。

Reference :

Christopher M. Bishop. Pattern Recognition and Machine Learning

2. 李航.《統計學習方法》

3. CS838. The EM Algorithm