大家在學(xué)習(xí)算法的時候會學(xué)習(xí)到關(guān)于Kmeans的算法，但是網(wǎng)絡(luò)和很多機器學(xué)習(xí)算法書中關(guān)于Kmeans的算法理論核心一樣，但是代碼實現(xiàn)過于復(fù)雜，效率不高，不方便閱讀。這篇文章首先列舉出Kmeans核心的算法過程，并且會給出如何最大限度的在不用for循環(huán)的前提下，利用numpy, pandas的高效的功能來完成Kmeans算法。這里會用到列表解析，它是相當(dāng)于速度更快的for循環(huán)，標(biāo)題里指出的無for loop指的是除了列表解析解析以外不用for循環(huán)，來完成Kmeans算法。

一般在python數(shù)據(jù)清洗中，數(shù)據(jù)量大的情況下，for循環(huán)的方法會使的數(shù)據(jù)處理的過程特別慢，效率特別低。一個很好的解決方法就是使用numpy，pandas自帶的高級功能，不僅可以使得代碼效率大大提高，還可以使得代碼方便理解閱讀。這里在介紹用numpy，pandas來進行Kmeans算法的同時，也是帶大家復(fù)習(xí)一遍numpy，pandas用法。

1 Kmeans的算法原理

創(chuàng)建k個點作為初始質(zhì)?心（通常是隨機選擇）

當(dāng)任意一個點的簇分配結(jié)果發(fā)生改變時:

對數(shù)據(jù)集中的每個點:

對每個質(zhì)?:

計算質(zhì)?與數(shù)據(jù)點之間的距離

將數(shù)據(jù)點分配到據(jù)其最近的簇

對每個簇，計算簇中所有點的均值并將均值作為新的質(zhì)?點

直到簇不再發(fā)?變化或者達到最大迭代次數(shù)

2 聚類損失函數(shù)

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個簇的質(zhì)心Ci指的是第i個簇,x是i個簇中的點，ci是第i個簇的質(zhì)心

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維的數(shù)據(jù)，可以用來嘗試完成Kmeans的代碼。

實際過程中，Kmeans需要能運行在多維的數(shù)據(jù)上，所以下面的代碼部分，會考慮多維的數(shù)據(jù)集，而不是僅僅2維的數(shù)據(jù)。

3 隨機生成數(shù)據(jù)點

這里的嚴(yán)格意義上不是隨機的生成k個質(zhì)心點，而是取出每個特征的最大值最小值，在最大值和最小值中取出一個隨機數(shù)作為質(zhì)心點的一個維度

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個隨機的質(zhì)心點
    #np.random.uniform(a, b, c) 隨機生成在[a,b)區(qū)間里的3個數(shù)
    #對每個特征都做此操作
    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的廣播的功能，可以直接將所有的實例和其中一個質(zhì)心點相減。如下圖，下圖中是給出相加的例子，而我們的例子是減法。

所以對于一個DataFrame來說，比如說這里的只包含x和y的data，假設(shè)我們的質(zhì)心是c = [1,1]，可以用以下的方式來給出所有的實例點的x和y和點(1,1)之間的差值。注意，這里的c可以是list，也可以是numpy array，甚至可以是元組。

$$

$$

算出每個實例的每個特征和質(zhì)心點的差距之后，則需要將所有的數(shù)平方一下，然后按每一行加起來則給出了每一個實例點到質(zhì)心的距離了

$$

$$

用的方法就是使用np.power(data - c, 2).sum(axis = 1)

def cal_distant(dataset, centers):
    #選出不是label的那些特征列
    data = dataset.loc[:, dataset.columns != 'label']
    
    #使用列表解析式的格式，對centers表里的每一行也就是每一個隨機的質(zhì)心點，都算一遍所有的點到該質(zhì)心點的距離，并且存入一個list中
    d_to_centers = [np.power(data - centers.loc[i], 2).sum(axis = 1)
                    for i in centers.index]
    
    #所有的實例點到質(zhì)心點的距離都已經(jīng)存在了list中，則可以直接帶入pd.concat里面將數(shù)據(jù)拼起來
    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 找出最近的質(zhì)心點

當(dāng)每個實例點都和中心點計算好距離后，對于每個實例點找出最近的那個中心點,可以用np.where的方法，但是pandas已經(jīng)提供更加方便的方法，用idxmin和idxmax,這2個函數(shù)可以直接給出DataFrame每行或者每列的最小值和最大值的索引，設(shè)置axis = 1則是想找出對每個實例點來說，哪個質(zhì)心點離得最近。

curr_group = d_to_centers.idxmin(axis=1)

這個時候，每個點都有了新的group，這里我們則需要開始更新我們的3個中心點了。對每一個臨時的簇來說，算出X的平均， 和Y的平均，就是這個臨時的簇的中心點。

6 重新計算新的質(zhì)心點

centers = data.loc[:, data.columns != 'label'].groupby(curr_group).mean()
centers

xy00.5484680.5234741-1.0036801.0449552-0.125490-0.475373

7 迭代

這樣我們新的質(zhì)心點就得到了，只是這個時候的算法還是沒有收斂的，需要將上面的步驟重復(fù)多次。

Kmeans代碼迭代部分就完成了，將上面的步驟做成一個函數(shù)，做成函數(shù)后，方便展示Kmeans的中間過程。

def iterate(dataset, centers):
    #計算所有的實例點到所有的質(zhì)心點之間的距離
    d_to_centers = cal_distant(dataset, centers)
    
    #得出每個實例點新的類別
    curr_group = d_to_centers.idxmin(axis=1)
    
    #算出當(dāng)前新的類別下每個簇的組內(nèi)誤差
    SSE = d_to_centers.min(axis = 1).sum()
    
    #給出在新的實例點類別下，新的質(zhì)心點的位置
    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差值不變的時候，算法停止

#創(chuàng)建一個空的SSE_list,用來存SSE的，第一個位置的數(shù)為0，無意義，只是方便收斂時最后一個SSE和上一個SSE的對比
SSE_list = [0]

#初始化質(zhì)心點
centers = initial_centers(data, k = 3)

#開始迭代
while True:
    #每次迭代中得出新的組，組內(nèi)誤差，和新的質(zhì)心點，當(dāng)前的新的質(zhì)心點會被用于下一次迭代
    curr_group, SSE, centers = iterate(data,centers)
    
    #檢查這一次算出的SSE和上一次迭代的SSE是否相同，如果相同，則收斂結(jié)束
    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

上面的函數(shù)已經(jīng)完成，當(dāng)然這里推薦大家盡量寫成class的形式更好，這里為了方便觀看，則用簡單的函數(shù)完成。

最后的函數(shù)是Kmeans_regular函數(shù)，這個函數(shù)里面包含了上面所有的函數(shù)。現(xiàn)在需要測試Kmeans_regular代碼對于多特征的數(shù)據(jù)集鳶尾花數(shù)據(jù)集，是否也能進行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

最后可以看出我們的代碼是可以適用于多特征變量的數(shù)據(jù)集，并且對于鳶尾花數(shù)據(jù)集來說，對角線上的數(shù)是預(yù)測正確的個數(shù)，準(zhǔn)確率大約為90%。

9 Kmeans中間過程以及可視化展現(xiàn)

在完成代碼后，還是需要討論一下，為什么我們的代碼的算法是那樣的，這個算法雖然看起來很有邏輯，但是它到底是從哪里來的。

這個時候，我們就需要從Kmeans的損失函數(shù)出發(fā)來解釋剛才提出的問題。對于無監(jiān)督學(xué)習(xí)算法來說，也是有一個損失函數(shù)。而我們的Kmeans的中間過程的邏輯，就是從最小化Kmeans的損失函數(shù)的過程。

假設(shè)我們有一個數(shù)據(jù)集{x_1, x_2, ..., x_N}x1,x2,...,xN， 每個樣本實例點x有多個特征。我們的目標(biāo)是將這個數(shù)據(jù)集通過某種方式切分成K份，或者說我們最后想將每個樣本點標(biāo)上一個類別(簇)，且總共有K個類別，使得每個樣本點到各自的簇中心點的距離最小，并且u_kuk來表示各個簇的中心點。

我們還需要一些其他的符號，比如說r_{nk}rnk, 它的值是0或者1。下標(biāo)k代表的是第k個簇，下標(biāo)n表示的是第n個樣本點。

舉例說明，加入當(dāng)前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

這個也可以把想象成獨熱編碼。

將上面的符號解釋完了后，以下就是損失函數(shù)。這里是使用了求和嵌套了求和的公式，并且也引入了剛才所提到了r_{nk}rnk。這個損失函數(shù)其實很好理解，在給定的k個中心點u_kuk以及分配好了各個實例點屬于哪一個簇之后，計算各個實例點到各自的簇中心點的距離，距離平方以下并且相加起來，就是損失函數(shù)。這個公式其實也就是在算簇內(nèi)誤差和。

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

那怎么來最小化這個損失函數(shù)呢，用的就是EM算法，這個算法總體來說分2個步驟，Expectation和Maximization，對Kmeans來說M應(yīng)該說是Minimization

Expection:

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

Minimization:

保持當(dāng)前實例點的簇的類別不變，為了整體降低損失函數(shù)，可以對每個簇內(nèi)的損失函數(shù)公式做微分。由于當(dāng)前我們的各個點的類別是不變的，變的是u_{k}uk，所以做的微分是基于u_{k}uk的

\fracgeybsqlxm7mc{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其實就是在算各個簇內(nèi)的新的中心點，也就是對各個簇內(nèi)所有的實例點的各個特征做平均數(shù)。

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

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

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

假如說現(xiàn)在有4個簇的話，Kmeans算法不一定能使最后的SSE最小。對于2列的數(shù)據(jù)集來說，我們?nèi)?組隨機的質(zhì)心點來做對比。

第一組為設(shè)置seed為5的時候，以下為演示的結(jié)果。

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

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

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

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