IA9b - K-means分群與EM算法: 實作

這一篇文章延續了上一篇IA9a的內容，但側重在python的實作，所以所有的節數、圖案以及公式編號都會依照IA9a的編號接序下來，所以重複的內容就直接引用，不在這裡重新贅述。

四、實作EM

4.1 高斯混合分佈

讓我們從1.1節的協方差矩陣開始，生成的圖形就如圖4 (腳本內加了seed random方便我們每次追蹤)，另外也看到了執行5次K-means算法後 (再多次也不會變了) 得到的分群結果如圖5，這顯然與我們的預期有極大的落差！
　但我們將按照3.3結給出EM算法的偽代碼，假設資料群是高斯混合分佈的形式，重新執行一遍分群，我們將看到EM算法能很好的分出我們預期的兩群的結果。我們預計會有兩群，所以$k=1,2$，且第一步需要初始化參數組$\boldsymbol{\theta}=\{\mu_k,\Sigma_k,\pi_k\}$。初始化的過程是給定一個起始點，原則上隨機選擇的方式是合理的，不過為了計算上方便，可以讓初始的協方差矩陣都為單位矩陣，$\pi_k$則為每個群均分 (這裡就是每個都是0.5)，如下代碼

# Initializing parameters
sigma1,sigma2 = np.eye(2),np.eye(2)
pi1 = 0.5
pi2 = 1-pi1
N = len(cluster)

完成初始的參數選擇後，我們用偽代碼的第二步，也就是E-step計算置信度$\gamma (z_{nk})$。這裡我將Eq. (10)列出來方便解說 $$\gamma (z_{nk})=\frac{\pi_k \mathcal{N}(\mathbf{x}|\mu_k,\Sigma_k)}{\sum_{k=1}^2 \pi_k \mathcal{N}(\mathbf{x}|\mu_k,\Sigma_k)}\tag{10}.$$ 因為我們是假設兩群都是高斯分佈，所以對第一個群而言的分子就是

# Starting E-step: Calculating responsibilities
# Numerator for k=1
k1 = pi1*multivariate_normal.pdf(cluster,p1,sigma1)

直接透過python腳本我們很容易理解了分子的意義便是計算每個資料點分別對應到$\{\mu_1,\Sigma_1\}$的高斯分佈的機率密度，所以愈靠近$k=1$群的點，它分子的值就愈大。接著計算對$k=2$以及分母並完整的寫下Eq. (10)的結果

k2 = pi2*multivariate_normal.pdf(cluster,p2,sigma2)
gamma1 = k1/(k1+k2)
gamma2 = k2/(k1+k2)

不難理解Eq. (10)的分母就是將兩個分子結果相加即可，最後的gamma1和gamma2分別表示$\gamma(z_{n1})$和$\gamma(z_{n2})$。所以我們可以看下圖7，它是透過初始值所計算出的置信度。顏色代表它們有大於0.5的置信度是屬於同一群。

圖7. 初始條件下所計算的置信度，實線圓圈代表confidence level，內圈是$1\sigma$而外圈是$2\sigma$

　既然我們已經分出兩群了，我們變透過這兩群的點來計算M-step所需要的值，基本上在M-step裡就是要透過置信度來更新參數組$\boldsymbol{\theta}$，計算所需的公式是Eqs. (12-15)，實現的代碼如下

# Starting M-step: Re-estimate parameters that maximizes expectation
# Re-evaluating N_k, Eq. (13) 
N1,N2 = np.sum(gamma1),np.sum(gamma2)

# Re-evaluating mean mu_k, Eq. (12)
p1 = np.sum(((cluster.T)*gamma1).T,axis=0)/N1
p2 = np.sum(((cluster.T)*gamma2).T,axis=0)/N2

# Re-evaluating pi_k, Eq. (15)
pi1 = N1/N
pi2 = N2/N

上面的部份更新了$\{\mu_k,\pi_k\}$，因這兩個計算比較簡單。而協方差矩陣的計算按照Eq. (14)是 $$\Sigma_k = \frac{1}{N_k}\sum_{n=1}^N \gamma(z_{nk})(x_n-\mu_k)\otimes (x_n-\mu_k).\tag{14}$$ 假使我要計算新的$\Sigma_1$，那麼它必須逐一計算$\mathbf{x}$內每個資料點對$\mu_1$的差外積後乘上該點對應的$\gamma (z_{n1})$，當計算完$\mathbf{x}$內所有的點後，我們會得到$N$個$2\times 2$的矩陣 (因為我們考慮的是2維的資料)，將這$N$個矩陣加總並除以$N_1$即得$\Sigma_1$，對$\Sigma_2$的情形一樣，代碼如下

# Re-evaluating covariance matrix
sigma1,sigma2 = [],[]
for i in range(N):
    sigma1.append(gamma1[i]*np.outer((cluster[i]-p1),(cluster[i]-p1)))
    sigma2.append(gamma2[i]*np.outer((cluster[i]-p2),(cluster[i]-p2)))
sigma1 = np.sum(sigma1,axis=0)/N1
sigma2 = np.sum(sigma2,axis=0)/N2

OK, very well! 從代碼註解# Starting E-step...開始一路到上面結束我們變完成一次迭代啦，結果如下圖8。

圖8. EM算法經過一次迭代後的結果，為了排版乾淨我們只秀出兩個群對應到$1\sigma$所包含的區域

　不難發現我們的$\mu_k$已經透過E-step的置信度分群之後更新位置了，但距離正確的分群還有一段距離，所以我們可以加上for迴圈來重複迭帶個幾次，如下

for iters in range(max_iters):
    ...

其中...就是上面E/M-steps交互計算的代碼部份，最大迭代次數由max_iters決定，當然在本文提供的腳本裏面，你會注意到我還提供了一個停止條件，若當$\mu_k$和上一次迭代的距離變化小於0.01時，不管有沒有達到max_iters的次數都會結束計算，這裡不贅述。而這份代碼經過23次迭帶後會碰到停止條件，我們把部份的結果羅列如下圖9。

圖9. 部份迭代後的結果，實線圓圈的區域表$1\sigma$ CL

　在圖9裡秀出了每一次迭代過程中$\mu_k$走過的軌跡，以及該次迭代後修正的$\Sigma_k$對應到$1\sigma$ CL區間，不像K-means (圖5和6) 最終錯誤的分群，EM很正確的分出來我們預期的結果。

4.2 Beyond Gaussians

在上一節中我們用高斯混合為例，發現了EM算法比傳統的K-means更正確的分出了結果，而事實上數學可以證明K-means就是將EM中置信度全取為1的特例，這種我們叫作hard K-means，而EM又叫作soft K-means (McKay 2003)，不過這只是歷史名詞，我們不必太拘泥。
　事實上，分佈$p(x|\theta)$不一定要是高斯，同理$p(z|x)$也不一定要是Eq. (6)的形式，在這種情況下我們就必須假定一個分佈$p(x,z)$並採用第三節的方案，而代碼也不能直接照抄4.1節的結果 (It's Gaussian mixtures ad hoc!)，腳本是死的，方法是活的，了解後便可以理出正確的思路，高斯以外的討論可以見Bishop 2008和Murphy 2012。

五、小結

在IA9裡我們看出EM算法將機率分佈的概念引入分群的方法中，使得它較hard K-means更有彈性，雖然它必須花費較多的計算效能，但卻能得到更接近我們預期的結果，並且這個算法在computer vision裡用來分割圖像或分離前後景也有著重要的應用 (Prince 2012)。但即使這裡引入了機率 (主要是每個資料點對應的pdf) 的計算，我們仍然離真正的隨機性算法八字沒一撇，EM仍是所謂的決定性算法，與IA8用Gibbs採樣去噪本質仍是不一樣的，我會在之後的文章繼續討論這個議題。啊對了，雖然在IA9a的引文裡提到要證明如何從高斯混合EM到K-means，不過就和費馬 (Pierre de Fermat) 說的一樣：空間不夠，而目前我是時間不夠😜，所以只能留待之後再補遺啦！


Source on github: IA9_EM.ipynb