浅谈K-means算法和实现(基于Python)

Kmeans可视化
https://www.naftaliharris.com/blog/visualizing-k-means-clustering/

K-means原理

K-means 有一个著名的解释：牧师—村民模型：

有四个牧师去郊区布道，一开始牧师们随意选了几个布道点，并且把这几个布道点的情况公告给了郊区所有的村民，于是每个村民到离自己家最近的布道点去听课。
听课之后，大家觉得距离太远了，于是每个牧师统计了一下自己的课上所有的村民的地址，搬到了所有地址的中心地带，并且在海报上更新了自己的布道点的位置。
牧师每一次移动不可能离所有人都更近，有的人发现A牧师移动以后自己还不如去B牧师处听课更近，于是每个村民就去了离自己最近的布道点……
就这样，牧师每个礼拜更新自己的位置，村民根据自己的情况选择布道点，最终稳定了下来
1
2
3
4

Kmeans，一种无监督聚类算法，也是最为经典的基于划分的聚类方法，思想是: 对于给定的样本集，按照样本之间的距离大小，将样本集划分为K个类。让类内的点尽量紧密的连在一起，而让类间的距离尽量的大
步骤:

1. 先确定数据集类的个数k
2. 数据预处理: 剔除离群点(剔除掉那些可能由人为因素造成的，并非真实数据的样本点)、数据归一化(消除量纲的影响，转换方法为: (X - min) / (max - min) )、数据标准化(将采集到的数据分布转换成标准正态分布，转换公式: X减去均值，再除以标准差)
3. 在数据集中随机选取k个数据，作为初始质心
4. 计算数据集中每个样本到每个质心的欧式距离，把样本划分到距离最小的质心所属的类别
5. 根据聚类结果，重新计算质心，使用的是求均值的方法进行更新: ${\mu }_i=\frac{1}{m_i}{\sum }_{x_j\in C_i}{x}_j$ 其中，C_i是第i个类，x_j是C_i中的样本点，μ_i是第i个类的类中心点，m_i是第i个类中的样本点个数 (注意新的质心不一定是实际的样本点)，当本次计算的质心与上一次质心完全一样（或者收敛）时，停止迭代；否则更新质心，重复执行3、4两步操作
6. 考虑到临近性度量是欧几里得距离，使用误差的平方和(SSE指标)，作为度量聚类质量的好坏，后面再讲

例子: https://www.bilibili.com/video/BV1py4y1r7DN?share_source=copy_web
有以下6个点，将A3和A4作为初始质心

1.计算每个点到质心距离，将距离近的点归为一类。从下图可以看到A1与A3的距离是2.24，与A4的距离是3.61，所以A1被归为A3这一类

2.将蓝色每个点，和紫色每个点的x,y值分别求平均。获得新的质心

3.计算每个点到质心的新距离，将距离近的点归为一类

4.由于关联点没有变化，所以之后的计算结果不会改变，停止计算
5.蓝色类: A1, A3, A5。紫色类：A2, A4, A6。

代码实现

# 不调库实现
class my_Kmeans:
    def __init__(self, k):
        self.k = k

    # 距离度量函数: 计算任意点之间的距离
    def calc_distance(self, vec1, vec2):
        return np.sqrt(np.sum(np.power(vec1 - vec2, 2)))

    def fit(self, data):
        numSamples, dim = data.shape # numSamples指样本总数，dim指特征维度
        # 随机并且不重复地选取k个数据作为初始聚类中心点
        self.centers_idx = np.random.choice(numSamples, self.k, replace=False) # 得到的是质心点的索引
        self.centers = data[self.centers_idx].astype(np.float32) # 初始化聚类中心

        # ClusterAssment记录每个样本的信息: 第一列记录样本所属聚类中心的索引，第2列记录样本与聚类中心的距离的平方(用于SSE指标)
        ClusterAssment = np.zeros((numSamples, 2))
        ClusterChanged = True

        while ClusterChanged:
            ClusterChanged = False
            # 遍历所有点求每个点与各个聚类中心点的距离
            for i in range(numSamples):
                mindist = self.calc_distance(data[i], self.centers[0])
                label = 0
                for j in range(1, self.k):
                    distance = self.calc_distance(data[i], self.centers[j])
                    # 把第i个数据分配到距离最近的聚类中心
                    if distance < mindist:
                        mindist = distance
                        label = j

                # print("mindist = {}，label = {}".format(mindist, label))

                # 判断样本所属聚类中心是否发生了变化，若发生了变化则更新数据
                if ClusterAssment[i, 0] != label:
                    ClusterChanged = True
                    ClusterAssment[i, :] = label, mindist ** 2

            # print(ClusterAssment, '\n')

            # 对每个类别的样本重新求聚类中心，并更新聚类中心
            for j in range(self.k):
                # 找到所有属于类中心j的样本数据
                pointsInCluster = data[ClusterAssment[:, 0] == j]
                # print("{}: {}".format(j, pointsInCluster))
                self.centers[j, :] = (np.mean(pointsInCluster, axis=0).tolist()) # 更新聚类中心的位置

    def predict(self, data):
        numSamples, dim = data.shape
        ClusterAssment = np.zeros((numSamples, 2))
        for i in range(numSamples):
            mindist = self.calc_distance(data[i], self.centers[0])
            label_pred = 0
            for j in range(1, self.k):
                distance = self.calc_distance(data[i], self.centers[j])
                if distance < mindist:
                    mindist = distance
                    label_pred = j

            ClusterAssment[i, :] = label_pred, mindist ** 2

        return self.centers, ClusterAssment

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64

生成一些数据集，来检验Kmeans聚类效果

import matplotlib.pyplot as plt
import numpy as np
from sklearn.datasets import make_blobs

data, true_labels = make_blobs(centers = 5, random_state = 20, cluster_std = 1) # data是坐标集合，ture_labels是正确标签
plt.scatter(data[:, 0], data[:, 1])
plt.show()
1
2
3
4
5
6
7

# 不调库实现
model = my_Kmeans(k=5)
model.fit(data)

'''
model.predict()函数返回值
  centers: 聚类中心坐标
  ClusterAssment: 每个样本的信息: 第一列记录样本所属聚类中心的索引，第2列记录样本与聚类中心的距离的平方(用于SSE指标)
'''
centers, ClusterAssment = model.predict(data)

plt.figure(figsize=(6, 5))
plt.scatter(data[:, 0], data[:, 1], c = ClusterAssment[:, 0])
plt.scatter(centers[:, 0], centers[:, 1], marker = '*', color = 'r', s = 100) # 绘制类中心点
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

下图便是代码跑出来的结果，多运行几次会发现不调库实现的Kmeans聚类结果会各种各样，即不调库实现的Kmeans结果大部分情况下都是局部最优，并非全局最优

调库实现的话，无论运行多少次，结果都很好

# 调库实现
from sklearn.cluster import KMeans

model_1 = KMeans(n_clusters=5)
label_pred_1 = model_1.fit_predict(data)

plt.figure(figsize=(6, 5))
plt.scatter(data[:, 0], data[:, 1], c = label_pred_1)
plt.show()
1
2
3
4
5
6
7
8
9

缺点:

1. 对初始的类中心敏感，步骤2随机选取数据作为初始质心，可以想象，当质心之间距离很近的时候，将需要迭代很多次，会影响收敛速度，而且有可能最终只能得到局部最小值
2. 步骤3中计算距离，我们应该想到，距离受到量纲的影响，在使用时，需要考虑一下归一化；而且因为是应用于所有样本，所以会对异常数据比较敏感，对孤立点数据非常敏感，少量的噪声数据就能对平均值造成极大的影响，所以稳定性较差
   Q: 前面数据预处理部分不是已经将离群点剔除了吗？为什么还会受离群点的影响？
   A: 有些离群点是很容易判断出来的，但有些离群点并不好判断，比如说离群点处于数据分布的边缘，我们并不知道它是真实数据还是离群点，所以还是有离群点会被加入到整个K-means算法过程中
3. 样本只能归为一类，不适合多分类任务
4. k值对K-means影响很大，因为在执行算法之前，并不知道要将样本分为几类，k值选择不恰当，可能会导致模型难以收敛

优点:

1. 对于大数据集，算法时间复杂度为线性 O(NKT) 其中，N: 样本点个数、K: 聚类中心个数、T: 迭代次数
   因为每次迭代更新质心时，需要计算所有样本点和所有类中心的距离，所以是N·K，要迭代T次，再乘以T
2. 局部最优解通常已经可以满足问题需要

---

Kmeans++

前言

前面我们提到过K-means存在一些问题，其中一个问题就是K-means对初始的类中心非常敏感，Kmeans++就是专门对初始值的选择进行优化的

问题点: 随机选取k个数据，若数据非常靠近，会导致收敛很慢，甚至收敛到局部最小值
可见聚类的结果高度依赖质心的初始化，所以为了解决初始时如何选取k个类中心点的问题，引入了K-means++算法，它的核心思想是，距离已选中心点越远，被选为下一个中心点的概率越大
那么问题来了，那是不是选择距离所有已选中心点最远的那个点作为下一个中心点呢？并不是，如果这样做很容易选到异常的离群点，从而影响最终的聚类效果

轮盘法

在讲Kmeans之前，先讲一下轮盘赌法，基本上很多文章谈到轮盘赌算法，都会和遗传算法，蚁群算法一起解释
步骤:

1. 计算出每个群体照顾你每个个体的适应度 $f(i=1,2,...,m)$，m为群体大小
2. 计算每个个体被遗传到下一代群体中的概率 $P(x_i)=\frac{f(x_i)}{{\sum }_{j=1}^{m}f(x_j)}$
3. 计算每个个体的累计概率 $q(x_i)={\sum }_{k=1}^{i}P(x_k)$，这样，如果某个个体的适应度高，则它所对应的个体选择概率就会越大，通过累计概率转换后对应的线段就会越长
   
4. 选择某个个体策略为在区间[0 1]中随机产生一个数，看看该数字落在那个区间，很明显，对于适应度值较大的个体，对应的线段长度会长，这样随机产生的数字落在此区间的概率就大，该个体被选中的概率也大。同时，对于适应度较小的个体，线段长度会相对较短，随机数字在该区间的概率相对较小，但是也有被选中的可能，避免了适应度数值较小的个体被直接淘汰的问题

举例
初始种群适应度为 [169, 576, 64, 361]，那么每个个体被遗传到下一代群体的概率为:
$$\begin{gathered} P(x_1)=\frac{f(x_1)}{{\sum }_{j=1}^{m}f(x_j)}=\frac{169}{169+576+64+361}=0.14 \\ P(x_2)=\frac{f(x_2)}{{\sum }_{j=1}^{m}f(x_j)}=\frac{576}{169+576+64+361}=0.49 \\ P(x_3)=\frac{f(x_3)}{{\sum }_{j=1}^{m}f(x_j)}=\frac{64}{169+576+64+361}=0.06 \\ P(x_4)=\frac{f(x_4)}{{\sum }_{j=1}^{m}f(x_j)}=\frac{361}{169+576+64+361}=0.31 \end{gathered}$$
每个个体的累计概率为:
$$\begin{gathered} q(x_1)=\sum _{j=1}^{1}P(x_j)=0.14 \\ q(x_2)=\sum _{j=1}^{2}P(x_j)=0.14+0.49=0.63 \\ q(x_3)=\sum _{j=1}^{3}P(x_j)=0.14+0.49+0.06=0.69 \\ q(x_4)=\sum _{j=1}^{1}P(x_j)=0.14+0.9+0.06+0.31=1 \end{gathered}$$

随机产生一个数r=0.672452，发现是落在[q₂, q₃]之间，则x₃被选中

Kmeans++代码实现

步骤:

1. 从输入的数据点集合中随机选择一个点作为第一个聚类中心
2. 对于数据集中的每一个点x_i，计算它与当前已有聚类中心之间的最短距离(即与最近的一个聚类中心的距离)，用D(x)表示
3. 选择一个新的数据点作为新的聚类中心，选择的原则是：D(x)较大的点，被选取作为聚类中心的概率较大，方法是计算每个样本被选为下一个聚类中心的概率$\frac{D(x{)}^2}{{\sum }_{x\in X}D(x{)}^2}$，然后按照轮盘算法选择出下一个聚类
4. 重复2、3步骤直到选择出k个聚类质心
5. 利用这k个质心来作为初始化质心去运行标准的K-Means算法

代码

class my_Kmeans_plus:
    def __init__(self, k, init):
        self.k = k
        self.init = init # init='random'就是标准Kmeans，init='Kmeans++'就是Kmeans++

    # 距离度量函数: 计算任意点之间的距离
    def calc_distance(self, vec1, vec2):
        return np.sqrt(np.sum(np.power(vec1 - vec2, 2)))
    
    '''
    利用轮盘法选择下一个聚类中心
    参数: 
      P是各样本点被选为下一个聚类中心的概率集合
      r是区间[0, 1]之间随机生成的一个数，判断该数落在哪个区间，则这个区间被选中
    返回的是被选中的样本索引
    '''
    def RWS(self, P, r):
        q = 0 # 累计概率
        for i in range(len(P)):
            q += P[i] # P[i]表示第i个样本被选中的概率
            if r <= q: # 产生的随机数在该区间内
                return i

    def fit(self, data):
        numSamples, dim = data.shape # numSamples指样本总数，dim指特征维度

        if self.init == 'random':
            # 随机并且不重复地选取k个数据作为初始聚类中心点
            self.centers_idx = np.random.choice(numSamples, self.k, replace = False) # 得到的是质心点的索引
            # 默认类别是从0到k-1
            self.centers = data[self.centers_idx].astype(np.float32) # 初始化聚类中心
            # print(self.centers)

        elif self.init == 'Kmeans++':
            first = np.random.choice(numSamples) # 随机选出一个样本点作为第一个聚类的中心
            index_select = [first] # 将聚类中心索引存储到一个列表中
            D = np.zeros((numSamples, 1)) # 初始化D(x)
            # 继续选取k-1个点
            for j in range(1, self.k):
                for i in range(numSamples):
                    D[i] = float('inf')
                    for idx in index_select:
                        distance = self.calc_distance(data[i], data[idx])
                        D[i] = min(D[i], distance)

                # 计算概率
                P = np.square(D) / np.sum(np.square(D))
                r = np.random.random() # r为0到1的随机数
                choiced_index = self.RWS(P, r)
                index_select.append(choiced_index) # 利用轮盘法选择下一个聚类中心

            # print(index_select)
            self.centers = data[index_select].astype(np.float32)
            # print(self.centers)

        # ClusterAssment记录每个样本的信息: 第一列记录样本所属聚类中心的索引，第2列记录样本与聚类中心的距离
        ClusterAssment = np.zeros((numSamples, 2))
        ClusterChanged = True

        while ClusterChanged:
            ClusterChanged = False
            # 遍历所有点求每个点与各个聚类中心点的距离
            for i in range(numSamples):
                mindist = self.calc_distance(data[i], self.centers[0])
                label = 0
                for j in range(1, self.k):
                    distance = self.calc_distance(data[i], self.centers[j])
                    # 把第i个数据分配到距离最近的聚类中心
                    if distance < mindist:
                        mindist = distance
                        label = j

                # print("mindist = {}，label = {}".format(mindist, label))

                # 判断样本所属聚类中心是否发生了变化，若发生了变化则更新数据
                if ClusterAssment[i, 0] != label:
                    ClusterChanged = True
                    ClusterAssment[i, :] = label, mindist

            # print(ClusterAssment, '\n')

            # 对每个类别的样本重新求聚类中心，并更新聚类中心
            for j in range(self.k):
                # 找到所有属于类中心j的样本数据
                pointsInCluster = data[ClusterAssment[:, 0] == j]
                # print("{}: {}".format(j, pointsInCluster))
                self.centers[j, :] = (np.mean(pointsInCluster, axis=0).tolist()) # 更新聚类中心的位置

    def predict(self, data):
        numSamples, dim = data.shape
        ClusterAssment = np.zeros((numSamples, 2))
        for i in range(numSamples):
            mindist = self.calc_distance(data[i], self.centers[0])
            label_pred = 0
            for j in range(1, self.k):
                distance = self.calc_distance(data[i], self.centers[j])
                if distance < mindist:
                    mindist = distance
                    label_pred = j

            ClusterAssment[i, :] = label_pred, mindist ** 2

        return self.centers, ClusterAssment

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104

调用自己写的类

# 不调库实现Kmeans++
model = my_Kmeans_plus(k = 5, init = 'Kmeans++')
model.fit(data)
centers, ClusterAssment = model.predict(data)

'''
# 不调库实现Kmeans
model = my_Kmeans_plus(k = 5, init = 'random')
model.fit(data)
centers, ClusterAssment = model.predict(data)
'''

plt.figure(figsize=(6, 5))
plt.scatter(data[:, 0], data[:, 1], c = ClusterAssment[:, 0])
plt.scatter(centers[:, 0], centers[:, 1], marker = '*', color = 'r', s = 100) # 绘制类中心点
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

模型的稳定性提高了一些

---

SSE指标

SSE (sum of the squared errors，误差平方和)
$$SSE=\sum _{i=1}^k\sum _{x_j\in C_i}dist(x_j-{\mu }_i{)}^2$$
其中，C_i是第i个类，x_j是C_i中的样本点，μ_i是C_i的类中心(C_i中所有样本的均值)，dist是欧几里得距离，SSE是所有样本的聚类误差，用来评估聚类效果的好坏，可以理解为损失函数

因为我们要考虑如何更好地更新类中心，使得对于所有样本点，到其所属类的类中心距离之和最小
m i n S S E = ∑ i = 1 k ∑ x j ∈ C i ( x j − μ i ) 2 = ∑ i = 1 k ∑ x j ∈ C i ( x j T x j − x j T μ i − μ i T x j + μ i T μ i ) = ∑ i = 1 k ( ∑ x j ∈ C i x j T x j − ( ∑ x j ∈ C i x j T ) μ i − μ i T ( ∑ x j ∈ C i x j ) + m i μ i T μ i )

minSSE===​i=1∑k​xj​∈Ci​∑​(xj​−μi​)2i=1∑k​xj​∈Ci​∑​(xjT​xj​−xjT​μi​−μiT​xj​+μiT​μi​)i=1∑k​(xj​∈Ci​∑​xjT​xj​−(xj​∈Ci​∑​xjT​)μi​−μiT​(xj​∈Ci​∑​xj​)+mi​μiT​μi​)​
其中，m_i是第i个类中样本点个数
对类中心求导:
$$\frac{\partial SSE}{\partial {\mu }_i}=-(\sum _{x_j\in C_i}{x}_j)-(\sum _{x_j\in C_i}{x}_j)+2m_i{\mu }_i$$
令$\frac{\partial SSE}{\partial {\mu }_i}=0$
$${\mu }_i=\frac{1}{m_i}\sum _{x_j\in C_i}{x}_j$$
由上式可知，类的最小化SSE的最佳类中心点是类中所有样本点的均值

---

手肘法

前面提到Kmeans的缺点中，有提到过k值的选取对K-means影响很大。引用的这个数据集可以很容易看出来数据应该划分为5类，但不是所有数据集都能这么容易得出k值的大小，一种选取k值的办法就是通过不断尝试: 手肘法

手肘法的核心指标便是SSE

如图所示便是该数据的肘部曲线图，横坐标是k值的选取，纵坐标是误差平方和，曲线呈现递降的趋势，为什么？
当k=1时，表示不对样本点进行分类，直接认为所有样本点为同一类。当 k = numSamples，即k为样本点的个数，就是一个样本点分为一类，每一个类的类中心就是样本点自己本身，所以样本点到自己的类中心的距离都为0，SSE就为0

手肘法的核心思想是：随着聚类数k的增大，样本划分会更加精细，每个簇的聚合程度会逐渐提高，那么误差平方和SSE自然会逐渐变小。并且，当k小于真实聚类数时，由于k的增大会大幅增加每个簇的聚合程度，故SSE的下降幅度会很大，而当k到达真实聚类数时，再增加k所得到的聚合程度回报会迅速变小，所以SSE的下降幅度会骤减，然后随着k值的继续增大而趋于平缓，也就是说SSE和k的关系图是一个手肘的形状，而这个肘部对应的k值就是数据的真实聚类数。当然，这也是该方法被称为手肘法的原因

手肘法可视化代码实现

from pylab import *
mpl.rcParams['font.sans-serif'] = ['SimHei']

# 画肘部法曲线
K = []
inertia = []
for k_num in range(20): # 这里按理说应该是data.shape[0]，但数据的样本点太多，而且也没必要
    model = my_Kmeans_plus(k = k_num + 1, init = 'Kmeans++')
    model.fit(data)
    centers, ClusterAssment = model.predict(data)
    K.append(k_num + 1)
    inertia.append(np.sum(ClusterAssment[:, 1]))

plt.title('肘部曲线')
plt.xlabel('k')
plt.ylabel('SSE')
plt.plot(K, inertia)
plt.scatter(K, inertia, c='r', s=4)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19

调库实现Kmeans并绘制手肘法

from sklearn.cluster import KMeans
from pylab import *
mpl.rcParams['font.sans-serif'] = ['SimHei']

sse = []
for k_num in range(1, 20):
    model = KMeans(n_clusters=k_num, init="random", n_init=10, max_iter=200)
    model.fit(data)
    sse.append(model.inertia_)

plt.title('肘部曲线')
plt.xlabel('k')
plt.ylabel('SSE')
plt.plot(range(1, 20), sse)
plt.scatter(range(1, 20), sse, c='r', s=4)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

参考资源

• 手写算法-python代码实现Kmeans
• 手写算法-python代码实现Kmeans++以及优化
• 聚类算法之——K-Means++聚类算法