2018/12/6

k-means clustering(k-均值叢集化)

10 k-means clustering (k-均值叢集化)

k-means clustering 是指將給定的資料分成k個聚群 (group),每一群又稱為一個叢集 (cluster),此種聚類是一種unsupervised learning,也就是沒有訓練的過程,直接以資料進行處理。演算法決定如何以給定的度量將資料集聚類於 k 個叢集,也就是找出 k 個重心 (centroids) 以這些重心為叢集的中心將資料收納,所以是一種聚類的演算法,一般 k 是一個系統的超參數,可由 elbow method (手肘方法)決定最佳值。而度量一般可採用歐氏距離或距離平方,根據距離於多維的空間將資料點指配於最近的重心。此聚類常用於如社群網路將會員指配到適當的群組或是推薦系統用於產品行銷。

下式可用為 k-means 演算法的 cost function,在一個固定的k值下計算
$J=\sum_{i=1}^{K}\sum \nolimits_{i\in C_k}\left \| x_i-\mu _k \right \|^2$
其中$\mu_k$是叢集 k 的重心,此式計算每個叢集裡的成員到重心的距離平方。

∎ k-means algorithm
  • 事先指定 k值
  • 隨機選 k 個點為重心,對每一個點算出到這 k 個重心的矩離,並依距離將所有點分配到最近的重心,以每一重心為中心形成叢集,共形成 k 個叢集.
  • 再以每個叢集的平均重心為新的重心重新計算每個點到每個新重心的距離,算過後再把點重新指配給最近的重心,形成新的 k 個叢集
  • 重覆上述動作,直至沒有任何點改變重心歸屬。
∎ k-means 實作練習 (from scratch)
以文獻裡的一組資料,用 python 程式練習 k-means 演算法,資料如下
import numpy as np
import matplotlib.pyplot as plt
x0=np.array([[7,5,7,3,4,1,0,2,8,6,5,3]])
x1=np.array([[5,7,7,3,6,4,0,2,7,8,5,7]])
plt.plot(x0,x1, 'r.')
plt.xlabel('x0')
plt.ylabel('x1')
plt.show()
XXX = np.vstack((x0, x1)).T
練習的資料
把圖中的點分成兩個叢集,一開始隨機抓出兩點當重心,如抓出第5和第11點 亦即$(4,6)$和$(5,5)$,以下列程式計算
C1 = np.array([4,6])
C2 = np.array([5,5])
dist1 = [0,0,0,0,0,0,0,0,0,0,0,0]
dist2 = [0,0,0,0,0,0,0,0,0,0,0,0]
lastC = [0,0,0,0,0,0,0,0,0,0,0,0]
newC = [0,0,0,0,0,0,0,0,0,0,0,0]
change = [0,0,0,0,0,0,0,0,0,0,0,0]
for i in range(12):
    dist = np.linalg.norm(XXX[i] - C1)
    dist1[i] = dist
    dist = np.linalg.norm(XXX[i] - C2)
    dist2[i] = dist
    if dist1[i] < dist2[i]:
        newC[i] = 1
    else:
        newC[i] = 2
    if newC[i] != lastC[i]:
        change[i] = 1
count1 = 0
count2 = 0
C1_new = np.array([0,0])
C2_new = np.array([0,0])
for i in range(12):
    if newC[i] == 1:
        count1 += 1;
        C1_new = C1_new + XXX[i]
    else:
        count2 += 1;
        C2_new = C2_new + XXX[i]
C1_new = C1_new / count1
C2_new = C2_new / count2
for i in range(12):
    print('To C1 %.4f, C2 %.4f, newC %d, change: %d' % 
          (dist1[i], dist2[i], newC[i],change[i]))
print(C1_new,C2_new)
To C1 3.1623, C2 2.0000, newC 2, change: 1
To C1 1.4142, C2 2.0000, newC 1, change: 1
To C1 3.1623, C2 2.8284, newC 2, change: 1
To C1 3.1623, C2 2.8284, newC 2, change: 1
To C1 0.0000, C2 1.4142, newC 1, change: 1
To C1 3.6056, C2 4.1231, newC 1, change: 1
To C1 7.2111, C2 7.0711, newC 2, change: 1
To C1 4.4721, C2 4.2426, newC 2, change: 1
To C1 4.1231, C2 3.6056, newC 2, change: 1
To C1 2.8284, C2 3.1623, newC 1, change: 1
To C1 1.4142, C2 0.0000, newC 2, change: 1
To C1 1.4142, C2 2.8284, newC 1, change: 1
[3.8 6.4] [4.57142857 4.14285714]
輸出列印的是點到C1 和C2 的距離,還有新指配的重心,因一開始都還沒指配,所以執行後指狀況都是True。結果顯示第一輪計算後有5個指配到C1,其他的指配到C2,而最一列 是在兩叢集裡的點的新重心。將新重心值帶入程式的第1,2行,也把新重心帶入第6行,再執行一次計算,得
To C1 3.4928, C2 2.5754, newC 2, change: 0
To C1 1.3416, C2 2.8891, newC 1, change: 0
To C1 3.2558, C2 3.7498, newC 1, change: 1
To C1 3.4928, C2 1.9431, newC 2, change: 0
To C1 0.4472, C2 1.9431, newC 1, change: 0
To C1 3.6878, C2 3.5743, newC 2, change: 1
To C1 7.4431, C2 6.1694, newC 2, change: 0
To C1 4.7539, C2 3.3472, newC 2, change: 0
To C1 4.2426, C2 4.4630, newC 1, change: 1
To C1 2.7203, C2 4.1132, newC 1, change: 0
To C1 1.8439, C2 0.9583, newC 2, change: 0
To C1 1.0000, C2 3.2608, newC 1, change: 0
[5.5 7. ] [3.         3.16666667]
可見有三個點歸屬被改變,而新叢集的位置也改變了。再依上述方法分別再把新重心及新指配值代入程式中,再執行第三輪計算,得
To C1 2.5000, C2 4.4001, newC 1, change: 1
To C1 0.5000, C2 4.3237, newC 1, change: 0
To C1 1.5000, C2 5.5403, newC 1, change: 0
To C1 4.7170, C2 0.1667, newC 2, change: 0
To C1 1.8028, C2 3.0046, newC 1, change: 0
To C1 5.4083, C2 2.1667, newC 2, change: 0
To C1 8.9022, C2 4.3621, newC 2, change: 0
To C1 6.1033, C2 1.5366, newC 2, change: 0
To C1 2.5000, C2 6.3004, newC 1, change: 0
To C1 1.1180, C2 5.6887, newC 1, change: 0
To C1 2.0616, C2 2.7131, newC 1, change: 1
To C1 2.5000, C2 3.8333, newC 1, change: 0
[5.625 6.5  ] [1.5  2.25]
可見依還有兩個點歸屬被改變,再帶入新的重心及指配,再執行第四輪,見
To C1 2.0349, C2 6.1492, newC 1, change: 0
To C1 0.8004, C2 5.9002, newC 1, change: 0
To C1 1.4631, C2 7.2672, newC 1, change: 0
To C1 4.3750, C2 1.6771, newC 2, change: 0
To C1 1.7002, C2 4.5069, newC 1, change: 0
To C1 5.2574, C2 1.8200, newC 2, change: 0
To C1 8.5960, C2 2.7042, newC 2, change: 0
To C1 5.7785, C2 0.5590, newC 2, change: 0
To C1 2.4271, C2 8.0506, newC 1, change: 0
To C1 1.5462, C2 7.3015, newC 1, change: 0
To C1 1.6250, C2 4.4511, newC 1, change: 0
To C1 2.6722, C2 4.9812, newC 1, change: 0
[5.625 6.5  ] [1.5  2.25]
這輪的計算中沒有點的歸屬被改變,所以 k-means 聚類已完成了,最後的重心列在最後一行,兩個重心是$(5.625, 6.5)$ 和 $(1.5, 2.25)$。接著就可依此兩個叢集求 cost function,如此完成$k=2$ 叢集化。後續可用相同的方法求$k=3, \cdots n$ 的 cost function,再決定最佳的 $k$ 值。

上列程式中我們採用的 cost function 是距離而不是距離平方,但因是比較大小,所以有沒有平方都一樣。因 sklearn 的函數有屬性可求出距離的平方,所以我們再以下列程式計算分配完成後所有的點到所屬重心的距離平方的總和,這屬性在sklearn 叫 inertia (慣量),
import math
tot = 0
for i in range(12):
    if newC[i] == 1:
        tot = tot + math.pow(dist1[i],2)
    else:
        tot = tot + math.pow(dist2[i],2)
print(tot)
得41.625

∎ k-means with sklearn
sklearn 有個 KMeans class 可叫用。此函式的 fit 方法接受以 $m \times n$  array 形式儲存的資料,  $m$ 是實例的數目, $n$ 是 features 的數目。 如有12個實例,每實例有兩個 features,則輸入的矩陣形式是 $12\times 2$。 KMeans 的constructer 裡有很多參數可以 宣告,重要的如下:
  • n_clusters: 指定 k 值
  • max_iter: 對於單次初始值計算的最大迭代次數
  • n_init: 重新選擇初始值的次數
  • init: 制定初始值選擇的算法,如 k-means++
  • n_jobs: CPU進程的數目 ,設為 $-1$ 指跑滿 CPU
以上參數中以 n_clusters最重要,其他的參數初學可用預設值就好了。另此 class 的幾個重要的屬性如下:
  • labels_ : 取得聚類後標籤,也就是成員歸屬 。
  • cluster_centers_: 取得聚類的中心(重心) 
  • inertia_: 取得聚類完成後所有點到中心距離平方的總和。
下列程式以sklearn 的 KMeans class 重練上面的例子
from sklearn.cluster import KMeans
estimator = KMeans(n_clusters=2)
estimator.fit(XXX)
label_pred = estimator.labels_
centroids = estimator.cluster_centers_ 
inertia = estimator.inertia_
print(label_pred)
print(centroids)
print(inertia)
[0 0 0 1 0 1 1 1 0 0 0 0]
[[5.625 6.5  ]
 [1.5   2.25 ]]
41.625
上列程式 estimator 是 KMeans 的實例,之後用fit method 帶入data進行聚類,並用 label_pred 取出聚類結果,centroids 取出兩個重心,也用 inertia 取出距離平方的總和。結果上我們沒用KMeans class 的程式所求的結果是相同的。注意 沒用 KMeans class 我們的 C1 是[5.625 6.5  ] 而 C2 是 [1.5  2.25],使用 Kmeans class 所得的是 C0 [[5.625 6.5  ] 而 C1是 [1.5   2.25 ]],且兩種求法的慣量都是 41.625,這指出兩種方法的結果是相同的。

∎ elbow method 決定最佳 k 值
以下列 的程式試圖以平的距離平方來決定系統的最佳k值,程式中$k=1\sim  12$
from scipy.spatial.distance import cdist
K = range(1,13)
meanDispersions = []
for k in K:
    kmeans = KMeans(n_clusters=k)
    kmeans.fit(XXX)
    meanDispersions.append(sum(np.min(cdist(XXX,
        kmeans.cluster_centers_, 'euclidean'), axis=1)) / XXX.shape[0])
plt.plot(K, meanDispersions, 'bx-')
plt.xlabel('k')
plt.ylabel('Average Dispersion')
plt.title('Selecting k with the Elbow Method')
plt.show()
不同 k 值
此結果顯示很難斷定一個較佳的 k值。k-mean 演算法的兩大問題中一個是若初始的重心選得不好,最後結果可能收斂到一個區域的最佳而不是整體最佳,另一個問題就是有時很難找到一個較佳的k值。像上圖,平均距離由最(只有一個重心)到最小(每個點都是重心,當然是0)遞減,則並沒有很明顯的較佳值,最主要的原因是這組資料點分散的很平均,本身並沒有形成簇狀,所以也就很難說應聚成幾個群。 再看下列例子
import matplotlib.pyplot as plt
c1x = np.random.uniform(0.5, 1.5, (1, 10))
c1y = np.random.uniform(0.5, 1.5, (1, 10))
c2x = np.random.uniform(3.5, 4.5, (1, 10))
c2y = np.random.uniform(3.5, 4.5, (1, 10))
x = np.hstack((c1x, c2x))
y = np.hstack((c1y, c2y))
X = np.vstack((x, y)).T
K = range(1, 10)
meanDispersions = []
for k in K:
    kmeans = KMeans(n_clusters=k)
    kmeans.fit(X)
    meanDispersions.append(sum(np.min(cdist(X,
        kmeans.cluster_centers_, 'euclidean'), axis=1)) / X.shape[0])
plt.plot(K, meanDispersions, 'bx-')
plt.xlabel('k')
plt.ylabel('Average Dispersion')
plt.title('elbow method')
plt.show()
elbow method
此例中很明顯$k=2$應該是最好的,因這組資料是以 $(1,1)$ 和 $(4,4)$ 為均勻分佈平均值所產生的,所以本身具有兩個簇狀。

參考文獻 Gavin Hackeling, Mastering Machine Learning with scikit-learn, 2nd, Packt Publishing, 2017

沒有留言:

張貼留言