EM演算法也稱期望最大化（Expectation-Maximum,簡稱EM）演算法，它是一個基礎演算法，是很多機器學習領域演算法的基礎，比如隱式馬爾科夫演算法（HMM）， LDA主題模型的變分推斷等等。本文就對EM演算法的原理做一個總結。
1. EM演算法要解決的問題
我們經常會從樣本觀察資料中，找出樣本的模型引數。 最常用的方法就是極大化模型分佈的對數似然函式。
但是在一些情況下，我們得到的觀察資料有未觀察到的隱含資料，此時我們未知的有隱含資料和模型引數，因而無法直接用極大化對數似然函式得到模型分佈的引數。怎麼辦呢？這就是EM演算法可以派上用場的地方了。
EM演算法解決這個的思路是使用啟發式的迭代方法，既然我們無法直接求出模型分佈引數，那麼我們可以先猜想隱含資料（EM演算法的E步），接著基於觀察資料和猜測的隱含資料一起來極大化對數似然，求解我們的模型引數（EM演算法的M步)。由於我們之前的隱藏資料是猜測的，所以此時得到的模型引數一般還不是我們想要的結果。不過沒關係，我們基於當前得到的模型引數，繼續猜測隱含資料（EM演算法的E步），然後繼續極大化對數似然，求解我們的模型引數（EM演算法的M步)。以此類推，不斷的迭代下去，直到模型分佈引數基本無變化，演算法收斂，找到合適的模型引數。
從上面的描述可以看出，EM演算法是迭代求解最大值的演算法，同時演算法在每一次迭代時分為兩步，E步和M步。一輪輪迭代更新隱含資料和模型分佈引數，直到收斂，即得到我們需要的模型引數。
一個最直觀瞭解EM演算法思路的是K-Means演算法，見之前寫的K-Means聚類演算法原理。在K-Means聚類時，每個聚類簇的質心是隱含資料。我們會假設$K$個初始化質心，即EM演算法的E步；然後計算得到每個樣本最近的質心，並把樣本聚類到最近的這個質心，即EM演算法的M步。重複這個E步和M步，直到質心不再變化為止，這樣就完成了K-Means聚類。
當然，K-Means演算法是比較簡單的，實際中的問題往往沒有這麼簡單。上面對EM演算法的描述還很粗糙，我們需要用數學的語言精準描述。
2. EM演算法的推導
對於$m$個樣本觀察資料$x=(x^{(1)},x^{(2)},...x^{(m)})$中，找出樣本的模型引數$\theta$, 極大化模型分佈的對數似然函式如下：$$\theta = arg \max \limits_{\theta}\sum\limits_{i=1}^m logP(x^{(i)};\theta)$$
如果我們得到的觀察資料有未觀察到的隱含資料$z=(z^{(1)},z^{(2)},...z^{(m)})$，此時我們的極大化模型分佈的對數似然函式如下：$$\theta = arg \max \limits_{\theta}\sum\limits_{i=1}^m logP(x^{(i)};\theta) = arg \max \limits_{\theta}\sum\limits_{i=1}^m log\sum\limits_{z^{(i)}}P(x^{(i)}， z^{(i)};\theta) $$
上面這個式子是沒有 辦法直接求出$\theta$的。因此需要一些特殊的技巧，我們首先對這個式子進行縮放如下：$$\begin{align} \sum\limits_{i=1}^m log\sum\limits_{z^{(i)}}P(x^{(i)}， z^{(i)};\theta) & = \sum\limits_{i=1}^m log\sum\limits_{z^{(i)}}Q_i(z^{(i)})\frac{P(x^{(i)}， z^{(i)};\theta)}{Q_i(z^{(i)})} \\ & \geq \sum\limits_{i=1}^m \sum\limits_{z^{(i)}}Q_i(z^{(i)})log\frac{P(x^{(i)}， z^{(i)};\theta)}{Q_i(z^{(i)})} \end{align}$$
上面第(1)式引入了一個未知的新的分佈$Q_i(z^{(i)})$，第(2)式用到了Jensen不等式：$$log\sum\limits_j\lambda_jy_j \geq \sum\limits_j\lambda_jlogy_j\;\;, \lambda_j \geq 0, \sum\limits_j\lambda_j =1$$
或者說由於對數函式是凹函式，所以有:$$f(E(x)) \geq E(f(x))\;\; 如果f(x) 是凹函式$$
此時如果要滿足Jensen不等式的等號，則有：$$\frac{P(x^{(i)}， z^{(i)};\theta)}{Q_i(z^{(i)})} =c, c為常數$$
由於$Q_i(z^{(i)})$是一個分佈，所以滿足：$$ \sum\limits_{z}Q_i(z^{(i)}) =1 $$
從上面兩式，我們可以得到：$$Q_i(z^{(i)}) = \frac{P(x^{(i)}， z^{(i)};\theta)}{\sum\limits_{z}P(x^{(i)}， z^{(i)};\theta)} = \frac{P(x^{(i)}， z^{(i)};\theta)}{P(x^{(i)};\theta)} = P( z^{(i)}|x^{(i)};\theta))$$
如果$Q_i(z^{(i)}) = P( z^{(i)}|x^{(i)};\theta)) $, 則第(2)式是我們的包含隱藏資料的對數似然的一個下界。如果我們能極大化這個下界，則也在嘗試極大化我們的對數似然。即我們需要最大化下式：$$ arg \max \limits_{\theta} \sum\limits_{i=1}^m \sum\limits_{z^{(i)}}Q_i(z^{(i)})log\frac{P(x^{(i)}， z^{(i)};\theta)}{Q_i(z^{(i)})}$$
去掉上式中為常數的部分，則我們需要極大化的對數似然下界為：$$ arg \max \limits_{\theta} \sum\limits_{i=1}^m \sum\limits_{z^{(i)}}Q_i(z^{(i)})log{P(x^{(i)}, z^{(i)};\theta)}$$
上式也就是我們的EM演算法的M步，那E步呢？注意到上式中$Q_i(z^{(i)})$是一個分佈，因此$\sum\limits_{z^{(i)}}Q_i(z^{(i)})log{P(x^{(i)}， z^{(i)};\theta)}$可以理解為$logP(x^{(i)}， z^{(i)};\theta)$基於條件概率分佈$Q_i(z^{(i)})$的期望。
至此，我們理解了EM演算法中E步和M步的具體數學含義。
3. EM演算法流程
現在我們總結下EM演算法的流程。
輸入：觀察資料$x=(x^{(1)},x^{(2)},...x^{(m)})$，聯合分佈$p(x,z ;\theta)$, 條件分佈$p(z|x; \theta)$, 最大迭代次數$J$。
1) 隨機初始化模型引數$\theta$的初值$\theta^{0}$。
2） for j from 1 to J開始EM演算法迭代：
a) E步：計算聯合分佈的條件概率期望：$$Q_i(z^{(i)}) = P( z^{(i)}|x^{(i)}，\theta^{j})) $$$$L(\theta, \theta^{j}) = \sum\limits_{i=1}^m\sum\limits_{z^{(i)}}Q_i(z^{(i)})log{P(x^{(i)}， z^{(i)};\theta)}$$
b) M步：極大化$L(\theta, \theta^{j})$,得到$\theta^{j+1}$: $$\theta^{j+1} = arg \max \limits_{\theta}L(\theta, \theta^{j}) $$
c) 如果$\theta^{j+1}$已收斂，則演算法結束。否則繼續回到步驟a)進行E步迭代。
輸出：模型引數$\theta$。
4. EM演算法的收斂性思考
EM演算法的流程並不複雜，但是還有兩個問題需要我們思考：
1） EM演算法能保證收斂嗎？
2） EM演算法如果收斂，那麼能保證收斂到全域性最大值嗎？
首先我們來看第一個問題, EM演算法的收斂性。要證明EM演算法收斂，則我們需要證明我們的對數似然函式的值在迭代的過程中一直在增大。即：$$\sum\limits_{i=1}^m logP(x^{(i)};\theta^{j+1}) \geq \sum\limits_{i=1}^m logP(x^{(i)};\theta^{j})$$
由於$$L(\theta, \theta^{j}) = \sum\limits_{i=1}^m\sum\limits_{z^{(i)}}P( z^{(i)}|x^{(i)};\theta^{j}))log{P(x^{(i)}， z^{(i)};\theta)}$$
令：$$H(\theta, \theta^{j}) = \sum\limits_{i=1}^m\sum\limits_{z^{(i)}}P( z^{(i)}|x^{(i)};\theta^{j}))log{P( z^{(i)}|x^{(i)};\theta)}$$
上兩式相減得到：$$ \sum\limits_{i=1}^m logP(x^{(i)};\theta) = L(\theta, \theta^{j}) - H(\theta, \theta^{j})$$
在上式中分別取$\theta$為$\theta^j$和$\theta^{j+1}$，並相減得到：$$\sum\limits_{i=1}^m logP(x^{(i)};\theta^{j+1}) - \sum\limits_{i=1}^m logP(x^{(i)};\theta^{j}) = [L(\theta^{j+1}, \theta^{j}) - L(\theta^{j}, \theta^{j}) ] -[H(\theta^{j+1}, \theta^{j}) - H(\theta^{j}, \theta^{j}) ] $$
要證明EM演算法的收斂性，我們只需要證明上式的右邊是非負的即可。
由於$\theta^{j+1}$使得$L(\theta, \theta^{j})$極大，因此有:$$L(\theta^{j+1}, \theta^{j}) - L(\theta^{j}, \theta^{j}) \geq 0$$
而對於第二部分，我們有：$$ \begin{align} H(\theta^{j+1}, \theta^{j}) - H(\theta^{j}, \theta^{j}) & = \sum\limits_{i=1}^m\sum\limits_{z^{(i)}}P( z^{(i)}|x^{(i)};\theta^{j})log\frac{P( z^{(i)}|x^{(i)};\theta^{j+1})}{P( z^{(i)}|x^{(i)};\theta^j)} \\ & \leq \sum\limits_{i=1}^mlog(\sum\limits_{z^{(i)}}P( z^{(i)}|x^{(i)};\theta^{j})\frac{P( z^{(i)}|x^{(i)};\theta^{j+1})}{P( z^{(i)}|x^{(i)};\theta^j)}) \\ & = \sum\limits_{i=1}^mlog(\sum\limits_{z^{(i)}}P( z^{(i)}|x^{(i)};\theta^{j+1})) = 0 \end{align}$$
其中第（4）式用到了Jensen不等式，只不過和第二節的使用相反而已，第（5）式用到了概率分佈累積為1的性質。
至此，我們得到了：$\sum\limits_{i=1}^m logP(x^{(i)};\theta^{j+1}) - \sum\limits_{i=1}^m logP(x^{(i)};\theta^{j}) \geq 0$, 證明了EM演算法的收斂性。
從上面的推導可以看出，EM演算法可以保證收斂到一個穩定點，但是卻不能保證收斂到全域性的極大值點，因此它是區域性最優的演算法，當然，如果我們的優化目標$L(\theta, \theta^{j}) $是凸的，則EM演算法可以保證收斂到全域性最大值，這點和梯度下降法這樣的迭代演算法相同。至此我們也回答了上面提到的第二個問題。
5. EM演算法的一些思考
如果我們從演算法思想的角度來思考EM演算法，我們可以發現我們的演算法裡已知的是觀察資料，未知的是隱含資料和模型引數，在E步，我們所做的事情是固定模型引數的值，優化隱含資料的分佈，而在M步，我們所做的事情是固定隱含資料分佈，優化模型引數的值。比較下其他的機器學習演算法，其實很多演算法都有類似的思想。比如SMO演算法（支援向量機原理(四)SMO演算法原理），座標軸下降法(Lasso迴歸演算法： 座標軸下降法與最小角迴歸法小結), 都使用了類似的思想來求解問題。
大家也可以去比較下這些演算法的優化方法，看思路上是不是有共同之處。
（歡迎轉載，轉載請註明出處。歡迎溝通交流： liujianping-ok@163.com）