大家在学习算法的时候会学习到关于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函数，这个函数里面包含了上面所有的函数。现在需要测试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算法。