斯坦福大學機器學習——EM演算法求解高斯混合模型

阿新 • • 發佈：2019-01-26

EM演算法（Expection-Maximizationalgorithm，EM）是一種迭代演算法，通過E步和M步兩大迭代步驟，每次迭代都使極大似然函式增加。但是，由於初始值的不同，可能會使似然函式陷入區域性最優。辜麗川老師和其夫人發表的論文：基於分裂EM演算法的GMM引數估計（提取碼：77c0）改進了這一缺陷。下面來談談EM演算法以及其在求解高斯混合模型中的作用。

一、高斯混合模型（Gaussian MixtureModel， GMM）

之前寫過高斯判別分析模型，利用引數估計的方法用於解決二分類問題。下面介紹GMM，它是對高斯判別模型的一個推廣，也能借此引入EM演算法。

假設樣本集為 $\{x^{(1)},x^{(2)},...,x^{(m)}\}$ 並且樣本和標籤滿足聯合分佈

$p(x^{(i)},z^{(i)})=p(x^{(i)}|z^{(i)})p(z^{(i)})$ 。這裡： $z^{(i)}$ 服從多項式分佈，即 $z^{(i)}\sim Multinomial(\phi)$ （ $\phi_{j}=p(z^{(i)}=j)$ ， $\phi_{j}\ge 0$ ， $\underset{j=1}{\overset{k}{\sum}}\phi_{j}=1$ ），且 $z\in\{1,2,...,k\}$ ；在 $z^{(i)}$ 給定的情況下， $x^{(i)}$ 服從正態分佈，即 $x^{(i)}|z^{(i)}=j\sim N(\mu_{j},\Sigma_{j})$ 。這樣的模型稱為高斯混合模型。

該模型的似然函式為：

$\begin{aligned}l(\phi,\mu,\Sigma)&=\underset{i=1}{\overset{m}{\sum}}log\; p(x^{(i)};\phi,\mu,\Sigma)\\ &= \underset{i=1}{\overset{m}{\sum}}log \underset{z^{(i)}=1}{\overset{k}{\sum}}p(x^{(i)}|z^{(i)};\mu,\Sigma)p(z^{(i)};\phi) \end{aligned}$

如果直接令 $l(\phi,\mu,\Sigma)$ 的各變數偏導為0，試圖分別求出各引數，我們會發現根本無法求解。但如果變數 $z^{(i)}$ 是已知的，求解便容易許多，上面的似然函式可以表示為：

$l(\phi,\mu,\Sigma)=\underset{i=1}{\overset{m}{\sum}}log\;p(x^{(i)}|z^{(i)};\mu,\Sigma)+log\;p(z^{(i)};\phi)$

利用偏導求解上述式，可分別得到引數 $\phi,\mu,\Sigma$ 的值：

$\phi_{j}=\frac{1}{m}\underset{i=1}{\overset{m}{\sum}}\#\{z^{(i)}=j\}$

$\mu_{j}=\frac{\underset{i=1}{\overset{m}{\sum}}\#\{z^{(i)}=j\}x^{(i)}}{\underset{i=1}{\overset{m}{\sum}}\#\{z^{(i)}=j\}}$

$\Sigma_{j}=\frac{\underset{i=1}{\overset{m}{\sum}}\#\{z^{(i)}=j\}(x^{(i)}-\mu_{j})(x^{(i)}-\mu_{j})^T}{\underset{i=1}{\overset{m}{\sum}}\#\{z^{(i)}=j\}}$

其中,#{ }為指示函式，表示滿足括號內條件的數目。

那麼，變數 $z^{(i)}$ 無法通過觀察直接得到， $z^{(i)}$ 就稱為隱變數，就需要通過EM演算法，求解GMM了。下面從Jensen不等式開始，介紹下EM演算法：

二、 Jensen不等式（Jensen’s inequality）

引理：如果函式f的定義域為整個實數集，並且對於任意x或

$\overset{\rightarrow}{x}$ 存在 $f''(x)\ge 0$ 或函式的Hessian矩陣 $H\ge0$ ，那麼函式f稱為凹函式。 $f''(x)$ 或函式的Hessian矩陣H>0，那麼函式f為嚴格凹函式。

（存在 $f''(x)\le 0$ 或函式的Hessian矩陣 $H\le0$ ，那麼函式f稱為凸函式；如果 $f''(x)$ 或函式的Hessian矩陣 H<0，那麼函式f為嚴格凸函式。）

定理：如果函式f是凹函式，X為隨機變數，那麼：

$E[f(x)]\ge f(EX)$

不幸的是很多人都會講Jensen不等式記混，我們可以通過圖形的方式幫助記憶。下圖中，橫縱座標軸分別為X和f(X)，f(x)為一個凹函式，a、b分別為變數X的定義域，E[X]為定義域X的期望。圖中清楚的看到各個量的位置和他們間的大小關係。反之，如果函式f是凸函式，X為隨機變數，那麼：

$E[f(x)]\le f(EX)$

Jensen不等式等號成立的條件為：

$E[X]=X$ ，即X為一常數。

三、 EM演算法

假設訓練集 $\{x^{(1)},x^{(2)},...,x^{(m)}\}$ 是由m個獨立的樣本構成。我們的目的是要對 $p(x,z)$ 概率密度函式進行引數估計。它的似然函式為：

$\begin {aligned}l(\theta)&=\underset{i=1}{\overset{m}{\sum}}log\;p(x,\theta)\\ &=\underset{i=1}{\overset{m}{\sum}}log\underset{z}{\sum}p(x,z;\theta) \end{aligned}$

然而僅僅憑藉似然函式，無法對引數進行求解。因為這裡的隨機變數 $z^{(i)}$ 是未知的。

EM演算法提供了一種巧妙的方式，可以通過逐步迭代逼近最大似然值。下面就來介紹下EM演算法：

假設對於所有i， $Q_{i}$ 皆為隨機變數 $z^{(i)}$ 的分佈函式。即： $\sum_{Z}Q_{i}(z)=1,Q_{i}(z)\ge0$ 。那麼：

$\begin{aligned}l(\theta)&=\underset{i}{\sum}log \; p(x^{(i)};\theta)\\ &=\underset{i}{\sum}log\underset{z^{(i)}}{\sum}p(x^{(i)},z^{(i)};\theta)\quad\quad\quad\quad &(1)\\&=\underset{i}{\sum}log\underset{z^{(i)}}{\sum}Q_{i}(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}&(2)\\ &\ge \underset{i}{\sum}\underset{z^{(i)}}{\sum}Q_{i}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}&(3) \end{aligned}$

其中第（2）步至第（3）步的推導就使用了Jensen不等式。其中：f(x)=log x， $f''(x)=-\frac{1}{x^{2}}<0$ ，因此為凸函式； $\underset{z^{(i)}}{\sum}Q_{i}(z^{(i)})\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}$ 表示隨機變數為 $Q_{i}(z^{(i)})$ 概率分佈函式為 $\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}$ 的期望。因此有：

$f(\underset{z^{(i)}\sim Q_{i}}{E}[\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}])\ge \underset{z^{(i)}\sim Q_{i}}{E}[f(\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})})]$

這樣，對於任意分佈 $Q_{i}$ ，（3）都給出了 $l(\theta)$ 的一個下界。如果我們現在通過猜測初始化了一個 $\theta$ 的值，我們希望得到在這個特定的 $\theta$ 下，更緊密的下界，也就是使等號成立。根據Jensen不等式等號成立的條件，當 $\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}$ 為一常數時，等號成立。即：

$\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})} =c$

由上式可得 $p(x^{(i)},z^{(i)};\theta)=cQ_{i}(z^{(i)})$ ，又 $\sum_{Z}Q_{i}(z)=1$ ，因此 $p(x^{(i)},z^{(i)};\theta)=c$ 。再由上式可得：

$\begin {aligned}Q_{i}(z^{(i)})&=\frac{p(x^{(i)},z^{(i)};\theta)}{\sum_{z}{p(x^{(i)},z;\theta)}}\\&=\frac{p(x^{(i)},z^{(i)};\theta)}{{p(x^{(i)};\theta)}}\\&=p(z^{(i)}|x^{(i)};\theta) \end{aligned}$

上述等式最後一步使用了貝葉斯公示。

EM演算法有兩個步驟：

（1）通過設定初始化 $\theta$ 值，求出使似然方程最大的 $Q_{i}$ 值，此步驟稱為E-步（E-step）

（2）利用求出的 $Q_{i}$ 值，更新 $\theta$ 。此步驟稱為M-步（M-step）。過程如下：

repeat until convergence{

(E-step) for each i, set

$Q_{i}(z^{(i)}):=p(z^{(i)}|x^{(i)};\theta)$ $Q_{i}(z^{(i)})=p(z^{(i)}|x^{(i)};\theta)$

(M-step) set

$\theta:=arg\,\underset{\theta}{max}\underset{i}{\sum}\underset{z^{(i)}}{\sum}Q_{i}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}$

｝

那麼，如何保證EM演算法是收斂的呢？下面給予證明：

假設 $\theta^{(t)}$ 和 $\theta^{(t+1)}$ 是EM演算法第t次和第t+1次迭代所得到的引數 $\theta$ 的值，如果有 $l(\theta^{(t)})\le l(\theta^{(t+1)})$ ，即每次迭代後似然方程的值都會增大，通過逐步迭代，最終達到最大值。以下是證明：

$\begin{aligned} l(\theta^{(t+1)}) &\ge \underset{i}{\sum}\underset{z^{(i)}}{\sum}Q_{i}^{(t)}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta^{(t+1)})}{Q_{i}^{(t)}(z^{(i)})}\quad\quad&(4)\\ &\ge \underset{i}{\sum}\underset{z^{(i)}}{\sum}Q_{i}^{(t)}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta^{(t)})}{Q_{i}^{(t)}(z^{(i)})}&(5)\\ &=l(\theta^{(t)})&(6) \end {aligned}$

不等式（4）是由不等式（3）得到，對於任意 $Q_{i}$ 和 $\theta$ 值都成立；得到不等式（5）是因為我們需要選擇特定的 $\theta^{(t+1)}$ 使得方程 $arg\,\underset{\theta}{max}\underset{i}{\sum}\underset{z^{(i)}}{\sum}Q_{i}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}$ 在 $\theta^{(t+1)}$ 處的值大於在 $\theta^{(t)}$ 處的值；等式（6）是找到特定的 $\theta^{(t)}$ 的值，使得等號成立。

最後我們通過圖形的方式再更加深入細緻的理解EM演算法的特點：

$J(Q,\theta)=\underset{i}{\sum}\underset{z^{(i)}}{\sum}Q_{i}(z^{(i)})log\frac{p(x^{(i)},z^{(i)};\theta)}{Q_{i}(z^{(i)})}$

由上文我們知道有這樣的關係： $l(\theta)\ge J(Q,\theta)$ ，EM演算法就是不斷最大化這個下界，逐步得到似然函式的最大值。如下圖所示：

首先，初始化 $\theta^{(0)}$ ，調整 $Q(z)$ 使得 $J(Q,\theta^{(0)})$ 與 $l(\theta^{(0)})$ 相等，然後求出 $J(Q,\theta^{(0)})$ 使得到最大值的 $\theta^{(1)}$ ；固定 $\theta^{(1)}$ ，調整 $Q(z)$ ，使得 $J(Q,\theta^{(1)})$ 與 $l(\theta^{(1)})$ 相等，然後求出使 $J(Q,\theta^{(1)})$ 得到最大值的 $\theta^{(2)}$ ；……；如此迴圈，使得 $l(\theta)$ 的值不斷上升，直到k次迴圈後，求出了 $l(\theta)$ 的最大值 $l(\theta^{(k)})$ 。

四、 EM演算法應用於混合高斯模型（GMM）

再回到GMM：上文提到由於隱變數 $z^{(i)}$ 的存在，無法直接求解引數，但可以通過EM演算法進行求解：

E-Step:

$\omega_{j}^{(i)}=Q_{i}(z^{(i)}=j) =P(z^{(i)}=j|x^{(i)};\phi,\mu,\Sigma)$ M-Step:
$\begin{aligned} &\underset{i=1}{\overset{m}{\sum}}\underset{z^{(i)}}{\sum}Q_{i}(z^{(i)})log \frac{p(x^{(i)},z^{(i)};\phi,\mu,\Sigma)}{Q_{i}(z^{(i)})}\\ = &\underset{i=1}{\overset{m}{\sum}}\underset{j=1}{\overset{k}{\sum}}Q_{i}(z^{(i)})log \frac{p(x^{(i)}|z^{(i)}=j;\mu,\Sigma)p(z^{(i)}=j;\phi)}{Q_{i}(z^{(i)})}\\ = &\underset{i=1}{\overset{m}{\sum}}\underset{j=1}{\overset{k}{\sum}}\omega_{j}^{(i)}log\frac{\frac{1}{(2\pi)^{\frac{n}{2}}\left|\Sigma_{j}\right|^\frac{1}{2}}exp(-\frac{1}{2}(x^{(i)}-\mu_{j})^{T}\Sigma_{j}^{-1}(x^{(i)}-\mu_{j}))\cdot\phi_{j}}{\omega_{j}^{(i)}} \end {aligned}$

（1）引數 $\mu$
對期望 $\mu$ 的每個分量 $\mu_{l}$ 求偏導：

$\begin{aligned}&\triangledown_{\mu_{l}}\underset{i=1}{\overset{m}{\sum}}\underset{j=1}{\overset{k}{\sum}}\omega_{j}^{(i)}log\frac{\frac{1}{(2\pi)^{\frac{n}{2}}\left|\Sigma_{j}\right|^\frac{1}{2}}exp(-\frac{1}{2}(x^{(i)}-\mu_{j})^{T}\Sigma_{j}^{-1}(x^{(i)}-\mu_{j}))\cdot\phi_{j}}{\omega_{j}^{(i)}}\\ =&-\triangledown_{\mu_{l}}\underset{i=1}{\overset{m}{\sum}}\underset{j=1}{\overset{k}{\sum}}\omega_{j}^{(i)}\frac{1}{2}(x^{(i)}-\mu_{j})^{T}\Sigma_{j}^{-1}(x^{(i)}-\mu_{j})\\= &\frac{1}{2}\underset{i=1}{\overset{m}{\sum}}\omega_{l}^{(i)}\Delta_{\mu_{l}}(2\mu_{l}^{T}\Sigma_{l}^{-1}x^{(i)}-\mu_{l}^{T}\Sigma_{l}^{-1}\mu_{l})\\=&\underset{i=1}{\overset{m}{\sum}}\omega_{l}^{(i)}\Sigma_{l}^{-1}(x^{(i)}-\mu_{l}) \end{aligned}$

令上式為0，得：

$\mu_{l}:=\frac{\sum_{i=1}^{m}\omega_{l}^{(i)}x^{(i)}}{\sum_{i=1}^{m}{\omega_{l}^{(i)}}}$

（2）引數

觀察M-Step，可以看到，跟 $\phi_{i}$ 相關的變數僅僅有 $\omega_{j}^{(i)}$ 。因此，我們僅僅需要最大化下面的目標函式：

$max\,\underset{i=1}{\overset{m}{\sum}}\underset{j=1}{\overset{k}{\sum}}\omega_{j}^{(i)}log\phi_{j}$

又由於 $\sum_{j=1}^{k}\phi_{j}=1$ ，為約束條件。因此，可以構造拉格朗日運算元求目標函式：

$\mathcal L(\phi)=\underset{i=1}{\overset{m}{\sum}}\underset{j=1}{\overset{k}{\sum}}\omega_{j}^{(i)}log\phi_{j}+\beta(\underset{j=1}{\overset{k}{\sum}}\phi_{j}-1)$

求偏導：

$\frac{\partial}{\partial \phi_{j}} \mathcal L(\phi)=\underset{i=1}{\overset{m}{\sum}}\frac{\omega_{j}{(i)}}{\phi_{j}}+1$

令 $\frac{\partial}{\partial \phi_{j}} \mathcal L(\phi)=0$ 得：

$\phi_{j}:=\frac{\underset{i=1}{\overset{m}{\sum}}\omega_{j}^{(i)}}{-\beta}$

將 $\sum_{j=1}^{k}\phi_{j}=1$ 帶入上式得：

$\begin{aligned} -\beta=\sum_{i=1}^{m}\sum_{j=1}^{k}\omega_{j}^{(i)} &= \sum_{i=1}^{m}1 &=m\end{aligned}$

最後將 $\beta=-m$ 帶入得：

$\phi_{j}:=\frac{\underset{i=1}{\overset{m}{\sum}}\omega_{j}^{(i)}}{m}$

（3）引數 $\Sigma$

$\begin{aligned}&\triangledown _{\Sigma_{j}} \underset{i=1}{\overset{m}{\sum}}\underset{j=1}{\overset{k}{\sum}}\omega_{j}^{(i)}log\frac{\frac{1}{(2\pi)^{\frac{n}{2}}\left|\Sigma_{j}\right|^\frac{1}{2}}exp(-\frac{1}{2}(x^{(i)}-\mu_{j})^{T}\Sigma_{j}^{-1}(x^{(i)}-\mu_{j}))\cdot\phi_{j}}{\omega_{j}^{(i)}}\\ =&\triangledown_{\Sigma_{j}}\underset{i=1}{\overset{m}{\sum}}\underset{j=1}{\overset{k}{\sum}} \omega_{j}^{(i)}[-\frac{n}{2}log2\pi-\frac{1}{2}log\left|\Sigma_{j}\right|-\frac{1}{2}(x^{(i)}-\mu_{j})^{T}\Sigma_{j}^{-1}(x^{(i)}-\mu_{j})+log\phi_{j}-log\omega_{j}^{(i)}]\\=&\underset{i=1}{\overset{m}{\sum}}\underset{j=1}{\overset{m}{\sum}}\omega_{j}^{(i)}\frac{1}{2}\Sigma_{j}^{-T}+\frac{1}{2}(x^{(i)}-\mu_{j})^{T}\Sigma^{-2}(x^{(i)}-\mu_{j}) \end{aligned}$

令上式為零，解得：

$\Sigma_{j}:=\frac{\sum_{i=1}^{m}\omega_{j}^{(i)} (x^{(i)}-\mu_{j})(x^{(i)}-\mu_{j})^{T} }{ \sum_{i=1}^{m}\omega_{j}^{(i)} }$

五、總結

EM演算法利用不完全的資料，進行極大似然估計。通過兩步迭代，逐漸逼近最大似然值。而GMM可以利用EM演算法進行引數估計。

最後提下辜老師論文的思路：EM模型容易收斂到區域性最大值，並且嚴重依賴初試值。傳統的方法即上文中使用的方法是每次迭代過程中，同時更新高斯分佈中所有引數，而辜老師的方法是把K個高斯分佈中的一個分量，利用奇異值分解的方法將其分裂為兩個高斯分佈，並保持其他分量不變的情況下，對共這K+1個高斯分佈的權值進行更新，直到符合一定的收斂條件。這樣一來，雖然演算法複雜度沒有降低，但每輪只需要更新兩個引數，大大降低了每輪迭代的計算量。

斯坦福大學機器學習——EM演算法求解高斯混合模型

斯坦福大學機器學習——EM演算法求解高斯混合模型

【機器學習】EM演算法在高斯混合模型學習中的應用

機器學習讀書筆記（高斯混合模型GMM與EM）（改）

EM演算法在高斯混合模型中的應用（詳細解釋與求解）

EM演算法與高斯混合模型

EM演算法之高斯混合模型(二)

EM演算法之高斯混合模型（一）

R語言：EM演算法和高斯混合模型的R語言實現

EM演算法和高斯混合模型GMM介紹

【機器學習】機器學習（十二、十三）：K-means演算法、高斯混合模型

EM(期望最大演算法)在高斯混合模型中的python實現

斯坦福大學機器學習筆記——異常檢測演算法（高斯分佈、多元高斯分佈、異常檢測演算法）

斯坦福大學機器學習筆記——聚類（k-均值聚類演算法、損失函式、初始化、聚類數目的選擇）

斯坦福大學機器學習第一課“引言(Introduction)”

吳恩達斯坦福大學機器學習 CS229 課程學習筆記（一）

斯坦福大學機器學習——誤差理論（Error Theory）

斯坦福大學機器學習——交叉驗證（Cross Validation）

高斯混合模型（GMM）及其求解（期望最大化（EM）演算法）

斯坦福大學機器學習筆記——推薦系統（協同過濾、低秩分解、推薦系統）

Machine Learning（Stanford）| 斯坦福大學機器學習筆記--第一週（1.監督學習與無監督學習的介紹）

斯坦福大學機器學習——EM演算法求解高斯混合模型

相關推薦