混合模型與EM演算法

02-24

一、前言

在用機器學習方法對數據進行建模時，我們一方面總是希望所用的模型盡量簡單，這樣對模型相關的各類問題能夠直接求出解析解。但簡單的模型表達能力有限，無法對複雜的數據特性進行建模。通過將簡單模型升級改造為複雜模型可以大大提高模型的表達能力，但往往要面臨維度災難問題，很有可能當你將模型升級為可以表達期望的數據關係時，因為維度災難問題模型已經變得不可求解了。因此，我們希望找到一種方法，既能保持簡單模型的簡單求解性，又能使模型具有豐富的表達能力。混合模型即是這樣一種方法，它以簡單模型作為基礎模塊，通過合理的引入隱藏變數，將基礎模塊組合在一起形成具有豐富表達能力的「組合式」模型——混合模型（Mixture Models）。這是進行數據建模非常常用的一種方法，《隱馬爾科夫模型》一文即詳細介紹了此類模型的一個實例，它展示了如何利用混合模型的思想在控制模型複雜度的同時表達變數之間長距離的相關性。除了隱馬爾科夫模型之外，我們平常聽到的很多模型其本質都是混合模型，例如K均值演算法，混合高斯模型，以及線性動態系統（卡爾曼濾波、粒子濾波）等等。因為混合模型是一種引入隱變數的建模方法，因此也叫隱變數模型，根據引入的隱變數是離散的還是連續的，隱變數模型又分離散隱變數模型和連續隱變數模型，前面說的線性動態系統就是連續隱變數的例子。這篇文章主要探討離散隱變數模型的思想和本質，以及離散隱變數模型的最大似然問題的求解方法——EM演算法的演算法原理。

本文先從傳統視角介紹K均值聚類演算法和混合高斯模型的最大似然問題。隨後通過引入隱藏變數，解釋K均值聚類演算法和混合高斯模型為何是離散隱變數模型的實例，並藉助求解隱變數模型的最大似然問題的有力工具——EM演算法重新推導K均值聚類演算法和混合高斯模型的最大似然解，從而揭示EM演算法和傳統方法的關聯。最後，從更一般化的角度證明EM演算法為何能夠收斂到最大似然問題的極大值。增加K均值聚類相關的章節是為了展示離散隱變數模型的廣泛應用，以及用更多的實例解釋EM演算法，從而加深讀者對EM演算法的理解。如果您對K均值聚類不感興趣，完全可以跳過相關章節（第二、七節）。雖然您也可以跳過混合高斯模型的相關章節（第三、四、六節），直接閱讀有關EM演算法的理論介紹（第五、八節），但我還是建議您不要這樣做，因為混合高斯模型是引入EM演算法非常合適的例子，有具體的例子做輔助，會更容易理解EM演算法理論。如果您已經對混合高斯模型相關章節的內容非常熟悉了，此時便可直接參考EM演算法理論介紹的章節。

二、K均值聚類

作為混合模型最簡單的一個實例，K均值聚類在很多領域都有廣泛的應用。本節先從傳統的機器學習視角介紹K均值演算法，然後在深入介紹混合模型和EM演算法之後，再解釋K均值演算法如何是混合模型的一個實例，以及如何用EM演算法的理論求解K均值聚類問題。

K均值聚類是要解決這樣一類問題：給定 $N$ 個數據點 $left { m{x}_1, m{x}_2, cdots, m{x}_N ight }$ ，這些點的維數不重要，可以是一維的也可以是多維的，我們的目標是要將這些點分成 $K$ 組，使得同一組內的點之間更「相近」或「相似」，不同組之間的點「距離」更遠，一般假設 $K$ 是事先給定的值。為了用更形式化的語言描述這個問題，我們需要引入一些記號。為了表達 $K$ 個分組，我們引入 $K$ 個與數據點維數相同的向量 $m{mu}_k$ ，叫做分組的「原型（prototype）」，我們可以暫且將其理解為分組的聚類中心。為了表達將每個數據點分給哪個分組，對每一個數據點 $m{x}_n$ ，我們引入一個對應的二值指示變數 $m{gamma}_n$ ，該變數為 $K$ 維的，每個元素的取值只能是0或1，且同時只能有一個元素取1，其它都取0。如果 $gamma_{n,k} = 1$ ，即表示將 $m{x}_n$ 分配給第 $k$ 個分組。

如果用歐式距離表示數據點之間的「距離」或「不相似性」的話，那麼K均值聚類問題的目標就是，找到 $K$ 個分組的聚類中心 $m{mu}_k$ ，並將每個數據點分配到其中一個分組，使得所有數據點到其所屬的聚類中心的距離之和最小，目標函數為

$J = sum_{n=1}^{N} { sum_{k=1}^{K} {gamma_{n,k} ||m{x}_n - m{mu}_k || ^ 2}}$ （1）

該問題可用迭代方法求解：首先，假定 $m{mu}_k$ 已知且固定，讓 $J$ 對 $gamma_{n,k}$ 求最小。因為 $J$ 是 $gamma_{n,k}$ 的線性函數，且不同 $n$ 之間相互獨立，因此固定 $n$ 為某一個值，則有由 $||m{x}_n - m{mu}_k||^2$ 給定的 $K$ 個距離值，且每個距離值有一個二值變數 $gamma_{n,k}$ 作為係數。因為這 $K$ 個二值變數係數同時只能有一個取1，因此，為了使得最終的加權距離和最小，只能讓距離最小者的係數取1。也就是說，對每一個數據點 $m{x}_n$ ，我們將其分配給與之距離最小的那個分組：

$gamma_{n,k} = { egin{align} & 1, quad if k = arg min_{j} || m{x}_n - m{mu}_j||^2 \ &0, quad otherwise. end{align}$ （2）

然後，假定 $gamma_{n,k}$ 已知且固定，讓 $J$ 對 $m{mu}_k$ 求最小。因為 $J$ 是 $m{mu}_k$ 的二次函數，因此對其求導並令導數為0即可：

$frac{partial{J}}{partial{m{mu}_k}} = sum_{n}{2 gamma_{n,k} (m{x}_n - m{mu}_k)} = 0$ （3）

容易求解得到

$m{mu}_k = frac{sum_{n}{gamma_{n,k}m{x}_n}}{sum_{n}{gamma_{n,k}}}$ （4）

該表達式分母表示分配給第 $k$ 個分組的數據點個數，因此整個表達式表示將聚類中心設置為分配給該分組的所有數據點的平均值。正因為這個原因，該演算法被叫做K均值聚類演算法。

以上兩個步驟交互迭代進行，因為每一步都會減小目標值，因此演算法最終收斂到極小值。因為以上介紹中我們用歐式距離表達兩點間距離，因此該極小值一定為全局最小值。但如果使用更複雜的距離表達式，該演算法有可能收斂到局部極小值。有研究者專門研究了K均值演算法的收斂問題（MacQueen 1967），有興趣的讀者可以參考相關研究成果。

如果我們用更一般的形式 $V(m{x}, m{x^{}})$ 表示兩個數據點之間的距離，則目標函數將變為

$ilde{J} = sum_{n=1}^{N} {sum_{k=1}^{K}{gamma_{n,k} V(m{x}_n, m{mu}_k)}}$ （5）

此即K中心點演算法（K-medoids algorithm）。第一步與K均值演算法一樣，也是將每個數據點分配到與之距離最近的分組。而第二步稍微複雜一點，尤其當距離函數過於複雜時往往無法直接求解。此時一般限制聚類中心必須為該分組中的某一個數據點，因此求聚類中心的問題就簡化為選擇該分組中的某一個點，使得當該點為聚類中心時類內距離和最小。該簡化版本只要求已知 $m{x}$ 和 $m{x^{}}$ 時，式 $V(m{x}, m{x^{}})$ 可以求值即可。

三、混合高斯模型：表示

本節介紹混合模型另一個應用非常廣泛的實例——混合高斯模型，它比K均值演算法稍微複雜一點，但能帶領我們更接近EM演算法的神秘面貌。混合高斯模型一般用來對數據進行概率建模，跟上一節一樣，我們還是先從傳統的機器學習視角，建立最大似然問題，通過純數學的公式推導，得到該最大似然問題的迭代求解方法。隨後逐步引入EM演算法，在深入介紹EM演算法之後，將解釋該迭代演算法如何被裝進EM演算法的框架，即如何從EM演算法的角度推導出相同的結論。

高斯分布因其簡單、常見（大部分自然規律都與高斯分布有關）的特性，在機器學習領域有著非常廣泛的應用。然而單獨的高斯分布由於太簡單，表達能力有限，面對複雜的實際問題往往無能為力。如果我們取多個高斯分布，並以合適的係數加權平均，就能得到表達能力更豐富一類模型——高斯混合模型：

$p(m{x}) = sum_{k = 1}^{K} {pi_k N(m{x} | m{mu}_k, Sigma_k)}$ （6）

其中 $K$ 為高斯分布的個數， $m{mu}_k$ 、 $Sigma_k$ 為第 $k$ 個高斯分布的均值向量和方差， $pi_k$ 為第 $k$ 個高斯分布的加權係數，也叫混合係數。同樣的，我們不關心 $m{x}$ 的維數，可以是一維的，也可以是多維的，姑且用 $D$ 表示 $m{x}$ 的維數。因為概率密度函數要滿足非負和歸一化的條件，因此混合係數滿足

$egin{align} 0 le pi_k le 1 \ sum_{k = 1}^{K}{pi_k} = 1end{align}$ （7）

接下來，我們通過引入隱藏變數，構造混合高斯模型的另一種表示方法。我們可以這樣想像一個符合混合高斯模型分布的數據點的產生過程：首先以一定的概率從 $K$ 個高斯分布中選擇一個，其中第 $k$ 個高斯分布被選中的概率為 $pi_k$ ，然後依據被選中的高斯分布生成一個數據點。重複以上過程生成多個數據點，這些數據點將符合（6）式的混合高斯分布。為了表示選中哪一個高斯分布，我們引入一個二值指示變數 $m{z}$ ，該變數為 $K$ 維的，每個元素的取值只能是0或1，且同時只能有一個元素取1，其它都取0。如果 $z_{k} = 1$ ，即表示選擇第 $k$ 個高斯分布。由於第 $k$ 個高斯分布被選中的概率為 $pi_k$ ，即

$p(z_k=1) = pi_k$ （8）

因此，隨機變數 $m{z}$ 的分布可表示為

$p(m{z}) = prod_{k = 1} ^ {K} {pi_k^{z_k}}$ （9）

所以，當已知選中第 $k$ 個高斯分布時， $x$ 的條件分布為

$p(m{x} | z_k = 1) = N(m{x} | m{mu}_k, Sigma_k)$ （10）

或表示為

$p(m{x} | m{z}) = prod_{k = 1} ^ {K} {N(m{x} | m{mu}_k, Sigma_k) ^ {z_k}}$ （11）

從而， $m{x}$ 的邊緣分布為

$p(m{x}) = sum_{m{z}} {p(m{x} | m{z}) p(m{z})} = sum_{m{z}} {prod_{k = 1} ^ {K} { left[ pi_k N(m{x} | m{mu}_k, Sigma_k) ight] ^ {z_k}}}$ （12）

上式等號右邊又是求和、又是連乘積，看似很複雜，其實非常簡單。因為 $m{z}$ 總共只有 $K$ 種取值，對每一個取值，求和符號裡面的連乘積的結果只剩下一項，因為 $m{z}$ 的每一個取值（ $K$ 維二值向量）只有一個元素值為1。例如，當 $m{z}$ 取第 $k$ 種取值（ $z_k = 1$ ）時，連乘積只剩下 $pi_k N(m{x} | m{mu}_k, Sigma_k)$ 一項，因此（12）式簡化為

$p(m{x}) = sum_{m{z}} {p(m{x} | m{z}) p(m{z})} = sum_{k = 1} ^ {K} {pi_k N(m{x} | m{mu}_k, Sigma_k)}$ （13）

可見，用隱藏變數重新表達混合高斯模型後，最終得到的 $m{x}$ 的邊緣分布與（6）式一模一樣，這是必要的。因為無論你用什麼方式重新表達混合高斯模型，最後 $m{x}$ 的邊緣分布一定要「回歸」到（6）式，否則你的表達是毫無意義的，或者說是錯誤的。

如果將（6）式看做混合高斯模型的「判別式」表示方法，則（8）~（13）式可以看做混合高斯模型的「產生式」表示方法，因為這是從數據產生的角度對模型進行的表示。

從目前來看，我們用隱藏變數的方法重新表示混合高斯模型，似乎沒什麼受益。事實上，從後面的推導將會看出，這種表示方法能夠大大簡化最大似然問題的求解過程。因為如果直接用（6）式建模，似然函數將會是高斯函數和的連乘積形式，求解此種形式的最優化問題將會是非常困難的。但是引入隱藏變數後，我們就能夠使用聯合分布 $p(m{x}, m{z})$ 簡化問題，因為聯合分布直接是高斯函數的連乘積形式，這種形式便於求解析解，我們從接下來的推導中將會看到。在此之前，我們先看看變數 $m{z}$ 的條件概率的形式。如果用 $gamma(z_k)$ 表示條件概率 $p(z_k = 1 | m{x})$ ，利用貝葉斯定理，

$egin{align} gamma(z_k) = p(z_k = 1 | m{x}) &= frac{p(m{x} | z_k = 1) p(z_k = 1)}{sum_{j = 1} ^{K} {p(m{x} | z_j = 1) p(z_j = 1)}} \ &= frac{pi_k N(m{x} | m{mu}_k, Sigma_k)}{sum_{j = 1} ^{K} {pi_j N(m{x} | m{mu}_j, Sigma_k)}} end{align}$ （14）

上式是在已知觀測變數 $m{x}$ 的條件下 $z_k = 1$ 的後驗概率，而其先驗概率已經知道為 $pi_k$ 。

接下來我們就可以形式化地表示混合高斯模型的最大似然問題了。已知 $N$ 個數據點 $left { m{x}_1, m{x}_2, cdots, m{x}_N ight }$ ，我們希望用混合高斯模型對其進行概率密度建模。對每一個數據點 $m{x}_n$ ，都有一個與之對應的二值指示變數 $m{z}_n$ 。為了記號的緊湊，我們用一個 $N imes D$ 的矩陣 $X$ 表示所有的數據點，其中矩陣的第 $n$ 行為 $m{x}_n^T$ 。類似地，我們用一個 $N imes K$ 的矩陣 $Z$ 表示所有的二值指示變數，其中矩陣的第 $n$ 行為 $m{z}_n^T$ 。如果我們假設給定的數據集獨立同分布，那麼根據（6）式，對數似然函數為

$p(X | m{pi}, m{mu}, Sigma) = sum_{n = 1} ^ {N} {ln { { sum_{k = 1} ^ {K} {pi_k N (m{x} | m{mu}_k, Sigma_k) }}}}$ （15）

在求解該最大似然問題之前，我們需要首先討論一下應用最大似然方法求解混合高斯模型參數時將要面臨的兩個問題。第一個是奇點問題。簡單起見，假設某一個高斯分布的協方差矩陣為對角陣，即 $Sigma_k = sigma^2 I$ ，其中 $I$ 為單位陣。需要說明的是，不僅僅是對角陣，奇點問題對任意形式的協方差矩陣都是存在的。假設這個高斯分布的均值向量正好等於某一個數據點，即 $m{mu}_k = m{x}_n$ ，則該數據點在這個高斯分布下的似然值為

$N(m{x}_n | m{mu}_k, sigma^2 I) = frac{1}{(2 pi) ^{D / 2}} frac{1}{ { sigma}}$ （16）

當 $sigma o 0$ 時，上式將趨近於無窮大，因此似然函數也將趨近於無窮大。因此混合高斯模型的最大似然問題不是一個良好定義（well-posed）的問題，因為一旦其中某一個高斯分布「坍縮」到某個數據點上，似然值將變為無窮大，此時求得的參數將是不可信的。這可以看做是所有最大似然法求解面臨的過擬合問題。這個問題並不是無法解決的，以後我們將會看到，通過採用貝葉斯方法可以解決最大似然求解的過擬合問題，不過暫時我們不打算這麼做。因此，在應用最大似然法求解混合高斯模型的參數時，需要注意避免奇點問題帶來的病態解，可以在迭代求解的過程中，不斷檢測是否有高斯分布坍縮到某個數據點，倘若如此，則將其均值重新初始化為隨機值，並將其協方差重置為某個較大的值，然後繼續迭代。

第二個問題是唯一性問題。這個問題很好理解，假定我們已經求得了混合高斯模型的一組解，然後不改變這組解的值，僅僅將 $K$ 個參數重新排列組合，則可以得到另一組等價的解。這樣的組合共有 $K!$ 種。以 $K = 3$ 的一維模型為例，假設我們已經求得一組解： $m{pi} = [0.5, 0.3, 0.2] ^ T$ ， $mu_1 = 1$ ， $mu_2 = 2$ ， $mu_3 = 3$ ， $Sigma_1 = 0.1 ^ 2$ ， $Sigma_2 = 0.2 ^ 2$ ， $Sigma_3 = 0.3 ^ 2$ 。如果交換前兩組參數的順序，可以得到另一組等價的解： $m{pi} = [0.3, 0.5, 0.2] ^ T$ ， $mu_1 = 2$ ， $mu_2 = 1$ ， $mu_3 = 3$ ， $Sigma_1 = 0.2 ^ 2$ ， $Sigma_2 = 0.1 ^ 2$ ， $Sigma_3 = 0.3 ^ 2$ 。這些不同組合的等價解共有 $3! = 6$ 種。唯一性問題在需要對求解的參數進行解釋時需要認真對待，但如果僅僅是為了找到數據集的概率密度模型，那麼最終得到哪一組解並不重要，因為任意一組解都是等價的。

四、混合高斯模型：求解

求解的目的是要對（15）式求最大值，並求得取最大值時所有參數的取值。為了簡化求導過程，一般將原始問題轉化為對數似然的最大化問題，對（15）式兩邊取對數

$ln {p(X | m{pi}, m{mu}, Sigma)} = sum _ {n = 1} ^ {N} {ln { sum _ {k = 1} ^ {K} {pi_k N(m{x}_n | m{mu}_k, Sigma_k)} }}$ （17）

我們不妨先看看（17）式取極大值時參數應該滿足什麼條件。

先讓（17）式對 $m{mu}_k$ 求導並令導數為0：

$frac{partial P(X | m{pi}, m{mu}, Sigma)}{partial {m{mu}_k}} = -sum_{n = 1} ^ {N} {frac{pi_k N(m{x}_n | m{mu}_k, Sigma_k)}{sum_{j = 1} ^ {K} {pi_j N(m{x}_n | m{mu}_j, Sigma_j)}} Sigma_k^{-1} (m{x}_n - m{mu}_k)} = 0$ （18）

上式最外層求和符號裡面的第一項分式恰好就是由（14）式定義的後驗概率，僅僅多一個下標 $n$ ，可以簡記為 $gamma(z_{nk})$ ，帶入上式容易求解得到

$m{mu}_k = frac{sum_n{gamma(z_{nk})m{x}_n}}{sum_n{gamma(z_{nk})}} = frac{1}{N_k}sum_n{gamma(z_{nk})m{x}_n}$ （19）

其中，定義了

$N_k = sum_n{gamma(z_{nk})}$ （20）

由於 $gamma(z_{nk})$ 可以理解為數據點 $m{x}_n$ 來自第 $k$ 個高斯分布的可能性，或者第 $k$ 個高斯分布對生成該數據點承擔的「責任（responsibility）」，因此 $N_k$ 可以理解為所有數據點中來自第 $k$ 個高斯分布的數據點的「有效個數」。從而（19）式擁有良好的解釋性，它表示，第 $k$ 個高斯分布的均值 $m{mu}_k$ 等於所有數據點的加權平均，其中加權係數為該高斯分布應當為產生該數據點承擔的責任。

注意到，

$egin{align} & sum_k{gamma(z_{nk})} = 1 \ & sum_n{sum_{k} gamma(z_{nk})} = N end{align}$ （21）

上式的結論隨後將會用到。

然後讓（17）式對 $Sigma_k$ 求導並令導數為0，

$frac{partial P(X | m{pi}, m{mu}, Sigma)}{partial {Sigma_k}} = - frac{1}{2} sum _ {n = 1} ^ {N} {gamma(z_{nk}) left[ |Sigma_k| ^ {-1} frac{partial{|Sigma_k|}}{partial{Sigma_k}} + frac{partial{(m{x}_n - m{mu}_k) ^ T Sigma_k^{-1} (m{x}_n - m{mu}_k)}}{partial{Sigma_k}} ight]} = 0$ （22）

要求解該式，需要一些矩陣求導相關的知識。由於矩陣求導不是本文的重點，此處直接給出相關結論，對詳細推導過程感興趣的讀者可以參考相關資料。

$frac{partial{|A|}}{partial{A}} = |A| (A ^ {-1})T$ （23）

$frac{partial{A ^ {-1}}}{partial{x}} = -A^{-1} frac{partial{A}}{partial{x}} A^{-1}$ （24）

將（23）、（24）兩式帶入（22）式將其進一步簡化為

$- frac{1}{2} sum_{n = 1} ^ {N} {gamma(z_{nk}) Sigma_k^{-1} left[ I - (m{x}_n - m{mu}_k) (m{x}_n - m{mu}_k)^T Sigma_k^{-1} ight]} = 0$ （25）

求解得到

$Sigma_k = frac{1}{N_k} sum_n {gamma(z_{nk}) (m{x}_n - m{mu}_k) (m{x}_n - m{mu}_k)^T}$ （26）

與（19）式類似，上式可解釋為，第 $k$ 個高斯分布的協方差矩陣 $Sigma_k$ 等於所有數據點以 $m{mu}_k$ 為中心的協方差矩陣的加權平均，其中加權係數為該高斯分布應當為產生該數據點承擔的責任。

最後，讓（17）式對 $pi_k$ 求極大值。因為 $pi_k$ 滿足（7）式的歸一化約束，因此，在原目標函數的基礎上增加一個拉格朗日乘子，

$ln{p(X | m{pi}, m{mu}, Sigma) } + lambda (sum_k{pi_k} - 1)$ （27）

讓上式對 $pi_k$ 求導並令導數為0：

$sum_{n = 1} ^ {N} {frac{N(m{x}_n | m{mu}_k, Sigma_k)}{sum_{j = 1} ^ {K} {pi_j N(m{x}_n | m{mu}_j, Sigma_k)}}} + lambda = 0$ （28）

等式兩邊同時乘以 $pi_k$ 得

$sum_n{gamma(z_{nk})} + pi_k lambda = 0$ （29）

上式兩邊對 $k$ 求和，並利用（7）式和（21）式的結論可得

$lambda = -N$ （30）

帶入（29）式得

$pi_k = frac{N_k}{N}$ （31）

因此上式可以解釋為，混合係數 $pi_k$ 等於第 $k$ 個高斯分布為產生所有數據點承擔的責任的平均值。

注意到，（19）、（26）、（31）式並不構成該最大似然問題的解析解，因為 $m{mu}_k$ 、 $Sigma_k$ 、 $pi_k$ 的求解依賴 $gamma(z_{nk})$ ，而 $gamma(z_{nk})$ 的求解又反過來依賴 $m{mu}_k$ 、 $Sigma_k$ 、 $pi_k$ ，詳見（14）式。不過這很容易讓我們想到用迭代的方法求得參數值。第一步，隨機選取 $m{mu}_k$ 、 $Sigma_k$ 、 $pi_k$ 的初始值，依據（14）式求得 $gamma(z_{nk})$ 。第二步，將第一步求得的 $gamma(z_{nk})$ 依次帶入（19）、（26）、（31）式更新 $m{mu}_k$ 、 $Sigma_k$ 、 $pi_k$ 的值。如此往複迭代，直至達到收斂條件（似然值或參數值的變化在某個閾值以下）。因為每一步的迭代都會增加似然值，因此該迭代演算法可以保證找到似然函數的一個極大值，但不能保證找到全局最大值。迭代演算法的收斂速度很大程度上受初始值選取的影響，為了加快演算法的收斂速度，一般先用K均值演算法求得 $m{mu}_k$ 、 $Sigma_k$ 、 $pi_k$ 的初始值，然後以此初始值開始迭代。

下圖是用包含兩個高斯分布的混合高斯模型對old faithful數據集進行概率密度建模的示例圖。圖中，紅色（R）圓圈和藍色（B）圓圈分別為兩個高斯分布在一個標準差處的等高線，數據點用紅色和藍色的混合表示，紅色的多少用該數據點屬於紅色高斯分布的概率表示，藍色的多少用該數據點屬於藍色高斯分布的概率表示。因此，如果一個點完全屬於藍色，則該點將用藍色畫出，如果完全屬於紅色，則用紅色畫出，如果屬於紅色和藍色的概率恰好相等，則用紫色畫出。如圖（1），初始時，隨機設置兩個高斯分布的均值和方差，由於此時還未計算數據點屬於哪個高斯分布的概率，因此所有數據點用藍色（G）畫出（紅色和藍色通道為0）。如圖（2），經過一輪EM迭代後，兩個高斯分布的均值和方差均得到了更新，且每個數據點屬於兩個高斯分布的概率也計算出來了。可以看到，靠近藍色高斯分布中心的數據點幾乎完全為藍色，靠近紅色高斯分布中心的數據點幾乎完全為紅色，而處在兩個高斯分布中間的那些點帶著紫色，說明這些點屬於哪個高斯分布還不好確定。此時兩個高斯分布的方差還比較大，對數據點的概率密度的擬合還比較模糊。如圖（3），經過5輪EM迭代後，兩個高斯分布的均值和方差開始逐漸收斂。本例設置的收斂條件為數據集的對數似然值的變化量小於1e-4，經過11輪EM迭代後演算法收斂，最終收斂結果如圖（4）所示。圖（5）給出了數據集的對數似然隨迭代次數的變化趨勢，可見對數似然值在開始幾輪迭代中迅速增加，隨後到達平台期逐步上升最終達到收斂條件。本例中參數的初始值是隨機選取的，因此收斂速度較慢。如果先用K均值演算法找到一組「差不多」的參數，並以此作為初始參數開始迭代，可以大大加快收斂速度。文末給出了本例所使用的數據集和matlab代碼。

圖（1）初始化狀態

圖（2）經過第1論EM迭代後，兩個高斯分布的均值和方差均得到了更新。

圖（3）經過5輪EM迭代後，模型已經能較好地擬合觀測數據。

圖（4）經過11輪EM迭代後達到收斂。

圖（5）數據集的對數似然隨迭代次數的變化趨勢。

至此，我們還沒有引入EM演算法求解混合高斯模型的最大似然問題，不過已經很接近了。EM演算法也是一種迭代演算法，每一步迭代分E步和M步兩步。從接下來的介紹將會看到，上述迭代演算法的第一步和第二步恰好對應EM演算法的E步和M步。

五、EM演算法

此節將會脫離具體問題背景，介紹一般化的EM演算法。EM演算法用於求解帶隱藏變數的模型的最大似然問題。EM的全拼為，Expectation Maximization，即期望最大化演算法。

為此首先按照上文類似的方式引入幾個記號：用 $N imes D$ 的矩陣 $X$ 表示所有的觀測變數，其中矩陣的第 $n$ 行為 $m{x}_n^T$ ，數據點的維數為 $D$ ；用 $N imes K$ 的矩陣 $Z$ 表示所有的隱藏變數，其中矩陣的第 $n$ 行為 $m{z}_n^T$ ；用 $heta$ 表示所有的模型參數。因此，對數似然函數為

$ln {p(X | heta)} = ln {{sum_{Z}{p(X, Z | heta)}}}$ （32）

觀察上式可以發現，對數符號直接作用在聯合分布的和上，而不是直接作用於聯合分布。這使得要求得上式最優化問題的解析解幾乎是不可能的。

假設給定數據集 $X$ ，對應的 $Z$ 是知道的，那麼 $X$ 和 $Z$ 的聯合分布的對數似然函數為 $ln {p(X, Z | heta)}$ ，一般來說，最大化該似然函數的解析解是比較好求的。我們把 ${X, Z}$ 叫做完全數據集，而把單獨的 ${ X}$ 的叫做不完全數據集。

但實際中，我們只知道 $X$ ， $Z$ 是不知道的，所有關於 $Z$ 的信息只能從其後驗分布 $p(Z | X, heta)$ 中得知，因此我們無法求得全數據的似然。不過我們可以考慮全數據的似然在後驗分布 $p(Z | X, heta)$ 下的期望，此即期望最大化演算法中「期望」二字的來歷，對應EM演算法的E步。然後，我們將此期望值最大化，就可以求得一組參數值，此即期望最大化演算法中「最大化」一詞的來歷，對應EM演算法的M步。從而，EM演算法可被描述為包含如下兩步的迭代過程。在E步，用參數的當前值 $heta^{old}$ 求得隱藏變數的後驗分布 $p(Z | X, heta^{old})$ ；在M步，用此後驗分布計算全數據的似然的期望，將此期望記作 $Q( heta, heta^{old})$ ，則

$Q( heta, heta^{old}) = sum_{Z} {p(Z | X, heta^{old}) ln p(X, Z | heta)}$ （33）

然後求得此期望取最大值時的參數

$heta^{new} = argmax_{ heta} {Q( heta, heta^{old})}$ （34）

注意到，（33）式中，對數符號直接作用於聯合分布 $p(X, Z | heta)$ 上，因此上述最大化問題一般比較容易求得解析解。

將新求得的 $heta^{new}$ 作為下一步的 $heta^{old}$ 重複以上步驟，直至演算法收斂。

讀到此處，相信的大家心中都有一個疑問，為什麼取全數據對數似然的期望並最大化就可以最大化（32）式的似然呢？這個操作看上去似乎毫無依據啊？難道隨便蒙了一個恰好就蒙對了？當然不是蒙的，而且這樣選擇是有明確的理論依據的，問題的答案要在下文證明EM演算法時給出。此時，我們先給出EM演算法一般化的演算法框架。

EM演算法： 給定全數據的似然 $p(X, Z | heta)$ 表達式，目的是求得參數 $heta$ 使得觀測數據的似然 $p(X | heta)$ 最大。 1. 選擇參數的初始值，記為 $heta^{old}$ ； 2. E步：計算隱藏變數的後驗分布 $p(Z | X, heta^{old})$ ； 3. M步：最大化全數據的對數似然在 $p(Z | X, heta^{old})$ 下的期望，從而得到新的參數值 $heta^{new}$ : $heta^{new} = argmax_{ heta} {Q( heta, heta^{old})}$ 其中， $Q( heta, heta^{old})$ 為 $Q( heta, heta^{old}) = sum_{Z} {p(Z | X, heta^{old}) ln p(X, Z | heta)}$ 4. 檢查迭代是否收斂（似然值或參數值的變化量是否在給定閾值以下）；如果收斂，則停止，演算法結束；否則，令 $heta^{old} leftarrow heta^{new}$ ，返回第2步。

六、EM視角下的混合高斯模型求解

有了上一節的知識基礎，本節我們應用EM演算法框架重新求解混合高斯模型的最大似然問題，從而加深對EM演算法的理解。可以看到，採用EM演算法將會得到與第四節相同的結論。

EM演算法要求計算全數據的對數似然。由式（9）和（11）可以得到，全數據的似然函數

$egin{align} p(X, Z | m{pi}, m{mu}, Sigma) &= prod_{n = 1} ^ {N} {p(m{x}_n, m{z}_n | m{pi}, m{mu}, Sigma)} \ &= prod_{n = 1} ^ {N} {p(m{x}_n | m{z}_n, m{mu}, Sigma) p(m{z}_n | m{pi})} \ &= prod_{n = 1} ^ {N} {prod_{k = 1} ^ {K} {left[ pi_k N(m{x}_n | m{mu}_k, Sigma_k) ight] ^ {z_{nk}}}} end{align}$ （35）

對兩邊取對數得

$ln {p(X, Z | m{pi}, m{mu}, Sigma)} = sum_{n = 1} ^ {N} {sum_{k = 1} ^ {K} { z_{nk} left[ lnpi_k + ln N(m{x}_n | m{mu}_k, Sigma_k) ight]}}$ （36）

因此，全數據的對數似然在後驗分布 $p(Z | X, heta^{old})$ 下的期望為

$egin{align} Q( heta, heta^{old}) &= sum_{Z} {p(Z | X, heta^{old}) sum_{n = 1} ^ {N} {sum_{k = 1} ^ {K} { z_{nk} left[ lnpi_k + ln N(m{x}_n | m{mu}_k, Sigma_k) ight]}}} \ &= sum_{n = 1} ^ {N} {sum_{k = 1} ^ {K} { sum_{Z} {p(Z | X, heta^{old})} z_{nk} left[ lnpi_k + ln N(m{x}_n | m{mu}_k, Sigma_k) ight]}} end{align}$ （37）

由於給定 $X$ 時各個隱藏變數之間相互獨立，因此

$egin{align} sum_{Z} {p(Z | X, heta^{old})} z_{nk} &= sum_{Z} {prod_{n = 1} ^ {N} {p(m{z}_n | m{x}_n, heta^{old}) z_{nk}}} \ &= sum_{m{z}_n} {p(m{z}_n | m{x}_n, heta^{old}) z_{nk}} \ &= p(m{z}_{nk} = 1 | m{x}_n, heta^{old}) \ &= gamma(z_{nk}) end{align}$ （38）

從而，

$Q( heta, heta^{old}) = sum_{n = 1} ^ {N} {sum_{k = 1} ^ {K} { gamma(z_{nk}) left[ lnpi_k + ln N(m{x}_n | m{mu}_k, Sigma_k) ight]}}$ （39）

對上式求極大值就比較容易了，可以直接得到解析解。首先，讓上式對 $m{mu}_k$ 求導並令導數為0，可以得到與式（19）相同的結論。然後，讓上式對 $Sigma_k$ 求導並令導數為0，可以得到與式（26）相同的結論。最後，考慮 $pi_k$ 的歸一化約束條件，增加相應的拉格朗日乘子後令導數為0，可以得到與式（31）相同的結論。

七、EM視角下的K均值聚類演算法

前文分別介紹了K均值演算法和混合高斯模型，本節將會看到K均值演算法其實是混合高斯模型的一個特例，只需要對混合高斯模型增加一個簡單的約束即可。因此，也可以將EM演算法應用到混合高斯模型的結論應用於K均值演算法，可以看到，用EM演算法求解K均值聚類問題的結果與第二節介紹的直接最小化聚類誤差得到的結果一模一樣。

從前文的介紹可以看到，K均值演算法和混合高斯模型其實有一些相似之處。例如，都假設模型由 $K$ 個「組件」構成，K均值聚類假設有 $K$ 個聚類，混合高斯模型假設有 $K$ 個高斯分布；都從產生式的角度將數據點分配給各個組件，只不過K均值聚類將每個數據點唯一分配給其中一個組件，具體分配給哪個組件由 $gamma_{nk}$ 指定（（2）式），而混合高斯模型認為每個數據點可以從屬於每個高斯分布，而從屬的程度由後驗分布 $gamma(z_{nk})$ 給定。如果說K均值聚類是「硬」分配（hard assignment）的話，那麼混合高斯模型則是「軟」分配（soft assignment）。

在混合高斯模型中，如果假設所有高斯分布的協方差矩陣相同，且都為 $epsilon{I}$ ，其中 $epsilon$ 為已知的常數，則此時

$gamma(z_{nk}) = frac{pi_kexp{-||m{x}_n-m{mu}_k||^2/(2epsilon)}}{sum_{j = 1} ^ {K} {pi_j exp{ -||m{x}_n - m{mu}_j|| ^ 2 / (2epsilon)}}}$ （40）

現在讓 $epsilon o 0$ ，則上式分子分母中的指數項都將趨近於0，但趨近的速度不一樣，其中 $||m{x}_n - m{mu}_j||^2$ 最小的那一項趨近的速度最慢。因此，只有當 $k = arg min_{j} {||m{x}_n - m{mu}_j||^2}$ 時， $gamma(z_{nk}) = 1$ ， $k$ 為其他值時 $gamma(z_{nk})$ 均為0。即

$gamma(z_{nk}) = { egin{align} &1, quad if k = arg min _{j} {||m{x}_n - m{mu}_j||^2} \ &0, quad otherwise end{align}$ （41）

此式與（2）式一模一樣。可見，此時 $gamma(z_{nk}) o gamma_{nk}$ ，即混合高斯模型的軟分配退化成了K均值聚類的硬分配。同時，（19）式的均值更新方程退化成（4）式的聚類中心更新方程。從而可以看到，混合高斯模型退化成了K均值聚類演算法。

八、EM演算法的證明

前文一直在講EM演算法就是先求隱藏變數的後驗分布，然後最大化全數據的對數似然在此後驗分布下的期望，這樣不斷迭代就可以得到EM演算法的解。但是相信讀者心中還帶著疑問，為什麼要選擇全數據的對數似然？為什麼要最大化全數據的對數似然在後驗分布下的期望？為什麼這樣迭代後卻會收斂到觀測數據的似然 $p(X | heta)$ 的極大值？

要回答以上問題很簡單，不過需要對觀測數據的似然做一個非常巧妙的變換：

$egin{align} ln p(X | heta) &= sum_{Z} {q(Z) ln p(X | heta)} \ &= sum_{Z}{q(Z) ln frac{p(X,Z | heta)}{p(Z|X, heta)}} \ &= sum_Z{q(Z) ln frac{p(X, Z | heta)}{q(Z)}} + sum_{Z} {q(Z)lnfrac{q(Z)}{p(Z | X, heta)}} \ &= L(q, heta) + KL(q || p) end{align}$ （42）

其中，

$L(q, heta) = sum_Z{q(Z) ln frac{p(X, Z | heta)}{q(Z)}}$ （43）

$KL(q || p) = sum_{Z} {q(Z)lnfrac{q(Z)}{p(Z | X, heta)}}$ （44）

因此，將目標函數 $ln p(X | heta)$ 分解成了 $L(q, heta)$ 和 $KL(q||p)$ 兩項。注意到， $L(q, heta)$ 是 $heta$ 的函數，是 $q$ 的泛函，因為 $q$ 本身是一個函數。 $KL(q||p)$ 是分布 $q$ 和 $p$ 之間的KL距離，也是一個泛函，它具有非負的特性，即 $KL(q||p) ge 0$ ，當且僅當 $q = p$ 時， $KL(q||p) = 0$ 。此分解的示意圖如下：

因此，根據（42）式的分解，可以推導出最大化 $ln p(X | heta)$ 的迭代演算法，即EM演算法。

在E步，假設參數為已知的常數，記為 $heta^{old}$ ，因此（42）式的左邊為一個常數，而等式右邊分解成兩項，這兩項的大小根據 $q$ 的取值不同而不同。因為 $KL(q||p) ge 0$ ，所以 $L(q, heta)$ 是 $ln p(X | heta)$ 的下界，我們要讓這個下界盡量大。顯然，只有當 $KL(q||p) = 0$ 即 $q(Z) = p(Z | X, heta ^ {old})$ 時，該下界最大，且正好等於 $ln p(X| heta^{old})$ 。如下圖，紅色線的高度表示等式左邊的值 $ln p(X | heta^{old})$ ，此時它是一個常量。藍色虛線為分解的位置，虛線下面是 $L(q, heta^{old})$ 的值，虛線上面是 $KL(q||p)$ 的值。根據 $q$ 的取值不同藍色虛線可以上下移動，當 $q(Z) = p(Z | X, heta ^ {old})$ 時，藍色虛線移動到與紅色線等高，此時下界 $L(q, heta^{old})$ 達到最大。

在接下來的M步，固定 $q(Z)$ 的取值，以 $heta$ 為參數最大化 $L(q, heta)$ ，得到新的參數 $heta^{new}$ 。因為 $L(q, heta)$ 是 $ln p(X | heta)$ 的下界，下界增大了， $ln p(X | heta)$ 也會跟著增大。且此時因為 $heta^{new} eq heta^{old}$ ，從而 $q eq p$ ，所以 $KL(q||p) gt 0$ 。因此 $ln p(X | heta)$ 增大的量比下界增大的量更大。如下圖，紅色虛線和藍色虛線分別是經過E步之後目標函數及其下界的值，此時二者相等。在經過M步之後，下界增大到藍色實線的位置，同時 $KL(q||p)$ 不再等於0，所以目標函數從紅色虛線增大到紅色實線的位置，增大的量等於兩項分解式增大的量的和。

將 $q(Z) = p(Z | X, heta ^ {old})$ 帶入（43）式，可以看到我們在M步要最大化的下界變為

$egin{align} L(q, heta) &= sum_{Z} {p(Z | X, heta^{old}) {p(X, Z | heta)}} - sum_{Z} {p(Z | X, heta^{old}) {p(Z | X, heta^{old})}} \ &= Q( heta, heta^{old}) + const end{align}$ （45）

因此，在M步要最大化的實際是全數據的對數似然在後驗分布 $p(Z | X, heta ^ {old})$ 下的期望，這正是前文在用EM演算法求解混合高斯模型時給出的結論。

因為每經過一輪E步和M步，都會增大目標函數的值，因此當迭代收斂時，找到的一定是目標函數的極大值。

EM演算法的核心步驟在於求解M步的最大化問題，當 $x$ 和 $z$ 的聯合分布屬於指數分布族時，對數正好和指數相消，可以較容易求得解析解。但當聯合分布是其他複雜的分布時，M步往往不好直接求解，此時需要對標準EM演算法進行擴展，其中GEM（Generalized EM）和ECM（Expectation conditional maximization）演算法即是其中的例子。在數據量比較大無法一次性載入內存進行計算時，還可以將標準EM演算法轉化為序列化演算法，一次只需要載入一個數據點對參數進行序列化地更新。除此之外，EM演算法還可以用來求解參數 $heta$ 的MAP問題。這些內容由於不在本文的範圍之內，就不再贅述了。此文完。

附matlab源碼

用包含兩個高斯分布的混合高斯模型對old faithful數據集進行概率密度建模的matlab源碼如下：

% main.mclearclcload old_faithful[sof, z_mu, z_sigma] = zscore(old_faithful); % scaled old faithful data[N, D] = size(sof);K = 2;X = sof;% initialize. You can also initialize parameters using K-means algorithm.pie = [0.5, 0.5];mu = [-1.0, 0; 1 1.5];sigma(:, :, 1) = [0.5, 0; 0, 0.5];sigma(:, :, 2) = [0.5, 0; 0, 0.5];index = 1;% EM iterationswhile 1 % plot before E step (the label of each point is unknown and so ploted % in green color. figure plot(X(:, 1), X(:, 2), g.); hold on h = ezplot(@(x,y)myfunc(x, y, mu(1, :), inv(sigma(:, :, 1)))); set(h, color, b, lineWidth, 2) h = ezplot(@(x,y)myfunc(x, y, mu(2, :), inv(sigma(:, :, 2)))); set(h, color, r, lineWidth, 2) title(sprintf(iteration %d, before E step., index)); axis equal axis([-2 2 -3 3]) hold off % E step -- calculate responsibilities R = zeros(N, K); for k = 1 : K R(:, k) = mvnpdf(X, mu(k, :), sigma(:, :, k)) * pie(k); end LL(index) = sum(log(sum(R, 2))); % log likelihood R = R ./ repmat(sum(R, 2), 1, K); % M step -- update parameters Nk = sum(R); for k = 1 : K mu(k, = sum(diag(R(:, k)) * X) / Nk(k); MU = repmat(mu(k, :), N, 1); sigma(:, :, k) = (X - MU) * diag(R(:, k)) * (X - MU) / Nk(k); end pie = Nk / N; % plot after M step figure hold on for i = 1 : N plot(X(i, 1), X(i, 2), ., color, [R(i, 2) 0 R(i, 1)]); end h = ezplot(@(x,y)myfunc(x, y, mu(1, :), inv(sigma(:, :, 1)))); set(h, color, b, lineWidth, 2) h = ezplot(@(x,y)myfunc(x, y, mu(2, :), inv(sigma(:, :, 2)))); set(h, color, r, lineWidth, 2) title(sprintf(iteration %d, after M step, index)); axis equal axis([-2 2 -3 3]) hold off mu sigma pie LL(index) if index > 1 && LL(index) - LL(index - 1) <= 1e-4 break; end index = index + 1;end

function [ z ] = myfunc( x, y, mu, inv_sigma )%UNTITLED Summary of this function goes here% Detailed explanation goes herez = (x - mu(1)) .* inv_sigma(1, 1) .* (x - mu(1)) ... + (x - mu(1)) .* inv_sigma(1, 2) .* (y - mu(2)) ... + (y - mu(2)) .* inv_sigma(2, 1) .* (x - mu(1)) ... + (y - mu(2)) .* inv_sigma(2, 2) .* (y - mu(2)) - 1;end

old faithful數據集可以從網上公開的地址下載，或者從本人以下網盤地址下載：

鏈接: https://pan.baidu.com/s/1nw0kNFV 密碼: ujnb