k 均值算法（k-means）

k-means 聚类算法，属于无监督学习算法。也就是说样本中却没有给定y，只有特征x。聚类的目的是找到每个样本x潜在的类别y，并将同类别y的样本x放在一起。k-means 算法实际上是一种最大期望算法（EM 算法）。

1. k-means算法步骤

在k-means算法中，用质心（样本均值）来表示簇类（cluster）；且容易证明k-means算法收敛等同于所有质心不再发生变化。基本的k-means算法流程如下：

---

设训练样本集是 D = { x 1 , x 2 , . . . , x m } D = \{x_1, x_2,..., x_m\} D={x1​,x2​,...,xm​}，每个$x_i\in R$，

1. 随机选取k个初始质心（作为初始cluster） U = { u 1 , u 2 , . . . , u k } U=\{u_1, u_2, ..., u_k\} U={u1​,u2​,...,uk​}，对应的分类（簇）为 C = { C 1 , C 2 , . . . , C k } C = \{C_1, C_2,..., C_k\} C={C1​,C2​,...,Ck​};
2. repeat：
   对于每个样本点，遍历所有k个质心，并计算得到距离该样本点 $x_i$ 最近的质心 $u_j$ ，将样本点 $x_i$ 的类别设置为质心 $u_j$ 所对应的cluster（$c_i$也就是$x_i$的分类）:
   $$c_i:=argmin||x_i-u_j|{|}^2$$
   更新完所有样本点的类别之后，遍历每个cluster中的所有样本点，计算该类下所有样本的均值向量，然后就得到了k个cluster对应的新的质心：
   u j : = ∑ i = 1 m 1 { x i = j } ⋅ x i ∑ i = 1 m 1 { x i = j } u_j := \frac{\sum_{i=1}^{m}1\{x_i=j\}\cdot x_i}{\sum_{i=1}^{m}1\{x_i=j\}} uj​:=∑i=1m​1{xi​=j}∑i=1m​1{xi​=j}⋅xi​​
   until 质心不再发生变化。

注：这里的 1 { x i = j } 1\{x_i=j\} 1{xi​=j}表示$x_i=j$的时候为1（也就是要加上$x_i$），否则为0.

---

注：划分所得簇的误差平方和（Sum of the Squared Error,SSE）为：

$$E=\sum _{j=1}^k\sum _{x\in C_j}{||x-u_j||}^2$$

这也是k-means的代价函数。算法找到最近质心点其实就是使得SSE最小化。

2. k-means 的优缺点

缺点

• k-means是局部最优的，容易受到初始质心的影响；
• 同时，k值的选取也会直接影响聚类结果，最优聚类的k值应与样本数据本身的结构信息相吻合，而这种结构信息是很难去掌握，因此选取最优k值是非常困难的。除了基于经验以及多次试验交叉验证，还可以是使用Gap Statistic($Gap=E(log{D}_k)-log{D}_k$,取Gap最大的k值)的方法。
• K均值算法并不适合所有的数据类型。它不能处理非球形簇、不同尺寸和不同密度的簇。面对非凸的数据分布形状时，可能新需要引入核函数来优化。
• 对离群点的数据进行聚类时，K均值也有问题，也就是对异常数据敏感。这种情况下，离群点检测和删除有很大的帮助。一般使用k-means之前需要对数据做预处理。

优点

• 实现简单
• 收敛快
• 虽然是局部最优结果，但一般情况就已经可以满足聚类需求了

3. K值的评估标准

不像监督学习的分类问题和回归问题，我们的无监督聚类没有样本输出，也就没有比较直接的聚类评估方法。但是我们可以从簇内的稠密程度和簇间的离散程度来评估聚类的效果。常见的方法有轮廓系数Silhouette Coefficient和Calinski-Harabasz Index。个人比较喜欢Calinski-Harabasz Index，这个计算简单直接，得到的Calinski-Harabasz分数值s越大则聚类效果越好。Calinski-Harabasz分数值s的数学计算公式是：
$$s(k)=\frac{tr(B_k)}{tr(W_k)}\frac{m-k}{k-1}$$

其中m为训练集样本数，k为类别数。Bk为类别之间的协方差矩阵，Wk为类别内部数据的协方差矩阵。tr为矩阵的迹。

也就是说，类别内部数据的协方差越小越好，类别之间的协方差越大越好，这样的Calinski-Harabasz分数会高。在scikit-learn中， Calinski-Harabasz Index对应的方法是metrics.calinski_harabaz_score.

4. 优化

3.1 二分k均值：bisecting k-Means

由于传统的KMeans算法的聚类结果易受到初始聚类中心点选择的影响，在基本k-means的基础上发展而来二分 (bisecting) k-means，其主要思想：一个大cluster进行分裂后可以得到两个小的cluster；为了得到k个cluster，可进行k-1次分裂。算法流程如下：

初始只有一个cluster包含所有样本点；
repeat:
    从待分裂的clusters中选择一个进行二元分裂，所选的cluster应使得SSE最小（可以按照原始k-means的方法进行质心选择）；
until 有k个cluster

二分k-means算法对初始质心的选择不太敏感，因为初始时只选择一个质心。

3.2 k-means++ 算法

针对原始k-means的初始值选择的改进是很重要的一部分，其中最具影响力的就是k-means++算法。k-means++选取k个点组委聚类中心的方法如下：

假设已经选取了n(0<n<k)个初始聚类中心，则在选取第n+1个聚类中心时，距离当前n个聚类中心越远的点会有更高概率被选为第n+1个聚类中心。此外，还是通过随机的方法选取第一个聚类中心。

选择完初始点后，k-means++的后续执行和原始k-means算法相同。

k-means+ 和 bisecting k-means 的不同在于：k-means++是先选取好k个聚类点，再进行聚类，而bisecting k-means是不断扩大分类的数量。

3.3 Canopy 算法

原始k-means的k值必须由用户给定。为了解决这个不足，可以先由Canopy（遮篷）算法粗略得到初始K值，然后再由k-means算法进行聚类。Canopy算法流程：

• (1) 将数据集向量化得到一个样本列表 $L=x_1,x_2,...,x_m$ ，（根据经验）选择两个距离阈值：T1和T2，其中T1 > T2，对应上图，外圈圈为T1，内圈为T2，T1和T2的值可以用交叉校验来确定；
• (2) 从list中任取一点P，用低计算成本方法快速计算点P与所有Canopy(聚簇中心)之间的最小距离$D(P,c_j)$（如果当前不存在Canopy，则把点P作为一个Canopy，这样就有了第一个聚簇中心，也就是一个分类），如果距离D<T1，则将点P加入到这个$c_j$对应的这个Canopy，此时暂不删除点P；
• 如果距离D<T2，则需要把点P从list中删除，这一步是认为点P此时与这个Canopy已经够近了，因此它不可能再属于别的聚簇了；
• 重复步骤2、3，再依次选出其他聚簇中心，直到list为空结束。

于是最终达到的效果如下：

3.4 ISODATA 算法

ISODATA(Iterative Selforganizing Data Analysis Techniques Algorithm) 的全称是迭代自组织数据分析法。ISODATA使用归并与分裂的机制，当某两类聚类中心距离小于某一阈值时，将它们合并为一类，当某类标准差大于某一阈值或其样本数目超过某一阈值时，将其分为两类。在某类样本数目少于某阈值时，需将其取消。如此，根据初始聚类中心和设定的类别数目等参数迭代，最终得到一个比较理想的分类结果

3.5 Mini Batch k-Means

在原始的K-means算法中，每一次的划分所有的样本都要参与运算，如果数据量非常大的话，这个时间是非常高的，因此有了一种分批处理的改进算法。 使用Mini Batch（分批处理）的方法对数据点之间的距离进行计算。

Mini Batch的好处：不必使用所有的数据样本，而是从不同类别的样本中抽取一部分样本来代表各自类型进行计算。n 由于计算样本量少，所以会相应的减少运行时间n 但另一方面抽样也必然会带来准确度的下降。

4. k-means 实践

先随机生成一个样本集，然后分别使用k-means和mini batch k-means进行聚类。观察k值取不同大小的聚类情况，并计算Calinski-Harabasz分数。

Python 代码实现

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from sklearn.cluster import KMeans
from sklearn.cluster import MiniBatchKMeans
from sklearn import metrics
from sklearn.datasets.samples_generator import make_blobs
# X为样本特征，Y为样本簇类别， 共1000个样本，每个样本2个特征，共4个簇，簇中心在[-1,-1], [0,0],[1,1], [2,2]， 簇方差分别为[0.4, 0.2, 0.2]
X, y = make_blobs(n_samples=1000, n_features=2, centers=[[-1,-1], [0,0], [1,1], [2,2]], cluster_std=[0.4, 0.2, 0.2, 0.2], 
                  random_state =9)
plt.scatter(X[:, 0], X[:, 1], marker='o')
plt.show()

for index, k in enumerate((2,3,4,5)):
    plt.subplot(2,2,index+1)
    # 执行聚类操作
    y_pred = KMeans(n_clusters=k, random_state=9).fit_predict(X)
    # 计算得分
    score= metrics.calinski_harabaz_score(X, y_pred)  
    # 画出不同簇，并用不同颜色标记
    plt.scatter(X[:, 0], X[:, 1], c=y_pred)
    plt.text(.99, .01, ('k=%d, score: %.2f' % (k,score)),
                 transform=plt.gca().transAxes, size=10,
                 horizontalalignment='right')
plt.show()

for index, k in enumerate((2,3,4,5)):
    plt.subplot(2,2,index+1)
    y_pred = MiniBatchKMeans(n_clusters=k, batch_size = 200, random_state=9).fit_predict(X)
    score= metrics.calinski_harabaz_score(X, y_pred)  
    plt.scatter(X[:, 0], X[:, 1], c=y_pred)
    plt.text(.99, .01, ('k=%d, score: %.2f' % (k,score)),
                 transform=plt.gca().transAxes, size=10,
                 horizontalalignment='right')
plt.show()

看一下实际的聚类效果以及不同k值的得分情况：

可以看到k取4的时候，两种聚类方的Calinski-Harabasz分数都是最大的。

---

参考：

1. 《机器学习（周志华）》：9.4.1 k均值算法
2. 《百面机器学习》：5.1 k均值聚类
3. 用scikit-learn学习K-Means聚类

---

THE END.