本文內(nèi)容

閱讀須知：

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

本文將介紹：

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

K-means簡介

聚類是一個將數(shù)據(jù)集中在某些方面相似的數(shù)據(jù)成員進(jìn)行分類組織的過程，聚類就是一種發(fā)現(xiàn)這種內(nèi)在結(jié)構(gòu)的技術(shù)，聚類技術(shù)經(jīng)常被稱為無監(jiān)督學(xué)習(xí)。

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

K-means原理

設(shè)有樣本$(x_1,x_2,...,x_n)$，其中每個樣本有d個特征（由d維實向量組成），K-means聚類的目標(biāo)是將$n$個樣本聚類至$k(\le n)$個集合 $S=\{S_1,S_2,...,S_k\}$ 中，使簇中樣本的距離和達(dá)到最小。
即用公式表示為：
$$\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$中所有點的質(zhì)心，因此上式相當(dāng)于最小化簇中每兩點的平方距離：
$$\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算法根據(jù)此從距離入手，迭代找出（局部）最小時的聚類中心$S_i$。具體步驟如下：

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

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

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

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

Python實現(xiàn)

基本實現(xiàn)

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

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數(shù)組，約定shape為：(數(shù)據(jù)數(shù)量，數(shù)據(jù)維度):type data: numpy.ndarray"""self.features_count = data.shape[1]# 初始化聚類中心（維度：k個 * features種數(shù)）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)]# 對每個點與聚類中心進(jìn)行距離計算for feature_set in data:# 預(yù)測分類classification = self.predict(feature_set)# 加入類別self.classifications[classification].append(feature_set)# 記錄前一次的結(jié)果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# 如果都滿足條件（上面循環(huán)沒break），則返回else:returndef predict(self, data):# 距離distances = np.linalg.norm(data - self.centroids, axis=1)# 最小距離索引return distances.argmin()

測試代碼

blobs函數(shù)產(chǎn)生若干個數(shù)據(jù)點云。n_features是維度，centers是中心數(shù)目，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

輸出聚類圖象的函數(shù)，支持2D和3D，8類別以內(nèi)不同顏色區(qū)分。

def kmeans_plot(kmeans_model):"""又長又臭的函數(shù)，簡單可視化2d或3d kmeans聚類結(jié)果，不是算法必須的，直接使用即可。:param kmeans_model: 訓(xùn)練的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()

測試調(diào)用

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)

測試結(jié)果

計算效率問題

問題描述

該實現(xiàn)雖然可以使用，但是在進(jìn)行大規(guī)模數(shù)據(jù)處理時非常慢，通過性能監(jiān)控找出耗時最大的過程：

問題分析

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

解決辦法

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

計算優(yōu)化實現(xiàn)

numpy廣播

計算優(yōu)化借助了numpy庫中向量廣播的特性。

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

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

實現(xiàn)代碼

• 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數(shù)組，約定shape為：(數(shù)據(jù)數(shù)量，數(shù)據(jù)維度):type data: numpy.ndarray"""self.features_count = data.shape[1]# 初始化聚類中心（維度：k個 * features種數(shù)）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個數(shù)據(jù)進(jìn)行從1到20類的聚類。記錄兩種做法分別的耗時。

500100050001000020000

普通實現(xiàn)	0.86	2.91	23.72	59.17	129.69
計算優(yōu)化實現(xiàn)	0.07	0.98	3.84	6.41	28.34
比率	12.2857	2.9694	6.1771	9.2309	4.5762

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

自動K值選取

常用方法

手肘法

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

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

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

Gap statistic

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

簡單解釋

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

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

詳細(xì)解釋

• 原理性較強，非必須理解
• 本節(jié)內(nèi)容翻譯自原論文

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

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

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

使用Gap Statistic選取最優(yōu)的k

首先，該算法需要評估聚類誤差$D_k$，由于最后的標(biāo)準(zhǔn)是相對的，因此$D_k$的求法并無過多約束，只要能體現(xiàn)其誤差即可，因此還是選擇最簡便的做法：各點到聚類中心的距離為單個誤差，將其加和作為最終的誤差，由sum_distance處理這一過程。
雖然聚類中心K的最優(yōu)解是不依賴算法客觀存在的，但由于不同K-means實現(xiàn)會得出不同的$D_k$，因此為了$Gap(K)$具有可比性，需要借助同一種K-means實現(xiàn)求解誤差$D_k$，因此sum_distance中統(tǒng)一使用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函數(shù)需要給定k的測試范圍ks以及蒙特卡洛模擬次數(shù)nrefs，負(fù)責(zé)產(chǎn)生隨機樣品估計$E(log{D}_k)$、計算$log{D}_k$，然后返回范圍內(nèi)各個$k$值對應(yīng)的$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

調(diào)用代碼

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

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

結(jié)果輸出

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

[-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

結(jié)果分析

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

Gap Statistic總結(jié)

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

參考文獻(xiàn)

[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.

非文獻(xiàn)參考

【機器學(xué)習(xí)】K-means（非常詳細(xì)）
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/

總結(jié)

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

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。