Kmeans算法精簡版(無for loop循環)

​ 你們在學習算法的時候會學習到關於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 = \sum_{i=1}^k\sum_{x\in C_{i}}(c_{i} - x)^2SSE=i=1∑k​x∈Ci​∑​(ci​−x)2

C_{i}指的是第i個簇, x是i個簇中的點，c_{i}是第i個簇的質心Ci​指的是第i個簇,x是i個簇中的點，ci​是第i個簇的質心

import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt
from sklearn.datasets import make_blobs
import seaborn as sns

#r = np.random.randint(1,100)
r = 4
#print(r)
k = 3
x , 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'] = y
sim_data.head(5)

data = sim_data.copy()

plt.scatter(sim_data['x'], sim_data['y'], c = y)

上圖是一個隨機生成的2維的數據，能夠用來嘗試完成Kmeans的代碼。

實際過程當中，Kmeans須要能運行在多維的數據上，因此下面的代碼部分，會考慮多維的數據集，而不是僅僅2維的數據。

3 隨機生成數據點

​ 這裏的嚴格意義上不是隨機的生成k個質心點，而是取出每一個特徵的最大值最小值，在最大值和最小值中取出一個隨機數做爲質心點的一個維度

def initial_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

x

y

0

0.122575

0.021762

1

-0.922596

1.367504

2

-0.677202

-0.411821

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)

def cal_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)

0

1

2

0

0.153365

3.935546

0.528286

1

1.987879

0.088006

2.462444

2

0.027977

2.361753

0.795004

3

0.543410

5.183283

0.565696

4

1.505514

2.248264

4.031165

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

x

y

0

0.548468

0.523474

1

-1.003680

1.044955

2

-0.125490

-0.475373

7 迭代

這樣咱們新的質心點就獲得了，只是這個時候的算法仍是沒有收斂的，須要將上面的步驟重複屢次。

Kmeans代碼迭代部分就完成了，將上面的步驟作成一個函數，作成函數後，方便展現Kmeans的中間過程。

def iterate(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)

#開始迭代
while True:
    #每次迭代中得出新的組，組內偏差，和新的質心點，當前的新的質心點會被用於下一次迭代
    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 代碼整合

算法完成了，將全部的代碼整合在一塊兒

def initial_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.T

def cal_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)

def iterate(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

def Kmeans_regular(data, k = 3):
    SSE_list = [0]
    centers = initial_centers(data, k = k)

    while True:
        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_iris
data_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)

col_0

0

1

2

label

0

50

0

0

1

0

48

2

2

0

14

36

np.diag(pd.crosstab(iris['label'], curr_group)).sum() /  iris.shape[0]

0.8933333333333333

最後能夠看出咱們的代碼是能夠適用於多特徵變量的數據集，而且對於鳶尾花數據集來講，對角線上的數是預測正確的個數，準確率大約爲90%。

9 Kmeans中間過程以及可視化展示

​ 在完成代碼後，仍是須要討論一下，爲何咱們的代碼的算法是那樣的，這個算法雖然看起來頗有邏輯，可是它究竟是從哪裏來的。

​ 這個時候，咱們就須要從Kmeans的損失函數出發來解釋剛纔提出的問題。對於無監督學習算法來講，也是有一個損失函數。而咱們的Kmeans的中間過程的邏輯，就是從最小化Kmeans的損失函數的過程。

​ 假設咱們有一個數據集{x_1, x_2, ..., x_N}x1​,x2​,...,xN​， 每一個樣本實例點x有多個特徵。咱們的目標是將這個數據集經過某種方式切分紅K份，或者說咱們最後想將每一個樣本點標上一個類別(簇)，且總共有K個類別，使得每一個樣本點到各自的簇中心點的距離最小，而且u_kuk​來表示各個簇的中心點。

咱們還須要一些其餘的符號，好比說r_{nk}rnk​, 它的值是0或者1。下標k表明的是第k個簇，下標n表示的是第n個樣本點。

舉例說明，加入當前K=3，k的可取1，2，3。對於第一個實例點n = 1來講它屬於第3個簇，因此

r_{n=1, k = 1} = 0rn=1,k=1​=0

r_{n=1, k = 2} = 0rn=1,k=2​=0

r_{n=1, k = 3} = 1rn=1,k=3​=1

這個也能夠把想象成獨熱編碼。

將上面的符號解釋完了後，如下就是損失函數。這裏是使用了求和嵌套了求和的公式，而且也引入了剛纔所提到了r_{nk}rnk​。這個損失函數其實很好理解，在給定的k箇中心點u_kuk​以及分配好了各個實例點屬於哪個簇以後，計算各個實例點到各自的簇中心點的距離，距離平方如下而且相加起來，就是損失函數。這個公式其實也就是在算簇內偏差和。

C = \sum_{n=1}^N\sum_{k=1}^K r_{nk} (x_n - u_k)^2C=n=1∑N​k=1∑K​rnk​(xn​−uk​)2

那怎麼來最小化這個損失函數呢，用的就是EM算法，這個算法整體來講分2個步驟，Expectation和Maximization，對Kmeans來講M應該說是Minimization

Expection:

保持u_{k}uk​不變，也就是各個簇的中心點的位置不變，計算各個實例點到哪一個u_{k}uk​最近，將各個實例點劃分到離各自最近的那個簇裏面去，從而使得總體SSE下降。

Minimization:

保持當前實例點的簇的類別不變，爲了總體下降損失函數，能夠對每一個簇內的損失函數公式作微分。因爲當前咱們的各個點的類別是不變的，變的是u_{k}uk​，因此作的微分是基於u_{k}uk​的

\frac{d}{du_{k}}\sum_{k=1}^K r_{nk} (x_n - u_k)^2 = 0duk​d​k=1∑K​rnk​(xn​−uk​)2=0

-2\sum_{k=1}^K r_{nk} (x_n - u_k) = 0−2k=1∑K​rnk​(xn​−uk​)=0

u_{k} = \frac{\sum_{n} r_{nk} x_{n}}{\sum_{n} r_{nk}}uk​=∑n​rnk​∑n​rnk​xn​​

得出來的u_{k}uk​其實就是在算各個簇內的新的中心點，也就是對各個簇內全部的實例點的各個特徵作平均數。

這時候獲得新的中心點u_{k}uk​, 緊接着再到E階段，保持u_{k}uk​更新簇類別，再到M階段，保持簇類別不變動新u_{k}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算法。