IA9a - K-means分群與EM算法: 理論

市場內分門別類的產品

原則上K-means分群是EM算法 (expectation-maximization algorithm) 的一個特例，相較於EM從似然函數 (likelihood) 的方式出發並導出不同資料點對不同群 (或稱聚類、簇，cluster) 的貢獻度，K-means單純用歐氏距離來分類，雖然這加快了K-means算法的計算效率，但顯然的這麼naïve的方法會產生明顯的缺陷。但仍讓我們會先從簡單的K-means開始，並探討K-means算法的結果在某些情況下和人類肉眼辨識的結論會產生明顯的歧義，接著再推進到通用的EM算法，除了python的demo外，我也會證明在數學上如何從EM得到K-means的極限 (僅限高斯混合的情況)，而在這個part內，讓我們先從數學理論的角度來面。

一、K-means分群

今天如果我有一大坨資料共120個點如圖1所示，從資料的散佈趨勢來看，很明顯的分成兩群，那我們是不是可以有一個演算法來自動化將資料點如同我們臆測的一樣分成兩袋？一個最簡單的發法便是使用K-means分群，而這個大寫K表示群的數量，以圖1而言，我們猜測K=2，所以是個2-means分群問題。

圖1. 資料散佈圖 (普魯士藍) 與初始化的紅綠兩點

　所以為了達成分群，我們需要先初始化兩個點，這兩個點代表兩群的中心。事實上，初始化的位置不必要一定剛好在兩群的中心，可以是任意挑選的地方，但為了計算方便，我們會將起始化的兩點開始在靠近資料點的平均值附近，如同圖1中紅綠兩點。我們將這兩個初始化的點叫作$p_k$其中$k=1,2$表示第$k$群。另外將資料點的集合用$\mathbf{D}=\{d_1,d_2\cdots d_n \cdots d_{120}\}$表示，所以我們計算每個點$d_n=(x_{d_n},y_{d_n})$和$p_k=(x_{p_k},y_{p_k})$的歐氏距離 $$\ell_{n,k}=\sqrt{(x_{d_n}-x_{p_k})^2+(y_{d_n}-y_{p_k})^2}. \tag{1}$$ 若該點和$p_1$的距離小於和$p_2$的距離 ($\ell_{n,1}<\ell_{n,2}$)，則我們說$d_n$比較靠近$p_1$且應該和$p_1$分在一起。當所有的點都輪過一次Eq. (1)的計算後，應該會分出兩群了，我們在分別計算出這兩群的平均位置並將之重新賦與給$p_{1,2}$並完成一次迭帶的結果，如圖2。

圖2. 經過一次迭代後的結果

　不難發現淺綠色和淡橘色的兩群資料點便是經過Eq. (1)後認定應該是分在一起的群組，然後透過這兩群我們重新計算這兩群的平均點，新的$p_k$已經和圖1的位置不一樣了。我們將新的$p_k$代入Eq. (1)同樣的重複迭代過程，直到紅綠兩點都不再變化位置為止，最終的結果如圖3，這是經過3次迭代後的結果，現在每個資料點都已經找出和它最靠近的群了。由於這個過程非常簡單，我們就不列出代碼細節，可以參考隨附的python腳本 (會在part 2和EM算法一起給)，但我們接下來會看到K-means算法會產生一個和人類視覺相抵觸的結果，而這顯然是K-means分群的謬誤。

圖3. 經歷三次迭代後紅綠兩個資料群的平均點皆不再移動位置

1.1 錯的離譜！

為了容易比較，我們生成兩群一樣120個資料點，分別從均值為 $$\mu_{1} =(0,1),\\ \mu_{2} =(2.7,3.8)$$ 和協方差矩陣 (covariance matrix) 為 $$\Sigma_{1} =\left(\begin{array}{cc} 0.1 & 0.084\\ 0.084 & 7 \end{array}\right),\\ \Sigma_{2} =\left(\begin{array}{cc} 0.1 & -0.142\\ -0.142 & 9 \end{array}\right)$$ 的高斯分佈中所生成，當然電腦只會看到它們merge在一起而不是兩個分開的群。接著我們一樣給定兩個初始點$p_k$如圖4。

圖4. 由高斯分佈生成的兩群資料點與初始化點$p_k$

　很合理的認為，通過Eq. (1)的K-means算法，最終我們應該可以分出左右兩群，況且初始化點也離我們的視覺所告訴我們群的位置，這應該不是個難事。我們通過3次迭代後如圖5所示，這個結果非常令人跌破眼鏡，一個對人類來說如此容易辨別的分群，電腦居然出現除此離譜的謬誤。K-means並沒有將資料點分成左右兩群，而是將它分成上下兩群，如果你以為在多幾次迭代結果會自動修正，那麼很遺憾的只是徒勞而已，K-means truncated在這個莫名其妙的地方了！

圖5. 經過5次迭代K-means分群的結果

　但如果我們真的仔細檢視每一次迭代得出對每個$d_n$的$\ell_{n,k}$，我們會發現K-means其實非常正確的分出了該$d_n$比較靠近$p_1$或$p_2$，所以數學上如圖5的結果並沒有違背我們對K-means的認知，所以可見K-means並不善於分類這種狹長形不對稱結構的資料分佈。也許你會從圖3的比較猜484這個K-means比較能分上下而不分左右，答案是no-no，見圖6，我們將$\Sigma_k$的對角元素對調 (原本朝$x$方向伸展的變成朝$y$方向伸展)，並重新run一次K-means，分類結果一樣不如預期啊！

圖6. 對調$\Sigma_k$的對角元素生成的資料點使用K-means分群經過5次迭代後的結果，一樣和我們的分群預測相抵觸

二、是你的還是我的？EM演算法

我們在上一節中明顯的看到，K-means透過歐氏距離$\ell$來區分每個資料點$d_n$到底比較靠近哪個$p_k$，如果在這次迭代中我發現某個點比較靠近第$p_1$ (即$\ell_{n,1}<\ell_{n,2}$，則它強迫$d_n$被分到第1群去，而最後再用所有被強迫分到第1群去的點的平均算出下一個$p_1$。這種方式就像IA8裡提到的ICM一樣，到底該取什麼值 (分到哪群去) 都是絕對的，所以只要assign的時候分錯 (這種錯不是電腦的錯，而是我們在implant算法的時候思考就錯了)，我又用錯誤的分群在計算下一次分群所需的$p_k$，最後K-means不斷使用錯誤的分群來做迭代，最終得到的分類結果當然也不會是正確的了，就像是被truncated在局域而非全域的最佳解上！
　所以我們希望每個點被分到哪個群不再是一個yes/no (binary系統) 的問題，而是一個機率決定的問題。比如一樣以一個K=2的情況，我們讓每個點都得到一個值$0\leq\gamma_n\leq 1$，愈小則表示離群1的中心值愈靠攏，愈大則離群2的中心愈靠攏，所以我們不再像之前一樣強迫分兩類，而是以比例原則來調控每個點到底比較有可能屬於哪個群。如果群與群之間有著明確的楚河漢界，那群的界線就會非常明確，若群與群彼此交融，那麼這條界線就會是一個漸進的過程，我們將在下面看到。在正式進入EM算法以前，讓我們先談一談高斯混合模型 (mixtures of Gaussians) 。

2.1 高斯混合模型

首先將資料集寫成$\mathbf{x}$ (其實就是前面的$\mathbf{d}$)，因為每個資料點都可以看成是從某個分佈中採樣的點所組成的隨機變量 (這有個專門稱呼讀術語叫獨立同分佈iid)。在繼續以前，讓我們先介紹一個有$K$個元素的向量$\mathbf{z}$，特別的是在這個向量裡面總是只有一個元素不為0，其值為1，考慮一個$K=6$ (或稱6維) 的例子我們有$\mathbf{z}=(0,0,1,0,0,0)$，所以對$\mathbf{z}$裡的每個元素$z_k \in \{0,1\}$，但又限製了只有一個不為0，我們又稱這種向量叫one-hot向量 (one-hot vector) 或1-of-K表示 (1-of-K representation)。

圖7. z和x之間的關係圖

　所以對一個$\mathbf{z}$，假設第$k$元素為1 (也透露了其它都是0) 的機率是 $$p(z_k) = \pi_k,\tag{2}$$ 再者，我們假設$\mathbf{z}$和隨機變量$\mathbf{x}$之間有一個如圖7所描述的關係。所以我們稱$\mathbf{z}$是$\mathbf{x}$的隱變數 (hidden/latent variable)，而聯合機率為 $$p(\mathbf{x},\mathbf{z})=p(\mathbf{x}|\mathbf{z})p(\mathbf{z})\tag{3}.$$ 而對一個$\mathbf{z}$而言，我們有 $$p(\mathbf{z})=\prod_k^{K}\pi_k ^{z_k}, \tag{4}$$ 強調一下Eq. (4)並不是likelihood function而是和Eq. (2)等價的描述。假定$\mathbf{z}$的第$k$個元素$z_k=1$且滿足$1\leq k\leq K$，那麼我們有Eq. (4)變為 $$p(\mathbf{z})=\pi_1^0 \pi_2^0 \cdots \pi_k^{1} \cdots \pi_{K-1}^0\pi_K^0=1\times 1\times\cdots\times \pi_k\times\cdots\times 1\times 1 =\pi_k$$ 正是Eq. (2)！假使給定$z_k=1$的條件生成$\mathbf{x}$的機率是 $$p(\mathbf{x}|z_k=1)=\mathcal{N}(\mathbf{x}|\mu_k,\Sigma_k) \tag{5}$$ 則等價為 $$p(\mathbf{x}|\mathbf{z})=\prod_{k=1}^K \mathcal{N}(\mathbf{x}|\mu_k,\Sigma_k)^{\pi_k}.\tag{6}$$ 這樣便完成大部分的數學工作了。
　由於我們只在意$\mathbf{x}$，所以我們要計算Eq. (3)的邊緣機率 $$p(\mathbf{x})=\sum_{\mathbf{z}}p(\mathbf{x}|\mathbf{z})p(\mathbf{z})\tag{7}$$ 來得到關於$\mathbf{x}$的所有訊息，通過Eqs. (4,6)展開後我們有 $$\begin{align*} \sum_{\mathbf{z}}p(\mathbf{x}|\mathbf{z})p(\mathbf{x}) & =\sum_{\mathbf{z}}\prod_{k}\mathcal{N}(\mathbf{x}|\mu_{k},\Sigma_{k})^{z_{k}}\prod_{k}\pi_{k}^{z_{k}}\\ & =\left.\mathcal{N}(\mathbf{x}|\mu_{1},\Sigma_{1})\pi_{1}\right|_{\mathbf{z}=(1,0,0....)}+\left.\mathcal{N}(\mathbf{x}|\mu_{2},\Sigma_{2})\pi_{2}\right|_{\mathbf{z}=(0,1,0....)}\\ & \quad+\left.\mathcal{N}(\mathbf{x}|\mu_{3},\Sigma_{3})\pi_{3}\right|_{\mathbf{z}=(0,0,1....)}+\cdots+\left.\mathcal{N}(\mathbf{x}|\mu_{K},\Sigma_{K})\pi_{K}\right|_{\mathbf{z}=(0...0,1)}\\ & =\sum_{k=1}^{K}\pi_{k}\mathcal{N}(\mathbf{x}|\mu_{k},\Sigma_{k}),\tag{8} \end{align*}$$ 故得 (上面那些$\prod_k$的求積符號當求和給出一個特定的比如說$z_2=1$的向量$\mathbf{z}=(0,1,0...0)$時就可以直接拿掉了，因為這些在相乘的因子除了$k=2$外都有一個0的次方，所以任何數的0次方都變成1了) $$p(\mathbf{x})=\sum_{k=1}^{K}\pi_{k}\mathcal{N}(\mathbf{x}|\mu_{k},\Sigma_{k})\tag{9}$$ 正是我們所要的高斯混合模型的表達式。Eq. (8)顯示了資料集$\mathbf{x}$是由$K$個均值$\mu_k$和斜方差矩陣$\Sigma_k$的高斯分佈所構成，每個分佈的權重是$\pi_k$的機率分佈$p(\mathbf{x})$。
　另一個有趣的分佈是若我們把圖7的箭頭反向，若已知$\mathbf{x}$可不可以逆推回去對應的哪個$k$有$z_k=1$的機率，實際上就是將Eq. (5)倒施逆行的意思，故我們通過貝氏定理有 $$p(z_k=1|\mathbf{x})=\frac{p(\mathbf{x}|z_k=1)p(z_k=1)}{p(\mathbf{x})}=\frac{\pi_k \mathcal{N}(\mathbf{x}|\mu_k,\Sigma_k)}{\sum_{k=1}^{K}\pi_{k}\mathcal{N}(\mathbf{x}|\mu_{k},\Sigma_{k})}\equiv \gamma(z_k)\tag{10}$$ 我們在當中使用了Eqs. (2,5,9)並定義該結果叫做$\gamma(z_k)$。在上式中$\pi_k$為$z_k=1$的先驗機率，而$\gamma(z_k)$則視為是當$\mathbf{x}$給定 (被觀測到) 時$z_k=1$的後驗機率。我們也稱$\gamma(z_k)$為成分$k$作為解釋$\mathbf{x}$的置信度 (responsibilities)。

2.2 EM演算法

OK，讓我們從Eq. (9)出發，整個式子表示了當給定參數$\mu_k$和$\Sigma_k$後它們生成$\mathbf{x}$的機率，若使$\mathbf{x}=\{x_n\}_{n=1}^N$且每筆資料都是從某個高斯混合中抽出的獨立同分佈樣本，我們有 $$p(x_n)=\sum_k \pi_k \mathcal{N}(x_n|\mu_k,\Sigma_k),$$ 故資料集$\mathbf{x}$的似然函數顯然是 $$L(\mathbf{x})=\prod_{n=1}^N p(x_n)=\prod_{n=1}^N \left[ \sum_k \pi_k \mathcal{N}(x_n|\mu_k,\Sigma_k)\right].$$ 將其取對數後得到 $$\ln L(\mathbf{x})=\sum_{n=1}^N \ln\left[ \sum_k \pi_k \mathcal{N}(x_n|\mu_k,\Sigma_k)\right].\tag{11}$$ 而最佳的參數$\boldsymbol{\theta}=\{\boldsymbol{\pi},\boldsymbol{\mu},\boldsymbol{\Sigma}\}$便是使得Eq. (11)有最大的組合，我們要使用maximum likelihood法 (可以參考PS6)。
　首先使Eq. (11)最大的$\mu_k$可計算 $$\frac{\partial \ln L}{\partial \mu_k}=-\sum_{n=1}^{N}\underbrace{\frac{\pi_{k}\mathcal{N}(x_{n}|\mu_{k},\Sigma_{k})}{\sum_{j}\pi_{j}\mathcal{N}(x_{n}|\mu_{j},\Sigma_{j})}}_{\gamma(z_{nk})}\Sigma_{k}(x_{n}-\mu_{k})=0.$$ 預設$\Sigma_k$是non-singular，我們可以將上式乘上$\Sigma_k^{-1}$得到 $$\mu_k=\frac{1}{N_k}\sum_{n=1}^N\gamma(z_{nk})x_n\tag{12}$$ 和 $$N_k=\sum_{n=1}^N \gamma(z_{nk})\tag{13}$$ 其中$\gamma(z_{nk})$解釋同Eq. (10)，多了個下標n表示結果是對特定某個資料點$x_n$最有可能的解釋。同理$\Sigma_k$和$\pi_k$都可以用相同的方法導出，再繼續給出答案以前，讓我們先回顧一下Eq. (12)。
　還記得一開始是假設資料點可以分成$K$個高斯分佈的群體，而$\mu_k$表示的是第$k$個高斯的均值點，但與K-means不同的是，這個$\mu_k$並不像K-means只考慮那些被assigned到第$k$個群裡的點才有貢獻，而是整個資料集的$N$個點都納入計算當中了。所以代表我們並沒有強迫哪個點被分類到哪個群，而是它們都是屬於整體分佈$p(\mathbf{x})$的一份子，我們改用在Eq. (12)裡的置信度$\gamma(z_{nk})$來決定資料點$x_n$有多大的貢獻會到第$k$個群裡，離$k$群近的置信度會比較遠的來得大以調控貢獻的程度。這樣折衷的方式較K-means來得有彈性，也避免落入將某次分錯結果的影響不斷傳播到往後的迭代中，我們稍後便會看到結果。
　接著對$\Sigma_k$和$\pi_k$的計算方法同上，只是將偏微分中的$\mu_k$換掉罷了，結果如下 $$\Sigma_k = \frac{1}{N_k}\sum_{n=1}^N \gamma(z_{nk})(x_n-\mu_k)\otimes (x_n-\mu_k)\tag{14}$$ 和 $$\pi_k = \frac{N_k}{N}.\tag{15}$$ 在Eq. (14)裡的符號$\otimes$表示外積 (outer product)，以2維平面空間為例，兩個向量A和B外積變成 $$\left(\begin{array}{c} a_{1}\\ a_{2} \end{array}\right)\otimes\left(\begin{array}{c} b_{1}\\ b_{2} \end{array}\right)=\left(\begin{array}{c} a_{1}\\ a_{2} \end{array}\right)(b_{1},b_{2})=\left(\begin{array}{cc} a_{1}b_{1} & a_{1}b_{2}\\ a_{2}b_{1} & a_{2}b_{2} \end{array}\right),$$ 而因為上式第2步的關係，也有人會寫成$\mathbf{A}\otimes\mathbf{B}=\mathbf{A}\mathbf{B}^T$，其中T表轉置 (transpose)。
　所以我們便完成了所有的計算啦！不難發現在求參數$\boldsymbol{\theta}$前我們都需要先知道置信度$\gamma(z_{nk})$，置信度告訴我們每個資料點$x_n$對第$k$個群的預期 (expectation) 貢獻是多少，所以首先求得$\gamma(z_{nk})$便是EM算法中的第E步 (E-step)，而有了置信度後，我們接著便可以導出使likelihood function最大化 (maximization) 的參數$\boldsymbol{\theta}$，這便是第M步 (M-step)，故EM算法便是在這兩個步驟之間交替而成，直到達成特定的停止條件後便回傳數值結果 (EM-EM-EM-... until met stopping criterion)。

2.3 高斯混合的EM偽代碼

所以我們現在就可以來給出計算步驟的偽代碼了，當然最初的高斯混合分佈的所有參數都是隨機的，透過EM步驟不斷交替逐漸收斂到最佳的解上：

1. 初始化參數群$\boldsymbol{\theta}=\{\mu_k,\Sigma_k,\pi_k\}$
2. E-step: 利用Eq. (10)計算每個資料點對第$k$個高斯分佈的置信度$\gamma(z_{nk})$
3. M-step: 利用$\gamma(z_{nk})$和Eqs. (12-15)重新計算最佳參數群$\boldsymbol{\theta}$
4. 檢查M-step所計算得到的參數是否收斂或達成停止條件，若有則停止計算，若否則重回第2步

明顯地EM演算法在計算是單純是演算Eqs. (10,12-15)這些數學公式罷了，並沒有用到複雜的邏輯判別或是樣本生成，僅僅是讓電腦幫我們做這些拉雜的計算。

三、超越高斯分佈: 通用EM

我們這節簡要的描述一個通用EM算法的概念，且不再要求資料的生成必須來自高斯分佈。假定如圖7的圖模型，我們知道資料集$\mathbf{x}$和隱變量$\mathbf{z}_k$之間存在一定的連結 (兩個其實都是隨機變量)，而它們之間又條件相依於參數群$\boldsymbol{\theta}$，一個顯然的聯合分佈變可以寫成$p(\mathbf{x},\mathbf{z}_k|\boldsymbol{\theta})$。原則上如果$\{\mathbf{x},\mathbf{z}_k\}$構成一個完整資料集 (complete data set) 那麼直接計算這個聯合分佈就得到$\mathbf{x}$的邊緣機率 $$p(\mathbf{x}|\boldsymbol{\theta}) =\sum_{k=1}^K p(\mathbf{x},\mathbf{z}_k|\boldsymbol{\theta}).\tag{16}$$ 然而隱變量$\mathbf{z}_k$不能直接觀測到，我們只獲得了$\mathbf{x}$的資訊，這稱做不完整資料集 (incomplete data set)，所以直接計算Eq. (16)是不可能的，因為我們缺了資料輸入$\mathbf{z}_k$。為了能夠計算Eq. (16)我們必須先得到關於$\mathbf{z}$的間接資訊，而這僅能從後驗機率$p(\mathbf{z}_k|\mathbf{x},\boldsymbol{\theta})$中獲得 (如上一節的Eq. (10))。
　若現在我們有一組參數群的結果叫$\boldsymbol{\theta}^{\rm old}$，它可能來自初始化的猜測，或是上一步迭代的結果，我們用它和$\mathbf{x}$一起計算後驗機率$p(\mathbf{z}_k|\mathbf{x},\boldsymbol{\theta}^{\rm old})$，這便是某個$\mathbf{z}_k$出現的機率，接著在將它乘上聯合機率分佈得到了 $$\sum_{k=1}^{K}\underbrace{p(\mathbf{z}_{k}|\mathbf{x},\boldsymbol{\theta}^{{\rm old}})}_{上一步得到的\boldsymbol{\theta}^{{\rm old}}} p(\mathbf{x},\mathbf{z}_{k}|\boldsymbol{\theta})\equiv\mathcal{Q}(\boldsymbol{\theta},\boldsymbol{\theta}^{{\rm old}}).\tag{17}$$ 而我們必須找到一個新的參數$\boldsymbol{\theta}$使得$\mathcal{Q}$有最大，即滿足 $$\boldsymbol{\theta}^{\rm new} = {\rm argmax}~\mathcal{Q}(\boldsymbol{\theta},\boldsymbol{\theta}^{{\rm old}}),\tag{18}$$ 其中稱$\mathcal{Q}$作輔助函數 (auxiliary function)。而事實上Eq. (17)並不是很rigorously的從數學中推導出來的，它是透過一個紆迴的過程，期望找到參數$\boldsymbol{\theta}$能讓data likelihood function有最大 $$\hat{\boldsymbol{\theta}}={\rm argmax}~L(\boldsymbol{\theta}) = {\rm argmax}~\prod_{n=1}^N \left[ \sum_{k=1}^K p(x_n,z_n^k|\boldsymbol{\theta})\right]\tag{19},$$ 這其實也是之前提過的MAP方法，我們假定Eq. (19)和Eq. (18)有正相關，能使Eq. (18)最大的參數也能使真正的likelihood function Eq. (19)有最大。

3.1 通用EM偽代碼

所以這裡給出通用EM的pseudocode，當給定觀測量$\mathbf{x}$和隱變數$\mathbf{z}$聯合分佈的形式 $p(\mathbf{x},\mathbf{z}|\boldsymbol{\theta})$ (就忽略下標$k$了)，我們就是要找出使$p(\mathbf{x}|\boldsymbol{\theta})$最大的參數$\boldsymbol{\theta}$： 

1. 初始化參數$\boldsymbol{\theta}^{\rm old}$
2. E-step: 計算後驗機率 (或稱置信度) $p(\mathbf{z}_{k}|\mathbf{x},\boldsymbol{\theta}^{{\rm old}})$
3. M-step: 通過Eqs. (17,18)計算新的$\boldsymbol{\theta}^{{\rm new}}$
4. 檢查是否收斂或達成停止條件，若是則停止計算，若否則用$\boldsymbol{\theta}^{{\rm new}}$取代$\boldsymbol{\theta}^{\rm old}$並重回步驟2

好了，那我們就完成所有關於EM演算法的理論了，由於這次的內容有點長，我們將在這裡暫停，在下一個part裡會直接實作EM算法，用它來修正K-means所犯的分群錯誤！