龔運(yùn)鴻，周新志，雷印杰

(四川大學(xué) 電子信息學(xué)院，成都 610065 )

0 引言

隨著程序需要處理的數(shù)據(jù)量越來越龐大，現(xiàn)如今以GB或TB為單位的數(shù)據(jù)集已經(jīng)十分普遍了，數(shù)據(jù)挖掘中必須重視的一個(gè)問題就是如何高效得處理如此龐大的數(shù)據(jù)。即使算法的復(fù)雜程度是線性增長的，時(shí)間和空間的消耗也不容忽視。K均值聚類算法由Stuart Lloyd等人在1957年第一次提出[1]，K均值聚類算法源于信號(hào)處理的一種矢量量化方法，由于其概念簡(jiǎn)單、收斂速度快、易于實(shí)現(xiàn)等特點(diǎn)，現(xiàn)如今在數(shù)據(jù)挖掘領(lǐng)域聚類分析中十分流行。然而K均值聚類算法的復(fù)雜度比較高，如何高效進(jìn)行算法計(jì)算是一個(gè)急切的研究方向。目前，陶冶[2]等人證明并實(shí)現(xiàn)了并行K均值聚類算法。喻金平[3]等人提出一種改進(jìn)的人工蜂群算法，解決了K均值算法的搜索能力差的問題。霍瑩秋[4]等人提出分塊、并行的K均值聚算法，采用“合并訪問”、“多級(jí)規(guī)約求和”和“負(fù)載均衡”等優(yōu)化策略優(yōu)化并行算法，提高了算法的運(yùn)行速度。對(duì)于K均值聚類并行計(jì)算的研究還存在許多缺陷，比如沒有針對(duì)CUDA并行計(jì)算平臺(tái)進(jìn)行優(yōu)化，也沒有針對(duì)中心點(diǎn)更新效率問題提出解釋等。根據(jù)以上研究的缺陷，本文利用NVIDIA的CUDA并行平臺(tái)，在傳統(tǒng)K均值算法基礎(chǔ)上，采用了一種滑動(dòng)門并行計(jì)算中心點(diǎn)算法優(yōu)化K均值算法更新中心點(diǎn)的耗時(shí)問題。通過與規(guī)約法計(jì)算中心點(diǎn)算法相比，獲得了很好的加速比。

1 CUDA并行計(jì)算平臺(tái)

隨著社會(huì)的發(fā)展，CPU逐漸達(dá)到了速度極限并且購置成本也在急速上升。通過對(duì)顯示圖像的優(yōu)化，GPU技術(shù)卻逐漸完善了起來，在通用計(jì)算領(lǐng)域，GPU的能力逐漸超過了傳統(tǒng)CPU。通過近幾年的發(fā)展，計(jì)算技術(shù)正在從單一CPU串行計(jì)算方式向GPGPU并行協(xié)同計(jì)算發(fā)展。為了讓開發(fā)者無需學(xué)習(xí)復(fù)雜的著色語言和圖像處理原語，方便更多的開發(fā)者能夠進(jìn)行GPU編程，顯卡廠家NVIDIA在2007年提供了一個(gè)方便開發(fā)者使用的接口——Compute Unified Device Architecture,即CUDA。CUDA可以讓開發(fā)者直接訪問GPU的虛擬指令設(shè)定和并行計(jì)算元素。與在GPU上使用圖形API進(jìn)行計(jì)算的傳統(tǒng)方式相比，CUDA具有十分巨大的優(yōu)勢(shì)：1)程序可以在內(nèi)存的任意位置進(jìn)行讀??；2)在CUDA4.0以上的版本擁有統(tǒng)一虛擬內(nèi)存；3)在CUDA6.0以上的版本擁有統(tǒng)一內(nèi)存；4)CUDA為開發(fā)者開辟一個(gè)快速共享內(nèi)存區(qū)域，其數(shù)據(jù)可以在線程中共享。這可以幫助開發(fā)者管理緩存、建立更高的帶寬；5)可以更快地從GPU上下載和回寫數(shù)據(jù)。在科技研究方面，CUDA技術(shù)可謂炙手可熱，越來越多的研發(fā)人員投入到了CUDA并行技術(shù)研究當(dāng)中，所以在科研領(lǐng)域，CUDA也擁有了許多研究成果：在醫(yī)療領(lǐng)域，可以將CUDA加速技術(shù)運(yùn)用在CT或者M(jìn)RI掃描圖像的VR技術(shù)上；在物理建模上，特別是流體動(dòng)力學(xué)領(lǐng)域有很好的成果；在神經(jīng)網(wǎng)絡(luò)訓(xùn)練領(lǐng)域中，CUDA技術(shù)可以很好地解決機(jī)器學(xué)習(xí)遇到的瓶頸。

CUDA將C語言進(jìn)行擴(kuò)展，這樣可以使開發(fā)者使用C語言等高級(jí)語言在GPU設(shè)備上進(jìn)行并行程序編程，方便開發(fā)者使用。使用CUDA進(jìn)行編寫的代碼，既可以在CPU上使用，也可以在GPU上使用。使用CUDA并行計(jì)算平臺(tái)時(shí)，主機(jī)端為CPU，設(shè)備端為GPU。在GPU運(yùn)行時(shí)，基于NVIDIA自身底層硬件特點(diǎn)，采用了單程序多數(shù)據(jù)(SPMD)的并行計(jì)算方式來處理數(shù)據(jù)，屬于單指令多數(shù)據(jù)(SIMD)的一種變體。并行編程的核心是線程的概念，運(yùn)行的程序稱為內(nèi)核函數(shù)(kernel)，并行計(jì)算時(shí)，每一個(gè)線程都會(huì)同時(shí)處理一個(gè)kernel，CUDA中數(shù)據(jù)通過線程-塊-網(wǎng)格計(jì)算方式進(jìn)行分配。每個(gè)線程塊都有自己唯一的識(shí)別id，一個(gè)或者多個(gè)線程塊會(huì)由流多處理器SM來并行處理，每一個(gè)SM由很多個(gè)32位寄存器組成。每個(gè)SM中有8個(gè)流處理器(SP)構(gòu)成。每一個(gè)流處理器都有一個(gè)共享存儲(chǔ)器，所有SP都可以訪問共享存儲(chǔ)器中的內(nèi)容。如果線程塊是一維的，那么blockId.x就可以確定線程塊的id了。如果線程塊是二維的，那么需要blockId.x和blockId.y才能確定線程塊的id。每一個(gè)線程也有自己唯一識(shí)別id，通過這個(gè)id查找需要處理數(shù)據(jù)的位置。如果線程是一維的，那么threadIdx.x就可以確定線程的id了。如果線程是二維的，那么需要threadIdx.x和threadIdx.y才能確定線程的id。如果線程是三維的，那么需要threadIdx.x，threadIdx.y，threadIdx.z這3個(gè)參數(shù)才能確定線程的id。同一線程塊中的線程如果需要進(jìn)行數(shù)據(jù)傳輸，一般是在共享存儲(chǔ)器中進(jìn)行的。

2 K均值聚類算法

2.1 K均值聚類算法介紹

聚類是一種無監(jiān)督的學(xué)習(xí)，它將相似的對(duì)象歸到同一簇中，簇內(nèi)的對(duì)象越相似，聚類的效果越好。K均值聚類算法可以發(fā)現(xiàn)k個(gè)不同的簇，且每個(gè)簇的中心采用簇中所含值的均值計(jì)算而成。

K均值聚類算法是將含有n個(gè)數(shù)據(jù)點(diǎn)劃分為k個(gè)簇Cj(j=1,2,...,k;k≤n)。首先在n個(gè)數(shù)據(jù)點(diǎn)中隨機(jī)選取k個(gè)數(shù)據(jù)點(diǎn)代表k個(gè)簇的質(zhì)心，再根據(jù)剩余數(shù)據(jù)點(diǎn)與k個(gè)質(zhì)心的距離，將剩余數(shù)據(jù)點(diǎn)劃分到與其距離最近的簇Cj中。然后重新計(jì)算每個(gè)簇中所有值的均值，并將這個(gè)均值作為新的質(zhì)心。該過程不斷重復(fù)，直到誤差平方和函數(shù)SSE收斂，其定義如公式(1)所示：

(1)

(2)

式中,k代表選擇的簇的數(shù)量；Cj代表第j(j=1,2,...,k;k≤n)個(gè)簇；x代表簇Cj中的任意一個(gè)數(shù)據(jù)點(diǎn)；cj是簇Cj的均值；mj代表簇Cj中數(shù)據(jù)點(diǎn)的個(gè)數(shù)。

2.2 K均值聚類流程

1)初始化聚類中心。在n個(gè)數(shù)據(jù)點(diǎn)鐘隨機(jī)選取k個(gè)數(shù)據(jù)點(diǎn)作為初始聚類中心C1,C2,...,Ck。

2)計(jì)算剩余數(shù)據(jù)點(diǎn)到每個(gè)初始聚類中心的距離，將剩余數(shù)據(jù)點(diǎn)幾個(gè)劃分到與其最近的簇中心代表的簇中，根據(jù)公式(1)計(jì)算SSE的值。

3)分別算出各個(gè)簇中所有數(shù)據(jù)點(diǎn)的均值，用這些均值替換初始聚類中心，用新的聚類中心重復(fù)步驟2)。

4)將上一次SSE的值和本次SSE的值進(jìn)行比較，差值絕對(duì)值大于閥值，代表SSE收斂，則進(jìn)行步驟5)，否則進(jìn)行步驟5)。

5)算法結(jié)束。

排除其他因素只考慮算法核心部分，把每一次的比較、乘法、加法操作都當(dāng)做一次浮點(diǎn)運(yùn)算，用Tflop來代表一次浮點(diǎn)運(yùn)算花費(fèi)的時(shí)間，則算法核心部分每次迭代所用的時(shí)間見表1。

表1 K均值聚類算法核心部分所用運(yùn)算時(shí)間

3 數(shù)據(jù)準(zhǔn)備階段并行設(shè)計(jì)

3.1 初始化數(shù)據(jù)的并行設(shè)計(jì)

選取初始化聚類中心的過程，不是完全隨機(jī)產(chǎn)生的，需要先確定數(shù)據(jù)點(diǎn)在所有維度的最大最小值，如圖1所示。

圖1 二維數(shù)據(jù)點(diǎn)分布圖

在選取初值時(shí)，求最大最小值是典型的規(guī)約運(yùn)算，而規(guī)約運(yùn)算是可以并行化的。對(duì)于n個(gè)輸入數(shù)據(jù)和操作⊕，規(guī)約可表示為：

⊕a1⊕a2⊕...⊕an

(3)

圖2展示了處理N個(gè)元素規(guī)約操作的實(shí)現(xiàn)。

圖2 3種不同規(guī)約操作的實(shí)現(xiàn)

由圖2可以發(fā)現(xiàn)，不同的實(shí)現(xiàn)方式，其時(shí)間復(fù)雜度也是不同的。其中，串行實(shí)現(xiàn)完成規(guī)約操作需要n-1步，兩種對(duì)數(shù)步長的規(guī)約操作只需要lgn步，大大減少了時(shí)間復(fù)雜度。但是，由于成對(duì)方式不能合并內(nèi)存事務(wù)，成對(duì)方式在CUDA實(shí)現(xiàn)中性能比較差。在CUDA中，交替方式在全局內(nèi)存和共享內(nèi)存都有很好的效果。在全局內(nèi)存中，將blockDim.x*gridDim.x的倍數(shù)作為作為交替因子，這樣可以獲得比較好的性能。在共享內(nèi)存中，需要保持線程塊相鄰的線程保持活躍狀態(tài)，同時(shí)，為了避免內(nèi)存的沖突，需按確定的交替因子來累計(jì)結(jié)果，這樣可以獲得良好的效果。

3.2 分配數(shù)據(jù)的并行設(shè)計(jì)

基于GPU的K均值聚類算法在進(jìn)行數(shù)據(jù)對(duì)象分配這一步驟的時(shí)候有兩種策略。第一種策略是面向每個(gè)簇的中心，通過計(jì)算出該簇的中心與每個(gè)數(shù)據(jù)對(duì)象的距離，然后將每個(gè)數(shù)據(jù)對(duì)象歸并到距離簇中心最近的那個(gè)簇中。這種策略適應(yīng)于GPU的處理核心數(shù)量比較少的情況，此時(shí)，GPU中的每個(gè)處理核心可以對(duì)應(yīng)一個(gè)簇的中心，并且能夠連續(xù)的處理所有的數(shù)據(jù)對(duì)象。另一種策略是面向每個(gè)數(shù)據(jù)對(duì)象，每個(gè)數(shù)據(jù)對(duì)象計(jì)算與所有簇中心的距離，然后數(shù)據(jù)對(duì)象將會(huì)被分配到距離簇中心最近的那個(gè)簇當(dāng)中。該策略適應(yīng)于GPU的處理核心比較大的情況，目前主流的GPU一般都有超過100個(gè)處理核心，因此第二種處理策略比較合適。

圖3 數(shù)據(jù)分配過程不同策略對(duì)比圖

4 滑動(dòng)門中心點(diǎn)并行算法設(shè)計(jì)

在目前的CUDA并行計(jì)算算法優(yōu)化中,一般使用帶狀劃分的并行劃分方法[5-8]，如圖4所示，K均值聚類算法的中心點(diǎn)存儲(chǔ)在共享內(nèi)存中，而樣本保存在全局存儲(chǔ)器中，每一行數(shù)據(jù)由一個(gè)線程進(jìn)行處理。帶狀劃分方法能最大化地利用GPU的計(jì)算內(nèi)核。

采用帶狀劃方法的方法，線程1～m將在同一時(shí)間去讀取共享儲(chǔ)存器索引位置0的數(shù)據(jù)。在CUDA中，寄存器的訪問速度最快。所以為了減少重復(fù)訪問的時(shí)間，可以把共享存儲(chǔ)器中的數(shù)據(jù)轉(zhuǎn)化到寄存器當(dāng)中進(jìn)行存儲(chǔ)。

圖4 帶狀劃分的并行聚類算法

采用滑動(dòng)門并行算法時(shí)，也是采用帶劃分的劃分方式，在每個(gè)塊上開啟m個(gè)線程，相同簇內(nèi)的樣本都會(huì)在這些線程中進(jìn)行計(jì)算。每個(gè)塊中的樣本點(diǎn)值最后合并至臨時(shí)中心點(diǎn)存儲(chǔ)區(qū)s_cust中。m的最佳值由式(4)計(jì)算。

(4)

式中，L代表樣本向量的維度；wsize代表wrap的大小；Ci代表簇i中的樣本總數(shù)；tNum代表每個(gè)block中可以開啟的最大線程數(shù)量；SMsize代表共享內(nèi)存的大小。

在長度為L的s_cust上分配大小為m的滑動(dòng)門，block內(nèi)的線程按照線程號(hào)分配其計(jì)算的數(shù)據(jù)。在同一時(shí)刻，block內(nèi)所有線程的存儲(chǔ)數(shù)據(jù)均或落在滑動(dòng)門范圍內(nèi)，且每個(gè)線程會(huì)計(jì)算獨(dú)立維度的數(shù)據(jù)，這樣就避免了存儲(chǔ)數(shù)據(jù)時(shí)的線程同步。并行計(jì)算完m個(gè)數(shù)據(jù)之后，滑動(dòng)門會(huì)向前移動(dòng)，同時(shí)線程也會(huì)向前啟動(dòng)，計(jì)算下一批m個(gè)數(shù)據(jù)，直到滑動(dòng)門回到初始位置為止?；瑒?dòng)門并行計(jì)算中心點(diǎn)的算法過程偽代碼如下：

圖4 滑動(dòng)門并行計(jì)算中心點(diǎn)

輸入：所有樣本的聚類結(jié)果；總數(shù)為n的樣本集D，其中樣本為L維向量；簇?cái)?shù)目k。

輸出：所有簇的中心點(diǎn)集合U

for每個(gè)簇do

1)根據(jù)(4)選取合適的m值。

2)for j=ni,do:

(1) j = j/m

(2) for i=0 to L

for 每個(gè)線程 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 每個(gè)線程 paralle do:

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

end loop

根據(jù)實(shí)驗(yàn)環(huán)境的GPU配置，由(1)可以計(jì)算出最合適的m值。通過輸入的總數(shù)為n的樣本集，可以計(jì)算出程序所需的迭代次數(shù)為L×logmn。其中，步驟(2)的迭代次數(shù)是L，m個(gè)維的值將會(huì)通過滑動(dòng)門的方式被并行計(jì)算。同一個(gè)塊中所有

的線程都會(huì)對(duì)存儲(chǔ)器s_cust中的值進(jìn)行訪問，為了降低訪問數(shù)據(jù)時(shí)的延遲，可以采用共享存儲(chǔ)器來存儲(chǔ)s_cust中的所有的值。通過步驟3)，共享存儲(chǔ)器中的所有的數(shù)據(jù)都會(huì)按照塊的編號(hào)轉(zhuǎn)移到全局存儲(chǔ)器當(dāng)中，這樣就不用重新分配其他的存儲(chǔ)空間，可以降低運(yùn)算的速度。當(dāng)步驟2)運(yùn)行結(jié)束之后，該簇內(nèi)的樣本的和將會(huì)存儲(chǔ)在全局存儲(chǔ)器的的第一行中，中心點(diǎn)存儲(chǔ)區(qū)將會(huì)保存根據(jù)步驟3)得到的均值。通過分析，滑動(dòng)門中心點(diǎn)并行算法的時(shí)間復(fù)雜度為logmn，相比傳統(tǒng)規(guī)約法的計(jì)算中心點(diǎn)的時(shí)間復(fù)雜度相比，得到了明顯的提高。

4 實(shí)驗(yàn)結(jié)果

為了驗(yàn)證滑動(dòng)門并行算法的有效性，本實(shí)驗(yàn)對(duì)一個(gè)含有25789個(gè)數(shù)據(jù)點(diǎn)的樣本集進(jìn)行K均值聚類算法并行計(jì)算，采用單個(gè)簇一次迭代中分別使用滑動(dòng)門法與規(guī)約法在不同m值下并行計(jì)算中心點(diǎn)的運(yùn)行時(shí)間比。

表2 實(shí)驗(yàn)環(huán)境

表3 滑動(dòng)門法與規(guī)約法計(jì)算中心點(diǎn)時(shí)間加速比

通過表(3)分析，可以看出當(dāng)m值分配較小值時(shí)，傳統(tǒng)規(guī)約并行算法可以得到比較好的加速比，但是與滑動(dòng)門法相比，優(yōu)勢(shì)并不十分明顯，只有0.05的優(yōu)勢(shì)。但是，當(dāng)分配的m達(dá)到相當(dāng)規(guī)模之后，滑動(dòng)門中心點(diǎn)并行算法的優(yōu)勢(shì)逐漸顯示了出來，并且優(yōu)勢(shì)逐漸擴(kuò)大，從m=64的0.09到m=512時(shí)的0.47。

5 結(jié)語

本文基于CUDA平臺(tái)，對(duì)K均值聚類算法的并行加速計(jì)算研究，討論了CUDA在實(shí)現(xiàn)并行計(jì)算的關(guān)鍵技術(shù)。本文介紹了CUDA平臺(tái)的基本內(nèi)容和K均值聚類算法的基本原理。在并行步驟中，基于以往研究的不足之處，提出了對(duì)初始值規(guī)約計(jì)算的并行方法進(jìn)行的討論和選擇、數(shù)據(jù)分配時(shí)的并行策略選擇和計(jì)算中心點(diǎn)的并行滑動(dòng)門方式，這些方法提高了GPU的利用率，降低了數(shù)據(jù)獲取延遲，充分發(fā)揮了CUDA架構(gòu)的并行計(jì)算優(yōu)勢(shì)。通過與傳統(tǒng)規(guī)約并行方法相比，滑動(dòng)門并行計(jì)算算法擁有更好的加速比。通過實(shí)驗(yàn)證明了CUDA在并行計(jì)算方面的有效性，結(jié)果表明了在GeForce730上實(shí)現(xiàn)了算法2.8倍的加速比，有效地提高了K均值算法在計(jì)算機(jī)上的運(yùn)行效率。

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

[2] 陶 冶,曾志勇,余建坤,等.并行k均值聚類算法的完備性證明與實(shí)現(xiàn)[J].計(jì)算機(jī)工程, 2010, 36 (22) :72-74.

[3] 喻金平,鄭 杰,梅宏標(biāo),等. 基于改進(jìn)人工蜂群算法的K均值聚類算法[J]. 計(jì)算機(jī)應(yīng)用,2014,34(4):1065-1069,1088.

[4] 霍迎秋，秦仁波，刑彩燕，等.基于CUDA的并行K-means聚類圖像分割算法優(yōu)化[J].農(nóng)業(yè)機(jī)械學(xué)報(bào),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.