Science14年的聚类论文——Clustering by fast search and find of density peaks

这是一个比较新的聚类方法（文章中没看见作者对其取名，在这里我姑且称该方法为local density clustering，LDC），在聚类这个古老的主题上似乎最近一些年的突破不大，这篇文章算是很好的了，方法让人很有启发（醍醐灌顶），并且是发表在Science上，受到的关注自然非常大。

本文的核心亮点：1是用比较新颖的方法来确定聚类中心，2是采用距离的local density来进行聚类的划分。在这两点中，常见的Kmeans算法采用的方法是：用每一类的均值作为中点，用距离的最近的点来确定聚类划分。下面对LDC方法进行描述。

首先来讲一下聚类中心，什么样的点才是好的聚类中心呢？作者认为需要满足下面两点：
1、样本点本身的密度大，即它周围点的密度应该小于它；
2、和比它密度更大的点之间的距离应该尽可能的大；

定义每个样本点的（1）局部密度函数，采用cut-off kernel：

$$\begin{equation} \rho_i = \sum_{j\in I_{S}\backslash\{i\}}\chi(d_{ij} - d_c) \end{equation}$$
其中
$$\begin{equation} \chi(a) =\left\{\begin{matrix} 1, &a<0 \\ 0, &a>=0\\end{matrix}\right. \end{equation}$$
$d_c$是截断距离（cutoff distance），是预先设定的一个值，$d_{ij}$表示点i和j之间的距离，可以用任意距离计算公式计算；$I_S$表示数据集中所有样本点的标号集合；很明显，这个局部密度函数是一个离散的函数，它的含义是所有样本点中（不含i点）和i点距离在$d_c$之内的样本点的数量。当然也可以采用连续的函数来定义局部密度，如采用Gaussian Kernel：
$$\begin{equation} \rho_i = \sum_{j\in I_{S}\backslash\{i\}}\exp(-\frac{d_{ij}}{d_c})^2 \end{equation}$$

Ok，不管采用离散的还是连续的，我们已经得到了局部密度；再定义每个样本点的距离$\delta_i$，为了避免和距离计算的混淆，我这里称之为（2）距离偏量：

$$\begin{equation} \delta_i=\left\{\begin{matrix} \min_{j:\rho_j>\rho_i}(d_{ij}) &i不是密度最大的点 \\ \max_{j}(d_{ij}) &i是密度最大的点 \end{matrix}\right. \end{equation}$$
很明显，i点的距离偏量的含义是，所有局部密度大于i点的样本点到i点的距离的最小值；单是对密度最大的那个点单独处理一下，其距离偏量为所有两两点之间距离的最大值（其实这个最大值的数值本身并不重要，只是为了使得这个密度最大点的距离偏量也是最大的，使得它会被选为聚类中心之一）。上面的定义有一个小处定义的不是很周到，就是对于密度相同的点，距离到底是计算还是不计算？因此可以在计算完局部密度以后，可以对所有点的局部密度进行降序排序，这样相同密度的点也会有个先后顺序，而对于i点来说，密度比它大的点就是那些排序在它前面的点，包括了和它密度相同的点[2]。

好，接下来看paper中的这个例子，下面左图是数据样本的分布，其中样本点的标号是根据局部密度降序排序的。右图是将每个点的局部密度和距离偏移绘制出来的结果，可以看到1和10号点脱颖而出，因此聚类的时候把他们作为两类聚类的中心。

这种聚类中心的判断主要还是通过人眼来看的，这一点有一些麻烦，也不是很普适。在文献[1]作者提到，可以用一种综合性的指标$\gamma_i=\rho_i\times\delta_i$来判断，$\gamma_i$大的点被选为中心点，因此对所有的$\gamma_i$降序排序取前K个点作为聚类中心点，具体取多大的K，可以观察下图找到分割点（具体内容需要参看文献[1]的补充资料）：

OK，介绍完如何选取聚类中心，接下来记录一下如何做聚类的过程。前面我们提到，已经有一个关于样本点局部密度的降序排列，在确定完中心点后，对余下的点，我们依次做下面的操作：
对于未标记聚类类别的点a，找到密度比它大的点中和它距离最近的那个点c，把c点的类别标号赋给a点。根据局部密度的降序依次给所有样本点进行类别赋值，就可以一遍得到所有样本点的聚类结果。下图（来自于[3]）简单说明了这个过程，先确定了点1和2为中心，然后点3（找到点1）标记为聚类1,点4（找到点3）也标记为聚类1。可以看到，聚类这个过程是一个扫描一遍的过程，是线性的。

当然，本方法还有一点是有区别于其他聚类方法的，它还做了outlier的过滤。在前面第一幅图中，26-28三个点被认为是outlier。具体的思路是这样的，首先对每一个cluster中的点，分为cluster core点和cluster halo点，halo点是最终要被排除的。需要先确定每个cluster的边界区域，边界区域由cluster中这样性质的数据点构成——它们的距离为$d_c$的邻域内，存在着其他cluster的点。计算这样跨cluster的点对(记为i和j，i为本cluster内部的点)的平均局部密度$\bar{\rho}=\frac{1}{2}(\rho_i+\rho_j)$，并取一个cluster的所有存在的跨cluster点对的平均局部密度的最大值作为该cluster的平均局部密度上界；然后对所有样本点进行一次标记，如果一个样本点的局部密度小于它所在类别的平均局部密度上界，则被标记为halo点，否则标记为core点；halo点是被抛弃的离群点。看下面这个示意图[3]：

一开始点1,3，6,7都被标记为cluster 1，点6确定了cluster 1的边界区域，并计算出了该类的平均局部密度上界（通过6-5和6-4），由于点7的局部密度小于该平均局部密度上界，被标记为halo点，排除在cluster 1之外。

OK，到这里就讲完了聚类过程的核心过程思想了，但是整个程序的执行还有很多小的细节点，比如如何确定$d_c$？作者提出的方法是经验性的，设定为所有点的相互距离中由小到大排列占总数2%位置的距离数值。当然，这个2%也是经验性质的，事实上该数值的选取是非常影响算法最终的结果的，太大太小都不行，大概2%在作者测试的data set里面还算不错。这样通过data depended的方法来选取$d_c$相对来说比直接设定一个具体的数值来的鲁棒性更好一些。
下面来看一下实验结果（看起来肯定是很棒的啦，所有的聚类paper都号称自己方法可以做到很好的效果，不过实际上最重要的是鲁棒性以及参数调节难度低，对用户来说参数调节最好自动或者没有/很少参数，这也是为什么kmeans算法依然在实际应用中使用的最多的原因——使用方便）

下图是一个非球形类分布图，同时加入黑色噪音点后，A图为类的概率分布，B、C图为4000个和1000个样本点，E和F生成的每个点对应的$\rho$和$\delta$的函数图，可以明显看出类别中心及个数。F图为随着样本点的增加，错误指派点的比率[3]。

以下是使用该算法在其他数据集上进行聚类的效果图

下面是算法执行流程，在这里收录一下[2]，：

参考资料
[1] Alex Rodriguez and Alessandro Laio. Clustering by fast search and find of density peaks, Science 344, 1492 (2014);
[2] http://blog.csdn.net/itplus/article/details/38926837
[3] http://blog.csdn.net/lvxiong1990/article/details/40540065