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

EM算法及其應用（一）

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

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

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}\}\)：

1. 固定\(\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. \]

2. 固定$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算法的優缺點：

優點：

1. 可解釋性比較強。
2. 調參的參數僅為簇數k。
3. 相對於高斯混合模型而言收斂速度快，因而常用於高斯混合模型的初始值選擇。K-means 的時間復雜度為 \(\mathcal{O}(N\cdot K \cdot I)\) ，簇數 \(K\) 和 迭代次數 \(I\) 通常遠小於\(N\)，所以可優化為 \(\mathcal{O}(N)\) ，效率較高。

缺點：

1. 對離群點敏感。

2. K值難以事先選取，交叉驗證不大適合，因為簇越多，目標函數 \(\sum\limits_{n=1}^N r_{nk} ||\mathbf{x}_n - \boldsymbol{\mu}_k||^2\) 就越小。常采用的方法有：一、“拐點法”，如下圖 K=3 就是一個拐點。

二、 加上正則化系數 \(\lambda\) ，使得 \(\sum\limits_{n=1}^N \left(r_{nk} ||\mathbf{x}_n - \boldsymbol{\mu}_k||^2\right) + \lambda K\) 最小。
3. 無法保證收斂到全局最優值，常使用不同的初始值進行多次試驗。也可以通過 K-means++ 算法優化，核心思想是選取與已有質心距離較遠的點作為初始值。

4. 只能發現球狀的簇。

5. 由於采用歐氏距離，無法直接計算類別型變量。

高斯混合模型

---

高斯混合模型同樣多用於聚類，與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}\)

1. 初始化各個分模型參數 \(\pi_k, \boldsymbol{\mu}_k, \boldsymbol{\Sigma}_k\)
2. 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*} \]

3. 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} $$
4. 重復2. 和3. 步直至收斂。

最后總結一下高斯混合模型的優缺點：

優點：

1. 相比於K-means更具一般性，能形成各種不同大小和形狀的簇。K-means可視為高斯混合聚類中每個樣本僅指派給一個混合成分的特例，且各混合成分協方差相等，均為對角矩陣\(\sigma^2\mathbf{I}\)。
2. 僅使用少量的參數就能較好地描述數據的特性。

缺點：

1. 高斯混合模型的計算量較大收斂慢。因此常先對樣本集跑k-means，依據得到的各個簇來定高斯混合模型的初始值。其中質心即為均值向量，協方差矩陣為每個簇中樣本的協方差矩陣，混合系數為每個簇中樣本占總體樣本的比例。
2. 分模型數量難以預先選擇，但可以通過划分驗證集來比較。
3. 對異常點敏感。
4. 數據量少時效果不好。

EM算法在半監督學習上的應用

在現實的分類問題中常遇到數據集中只有少量樣本是帶有標簽的，而大部分樣本的標簽缺失。如果直接將無標簽的樣本丟棄，則會容易造成大量的有用信息被浪費。而如果使用EM算法，將缺失標簽視為隱變量，很多時候能達到利用所有數據進行訓練的目的。流程如下：

1. 僅用帶標簽樣本訓練學習器，得到初始參數\(\theta\)。
2. E步： 利用訓練好的學習器預測無標簽樣本，將其分類到概率最大的類別，這樣所有樣本就都有標簽了。
3. M步： 用所有樣本重新訓練學習器，得到參數\(\theta^t\)。
4. 重復2. 和 3. 步直至收斂，得到最終模型。

/