EM演算法總結

05-06

最近因工作需要，要學習概率圖模型，來對一些序列對象進行一些建模。概率圖模型根據不同情況的不同概率假設，能夠有效解釋一些現實問題，使得我們的模型能夠在一定程度上擬合現實的問題，但是存在模型不確定、參數不確定，建模比較難的特點。因此有很多思想方法幫助我們去進行參數估計和建模，本文準備將EM演算法總結一下。

什麼是EM演算法

EM演算法全稱Expectation-Maximization演算法，即最大化期望演算法。根據這個名字就能很清楚知道這個演算法的核心元素，一個是期望，另一個是期望的最大化。而這兩個元素代表了EM演算法的兩個流程。首先EM演算法是一個迭代性質的演算法，即每個迭代都能得到新的經過學習後的參數，迭代的結束標誌一般是人為定義閾值或者參數收斂。其次EM演算法每個迭代主要執行兩個步驟：

求期望，得到一個關於參數的期望函數式。類似於我將樣本數據輸入到
對期望函數式求使其最大化的參數。

上述描述比較簡略，下面會詳細將這兩個步驟說明。

為什麼用EM演算法

EM演算法的好處，或者說它被需要的原因是什麼？其實之前我已經在HMM模型的參數估計章節中，對EM演算法的具體應用已經有了一個專門的描述，感興趣的同學可以去看我之前的文章：《隱馬爾科夫模型學習總結之三》。從HMM的例子中，我們知道當一個實際問題中，除了可觀測的變數外，還需要一些隱變數幫助我們去建模，此時EM演算法可以幫助對這種包含隱變數的問題進行參數估計，除了HMM外，還有k-means聚類演算法也運用了EM演算法，它的隱變數可以說是具體的那些聚類本身。

除此以外，我們會存在無法確定問題對應的具體的模型的樣式的時候。說的有點拗口，我舉個例子，大家感受一下。

logistics regression的二分類例子

首先，我舉個能確定具體模型的例子，即使用logitstic regression對二分類數據進行建模，其中樣本數據之間是相互獨立的，沒有關聯，且當前數據標籤只有兩種0和1。其實這種數據分布可以看成是一個二項分布，即重複n次獨立的伯努利試驗。在每次試驗中只有兩種可能的結果，而且兩種結果發生與否互相對立，並且相互獨立，與其它各次試驗結果無關。對應於我們的數據就是每個樣本相當於一次拋硬幣，要麼正面要麼反面，且兩次拋硬幣（兩個數據樣本）之間沒有任何關聯，互不影響。那麼我們假設模型參數為 $heta$ ,模型為 $h(x, heta)$ 樣本數據集為 ${(x_1,y_1),(x_2,y_2),...,(x_m,y_m)}$ ，m表示m個樣本，那麼我們給定參數的情況下，每個樣本數據標籤的條件概率為： $P(y_i| heta)=h(x_i, heta)^{y_i}(1-h(x_i, heta))^{(1-y_i)}$ ,這個式子其實是兩個式子的合併，我們知道logistics regression是生成樣本屬於標籤1的概率，因此有如下兩個式子： $P(y_i=1| heta)=h(x_i, heta)$ ， $P(y_i=0| heta)=1-h(x_i, heta)$ ，兩個式子合併，就成上面的式子。

由於樣本數據之間相互獨立，因此整個樣本集的預測目標概率為：

$P(y| heta)=P(y_1| heta)P(y_2| heta)...P(y_m| heta)$

我們訓練模型的利用已知的樣本數據，反推最大概率能夠得到樣本數據的模型參數。因此對於這種已經確定模型形式，需要進行參數估計的問題，我們很自然想到使用極大似然估計方法MLE，即 $argmax_{ heta}P(y|x, heta)$ 。由於連乘式的求極值很好不做，且概率的多次連乘，尤其是在數據樣本很大的情況下，很容易造成最後結果幾乎為0，因此考慮在概率式子之前增加log函數，由於log函數的單調性，保證了原來概率式的性質，且可以將連乘式化簡為求和式。

然後就是以 $heta$ 為參數，對 $L( heta)$ 求偏導數，並令偏導為0，求出 $heta$ 的表達式。

上述的例子說明只要是數據分布確定後，求具體數據分布的參數是比較容易的，一般是使用MLE方法。然而現實生活中存在很多數據分布無法確定的情形，如之前講的HMM中的問題，因為它引入了一個隱變數，帶來了不確定性，因此不能簡單地使用MLE方法，一步到位得將我們的參數估計出來，因此需要EM演算法，迭代得去估計我們的參數。

高斯混合模型應用EM演算法舉例

EM演算法問題描述

下面準備拿很多人都舉過的高斯混合模型的例子來講解EM演算法，參考了cs229以及徐亦達老師的講義。

假設當前數據符合上述分布，很明顯，該數據不是一種概率分布就能簡單得描述的，大致可以看出該數據有三個不同的概率分布，每個概率分布可先用高斯分布來解釋（分布中心數據密度大，周圍數據密度小，與高斯分布特徵相符）。高斯分布的概率密度公式如下：

其中， $mu$ 表示期望， $Sigma$ 表示方差。在混合高斯模型中，假設有k個高斯分布混合，概率密度公式為：

其中，有約束條件：

因此參數集合為：

若使用MLE方法來進行參數估計，則極大似然估計函數為：

如果單純使用求偏導的方式來求極值的參數是很難的，無法簡單做到。

此時，由於我們不知道這些數據具體從屬於哪個概率分布，因此引入一個隱變數Z，表示數據屬於哪個高斯分布。原始的MLE函數式為：

這個函數式無法一次使用求偏導的方法求出極值，因此需要迭代的方法，每個迭代產生參數估計。然後我們引入隱變數Z，在不改變原式子的情況下，我們可以採用貝葉斯公式以及期望積分的方式，將隱變數Z引入上式（之前講HMM的時候已經講過這種方法，感興趣的同學可以回頭去看一下），因此得到EM演算法每步迭代的問題描述式：

首先這是某步迭代的函數式，假設之前的迭代已經估計出了參數 $heta^{(g)}$ ,因此在此式中 $heta^{(g)}$ 是一個已知的常數值，這個在後面化簡推導的時候很有用。其次在外層增加了對隱變數z的積分，結合內部的貝葉斯條件概率公式，可以將z的影響消除。

備註：隱變數引入不是隨便的，有一些要求，比如需要引入後簡化模型的求解，另外不能改變原邊緣分布，即 $P(x)=int_{z}^{}P(x|z)P(z)dz$ ,在高斯混合模型中，引入的隱變數z表示對應數據屬於哪個高斯，因此P(z)屬於先驗知識，來自於上面提到的 $alpha_l$ ,表示高斯分布的權重。因此可以很容易證明 $P(x)=int_{z}^{}P(x|z)P(z)dz$ ，這邊就不寫了，有興趣的同學自己推導，很簡單。

EM演算法迭代收斂性證明

既然是迭代的方法，那麼勢必要保證這種方法最後能夠收斂，事實上EM演算法的收斂性是可以證明的，下面簡單對這個收斂性做一個證明。即證明：

首先根據條件概率公式，有：

等式兩邊同時加log，得到：

等式左邊引入隱變數z，並對z做積分的期望，右邊同理，可得：

（上述使用了這個等式： $int_{z}^{}P(z|x , heta^{(g)})=1$ ）

我們將上面等號右邊的被減子式定義為 $Q( heta, heta^{(g)})$ ,減子式定義為 $H( heta, heta^{(g)})$ ，因此 $log(P(X| heta))=Q( heta, heta^{(g)})-H( heta, heta^{(g)})$ ,我們要證明 $log(P(X| heta^{(g+1)}))geq log(P(X| heta^{(g)}))$ ，可以看到我們的Q項其實與我們一開始提出的EM演算法問題描述式是一樣的，即每個迭代， $heta^{(g+1)}$ 都是在Q項上求極大值後的參數，因此 $Q( heta^{(g+1)}, heta^{(g)})geq Q( heta^{(g)}, heta^{(g)})$ 是必然成立的（因為本來就是求極大值，因此極大值只可能比上個迭代大或者一樣）。那麼我們接下來只要證明H項在每個迭代都在逐漸減小或者不變，即： $H( heta^{(g+1)}, heta^{(g)})leq H( heta^{(g)}, heta^{(g)})$

證明：我們將要證明的不等式做一個轉換，問題可以轉換為：

對任意的 $heta$ ，我們要證明 $H( heta, heta^{(g)})leq H( heta^{(g)}, heta^{(g)})$ ，有如下化簡：

其中-log()是一個凸函數，函數圖如下：

對於凸函數我先介紹一個性質：Jensen不等式：

我們引入一個變數 $tin(0,1)$ ，假設當前數據區間為 $[x_1,x_2]$ ，那麼該區間內任意一個數都可以用如下式子表示： $(1-t)x_1+tx_2$ 。而將函數曲線的兩端的點連接構成一條直線，該直線上的點可以表示為：(1-t)f(x_1)+tf(x_2)。那麼在函數為凸的情況下，根據圖上的示例（雖然比較難看），可以看出 $(1-t)f(x_1)+tf(x_2)geq f((1-t)x_1+tx_2)$ ，這個就是Jensen不等式，應用到我們這邊證明，其實可以泛化為不等式左邊的式子可以看成是函數的期望，不等式右邊看成是期望的函數，其中t，1-t都是某個函數值的概率，因此可以寫成 $期望的函數leq 函數的期望$ ，即

而

是函數的期望，因此它比大於期望的函數：

因此，證畢。

使用EM演算法解決高斯混合模型問題

下面定義問題描述式中的具體項：

其中 $N(X|mu_l,Sigma_l)$ 表示一個高斯分布， $alpha_l$ 表示這個高斯分布的權重。k表示有k個不同的高斯分布，n表示數據樣本數。

定義 $P(X,z| heta)$ ，因為數據間獨立，且有貝葉斯公式：

因此有：

定義 $P(z|X, heta)$ ，根據數據獨立性，貝葉斯公式和積分期望，有如下式子：

因此可定義 $P(z|X, heta)$ 為：

將上述兩個定義式代入之前的Q項中（H項已不需要），開始EM演算法中的E-step：

E-step：

將式子

定義為 $f_i(z_i)$ ,將式子

定義為 $P(z_1,....z_n)$ 。

對上述式子肯定是要進行化簡的：

把上述式子展開，其實是一個包含N個sum項的式子，我們把第一個sum式拿出來單獨分析：

上式中，後面對z2...zn的積分式，其實可以看做是z1變數的邊緣概率分布，因為其他z都被積分積掉了（原諒我用這麼口語的方式表達），因此上式可化簡為：

很明顯，其他sum項也可以做這樣的化簡，因此對所有N個sum項都做如此化簡後，可得如下式子：

最後得到的式子中，加號前面的只包含參數 $alpha$ ，加號後面的只包含參數 $mu,Sigma$ ，因此可以單獨對兩個式子分別求極大值，然後求得本次迭代的參數。下面就是我們的M-step。

M-Step：

主要是開始對E-Step構造的式子進行參數最大化，按照上面說的，將式子根據加號一拆為二，首先估計參數 $alpha$ 。對式子求偏導，並令偏導為0，得到：

說明一下，為了描述方便，這裡定義 $z_i=l$ 。

這個式子求解方法已經很熟悉了，之前SVM以及HMM中都講過這個，就不多說了，直接使用拉格朗日乘子法，得到：

接下來就是求參數 $mu,Sigma$ 。這個推導比較複雜，因為牽扯到一些矩陣向量的求導的trick和小公式，先把公式列出來，有同學感興趣，我再把這些公式的推導過程列出來：

我們對Q項加號的第二項進行分析，即：

首先，我們對參數 $mu$ 進行估計，我們將下面的高斯分布的概率密度公式代入。

其中， $Sigma$ 的多項式在此式可以視作一個const最後，對式子求偏導，令偏導為0，得到：

接下來是對 $Sigma$ 的推導，為了推導簡便，我們改變求偏導的變數對象，轉而對 $Sigma^{-1}$ 求偏導，這樣利用上面的一些facts（fact2），然後得到：

這樣我們就將所有本迭代的參數求出來了，細心的同學可以發現所有參數裡面都包含一個式子： $P(l|X_i, heta)$ ，這個式子我們在之前的推導中已經得到，為：

因此，使用這個式子，加上迭代的方式，就可以一步一步逼近最優解。

備註：

放到最後一個講，但是非常重要的一個地方就是EM演算法不能保證對每個應用場景的問題都能都到全局最優解，如果我們的目標函數是凸函數（本文我們討論到的-log()），那麼可以保證得到全局最優解。但是，若我們的目標函數不是凸的，比如在聚類問題中，若對文本分類中的餘弦相似度計算函數就不能保證是凸的，此時EM演算法就有可能給出局部最優。

reference:

cs229 notes

徐亦達 notes

《數學之美》吳軍