前言

本文主要对KMeans算法中的数学原理进行介绍，将不涉及或很少涉及代码实现，阅读前请确保已了解基本的机器学习相关内容和数学基础。

文章概述

本文从数学角度出发介绍KMeans算法及相关知识，主要内容包括：

• KMeans算法原理：算法工作逻辑及原理，簇内平方和的计算，以及总体平方和、KMeans算法的收敛性证明
• 聚类算法评估指标：重点介绍真实标签已知时的评估指标 V-Measure，真实标签未知的评估指标 Silhouette Coeffcient（轮廓系数）和 Calinski-Harabasz Index（卡林斯基-哈巴内斯指数）
• KMeans算法优化：KMeans++、WKMeans、elkan KMeans、Mini-Batch KMeans、Bisecting KMeans

KMeans算法本身原理较为简单，因此本文重点介绍其中包含的数学原理以及相关方法的推导和原理。

KMeans算法基本原理

基本工作流程

KMeans算法是处理聚类任务的一种，同时也是最简单的聚类算法，原理相对很容易就能理解，在讲解KMeans算法之前，先要了解一下 “簇” 和 “质心” 的概念。

簇：直观上来看簇就是一组以某种形式”聚集”在一起的数据，在一个簇中的数据就认为是同一类；簇就是聚类的结果表现。

质心：簇中所有数据的均值位置通常被称为这个簇的”质心”（centroid）；例如在一个二维平面中，某个簇的质心的横坐标就是这该簇中所有数据点的横坐标的均值，纵坐标就是该簇中所有数据点的纵坐标的均值，其它任意维度空间也同理。

【含有两个特征的某数据集使用K=4的KMeans聚类的一个迭代过程】

KMeans工作流程：在KMeans算法中，簇的个数K是一个超参数，需要我们人为输入来确定；KMeans的核心任务就是根据我们设定好的簇数K，找出K个最优的质心，并将离这些质心最近的数据分别分配到这些质心代表的簇中去，具体工作过程如下：

1. 对于某数据集，随机抽取K个样本作为最初的质心
2. 进入循环
   1. 将所有样本点分别分配到各自距离最短的质心
   2. 对于分好的每个簇，分别重新计算每个簇对应的新的质心
3. 当质心位置不再变化（或变化小于某一阈值）时，循环停止，聚类完成

聚类完成条件：当我们找到一个质心，在每次迭代中被分配到这个质心上的样本都是一致的，即每次新生成的簇都是一致的，所有的样本点都不会再从一个簇转移到另一个簇，质心就不会变化了。

基本工作原理

在生活中，我们把相似的东西归为一类，例如将 “青苹果” 分为 “苹果” 这一类而不是 “梨” 这一类，数据点也是如此，同一类别下的数据都在某方面有着一定的相似性，我们基于数据点之间的相似程度将不同数据点进行归类，这就是聚类算法的依据；那怎么衡量数据点之间的相似程度呢？

“距离” 衡量差异：不难想到，特征空间中两数据点之间的距离就代表着两个数据点之间的差异程度，对于某数据集而言，每个特征都是衡量数据的一个方向，而两数据点之间的距离是对于所有方向组成的特征空间而言的，其之间的距离当然可以代表它们的差异或相似程度；例如对于有长和宽两个特征的一组矩形，对于长和宽这两个值相差越小，我们则认为两个矩形差异越小；同理，对于拥有相似特征的数据点，在特征空间中 “挨得越近”，我们就认为差异越小，反之差异越大。

和人一样，把具有相似特征的东西划为一类；特征空间中的距离通过衡量不同数据点的各特征之间的差异来代表数据点之间的差异。

基于这种原理再回头看KMeans算法的工作流程：对于质心，质心是某组簇中各数据点的均值，在特征空间中处于该簇数据的正中心位置，用来代表每个簇；对于数据点，我们基于每个数据点到各个质心的距离，将其归类于距离最短的质心所代表的的簇。对于每次迭代，在假设每个质心的划分合理下，保证每个簇中的数据点都是相对最相似、最合理的，这时我们再回过头来检查质心的划分是否合理，也就是让每一类中新的质心替换掉原来的质心，进行更替的原因也很简单，因为对于每一个簇而言，旧的质心与该簇中所有样本的平均距离一定是大于新的质心与该簇中所有样本的平均距离的，新的质心更能代表这一个簇，因此当然要进行替换；随着质心的替换，簇中样本也要随之发生变化。

到此不妨把眼界放地更广一些，对于每次迭代的所有簇而言，不难发现各个簇中质心与样本的平均距离的总和一定是在减少的（下方会进行证明），当其达到最小值时，也就是模型收敛、聚类结束的时候。

簇内平方和 / 整体平方和

KMeans算法的目的：对于聚类的结果，我们追求 “簇内差异小，簇外差异大”，也就是同一簇下的数据点差异小，差异小，而不同簇下的数据点则希望差异大，差异尽量大。

• 对于一个簇来说，所有样本点到该簇质心的距离之和越小，我们就认为这个簇中的样本越相似，簇内差异就越小，这个簇聚得就越好。

上方提及过，特征空间中两数据点之间的距离越大，相似程度越低，差异越大，距离越小，则相似程度越高，差异越小；这种距离的衡量方法有多种，对于某特征空间，样本点到质心的距离可以由以下常用距离来度量：

$$\begin{array}{rl}& EuclideanDistance（欧氏距离）：d(x,\mu )=\sqrt{\sum _{i=1}^n(x_i-{\mu }_i{)}^2}\\ & ManhattanDistance（曼哈顿距离）：d(x,\mu )=\sum _{i=1}^n|x_i-{\mu }_i|\\ & CosineDistance（余弦距离）：cos\theta =\frac{{\sum }_{i=1}^n(x_i\ast {\mu }_i)}{\sqrt{{\sum }_{i=1}^n{x}_i^2}\ast \sqrt{{\sum }_{i=1}^n{\mu }_i^2}}\end{array}$$

• $n$：表示每个样本点中的特征数目
• $x$：表示簇中的一个样本点
• $\mu$：表示该簇中的质心
• $x_i$：表示该样本点中第 $i$ 个特征的样本值
• $\theta$：以向量角度来看的两样本向量之间的夹角

簇内平方和/整体平方和：对于每个样本点来说计算每个样本点与质心的距离就是两者中各个特征之间距离的累加；对于每个簇而言，某簇的簇内平方和就是该簇中各个样本点与质心之间距离的总和；对于整个聚类而言，整体平方和就是所有簇内平方和的累加；以欧式距离为例，计算簇内平方和与整体平方和的公式如下：

$$\begin{array}{rl}& 整体平方和：TotalInertia=\sum _{k=1}^{K}CS{S}_k\\ & 簇内平方和：ClusterSumofSquare(CSS)=\sum _{j=1}^m\sum _{i=1}^n(x_{ji}-{\mu }_i{)}^2\end{array}$$

• $K$：簇数
• $m$：簇中样本数
• $n$：表示每个样本点中的特征数目
• $x_{ji}$：某簇中第 $j$ 个样本点的第 $i$ 个特征对应的特征值
• ${\mu }_i$：某簇中质心的第 $i$ 个特征对应的特征值
• $CS{S}_k$：第 $k$ 个簇对应的CSS（簇内平方和）

Total Inertia（整体平方和）越小，代表着每个簇内样本越相似，说明在当前质心下聚类的效果就越好，因此KMeans算法的本质就是求解能够让 Inertia 最小化的质心。

在质心不断变化不断迭代的过程中，总体平方和是越来越小的（下方会进行证明），当整体平方和最小的时候，质心就不再发生变化了，该求解过程也就变成了一个最优化问题。在KMeans中，在确定簇数K后，最小化 Inertia 来求解最佳质心，并基于质心的存在重新进行聚类，直到模型收敛。

KMeans的损失函数：可以说 Total Inertia 是KMeans的损失函数，也可以说KMeans没有损失函数，业界什么说法都有，不必在这一点纠结。

损失函数的设定是为了求解某些参数，但是从上方介绍的过程和原理来看，KMeans并不是在进行像线性回归或逻辑回归那样的参数求解，它只是在不断迭代和更新质心位置，且不需要基于损失函数来更新质心，而是通过每簇的样本点来计算，因此当然不存在什么损失函数；但若硬要使用优化算法来求解KMeans的话，也不是不可以，例如可以使用梯度下降的方法来不断更新质心在各个方向上的位置，最后得到收敛，也就是说 Total Inertia 也可以是KMeans的损失函数；各有各的理，不必在此纠结过多。

质心与距离的选择：除了欧式距离，我们也可以使用其它距离度量方法，每种距离度量方法都有其对应的簇内平方和整体平方和，同时质心的选择也不止平均数这一种方法，在过去的经验中，我们总结出不同距离度量方法所对应的质心选择方法组合：

距离度量方法	质心选择方法
欧氏距离	均值
曼哈顿距离	中位数
余弦距离	均值

使用这些距离度量方法和对应的质心选择方式一般都可以达到不错的效果。

KMeans 收敛性证明

对于 Total Inertia (整体平方和)，我们在上方提到会随着迭代次数的增加逐渐减小，最终趋于收敛，但显然无法仅凭直觉来理解 Total Inertia 的收敛性，本文在此进行收敛性证明以供参考。

Total Inertia 的计算公式：

$$\begin{array}{rl}TotalInertia=\sum _{k=1}^{K}CS{S}_k& =\sum _{k=1}^K\sum _{j=1}^{m_k}\sum _{i=1}^n(x_{kji}-{\mu }_{ki}{)}^2\\ & =\sum _{k=1}^K\sum _{j=1}^{m_k}||x_{kj}-{\mu }_k|{|}_2^2\end{array}$$

• $K$：簇数
• $m_k$：第 k 个簇下的样本数
• $n$：样本对应的特征数
• $x_{kji}$：第 k 个簇下的第 j 个样本下的第 i 个特征对应的特征值
• ${\mu }_{ki}$：第 k 个簇对应质心的第 i 个特征对应的特征值
• $||x_{kj}-{\mu }_k|{|}_2$：第 k 个簇下的第 j 个样本与第 k 个簇对应质心的欧几里得范数

Total Inertia的计算时机：这里需要明确的是，Total Inertia 的每次计算发生在样本划分更新后，而不是计算出质心后进行计算，在计算出新的质心并对样本进行划分后，在计算下一次质心之前计算当前质心的 Total Inertia.

Total Inertia的单调性证明：KMeans算法共分为两步，第一步为对各个样本点基于就近原则分配给各距离最短的簇，第二步为基于各个簇的样本点重新计算每个簇的质心，我们在此用 $J$ 表示 Total Inertia，对于单个簇来讲，则有：

• 若簇的质心位置不变，则 Total Inertia 不增.
• 若簇的质心发生改变，对于质心更新前、样本划分更新前，旧质心所在簇的簇内平方和为最小，对于质心更新前、样本划分更新后，由于旧质心已经在样本划分更新后已经不是各样本距离的均值了，因此以旧质心计算的簇内平方和也已经不是最小值了，对于质心更新后，样本划分更新后，由于新质心为样本划分后簇内各样本的距离均值，因此以旧质心计算的簇内平方和一定小于新质心计算的簇内平方和.

基于以上理论，因此对于每个簇的质心的每次更新有：

$$\{\begin{array}{l}\hat{J}=J-{\sum }_{j=1}^{m_k}||x_{kj}-{\mu }_k|{|}_2^2+{\sum }_{j=1}^{m_k}||x_{kj}-{\mu }_k^{\prime }|{|}_2^2\\ {\mu }_k|{|}_2^2+{\sum }_{j=1}^{m_k}||x_{kj}<={\sum }_{j=1}^{m_k}||x_{kj}-{\mu }_k^{\prime }|{|}_2^2\end{array}⟶\hat{J}\le J$$

• $m_k$：第 k 个簇下的样本数
• $x_{jk}$：第 k 个簇下的第 j 个样本
• ${\mu }_k$：第 k 个簇对应质心
• $||x_{kj}-{\mu }_k|{|}_2$：第 k 个簇下的第 j 个样本与第 k 个簇对应质心的欧几里得范数

Total Inertia存在极小值证明：至此，单调性证明完毕，可以发现 Total Inertia 单调递减。

将 ${\mu }_k$ 作为自变量，Total Inertia 作为因变量，对该函数求导得：

$$J^{\prime }=-2\sum _{k=1}^K\sum _{j=1}^{m_k}\sum _{i=1}^n(x_{kji}-{\mu }_{ki})$$

• $K$：簇数
• $m_k$：第 k 个簇下的样本数
• $n$：样本对应的特征数
• $x_{kji}$：第 k 个簇下的第 j 个样本下的第 i 个特征对应的特征值
• ${\mu }_{ki}$：第 k 个簇对应质心的第 i 个特征对应的特征值
• $||x_{kj}-{\mu }_k|{|}_2$：第 k 个簇下的第 j 个样本与第 k 个簇对应质心的欧几里得范数

注意：参考上方说明的 Total Inertia 的计算时机，这里的 ${\mu }_{k_i}$ 并不是当前簇中所有样本距离的均值，不存在 $m_k{\mu }_{ki}={\sum }_{j=1}^{m_k}{x}_{kji}$，也就是 $J^{\prime }$ 并不是恒为0.

不难发现，对于函数 $J^{\prime }$，在样本更新前，在某簇中对于各个样本点一定有各个特征方向距离的均值与质心的各个特征方向距离相等，也就满足 $J^{\prime }=0$，在本轮样本划分前，新的质心更新前的样本和质心分布下，对于单个簇中有：

$$\begin{array}{r}\frac{1}{m}\sum _{j=1}^m{x}_{ji}={\mu }_i\end{array}$$

• $x_{ji}$：某簇中的第 j 个样本对应的第 i 个特征的特征值
• ${\mu }_i$：某簇中的质心的第 i 个特征的特征值

注意这里是样本划分前，样本划分前 ${\mu }_i$ 仍然是该簇中各样本距离的均值

因此就存在该值满足 $J^{\prime }=0$，即函数 $J$ 存在极小值。

收敛性证明：因为该函数单调递减且存在极小值，又因为该函数的二阶导数为一个常数，则函数 $J$ 当然为一个严格意义上的凸函数，最终当然会趋于收敛。

聚类算法的评估指标

聚类模型的结果不是某种标签输出，并且聚类的结果是不确定的，其优劣由业务需求或者算法需求来决定，并且不存在一定的正确答案。

不能使用 Total Inertia 作为评估指标：观察 Total Inertia 的计算公式不难看出，它是用来衡量簇内差异的指标，会严重受到超参数K的影响，随着超参数K的增大，也就是随着簇的数量的增大，Total Inertia 的值必然会越来越小，但这并不意味着模型的效果越来越好了，因此不能用于作为评估指标。

聚类好坏的评估标准：簇内差异小，簇外差异大。

Total Inertia 只衡量簇内差异，不衡量簇外差异，因此随着簇数 K 的增大其必然越来越小，但显然随着簇数 K 的增大簇外差异（簇与簇之间的差异）其实也在不断减小，我们要找的评估指标要同时涵盖簇内和簇外两方面的差异才行。

虽然我们在聚类中不输入真实标签，但这不代表我们拥有的数据中一定不具有真实标签，或者一定没有任何参考信息；因此考虑到含不含有真实标签，本文在此主要介绍以下不同情况的常用评估指标：

• 真实标签已知：
  • V-Measure：同质性度量和完整性度量的调和平均
• 真实标签未知：
  • 轮廓系数
  • 卡林斯基-哈巴内斯系数

在实际应用时，拥有真实标签的情况非常少见，若有真实标签，我们更倾向于使用分类算法而不是聚类算法，但不排除我们依然可能使用聚类算法的可能性，但通常来讲，真实标签未知时的评估指标才是最常用的。

真实标签已知时的评估指标

V-Measure 对样本分布没有假设，在任何分布上都有不错的表现；除了以下介绍的V-Measure外，互信息分、调整兰德系数等其它相关指标也很常见，本文在此主要介绍以下两种指标，在绝大部分场景下，V-Measure已经可以满足需求。

聚类后的各个簇与真实标签：对于包含某真实标签样本最多的某个簇，我们一般选择该簇作为该真实标签对应的结果。

注：本文下方对于类别和簇，若不指明，则为以下含义：

• 类别：数据样本对应的真实标签
• 簇：聚类后各个数据样本对应的簇（标签）

V-Measure

计算不纯度

对于信息熵我们都很熟悉，信息熵是可能发生的事件的对应信息量的期望，因此信息熵可以被用来衡量不确定度，当信息熵越大时，不确定度或不纯度越大，信息熵公式如下：

$$Entropy=-\sum _{c=1}^C{P}_{c}lo{g}_{2}P(c)$$

• $C$：类别数
• $P_c$：第 c 个可能发生的事件的概率或第 c 组样本占样本总数的比例

若不熟悉信息量与信息熵，可以参考：信息量与信息熵_信息量和信息熵-CSDN博客

已知共有 C 个类别，则我们可以计算出在原数据 C 个类别下的不纯度（熵） $H(C)$，公式如下：

$$\begin{array}{r}H(C)=-\sum _{c=1}^C\frac{n_c}{n}log(\frac{n_c}{n})\end{array}$$

• $n$：样本总数
• $C$：类别数
• $n_c$：类别 c 下的样本数

得到原数据样本整体的不纯度后，接下来我们需要求出对于聚类后、在 K 个簇划分条件下的原类别划分对应的不纯度，这么讲有点绕口，不过没有关系，这里将会一步一步来进行推导。首先我们要先得到聚类后对于单个簇的不纯度（聚类后的单个簇可能包含多个类别的样本），公式如下（这里简称CI，Cluster Impurity）：

$$CI=-\sum _{c=1}^{C}P(c|k)logP(c|k)=-\sum _{C=1}^C\frac{n_{c,k}}{n_k}log(\frac{n_{c,k}}{n_k})$$

• $C$：类别数
• $k$：聚类后的第 k 个簇
• $n_k$: 簇 k 下的样本数
• $n_{c,k}$：类别 c 中被划分到簇 k 中的样本数

对于聚类后的所有簇，我们对每一个簇对应的不纯度取期望即可，得到整体的在 K 个簇划分条件下的原类别划分对应的不纯度，公式如下：

$$H(C|K)=\sum _{k=1}^K\frac{n_k}{n}C{I}_k=-\sum _{k=1}^K\sum _{c=1}^C\frac{n_k}{n}\ast \frac{n_{c,k}}{n_k}log(\frac{n_{c,k}}{n_k})=-\sum _{k=1}^K\sum _{c=1}^C\frac{n_{c,k}}{n}log(\frac{n_{c,k}}{n_k})$$

• $n$：样本总数
• $C$：类别数
• $K$：簇数
• $n_c$：类别 c 下的样本数
• $n_k$: 簇 k 下的样本数
• $n_{c,k}$：类别 c 中被划分到簇 k 中的样本数

显然，$H(C)$ 为某数据集下的特定值；由于我们追求聚类后每个簇的不纯度越小越好，也就是每个簇中对于所包含类别的混乱程度越小越好，放到整体来看，也就是所有簇的不纯度的期望、也就是 $H(C|K)$ 越小越好。

此外，对于聚类后的每个簇，当每个簇中包含的样本都是相同类别的，这时候每个簇的不纯度为 0，$H(C|K)=0$，当每个簇中包含各个类别的样本，且占比相同（最混乱的情况），此时不纯度最高，有如下推导结论：

$$\begin{array}{rl}H(C|K)& =-\sum _{c=1}^C[\frac{n_{c,1}}{n}log(\frac{n_{c,1}}{n_1})+\frac{n_{c,2}}{n}log(\frac{n_{c,2}}{n_2})+\cdots +\frac{n_{c,K}}{n}log(\frac{n_{c,K}}{n_K})]\\ & =-\sum _{c=1}^C[\frac{n_{c,1}}{n}log(\frac{n_c}{n})+\frac{n_{c,2}}{n}log(\frac{n_c}{n})+\cdots +\frac{n_{c,K}}{n}log(\frac{n_c}{n})]\\ & =-\sum _{c=1}^C\frac{n_c}{n}log(\frac{n_c}{n})=H(C)\end{array}$$

要保证整体”最混乱”，即不纯度最高，则各簇中的各类别样本应占比相同且等于每个类别占总体的比例（这里就不详细证明了），即：

$$\frac{n_c}{n}=\frac{n_{c,1}}{n_1}=\frac{n_{c,2}}{n_2}=\cdots =\frac{n_{c,K}}{n_K}$$

即 $H(C|K)$ 的取值范围为 $[0,H(C)]$。

计算思路与以上同理，我们也可以求出 $H(K)，H(K|C)$，公式总结如下：

$$\begin{array}{r}H(C)=-\sum _{c=1}^C\frac{n_c}{n}log(\frac{n_c}{n})\\ H(K)=-\sum _{k=1}^K\frac{n_k}{n}log(\frac{n_k}{n})\end{array}，\begin{array}{r}H(C|K)=-\sum _{k=1}^K\sum _{c=1}^C\frac{n_{c,k}}{n}log(\frac{n_{c,k}}{n_k})\\ H(K|C)=-\sum _{c=1}^C\sum _{k=1}^K\frac{n_{k,c}}{n}log(\frac{n_{k,c}}{n_c})\end{array}$$

• $n$：样本总数
• $C$：类别数
• $K$：簇数
• $n_c$：类别 c 下的样本数
• $n_k$: 簇 k 下的样本数
• $n_{c,k}$：类别 c 中被划分到簇 k 中的样本数

四种不纯度的解释总结如下（这里的名字是本文随便取的，不权威）：

• $H(C)$：类别不纯度（类别划分熵）；也就是原数据集中按类别划分对应的不纯度.
• $H(K)$：聚类不纯度（聚类划分熵）；也就是原数据集中按聚类后的簇划分对应的不纯度.
• $H(C|K)$：给定簇划分条件下的类别划分熵；也就是按聚类后的簇对样本数据进行划分后，按类别计算每个簇的不纯度，然后取所有簇的不纯度的期望，取值范围为 $[0,H(C)]$.
• $H(K|C)$：给定类别划分条件下的聚类划分熵；也就是按类别对样本数据进行划分后，按聚类后的簇分贝计算每个类别的不纯度，然后取所有类别的不纯度的期望，取值范围为 $[0,H(K)]$.

Homogeneity & Completeness & V-Measure

Homogeneity（同质性）度量：考虑聚类后每个簇中只包含某单个类别的样本的程度；取值范围为 $[0,1]$，其值越接近1说明效果越好，越接近0说明效果越差；本文下方称其为 h 值，计算公式如下：

$$h=1-\frac{H(C|K)}{H(C)}$$

Completeness（完整性）度量：考虑给定某类别的所有样本，其都被划分给同一个簇的程度；取值范围为 $[0,1]$，其值越接近1说明效果越好，越接近0说明效果越差；本文下方称其为 c 值，计算公式如下：

$$c=1-\frac{H(K|C)}{H(K)}$$

同质性越强，则聚类后每个簇越倾向于仅包含单个类别；完整性越强，则相同类别下的样本越倾向于划分到同一个簇；我们要做的就是从这两个度量方式中做出权衡。

试想一下两种极端情况，假设共有 $n$ 个样本，10个类别，则将这 $n$ 个样本聚成 $n$ 个簇得到的Homogeneity（同质性）是最强的，因为每个簇都仅包含1个样本，也就是每个簇中仅包含同一个类别的样本；同理，将这 $n$ 个样本聚成 1 个簇得到的Completeness（完整性）是最强的，因为这样的话无论给定哪个类别的样本，都会被划分到同一个簇中；因此我们要做的是从这两个度量方式中取得权衡，由此得到我们的聚类评估指标 V-Measure。

V-Measure：Validity Measure，用于权衡 h 值和 c 值，是两者的调和平均数，取值范围为 $[0,1]$，越接近1说明聚类效果越好，反之越差；计算公式如下：

$$V=\frac{1}{\frac{1}{2}\ast (\frac{1}{h}+\frac{1}{c})}=\frac{2hc}{h+c}$$

• 调和平均数：所有数字的倒数的算术平均数的倒数

为什么使用调和平均数：调和平均数对异常值更不敏感，当 h 值与 c 值大小差不多时，调和平均数与算数平均数之间差异不大；但当 h 值与 c 值大小相差较大时，这时若使用算数平均数，V-Measure的值将会被其中较大的一方主导而导致V-Measure较大，但实际上我们要做的是权衡两者，当其中一方较大时，平均数应该更偏向值比较小的那一方，因此我们在这里选择调和平均数。

其它常见指标

对于其它常见指标本文在此只做简单介绍。

互信息分：取值范围为 $[0,1]$，其中越接近 1 表示聚类效果越好，反之聚类效果越差

• MI：Mutual Info，普通互信息分
• AMI：Adjusted Mutual Info，调整的互信息分
• NMI：Normalized Mutual Info，标准化互信息分

调整兰德系数：ARI（Adjusted Rand Index），取值范围为 $[-1,1]$，越接近 1 表示聚类效果越好，越接近-1表示聚类效果越差；对样本分布没有假设，在任何分布上都可以有不错的表现，尤其是在具有”折叠”形状的数据上表现优秀。

真实标签未知时的评估指标

对于聚类算法，一般真实标签都是未知的，本文在此重点介绍 SC（轮廓系数）和CHI（卡林斯基-哈拉巴斯指数），除了这两个评估指标外，戴维斯-布尔丁指数（Davies-Bouldin）和权变矩阵（Contingency Martix）等也是常用的指标，本文在此重点介绍前两者，SC和CHI一般已经可以应对绝大部分场景。

对于真实标签未知的聚类效果评估，完全依赖于 评价簇内的稠密程度和簇外的离散程度，也就是 “簇内差异小，簇外差异大”.

Silhouette Coefficient

Silhouette Coefficient：轮廓系数；轮廓系数是对于单个样本而言的，对于某样本的轮廓系数公式如下：

$$SC=\frac{b-a}{max(a,b)}=\{\begin{array}{l}1-\frac{a}{b},a<b\\ 0,a=b\\ \frac{b}{a}-1,a>b\end{array}$$

• $a$：某样本点与该样本所在簇中其它样本点之间的平均距离，被称为该样本与其所在簇中其它样本之间的相似度
• $b$：某样本点与和该样本相距最短的其它簇中的样本点之间的平均距离，被称为该样本与其它簇中样本之间的相似度

为了衡量簇内差异小，簇间差异大，则相似度 $a$ 越小越好，相似度 $b$ 越大越好，也就是轮廓系数的值越大越好；此时不难发现，轮廓系数的取值范围为 $[-1,1]$，对于某单个样本而言，当轮廓系数越趋近于1，代表效果越好，反之若越趋于-1，则结果越差，其中若某样本的轮廓系数为负数，则说明该样本的聚类结果非常差。

若一个簇中的大多数样本整体都具有比较高的轮廓系数，则该簇会有较高的总轮廓系数，则整个数据集中所有样本的平均轮廓系数越高，说明聚类效果越好；同理，若存在具有低轮廓系数甚至负值的样本点个数越多，则说明聚类效果越差。

轮廓系数对数据分布没有假设，在很多数据集上都表现良好，其中在每个簇的分割比较清楚时表现最好；但其对于簇结构为凸的数据或使用例如DBSCAN的基于密度的聚类算法，轮廓系数较高（存在一部分虚高）。

凸数据集：在某数据集构成的特征空间中，任何两样本点之间连线也在该数据集所构成的集合内；最常见的例子就是环形数据集，例如一个环包着另一个环，对于这种数据集，考虑到轮廓系数中的每个样本天生地 $a$ 值就相对比较小，$b$ 值相对比较大，因此当然会虚高。

Calinski-Harabasz Index

Calinski-Harabasz Index：卡林斯基-哈拉巴斯指数，也被称为方差比标准；本质是簇间距离与簇内距离的比值，对于有 K 个簇的聚类而言，计算公式如下：

$$CHI==\frac{Tr(B_k)}{K-1}/\frac{Tr(W_k)}{N-K}=\frac{Tr(B_k)}{Tr(W_k)}\ast \frac{N-K}{K-1}$$

• $B_k$：簇间离散矩阵，即不同簇之间的协方差矩阵
• $W_k$：簇内离散矩阵，即某簇内数据的协方差矩阵
• $Tr(B_k)$：矩阵 $B_k$ 的迹，也就是正对角线元素相加
• $N$：数据集中的样本总数
• $K$：簇数

对于 $\frac{Tr(B_k)}{K-1}$：簇间离散矩阵 $B_k$，已知各个簇的质心（中心点），则可以通过各个簇的质心计算出整个数据集的中心点，$Tr(B_k)$ 其实就是各个簇的质心与数据集中心点之间距离的平方和，除以自由度 $K-1$，得到各个簇的质心与数据集中心点之间的平均距离。

对于 $\frac{Tr(W_k)}{N-K}$：对于簇内离散矩阵 $W_k$，中心点也就是质心，$Tr(B_k)$ 其实就是某簇中各个样本点与质心之间距离的平方和，除以自由度 $N-K$，得到某簇中样本点与其质心之间的平均距离。

相比轮廓系数，CHI没有界，且对于簇结构为凸的数据或使用例如DBSCAN的基于密度的聚类算法也会存在一部分虚高，但它相对轮廓系数在计算速度上有着巨大的优势，计算地很快。

KMeans原理拓展与算法优化

KMeans++

若有足够的时间，Kmeans算法最终一定会收敛，但却可能收敛到局部最小值；KMeans是否能够收敛到全局最小值很大程度上取决于质心位置的初始化。

KMeans++：KMeans++由 Arthur, David, Sergei Vassilvitskii 三人在2007年发表的论文 “k-means++: The advantages of careful seeding” 提出；KMeans++ 为 KMeans 算法的改进版本，不过相对KMeans，KMeans++只是在质心位置的初始化上进行了改进，其它地方同KMeans算法一样，算是对质心初始化问题的一个解决方案。

原始质心初始化方法：初始时随机选择 $K$ 个样本点作为初始质心位置。

KMeans++流程：由于逐个选取 $K$ 个质心，且离其它质心越远的样本点越有可能被选为下一个质心，流程如下：

1. 从原数据集中随机选择一个样本点作为第一个质心
2. 计算数据集中每个样本点与其最近的质心的距离，其中第 i 个样本点与其最近的质心的距离记为 $D(x_i)$
3. 开始循环：
   1. 计算数据集中每个样本点被选为下一个质心的概率：
      $$P(x_i)=\frac{D(x_i{)}^2}{{\sum }_{i=1}^{n}D(x_i{)}^2}$$
      • $n$：数据集中的样本总数
      • $x_i$：数据集中第 i 个样本点
      • $D(x_i)$：第 i 个样本点与其最近的质心的距离
   2. 已知每个样本点被选为下一个质心的概率后，按照算出来的概率进行随机选择，确定下一个质心
4. 直到找出 $K$ 个质心，循环结束.

KMeans++原理：由于篇幅受限且无必要，本文在此抛开数学证明及数学原理，以理解的形式回答以下两个问题：

• 为什么要使得样本点离其它质心越远，其被选择为下一个质心的概率越大？
  • 簇是有向心性的：同一个质心附近的点都会被纳入这个簇的范围内，反过来说也就是离最近的质心比较远的点属于不同簇的可能性比离得近的点大.
• 以什么方法实现下一个质心的选择？
  • 轮盘法：也就是用程序算法模拟一个按各个概率划分对应占比面积的轮盘，实现以各个概率进行抽样的效果.

若想了解严格的理论证明和原理，可以直接阅读KMeans++提出的原论文文章

WKMeans

WKMeans：首次在2005年的论文 “Automated Variable Weighting in k-Means Type Clustering” 被提出，相对一般的 KMeans算法，由于各个特征的重要性不同，它对每个特征都考虑相应的权重，在已知各个特征权重的情况下，在计算样本点间距离的时候以相应权重计算加权距离，因此得到的效果理论上要比一般的KMeans算法更为精确。

例：假设有三个特征，两个样本点坐标分别为$\{\begin{array}{l}{x}_1=[1,2,3]\\ x_2=[4,4,2]\end{array}$，已知各个特征对应的权重为 $W=[w_1,w_2,w_3]=[0.50,0.40,0.10]$，以欧式距离（平方）为例，则不考虑加权的距离 $d$ 和 考虑加权的距离 $d_w$ 计算如下：

$$\begin{array}{rl}d& =(1-4{)}^2+(2-4{)}^2+(3-2{)}^2=14\\ d_w& =0.5\ast (1-4{)}^2+0.4\ast (2-4{)}^2+0.1\ast (3-2{)}^2=6.2\end{array}$$

簇内平方和 / 整体平方和：WKMeans相对KMeans仅仅是在计算每个特征对应距离时加上了相应权重，计算公式如下：

Total Inertia=∑k=1KCSSkCluster Sum of Square(CSS)=∑j=1m∑i=1nwi(xji−μi)2约束条件：∑i=1nwi=1begin{align}
& Total space Inertia = sum^K_{k=1}CSS_k \

& Cluster space Sum space of space Square(CSS) = sum^m_{j=1} sum^n_{i=1}w_i(x_{ji}-mu_i)^2 \

& 约束条件：sum^n_{i=1}w_i = 1
end{align}​Total Inertia=k=1∑K​CSSk​Cluster Sum of Square(CSS)=j=1∑m​i=1∑n​wi​(xji​−μi​)2约束条件：i=1∑n​wi​=1​​

• $K$：簇数
• $m$：簇中样本数
• $n$：表示每个样本点对应的特征数目
• $x_{ji}$：某簇中第 $j$ 个样本点的第 $i$ 个特征对应的特征值
• ${\mu }_i$：某簇中质心的第 $i$ 个特征对应的特征值
• $CS{S}_k$：第 $k$ 个簇对应的CSS（簇内平方和）
• 约束条件：各个特征对应权重的和为1

WKMeans只是对各个特征设置了对应特征权重，对于权重的确定则需要基于其它方法来确定，例如使用相关树模型得到的特征重要性、基于熵权法得到的权重等等，在确定了权重后才能使用WKMeans。

elkan KMeans

elkan KMeans：在传统的K-Means算法中，我们在每轮迭代时，要计算所有的样本点到所有的质心的距离，这样会比较的耗时，elkan KMeans算法基于该过程入手加以改进，以减少不必要的距离的计算，加速模型迭代。

elkan KMeans原理：利用三角形的性质，即两边之和大于第三边，两边之差小于第三边来减少不必要的距离计算，例如对于两个质心 $u_1,u_2$ 和 某样本点 $x$，对于传统KMeans算法，则是先分别计算样本点 $x$ 与 $u_1,u_2$ 的距离，然后通过比较两距离来选择距离较短的那一个；而对于elkan KMeans，则是先计算 $u_1,u_2$ 之间的距离 $d(u_1,u_2)$，然后计算 $x$ 与其中一个质心的距离，这里以计算 $x$ 与 $u_1$ 的距离为例，记为 $d(x,u_1)$，则由三角形两边之差小于第三边可得：

$$d(u_1,u_2)-d(x,u_1)<d(x,u_2)$$

设 $d(u_1,u_2)=nd(x,u_1)$，则有：

$$(n-1)d(x,u_1)<d(x,u_2)$$

易得当 $n-1>1$ 时，即 $d(u_1,u_2)>2d(x,u_1)$时，有：

$$d(x,u_2)>d(x,u_1)$$

此外，由两边之差小于第三边可得，即$d(x,u_2)>|d(u_1,u_2)-d(x,u_1)|$，因此也有以下推论：

$$\{\begin{array}{l}d(x,u_2)>|d(u_1,u_2)-d(x,u_1)|\\ d(x,u_1)\le |d(u_1,u_2)-d(x,u_1)|\end{array}⟶d(x,u_1)<d(x,u_2)$$

总结一下，有如下两个推论：

$$\{\begin{array}{l}d(x,u_1)>2d(u_1,u_2)\to d(x,u_2)>d(x,u_1)\\ d(x,u_1)\le |d(u_1,u_2-d(x,u_1))|\to d(x,u_2)>d(x,u_1)\end{array}$$

将两个推论合并在一起表示为：

$$d(x,u_1)>2d(u_1,u_2)\mathbf{or}d(x,u_1)\le |d(u_1,u_2-d(x,u_1))|⟶d(x,u_2)>d(x,u_1)$$

基于以上推论，我们在已知两质心之间距离后，只需要计算一次样本点与质心距离后我们就可以直接通过比较该距离与两倍的质心之间距离就可以间接得到哪个距离更近了，而各质心之间距离只需要计算一次，不需要重复去计算，因此可以节省大量算力。

elkan KMeans算法具体流程：对于计算每个样本的所属质心(簇)的算法部分逻辑，对于每次迭代则有：

1. 已知有 $k$ 个质心，首先要计算各质心中两两之间的距离，要得到共需要计算多少次，将其看作一个组合数问题，也就是共有 $k$ 个点，两个为一组，判断共有多少个组合，也就是 $C_k^2$.
2. 将每个样本随机分配给一个质心，并分别计算每个样本与所分配到的质心之间的距离
   1. 对每个样本进行遍历，以某个样本为例，计算其与所分配到的质心的距离并记为 $d$
   2. 已知各质心之间距离，对于剩余质心进行遍历，这里记为 d_i$, 基于上方推论我们可以不计算该样本与其它质心的距离而直接比较两者大小
   3. 将遍历的每个质心与该样本间的距离记为 $d_i$，对于遍历的每个质心则有：
      • 若已知 $d_i<d$，则计算得到 $d_i$ 后，令 $d=d_i$，将该质心对原来的质心进行替换
      • 若已知 $d_i>d$，则不做处理
3. 确定每个样本与其所属质心(簇)后，剩下的操作与原始KMeans一致

elkan KMeans 与传统KMeans的计算量比较：设某数据集下的样本总数为 $n$，确定的簇数（质心数）为 $k$，共需要进行 $m$ 此迭代，则有：

• 传统KMeans算法：在每轮迭代中，需要对每个样本分别计算与 $k$ 个质心间的 $k$ 个距离，因此总共需要计算的距离次数为 $mnk$.
• elkan KMeans算法：对于每轮迭代，有：
  • 质心计算：已知有 $k$ 个质心，首先要计算各质心中两两之间的距离，要得到共需要计算多少次，将其看作一个组合数问题，也就是共有 $k$ 个点，两个为一组，判断共有多少个组合，也就是 $C_k^2$.
  • 最坏的情况：对于每个样本，每个质心在计算时都比上一个要小，则共需计算的距离次数为 $m(C_k^2+nk)$
  • 最好的情况：对于每个样本，每个质心在计算时第一个就是最小的，则共需计算的距离次数为 $m(C_k^2+n)$
  • 平均情况：$\frac{m(C_k^2+nk)+m(C_k^2+n)}{2}=m(C_k^2+\frac{k+1}{2}n)$

为了简单（图方便），这里的平均情况直接用 (最坏情况+最好情况)/2 这样的方式计算了，与实际平均可能有误差（实际上是所有情况加起来/情况总数），本文在这里主要为了表达elkan KMeans相对能节省大量计算资源。

如图所示，以上三图从左到右：

• 图左：样本总量为1000、簇数为5时，随着迭代次数的增加，两算法中距离计算次数的变化折线图
• 图中：迭代次数为100、簇数为5时，随着样本总量的增加，两算法中距离计算次数的变化折线图
• 图右：迭代次数为100、样本总量为1000时，随着簇数的增加，两算法中距离计算次数的变化折线图

从图中可以很容易看出（注意y轴单位为 $10^6$），随着迭代次数、样本总量、簇数的逐渐增大，elkan算法相对原始KMeans可以节省大量的计算资源，特别对于数据量大的时候非常有用，即使数据量较小时也仍然有一定的提升。

（补充说明）Lloyd算法：该算法与一般KMeans算法的计算思路如出一辙，也是通过逐点计算距离的方式，计算流程也基本相同；但最主要的区别在于KMeans在每次迭代中根据划分后的样本来更新质心，而Lloyds则事先根据划分后的样本对各个区域进行划分，前者是一个由许多离散的样本点构成的区域，而后者则是一个连续的区域，Lloyds将这些连续区域看作Voronoi单元，在更新时需要计算每个单元的质心，而不是像KMeans那样简单计算均值。

很多时候会将Lloyd算法替代原始KMeans的使用，在算法比较的时候不妨也可以将Lloyd看作原始的KMeans算法；此外，由于逐点的计算方式，Lloyd算法相比elkan算法在内存上的使用会更少，但计算效率要弱于elkan算法，当内存受限时，可以选择Lloyd算法，但无论是在大样本数据还是小样本数据上，elkan算法都在有着相同程度的准确性下有着更快的速度。

Mini-Batch KMeans

Mini-Batch KMeans：Mini-Batch KMeans 是用来应对数据量太大的情况的，其原理非常简单，当数据集样本量过大时，就从整个数据集中随机抽样得到样本量更小的一个数据子集（Batch），然后使用该数据子集对模型进行训练和迭代，在大数据量时，这样可以很大程度上节省时间，模型效果最终会略差于使用全部数据训练的KMeans模型，但大多情况下相差不大。

Mini-Batch KMeans算法一般流程：基于上方选择原数据集中的一部分样本作为数据子集来对模型进行训练的思想，算法一般流程如下：

1. 从原数据集中抽取部分样本组成一个batch，使用该batch样本数据训练KMeans模型，得到 $K$ 个簇和其对应质心
2. 开始循环
   1. 对原数据中抽取一个batch的样本数据，将其添加到上方的KMeans模型中，将这些样本分配给距离其最近的簇(质心)
   2. 更新质心的位置
3. 直到各质心位置稳定或达到最大迭代次数，循环结束

Mini-Batch KMeans算法的原理非常简单，但在大数据集的场景下却十分有效，KMeans对计算资源的要求很高，当数据集过大时，使用Mini-Batch的方法能够节省大量时间和计算资源。

Bisecting KMeans

Bisecting KMeans (常被称为二分KMeans算法)：属于层次聚类中分裂法的一种，主要是为了克服KMeans算法易收敛于局部最小值的问题。

Bisecting KMeans核心思想：若某簇分地很差，则说明该簇很大可能包含多个类别，因此我们在最开始将所有数据看成一个簇，然后不断地将每轮中分地最差的簇进行分裂（一分为2），直到达到目标条件（如目标簇数）即可。

可以使用簇内平方和来反映某簇分得好坏，由于簇内平方和无界，因此在每轮中可以通过对所有簇进行横向对比来找出分得最差的簇，对该簇使用KMeans算法进行分裂。

Bisecting KMeans一般流程：

1. 把所有数据初始化为一个簇，并计算该簇质心位置.
2. 开始循环
   1. 计算当前所有簇的簇内平方和，选择簇内平方和（或SSE）最大的簇
   2. 通过KMeans算法将选择的簇进行多次分类（目标簇数为2的距离），得到多组分裂结果
   3. 分别计算每组分裂结果的总簇内平方和（或总SSE）；每组的分裂结果为两个簇，对于每组的总簇内平方和也就是分裂后得到的两个簇的簇内平方和的和
   4. 选择多组分裂中总簇内平方和最小的结果作为本次分裂结果
3. 通过分裂达到事先设置的簇数或其它限制条件，循环结束

过程不唯一，只要符合Bisecting KMeans的核心思想，能达到优化目的即可。

待分裂簇的选择：除了上方所介绍的结果的总簇内平方和，一般也可直接计算当前所有簇对应的簇内平方和（或SSE），选择拥有最大簇内平方和的簇作为待分裂簇；此外，因为所有样本都要参与计算，会占用较多计算资源，模型一般也提供另一种策略，也就是直接选择所含样本量最大的簇作为待分裂簇，因为往往样本量大的簇对应的簇内平方和（或SSE）也大，这种策略会丢失一部分精度，但可节省许多计算资源，且在实践中也有着相对不错的效果。

写在最后

① 在实践中，KMeans算法相对其它聚类算法有着相对较快的速度，但仍然需要不小计算资源。

② KMeans算法在拟合过程中容易陷入局部最小值，因此多次重启可能会得到更好的效果。

③ 真实标签已知时一般我们更倾向于使用有监督的分类算法，因此V-Measure在聚类任务中实际并不常用，大部分情况都使用轮廓系数和卡林斯基-哈拉巴斯指数，但使用这两个评价指标在例如DBSCAN以及凸簇上有着部分虚高的分数，因此一般不用于算法之间的效果比较，或进行比较的时候要考虑这部分误差。

④ 对于现在的KMeans算法，一般使用KMeans++提供的方法进行质心初始化处理，使用elkan KMeans算法进行样本划分处理，Mini-Batch KMeans在数据量较多时也常被用到，但WKMeans使用地相对较少。

Reference

• Elkan, Charles. “Elkan’s k-Means Algorithm for Graphs.” Journal of Machine Learning Research, vol. 4, no. 3, 2003, pp. 213-224.
• Andrew Rosenberg and Julia Hirschberg, 2007. V-Measure: A conditional entropy-based external cluster evaluation measure
• Rosenberg, Andrew and Hirschberg, Julia. “V-Measure: A Conditional Entropy-Based External Cluster Evaluation Measure.” Journal of Machine Learning Research, vol. 8, 2007, pp. 283-323.
• 聚类性能评估-V-Measure – 知乎 (zhihu.com)
• 兰德系数、调整兰德系数-CSDN博客
• 信息量与信息熵_信息量和信息熵-CSDN博客
• 菜菜的sklearn机器学习
• 详解Kmeans两大优化——mini-batch和Kmeans++ – 知乎 (zhihu.com)
• elkan K-Means-CSDN博客
• 【聚类算法】MiniBatchKMeans算法-CSDN博客