简明聚类分析入门

摘要 : 以“为什么需要聚类分析这一问题”作为引入，逐步阐述聚类分析领域是如何发展的。这篇文章主要阐述聚类分析的四类方法：划分方法、层次方法、基于密度的方法和基于网格的方法的基本原理以及它们中的代表算法和实现方式。将聚类算法的设计总结为两大核心：划分过程和相似度量的设计。

1 引言

聚类分析技术已经发展了近60年，至今该领域依旧非常活跃 $^{[1]}$ 。聚类分析的地位与其他的机器学习理论，如分类，SVM等，有所不同。首先，聚类分析是一个多学科交织的研究领域，大部分的聚类算法都需要跨域知识(领域知识)的参与，不同角度会产生不一样的聚类结果，因此聚类结果没有统一的衡量标准；其次，聚类分析领域中的算法多如牛毛，但很难对这些算法提出一个总的划分(概念的划分)，算法之间存在着重叠概念。以上原因导致了聚类分析理论难以形成系统化的理论体系，以至于有些教科书只会轻描淡写的阐述聚类分析的一些理论。但是不可否认的一点是，聚类分析在现实任务中占据着非常重要的地位。

本文从聚类分析本质的角度阐述聚类分析的来龙去脉，让读者们能够简单明了的掌握聚类分析的基本概念。在第二部分将会围绕着4个聚类分析的本质问题，即：

为什么需要聚类分析？
什么是聚类分析？
如何进行聚类分析？
如何评估聚类分析结果？

展开详细的论述；第三部分阐述聚类分析中一些具有代表性算法以及常用的聚类结果评估标准；希望通过本文能让读者对聚类分析有一个全方位的理解。

2 聚类分析的因与果

不同的研究者对聚类分析这个领域的观察角度是不同的，例如周志华教授的《机器学习》 $^{[2]}$ 就会从性能度量和距离度量着手，而韩家炜教授 $^{[3]}$ 则会以数据挖掘的一个工具的角度讲聚类定义和要求。尽管不论从何种角度看待聚类最终我们学到的东西都是一致的，但对于入门者来说这些就有点摸不着头脑了。就比如我当初学习聚类时我更希望教科书能告诉我为什么需要聚类技术？，而后才是如何定义聚类这个概念等等。因此本节从聚类最本质的问题着手，回答聚类分析的因与果。

2.1 聚类分析的因

问题一：为什么需要聚类分析？

在我们接触到的许多常见的机器学习任务(算法)大部分都是带有标记的(有监督的)，如分类与回归。分类与回归问题恐怕是绝大部分学习机器学习理论的人最开始接触到的内容。有监督学习(或者称有先验结果的学习)自然是从业者最愿意看到的一种学习任务，然而在真实的世界中，有人工标记的数据仅仅占很少的一部分，对于大数据时代而言甚至可以忽略不计（可以说无监督(不带标记)的学习才是机器学习的终极目标）。但我们仍然希望探究那些无法或者还未标记的数据的内在数据结构或规律。 聚类分析就此诞生。

问题二：什么是聚类分析？

从问题一中我们知道了聚类分析的目的，就是要寻找非预先得知的数据分类信息，其中分类信息值得就是我们现在所称呼的“簇”的概念，物以类聚。韩家炜教授 $^{[2]}$ 用了更为精炼的语言阐述这一概念，即：

Cluster analysis or simply clustering is the process of partitioning a set of data objects (or observations) into subsets. Each subset is a cluster, such that objects in a cluster are similar to one another, yet dissimilar to objects in other clusters. The set of clusters resulting from a cluster analysis can be referred to as a clustering.

聚类分析，是一个把数据对象(或观测)划分成子集的过程。每个子集是一个簇，使得簇中的对象彼此相似，但与其他簇中的对象不相似。由这个过程产生的一系列簇成为一个聚类。

2.1 聚类分析的果

问题三：如何进行聚类分析？

从问题二聚类分析的定义定义中我们可以提取两个关键的信息点，划分过程与(对象之间)相似与否。实际上聚类分析所有的工作都关注在这两者上，划分过程一般的有：

依距离划分
依层次划分
依密度估计划分
依概率估计划分
依网格划分

不同的划分方式有不同的优缺点（如簇的形状，数据集的大小等），没有包打一切的算法 $^{[2]}$ 。而对于相似度（一般是特征空间上的相似度，因此也称距离度量），一般有：

有序属性的相似度，常用的欧氏距离(Euclidean distance)，闵氏距离(Minkovski distance)等。
无序属性的相似度，如VDM(Value Difference Metric)。

具体可查阅周志华教授的机器学习 $^{[2]}$ ，这里就不再赘述。

问题四：如何评估聚类分析结果？

在问题三我们知道聚类分析主要关注划分方式与相似度量，自然地评价一个聚类结果的好与坏也主要关注在这两者上，即评估划分结果与聚类的质量。对于估算簇的个数，一般地常用以下两种：

经验法，簇数取 $\sqrt{\frac{n}{2}}$ ，其中 $n$ 为样本数。
肘方法(elbow method)，其中常用的为Gap Statistic $^{[4]}$ 的方法，在第三部分将会详细讲解。

对于评估聚类质量，一般地也有两种方法：

有监督方法(外在方法)，即我们有对聚类结果的期望结果(真实的聚类结果)，然后对比期望结果与聚类器的结果的差异，例如BCubed $^{[5]}$ 算法，该算法在scikit-learn库中有相应的实现。
无监督方法(内在方法)，此时我们并没有对聚类结果的期望结果，因此我们仅能从聚类的本质出发评估聚类结果，即簇内距离是否足够近且簇间距离是否足够远，其中轮廓系数(silhouette coefficient)较为常用，同样在被实现于scikit-learn中。

但从上文中我们了解一个事实，即使聚类技术发展到如今也没有产生一个绝对统一的衡量标准，因此所有的评估方式都只能充当一种参考，一种从数据角度(统计理论)出发的参考。

至此，我们已经理解了聚类分析的由来以及功用。接下来我们将对问题三与问题四展开详细的讨论，但由于聚类分析这棵大树的果实非常多，因此我们仅能阐述有代表性的部分，关于它们的变种或改进我们也会稍微提及。

3 聚类算法与评估方法

至今为止，都没有一个很准确的方法去归类聚类算法，有的研究者将聚类算法大致归并为两类 $^{[1]}$ ：划分方法(Partition methods)和层次方法(Hierarchical methods)，有的会归并为四类 $^{[3]}$ ：划分方法、层次方法、基于密度的方法(Density-based methods)和基于网格的方法(Grid-based methods)，而在最近出现了一种更为细致的 $^{[6]}$ 划分方式，如图3.1所示：

图3.1: 聚类算法归类方式，取自Rehioui, Idrissi, Abourezq和Zegrari[6]

新增了基于模型的方法(Model-based methods/clustering)的划分方式，但纵使有这些总结性的划分方式，也不见得所代表的算法就是“纯”的(即特性完全属于该分类)，例如CLIQUE(CLustering In QUEst) $^{[21]}$ 就是同时基于网格和密度的算法，而且算法间也经常互相借鉴思想。本文将会以图1中的归并方式为下文组织结构阐述若干个有代表性的算法(其中基于模型的方法属于高级内容，不在本文的进行阐述)，本文所有的代码均可以在附录中找到。

3.1 划分方法

划分方法，是聚类分析中两大经典的划分过程之一。根据距离度量将 $n$ 个观测对象划分成 $k$ 个簇(分区)，其中 $k \leq n$ 。显然当 $k=n$ 时是划分方法的一种极端情况，此时代价函数最小(为0)。其中以KMeans为代表的算法最为常用，由于划分方法基于距离度量，因此非常擅长发现球状的互斥簇，并且能发现每个簇的代表对象，也常用于压缩。但由于划分方法在计算时通常需要计算每一对对象之间的距离，因此通常在小数据集上表现得很好，但近几年不断有研究者创造出了高伸缩的适合大规模数据集的划分方法，如KMeans的变种MiniBatchKMeans $^{[7]}$ 。本节主要介绍KMeans以及另一种围绕中心点划分算法PAM(Partitioning Around Medoids)。

3.1.1 形心划分方法KMeans

KMeans以簇内变差度量为划分依据，即：

E = \sum^{k}_{i=1} \sum_{p \in C_i} dist(p, c_i)^2 \tag*{(3.1)}

将 $n$ 个观测对象划分成以 $k$ 个中心点(均值点)为代表的簇 $C_1, C_2 , \cdots , C_k$ ，其中 $k \geq 2$ 。KMeans算法通常有两个步骤组成：分配(Assigment)和更新(Update)，一开始 (1) 算法会随机在数据集 $D$ 中选择 $k$ 个中心点，然后 (2) 将 $n$ 个观测对象依据距离(一般为欧氏距离)最近原则分配给 $k$ 个簇，最后 (3) 以启发式的方式迭代计算(式3.1)得到新的 $k$ 个中心点，重复步骤(2)(3)直到中心点不再发生改变或者已达到设定的迭代次数，关于KMeans的详细内容可以查阅《KMeans原理、实现及分析》。

KMeans因此简便易操作的特性，使其逐渐成为了聚类分析任务中的首选聚类器。但KMeans也存在许多缺陷，如无法发现任意形状的簇，在图3.2中，由于KMeans基于欧氏距离导致其不能发现类似同心圆簇或是U型簇。不仅如此，KMeans使用了 平均值(Mean) 这种较为脆弱的统计量还导致了其对离群点(噪点)敏感。而且由于在步骤(1)使用了随机初始化中心点的策略，无法保证KMeans算法能够收敛到全局最优解，因此需要多次运行KMeans，才能的到一个较为满意的聚类结果。

图3.2: KMeans在任意形状簇上的表现

到目前为止，已有许多KMeans的变种算法被开发出来，用于处理不同需求的聚类任务。当要处理的数据是类别属性(标称属性)时，使用平均值可能是无意义的，这时就可以使用众数(mode)来代替均值，即K-mode算法 $^{[8]}$ ，在Python上可以使用实现了K-mode的第三方库；当需要处理较大规模的数据集时，KMeans就显得抓襟见肘，因此一种行之有效的方法就是选取适当规模的数据集进行多次训练，如MiniBatchKMeans，该算法被scikit-learn作为标准工具实现了；当数据存在噪声或离群数据时，KMeans算法产生的聚类结果与真实结果存在较大的偏差，一种抗噪方案是使用中心划分的技术，使用真实观测点代替均值点，如3.1.2节提到的K-medoids算法。

即使KMeans本身已经非常成熟，但仍然有许多新的变种相继被开发出来，甚至可以KMeans的历史是聚类分析发展史的缩影。更多的KMeans相关信息可以查阅Jain $^{[1]}$ 的文章。

3.1.2 中心划分方法PAM

为了减少离群点对聚类结果的影响，Kaufman和Rousseeuw $^{[9]}$ 于1987年提出了围绕中心划分的方法，一种流行的实现是K-medoids，如图3.3。K-medoids算法与KMeans一样可以分为两大步骤：分配与更新，在开始时 (1) 会随机地在数据集 $D$ 中选择 $k$ 个中心点 $o$ ，然后 (2) 将剩余的所有观测对象依据最近原则分配给每个中心点，形成 $k$ 个簇，最后 (3) 随机地选取 $k$ 个非中心对象 $o_{random}$ 尝试代替最原来的中心对象 $o$ ，计算 $o_{random}$ 的平均距离代价差 $S$ ，若 $S<0$ ，则用 $o_{random}$ 替换原中心对象 $o$ ，否则继续寻找，直到遍历完所有对象为止。算法迭代地执行步骤(2)(3)直到中心对象 $o$ 不再发生改变。该算法的代码实现可在附录文件中找到。

图3.3: K-medoids算法，取自Han J[3]

显然，步骤(3)的计算是非常消耗时间的，当 $n,k$ 都比较大时，其时间计算开销远大于KMeans，虽然K-medoids算法能够有较好的抗噪能力，但很难在大型数据集上有较好的表现。因此Kaufman和Rousseeuw于1990年又提出一种基于抽样的中心划分方法CLARA(Clustering LARge Applications) $^{[10]}$ 拟补K-medoids算法处理大型数据的不足，理论上只要抽样数据分布与原数据分布足够接近，是可以找到最佳中心点的，但这就会导致CLARA的聚类结果的好坏依赖抽样的好坏。于是乎Raymond和Jiawei Han提出了更为高效的CLARANS(Clustering Large Applications based on RANdomized Search) $^{[11]}$ 算法，基于随机搜索查找全局最优解，该算法在Python上被pyclustering库实现。此外Schubert和Rousseeuw ${^{[22]}}$ 在2019年提出了更快速的K-Medoids算法，有兴趣的读者可以自行查阅！

3.1.3 小结

划分方法的算法步骤大多都是相似的：随机的选择初始点，然后划分簇，不断尝试将一个簇中对象放至另一个簇中来降低代价函数。划分方法通常的划分依据是距离，常用的欧氏距离、马氏距离和闵氏距离等，通过更换距离度量有可能会提高聚类质量。划分方法通常在高斯簇上表现得很好，但不能发现任意形状的簇。

3.2 层次方法

与划分方法要求聚类结果为互斥簇不同，层次方法更加强调对象间层次关系。如在一个公司的管理体系内，可以分成董事局-总经理-部门经理-业务员，利用“树状”的图形表示这些层次关系，是层次聚类可视化的常用方法。一般地，使用两种划分策略来完成层次聚类：凝聚(Agglomerative) 和 分裂(Divisive)，用数据结构的语言来描述，即自底向上和自顶向下。这两者不仅含义相仿，甚至层次聚类里面的许多算法都与数据结构中的构建各种“树”的算法相似。而层次聚类常用的距离度量是有：

\begin{aligned}
& 最小距离: dist_{min}(C_i, C_j) = \min_{p \in C_i, p^{'} \in C_j} \{ | p - p^{'}| \} \qquad & (3.2) \\
& 最大距离: dist_{max}(C_i, C_j) = \max_{p \in C_i, p^{'} \in C_j} \{ | p - p^{'}| \} \qquad & (3.3) \\
& 均值距离: dist_{mean}(C_i, C_j) = | m_i - m_j | \qquad & (3.4) \\
& 平均距离: dist_{avg}(C_i, C_j) = \frac{1}{n_i n_j} \sum_{p \in C_i, p^{'} \in C_j} | p - p^{'}| \qquad & (3.5)\\
\end{aligned}

其中如果算法是以最小距离作为划分依据，那么此时算法就会在不同簇之间添加一条线段作为边，形成最大连通图，也称单链接算法(Single-linkage algorithm)，形象的理解为就是只有一条连通分量(线)链接所有的簇，如图3.4(b)中的情况所示，数据集被一条线全部链接。而如果是以最大距离作为划分依据，则此时形成的是最大完全图，也称全链接算法(Complete-linkage algorithm)，形象的理解为是由多个完全连通的子图构成的簇，如图3.5(c)所示。本文将在3.2.1阐述凝聚与分裂的层次聚类的基本原理，在3.2.2节探讨其他一些流行的大型层次聚类算法。

图3.4: 单链接与全链接，取自Han J[3]

3.2.1 凝聚与分裂的层次聚类

凝聚层次聚类采取的是自底向上合并簇的方式，最终将n个节点合并成1个簇(树根root)，类似最小生成树算法。只需要 $n$ 次迭代即可完成聚类，但弊端在于无法更正以前的操作。凝聚聚类算法也是静态层次聚类算法中最常用的，在scikit-learn等机器学习库中实现了相应的算法类AgglomerativeClustering，周志华教授用了西瓜4.0数据 $^{[2]}$ 来展示基于全链接的凝聚层次聚类，我们也可以尝试使用其数据来复现结果，使用SciPy提供的dendrogram工具能够使用树状图展示聚类结果，如图3.5所示(除顺序不同外，分类结果一致)，横轴对应样本编号，纵轴表示样本距离。

图3.5: 西瓜数据4.0 全链接聚类算法

层次聚类的一个好处是可以根据需要筛选簇，如图3.5中，若簇的粒度为4，则可以得到：

\begin{aligned}
C_1 & = \{1, 2, 3, 4, 21, 22, 26, 29\}, \\
C_2 &= \{23, 24, 25, 27, 28, 30\}, \\
C_3 &= \{5, 7, 9, 13, 14, 16, 17\}, \\
C_4 &= \{6, 8, 10, 11, 12, 15, 18, 19, 20 \}。
\end{aligned}

而分裂的层次聚类要面对的挑战就比凝聚的方式要多，最大的一个挑战是如何将大簇变成一个小簇，拆分成互斥簇的方法是多种多样的，当样本数 $n$ 足够大时，问题是NP-Hard的。因此只能牺牲精度采用启发式的算法，这种策略通常无法修改以前的操作，因此分裂的算法在工业界并不常用，在学术界受到的关注也比凝聚型的要少，目前较新的研究成果是Tengke Xiong $^{[12]}$ 等人的DHCC(Divisive Hierarchical Clustering of Categorical data)算法用于处理标称属性数据，同时该算法也具有很好的伸缩性，适用于大型数据集。

3.2.2 现代流行的层次聚类

经典的层次聚类虽然拥有非常好的可视化潜力，但在面对大型数据时就会显得很无力，且每一层一旦划分结束后便不可修改，很有可能会错失最优解。为了针对层次聚类不同方面的局限性，现代的层次聚类技术都采用了复合聚类的策略，即与其他聚类方式结合：划分方法，如BIRCH(Balanced Iterative Reducing and Clustering using Hierarchies) $^{[13]}$ ；或密度方法，如Chameleon算法 $^{[14]}$ 等。

BIRCH算法(译为：利用层次结构的平衡迭代归约和聚类)，该算法使用聚类特征树(Clustering Feature Tree)来描述一个簇，聚类特征(CF)是一个三维向量：

CF=\langle n, LS, SS\rangle \tag*{(3.6)}

其中 $n$ 是一个簇 $C$ 的样本数， $LS$ (Linear Sum)是n个样本的线性和， $SS$ (Square Sum)是n个样本的平方和。使用CF能够很轻易地推导出簇的形心 $x_0$ ，半径 $R_0$ 和直径 $D_0$ 这些统计量：

\begin{aligned}
x_0 & = \frac{\sum^{n}_{i=1}x_i}{n} = \frac{LS}{n} \qquad \qquad &(3.7) \\
R_0 & = \sqrt{\frac{\sum^{n}_{i=1}(x_i - x_0)^2}{n}} & (3.8) \\
D_0 & = \sqrt{\frac{\sum^{n}_{i=1} \sum^{n}_{j=1} (x_i - x_j)^2}{n (n-1)}} = \sqrt{\frac{2nSS - 2LS^2}{n(n-1)}} &(3.9)
\end{aligned}

其中半径 $R_0$ 和直径 $D_0$ 都反映了簇 $C$ 的紧凑程度，且对于任意两个不相交的簇来说，它们的CF是可加的。如 $C_1, C_2$ 是两个不相交的簇，对应的CF为 $CF_1=\langle n_1, LS_1, SS_1 \rangle, CF_2=\langle n_2, LS_2, SS_2 \rangle$ ，合并 $C_1, C_2$ 成为 $C_3$ ，则 $C_3$ 的CF为： $CF_3 = \langle n_1 + n_2, LS_1 + LS_2, SS_1 + SS_2 \rangle$ 。该算法集成了层次聚类与划分聚类的思想，能够适应大规模数据的训练，拟补了凝聚方法的不足。CF树(CF-Tree)是一个高度平衡的树，CF树有两个重要的参数：分支因子(Branch factor) $B$ 和阈值(Threshold) $T$ ，其中 $B$ 定义的是非叶结点的子女的最大数目，即最大分支数，而 $T$ 定义的是叶节点存储的最大直径 $R$ 。两个参数超过任意一个都会产生分裂，类似B+树。BIRCH是一种动态的层次聚类方法，从CF的特征上看该方法能够很好的发现圆形的簇，且支持增量聚类，因此常用于处理流数据聚类任务。scikit-learn也实现了BIRCH算法。

Chameleon同样也是动态层次聚类算法，结合了密度方法的特点，引入了簇间互连度 $RI(C_i, C_j)$ 和簇间接近度 $RC(C_i, C_j)$ ，这两个度量一同描述了 $C_i, C_j$ 的相似程度。该方法由于结合了密度方法，因此能够发现任意形状的簇，同时在高维数据下比DBSCAN方法的表现要更好。需要深入了解的读者可以查阅Karypis等人 $^{[14]}$ 的文章。

此外还有一些基于概率估计的层次聚类方法，如贝叶斯层次聚类 $^{[15]}$ ，使用概率估计的层次聚类不需要寻找合适的连接度量，同时能够很好的处理数据缺失的问题，且基于概率估计的优化目标明确，能够直接求解全局最优解，而不会像前面提到的方法一样陷入局部陷阱。但这类算法的弊端是当数据量比较大时，参数估计的计算开销是巨大的。

3.2.3 小结

静态的层次聚类有凝聚和分裂两种策略进行聚类，它们都是不可修改的，其中凝聚类的算法比较多，效果也比较好。动态层次聚类如BIRCH等能够更好的适应大型数据集，且由于结合了其他聚类方法，因此性能表现要优于静态层次聚类。

3.3 基于密度的方法

前面我们阐述了划分方法和层次方法的基本原理，但它们都有一个共同的特点是：都是基于距离度量，因此它们仅能处理圆形簇的情况，类似于图3.6中展示的数据形状它们就很难处理好。如果数据是充满噪声的，那么它们就表现大打折扣。

为了克服它们无法发现任意簇的缺陷，研究者们重新考虑了簇的划分过程，一个事实是在一个足够大的数据空间中，噪点总是稀疏的，而真正的簇总是稠密的，因此就有了基于密度的划分方法。有了划分策略，那么现在的问题就是如何定义相似度呢？换句话来说到底如何才算是稠密的呢？这一个问题恐怕到目前为止也很难有统一的标准，就如同为划分方法和层次方法寻在一个绝对统一的距离度量一样难。如著名的DBSCAN(Density-Based Spatial Clustering of Application with Noise) $^{[16]}$ 算法就是需要用户手动调节密度参数；而为了克服DBSCAN的需要设置全局参数的弊端，Ankerst等人 $^{[17]}$ 发明了OPTICS(Ordering Points To Identify the Clustering Structure)算法，通过簇排序的方式来辨识稠密的簇，但其仍然要手工设置广泛密度参数；为了完全摆脱参数的限制，Hinneburg和Keim $^{[18]}$ 提出了基于密度函数估计的算法DENCLUE(DENsity-based CLUstEring)，该算法使用了高斯核函数作为相似度量，使用者不需要考虑经验参数的问题。本节在3.3.1详细介绍DBSCAN的基本原理以及在Python上的使用方法，在3.3.2介绍OPTICS算法，在3.3.3阐述DENCLUE的原理以及其变种算法的发展。

图3.6: 任意形状的簇

3.3.1 DBSCAN算法

DBSCAN的划分策略是：根据指定的邻域半径 $\varepsilon$ 以及邻域内至少需要包含的对象个数 $MinPts$ 来确定一块稠密区域，密度相连(density-connected)的最大稠密区域构成一个簇（其算法原理类似寻找最大连通块的问题）。从该算法的划分策略上，DBSCAN是可以区分噪声的。在DBSCAN算法中有3个非常重要的概念：

概念1: 核心对象(Core object) 核心对象指的是该对象 $o$ 在其邻域半径 $\varepsilon$ 内存在不少于 $MinPts$ 个对象。核心对象是一个稠密区域的支撑对象，它的作用在于不断的发现密度连通的对象。一个簇可能拥有许多个核心对象。
概念2: 直接密度可达(directly density-reachable) 对于一个核心对象 $o$ 和一个对象 $p$ ，如果 $p$ 在 $o$ 的 $\varepsilon$ -领域内，那么 $p$ 是可以从 $o$ 直接密度可达的。此时如果 $p$ 也是一个核心对象，则也可以说 $o$ 是从 $p$ 直接密度可达的，如果 $p$ 不是核心对象，则错误。直接密度可达的作用在于将可以从核心对象直接密度可达的对象全部纳入考虑范围内，发现稠密区域。
概念3: 密度相连对于一个核心对象 $o$ 以及两个对象 $p_1, p_2$ ，且 $p_1 \in A_1, p_2 \in A_2$ ，此时如果 $p_1, p_2$ 都是可以从 $o$ 直接密度可达的，那么 $p_1, p_2$ 就是密度相连的。密度相连与密度可达是有区别的，密度相连是等价关系，而密度可达除非两者都是核心对象否则不是等价的。密度相连的作用在于将若干个密度可达的稠密区域连接成一个大的稠密区域，即形成簇。

DBSCAN算法在开始时 (1) 将数据集 $D$ 中所有的对象都标记为 "unvisited"，然后 (2) 随机选择 $D$ 中一个 "unvisited" 的对象 $p$ ，将其标记为 "visited"，接着 (3) 检查 $p$ 的 $\varepsilon$ -邻域内是否包含 $MinPts$ 个对象，如果是则为 $p$ 创建一个新的簇，并且将该邻域内的所有对象加入待观测队列中，如果不是，则标记 $p$ 为噪声。最后 (4) 继续从待观测队列或者数据集 $D$ 中取出对象 $p$ 执行步骤(3)，直到队列为空且 $D$ 中没有 "unvisited" 的对象为止。伪代码如图3.7所示。

图3.7: DBSCAN伪代码，取自Han J[3]

在sckit-learn中有相应的DBSCAN算法实现，只需设定邻域半径 $\varepsilon$ 以及邻域阈值 $MinPts$ 即可，图3.8展示了DBSCAN在任意形状簇上的聚类能力。从结果上来看，DBSCAN与预期聚类结果几乎一致。

图3.8: DBSCAN处理任意形状簇

但是正如前面所说，DBSCAN的聚类能力取决于邻域半径 $\varepsilon$ 以及邻域阈值 $MinPts$ 的设定，而这些往往都是经验使然的结果，因此使用DBSCAN算法需要大量的领域知识注入，才能取得一个较为满意的聚类结果。

3.3.2 OPTICS算法

OPTICS算法，是为了克服DBSCAN的全局参数对于描述数据能力不够而提出来的(因为在高维的数据中，全的参数不一定能很好的描绘簇内结构)。但是它并非完全不使用参数，而是使用一种广泛参数策略，即通过簇排序来寻找合适的簇结构，从而达到聚类的目的。其原理是：较稠密的对象在簇排序中的位置是相对接近的。它与DBSCAN的算法结构大致相同，区别在于它重新定义了簇的相似依据：

概念1: 对象p的核心距离(core-distance of an object p)

对于任何一个对象 $p$ 都会寻找其最小的半径阈值 $\varepsilon^{'}$ ， $\varepsilon^{'}$ 称为 $p$ 的核心距离，如果 $p$ 不是核心对象(见3.3.1)，则 $p$ 的核心距离是未定义的。这个核心距离的作用在于尽最大可能去使得 $p$ 成为一个核心对象，即自动的发现簇，显然如果 $p$ 是噪点，那么 $p$ 一定会是为定义的。对于不同的核心对象 $p$ 其 $\varepsilon^{'}$ 值很可能会不同，这是与DBSCAN的区别之一。
概念2: 可达距离(reachability-distance)

对于对象 $p, o$ ，如果 $o$ 是核心对象，则 $p$ 是从 $o$ 密度可达的(这点与DBSCAN一致)，此时为 $p$ 设置可达距离为 $max\{core-distance(o), dist(p, o)\}$ ，如果 $o$ 不是核心对象，则 $o$ 到 $p$ 无定义。显然， $p$ 很有可能同时拥有对个核心对象的可达距离，OPTICS只关注最短可达距离，这样就能发现最适合的簇结构。可达距离的作用在于对其排序后，能使得距离接近的簇位置相对较近，也就是得出了簇的结构。

另外的一个区别是使用了一个OrderSeeds表的数据结构，这个表维护了各个对象p到其核心对象的最短可达距离。OPTICS还有另一个有趣的用途，就是用于簇结构可视化。对于高维数据来说，我们很难观察其全貌，但OPTICS根据可达距离对所有数据进行排序后将其展示就能看到簇的结构，图3.9中有4个簇 $C_1, C_2, C_3, C_4$ ,使用可达距离(纵轴)与排序结果(横轴)展示簇，可以看出有四块空白区域(凸区域)，这能反映簇的真实情况。

图3.9: OPTICS用于数据可视化取自Ankerst, Breunig和Kriegel[16]，有修改

3.3.3 DENCLUE算法基本思想

无论是DBSCAN或OPTICS都是对半径值 $\varepsilon$ 敏感的，过大的 $\varepsilon$ 会导致簇内不够紧凑，而过小的 $\varepsilon$ 会导致不完全聚类。为了解决这一问题，Hinneburg等人 $^{[18]}$ 引入了统计学中非参数估计方法——核密度估计，DENCLUE使用高斯核估计待聚类的对象集密度。其局部最大点(即一个解)称为密度吸引点(density-attractor) $x^*$ ，每对密度吸引点之间定义了一条不大于噪声阈值 $\xi$ 的连通路径，通过这样的密度划分方式发现任意形状的簇。由于篇幅限制，不在这里详细展开DENCLUE算法的讨论。

DENCLUE发展至今已经衍生出许多变种算法，如Rehioui等人 $^{[19]}$ 提出的DENCLUE-IM算法，它能够应对更加复杂的数据种类，更多的数据量的聚类任务，感兴趣的读者可以自行研读。

3.3.4 小结

基于密度的方法能够发现任意形状的簇，但是难免会先于经验困境，即依赖经验参数的设定。但即便如此DBSCAN仍然是任意形状簇聚类任务的首选，因为其过程更可控，而OPTICS常常用于发现簇内结构。

3.4 基于网格的方法

如果从划分的角度去区分的话，那么前面所提到的所有方法都可以算是以数据驱动的划分方式 $^{[3]}$ ，即关注点在数据本身，从数据上划分出一系列空间。而基于网格的方法则属于空间驱动的划分方式，即从输入空间的角度考虑，将输入空间划分成一系列子空间，再把数据嵌入其中。其划分方法是使用一种多分辨率的网络数据结构，将对象空间量化成有限数目的单元。不同的算法使用不同的相似度量，如果STING(STatistical INformation Grid) $^{[20]}$ 使用统计信息作为度量，进行层次聚类，而CLIQUE则是一种子空间的密度聚类。CLIQUE是网格聚类中常用的算法之一。

3.4.1 STING: 统计信息网格聚类方法

STING利用分层和递归对输入空间进行划分，每一个高层都被划分为多个低层单元，如图3.10中，第 $(i-1)$ 层的单元被划分成第 $i$ 层的四个单元，每个单元格会预先存储一些来自低层单元格的统计信息，如均值，标准差，最大值和最小值等，使用者来可以预先指定数据的分布情况。最底层可以近似的等于DBSCAN的聚类结果，在产生查询请求时，STING会依照自顶向下的方式依次查询相关的单元格，直到最底层为止。

图3.10: STING层次结构取自Wang, Yang和Muntz[20]

这种基于网格查询层次聚类的特点使得STING的执行效率非常高，通常在线性时间内就能完成一次聚类信息的查询，不仅如此，STING还支持并行处理与增量更新操作。但其弊端也是很显然易见的，就是依赖于最底层的粒度，如果最底层的粒度很细，则STING模型的层数就会很多，因此处理的代价就会增加，如果最底层的粒度很粗，那么聚类的质量就会大打折扣。这是个两难的问题，因此导致STING的使用场景并不是很多。

3.4.2 CLIQUE: 子空间密度聚类方法

划分子空间是属性约简的一种常见手段，而随着信息的爆炸式增长需要处理的数据维度越来越高，高维度的信息通常含有许多冗余(不相关)的属性，因此这对在整个空间中寻找合适的簇是一件很难完成的任务，CLIQUE算法基于划分子空间的策略，对若干个子空间进行密度聚类就可以发现稠密的区域，即合适的簇。

CLIQUE算法有两个运行阶段：

第一阶段(分割与链接): CLIQUE将 $d-$ 维数据空间按照每个维度划分成若干个互不重叠的子空间，然后识别密度阈值超过 $l$ 的区域，之后迭代地链接各个子空间中重叠的区域。
第二阶段(合成): CLIQUE将第一阶段得到的若干个稠密单元合成为一个最大的稠密区域。使用贪心算法，从任意稠密单元开始找出覆盖这些稠密单元的最大区域，直到这一区域所有的稠密单元都已被覆盖完。循环这一步骤直到所有数据都已被覆盖完全。

虽然CLIQUE能够自动地发现含有高密度簇的子空间，但是它的聚类结果十分依赖密度阈值 $l$ 的设定，如图3.11与图3.12分别为 $l=10$ 和 $l=20$ 的聚类情况，预期的聚类结果应该是2个同心圆的簇，而当 $l=10$ 时，在图3.11(a)中可见满足这一阈值的稠密区域是非常多的，因此两个簇之间就非常相似，导致最后聚类结果仅有一个簇，如图3.11(b)所示。当重新调整 $l=20$ 后，在图3.12(a)中可以发现簇间距离增大，CLIQUE能够发现两个最大的稠密区域，最有将它们合成为两个互斥的簇，如图3.12(b)所示。

图3.11: CLIQUE聚类结果，网格密度因子为10

图3.12: CLIQUE聚类结果，网格密度因子为20

合适密度阈值 $l$ 是要花费大量实践时间去得到的，因此这也成为了CLIQUE最大的弊端。

3.5 评估聚类质量

至此，我们已经了解了聚类分析领域几个有代表性的算法，聚类的总体步骤总是相似的，定义划分过程和相似度量。本文所使用的许多数据集实际上都是带有标记的数据，但对于无标记的数据，我们该如何判断这些聚类结果的好坏呢？我们从两个方面着手：(1) 划分的簇的个数是否合适？(2) 聚类质量是否够好？本节我们将讨论这两个问题。

3.5.1 确定簇的数量

如K-Means等算法需要用户设定预期的簇数，如网格方法需要用户定义粒度等，这些都要求我们对最后的簇数有正确的评估，否则聚类结果没有实际意义。通过前面的介绍我们知道簇的数量是有两个极端的：当簇数为1或为 $n$ 时。显然这两者都不可取，但如何在 $1～n$ 之间选定一个合适的数字呢？假如我们使用(式3.1)来衡量聚类的误差和，如图3.13所示，在(a)中我们判断实质上簇数为2是最好的选择，而在(b)中簇数越大，簇内变差就越小，这也解释了为什么当簇为 $n$ 时，误差最小。

图3.13: 簇内变差曲线，取自Tibshirani, Walther和Hastie[22]

那么显然误差不能作为衡量簇数的标准，常用的方法有二：(1) 经验法和 (2) 肘方法(elbow method)，经验法一般是取簇数为 $\sqrt{\frac{n}{2}}$ ，显然不适合所有场景，而肘方法在大多数情况下是适用的，如图3.13(b)肘显然出现在2的位置，但是当如果数据分布不近似高斯分布的时候呢？这时候肘的位置是不明显的，因为误差曲线会很平滑。因此本节要介绍的是Tibshirani, Walther和Hastie提出的Gap Statistic $^{[4]}$ 的方法，该方法使用分布相似来确定最合适的簇数。要想理解这一概念首先我们需要知道一个事实是：如果一个数据集的分布是均匀分布的，那么这个数据集的聚类结果没有任何实际意义，因为均匀分布产生的簇是随机的。Gap Statistic的思想是假设测试在每一个待定的簇数时，都对真实数据的分布与参考分布(reference distribution, 实质上是一个均匀分布)进行对比，随着簇数增加时，真实数据分布的误差函数的值与参考分布的误差函数的值一定是在同时下降的，而在肘方法我们得知在最合适的簇数出现的时候会出现一个拐点，但对于均匀分布对这一簇数并不敏感，因此在此时两个分布的误差函数值就会出现一个全局最大的空隙(Gap)，使得最大间隙出现的簇数就是最优解。因此我们就有了间隙函数(Gap)定义：

\tag*{(3.10)} Gap(k) = \frac{1}{B}\sum_{n}log(W^{*}_{kb}) - log(W_k)

其中 $k$ 是簇的数量， $B$ 是参考分布的个数， $W_k$ 是在簇数 $k$ 下的误差值， $W^{*}_{kb}$ 是在簇数 $k$ 下第 $b$ 个参考分布的误差值，这一理论的证明可以查阅Tibshirani等人的论文。而Gap Statistic的目标就是使得 $Gap(\cdot)$ 最大。具体代码实现见本文附录文件。

接下来我们可以尝试利用Gap Statistic来评估两个高斯簇合成的数据集和周志华老师的西瓜数据4.0。如图3.14对于两个高斯簇合成的数据集，其最大Gap出现在 $K=2$ 的位置，符合我们所预期的结果。在图3.15中使用的西瓜数据4.0，最大的Gap出现在 $K=4$ 的位置，也符合西瓜数据的预期结果 $^{[2]}$ 。

图3.14: 高斯簇与Gap Statistic

图3.15: 西瓜数据4.0与Gap Statistic

3.5.2 评估簇的质量

在得到了合适的簇数量之后，需要考虑的事情是评估不同聚类算法得出的结果。如对于图3.6的任意形状簇来说，假设我们未知数据的真实分布，那么到底如何选择使用KMEANS、DBSCAN和OPTICS进行聚类呢？一般地，有两种指标对簇的质量进行评估，有监督的评估指标和无监督的评估指标。这些在scikit-learn中的metric包中都有相应的工具。

(1) 有监督的评估指标

在我们已知簇的标记的时候，这时聚类的评估与分类的评估较为相似，即预测的对象是否都在相应的簇中。使用下面工具就可以轻松的得到聚类结果的得分：

from sklearn.metrics import v_measure_score

在同心圆数据上我们分别测试KMEANS、DBSCAN以及OPTICS的得分：

-------------------- Circle Datasets --------------------
KMeans spended time: 0.03s, and got score 0.00
DBSCAN spended time: 0.01s, and got score 1.00
OPTICS spended time: 0.90s, and got score 1.00
---------------------------------------------------------

分数结果会在 $[0, 1]$ 区间，得分为 $1$ 时说明聚类质量很高。

(2) 无监督的评估指标

与有监督相比，无监督的评估就较为复杂。我们没有可以参考的对象之后，就仅能从数据本身出发。一种朴素的策略是考虑簇间的分离情况与簇内的紧凑情况。如轮廓系数(silhouette coefficient)，轮廓系数有两个重要的指标 $a(o)$ 和 $b(o)$ ，其中：

\tag*{(3.11)} a(o) = \frac{\sum_{o^{'} \in C_i, o \neq o^{'}} dist(o, o^{'})}{|C_i| - 1}

表示对象o所属的簇的紧凑程度，

\tag*{(3.12)} b(o) = \min_{C_j:1 \leq j \leq k, j \neq i } \lbrace \frac{\sum_{o^{'} \in C_j} dist(o, o^{'})}{|C_j|} \rbrace

表示簇间的分离程度。轮廓系数:

\tag*{(3.13)} s(o) = \frac{b(o) - a(o)}{\max \{ a(o), b(o)\}}

显然当 $o$ 与所属的簇较为紧凑且远离其他簇，即 $b(o) > a(o)$ 时， $s(o)$ 为正数，越接近1说明聚类质量越高。当 $o$ 所属的簇各个对象距离较大且与其他簇相距较近，即 $b(o) < a(o)$ 时时， $s(o)$ 为正数，越小说明聚类结果越差。在scikit-learn中使用相应的工具:

from sklearn.metrics import silhouette_score

但轮廓系数的一个弊端是会偏向于圆形簇。同样的在同心圆数据集上使用轮廓系数对比KMeans、DBSCAN和OPTICS算法得以下结果：

-------------------- Circle Datasets (without real labels)--------------------
KMeans spended time: 0.04s, and got score 0.35
DBSCAN spended time: 0.01s, and got score 0.11
OPTICS spended time: 0.98s, and got score 0.11
------------------------------------------------------------------------------

此时KMeans的分数要高于我们已知的聚类结果较好的DBSCAN和OPTICS。原因就在于轮廓系数使用了基于距离的度量。

4 总结

我们已经了解了聚类的缘由，以及聚类分析中几个典型的方法以及它们所对应的代表算法。聚类分析目前仍然是一个富有活力的研究领域，尤其是如今数据具有3V的特点，聚类分析更加看重伸缩能力、处理不同类型数据的能力、增量能力、并行处理能力和抗噪能力等。但不管聚类分析算法如何变化，其核心思想是不会变的，即两大问题：划分过程和相似度量。

参考文献

[1] Jain, A. K. (2010). Data clustering: 50 years beyond K-means. Pattern Recognition Letters, 31(8), 651–666.

[2] 周志华. 机器学习[M]. Qing hua da xue chu ban she, (2016): 197-224.

[3] Han J, Pei J, Kamber M. Data mining: concepts and techniques[M]. Elsevier, (2011): 442-495.

[4] Tibshirani, R., Walther, G., & Hastie, T. (2001). Estimating the number of data clusters via the gap statistic. In Journal of the Royal Statistical Society: Series B (Vol. 63, Issue Part 2, pp. 411–423).

[5] Amigó, Enrique, et al.: A comparison of Extrinsic Clustering Evaluation Metrics based on Formal Constraints. In: Information Retrieval 12.4 (2009): 461-486.

[6] Rehioui, H., Idrissi, A., Abourezq, M., & Zegrari, F. (2016). DENCLUE-IM: A New Approach for Big Data Clustering. Procedia Computer Science, 83(Ant 2016), 560–567.

[7] Sculley D. Web-scale k-means clustering[C]//Proceedings of the 19th international conference on World wide web. 2010: 1177-1178.

[8] Huang, Z.: Clustering large data sets with mixed numeric and categorical values, Proceedings of the First Pacific Asia Knowledge Discovery and Data Mining Conference, Singapore, pp. 21-34, 1997.

[9] Kaufman, L. and Rousseeuw, P.J. (1987), Clustering by means of Medoids.

[10] Kaufman L, Rousseeuw P J. Clustering large applications (Program CLARA)[J]. Finding groups in data: an introduction to cluster analysis, 1990.

[11] Ng, Raymond T., and Jiawei Han. "CLARANS: A method for clustering objects for spatial data mining." IEEE transactions on knowledge and data engineering 14.5 (2002): 1003-1016.

[12] Xiong, T., Wang, S., Mayers, A. et al. DHCC: Divisive hierarchical clustering of categorical data. Data Min Knowl Disc 24, 103–135 (2012).

[13] Zhang T, Ramakrishnan R, Livny M. BIRCH: an efficient data clustering method for very large databases[J]. ACM Sigmod Record, 1996, 25(2): 103-114.

[14] Karypis G, Han E H, Kumar V. Chameleon: Hierarchical clustering using dynamic modeling[J]. Computer, 1999, 32(8): 68-75.

[15] Katherine A. Heller and Zoubin Ghahramani. 2005. Bayesian hierarchical clustering. In Proceedings of the 22nd international conference on Machine learning (ICML ’05). Association for Computing Machinery, New York, NY, USA, 297–304.

[16] Ester M, Kriegel H P, Sander J, et al. A density-based algorithm for discovering clusters in large spatial databases with noise[C]//Kdd. 1996, 96(34): 226-231.

[17] Ankerst M, Breunig M M, Kriegel H P, et al. OPTICS: ordering points to identify the clustering structure[J]. ACM Sigmod record, 1999, 28(2): 49-60.

[18] Hinneburg, A., & Keim, D. (1998). An Efficient Approach to Clustering in Large Multimedia Databases with Noise. In Proceedings of 4th International Conference in Knowledge Discovery and Data Mining (KDD 98), 58–65.

[19] Rehioui, H., Idrissi, A., Abourezq, M., & Zegrari, F. (2016). DENCLUE-IM: A New Approach for Big Data Clustering. Procedia Computer Science, 83(Ant 2016), 560–567.

[20] Wang, W., Yang, J., & Muntz, R. (1997). STING : A statistical information grid approach to spatial data mining. Proceedings of the 23rd International Conference on Very Large Databases, VLDB 1997, 186–195.

[21] Agrawal, R., Gehrke, J., Gunopulos, D. et al. Automatic Subspace Clustering of High Dimensional Data. Data Min Knowl Disc 11, 5–33 (2005).

[22] Schubert, E., & Rousseeuw, P. J. (2019). Faster k-Medoids Clustering: Improving the PAM, CLARA, and CLARANS Algorithms. Lecture Notes in Computer Science (Including Subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), 11807 LNCS, 171–187.

附录