热线电话:13121318867

登录
首页精彩阅读Kmeans算法精简版(无for loop循环)
Kmeans算法精简版(无for loop循环)
2020-05-27
收藏

大家在学习算法的时候会学习到关于Kmeans的算法,但是网络和很多机器学习算法书中关于Kmeans的算法理论核心一样,但是代码实现过于复杂,效率不高,不方便阅读。这篇文章首先列举出Kmeans核心的算法过程,并且会给出如何最大限度的在不用for循环的前提下,利用numpy, pandas的高效的功能来完成Kmeans算法。这里会用到列表解析,它是相当于速度更快的for循环,标题里指出的无for loop指的是除了列表解析解析以外不用for循环,来完成Kmeans算法

一般在pythonshujuqingxi/' style='color:#000;font-size:inherit;'>python数据清洗中,数据量大的情况下,for循环的方法会使的数据处理的过程特别慢,效率特别低。一个很好的解决方法就是使用numpypandas自带的高级功能,不仅可以使得代码效率大大提高,还可以使得代码方便理解阅读。这里在介绍用numpypandas来进行Kmeans算法的同时,也是带大家复习一遍numpypandas用法。

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


数据分析咨询请扫描二维码

最新资讯
更多
客服在线
立即咨询