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算法。