利用CUDA提高內(nèi)存數(shù)據(jù)聚類效能的研究

董麗麗，董 瑋，張 翔

西安建筑科技大學(xué)，西安 710055

1 引言

聚類分析通常用于將數(shù)據(jù)劃分成有用的或有意義的組或簇，在心理學(xué)、生物學(xué)、統(tǒng)計學(xué)、信息檢索、機(jī)器學(xué)習(xí)和數(shù)據(jù)挖掘中扮演著重要的作用。在心理學(xué)領(lǐng)域，特定的人群存在著相應(yīng)的心理學(xué)特征，將此類人群提取特征值并進(jìn)行聚類，可以幫助心理學(xué)專家對此類人群的心理特征進(jìn)行分析。K-Means算法作為傳統(tǒng)的數(shù)據(jù)挖掘方法，在心理學(xué)數(shù)據(jù)的聚類分析中具有重要作用。但K-Means算法存在初始聚類中心敏感問題，導(dǎo)致聚類有效性降低。針對K-Means算法的中心敏感問題，Bradley P S等學(xué)者提出[1]，在聚類分析的數(shù)據(jù)中隨機(jī)選取多個子樣本，對選取的子樣本重復(fù)執(zhí)行K-Means算法得到優(yōu)化的聚類中心，再對整個數(shù)據(jù)集進(jìn)行K-Means算法聚類。然而，此方法受到子樣本集選取方式的影響，只是一種次優(yōu)結(jié)果。將智能算法與傳統(tǒng)數(shù)據(jù)挖掘方法相結(jié)合的方式在近些年來廣泛被研究[2-6]，可有效提高聚類有效性，但在實(shí)際應(yīng)用當(dāng)中，執(zhí)行效率較低。針對傳統(tǒng)數(shù)據(jù)挖掘算法的效率問題，將Mapreduce編程模型運(yùn)用到提高大規(guī)模數(shù)據(jù)聚類效率方面已經(jīng)成為熱點(diǎn)，薛勝軍等學(xué)者利用開源組織apache基金會下的hadoop項(xiàng)目[7]，實(shí)現(xiàn)了K-Means算法的并行化，提高了效率。但在實(shí)際運(yùn)用當(dāng)中，心理學(xué)數(shù)據(jù)不具備網(wǎng)絡(luò)數(shù)據(jù)的即時性、高訪問性、交互性和大規(guī)模性，利用云平臺進(jìn)行此類數(shù)據(jù)挖掘并不適合，并且由于網(wǎng)絡(luò)的延遲性和不可靠性，心理學(xué)等內(nèi)存數(shù)據(jù)在利用hadoop平臺實(shí)際挖掘當(dāng)中，并不比傳統(tǒng)方法在執(zhí)行效率上好多少。

CUDA技術(shù)是NVIDIA公司針對其GPU（Graphic Processing Unit）產(chǎn)品所發(fā)明的一項(xiàng)高性能計算技術(shù)，此技術(shù)非常適合心理學(xué)等數(shù)據(jù)規(guī)模較小的內(nèi)存數(shù)據(jù)挖掘，此項(xiàng)技術(shù)為在有限時間內(nèi)多次聚類，并從多次聚類結(jié)果中觀察最優(yōu)結(jié)果提供了可能。Kumar等學(xué)者[8]利用CUDA實(shí)現(xiàn)了高斯混合模型期望最大化算法的并行化，加速比為165倍；張繁等學(xué)者[9]將CUDA應(yīng)用于計算蛋白質(zhì)分子長，對比CPU加速比為7.5倍；Takahashi等學(xué)者[10]利用CUDA加速了三維亥姆霍茲方程，加速比為23倍；Labaki等學(xué)者[11]把CUDA應(yīng)用到二維勢問題和二維彈性靜力學(xué)問題，加速比分別為56.8和39.6倍。

為了解決上述問題，本文提出了一種基于CUDA技術(shù)和FP-Tree關(guān)聯(lián)迭代相結(jié)合的K-Means聚類算法—AIK-Means（Associated IterationK-means）算法，利用CUDA解決多次聚類的執(zhí)行效率問題。利用Chameleon層次聚類算法找出K-Means算法的初始中心。Chameleon算法是一種基于層次的聚類圖算法，該算法利用基于圖的方式得到初始數(shù)據(jù)的劃分，使用簇間的接近性和互連性概念和簇的局部建模，發(fā)現(xiàn)具有高質(zhì)量的不同形狀、大小和密度的簇[12]。所以本文首先對內(nèi)存數(shù)據(jù)進(jìn)行初步聚類，得到K-Means算法的初始參數(shù)，從而利用Chameleon算法聚類緩解K-means初始聚類中心的隨機(jī)性和先驗(yàn)性導(dǎo)致聚類結(jié)果不穩(wěn)定的問題。利用FP-Tree對多次聚類結(jié)果進(jìn)行關(guān)聯(lián)分析，從K-Means算法的執(zhí)行結(jié)果入手，降低K-Means算法的中心敏感問題，提高聚類有效性。

2AIK-Means算法

2.1 AIK-Means算法框架

AIK-Means算法主要流程如圖1所示。

數(shù)據(jù)流在顯卡設(shè)備（GPU）中利用CUDA技術(shù)加速計算，在主機(jī)設(shè)備（CPU）中利用.NET Framework平臺進(jìn)行IO讀寫和任務(wù)處理。AIK-Means算法思想如下。第一步：將心理學(xué)數(shù)據(jù)預(yù)處理，以達(dá)到層次聚類算法Chameleon所要求的格式。第二步：利用Chameleon算法進(jìn)行初始聚類，為下一步K-Means算法找出K個初始簇。第三步：將初始簇利用K-Means算法進(jìn)一步深層聚類，找出K個初始簇的收斂中心點(diǎn)。第四步：重復(fù)第二步、第三步n次，得到n個聚類結(jié)果。第五步：利用FP-Tree對n次聚類結(jié)果進(jìn)行關(guān)聯(lián)分析，在多次聚類結(jié)果中找出最終的穩(wěn)定結(jié)果。

2.2 心理學(xué)數(shù)據(jù)的預(yù)處理

圖1 算法流程圖

本文利用行為心理學(xué)中的MMPI（Minnesota Multiphasic Personality Inventory）量表集，采集某集團(tuán)單位人群的行為心理學(xué)數(shù)據(jù)，對人群進(jìn)行聚類，幫助心理學(xué)研究者分析人群的多項(xiàng)人格。MMPI是行為心理學(xué)最為復(fù)雜和繁瑣的一種測試，但其準(zhǔn)確度與其他量表相比較為準(zhǔn)確[13]。MMPI的第一步為獲得測試數(shù)據(jù)，計算每個測試對象的25個測試項(xiàng)目的原始分，分為兩種量表，共566道測試題。由于每項(xiàng)測試量表題目數(shù)量不同，導(dǎo)致原始分總分分值不同，無法直接比較，需要將各個量表的原始分?jǐn)?shù)進(jìn)行歸一化操作。第二步就需要通過公式（1）：

計算出各個量表歸一化為標(biāo)準(zhǔn)分T，其中T為向量，向量中每一分量代表各個量表的標(biāo)準(zhǔn)分。X為各個測試量表的原始分組成的原始分向量，M與SD分別為各個測試量表正常組原始分?jǐn)?shù)的期望及標(biāo)準(zhǔn)差所組成的向量。MMPI數(shù)據(jù)預(yù)處理的重點(diǎn)在于將原始分轉(zhuǎn)換為標(biāo)準(zhǔn)T分，求其期望M和標(biāo)準(zhǔn)差SD。

在主機(jī)設(shè)備中計算期望M和標(biāo)準(zhǔn)差SD需要多次迭代，此時迭代規(guī)模隨著數(shù)據(jù)規(guī)模的增長而增加。可利用CUDA技術(shù)在同樣的數(shù)據(jù)規(guī)模中將迭代執(zhí)行轉(zhuǎn)換為并行執(zhí)行，從而減少迭代規(guī)模，提高執(zhí)行效率。在CUDA并行環(huán)境中，可利用縮減樹算法求出期望M和標(biāo)準(zhǔn)差SD，如圖2所示。

圖2 縮減樹算法

圖2中為前13個被測試者的某一項(xiàng)測試量表原始分通過縮減樹算法并行相加，通過5次迭代，將相加結(jié)果存入共享數(shù)組首元素中，此結(jié)果除以數(shù)組長度，即為此項(xiàng)測試量表的期望M。具體步驟如算法1所示。

算法1縮減樹算法

函數(shù)說明：函數(shù)Expectation為在GPU上運(yùn)行的核函數(shù)，本文算法函數(shù)前若帶Kernel function關(guān)鍵字皆為核函數(shù)。輸出形參R為在CUDA上各個線程塊之和的數(shù)組，將R拷貝到主機(jī)內(nèi)存中，利用CPU求公式（2）

在算法1中，threadIdx、blockIdx、gridDim、blockDim和shared均為CUDA內(nèi)置變量。threadIdx為線程塊內(nèi)線程標(biāo)號；blockIdx為線程格內(nèi)線程塊標(biāo)號；gridDim為線程格尺寸，通過此變量可計算出線程塊的數(shù)量；blockDim為線程塊尺寸，通過此變量可計算線程塊內(nèi)線程數(shù)量；shared為GPU塊內(nèi)共享內(nèi)存數(shù)組；CUDA內(nèi)置函數(shù)synchronize threads為同步函數(shù)，使所有線程在當(dāng)下時刻同步。對于求SD函數(shù)Kernel function StandardDeviation(T,M,R)，利用函數(shù)Expectation求得期望M，作為其第二實(shí)參傳入。將函數(shù)Expectation第4 行shared[tid]←shared[tid]+T[i]改 為shared[tid]←shared[tid]+(T[i]-M)2，即為函數(shù)StandardDeviation，利用縮減樹算法的結(jié)果開平方即可求得SD。

2.3 利用CUDA改進(jìn)Chameleon算法

Chameleon算法是一種兩階段聚類層次算法，第一階段采用K最鄰近法把點(diǎn)分成多個小簇，把一個點(diǎn)和其最近鄰的k個點(diǎn)連接起來；第二階段根據(jù)相近程度合并多個小簇，計算任意兩個簇的互聯(lián)性RI和緊密性RC，當(dāng)兩個指標(biāo)都較大時合并兩個簇。RI的公式（3）與RC的公式（4）如下所示。

計算出在下三角陣H中的位置(i,j)，其中，k為數(shù)組R下標(biāo)，m為當(dāng)前矩陣的行數(shù)，當(dāng)k及m使不等式成立時，i=m，j=k-m(m-1)/2。利用(i,j)找出向量數(shù)組R中的向量Ri與向量Rj，求出兩向量間的對稱度量距離，存入數(shù)組H下標(biāo)為k的元素中。在每個CUDA線程中，通過上述過程計算得到向量數(shù)組R中Ri與Rj向量的距離，存入數(shù)組元素Hk中。如圖3所示。

圖3 計算向量間距離

將數(shù)組H拷貝回主機(jī)存儲器，遍歷數(shù)組H，采用K最近鄰法構(gòu)造多個小簇，每個簇的元素數(shù)量不大于K，并計算小簇i的信息如點(diǎn)的個數(shù)|Ci|，簇內(nèi)所有邊權(quán)重和EC(Ci)，存入數(shù)組R中。具體步驟如算法2所示。

算法2計算向量間距離偽代碼

Chameleon算法的第二階段的第一步與核函數(shù)CalculateDistance執(zhí)行過程相似，R由存儲向量組改為存儲小簇組，H數(shù)組記錄小簇之間RI×RC的距離。第二步，將數(shù)組R拷貝回主機(jī)存儲器，設(shè)定閾值u，遍歷數(shù)組H，當(dāng)元素i中的RI×RC＞u時，合并i所代表的兩個小簇，將新簇組存入數(shù)組R中并且拷貝到設(shè)備存儲器，重復(fù)第一步、第二步過程，直到條件RI×RC＞u不再滿足條件為止。

2.4 利用CUDA改進(jìn)K-Means算法

圖4 計算向量組R與向量組K每對向量的距離

將矩陣H從設(shè)備存儲器拷貝回主機(jī)存儲器中，遍歷矩陣H，分別找出與K個中心向量距離最近的向量組R1～k，組成K個簇。具體步驟如算法3所示。

算法3計算中心點(diǎn)與數(shù)據(jù)點(diǎn)距離偽代碼

第二階段，將每個簇的向量組R拷貝到設(shè)備存儲器中，利用2.1節(jié)介紹的縮減樹算法計算新的中心向量，作為下一次迭代第一階段的K個簇的中心向量。重復(fù)第一階段和第二階段，直到中心向量收斂為止。

2.5 構(gòu)造FP-Tree數(shù)據(jù)結(jié)構(gòu)

將2.3節(jié)中的Chameleon算法與2.4節(jié)中的K-Means算法運(yùn)行多次，組成事務(wù)型數(shù)據(jù)，數(shù)據(jù)格式與內(nèi)容如圖5所示。

圖5 事務(wù)型數(shù)據(jù)格式（部分?jǐn)?shù)據(jù)）

圖6 FP-Tree示意圖

Root為根節(jié)點(diǎn)，Root的每一個子樹為一簇，簇中節(jié)點(diǎn)包含兩部分信息，冒號前部為被測試者編號，冒號后部為被測試者編號在某一簇的多次結(jié)果中的某個結(jié)果中出現(xiàn)次數(shù)。具體步驟如算法4所示。

算法4構(gòu)造FP-tree偽代碼

統(tǒng)計Root的每個子樹中每個被測試者出現(xiàn)的次數(shù)，當(dāng)出現(xiàn)次數(shù)大于其在整個FP-Tree樹中出現(xiàn)次數(shù)的1/2時，將每個子樹中符合條件的被測試者編號組成頻繁項(xiàng)集。具體步驟如算法5所示。

算法5構(gòu)造頻繁項(xiàng)集偽代碼

算法6 Traversal函數(shù)偽代碼

在進(jìn)行造林之后需要進(jìn)行及時的驗(yàn)收，但是在實(shí)際的驗(yàn)收工作中存在一定的問題。第一，對需要驗(yàn)收的項(xiàng)目沒有進(jìn)行及時的自查驗(yàn)收，沒有進(jìn)行建檔工作。第二，一些驗(yàn)收工作只是在表面上進(jìn)行，沒有真正落實(shí)到位。第三，驗(yàn)收圖紙和實(shí)際的情況存在一定才差距。

3 實(shí)驗(yàn)

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

實(shí)驗(yàn)環(huán)境：CPU：Intel i7-3770；CPU頻率：3.4 GHz；CPU核心數(shù)/線程數(shù)：4/8；主機(jī)內(nèi)存：8 GB；GPU：NVIDIA GTX 660；設(shè)備內(nèi)存（顯卡）：2 GB；CUDA版本：5.0；操作系統(tǒng)：Windows7。

3.2 實(shí)驗(yàn)數(shù)據(jù)

本文的實(shí)驗(yàn)數(shù)據(jù)為某集團(tuán)單位員工歷年來MMPI心理學(xué)測試數(shù)據(jù)，大小為353 MB，以行為單位，測試數(shù)據(jù)有907 689條，由5個.csv文件組成，數(shù)據(jù)格式如圖7所示。

圖7 某集團(tuán)單位員工心理測試數(shù)據(jù)

逗號將每一列分開，第一列為員工測試編號，后面每一列對應(yīng)于一道MMPI測試題，第一列為第一題，第二列為第二題，以此類推。答題中0代表肯定回答，1代表否定回答。

3.3 實(shí)驗(yàn)分析

（1）CUDA并行加速環(huán)境與純CPU串行執(zhí)行環(huán)境下AIK-Means算法的執(zhí)行效率比較

為比較算法性能，本文分別對AIK-Means算法的各個部分在CUDA執(zhí)行環(huán)境（以下表格中均簡稱并行加速環(huán)境）與純CPU執(zhí)行環(huán)境（以下表格中均簡稱串行環(huán)境）下進(jìn)行比較，利用加速比度量算法的效率。數(shù)據(jù)預(yù)處理的實(shí)驗(yàn)結(jié)果如表1所示。

表1 數(shù)據(jù)預(yù)處理實(shí)驗(yàn)結(jié)果

由表1可知，在串行環(huán)境下CPU計算MMPI中T分的時間主要用在計算其期望及標(biāo)準(zhǔn)差，在并行加速環(huán)境下把計算T分期望及標(biāo)準(zhǔn)差交給GPU，CPU主要工作變?yōu)镮O讀寫文件，將算法的預(yù)處理計算部分用于CUDA并行計算，對程序的加速比有很大提高。Chameleon部分的實(shí)驗(yàn)結(jié)果如圖8、圖9所示。

由圖8、圖9可以得出，當(dāng)RI×RC閾值遞增時，Chameleon算法的串行環(huán)境版本和并行加速版本都出現(xiàn)遞減情況。經(jīng)分析，Chameleon算法的串行環(huán)境版本不僅需要在算法的第一階段遍歷k最鄰近圖和在算法的第二階段遍歷最鄰近簇，還需要計算第一階段向量組成的下三角矩陣元素和第二階段簇組成的下三角矩陣元素。將Chameleon算法經(jīng)過CUDA并行化后，主要并行化的部分為計算其兩個階段的下三角矩陣，但遍歷k最鄰近圖及最鄰近簇依然在CPU上執(zhí)行，與串行環(huán)境算法的遍歷部分一致。隨著RI×RC閾值增大，滿足合并條件的簇越少，合并過程所付出的時間代價就越少，即串行環(huán)境版本和并行加速版本算法的遍歷部分時間代價都會減少。K-Means部分的實(shí)驗(yàn)結(jié)果如圖10、圖11所示。

圖8 Chameleon算法部分在并行與串行環(huán)境下的耗時對比

圖9 Chameleon算法部分的加速比

圖10 K-Means算法部分在并行與串行環(huán)境下的耗時對比

圖11 K-Means算法部分的加速比

由圖10、圖11可以得出結(jié)論，隨著初始簇的數(shù)量K越大，K-Means算法部分的加速比越高。經(jīng)分析，K-Means算法的兩個階段計算每個點(diǎn)與中心點(diǎn)距離和求中心點(diǎn)都進(jìn)行了CUDA并行化改進(jìn)，將算法中所有的數(shù)值計算都在GPU上運(yùn)行，CPU負(fù)責(zé)數(shù)據(jù)拷貝與調(diào)度，使得算法效率大大提高，尤其是在初始簇越多的情況下，加速比也隨之迅速提高。AIK-Means算法的總體比較如圖12、圖13所示。

從圖12、圖13上分析，RI×RC閾值的提高，增加了K-Means的初始簇個數(shù)K，造成了算法K-Means階段計算規(guī)模增加，加速比也隨之增加，由此可得出結(jié)論，當(dāng)AIK-Means算法計算規(guī)模和計算數(shù)據(jù)量增大時，算法效率更高。

（2）AIK-Means算法與K-Means算法聚類純度的比較

實(shí)驗(yàn)采用UCI機(jī)器學(xué)習(xí)數(shù)據(jù)庫的三個權(quán)威數(shù)據(jù)集進(jìn)行實(shí)驗(yàn)，分別為Iris、Wine和Glass。如表2所示。

本文采用純度標(biāo)準(zhǔn)作為聚類結(jié)果的有效性評估方法。設(shè)類簇Ci的大小為ni，則該類簇的純度定義為公式（6）：

圖12 AIK-Means算法總體在并行與串行環(huán)境下的耗時對比

圖13 AIK-Means算法總體的加速比

表2 數(shù)據(jù)及信息統(tǒng)計

其中表示Ci與第j類的交集大小。公式（7）：

為整個結(jié)果的聚類純度，k為最終形成的類簇數(shù)目，純度越高，聚類算法越有效。本文將AIK-Means算法的迭代參數(shù)設(shè)定為50次，利用FP-Tree在這50次聚類結(jié)果中進(jìn)行關(guān)聯(lián)分析，得出最終聚類結(jié)果。如圖14所示。

從圖14可以看出，AIK-Means算法比傳統(tǒng)的K-Means算法運(yùn)行結(jié)果要好。

4 總結(jié)

在內(nèi)存數(shù)據(jù)中，利用CUDA技術(shù)很好地加速了對內(nèi)存數(shù)據(jù)的分析處理效率。在對心理學(xué)MMPI數(shù)據(jù)集的聚類實(shí)驗(yàn)中，AIK-Means算法在可接受的時間內(nèi)執(zhí)行多次，隨著迭代次數(shù)的增加，算法的聚類純度也隨之提高，即執(zhí)行效率的提高有效地支持了聚類純度的提高。將AIK-Means算法應(yīng)用到某集團(tuán)心理學(xué)測試系統(tǒng)中，對單位員工的MMPI數(shù)據(jù)進(jìn)行了聚類分析，幫助心理學(xué)工作者對測試人員進(jìn)行準(zhǔn)確的劃分，達(dá)到了預(yù)期效果。

圖14 聚類純度的對比

[1]Bradley P S，F(xiàn)ayyad U M.Refining initial points forK-Means clustering[C]//Proceedings of the 15th International Conference on Machine Learning.San Francisco：Morgan Kaufmann Press，1998：91-99.

[2]Selim S Z，Alsultan K.A simulated annealing algorithms for the clustering problem pattern recognition[J].Pattern Recognition，1991，24（10）：1003-1008.

[3]馬超，練繼建.基于聚類遺傳算法的梯級水利樞紐短期電力調(diào)度優(yōu)化[J].天津大學(xué)學(xué)報，2010，43（1）：1-8.

[4]鞏敦衛(wèi)，蔣余慶，張勇，等.基于微粒群優(yōu)化聚類數(shù)目的K-均值算法[J].控制理論與應(yīng)用，2009，26（10）：1175-1179.

[5]劉靖明，韓麗川，侯立文.基于粒子群的K均值聚類算法[J].系統(tǒng)工程理論與實(shí)踐，2005，25（6）：54-58.

[6]王縱虎，劉志鏡，陳東輝，等.一種改進(jìn)的粒子群優(yōu)化快速聚類算法[J].西安電子科技大學(xué)學(xué)報：自然科學(xué)版，2012，39（5）：61-65.

[7]薛勝軍，潘吳斌.基于MapReduce的氣象數(shù)據(jù)并行PK-Means算法[J].武漢理工大學(xué)學(xué)報，2012，34（12）：139-141.

[8]Kumar N，Satoor S，Buck I.Fast parallel expectation maximization for Gaussian mixture models on GPUs using CUDA[C]//Proceedings of 11th IEEE International Conference on High Performance Computing and Communications.Piscataway：IEEE Computer Society Press，2009：103-109.

[9]Zhang Fan，Wang Zhangye，Yao Jian，et al.Accelerating the calculation of protein molecular field using GPU cluster[J].Journal of Computer-Aided Design&Computer Graphics，2010，22（3）：412-419.

[10]Takahashi T，Hamada T.GPU-accelerated boundary element method for Hemholtz’equation in three dimensions[J].International Journal for Numerical Methods in Engineering，2009，80（10）：1295-1321.

[11]Labaki J，Mesquita E，F(xiàn)erreira S L O.The BEM on General Purpose Graphics Processing Unit（GPGPU）：a study on three distinct implementations[C]//Proceedings of the 11th International Conference on Boundary Element Techniques XI.United Kingdom：EC Ltd Press，2010：316-323.

[12]杜新軍，王瑩.基于雙線性對的Chameleon簽名方案[J].軟件學(xué)報，2007，18（10）：2662-2668.

[13]聶小晶，邱昌建.焦慮癥和抑郁癥患者的MMPI對照研究[J].華西醫(yī)學(xué)，2009，24（6）：1356-1361.