大家在學習算法的時候會學習到關於Kmeans的算法，但是網絡和很多機器學習算法書中關於Kmeans的算法理論核心一樣，但是代碼實現過於複雜，效率不高，不方便閱讀。這篇文章首先列舉出Kmeans核心的算法過程，並且會給出如何最大限度的在不用for循環的前提下，利用numpy, pandas的高效的功能來完成Kmeans算法。這裡會用到列表解析，它是相當於速度更快的for循環，標題裡指出的無for loop指的是除了列表解析解析以外不用for循環，來完成Kmeans算法。
一般在python數據清洗中，數據量大的情況下，for循環的方法會使的數據處理的過程特別慢，效率特別低。一個很好的解決方法就是使用numpy，pandas自帶的高級功能，不僅可以使得代碼效率大大提高，還可以使得代碼方便理解閱讀。這裡在介紹用numpy，pandas來進行Kmeans算法的同時，也是帶大家複習一遍numpy，pandas用法。
1 Kmeans的算法原理
創建k個點作為初始質心（通常是隨機選擇）
當任意一個點的簇分配結果發生改變時:
對數據集中的每個點:
對每個質:
計算質與數據點之間的距離
將數據點分配到據其最近的簇
對每個簇，計算簇中所有點的均值並將均值作為新的質點
直到簇不再發變化或者達到最大迭代次數
2 聚類損失函數
SSE=i=1∑kx∈Ci∑(cix)2
Ci指的是第i個簇,x是i個簇中的點，ci是第i個簇的質心
import numpy as npimport pandas as pdimport matplotlib as mplimport matplotlib.pyplot as pltfrom sklearn.datasets import make_blobsimport seaborn as sns
#r = np.random.randint(1,100)r =4#print(r)k =3x , y = make_blobs(n_samples =51, cluster_std =[0.3,0.3,0.3], centers =[[0,0],[1,1],[-1,1]],random_state = r )sim_data = pd.DataFrame(x, columns =['x','y'])sim_data['label']= ysim_data.head(5)data = sim_data.copy()plt.scatter(sim_data['x'], sim_data['y'], c = y)
上圖是一個隨機生成的2維的數據，可以用來嘗試完成Kmeans的代碼。
實際過程中，Kmeans需要能運行在多維的數據上，所以下面的代碼部分，會考慮多維的數據集，而不是僅僅2維的數據。
3 隨機生成數據點
這裡的嚴格意義上不是隨機的生成k個質心點，而是取出每個特徵的最大值最小值，在最大值和最小值中取出一個隨機數作為質心點的一個維度
definitial_centers(datasets, k =3):#首先將datasets的特徵名取出來，這裡需要除去label那一列 cols = datasets.columns data_content = datasets.loc[:, cols !='label']#直接用describe的方法將每一列的最小值最大值取出來 range_info = data_content.describe().loc[['min','max']]#用列表解析的方法和np.random.uniform的方法生成k個隨機的質心點#np.random.uniform(a, b, c) 隨機生成在[a,b)區間裡的3個數#對每個特徵都做此操作 k_randoms =[np.random.uniform(range_info[i]['min'], range_info[i]['max'], k)for i in range_info.columns] centers = pd.DataFrame(k_randoms, index = range_info.columns)return centers.T
centers = initial_centers(data, k =3)centers
4 計算所有的點到所有中心點的距離
將每一個中心點取出來，然後使用pandas的廣播的功能，可以直接將所有的實例和其中一個質心點相減。如下圖，下圖中是給出相加的例子，而我們的例子是減法。
所以對於一個DataFrame來說，比如說這裡的只包含x和y的data，假設我們的質心是c = [1,1]，可以用以下的方式來給出所有的實例點的x和y和點(1,1)之間的差值。注意，這裡的c可以是list，也可以是numpy array，甚至可以是元組。
$$
$$
算出每個實例的每個特徵和質心點的差距之後，則需要將所有的數平方一下，然後按每一行加起來則給出了每一個實例點到質心的距離了
$$
$$
用的方法就是使用np.power(data - c, 2).sum(axis = 1)
defcal_distant(dataset, centers):#選出不是label的那些特徵列 data = dataset.loc[:, dataset.columns !='label']#使用列表解析式的格式，對centers表裡的每一行也就是每一個隨機的質心點，都算一遍所有的點到該質心點的距離，並且存入一個list中 d_to_centers =[np.power(data - centers.loc[i],2).sum(axis =1)for i in centers.index]#所有的實例點到質心點的距離都已經存在了list中，則可以直接帶入pd.concat裡面將數據拼起來return pd.concat(d_to_centers, axis =1)
d_to_centers = cal_distant(data, centers)d_to_centers.head(5)
5 找出最近的質心點
當每個實例點都和中心點計算好距離後，對於每個實例點找出最近的那個中心點,可以用np.where的方法，但是pandas已經提供更加方便的方法，用idxmin和idxmax,這2個函數可以直接給出DataFrame每行或者每列的最小值和最大值的索引，設置axis = 1則是想找出對每個實例點來說，哪個質心點離得最近。
curr_group = d_to_centers.idxmin(axis=1)
這個時候，每個點都有了新的group，這裡我們則需要開始更新我們的3個中心點了。對每一個臨時的簇來說，算出X的平均， 和Y的平均，就是這個臨時的簇的中心點。
6 重新計算新的質心點
centers = data.loc[:, data.columns !='label'].groupby(curr_group).mean()centers
7 迭代
這樣我們新的質心點就得到了，只是這個時候的算法還是沒有收斂的，需要將上面的步驟重複多次。
Kmeans代碼迭代部分就完成了，將上面的步驟做成一個函數，做成函數後，方便展示Kmeans的中間過程。
defiterate(dataset, centers):#計算所有的實例點到所有的質心點之間的距離 d_to_centers = cal_distant(dataset, centers)#得出每個實例點新的類別 curr_group = d_to_centers.idxmin(axis=1)#算出當前新的類別下每個簇的組內誤差 SSE = d_to_centers.min(axis =1).sum()#給出在新的實例點類別下，新的質心點的位置 centers = dataset.loc[:, dataset.columns !='label'].groupby(curr_group).mean()return curr_group, SSE, centers
curr_group, SSE, centers = iterate(data,centers)
centers, SSE
( x y 0 0.892579 0.931085 1 -1.003680 1.044955 2 0.008740 -0.130172, 19.041432436034352)
最後需要判斷什麼時候迭代停止，可以判斷SSE差值不變的時候，算法停止
#創建一個空的SSE_list,用來存SSE的，第一個位置的數為0，無意義，只是方便收斂時最後一個SSE和上一個SSE的對比SSE_list =[0]#初始化質心點centers = initial_centers(data, k =3)#開始迭代whileTrue:#每次迭代中得出新的組，組內誤差，和新的質心點，當前的新的質心點會被用於下一次迭代 curr_group, SSE, centers = iterate(data,centers)#檢查這一次算出的SSE和上一次迭代的SSE是否相同，如果相同，則收斂結束if SSE_list[-1]== SSE:break#如果不相同，則記錄SSE，進入下一次迭代 SSE_list.append(SSE)
SSE_list
[0, 37.86874675507244, 11.231524142566894, 8.419267088238051]
8 代碼整合
算法完成了，將所有的代碼整合在一起
definitial_centers(datasets, k =3): cols = datasets.columns data_content = datasets.loc[:, cols !='label'] range_info = data_content.describe().loc[['min','max']] k_randoms =[np.random.uniform(range_info[i]['min'], range_info[i]['max'], k)for i in range_info.columns] centers = pd.DataFrame(k_randoms, index = range_info.columns)return centers.Tdefcal_distant(dataset, centers): data = dataset.loc[:, dataset.columns !='label'] d_to_centers =[np.power(data - centers.loc[i],2).sum(axis =1)for i in centers.index]return pd.concat(d_to_centers, axis =1)defiterate(dataset, centers): d_to_centers = cal_distant(dataset, centers) curr_group = d_to_centers.idxmin(axis=1) SSE = d_to_centers.min(axis =1).sum() centers = dataset.loc[:, dataset.columns !='label'].groupby(curr_group).mean()return curr_group, SSE, centersdefKmeans_regular(data, k =3): SSE_list =[0] centers = initial_centers(data, k = k)whileTrue: curr_group, SSE, centers = iterate(data,centers)if SSE_list[-1]== SSE:break SSE_list.append(SSE)return curr_group, SSE_list, centers
上面的函數已經完成，當然這裡推薦大家儘量寫成class的形式更好，這裡為了方便觀看，則用簡單的函數完成。
最後的函數是Kmeans_regular函數，這個函數裡面包含了上面所有的函數。現在需要測試Kmeans_regular代碼對於多特徵的數據集鳶尾花數據集，是否也能進行Kmeans聚類算法
from sklearn.datasets import load_irisdata_dict = load_iris()iris = pd.DataFrame(data_dict.data, columns = data_dict.feature_names)iris['label']= data_dict.target
curr_group, SSE_list, centers = Kmeans_regular(iris.copy(), k =3)
np.array(SSE_list)
array([ 0. , 589.73485975, 115.8301874 , 83.29216169, 79.45325846, 78.91005674, 78.85144143])
pd.crosstab(iris['label'], curr_group)
np.diag(pd.crosstab(iris['label'], curr_group)).sum()/ iris.shape[0]
0.8933333333333333
最後可以看出我們的代碼是可以適用於多特徵變量的數據集，並且對於鳶尾花數據集來說，對角線上的數是預測正確的個數，準確率大約為90%。
9 Kmeans中間過程以及可視化展現
在完成代碼後，還是需要討論一下，為什麼我們的代碼的算法是那樣的，這個算法雖然看起來很有邏輯，但是它到底是從哪裡來的。
這個時候，我們就需要從Kmeans的損失函數出發來解釋剛才提出的問題。對於無監督學習算法來說，也是有一個損失函數。而我們的Kmeans的中間過程的邏輯，就是從最小化Kmeans的損失函數的過程。
假設我們有一個數據集x1,x2,...,xN， 每個樣本實例點x有多個特徵。我們的目標是將這個數據集通過某種方式切分成K份，或者說我們最後想將每個樣本點標上一個類別(簇)，且總共有K個類別，使得每個樣本點到各自的簇中心點的距離最小，並且uk來表示各個簇的中心點。
我們還需要一些其他的符號，比如說rnk, 它的值是0或者1。下標k代表的是第k個簇，下標n表示的是第n個樣本點。
舉例說明，加入當前K=3，k的可取1，2，3。對於第一個實例點n = 1來說它屬於第3個簇，所以
rn=1,k=1=0
rn=1,k=2=0
rn=1,k=3=1
這個也可以把想像成獨熱編碼。
將上面的符號解釋完了後，以下就是損失函數。這裡是使用了求和嵌套了求和的公式，並且也引入了剛才所提到了rnk。這個損失函數其實很好理解，在給定的k個中心點uk以及分配好了各個實例點屬於哪一個簇之後，計算各個實例點到各自的簇中心點的距離，距離平方以下並且相加起來，就是損失函數。這個公式其實也就是在算簇內誤差和。
C=n=1∑Nk=1∑Krnk(xnuk)2
那怎麼來最小化這個損失函數呢，用的就是EM算法，這個算法總體來說分2個步驟，Expectation和Maximization，對Kmeans來說M應該說是Minimization
Expection:
保持uk不變，也就是各個簇的中心點的位置不變，計算各個實例點到哪個uk最近，將各個實例點劃分到離各自最近的那個簇裡面去，從而使得整體SSE降低。
Minimization:
保持當前實例點的簇的類別不變，為了整體降低損失函數，可以對每個簇內的損失函數公式做微分。由於當前我們的各個點的類別是不變的，變的是uk，所以做的微分是基於uk的
dukdk=1∑Krnk(xnuk)2=0
2k=1∑Krnk(xnuk)=0
uk=∑nrnk∑nrnkxn
得出來的uk其實就是在算各個簇內的新的中心點，也就是對各個簇內所有的實例點的各個特徵做平均數。
這時候得到新的中心點uk, 緊接著再到E階段，保持uk更新簇類別，再到M階段，保持簇類別不變更新uk，不斷的迭代知道SSE不變位置。這個就是Kmeans的算法過程。下面將用plotly可視化，動態展示Kmeans的過程。
使用之前寫好的函數，然後將數據的中間過程通過plotly展示出來。因為代碼比較長，所以這裡就不展示代碼了。由於當前是一個markdown，這裡放入一個gif圖片用來展示最後的Kmeans中間過程。
[外鏈圖片轉存失敗,源站可能有防盜鏈機制,建議將圖片保存下來直接上傳(img-zaiSFuGB-1589536000909)(…/…/…/…/…/…/0 AI-work/B 部門/SEO/202001/rowdata/倪向陽_SEO_2020_01/Kmeans_Plotly中間過程/Kmeans_1.gif)]
對於這個數據集來看的話，我們的Kmeans算法可以使得每一個點最終可以找到各自的簇，但是這個算法也是有缺陷的，比如以下例子。
假如說現在有4個簇的話，Kmeans算法不一定能使最後的SSE最小。對於2列的數據集來說，我們取2組隨機的質心點來做對比。
第一組為設置seed為5的時候，以下為演示的結果。
[外鏈圖片轉存失敗,源站可能有防盜鏈機制,建議將圖片保存下來直接上傳(img-t7uRDf5U-1589536000912)(…/…/…/…/…/…/0 AI-work/B 部門/SEO/202001/rowdata/倪向陽_SEO_2020_01/Kmeans_Plotly中間過程/Kmeans_2.gif)]
從上面的動圖可以看出一共用了8次迭代，才收斂。那加入我們的seed為1的話，隨機的質心點的分布會變的很離譜，會導致下面的結果。這裡我們加快動畫的速度。
[外鏈圖片轉存失敗,源站可能有防盜鏈機制,建議將圖片保存下來直接上傳(img-mEk3AXPf-1589536000913)(…/…/…/…/…/…/0 AI-work/B 部門/SEO/202001/rowdata/倪向陽_SEO_2020_01/Kmeans_Plotly中間過程/Kmeans_3.gif)]
這裡用34次，數據才迭代收斂，並且可以看出，在迭代的過程中，差點陷入了一個局部最小的一個情況。所以對於複雜的數據來說的話，我們最後看到迭代的次數會明顯的增加。
假如說我們的數據集再變的集中一點，其中的2個簇，稍微近一點，我們會看到以下的結果。
[外鏈圖片轉存失敗,源站可能有防盜鏈機制,建議將圖片保存下來直接上傳(img-ztYH1wDQ-1589536000914)(…/…/…/…/…/…/0 AI-work/B 部門/SEO/202001/rowdata/倪向陽_SEO_2020_01/Kmeans_Plotly中間過程/Kmeans_4.gif)]
所以在這次迭代的過程中，我們明顯看到其中有個質心點消失了，原因就是因為由於點的分布的原因和初始質心點的原因，最開始隨機生成的一個離所有的點都最遠的質心點，由於它離所有的點都最遠，所以導致了在迭代的過程中，沒有任何一個點屬於這個質心點，最後導致這個點消失了。所以這個就是Kmeans算法的缺陷，那怎麼來優化這個算法了，我們可以利用BiKmeans算法。