大家在學習算法的時候會學習到關于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∑kx∈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

xy00.1225750.0217621-0.9225961.3675042-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)

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

xy00.5484680.5234741-1.0036801.0449552-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函數，這個函數里面包含了上面所有的函數?，F在需要測試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_0012label0500010482201436

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∑Nk=1∑Krnk(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 = 0dukdk=1∑Krnk(xn?uk)2=0

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

u_{k} = \frac{\sum_{n} r_{nk} x_{n}}{\sum_{n} r_{nk}}uk=∑nrnk∑nrnkxn

得出來的u_{k}uk其實就是在算各個簇內的新的中心點，也就是對各個簇內所有的實例點的各個特征做平均數。

這時候得到新的中心點u_{k}uk, 緊接著再到E階段，保持u_{k}uk更新簇類別，再到M階段，保持簇類別不變更新u_{k}uk，不斷的迭代知道SSE不變位置。這個就是Kmeans的算法過程。下面將用plotly可視化，動態展示Kmeans的過程。

使用之前寫好的函數，然后將數據的中間過程通過plotly展示出來。因為代碼比較長，所以這里就不展示代碼了。由于當前是一個markdown，這里放入一個gif圖片用來展示最后的Kmeans中間過程。

對于這個數據集來看的話，我們的Kmeans算法可以使得每一個點最終可以找到各自的簇，但是這個算法也是有缺陷的，比如以下例子。

假如說現在有4個簇的話，Kmeans算法不一定能使最后的SSE最小。對于2列的數據集來說，我們取2組隨機的質心點來做對比。

第一組為設置seed為5的時候，以下為演示的結果。

從上面的動圖可以看出一共用了8次迭代，才收斂。那加入我們的seed為1的話，隨機的質心點的分布會變的很離譜，會導致下面的結果。這里我們加快動畫的速度。

這里用34次，數據才迭代收斂，并且可以看出，在迭代的過程中，差點陷入了一個局部最小的一個情況。所以對于復雜的數據來說的話，我們最后看到迭代的次數會明顯的增加。

假如說我們的數據集再變的集中一點，其中的2個簇，稍微近一點，我們會看到以下的結果。

所以在這次迭代的過程中，我們明顯看到其中有個質心點消失了，原因就是因為由于點的分布的原因和初始質心點的原因，最開始隨機生成的一個離所有的點都最遠的質心點，由于它離所有的點都最遠，所以導致了在迭代的過程中，沒有任何一個點屬于這個質心點，最后導致這個點消失了。所以這個就是Kmeans算法的缺陷，那怎么來優化這個算法了，我們可以利用BiKmeans算法。