基于滑动门中心点计算的K均值聚类并行算法研究

龚运鸿，周新志，雷印杰

(四川大学 电子信息学院，成都 610065 )

0 引言

随着程序需要处理的数据量越来越庞大，现如今以GB或TB为单位的数据集已经十分普遍了，数据挖掘中必须重视的一个问题就是如何高效得处理如此庞大的数据。即使算法的复杂程度是线性增长的，时间和空间的消耗也不容忽视。K均值聚类算法由Stuart Lloyd等人在1957年第一次提出[1]，K均值聚类算法源于信号处理的一种矢量量化方法，由于其概念简单、收敛速度快、易于实现等特点，现如今在数据挖掘领域聚类分析中十分流行。然而K均值聚类算法的复杂度比较高，如何高效进行算法计算是一个急切的研究方向。目前，陶冶[2]等人证明并实现了并行K均值聚类算法。喻金平[3]等人提出一种改进的人工蜂群算法，解决了K均值算法的搜索能力差的问题。霍莹秋[4]等人提出分块、并行的K均值聚算法，采用“合并访问”、“多级规约求和”和“负载均衡”等优化策略优化并行算法，提高了算法的运行速度。对于K均值聚类并行计算的研究还存在许多缺陷，比如没有针对CUDA并行计算平台进行优化，也没有针对中心点更新效率问题提出解释等。根据以上研究的缺陷，本文利用NVIDIA的CUDA并行平台，在传统K均值算法基础上，采用了一种滑动门并行计算中心点算法优化K均值算法更新中心点的耗时问题。通过与规约法计算中心点算法相比，获得了很好的加速比。

1 CUDA并行计算平台

随着社会的发展，CPU逐渐达到了速度极限并且购置成本也在急速上升。通过对显示图像的优化，GPU技术却逐渐完善了起来，在通用计算领域，GPU的能力逐渐超过了传统CPU。通过近几年的发展，计算技术正在从单一CPU串行计算方式向GPGPU并行协同计算发展。为了让开发者无需学习复杂的着色语言和图像处理原语，方便更多的开发者能够进行GPU编程，显卡厂家NVIDIA在2007年提供了一个方便开发者使用的接口——Compute Unified Device Architecture,即CUDA。CUDA可以让开发者直接访问GPU的虚拟指令设定和并行计算元素。与在GPU上使用图形API进行计算的传统方式相比，CUDA具有十分巨大的优势：1)程序可以在内存的任意位置进行读取；2)在CUDA4.0以上的版本拥有统一虚拟内存；3)在CUDA6.0以上的版本拥有统一内存；4)CUDA为开发者开辟一个快速共享内存区域，其数据可以在线程中共享。这可以帮助开发者管理缓存、建立更高的带宽；5)可以更快地从GPU上下载和回写数据。在科技研究方面，CUDA技术可谓炙手可热，越来越多的研发人员投入到了CUDA并行技术研究当中，所以在科研领域，CUDA也拥有了许多研究成果：在医疗领域，可以将CUDA加速技术运用在CT或者MRI扫描图像的VR技术上；在物理建模上，特别是流体动力学领域有很好的成果；在神经网络训练领域中，CUDA技术可以很好地解决机器学习遇到的瓶颈。

CUDA将C语言进行扩展，这样可以使开发者使用C语言等高级语言在GPU设备上进行并行程序编程，方便开发者使用。使用CUDA进行编写的代码，既可以在CPU上使用，也可以在GPU上使用。使用CUDA并行计算平台时，主机端为CPU，设备端为GPU。在GPU运行时，基于NVIDIA自身底层硬件特点，采用了单程序多数据(SPMD)的并行计算方式来处理数据，属于单指令多数据(SIMD)的一种变体。并行编程的核心是线程的概念，运行的程序称为内核函数(kernel)，并行计算时，每一个线程都会同时处理一个kernel，CUDA中数据通过线程-块-网格计算方式进行分配。每个线程块都有自己唯一的识别id，一个或者多个线程块会由流多处理器SM来并行处理，每一个SM由很多个32位寄存器组成。每个SM中有8个流处理器(SP)构成。每一个流处理器都有一个共享存储器，所有SP都可以访问共享存储器中的内容。如果线程块是一维的，那么blockId.x就可以确定线程块的id了。如果线程块是二维的，那么需要blockId.x和blockId.y才能确定线程块的id。每一个线程也有自己唯一识别id，通过这个id查找需要处理数据的位置。如果线程是一维的，那么threadIdx.x就可以确定线程的id了。如果线程是二维的，那么需要threadIdx.x和threadIdx.y才能确定线程的id。如果线程是三维的，那么需要threadIdx.x，threadIdx.y，threadIdx.z这3个参数才能确定线程的id。同一线程块中的线程如果需要进行数据传输，一般是在共享存储器中进行的。

2 K均值聚类算法

2.1 K均值聚类算法介绍

聚类是一种无监督的学习，它将相似的对象归到同一簇中，簇内的对象越相似，聚类的效果越好。K均值聚类算法可以发现k个不同的簇，且每个簇的中心采用簇中所含值的均值计算而成。

K均值聚类算法是将含有n个数据点划分为k个簇Cj(j=1,2,...,k;k≤n)。首先在n个数据点中随机选取k个数据点代表k个簇的质心，再根据剩余数据点与k个质心的距离，将剩余数据点划分到与其距离最近的簇Cj中。然后重新计算每个簇中所有值的均值，并将这个均值作为新的质心。该过程不断重复，直到误差平方和函数SSE收敛，其定义如公式(1)所示：

(1)

(2)

式中,k代表选择的簇的数量；Cj代表第j(j=1,2,...,k;k≤n)个簇；x代表簇Cj中的任意一个数据点；cj是簇Cj的均值；mj代表簇Cj中数据点的个数。

2.2 K均值聚类流程

1)初始化聚类中心。在n个数据点钟随机选取k个数据点作为初始聚类中心C1,C2,...,Ck。

2)计算剩余数据点到每个初始聚类中心的距离，将剩余数据点几个划分到与其最近的簇中心代表的簇中，根据公式(1)计算SSE的值。

3)分别算出各个簇中所有数据点的均值，用这些均值替换初始聚类中心，用新的聚类中心重复步骤2)。

4)将上一次SSE的值和本次SSE的值进行比较，差值绝对值大于阀值，代表SSE收敛，则进行步骤5)，否则进行步骤5)。

5)算法结束。

排除其他因素只考虑算法核心部分，把每一次的比较、乘法、加法操作都当做一次浮点运算，用Tflop来代表一次浮点运算花费的时间，则算法核心部分每次迭代所用的时间见表1。

表1 K均值聚类算法核心部分所用运算时间

3 数据准备阶段并行设计

3.1 初始化数据的并行设计

选取初始化聚类中心的过程，不是完全随机产生的，需要先确定数据点在所有维度的最大最小值，如图1所示。

图1 二维数据点分布图

在选取初值时，求最大最小值是典型的规约运算，而规约运算是可以并行化的。对于n个输入数据和操作⊕，规约可表示为：

⊕a1⊕a2⊕...⊕an

(3)

图2展示了处理N个元素规约操作的实现。

图2 3种不同规约操作的实现

由图2可以发现，不同的实现方式，其时间复杂度也是不同的。其中，串行实现完成规约操作需要n-1步，两种对数步长的规约操作只需要lgn步，大大减少了时间复杂度。但是，由于成对方式不能合并内存事务，成对方式在CUDA实现中性能比较差。在CUDA中，交替方式在全局内存和共享内存都有很好的效果。在全局内存中，将blockDim.x*gridDim.x的倍数作为作为交替因子，这样可以获得比较好的性能。在共享内存中，需要保持线程块相邻的线程保持活跃状态，同时，为了避免内存的冲突，需按确定的交替因子来累计结果，这样可以获得良好的效果。

3.2 分配数据的并行设计

基于GPU的K均值聚类算法在进行数据对象分配这一步骤的时候有两种策略。第一种策略是面向每个簇的中心，通过计算出该簇的中心与每个数据对象的距离，然后将每个数据对象归并到距离簇中心最近的那个簇中。这种策略适应于GPU的处理核心数量比较少的情况，此时，GPU中的每个处理核心可以对应一个簇的中心，并且能够连续的处理所有的数据对象。另一种策略是面向每个数据对象，每个数据对象计算与所有簇中心的距离，然后数据对象将会被分配到距离簇中心最近的那个簇当中。该策略适应于GPU的处理核心比较大的情况，目前主流的GPU一般都有超过100个处理核心，因此第二种处理策略比较合适。

图3 数据分配过程不同策略对比图

4 滑动门中心点并行算法设计

在目前的CUDA并行计算算法优化中,一般使用带状划分的并行划分方法[5-8]，如图4所示，K均值聚类算法的中心点存储在共享内存中，而样本保存在全局存储器中，每一行数据由一个线程进行处理。带状划分方法能最大化地利用GPU的计算内核。

采用带状划方法的方法，线程1～m将在同一时间去读取共享储存器索引位置0的数据。在CUDA中，寄存器的访问速度最快。所以为了减少重复访问的时间，可以把共享存储器中的数据转化到寄存器当中进行存储。

图4 带状划分的并行聚类算法

采用滑动门并行算法时，也是采用带划分的划分方式，在每个块上开启m个线程，相同簇内的样本都会在这些线程中进行计算。每个块中的样本点值最后合并至临时中心点存储区s_cust中。m的最佳值由式(4)计算。

(4)

式中，L代表样本向量的维度；wsize代表wrap的大小；Ci代表簇i中的样本总数；tNum代表每个block中可以开启的最大线程数量；SMsize代表共享内存的大小。

在长度为L的s_cust上分配大小为m的滑动门，block内的线程按照线程号分配其计算的数据。在同一时刻，block内所有线程的存储数据均或落在滑动门范围内，且每个线程会计算独立维度的数据，这样就避免了存储数据时的线程同步。并行计算完m个数据之后，滑动门会向前移动，同时线程也会向前启动，计算下一批m个数据，直到滑动门回到初始位置为止。滑动门并行计算中心点的算法过程伪代码如下：

图4 滑动门并行计算中心点

输入：所有样本的聚类结果；总数为n的样本集D，其中样本为L维向量；簇数目k。

输出：所有簇的中心点集合U

for每个簇do

1)根据(4)选取合适的m值。

2)for j=ni,do:

(1) j = j/m

(2) for i=0 to L

for 每个线程 paraplle do：

s=blockIdx.x × block_Size + (theadIdx.x +i)%L;

s_cust[i] += data[s];

end loop;

(3) s=blockIdx.x × block_Size + i;

data[blockIdx.x × L + i]=s_cust[i];

end loop;

3)or 每个线程 paralle do:

U[threadIdx.x] = data[0][threadIdx.x] / ;

end loop

根据实验环境的GPU配置，由(1)可以计算出最合适的m值。通过输入的总数为n的样本集，可以计算出程序所需的迭代次数为L×logmn。其中，步骤(2)的迭代次数是L，m个维的值将会通过滑动门的方式被并行计算。同一个块中所有

的线程都会对存储器s_cust中的值进行访问，为了降低访问数据时的延迟，可以采用共享存储器来存储s_cust中的所有的值。通过步骤3)，共享存储器中的所有的数据都会按照块的编号转移到全局存储器当中，这样就不用重新分配其他的存储空间，可以降低运算的速度。当步骤2)运行结束之后，该簇内的样本的和将会存储在全局存储器的的第一行中，中心点存储区将会保存根据步骤3)得到的均值。通过分析，滑动门中心点并行算法的时间复杂度为logmn，相比传统规约法的计算中心点的时间复杂度相比，得到了明显的提高。

4 实验结果

为了验证滑动门并行算法的有效性，本实验对一个含有25789个数据点的样本集进行K均值聚类算法并行计算，采用单个簇一次迭代中分别使用滑动门法与规约法在不同m值下并行计算中心点的运行时间比。

表2 实验环境

表3 滑动门法与规约法计算中心点时间加速比

通过表(3)分析，可以看出当m值分配较小值时，传统规约并行算法可以得到比较好的加速比，但是与滑动门法相比，优势并不十分明显，只有0.05的优势。但是，当分配的m达到相当规模之后，滑动门中心点并行算法的优势逐渐显示了出来，并且优势逐渐扩大，从m=64的0.09到m=512时的0.47。

5 结语

本文基于CUDA平台，对K均值聚类算法的并行加速计算研究，讨论了CUDA在实现并行计算的关键技术。本文介绍了CUDA平台的基本内容和K均值聚类算法的基本原理。在并行步骤中，基于以往研究的不足之处，提出了对初始值规约计算的并行方法进行的讨论和选择、数据分配时的并行策略选择和计算中心点的并行滑动门方式，这些方法提高了GPU的利用率，降低了数据获取延迟，充分发挥了CUDA架构的并行计算优势。通过与传统规约并行方法相比，滑动门并行计算算法拥有更好的加速比。通过实验证明了CUDA在并行计算方面的有效性，结果表明了在GeForce730上实现了算法2.8倍的加速比，有效地提高了K均值算法在计算机上的运行效率。

[1]LloydS.LeastsquaresquantizationinPCM[J].IEEETransactionsonInformationTheory,1982,28 (2): 129-137.

[2] 陶 冶,曾志勇,余建坤,等.并行k均值聚类算法的完备性证明与实现[J].计算机工程, 2010, 36 (22) :72-74.

[3] 喻金平,郑 杰,梅宏标,等. 基于改进人工蜂群算法的K均值聚类算法[J]. 计算机应用,2014,34(4):1065-1069,1088.

[4] 霍迎秋，秦仁波，刑彩燕，等.基于CUDA的并行K-means聚类图像分割算法优化[J].农业机械学报,2014,45(11):72-74.

[5]NguyenH.GPUGems[M].Boston:Addison-WesleyProfession,2008.L

[6]RezaR,DanielR,EllickC,etal.Aparallelimplementationofk-meansclusteringonGPUs[A].ProcofInternationalConferenceonParallelandDistributedProcessingTechniquesandApplications[C].Springer-Verlag, 2008:340-345.

[7]MarioZ,MichaelG.AcceleratingK-meansonthegraphicsprocessorviaCUDA[A].Procofthe1stInternationalConferenceonIntensiveApplicationsandServices[C]:IEEE,2009:7-15.

[8]BaiHT,HeLL,OuyangDT,etal.K-meansoncommodityGPUswithCUDA[A].ProcofWRIWorldCongressonComputerScienceandInformationEngineering[C].ACMPress, 2009:651-655.