白話EM演算法

hearthougan發表於2018-09-04
演算法

EM演算法其實就是一種函式的近似求解方法，它是一種迭代演算法，在很多演算法中都有著重要的應用，比如HMM的學習問題，GMM，語音識別中的MM-GMM以及主題模型PLSA等等，所以掌握EM的求解方式還是很必要的。本文參考李航博士的《統計學習方法》，但前提是需要了解EM以及高斯混合分佈的基本概念，我不從最基礎說起，希望能說的明白。

EM演算法可以解決極大似然估計引數的求解方法，只不過當問題的分佈包含一個隱藏變數的時候，無法通過直接取對數，再求導得出極值，用這個極值作為模型的引數。所以要想用極大似然估計的思想，就必須想一種辦法來求解引數，這就是EM演算法。既然是極大似然估計的思想，那就先簡單說一下極大似然估計。

極大似然估計

極大似然估計是頻率學派的思想，頻率學派他們認為隨機事件的分佈都是一定的，只不過不知道而已，要想求出這個分佈（假設為p(\theta ,x)，就是解出分佈的引數），就需要從整體中抽出一些樣本x_{1},x_{2},\cdots x_{n}，由於頻率學派認為，分佈是固定的，所以可以直接連乘，L(\theta )=p(\theta ,x_{1})p(\theta ,x_{2})\cdots \cdots p(\theta ,x_{n})，然後取對數，再求導，得極值\theta ^{'}，\theta ^{'}就認為是最可能導致這個樣本出現的引數。這就是極大似然估計的思想。這裡沒有詳細敘述極大似然估計，想細究的可以參考極大似然估計。

極大似然估計就是估計最有可能導致這些結果（樣本）出現的原因（引數）。故：

上述我們直接取對數求導的方法，在有時候可行的，比如，我們的樣本服從的分佈是單一分佈，比如服從二項分佈，或者正態分佈等等，這種是可行的。那如果我門樣本服從的分佈不是單一分佈，比如高斯混合分佈，每一個樣本來自哪個高斯分佈你都不知道，你怎麼連乘？比如OSR中幀與音素的對應關係，或者大家常說的身高的例子等等，他們都是服從正態分佈的，以身高為例，因為常見，所以好理解。

首先你知道人的身高服從正態分佈，男女的正態分佈不一樣，也就是\mu _{1},\sigma _{1},\mu _{2},\sigma _{2}是未知的，如下圖。

然後給你一堆身高的資料，但是不知道這些資料對應的是男是女，那麼讓你估計一下GMM模型的引數。

GMM混合分佈

GMM模型引數估計，肯定不能再直接用某個高斯分佈直接相乘再求導，因為每個資料屬於哪個分佈你都不知道。因為這裡包含了一個隱藏的變數性別，所以要換一種辦法來求解了。GMM模型引數估計，是一個無監督的過程，非常類似於KMeans，比如要把資料分成k類，KMeans的估算過程：

1、給定k箇中心點，計算每個資料點到這k箇中心點的距離，找出距離最近的中心點m，那麼這個該資料就屬於類別m。

2、經過步驟1之後，現在就有k個類別，每個類別擁有一些資料。現在要重新計算中心點，就是把每個類別的資料求均值，作為新的中心點。

3、重複1、2，直至中心點不在移動，或者移動的距離小於某個閾值，則停止。

舉個KMeans工作的小例子，來源HMM-GMM。

KMeans過程

從KMeans的過程，如上圖，你可以看到，這是一步一步迭代的過程。GMM也是一樣，給定一些初值\mu_{1}^{0},\sigma_{1}^{0},\mu_{2}^{0},\sigma _{2}^{0}，根據這個初值，利用高斯分佈，求出每個資料屬於每個分佈的概率，當然，那個概率大，就屬於那個分佈。然後根據每個分佈的資料重新估算出\mu_{1}^{i},\sigma_{1}^{i},\mu_{2}^{i},\sigma _{2}^{i}，然後再重新估計，直至引數變化不大。所以兩者思路是一樣的，

當然估計GMM引數不是上述那麼簡單，但思想是一樣的。GMM引數估計的過程就是用EM來更新引數的。將GMM的目的，就是為了讓你明白，有了隱含變數，就不能直接求解，但是思想還是極大似然估計，只不過這個表示式，不用取對數求導等於0的方式來解，而換成了EM而已。接下來，就來講解EM演算法。將按照李航博士的《統計學習方法》符號表示方法，儘量我自己闡述。

儘量白話EM演算法

EM演算法本質上就是極大似然估計，它符合極大似然估計的思想，只不過它是一種迭代演算法，它是通過迭代來逼近極大似然估計值。下面就從極大似然估計來推匯出EM演算法。

假設你的觀測變數資料是Y，隱藏變數資料是Z，模型分佈為p(Y|\theta ),待估計的引數是\theta。

1、首先根據條件寫出似然估計的函式表示式：

p(Y|\theta )=\sum_{Z}p(Y,Z|\theta )

2、取對數

\begin{matrix} L(\theta )=logp(Y|\theta )=log\sum_{Z}p(Y,Z|\theta )\\\: \: \: =log\left ( \sum_{Z}p(Y,|Z,\theta )p(Z|\theta ) \right ) \end{matrix}

從這個式子，可以看出，求極大似然估計時，你無法利用求導等於等0的方式求出引數\theta，因為有隱藏變數Z存在，也求不出來引數\theta。我們說過EM演算法是一個迭代的方法，那麼它怎麼迭代去逼近極大似然估計呢？

一般地我們求argmax(f(\theta ))

1、假如在f(\theta )的某點\theta _{i}我們找到一個g(\theta )使得，f(\theta )\geq g(\theta )，且在某一點\theta _{i}處有f(\theta _{i})=g(\theta _{i})，那麼g(\theta )就是f(\theta )的一個最大下界，如下圖紅色曲線。

找到最大下界

2、然後我們找下界的極大值，以極大值對應的點\theta _{i+1}。

3、重複1、2，直至引數\theta不在變動或者小於某個閾值即可。

那麼，這一過程，就可以用影象描述為（自己畫的，比較醜，將就著看吧）：

一幅動圖說明EM的計算過程

上圖就是EM的整個過程，若你只想瞭解EM的工作的大概過程，到這裡即可。既然是演算法，就必須要有嚴謹的推理，那麼接下來，就是理論推導部分。

Jensen不等式

設f(x)是定義在實數上的函式，若對於所有的實數x，f(x)的二次導數大於等於0，則f(x)是凸函式（簡記：上凸下凹，上下是指的是函式的開口）。  

則Jensen不等式：若f(x)是凸函式，X是隨機變數，那麼：E(f(x))\geq f(E(x))。當且僅當X是常量時，等號成立。

（簡計：函式的期望大於期望的函式），Jensen不等式很簡單，比較容易理解，在這裡不過多闡述，只以下圖簡單地展示一下。

Jensen不等式的簡單示例

EM演算法通過不斷地極大化下界，從而來實現對極大似然的估計。假設通過第i次迭代之後得到的模型引數為\theta _{i}，既然我們是一個不斷極大化的過程，那麼我們肯定是希望，第i+1次迭代得到的值大於第i得到的值，即L(\theta _{i+1})> L(\theta _{i})，就是通過這樣不斷地迭代，來達到L(\theta )的極大值，那麼求第i+1次的引數\theta _{i+1}，就可以由差的形式來表述：

\begin{matrix} L(\theta )-L(\theta _{i})=log(\sum_{Z}p(Y|Z,\theta )p(Z|\theta ))-logp(Y|\theta _{i}) \\ \\ \, \, \, \, \, \, \: \: \: \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; =log(\sum_{Z}p(Z|Y,\theta _{i})\frac{p(Y|Z,\theta )p(Z|\theta )}{p(Z|Y,\theta _{i})})-logp(Y|\theta _{i}) \end{matrix}

這裡做一個恆等變換，這樣就是把上式的第一項log中的式子變為一個函式的期望，其中p(Y|Z,\theta _{i})表示概率，\frac{p(Y|Z,\theta )p(Z|\theta )}{p(Y|Z,\theta _{i})})就是對應的函式，由於log函式是一個凸函式，那麼上式用Jensen不等式求解下界：

\begin{matrix} L(\theta )-L(\theta _{i})=log(\sum_{Z}p(Z|Y,\theta _{i})\frac{p(Y|Z,\theta )p(Z|\theta )}{p(Z|Y,\theta _{i})})-logp(Y|\theta _{i})\\ \\ \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \geq \sum_{Z}p(Z|Y,\theta _{i})log\frac{p(Y|Z,\theta )p(Z|\theta )}{p(Z|Y,\theta _{i})})-logp(Y|\theta _{i})\\ \\= \sum_{Z}p(Z|Y,\theta _{i})log\frac{p(Y|Z,\theta )p(Z|\theta )}{p(Z|Y,\theta _{i})p(Y|\theta _{i})}) \end{matrix}

把上式的L(\theta _{i})移到不等式右側，那麼我們的下界就找到了B(\theta ,\theta _{i})，令

B(\theta ,\theta _{i})= L(\theta _{i})+\sum_{Z}p(Z|Y,\theta _{i})log\frac{p(Y|Z,\theta )p(Z|\theta )}{p(Z|Y,\theta _{i})p(Y|\theta _{i})})

則有 L(\theta)\geq B(\theta ,\theta _{i})且在\theta _{i}處的L(\theta_{i})\geq B(\theta_{i} ,\theta _{i})，故B(\theta ,\theta _{i})是L(\theta)的一個下界，那麼剩下的我們就是去極大化這個下界，因為使得B(\theta ,\theta _{i})增大的\theta同樣會使得L(\theta )增大，我們每次選擇的\theta _{i+1}就要使得L(\theta )儘可能的大，因此我們才說要極大化這個下界。我們在極大化B(\theta ,\theta _{i})的過程中，對於變數\theta來說是常量的項，都可以被省略。則：

\begin{matrix} \theta _{i+1}=argmax B(\theta ,\theta _{i})\\ \\ \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \: \: \: \: \: \: \: \:=argmax_{\theta }(L(\theta _{i})+\sum_{Z}p(Z|Y,\theta _{i})log\frac{p(Y|Z,\theta )p(Z|\theta )}{p(Z|Y,\theta _{i})p(Y|\theta _{i})}))\\ \\\; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \;=argmax_{\theta }\sum_{Z}[p(Z|Y,\theta _{i})log(p(Y|Z,\theta )p(Z|\theta ))]\\ \\ \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; =argmax_{\theta }\sum_{Z}p(Z|Y,\theta _{i})log(p(Y,Z|\theta)) \\ \\ \; \; \; \; \; \; \;=argmax_{\theta }Q(\theta ,\theta _{i})\end{matrix}

在求下界的過程中，我們得到Q函式，這是所有EM演算法都共有的步驟，求下界的目的就是為了得到這個Q函式，極大化下界就是極大化Q函式。Q函式是所有需要用EM演算法來求解問題的共有步驟，每個問題的Q函式都不同的，比如GMM和HMM學習問題的Q函式是不一樣的，但是形式是一樣的，所以求Q函式是關鍵，解Q函式不難，按照正常的方法極大化Q函式即可。

EM演算法中E步就是求下界，也就是等價於求Q函式；M步就是極大化這個Q函式

（Q函式）在給定觀測資料Y和當前引數\theta _{i}下，完全資料的對數似然函式logp(Y,Z|\theta )關於對未觀測資料Z的條件概率分佈p(Z|Y,\theta _{i})的期望稱之為Q函式

\begin{matrix} Q(\theta ,\theta _{i})=E_{Z}[log(p(Y,Z|\theta))|Y,\theta _{i}]\\ \\ \, \, \, \, \: \: \; \; \; \; \; \; \; \; \; \; \; \; \; \; \; =\sum_{Z}p(Z|Y,\theta _{i})log(p(Y,Z|\theta)) \end{matrix}

EM演算法對初值是敏感的，就是說有可能找到的不是L(\theta )的全域性極大值，有可能是區域性極小值，如下圖。

EM演算法初值敏感

 

至此EM演算法大概算是交代完了，對EM演算法的收斂性，這列沒有繼續證明，看懂了EM思想，這個就不是什麼難點。本文是對《統計學習方法》的個人理解。如有不正確之處，歡迎指正。

參考

李航《統計學習方法》

相關文章