本文內容

閱讀須知：

• 閱讀本文需要有一定的Python及Numpy基礎

本文將介紹：

• K-means算法實現步驟
• 使用Python實現K-means算法
• 借助Numpy的向量計算提升計算速度
• 使用Gap Statistic法自動選取合適的聚類中心數K

K-means簡介

聚類是一個將數據集中在某些方面相似的數據成員進行分類組織的過程，聚類就是一種發現這種內在結構的技術，聚類技術經常被稱為無監督學習。

k均值聚類是最著名的劃分聚類算法，由于簡潔和效率使得他成為所有聚類算法中最廣泛使用的。給定一個數據點集合和需要的聚類數目k，k由用戶指定，k均值算法根據某個距離函數反復把數據分入k個聚類中。

K-means原理

設有樣本$(x_1,x_2,...,x_n)$，其中每個樣本有d個特征（由d維實向量組成），K-means聚類的目標是將$n$個樣本聚類至$k(\le n)$個集合 $S=\{S_1,S_2,...,S_k\}$ 中，使簇中樣本的距離和達到最小。
即用公式表示為：
$$\underset{\mathbf{S}}{\mathrm{argmin}}\sum _{i=1}^k\sum _{\mathbf{x}\in S_i}{\Vert \mathbf{x}?{\mathbf{\mu }}_i\Vert }^2=\underset{\mathbf{S}}{\mathrm{argmin}}\sum _{i=1}^k|S_i|\mathrm{Var}{S}_i$$
其中${\mu }_i$ 是$S_i$中所有點的質心，因此上式相當于最小化簇中每兩點的平方距離：
$$\underset{\mathbf{S}}{\mathrm{argmin}}\sum _{i=1}^k\frac{1}{|S_i|}\sum _{\mathbf{x},\mathbf{y}\in S_i}{\Vert \mathbf{x}?\mathbf{y}\Vert }^2$$

K-means步驟

K-means算法根據此從距離入手，迭代找出（局部）最小時的聚類中心$S_i$。具體步驟如下：

選擇初始化的 $k$ 個樣本作為初始聚類中心 $a=\{a_1,a_2,...,a_k\}$；

針對數據集中每個樣本$x_i$計算它到 $k$ 個聚類中心的距離并將其分到距離最小的聚類中心所對應的類中；

針對每個類別$S_j$，重新計算它的聚類中心 $a_j=\frac{1}{|S_i|}{\sum }_{x\in S_i}x$；

重復上面2、3兩步操作，直到達到某個中止條件（迭代次數、最小誤差變化等）。

Python實現

基本實現

• 需要用聚類中心數量k進行初始化，此外可以給定迭代次數與最小誤差等終止條件。
• fit函數以若干n維特征的數據作為輸入。執行后，通過classifications獲取分類結果；centroids獲取聚類中心。
• Predict函數以一個n維特征的數據作為輸入，輸出歸屬的聚類中心索引。
• 每步操作的作用詳見代碼注釋

class Kmeans:def __init__(self, k=2, tolerance=0.01, max_iter=300):self.k = kself.tol = toleranceself.max_iter = max_iterself.features_count = -1self.classifications = Noneself.centroids = Nonedef fit(self, data):""":param data: numpy數組，約定shape為：(數據數量，數據維度):type data: numpy.ndarray"""self.features_count = data.shape[1]# 初始化聚類中心（維度：k個 * features種數）self.centroids = np.zeros([self.k, data.shape[1]])for i in range(self.k):self.centroids[i] = data[i]for i in range(self.max_iter):# 清空聚類列表self.classifications = [[] for i in range(self.k)]# 對每個點與聚類中心進行距離計算for feature_set in data:# 預測分類classification = self.predict(feature_set)# 加入類別self.classifications[classification].append(feature_set)# 記錄前一次的結果prev_centroids = np.ndarray.copy(self.centroids)# 更新中心for classification in range(self.k):self.centroids[classification] = np.average(self.classifications[classification], axis=0)# 檢測相鄰兩次中心的變化情況for c in range(self.k):if np.linalg.norm(prev_centroids[c] - self.centroids[c]) > self.tol:break# 如果都滿足條件（上面循環沒break），則返回else:returndef predict(self, data):# 距離distances = np.linalg.norm(data - self.centroids, axis=1)# 最小距離索引return distances.argmin()

測試代碼

blobs函數產生若干個數據點云。n_features是維度，centers是中心數目，random_state是隨機種子。

from sklearn.datasets import make_blobsdef blobs(n_samples=300, n_features=2, centers=1, cluster_std=0.60, random_state=0):points, _ = make_blobs(n_samples=n_samples,n_features=n_features,centers=centers,cluster_std=cluster_std,random_state=random_state)return points

輸出聚類圖象的函數，支持2D和3D，8類別以內不同顏色區分。

def kmeans_plot(kmeans_model):"""又長又臭的函數，簡單可視化2d或3d kmeans聚類結果，不是算法必須的，直接使用即可。:param kmeans_model: 訓練的kmeans模型:type kmeans_model: Kmeans | FastKmeans"""style.use('ggplot')colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'w']# 2Dif kmeans_model.features_count == 2:fig = plt.figure(figsize=(12, 12))ax = fig.add_subplot()for i in range(kmeans_model.k):color = colors[i%len(colors)]for feature_set in kmeans_model.classifications[i]:ax.scatter(feature_set[0], feature_set[1], marker="x", color=color, s=50, linewidths=1)for centroid in kmeans_model.centroids:ax.scatter(centroid[0], centroid[1], marker="o", color="k", s=50, linewidths=3)# 3Delif kmeans_model.features_count == 3:fig = plt.figure(figsize=(12, 12))ax = fig.add_subplot(projection='3d')for i in range(kmeans_model.k):color = colors[i%len(colors)]for feature_set in kmeans_model.classifications[i]:ax.scatter(feature_set[0], feature_set[1], feature_set[2], marker="x", color=color, s=50, linewidths=1)for centroid in kmeans_model.centroids:ax .scatter(centroid[0], centroid[1], centroid[2], marker="o", color="k", s=50, linewidths=3)plt.show()

測試調用

model = Kmeans(k=3) model.fit(blobs(centers=3, random_state=1, n_features=2)) kmeans_plot(model) model.fit(blobs(centers=3, random_state=3, n_features=3)) kmeans_plot(model)

測試結果

計算效率問題

問題描述

該實現雖然可以使用，但是在進行大規模數據處理時非常慢，通過性能監控找出耗時最大的過程：

問題分析

分析后發現：距離計算的函數被調用了很多次。進而說明，主要函數（如距離）的計算單位是小規模點集，因此這些函數在python層面上調用了多次。眾所周知python的執行效率比較慢，特別是耗時排第一的norm函數，對于小規模計算的開銷是極大的。

解決辦法

最好的解決辦法是計算向量化：設法一次計算一批、甚至全部數據。這樣就可以把計算交給python底層的c語言實現；而且向量（或矩陣）的計算本身也是被優化過的。如此速度上會提升許多。

計算優化實現

numpy廣播

計算優化借助了numpy庫中向量廣播的特性。

如下圖，被減數是(N, 1, D)維向量，表示N個D維數據的數據集；減數是(1, K, D)維向量，表示K個D維的聚類中心。相減時，由于向量維度不匹配，numpy會將向量廣播為矩陣（結果為虛線表示的矩陣）再相減。而該例子中，所得結果矩陣恰好為：每個點與各個聚類中心的距離信息。

如此一來，甚至不需要顯式寫一個循環就可以完成一整輪的距離計算。這樣做把計算函數交給了受優化的矩陣運算，交給了高效的底層c實現。

實現代碼

• compute_l2_distance利用了numpy廣播

class FastKmeans:def __init__(self, k=2, tolerance=0.01, max_iter=300):self.k = kself.tol = toleranceself.max_iter = max_iterself.features_count = -1self.classifications = Noneself.centroids = Nonedef fit(self, data):""":param data: numpy數組，約定shape為：(數據數量，數據維度):type data: numpy.ndarray"""self.features_count = data.shape[1]# 初始化聚類中心（維度：k個 * features種數）self.centroids = np.zeros([self.k, data.shape[1]])for i in range(self.k):self.centroids[i] = data[i]points_centroid = np.zeros(data.shape[0])for i in range(self.max_iter):prev_centroid = np.ndarray.copy(points_centroid)points_centroid = self.get_closest_centroid(data)for c in range(self.centroids.shape[0]):# 獲取該類的點classification = data[points_centroid == c]# 更新聚類中心self.centroids[c] = classification.mean(axis=0)if FastKmeans.compute_l2_distance(prev_centroid, points_centroid) < self.tol:breakself.classifications = []for c in range(self.centroids.shape[0]):# 獲取該類的點self.classifications.append(data[points_centroid == c])def predict(self, data):# 距離distances = np.linalg.norm(data - self.centroids, axis=1)# 最小距離索引return distances.argmin()@staticmethoddef compute_l2_distance(x, centroid):# 利用numpy廣播計算dist = ((x - centroid) ** 2).sum(axis=x.ndim - 1)return distdef get_closest_centroid(self, x):# 遍歷每個中心，并計算中心與該點間的距離dist = FastKmeans.compute_l2_distance(x[:, None, :], self.centroids[None, :, :])# 取得距離最小的中心的索引closest_centroid_index = np.argmin(dist, axis=1)return closest_centroid_index

速度對比

每組對比實驗分別使用兩種做法對n個數據進行從1到20類的聚類。記錄兩種做法分別的耗時。

500100050001000020000

普通實現	0.86	2.91	23.72	59.17	129.69
計算優化實現	0.07	0.98	3.84	6.41	28.34
比率	12.2857	2.9694	6.1771	9.2309	4.5762

可見在不同的數據規模下均有一定提升。兩種做法的速度差距在之后選擇k的過程中更加明顯。

自動K值選取

常用方法

手肘法

其思想是：
隨著聚類數k的增大，樣本劃分會更加精細，每個簇的聚合程度會逐漸提高，那么誤差平方和SSE自然會逐漸變小。

當k小于真實聚類數時，由于k的增大會大幅增加每個簇的聚合程度，故SSE的下降幅度會很大；而當k到達真實聚類數時，再增加k所得到的聚合程度回報會迅速變小，所以SSE的下降幅度會驟減，然后隨著k值的繼續增大而趨于平緩，也就是說SSE和k的關系圖是一個手肘的形狀，而這個肘部對應的k值就是數據的真實聚類數。

但該做法的缺點是需要人工判斷“手肘”位置到底在哪，在類別多時難以判斷。

Gap statistic

其定義為：
$$Gap(K)=E(log{D}_k)?log{D}_k$$
其中，$E(log{D}_k)$是$log{D}_k$的期望，一般使用蒙特卡洛模擬產生。

簡單解釋

模擬的基本過程是：首先在樣本所在區域內按照均勻分布隨機地產生和原始樣本數一樣多的隨機樣本，并用K-means模型對這個隨機樣本聚類，計算結果的誤差和，得到一個$D_k$。重復多次就可以近似計算出$E(log{D}_k)$。

實際上，計算$E(log{D}_k)$只是為了給評判實際誤差$log{D}_k$提供一個標準。當選取合適的$K$值時，$D_k$應處于偏離（遠低于）平均值的極端情況，因此通過與期望誤差做差，得到最大$Gap(K)$所對應的K即是最好的選擇。

詳細解釋

• 原理性較強，非必須理解
• 本節內容翻譯自原論文

設聚類簇$C_r$中兩點間距離和為：
$$D_r=\sum _{i,i^{\prime }\in C_r}{d}_{i{i}^{\prime }}$$
若$d_{i{i}^{\prime }}$為歐幾里得距離，則可定義“每個簇內平方和均值”的總加和$W_k$為：
$$W_k=\sum _{r=1}^k\frac{1}{2n_r}{D}_r$$
其中因子$2$恰好使上式成立；數據規模$n$被約分掉了。

該方法通過與一個合適的零分布比較來標準化$\log (W_k)$，而最優聚類數則是使$\log (W_k)$最偏離于上述參考分布時的值$k$。因此定義：
$$Ga{p}_n(k)=E_n^{?}\{\log (W_k)\}?\log (W_k)$$
其中$E_n^{?}$表示樣本規模$n$的數據在參考分布中的數學期望。考慮抽樣分布后，使$Ga{p}_n(k)$最大化的值就是所估計聚類數$k$的最優值。
這個方法的操作是一般化的，可以用在以任意形式計算距離$d_{i{i}^{\prime }}$的任意的聚類方法。

Gap Statistic的思路啟發是：假設有n個、k簇、均勻分布的p維數據點，設想它們在各自的簇中均勻分布，則$\log (W_k)$的期望近似是：
$$\log (\frac{pn}{12})?(\frac{2}{p})\log (k)+Constant$$
如果數據確實有K個相互分離的聚類，對于$k\le K$，$\log (W_k)$應比預期下降率$(\frac{2}{p})\log (k)$下降的更快；當$k>K$，事實上是在加入不必要的聚類中心，此時由簡單代數可得$\log (W_k)$應下降的比其預期速率更慢。因此$Ga{p}_n(K)$會在$k=K$時取得最大值。

使用Gap Statistic選取最優的k

首先，該算法需要評估聚類誤差$D_k$，由于最后的標準是相對的，因此$D_k$的求法并無過多約束，只要能體現其誤差即可，因此還是選擇最簡便的做法：各點到聚類中心的距離為單個誤差，將其加和作為最終的誤差，由sum_distance處理這一過程。
雖然聚類中心K的最優解是不依賴算法客觀存在的，但由于不同K-means實現會得出不同的$D_k$，因此為了$Gap(K)$具有可比性，需要借助同一種K-means實現求解誤差$D_k$，因此sum_distance中統一使用FastKmeans。

import scipydef sum_distance(data, k):model = FastKmeans(k=k)model.fit(data)disp = 0for m in range(len(model.classifications)):disp += sum(np.linalg.norm(model.classifications[m] - model.centroids[m], axis=1))return disp

gap函數需要給定k的測試范圍ks以及蒙特卡洛模擬次數nrefs，負責產生隨機樣品估計$E(log{D}_k)$、計算$log{D}_k$，然后返回范圍內各個$k$值對應的$Gap(K)$值。

def gap(data, refs=None, nrefs=20, ks=range(1, 11)):shape = data.shapeif refs == None:tops = data.max(axis=0)bots = data.min(axis=0)dists = scipy.matrix(np.diag(tops - bots))rands = scipy.random.random_sample(size=(shape[0], shape[1], nrefs))for i in range(nrefs):rands[:, :, i] = rands[:, :, i] * dists + botselse:rands = refsgaps = np.zeros((len(ks),))for (i, k) in enumerate(ks):disp = sum_distance(data, k)refdisps = np.zeros((rands.shape[2],))for j in range(rands.shape[2]):refdisps[j] = sum_distance(rands[:, :, j], k)gaps[i] = np.lib.scimath.log(np.mean(refdisps)) - np.lib.scimath.log(disp)return gaps

調用代碼

先生成隨機點云，此處生成了一組6個聚類中心的3維點云。
之后調用gap函數。

my_data = blobs(centers=6, random_state=121, n_features=3) gaps = gap(my_data, nrefs=100)

結果輸出

此處順便做了之前兩種實現的效率對比，下面兩個輸出分別對應K-means和Fast K-means實現。

[-0.05127935 -0.01656365 0.30313623 0.78059812 1.56439394 1.709803511.67235943 1.63859173 1.62538263 -0.84390957] best k: 6 58.349995613098145[-0.05141505 -0.01660111 0.30252902 0.78299894 1.56696863 1.706896031.67327937 1.63563849 1.62526949 -0.8487873 ] best k: 6 4.088269472122192

結果分析

可以看到兩種實現的gap值（兩個列表輸出）有略微不同，但都得到$k=6$的結果，而之前生成的數據正是$6$個中心，說明算法成功檢測出最優的$k$值。
此外，可以看到K-means版本的gap函數運行需要58秒，而Fast K-means只需要4秒，速度差距超過十倍，之前的優化還算是效果拔群的。

Gap Statistic總結

• 理論上可以自動化找出最優的K。
• 但由于實現上需要借助特定的K-means算法，而K-means具有局部最優的特點，因此該算法找出的K并不一定是最優的。
• 再加上期望使用蒙特卡洛方法，想得到穩定、靠譜的答案需要花費更多時間進行頻率測試。
• 有時會出現多個峰值（設想3類別聚類完全可以對半分成6類），一般根據經驗主義選取第一個。
• 無論如何，用于自動化估計較優的K，Gap statistic已經足夠了。

參考文獻

[1] Tibshirani R , Hastie W T . Estimating the number of clusters in a data set via the gap statistic[J]. Journal of the Royal Statistical Society B, 2001, 63(2):411-423.

非文獻參考

【機器學習】K-means（非常詳細）
https://zhuanlan.zhihu.com/p/78798251
A Python implementation of the Gap Statistic
https://gist.github.com/michiexile/5635273
Nuts and Bolts of NumPy Optimization Part 2: Speed Up K-Means Clustering by 70x
https://blog.paperspace.com/speed-up-kmeans-numpy-vectorization-broadcasting-profiling/

總結

以上是生活随笔為你收集整理的【机器学习】K-means算法Python实现教程的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。