EM算法及其應用（一）

EM算法及其應用（一）

EM算法及其應用（二）： K-means 與 高斯混合模型

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

\(\large{\S} \normalsize\mathrm{1}\) 先驗知識

---

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

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

\[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\)，則：

\[f\left(\sum\limits_{i=1}^n \lambda_i x_i\right) \leq \sum\limits_{i=1}^n \lambda_i f(x_i) \]

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

\[ln\left(\sum\limits_{i=1}^n \lambda_i x_i\right) \geq \sum\limits_{i=1}^n \lambda_i ln(x_i) \]

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

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

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)}\right] = \mathbb{E}_{\mathrm{x}\sim p}\left[log\,p(x) - log\,q(x) \right ] \]

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) \neq \mathbb{KL}(q||p)\)，因此不能作為一個距離指標。

極大似然估計和極大后驗估計

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

\[L(\theta) = ln \,P(X|\theta) \]

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

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

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

\[P(\theta|X) = \frac{P(X,\theta)}{P(X)} = \frac{P(\theta)P(X|\theta)}{P(X)} \propto P(\theta)P(X|\theta) \]

上式中\(P(X)\)不依賴於\(\theta\)因而為常數項可以舍去，則最終結果為 \(\theta = \mathop{\arg\max}\limits_{\theta}P(\theta)P(X|\theta)\)

\(\large{\S} \normalsize\mathrm{2}\) EM算法初探

---

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

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

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

\[L(\boldsymbol{\theta}) = ln\,P(\mathbf{X}|\boldsymbol{\theta}) = ln\,\sum\limits_{\mathbf{Z}}P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta}) \tag{1.1} \]

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

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

1. \(P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)}) \geqslant 0\)
2. \(\sum_{\mathbf{Z}}P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)}) = 1\)
3. \(ln(\cdot)\)為凹函數

因而可以利用Jensen不等式求出\(L(\boldsymbol{\theta})\)的下界：

\[\begin{align} L(\boldsymbol{\theta}) = ln\,\sum\limits_{\mathbf{Z}}P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta}) &= ln\,\sum\limits_{\mathbf{Z}}P(\mathbf{Z|\mathbf{X},\boldsymbol{\theta}^{(t)}})\frac{P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta}) }{P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)})} \tag{1.2}\\ & \geqslant \sum\limits_{\mathbf{Z}}P(\mathbf{Z|\mathbf{X},\boldsymbol{\theta}^{(t)}}) \,ln\frac{P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta}) }{P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)})} \tag{1.3} \\ & = \underbrace{\sum\limits_{\mathbf{Z}}P(\mathbf{Z|\mathbf{X},\boldsymbol{\theta}^{(t)}}) \,ln{P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta}) }}_{\mathcal{Q}(\boldsymbol{\theta},\boldsymbol{\theta}^{(t)})} \;\;\underbrace{- \sum\limits_{\mathbf{Z}}P(\mathbf{Z|\mathbf{X},\boldsymbol{\theta}^{(t)}}) \,ln{P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)})}}_{entropy} \tag{1.4} \end{align} \]

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

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

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

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

\[\begin{align*} & \sum\limits_{\mathbf{Z}}P(\mathbf{Z|\mathbf{X},\boldsymbol{\theta}^{(t)}}) \,ln\frac{P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta}^{(t)}) }{P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)})} \\ =\;\; & \sum\limits_{\mathbf{Z}}P(\mathbf{Z|\mathbf{X},\boldsymbol{\theta}^{(t)}}) \,ln\frac{P(\mathbf{Z|\mathbf{X},\boldsymbol{\theta}^{(t)}})P(\mathbf{X}|\boldsymbol{\theta}^{(t)})}{P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}^{(t)})} \\ = \;\; & \sum\limits_{\mathbf{Z}}P(\mathbf{Z|\mathbf{X},\boldsymbol{\theta}^{(t)}}) \,lnP(\mathbf{X}|\boldsymbol{\theta}^{(t)}) \\ = \;\; & lnP(\mathbf{X}|\boldsymbol{\theta}^{(t)}) \;\;=\;\; L(\boldsymbol{\theta}^{(t)}) \end{align*} \]

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

EM算法流程：

輸入： 觀測數據\(\mathbf{X}\)，隱變量\(\mathbf{Z}\)，聯合概率分布\(P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})\)

輸出：模型參數\(\boldsymbol{\theta}\)

1. 初始化參數\(\boldsymbol{\theta}^{(0)}\)
2. E步： 利用當前參數\(\boldsymbol{\theta}^{(t)}\)計算Q函數， \(\mathcal{Q}(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)}) = \sum\limits_{\mathbf{Z}}P(\mathbf{Z|\mathbf{X},\boldsymbol{\theta}^{(t)}}) \,ln{P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})}\)
3. M步： 極大化Q函數，求出相應的 \(\boldsymbol{\theta} = \mathop{argmax}\limits_{\boldsymbol{\theta}}\mathcal{Q}(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)})\)
4. 重復 2. 和3. 步直至收斂。

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

\[\boldsymbol{\theta} = \mathop{argmax}\limits_{\boldsymbol{\theta}} \left[\mathcal{Q}(\boldsymbol{\theta}, \boldsymbol{\theta}^{(t)}) + ln(p(\boldsymbol{\theta})\right] \]

\(\large{\S} \normalsize\mathrm{3}\) EM算法深入

---

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

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

\[\begin{align*} L(\boldsymbol{\theta}) = lnP(\mathbf{X}|\boldsymbol{\theta}) &= \sum\limits_{\mathbf{Z}}q(\mathbf{Z})\,lnP(\mathbf{X}|\boldsymbol{\theta}) \quad\qquad \text{上下同乘以 $q(\mathbf{Z}) \,P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})$}\\ & = \sum\limits_{\mathbf{Z}}q(\mathbf{Z}) ln\frac{P(\mathbf{X}|\boldsymbol{\theta})q(\mathbf{Z}) P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})}{q(\mathbf{Z}) P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})} \\ & = \sum\limits_{\mathbf{Z}}q(\mathbf{Z}) ln\frac{P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})}{q(\mathbf{Z})} + \sum\limits_{\mathbf{Z}}q(\mathbf{Z}) ln \frac{P(\mathbf{X}|\boldsymbol{\theta})q(\mathbf{Z}) }{P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})} \\ & = \sum\limits_{\mathbf{Z}}q(\mathbf{Z}) ln\frac{P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})}{q(\mathbf{Z})} + \sum\limits_{\mathbf{Z}}q(\mathbf{Z}) ln \frac{q(\mathbf{Z}) }{P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta})} \\ & = \underbrace{\sum\limits_{\mathbf{Z}}q(\mathbf{Z}) ln\frac{P(\mathbf{X},\mathbf{Z}|\boldsymbol{\theta})}{q(\mathbf{Z})}}_{L(q,\boldsymbol{\theta})} + \mathbb{KL}(q(\mathbf{Z})||P(\mathbf{Z}|\mathbf{X},\boldsymbol{\theta}))) \tag{2.1} \end{align*} \]

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

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


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

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

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

Reference:

1. Christopher M. Bishop. Pattern Recognition and Machine Learning
2. 李航. 《統計學習方法》
3. CS838. The EM Algorithm

/