EM算法及其應用： K-means 與高斯混合模型

時間 2019-12-11

原文原文鏈接

EM算法及其應用（一）

EM算法及其應用（二）： K-means 與高斯混合模型

上一篇闡述了EM算法的主要原理，這一篇來看其兩大應用 —— K-means 與高斯混合模型，主要由EM算法的觀點出發。html

K-means

K-means的目標是將樣本集劃分爲K個簇，使得同一個簇內樣本的距離儘量小，不一樣簇之間的樣本距離儘量大，即最小化每一個簇中樣本與質心的距離。K-means屬於原型聚類(prototype-based clustering)，原型聚類指聚類結構能經過一組原型刻畫，而原型即爲樣本空間中具備表明性的點。在K-means中，這個原型就是每一個簇的質心 $\boldsymbol{\mu}$ 。

從EM算法的觀點來看，K-means的參數就是每一個簇的質心 $\boldsymbol{\mu}$，隱變量爲每一個樣本的所屬簇。若是事先已知每一個樣本屬於哪一個簇，則直接求平均便可獲得 $\boldsymbol{\mu}$ 。但現實中不知道的狀況下，則須要運用EM的思想：

假設要k個簇，先隨機選定k個點做爲質心$\{\boldsymbol{\mu_1}, \boldsymbol{\mu_2} \cdots \boldsymbol{\mu_k}\}$：算法

固定$\boldsymbol{\mu_k}$，將樣本劃分到距離最近的$\boldsymbol{\mu_k}$所屬的簇中。若用$r_{nk}$表示第n個樣本所屬的簇，則
\[ r_{nk} = \left. \begin{cases} 1 \,\, & \text{if} \;\;\;k = \mathop{argmin}_j ||\mathbf{x}_n - \boldsymbol{\mu}_j||^2 \\ 0 \,\, & \text{otherwise} \end{cases} \right. \]
固定$r_{nk}$，根據上一步的劃分結果從新計算每一個簇的質心。因爲咱們的目標是最小化每一個簇中樣本與質心的距離，可將目標函數表示爲 $J = \sum\limits_{n=1}^N r_{nk} ||\mathbf{x}_n - \boldsymbol{\mu}_k||^2$，要最小化$J$則對$\boldsymbol{\mu}_k$求導得 $2\sum\limits_{n=1}^N r_{nk}(\mathbf{x}_n - \boldsymbol{\mu}_k) = 0$，則
\[ \boldsymbol{\mu}_k = \frac{\sum_nr_{nk} \mathbf{x}_n}{\sum_n r_{nk}} \]
即簇中每一個樣本的均值向量。函數

上面兩步分別更新$r_{nk}$和$\boldsymbol{\mu_k}$就對應了EM算法中的E步和M步。和EM算法同樣，K-means每一步都最小化目標函數 $J$，於是能夠保證收斂到局部最優值，但在非凸目標函數的狀況下不能保證收斂到全局最優值。另外，K-means對每一個樣本進行「硬分配(hard assignment)」，即只歸類爲一個簇，然而某些樣本可能處於簇與簇的邊界處，將這些樣本強行分到其中一個簇可能並不能反映確信度。後文的高斯混合模型可進行「軟分配(soft assignment)」，即對每一個樣本所劃分的簇進行機率估計。學習

最後總結一下K-means算法的優缺點：優化

優勢：spa

可解釋性比較強。
調參的參數僅爲簇數k。
相對於高斯混合模型而言收斂速度快，於是經常使用於高斯混合模型的初始值選擇。K-means 的時間複雜度爲 $\mathcal{O}(N\cdot K \cdot I)$ ，簇數 $K$ 和迭代次數 $I$ 一般遠小於$N$，因此可優化爲 $\mathcal{O}(N)$ ，效率較高。

缺點：prototype

對離羣點敏感。3d
K值難以事先選取，交叉驗證不大適合，由於簇越多，目標函數 $\sum\limits_{n=1}^N r_{nk} ||\mathbf{x}_n - \boldsymbol{\mu}_k||^2$ 就越小。常採用的方法有：1、「拐點法」，以下圖 K=3 就是一個拐點。htm

2、加上正則化係數 $\lambda$ ，使得 $\sum\limits_{n=1}^N \left(r_{nk} ||\mathbf{x}_n - \boldsymbol{\mu}_k||^2\right) + \lambda K$ 最小。
沒法保證收斂到全局最優值，常使用不一樣的初始值進行屢次試驗。也能夠經過 K-means++ 算法優化，核心思想是選取與已有質心距離較遠的點做爲初始值。blog
只能發現球狀的簇。
因爲採用歐氏距離，沒法直接計算類別型變量。

高斯混合模型

高斯混合模型一樣多用於聚類，與K-means不一樣的是其採用機率模型來表達聚類原型。
首先回顧一下高斯分佈(Gaussian Distribution)：對於隨機變量$x$，其機率密度函數可表示爲：
\[ \mathcal{N}(x|\mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x-\mu)^2}{2\sigma^2}} \]
若$\mathbf{x}$爲n維隨機向量，則多元高斯分佈(Multivariate Gaussian Distribution)爲：
\[ \mathcal{N}(\mathbf{x}| \boldsymbol{\mu},\mathbf{\Sigma}) = \frac{1}{(2\pi)^{\frac{n}{2}}|\Sigma|^\frac12}\,e^{-\frac12(\mathbf{x}-\boldsymbol{\mu})^T\mathbf{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})} \]
其中$\boldsymbol{\mu}$爲n維均值向量，$\mathbf{\Sigma}$爲$n\times n$的協方差矩陣，$|\mathbf{\Sigma}|$爲$\mathbf{\Sigma}$的行列式。

不少時候咱們發現單個高斯分佈沒法很好地描述數據的性質，以下圖數據分紅了兩個簇，若是使用兩個高斯分佈明顯能更好地描述其結構。

所以沿着這個思路就誕生了高斯混合模型(Mixture of Gaussians)，本質上是k個高斯分佈的線性組合，這樣靈活性大增，能描述更多樣的分佈：
\[ p(\mathbf{x}) = \sum\limits_{k=1}^{K}\pi_k\mathcal{N}(\mathbf{x}|\boldsymbol{\mu}_k,\mathbf{\Sigma}_k) \qquad \tag{1.1} \]
其中$0 \leqslant\pi_k\leqslant 1$爲混合係數(mixing coefficients)，知足$\sum\limits_{k=1}^{K} \pi_k= 1$ 。

因爲本文的主角是EM算法，因此下面以EM算法的角度來看待高斯混合模型。

回憶EM算法是一種對含有隱變量的機率模型參數的極大似然估計。一般隱變量$\mathbf{Z}$未知，而實際知曉的只有觀測數據$\mathbf{X}$。而對於高斯混合模型來講，觀測數據是這樣產生的：先依照$\pi_k$選擇第k個高斯分模型，而後依照其機率密度函數進行採樣生成相應的樣本。

能夠看到在通過上述過程後咱們拿到手的只有採樣得來的樣本，殊不知道每一個樣原本自於哪一個分模型，但這樣就表示咱們得到了高斯混合模型的隱變量。

這裏定義K維隨機向量 $\mathbf{z} = \begin{bmatrix} z_1 \\ z_2 \\ \vdots \\ z_k \\ \end{bmatrix}$，$\mathbf{z}$中只有一個$z_k$爲1，其他元素爲0，即$z_k = \left. \begin{cases}1\,, & \text{數據來自第k個分模型} \\ 0\,, & \text{不然} \end{cases} \right.$，

這樣$z_k$即爲高斯混合模型的隱變量，表示樣原本自於第k個分模型。

因爲$p(\mathbf{x}) = \sum\limits_zp(\mathbf{x},\mathbf{z}) = \sum\limits_{z}p(\mathbf{z})\,p(\mathbf{x}|\mathbf{z})$，對比$(1.1)$式中高斯混合模型的定義，可將$\pi_k$視爲選擇第k個分模型的先驗機率，即$\pi_k = p(z_k = 1)$；而對應的$\mathcal{N}(\mathbf{x}|\mathbf{\mu}_k, \mathbf{\Sigma}_k) = p(\mathbf{x}|z_k = 1)$。另外在獲得觀測數據$\{\mathbf{x}_1, \mathbf{x}_2 \cdots \mathbf{x}_n\}$後，每一個$\mathbf{x}_n$都對應一個隱變量$z_{nk}$，則運用貝葉斯定理，$\mathbf{x}_n$屬於第k個分模型的後驗機率 (記做$\gamma(z_{nk})$)爲：
\[ \begin{align*} \gamma(z_{nk}) = p(z_{nk} = 1|\mathbf{x}_n) & = \frac{p(z_{nk} = 1)\,p(\mathbf{x}_n|z_{nk} = 1)}{\sum\limits_{j=1}^Kp(z_{nj}=1)\,p(\mathbf{x}_n|z_{nj} = 1)} \\ & = \frac{\pi_k\,\mathcal{N}(\mathbf{x}_n|\mathbf{\mu}_k,\mathbf{\Sigma}_k)}{\sum\limits_{j=1}^K \pi_j\,\mathcal{N}(\mathbf{x}_n|\mathbf{\mu}_j,\mathbf{\Sigma}_j)} \tag{1.2} \end{align*} \]

下圖顯示了$\gamma(z_{nk})$的做用，圖a是依照徹底數據的聯合分佈$p(\mathbf{x},\mathbf{z})$做圖，每一個點依照其所屬的第k個分模型標記顏色，這相似於「硬分配」；圖b則不考慮隱變量$\mathbf{z}$，僅依照不徹底數據$\mathbf{x}$直接做圖；圖c則考慮了每一個樣原本自於各個分模型的後驗機率$\gamma(z_{nk})$，這樣一些在簇中間的點的顏色是紅藍綠三種顏色的混合，代表這些點來自於每一個簇的機率比較接近，這即爲「軟分配」。

爲了估計參數，若是直接採用極大似然法，則
\[ \begin{align*} L(\mathbf{\theta}) = ln\,p(\mathbf{X}|\mathbf{\theta}) = ln\,p(\mathbf{X}|\mathbf{\pi},\boldsymbol{\mu},\mathbf{\Sigma}) & = ln\left[\prod\limits_{n=1}^N\left(\sum\limits_{k=1}^K\,\pi_k \mathcal{N}(\mathbf{x}_n|\boldsymbol{\mu}_k,\mathbf{\Sigma}_k)\right) \right] \\[2ex] & = \sum\limits_{n=1}^N ln\left(\sum\limits_{k=1}^K \pi_k \frac{1}{(2\pi)^{\frac{n}{2}}|\Sigma|^\frac12}\,e^{-\frac12(\mathbf{x}-\boldsymbol{\mu})^T\mathbf{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})} \right) \end{align*} \]
上式直接求導比較複雜，由於存在「和的對數」 $ln (\sum_{k=1}^K\pi_k\,\mathcal{N} \cdot)$ ，而若是使用EM算法則會方便不少。

先依照上一篇EM算法的流程，寫出含Q函數：
\[ \mathcal{Q}(\mathbf{\theta}, \mathbf{\theta}^{(t)}) = \sum\limits_{n=1}^N\sum\limits_{\mathbf{z}}p(\mathbf{z}|\mathbf{x}_n,\mathbf{\pi}^{(t)},\boldsymbol{\mu}^{(t)},\mathbf{\Sigma}^{(t)})\,ln(\mathbf{x}_n,\mathbf{z}|\pi,\boldsymbol{\mu},\mathbf{\Sigma}) \tag{1.3} \]
由$(1.2)$式可知，$(1.3)$中第一項 $p(z_{nk} = 1|\mathbf{x}_n,\pi^{(t)},\boldsymbol{\mu}^{(t)},\mathbf{\Sigma}^{(t)}) = \gamma(z_{nk})$，表示當前參數下每一個樣本$\mathbf{x}_n$屬於第k個分模型的後驗機率。而第二項爲徹底數據的對數似然函數：
\[ \begin{align*} ln(\mathbf{x}_n,\mathbf{z}|\pi,\mathbf{\mu},\mathbf{\Sigma}) & = ln\prod\limits_{k=1}^K \left[\pi_k\mathcal{N}(\mathbf{x}_n|\boldsymbol{\mu}_k,\mathbf{\Sigma_k}) \right]^{z_{nk}} \\ & = \sum\limits_{k=1}^Kz_{nk} \left[ln \pi_k + ln\,\mathcal{N}(\mathbf{x_n}|\boldsymbol{\mu}_k,\mathbf{\Sigma}_k)\right] \end{align*} \]
由此$(1.3)$式變爲：
\[ \begin{align*} & \sum\limits_{n=1}^N\sum\limits_{k=1}^K \gamma (z_{nk}) \left[ln \pi_k + ln\,\mathcal{N}(\mathbf{x_n}|\boldsymbol{\mu}_k,\mathbf{\Sigma}_k)\right] \\ = \;& \sum\limits_{n=1}^N\sum\limits_{k=1}^K \gamma (z_{nk}) \left[ln\,\pi_k + ln\left(\frac{1}{(2\pi)^\frac{n}{2}|\mathbf{\Sigma|^\frac{1}{2}}}\,e^{-\frac12 (\mathbf{x}_n - \boldsymbol{\mu}_k)^T\mathbf{\Sigma}_k^{-1}(\mathbf{x}_n-\boldsymbol{\mu}_k)} \right) \right] \tag{1.4} \end{align*} \]
能夠看到上式括號中變成了「對數的和」，這使得求導方便了不少，且求解時式中$ln(\cdot)$和$e^{(\cdot)}$能夠相互抵消。

$(1.4)$式分別對$\boldsymbol{\mu}_k$，$\boldsymbol{\Sigma}_k$求導得：
\[ \boldsymbol{\mu}_k = \frac{\sum\limits_{n=1}^N\gamma(z_{nk}) \cdot \mathbf{x}_n}{\sum\limits_{n=1}^N\gamma(z_{nk})}\qquad,\qquad \boldsymbol{\Sigma}_k = \frac{\sum\limits_{n=1}^N \gamma(z_{nk})\cdot (\mathbf{x}_n - \boldsymbol{\mu}_k)\cdot(\mathbf{x}_n-\boldsymbol{\mu}_k)^T}{\sum\limits_{n=1}^N \gamma(z_{nk})} \]
能夠看到分模型k的$\boldsymbol{\mu}_k$和$\boldsymbol{\Sigma}_k$是全部樣本的加權平均，其中每一個樣本的權重爲該樣本屬於分模型k的後驗機率$\gamma(z_{nk})$ 。對比上文中K-means的 $\boldsymbol{\mu}_k = \frac{\sum_nr_{nk} \cdot \mathbf{x}_n}{\sum_nr_{nk}}$，兩者形式相似，區別爲K-means中的$r_{nk} \in \{0,1\}$，而高斯混合模型中$\gamma(z_{nk})$爲機率估計。

對於混合係數$\pi_k$，由於有限制條件$\sum_{k=1}^K \pi_k = 1$，所以須要使用拉格朗日乘子法轉化爲無約束優化問題：
\[ \sum\limits_{n=1}^N\sum\limits_{k=1}^K \gamma (z_{nk}) \left[ln\,\pi_k + ln\left(\frac{1}{(2\pi)^\frac{n}{2}|\mathbf{\Sigma|^\frac{1}{2}}}\,e^{-\frac12 (\mathbf{x}_n - \boldsymbol{\mu}_k)^T\mathbf{\Sigma}_k^{-1}(\mathbf{x}_n-\boldsymbol{\mu}_k)} \right) \right] + \lambda(\sum\limits_{k=1}^K \pi_k - 1) \]
對$\pi_k$求導得，
\[ \begin{align*} & \sum\limits_{n=1}^N\frac{\gamma(z_{nk})}{\pi_k} + \lambda = 0 \qquad \Longrightarrow \qquad \pi_k = \frac{\sum\limits_{n=1}^N \gamma(z_{nk})}{-\lambda} \tag{1.5} \\ & \\ & \text{因爲$\sum_{k=1}^K\pi_k = 1\;，\qquad$ 則} \;\;\sum_{k=1}^K \;\frac{\sum_{n=1}^N \gamma (z_{nk})}{-\lambda} = 1 \;, \;\lambda = -N \quad 代入(1.5) \;,\\ & \\ & \text{得}\; \pi_k = \frac{\sum_{n=1}^N\gamma(z_{nk})}{N} \end{align*} \]
即每一個分模型k的混合係數是屬於k的樣本的平均後驗機率，由此運用EM算法能大大簡化高斯混合模型的參數估計過程，在中間步只需計算$\gamma(z_{nk})$就好了。

高斯混合模型的算法流程

輸入：樣本集 $\mathbf{X} = \{\mathbf{x}_1, \mathbf{x}_2 \cdots \mathbf{x}_n\}$

輸出：高斯混合模型參數 $\pi, \boldsymbol{\mu}, \boldsymbol{\Sigma}$

初始化各個分模型參數 $\pi_k, \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k$
E步：依照當前模型參數，計算觀測數據$\mathbf{x}_n$屬於分模型k的後驗機率：
\[ \begin{align*} \gamma(z_{nk}) = p(z_{nk} = 1|\mathbf{x}_n) = \frac{\pi_k\,\mathcal{N}(\mathbf{x}_n|\mathbf{\mu}_k,\mathbf{\Sigma}_k)}{\sum\limits_{j=1}^K \pi_j\,\mathcal{N}(\mathbf{x}_n|\mathbf{\mu}_j,\mathbf{\Sigma}_j)} \end{align*} \]
M步：根據 $\gamma(z_{nk})$計算新一輪參數：
\[ \boldsymbol{\mu}_k^{new} = \frac{\sum\limits_{n=1}^N\gamma(z_{nk}) \cdot \mathbf{x}_n}{\sum\limits_{n=1}^N\gamma(z_{nk})} \]
\[ \boldsymbol{\Sigma}_k^{new} = \frac{\sum\limits_{n=1}^N \gamma(z_{nk})\cdot (\mathbf{x}_n - \boldsymbol{\mu}_k^{new})\cdot(\mathbf{x}_n-\boldsymbol{\mu}_k^{new})^T}{\sum\limits_{n=1}^N \gamma(z_{nk})} \]
\[ \pi_k^{new} = \frac{\sum_{n=1}^N\gamma(z_{nk})}{N} \]
重複2. 和3. 步直至收斂。