一種基于CUDA的K-Means多級并行優化方法

方玉玲，那麗春

(上海立信會計金融學院 信息管理學院，上海 201209)

1 引 言

作為機器學習中主要算法之一，聚類分析在多個領域得到充分應用，如數據挖掘，大數據分析，推薦系統等[1，2].通過聚類算法不僅能夠對用戶關注的類別進行區分，還能用于進行個性化推薦，如淘寶、微博的智能興趣推薦中最常使用的就是聚類技術.從使用效果上看，聚類分析也是數據分類的一種，但與分類技術還存在一定的差別，最大不同之處在于聚類處理數據的所屬類是未知的，它是一個無監督過程[3].是在沒有相關經驗的基礎上，對數據進行處理，分析出數據間內在關聯并找出規律，將樣本點間距離較近的數據分配到相同的聚類中，反之，將樣本點間距離較遠的數據分配到不同的聚類中.目前常用的聚類算法分別是基于密度、網格、層次、劃分和模型的聚類算法[4].

傳統K-Means算法主要是針對大型數據庫而生的，一般情況下計算執行得到局部最優解.目前串行K-Means算法及其優化方法的發展已相當成熟[5-7]，然而當需要處理的數據量越來越多時，硬件對K-Means聚類的計算性能的限制也愈加凸顯.因此，隨著各種并行框架的發展，對并行K-Means算法的研究和應用也越來越廣泛.文獻[4,7，8]都是基于MapReduce的K-Means并行算法，均使用了減少流量的思想加速程序執行.文獻[9]提出基于Spark框架的K-Means優化方法，利用網格來存儲樣本點的坐標信息降低冗余計算，并利用Spark框架對算法實現并行化.文獻[10-12]是基于CUDA的K-Means優化方法，對CUDA編程模型，OpenMP以及MPI三種優化方式進行了比較.文獻[10]指出使用GPU編程比CPU編程性能提升了68倍，證明了CUDA編程對于計算密集型算法的優勢.文獻[11]評估K-Means數據集在不同初始化情況下對性能的影響且側重于對初始化數據的選擇，其通過實驗對比證明CUDA編程模型對于大圖像或者大數據集處理時可獲得最佳效果.文獻[12]從GPU存儲資源限制和CPU-GPU數據傳輸時間對并行程序進行優化來提升并行K-Means算法執行性能.

上述文獻雖然使用了CUDA并行優化但均只停留在線程級并行的程度，并沒有探索繼續進行性能優化的可能.因此，本文提出一種新的基于CUDA的多級并行的K-Means性能優化方法(multi-level parallelism optimization method，MLPOM).該方法不僅考慮了K-Means算法執行時的存儲訪問時延，也充分考慮了程序執行時的算術操作時延，結合GPU的多線程，線程束調度機制分別從線程塊級，線程級，指令級，比特級進行分析以充分發揮GPU良好的運算能力.本文的優化方法與文獻[11，12]的兩種方法(方法1和方法2)進行比較，極大的改善了K-Means的運行性能并降低了資源占用率.

2 相關工作

2.1 K-Means算法

作為目前使用最廣泛的聚類算法之一，K-Means算法[13]被廣泛應用于多個領域.其在未知類別的情況下分析數據對象，且在一個聚類中的樣本點之間具有更多的相似性.現有需分類數據集C={c1，c2，…，cn}共有n個樣本點，其中每個點的數據維度為d，即cj=(cj1，cj2，…，cjd)，K-Means算法的目的是把上述數據集分為確定的k類.要找到以上問題的最優解需要遍歷所有可能的類劃分，具體步驟通過串行算法1進行描述.

算法1.基于CPU的串行K-Means聚類算法

輸入：劃分聚簇的數目k，包含n個樣本的數據集

輸出：滿足條件的均值表，k個聚簇中心

1.在樣本中隨機選取k個樣本點作為k個簇的中心點；

2.計算所有樣本點與各個簇中心之間的距離，然后把樣本點劃入最近的簇中；

3.根據簇中已有的樣本點，重新計算簇中心；

4.重復2、3直至樣本點到簇中心距離最小.

該算法的任務是將數據聚類成k個簇，其中ui為簇ci的中心點：

已知數據維度為d，因此，當數據集越來越大或者維度越來越高時，算法的計算復雜度也隨之升高，運用CUDA并行化算法通常能取得更好的加速效果.

2.2 CUDA編程模型及其架構

CUDA(Compute Unified Device Architecture)是由顯卡廠商NVIDIA推出的運算平臺，是一種通用的并行計算架構，其能充分結合CPU和GPU的優點.其中，CPU為主處理器(host)執行邏輯事務處理和串行計算，GPU作為協處理器(device)執行高度并行化的計算任務[14].GPU的計算核心均勻劃分到多個流多處理器(stream multiprocessor，SM)中.同時，GPU還有不同的存儲機制.其中，設備內存可接受來自CPU的數據；共享內存可供SM中的所有線程塊block公用；寄存器被分配的thread單獨使用；其他存儲器，如紋理內存，各級緩存在本文中不做詳細討論.

在CUDA程序執行過程中，以_global_為前綴運行在GPU上的部分稱為核函數(kernel)，且這段代碼在CPU上是全局可見的，核函數的類型，數量，功能是隨程序變化的.通過kernel_func<<>>(param1，param2，…)的形式來進行kernel的調用，圓括號()中表示kernel調用時的所有參數，num_blocks和num_threads分別代表執行kernel的block數和thread數，所有thread分布在不同的block中，所有的block又分布在不同的SM中.一個完整的并行程序可能包含多個并行kernel和串行處理步驟.

在CUDA架構中，SM中基本調度單元叫做線程束(warp)，每個warp包含同一block的32個thread.自Kepler架構開始，每個SM中有4個線程束調度器和8個指令調度單元，如圖1(a)和圖1(b)所示，可以同時issue/execute這些線程束.若SM中warp個數少于4個，則意味著一些指令調度單元處于閑置狀態，會對程序性能產生影響.

圖1 GTX 680 block diagramFig.1 GTX 680 block diagram

3 多級并行的Kmeans優化算法

由K-Means串行算法可知，程序中主要計算部分集中在樣本點到樣本聚類中心的距離以及中心坐標的更新兩部分.因此，本文的并行優化也將集中于上述計算過程.以下4個小節將對本文優化過程進行詳細介紹.

3.1 K-Means算法并行化

為更好地提升計算性能，K-Means算法中的距離計算可利用NVIDIA通用函數庫中矩陣乘的思想來進行.需要注意的是，與傳統矩陣按行優先存儲方式不同的是，并行算法中的矩陣都是依照列優先的順序存儲的.因此，首先要實現矩陣轉置的并行化，經轉置后的矩陣按照列優先的方式存儲，更方便匹配CUDA架構中的計算.矩陣轉置是通過對角元素的位置互換實現的，假設矩陣A為m行n列，滿足A=a(i，j)mxn，則A的轉置矩陣B可記為B=b(j，i)nxm，且為n行m列.其并行實現的主要代碼如下：

__global__void invert_mapping(int A[]，int B[]) {

int i=blockIdx.x*blockDim.x+threadIdx.x；

int j=blockIdx.y*blockDim.y+threadIdx.y；

if(i>=n‖j>=m) return；

B[j+i*m]=A[i+j*n]；

}

其中，blockIdx.x為線程塊X方向的索引，blockIdx.y為線程塊Y方向的索引，blockDim.x為該線程塊X方向上的線程數量，blockDim.y為該線程塊Y方向上的線程數量，threadIdx.x為線程塊X方向上的線程索引，threadIdx.y為線程塊Y方向上的線程索引.

并行的第2步需要計算各樣本點與當前類中心的距離，其主要代碼如下：

__host__…float cluster_dist_calc(int d，int n，int k，float*objects，float*clusters，int objectId，int clusterId) {

int i；

float distance=0；

for(i=0；i

distance+=(objects[n*i+objectId]-clusters[k*i+clusterId])*(objects[n*i+objectId]-clusters[k*i+clusterId])；

}

return distance；

}

上述代碼用于計算樣本點objects[objectId]與中心點clusters[clusterId]之間的距離.其中，d表示樣本點維度，n表示樣本點個數，k表示聚類中心數，objects表示所有樣本點，維度為n x d，clusters表示聚類中心.

并行第3步要比較并對聚類中心坐標進行更新，當第objectId個樣本點距某一簇類中心點clusterId的距離變小時就更新該樣本點的聚類中心，否則保持其聚類中心不變.對所有樣本點都要進行聚類中心的更新.其并行實現過程稍微復雜，故在此不再呈現.綜合上述過程，本文并行算法完整步驟如下：

算法2.基于GPU的并行K-Means算法

輸入：劃分聚簇的數目 ，包含n個樣本的數據集

輸出：滿足條件的均值表，k個聚簇中心

1.在樣本中隨機選取k個樣本點作為各個簇的中心點；

2.并行化實現矩陣轉置，核函數為invert_mapping；

3.并行化計算每個數據點到多個聚類中心的距離，并樣本點劃入最近的簇中，核函數為cluster_dist_calc；

4.并行化實現聚類中心的更新，核函數為update_cluster_center；

5.重復2、3至樣本點和簇中心距離最小.

上述算法主要描述了并行計算過程，對數據的初始化及CPU與GPU間傳輸并未涉及.

3.2 多線程并行

GPU憑借其多線程并行(multi-thread parallelism，MTP)技術盡可能的隱藏算術計算和內存訪問的處理時延來提高執行性能[15].若單個warp因所執行數據之間存在依賴關系或存儲時延而暫停，則warp調度器將在warp池發出另一個活躍warp，來減少時延.暫停覆蓋的可用性依賴于SM中符合條件的warp的數量，這也是GPU需要大量同步thread的主要原因[16].

傳統思想認為，block中的thread越多越好，因為越多的thread可以執行更多的算術運算，且設備的資源占用率(occupancy)也更[17].然而，在實際應用中并非如此，有時受活躍warp數量限制，GPU隱藏計算延遲的能力隨之降低.尤其是當warp仍在進行全局內存訪問時，反而并行度越低越易獲得高性能[18].已知較高的MTP并不總是意味著更高的性能，但是較低的MTP總是缺乏覆蓋內存延遲的能力.因此，在第一級并行優化中，為了量化SM中活躍warp的比例并分配最合適的MTP，我們應該遵循相關規則.

3.2.1 優化線程配置

本文主要以Maxwell GTX970為例進行敘述，其共有13個SM，每個SM中有128個CUDA cores (SPs)，具體參數如表1所示.

表1 基本參數Table 1 Basic parameters

在CUDA架構中，為方便SM進行基本的調度，每個block中的thread個數應該是warp能容納的thread數(32)的倍數.且每個線程網格grid中block的維度以及thread的維度并不會影響程序性能[16].因此，本文只需考慮基本的block及thread配置情況：

1.自定義block數Nblock和每個block中的實際thread數Nthread_block；

2.保證GPU中所有SM均處于活躍狀態(既不存在閑置SM)，并計算每個SM中的block數Nblock_SM，參見公式(1)；

3.每個SM中要有充足的活躍block以避免部分block因等待__syncthreads()指令而使硬件閑置的情況；

3.2.2 共享內存使用

在GPU中，共享內存是一塊可讀/寫的存儲資源，能夠被同一block中的所有thread進行訪問.使用共享內存能得到thread間最小的通信延遲.同時，在使用共享內存是要注意只有當指令載入寄存器后或thread間有共享數據時才會開始使用該存儲器.因此，在使用共享內存時要滿足如下條件：

每個block中的共享內存Mshared_block資源要按當前SM中block的數量進行分配，參見公式(3)；

通過對共享內存的合理分配，進一步細化了GPU中計算資源和存儲資源分配.

3.2.3 寄存器分布

寄存器文件是GPU上速度最快的存儲空間，也是程序性能提升的最佳選擇.為了更好地使用并行thread，可在kernel中將有很多操作的任務分成許多獨立的部分，以便使用一組寄存器來計算獨立的問題，并最大化寄存器使用.然而，受寄存器總量限制，每次可用的寄存器文件的數量也是有限的，因此當單個block需要使用過多的寄存器時SM上活躍的block數量將會減少.

在沒有寄存器溢出的情況下，每個SM中寄存器的最大數量為64K.適當減少每個thread中的私有變量可以減少單線程中寄存器使用量.因此，在實際使用中，需滿足如下條件：

1.每個block中的thread數Nthread_block與thread使用的register files數成反比Nreg_thread，參見公式(4)；

2.受register files數量限制，每個thread中register files數量越多意味著更少的thread數；

3.每個thread中最多可以擁有63個register files.

綜合前述對基本線程配置，共享內存及寄存器的使用情況，實際編程時核函數需要同時滿足的表達式如下：

(1)

(2)

(3)

Nreg_block=Nreg_thread×Nthread_block

(4)

(5)

Nreg_thread≤63

(6)

Minthread_SM≤Nthread_block≤Maxthread_SM

(7)

(8)

其中，Nblock_SM代表SM中block數，Nblock_SM在不同的SM中可能并不相同，因為它們是通過輪詢方式分配到SM中的.然而，無論block和thread的數量如何變化，都必須保證當前kernel中的總體并行度Pdegree前后一致.

當kernel滿足表1及公式(1)-公式(8)時，可得到一組或多組<<>>的配置組合，通過實驗對比保留性能較好的一組.若要繼續提升程序性能，則可以考慮引入多指令并行(multi-instruction parallelism，MIP)，即使用單個thread來處理數據集的多個元素.

3.3 多指令并行

在進行CUDA編程時，人們往往忽略了MIP的重要性，而MIP能進一步降低程序占用率occupancy而提升計算性能[19，20].MIP旨在增加每個活躍thread中的執行操作數和待操作數據，并以此來降低所需執行thread數以進一步提升程序性能.而且，它能夠很大程度上降低全局存儲器的數據采集操作.

因此，在任務量不變情況下(滿足公式(8))，使用MIP使thread中指令數增加則SM中block數減少，那么可以為每個block分配的共享內存，寄存器文件，緩存的比例都會相應增多.甚至在某些情況下，它可能會完全隱藏延遲.本文實現MIP的方法如下.

3.3.1 循環展開

循環展開能夠實現在運行一個循環的時間開銷里完成展開之前多個循環所執行的數據操作.如圖2所示.

圖2 循環展開對比圖Fig.2 Loop unrolling comparison chart

已知循環成本包括：循環變量增加，循環終止判斷，迭代分支判斷和每次計算內存地址訪問開銷.由上圖可知，使用循環展開后，kernel執行時間大大減少而程序的并行任務總量保持不變.

3.3.2 寄存器重分配

CUDA架構支持只當變量被使用時才對其進行訪問.因此，可以在變量被使用前才進行賦值減少寄存器用量.

一般地，我們在kernel運行前定義多個(此處僅以3個為例)變量i，j和k，但實際上，不同的變量在程序中并不同時被使用.因此，在實際編程中可以通過改變變量的定義順序來減少寄存器用量，即只當某個變量被使用前才對其進行聲明和賦值來減少同一時刻寄存器的使用個數.這樣，通過簡單的單個寄存器重分配就能實現增加寄存器總量的效果.

3.4 多比特并行

使用MTP和MIP技術帶來了新的性能提升并降低了GPU資源占用率(在實驗結果中體現)，但研究者對于提升性能的研究是永不止步的，經多次實驗發現使用多比特并行(multi-bit parallelism，MBP)能進一步通過降低延遲的方式來提升性能.

CUDA架構支持使用向量類型變量vector_N.其中，vector_N是包含N個基本類型元素的結構.例如int2包含了2個int類型的元素，float4則包含了4個float類型的元素.

根據實際需要，也可以自定義新的向量類型，如uint4，uchar4.通過使用vector_N可以增加單個變量所含元素的數量.這一過程也引進了一定量的比特級并行，如int2和int4 就分別引入了64比特和128位的隱式并行.因MBP的效果與MIP效果相似，在后續試驗中把MBP與MIP放在一起進行分析.

3.5 優化算法完整流程

綜合前述各小節內容，本文多級并行的K-Means優化算法完整流程如圖3所示.

圖3 K-Means優化算法流程圖Fig.3 Flow chart of K-Means optimization algorithm

其中，并行化程序包括3個核函數，分別用于實現矩陣轉置，樣本點到已知聚類中心距離，更新聚類中心.多線程并行包括調整線程配置，規劃共享內存，寄存器分配3個部分；多指令并行包括循環展開和寄存器重分配兩部分；多比特并行通過使用向量類型變量來實現.

4 程序執行模擬

根據優化后的程序流程圖模擬出的kernel在不同MIP時的線程結構，如圖4所示.其中，圖4(a)給出了MIP1(此時MIP=1)的thread間執行順序，為了方便敘述，我們假設每個block有4個warp.其按照輪詢的方式從每個warp中發出多個instructions后，warp命中global memory訪問停頓.圖4(b)給出了MIP2的執行仿真情況，其中MIP2=2*MIP1，因此warp中每個thread的指令數是圖4(a)的兩倍.為了保持不同MIP情況下并行度總數一致，圖4(b)中每個線程塊的thread數量自動減半.

圖4 線程間執行順序模擬Fig.4 Execution sequential simulation in block

使用MTPOM，通過減少每個block中的warp數量來減少長延遲，如圖4(b)所示，它與圖4(a)相比擁有更多的內存吞吐量，當執行相同數量的任務時，它可以利用較少的warp(thread)而使每個thread同時擁有更多的指令.因此，本文方法可以更好的隱藏延遲等待.線程內執行情況仿真如圖5所示，展示了每個時鐘周期cycle內單個warp在不同MIP的情況下發布指令數情況.

圖5 線程內執行順序模擬Fig.5 Execution sequence simulation in thread

當MIP=1時表示單個thread/cycle只能發出一條instruction，則一個warp中總共有32條instructions；MIP=2表示單個thread/cycle可以同時發出兩條instructions，則一個warp中共有64條instructions可執行；當MIP=4時表示一個warp/cycle發出128條instructions.如果有充足的計算資源和存儲資源，則每個cycle的instruction越多，執行相同任務所需的thread數就越少.當然，如果我們使用不同程度的MIP，它同時發布的instruction數以及需要的數據都是相互獨立的.

5 實驗與結果分析

5.1 實驗環境

本文主要實驗平臺的CPU為Intel Core i7-7700 CPU，主頻3.60GHz，內存16GB，GPU為GTX970.所有程序在Windows 10系統下的Visual Studio 2015和CUDA 8.0環境下運行.

5.2 結果分析

在實驗環境和樣本數據不變的情況下，分別對不同優化程度的并行算法進行性能驗證.不失一般性，本文所有結果均經過多次實驗取其平均統計值進行分析，實驗設計同時考慮了下列因素：

1)為了驗證本文算法帶來的性能優勢，1分別對方法1、方法2及本文的MLPOM方法展開多組數據的性能測試，測試結果如表2所示；

2)為了分析MLPOM不同優化條件下的性能變化，分別在對應的實驗環境中測試多組自變量.MTP的取值分別設置為32的倍數，包括32，64，128，256，…，1024，MIP的值分別為1，2，4，8，…，32.分析不同自變量情況下算法性能變化及可能的原因，同時對不同優化程度下的warp使用率進行分析.

5.2.1 不同并行方法結果對比

對文獻中的方法1、方法2和我們的并行方法MLPOM進行對比，實驗結果見表2.

表2 不同并行K-Means執行時間Table 2 Execution time of different parallel K-Means

本實驗所用數據集是包含10個聚類中心的二維樣本點，表2對3種不同優化方法在不同數據量下的運行結果進行比較，數據集單位為MB，運行時間單位為s(秒).從表2中可以得出本文的K-Means并行優化效果最好，且與前兩種方法中的較好方法(方法2)相比最高加速比達到了39.7%(數據量為2.38M時)，在數據量繼續增大時優化效果有所降低，但都保持在12.3%以上，表2中平均優化效果為22.15%.當數據量大于2.38M后本文優化方法的加速比雖受CPU-GPU數據傳輸的影響而下降，但多級并行的優勢仍然帶來較大的性能提升.

5.2.2 本文方法分析

對本文MLPOM不同優化程度的MTP和MIP組合進行了詳細的探索，并比較不同組合下的性能變化.為了方便顯示，本文把所有性能結果相對于基準性能(未使用本文方法優化前的性能)進行了歸一化.

圖6顯示了K-Means在不同MTP情況下對GPU的資源占用率occupancy及性能變化情況.首先可以看出occupancy并非隨著thread增多而增大，而是隨thread增多而降低.這是因為當K-Means在單個block中的thread增多時，每個block對寄存器的使用也相應增多，會導致block數減少，根據公式(4)和公式(10)可知SM中的總thread數反而下降，此時，GPU的occupancy變小.同時，圖6表明本文K-Means算法在thread數為128時取得性能峰值.圖7進一步對不同MIP時的性能和資源占用率進行分析.

圖6 不同MTP的占用率和性能Fig.6 Occupancy and performance of different MTPs

圖7 不同MIP的占用率和性能圖Fig.7 Occupancy and performance of different MIPs

由圖7可知，GPU的occupancy隨MIP增加而降低.同時，在單thread指令數為8時取得性能峰值，其較指令數為1時性能提升了近40%，實驗結果表明并非MIP越高算法性能越好.

為了進一步探究本文MTPOM提高性能的原因，我們對最優化K-Means在GTX 680(計算能力3.0)上的運行過程進行深入分析.圖8顯示了每個時鐘周期發布的相關指令(IPC)及對應不同MTP和MIP的情況，其中，IPC表示每個SM的指令吞吐量，其理論峰值與設備的計算能力相關.

圖8 IPC分析Fig.8 IPC Analysis

隨著MIP(MBP)的改變，程序發布的IPC(IPC issued)和執行的IPC(IPC executed)也隨之變化.當發布越多IPC時，每個cycle能夠發出更多指令使時延變短.執行較高的IPC意味著資源使用效率更高程序性能更好.從圖8中可以看出，K-Means在單個線程中有8個獨立指令時，每個cycle可以發出最大指令，因此能夠達到最佳性能.這一結果也說明性能改進的另一個重要因素是warp發布效率(kernel執行期間在所有周期內具有合格warp的周期百分比)的變化.

圖9顯示了最佳性能情況下K-Means的warp發布情況及導致warp暫停的不同原因所占比例，并給出了K-Means中活躍warp不符合條件的原因，有多種原因導致warp調度程序沒有合適的warp可供選擇.很明顯warp暫停原因隨著MIP的變化而變化.首先，當MIP=8時，它具有最高的warp issue效率35.03%(從副坐標可以看出).其次，執行依賴性和內存依賴性的百分比在MIP=8時最低，此時，它們對性能的影響最小.

圖9 K-Means線程束指令分布Fig.9 Warp instruction distribution

6 結 語

為了提升并行K-Means程序的執行性能，本文提出了基于CUDA編程架構的多級并行優化模型(Multi-level parallelism optimization model，MLPOM).該模型充分挖掘GPU計算資源和存儲資源未被合理充分利用的事實，使用MTP，MIP和MBP多級組合并行來重新分配K-Means并行算法中核函數對各級存儲器和計算資源的使用情況.實驗結果表明，在多thread并行時，并非block中thread數越多程序性能越高，而且GPU占用率也不總是和thread數呈正比.在多指令和多比特級并行中，K-Means在單thread指令數為8時獲得最佳性能.這是因為本文優化模型不僅通過MTP改善了計算以及存儲資源的利用率，而且通過使用MIP(MBP)增加了單cycle內運行的指令數，使得相同數據量下執行時間減少，達到提升算法性能的目的.在未來的工作中，在繼續優化CPU-GPU間數據傳輸性能的同時，會把本文優化思想形成標準化模型推廣到更多的并行算法的性能優化中.