【機器學習】算法原理詳細推導與實現(六):k-means算法

之前幾個章節都是介紹有監督學習，這個章節介紹無監督學習，這是一個被稱為k-means的聚類算法，也叫做k均值聚類算法。

聚類算法

在講監督學習的時候，通常會畫這樣一張圖：

這時候需要用logistic回歸或者SVM將這些數據分成正負兩類，這個過程稱之為監督學習，是因為對於每一個訓練樣本都給出了正確的類標簽。

在無監督學習中，經常會研究一些不同的問題。假如給定若干個點組成的數據集合：

所有的點都沒有像監督學習那樣給出類標簽和所謂的學習樣本，這時候需要依靠算法本身來發現數據中的結構。在上面的這張圖中，可以很明顯的發現這些數據被分成了兩簇，所以一個無監督學習算法將會是聚類算法，算法會將這樣的數據聚集成幾個不同的類。

聚類算法很多應用場景，舉幾個最常用的：

1. 在生物學應用中，經常需要對不同的東西進行聚類，假設有很多基因的數據，你希望對它們進行聚類以便更好的理解不同種類的基因對應的生物功能
2. 在市場調查中，假設你有一個數據庫，里面保存了不同顧客的行為，如果對這些數據進行聚類，可以將市場分為幾個不同的部分從而可以對不同的部分指定相應的銷售策略
3. 在圖片的應用中，可以將一幅照片分成若干個一致的像素子集，去嘗試理解照片的內容
4. 等等...

聚類的基本思想是：給定一組數據集合，聚集成若干個屬性一致的類。

k-means聚類

這個算法被稱之為k-means聚類算法，用於尋找數據集合中的類，算法的輸入是一個無標記的數據集合\({x^{(1)},x^{(2)},...,x^{(m)}}\)，因為這是無監督學習算法，所以在集合中只能看到\(x\)，沒有類標記\(y\)。k-means聚類算法是將樣本聚類成\(k\)個簇(cluster)，具體算法步驟如下：

step 1 隨機選取k個聚類質心點(cluster centroids)，那么就等於存在了\(k\)個簇\(c^{(k)}\)：

\[\mu_1,\mu_2,...,\mu_k \in R^n \]

step 2 對於每一個\(x^{(i)}\)，需要計算與每個質心\(\mu_j\)的距離，\(x^{(i)}\)則屬於與他距離最近質心\(\mu_j\)的簇\(c^{(j)}\)：

\[c^{(j)}:=argmin_j||x^{(i)}-\mu_j||^2,j \in 1,2,...,k \]

step 3 對於每一個類\(c^{(j)}\)，重新計算該簇質心的值：

\[\mu_j:=\frac{\sum^m_{i=1}l\{c^{(i)}=j\}x^{(i)}}{\sum^m_{i=1}l\{c^{(i)}=j\}} \]

之后需要重復step 2和step 3直到算法收斂，下面圖中對上述步驟進行解釋，存在數據點如下所示：

假設我們取\(k=2\)，那么會在數據集合中隨機選取兩個點作為質心，即下圖中紅色的點\(\mu_1\)和藍色的點\(\mu_2\)：

分別計算每一個\(x^{(i)}\)和質心\(\mu_1\)、\(\mu_2\)的距離，\(x^{(i)}\)離哪個\(\mu_j\)更近，那么\(x^{(i)}\)就屬於哪個\(c^{(j)}\)，即哪些點離紅色\(\mu_1\)近則屬於\(c^{(1)}\)，離藍色近則屬於\(c^{(2)}\)。第一次將\(x^{(i)}\)分類后效果如下：

下一步是更新簇\(c^{(j)}\)的質心，計算所有紅色點的平均值，得到新的質心\(\mu_{1\_new}\)；計算所有藍色點的平均值，得到新的質心\(\mu_{2\_new}\)，如下圖所示：

再次重復計算每一個\(x^{(i)}\)和質心的距離，更新質心的值。多次迭代收斂后，即使進行更多次的迭代，\(x^{(i)}\)的類別和質心的值都不會再改變了：

這里涉及到一個問題，如何保證k-means是收斂的？前面算法強調的是結束條件是收斂，為了保證算法完全收斂，這里引入畸變函數(distortion function)：

\[J(c,\mu)=\sum^m_{i=1}||x^{(i)}-\mu_{c^{(i)}}||^2 \]

\(J(c,\mu)\)表示每個樣本點\(x^{(i)}\)到其質心距離的平方和，當\(J(c,\mu)\)沒有達到最小值，可以固定\(c^{(j)}\)更新每個簇的質心\(\mu_j\)，質心變化后固定質心的值\(\mu_j\)重新划分簇\(c^{(j)}\)也變化，不斷迭代。當\(J(c,\mu)\)達到最小值時，\(\mu_j\)和\(c^{(j)}\)也同時收斂。(這個過程和前面【機器學習】算法原理詳細推導與實現(五):支持向量機(下)中的SMO優化算法算法的過程很相似，都是固定一組值或者一個值，更新另外一組或者一個值，使其函數優化到極值，這個過程叫做坐標上升，這里不作推導)。實際上可能會有多組\(\mu\)和\(c\)能夠使得\(J(c,\mu)\)取得最小值，但是這種情況並不多見。

由於畸變函數\(J(c,\mu)\)是非凸函數，所以意味着不能保證取的最小值是全局最小值，也就是說k-means對隨機取的質心的初始位置比較敏感。一般達到局部最優已經滿足分類的需求了，如果比較介意的話，可以多跑幾次k-means算法，然后取\(J(c,\mu)\)最小值的\(\mu\)和\(c\)。

k值確定

很多人不知道怎么確定數據集需要分多少個類(簇)，因為數據是無監督學習算法，k值需要認為的去設定。所以這里會提供兩種方法去確定k值。

第一種方法：

觀察法，本文的例子可以看出，把數據集畫在圖中顯示，就能很明顯的看到應該划分2個類(簇)：

並非所有的數據都像上面的數據一樣，一眼可以看出來分2個類(簇)，所以介紹一般比較常用的第二種方法。

第二種方法：

輪廓系數（Silhouette Coefficient），是聚類效果好壞的一種評價方式。最早由 Peter J. Rousseeuw 在 1986 提出。它結合內聚度和分離度兩種因素。可以用來在相同原始數據的基礎上用來評價不同算法、或者算法不同運行方式對聚類結果所產生的影響。

按照上面計算k-means算法的步驟計算完后，假設\(x^{(i)}\)屬於簇\(c^{(i)}\)，那么需要計算如下兩個值：

1. \(a(i)\)，計算\(x^{(i)}\)與同一簇中其他點的平均距離，代表\(x^{(i)}\)與同一簇中其他點不相似的程度
2. \(b(i)\)，計算\(x^{(i)}\)與另一個與\(x^{(i)}\)最近的簇\(c\)，與簇\(c\)內所有點平均距離，代表\(x^{(i)}\)與最相鄰簇\(c\)的不相似程度

最終輪廓系數的計算公式為：

\[S(i)=\frac{b(i)-a(i)}{max(a(i),b(i))} \]

輪廓系數的范圍為\([-1, 1]\)，越趨近於1代表內聚度和分離度都相對較優。

所以可以在k-means算法開始的時候，先設置k值的范圍\(k \in [2, n]\)，從而計算k取每一個值的輪廓系數，輪廓系數最小的那個k值就是最優的分類總數。

實例

假設存在數據集為如下樣式：

import matplotlib.pyplot as plt
import numpy as np
from sklearn.cluster import KMeans
from sklearn import metrics

# figsize繪圖的寬度和高度，也就是像素
plt.figure(figsize=(8, 10))
x1 = np.array([1, 2, 3, 1, 5, 6, 5, 5, 6, 7, 8, 9, 7, 9])
x2 = np.array([1, 3, 2, 2, 8, 6, 7, 6, 7, 1, 2, 1, 1, 3])
X = np.array(list(zip(x1, x2))).reshape(len(x1), 2)
# print(X)

# x,y軸的繪圖范圍
plt.xlim([0, 10])
plt.ylim([0, 10])
plt.title('sample')
plt.scatter(x1, x2)
# 點的顏色
colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'b']
# 點的形狀
markers = ['o', 's', 'D', 'v', '^', 'p', '*', '+']
plt.show()

雖然觀察法可以知道這個數據集合只要設置\(k=3\)就好了，但是這里還是想用輪廓系數來搜索最佳的k值。假設不知道的情況下，這里取\(k \in [2, 3, 4, 5, 8]\)：

# 測試的k值
tests = [2, 3, 4, 5, 8]
subplot_counter = 1
for t in tests:
    subplot_counter += 1
    plt.subplot(3, 2, subplot_counter)
    kmeans_model = KMeans(n_clusters=t).fit(X)
    for i, l in enumerate(kmeans_model.labels_):
        plt.plot(x1[i], x2[i], color=colors[l], marker=markers[l], ls='None')
        # 每個點對應的標簽值
        # print(kmeans_model.labels_)
        plt.xlim([0, 10])
        plt.ylim([0, 10])
        plt.title('K = %s, Silhouette Coefficient = %.03f' % (t, metrics.silhouette_score(X, kmeans_model.labels_, metric='euclidean')))

得到的結果為：

可以看到當\(k=3\)時輪廓系數最大為0.722

數據和代碼下載請關注公眾號【 機器學習和大數據挖掘 】，后台回復【 機器學習 】即可獲取