关于密度聚类的一种GPU加速算法

论文来源 ：《G-DBSCAN: A GPU Accelerated Algorithm for Density-based Clustering》

这篇论文提出了一种基于GPU的密度聚类算法，简称G-DBSCAN。在摘要中，作者说G-DBSCAN的速度比使用cpu的序列化DBSCAN要快100倍。

那G-DBSCAN具体是怎么实现的呢？请看下文。

在introduction部分，作者说这个算法的实现策略主要分为两步。

第一步是构图，每一个object当成图的一个节点，当两个节点之间的相似度（例如欧几里得距离）小于等于某个阈值时，就在这两个节点之间构建一条边。

第二步是识别簇，通过传统的广度优先搜索遍历整个图。

那么具体细节是什么样呢？请继续往下看。

首先，在介绍G-DBSCAN之前，作者首先把DBSCAN的基本概念介绍了一下。

DBSCAN是密度聚类算法的一种，主要思想是通过定义邻域来对集群中的对象进行分组，每一组就是一簇。

这其中涉及到三个概念：core，border，noise。这三个概念都与邻域有关。

那什么是邻域呢？

假设给出一个阈值R，对于对象A来说，以A为中心，与A之间的距离在R以内的空间就是A的邻域。这个距离可以由不同的计算方法得出，比如欧几里得距离（欧氏距离）或者余弦相似度。

阈值T需要人为指定，是算法的一个输入参数。

如果一个对象A的邻域内包含超过M个其他对象，那么A就是一个核心对象，也就是core。被包含在A的邻域内的对象就是border，可以理解为成员对象。

M也是一个阈值，可以理解为邻域最小成员个数，M需要人为指定，同样的，它也是算法的一个输入参数。

对了，还得补充俩概念：密度直达和密度可达。

密度直达指的是，如果A是核心对象，B是A的成员对象，那么B对于A来说就是密度直达的关系，可称为A密度直达B。

密度可达指的是，如果A密度直达B，B密度直达C，那么A密度可达C。

如果一个节点和其他任何节点之间都不存在密度直达或密度可达关系，那么这个节点就被视为一个noise。

这篇论文里画了一张图，还是比较形象的，如下所示：

 和其他传统的聚类算法不一样，DBSCAN不需要预先指定簇的数量，但是需要预先指定和簇密度相关的参数，也就是上文提到的阈值R和阈值M（论文里用的是MinPts）。

那么G-DBSCAN是怎么实现的呢？

作者说了，首先呢，把所有数据用图的形式来表示，也就是G(V, E)，V表示所有对象，E表示对象之间的边。注意，只有两个对象之间的距离小于等于阈值R，它们之间才会连一条边。

考虑到G-DBSCAN的主要目的是为了实现并行计算，这对内存有一定要求，因此作者决定采用紧凑一点的邻接表来表示图。

这主要通过两个向量来实现。一个是向量Va，存储所有的节点。一个是向量Ea，存储所有的邻接表。

Va的索引表示节点，每个向量元素存储俩值：该节点的邻接节点的数目，以及在向量Ea中该节点的邻接表开始的位置。

如下图所示，对于左边小图中的节点0，我们暂且称它小玲。在O(V)中，第一列元素表示小玲对应的值，上面的2表示小玲的邻接节点的个数为2（从左边小图能看出小玲和2、3相连，2和3就是小玲的邻接节点），下面的0表示小玲的邻接表在O(E)中的起始位置是0。那么对于小玲，在O(E)中从0开始的位置，依次往后推两个位置，每个位置所对应的节点就是小玲的邻接节点，从O(E)中可以看出是2和3 ，这与左边的小图是一致的。

 构图过程如下面的伪代码所示。首先输入阈值M和R，以及数据集。对于数据集中的每一个对象，计算它与其他对象的距离（论文中用的是欧氏距离），如果小于阈值R，那么在这两个对象之间构建一条边。等到所有对象与其他对象的距离都计算完毕，边也构造完毕后，这时候就能知道一个对象有几条边连接着其他对象，这个边的数量就是这个对象的邻域节点的数量。这个时候，利用一个函数ClassifyObject()去判断每个对象是core还是noise。判断标准和前文介绍到的关于core、border、noise的概念略有差别，此处的判断标准是，如果一个对象的邻域节点数量大于阈值M，那么该对象就是core，否则它就是noise。那么有同学就会问了，只有core和noise，那border哪去了？别着急嘛，下文会讲。

 以上呢，是G-DBSCAN的第一部分。这部分的复杂度为O(V^2)，因为最坏的情况是每两个节点之间都需要计算一次距离。

接下来介绍G-DBSCAN的第二部分。如下图的伪代码所示，在构建了图，并且对节点做了分类后，对于图中的每一个core节点，通过广度优先搜索算法（BFS）去遍历所有与core节点相连的节点，直到遍历结束。所有被遍历到的节点就被分类成border节点，与出发时的core节点一起，成为一个簇。在一次BFS结束后，如果还有core节点没被遍历到，那么就再从这个从core节点出发，再进行BFS，每次BFS结束都会生成一个簇。直到数据中的所有的core节点都被遍历过。 

 以上就是G-DBSCAN的第二部分，这部分的计算复杂度主要取决于BFS算法，BFS的复杂度为O(V+E) 。

综合第一部分和第二部分，G-DBSCAN的复杂度为O(V^2)。

好了，以上呢，就是G-DBSCAN的基本原理介绍。

接下来就是重头戏了！

究竟 G-DBSCAN是怎么实现并行计算的呢？请君喝点水，继续往下看。

G-DBSCAN的并行化实现

一、首先是构图过程怎么实现并行化

如前文所述，构图的过程主要通过两个向量实现，一个是Va，一个是Ea。

这两个向量的计算过程分以下三步：

1. 计算Va的第一个值，也就是各个节点的度数（节点与其他节点连接的总边数）。

     通过GPU给每个节点分配一个线程，该线程会计算它所负责的节点的所有邻接节点数量，这个过程是可以并行实现的（互相之间并不存在依赖关系），即所有线程同时开始计算邻接节点数目，并将计算结果填入Va向量中第一个值的对应位置（Va向量中每个位置对应一个节点，每个位置有两个值）。

2. 计算Va的第二个值，也就是各个节点对应的邻接列表索引。

        Va第二个值对应的是节点的邻接列表在Ea中的起始位置，主要取决于前一个节点的起始位置和度数（邻接节点个数）。比如说节点0是第一个节点，其实位置是0，度数是2，那么对于下一个节点，也就是节点1来说，其起始位置是上一个节点的起始位置加度数，也就是0+2，即节点1的起始位置为2。依此类推，如果节点1的度数是1，那么对于节点2来说，其起始位置为2+1（2是节点1 的起始位置，1是节点1的度数）=3。

        由此可见，这Va第二个值的计算，互相之间是有依赖关系的。针对这个问题，论文作者说这个问题能通过使用独占扫描操作（exclusive scan）高效完成。具体点说，就是他们使用推力库（thrust library）实现了并行计算。然而，文中只说推力库是CUDA SDK分发的一部分，具体推力库是怎么回事没细说。笔者查了查推力库的概念，也没查到有效资料，暂且就留一个坑在这里，日后再补吧。

3. 构造邻接表的集合。

        经过步骤1和2，Va就计算完成啦。接下来Ea的计算就比较简单了，和之前一样，给每个节点分配一个GPU线程，用于构造该节点的邻接表，表里填充的是所有与该节点邻接的节点。相当于把Ea分布到每一个线程中，一个线程负责维护一个节点的邻接表。

二、其次是广度优先搜索（BFS）怎么实现并行化

        对于BFS的并行化，作者采用层级同步（level synchronization）去实现。也就是说，以层级的形式去遍历图，被遍历过的层级将不再被遍历。关于图的层级同步，这个博主写了博客，地址是CUDA编程——图的表示及BFS算法_NKU_Yang的博客-CSDN博客_bfs cuda。

        在具体实现时，执行以下过程：

        首先给出两个布尔向量，分别是Fa和Xa，向量大小和当前层级的节点个数一致。Fa和Xa的所有初始值都为False。Fa的每个值对应于该节点是否开始执行基于GPU-BFS，若为True，则开始执行，若为False，则不执行。Xa中的每个值则对应于该节点是否被访问过，若被访问过，则置为True，否则置为False。

        接下来给每个节点都分配一个线程。

        如下图中的伪代码所示，指定一个core节点作为起始节点，该节点对应的Fa值设为True，从该节点开始，按照层级同步的方式，同一层级的节点并行执行GPU-BFS。

        对于起始节点，搜索开始后，该位置对应的Xa值被置为True，表示该节点已被访问，其邻居节点对应的Fa值则置为True。搜索结束后，起始节点对应的Fa值从Fa中移除。而邻居节点则作为起始节点继续迭代上述过程。直到Fa中的值全部移除，即Fa为空。

        下图中的伪代码与论文中描述的过程似乎有些不符，论文中说起始在搜索结束后，会从Fa中移除，但是伪代码中只是将对应的Fa值置为False。

         经过以上过程，起始节点已经通过BFS搜索找到了其密度可达的所有对象，这些对象如果有的是core，那么更改为border，和起始节点构成一个簇。

        如果一次CPU-BFS结束后，还有core节点没被访问，那么从这个没被访问的core节点开始继续进行GPU-BFS，从而得到新的簇。

        该过程循环执行，直到所有core节点均被访问过。

        所有core节点按照上面这种方式即可找到自己密度可达的所有对象，从而聚成一个个簇。

        好了，以上就是BFS的并行化实现（实际上还是有一些迷惑的点没搞明白，对于层级的概念不是很了解，所以这部分写得连自己也觉得模棱两可，就先这样吧，日后想明白了再来补坑。）

三、实验部分

        论文里说实验部分是通过基于C语言的CUDA编程实现的。

        数据规模从5千一直增加到了70万，数据维度是2维。从规模和数据维度上看，相较于实际业务中动辄上千万的数据量，论文的数据集规模着实是有点小，而且数据达到60万以后，时间消耗的增速就不再提升了。

        从实验结果看，比起CPU实现的密度聚类，G-DBSCAN降低了时间消耗。

        在构图这一步，时间消耗降低了82倍，在BFS搜索中，时间消耗降低了21倍。综合全过程，时间消耗降低了112倍。

四、不足之处

        这篇论文看完，感觉看了个寂寞。

        从算法角度看，复杂度并没有得到降低，仍为O(V^2)，只是通过CUDA编程的方式，使得构图和广度优先搜索的一些中间过程能并行执行，虽然降低了时间开销，但是计算资源上的开销也相应增大了，而且时间上的增速存在瓶颈，到达某个值后就不再增加了。

        不过比起sklearn里的包来讲，这篇论文一定程度上提供了将密度聚类部署到GPU的思路，毕竟sklearn里的包只能用于CPU，且处理不了大规模数据，一旦维度大于10的数据量增加到10万以上，CPU就会负载过高，然后进程宕机。。。

        算法道阻且长，我辈仍需努力啊。。。