周天陽,王慶林*,李榮春,梅松竹,尹尚飛,郝若晨,劉 杰
(1. 國防科技大學 計算機學院, 湖南 長沙 410073; 2. 國防科技大學 并行與分布計算全國重點實驗室, 湖南 長沙 410073)
經典的K-means算法[1]是最重要的無監督學習算法之一,被廣泛應用于數據挖掘、文檔聚類、入侵檢測等領域中。該算法包括兩個步驟:樣本點分配和中心點重選。兩個步驟不斷重復,直到每個中心點的位置不再改變或者迭代達到一定的次數。
當樣本點數或中心點數較多時,經典的K-means算法往往運行緩慢。因此,學術界提出了許多優化方法來加速算法本身,例如中心點的初始位置優化[2]、算法結構優化[3-4]和三角不等式優化[5-7]。陰陽K-means[8]是基于三角形不等式優化的最流行實現之一,其通過三角形不等式有效地避免了大多數冗余的距離計算,從而性能得到大幅提升。
與此同時,許多學者在K-means算法的并行優化方面也開展了許多工作,以期通過并行計算的方式來減少其執行時間。例如:Wu等[9]在Intel MIC架構上實現了K-means算法的細粒度的單指令多數據(single instruction/multiple data, SIMD)并行;Kwedlo等[10]提出了Drake、Elkan、Annulus和陰陽四種K-means變種算法的MPI/OpenMP混合并行化;Zhao等[11]提出了一種基于MapReduce的并行K-means算法;Kumar等[12]實現了在多核超級計算機上使用大型數據集進行定量生態區域劃定的并行K-means聚類。此外,現有研究還嘗試面向GPU[13-15]、FPGA[16-17]等異構硬件平臺進行K-means 算法的并行優化。例如,Taylor等[18]優化了陰陽K-means算法在GPU上的實現。因此,K-means算法的并行優化是學術界關注的研究熱點。
目前,眾核處理器已經成為高性能計算領域中的主流平臺。考慮功耗與性能的平衡,基于SIMD方式運行的向量單元幾乎成為眾核處理器的標配,如x86 AVX-512[19]、ARMv8 Neon[20],從而也使得向量單元成為眾核處理器計算性能提升的主要來源之一。同時,眾核處理器還通常采用非一致內存訪問(non-unified memory access, NUMA)架構,以擴展片上計算核的數量以及訪存帶寬。單個處理器的所有核心組織成多個NUMA節點,如AMD EPYC系列[21]、飛騰FT-2000+[22-24];或者多個處理器以NUMA方式組織成為一個多socket節點。
盡管現有陰陽K-means算法的開源多線程實現可以直接運行于眾核處理器上,但其性能通常不是最優的。其原因主要在于這些開源實現既沒有利用眾核處理器中的向量單元,也沒有面向NUMA架構開展針對性的性能優化。
本文主要面向眾核處理器上的陰陽K-means 算法的高性能優化,通過向量化、NUMA親和性內存訪問優化以及數據布局優化等手段實現陰陽K-means算法的高效并行,設計了多組實驗來評估優化方法的性能,并將本文提出的陰陽K-means實現與開源多線程陰陽K-means[25]實現進行性能對比。
傳統陰陽K-means算法的偽代碼如算法1所示。該算法隨機選擇k個樣本點作為中心點,并將它們分為t個中心點組,然后初始化所有n個樣本點(第1~3行)。算法的核心是使用一個三級過濾器過濾不必要進行距離計算的樣本點或中心點。具體來說,算法每輪迭代使用全局過濾器(第9~11行)過濾不必要進行距離計算的樣本點,使用組過濾器(第12~13行)和局部過濾器(第14~15行)過濾不必要進行距離計算的中心點,最后剩下的每個樣本點與所有中心點進行距離計算,根據計算結果為每個樣本點找到最近的中心點,如此反復迭代,直到迭代達到一定的次數或者所有中心點的位置不再改變。

算法1 傳統陰陽K-means算法
樣本點和中心點之間的距離使用歐幾里得距離公式來計算:
(1)


圖1 陰陽K-means算法執行過程Fig.1 Execution process of Yinyang K-means algorithm
值得一提的是,本文提出的陰陽K-means算法的并行實現中沒有考慮局部過濾器,因為Newling等的研究[26]表明局部過濾器會降低性能。
向量化也稱為數據級并行化,它依賴于眾核處理器中的向量處理單元以及對應的向量指令。以Intel Xeon Gold 6240為例,該處理器配備512位向量處理單元,在AVX-512指令集的支持下,可以利用向量乘加指令同時處理8個雙精度浮點數或者16個單精度浮點數的乘加操作,大大提高了運算的并行性,從而可有效提高應用程序的性能。
NUMA是一種由多個節點組成的內存架構。 每個節點包含一個或多個處理器核心以及本地存儲器,如圖2所示。對于單個節點來說,每個節點內部的存儲器稱為本地內存,其他所有節點的存儲器被稱為遠程內存,多個NUMA節點通過互聯模塊進行通信。在NUMA架構中,處理器核心訪問內存的時間取決于內存相對于處理器核心的位置,任何處理器核心訪問其本地內存的速度都比遠程內存快得多。

圖2 典型NUMA體系結構Fig.2 Classic NUMA architecture
由于現有的研究主要使用多線程技術來對陰陽K-means算法進行并行優化,卻沒有利用眾核處理器的向量處理單元來進一步提高算法的性能,所以對陰陽K-means算法進行向量化并行優化是非常值得研究的。向量化是并行計算中最常見的優化技術之一,對陰陽K-means算法進行向量化優化之前需要確定基于樣本點數據還是中心點數據進行向量化。前者將所有樣本點數據視為一個長向量,一次計算多個樣本點到每個中心點的距離;后者則將所有中心點數據視為一個長向量,一次計算多個中心點到一個樣本點的距離。Wu等的研究[9]表明,如果需要獲得良好的性能,所構建的長向量必須足夠長才能有效隱藏相關單元的延遲,而樣本點數量要比中心點數量大得多,從而可以為SIMD指令流構建更長的向量。因此,決定基于樣本點數據進行向量化優化??紤]到陰陽K-means算法的執行特性,不僅可以在距離計算階段引入向量化,在其他任何可以批處理的地方都可以引入向量化。例如,可以將全局過濾器條件中的標量距離界限的比較轉換為向量距離界限的比較,從而一次檢查多個樣本點的全局過濾器條件。但是必須調整數據內存布局,使其更適合向量化并有更好的數據局部性。
NUMA架構良好的擴展性以及更高的訪存帶寬使其在超級計算機或服務器上得到了廣泛應用。但是NUMA架構具有其天然的缺點,即遠程內存訪問成本高。因此,傳統陰陽K-means算法的實現直接在NUMA架構處理器上運行時可能會因為數據沒有合理分配而導致遠程內存訪問次數過多,從而可能會性能較差,當NUMA節點較多時性能下降尤其明顯。所以如何減少遠程內存訪問量是提升算法性能的關鍵。 這個問題可以通過線程綁定、數據分配、工作負載分配三個步驟來解決。
根據上述分析,為了提升陰陽K-means算法的性能,可以利用眾核處理器的向量單元進行向量化優化、針對眾核處理器的NUMA架構進行NUMA親和性內存訪問優化。此外,為了支持這兩種優化方法,還需優化內存數據布局,以獲得更好的數據局部性,以及便于在各個NUMA節點進行數據分配。
傳統陰陽K-means算法由于沒有利用眾核處理器中的向量處理單元,所以處理樣本點時單周期最多只能處理一個樣本點,而基于樣本點數據進行向量化后則可以在單周期內同時處理s(s為向量長度,且大于1)個樣本點,進而算法的性能得到有效的提升。以ARMv8架構處理器為例,由于其配備了128位的向量處理單元,所以利用向量處理單元可在單周期內同時處理2個雙精度存儲的樣本點(即向量長度s=2)。算法2展示了同時采用向量化優化和NUMA親和性優化的陰陽K-means算法的執行過程,其中多處借助了向量處理單元以提升性能。如樣本點初始化階段需要計算所有樣本點與中心點之間的距離時(第5行),一次可以同時計算兩個樣本點到一個中心點的距離。具體做法是遍歷樣本點和中心點的所有d個維度,將兩個相鄰樣本點的同一個維度數據放入一個向量中,并將單個中心點的對應維度廣播到另一個向量,最后使用SIMD指令計算兩個向量的歐幾里得距離。得益于分組聚合內存數據布局所帶來的良好的數據局部性,計算多個樣本點與一個中心點的距離所需的內存訪問都是連續的,所以非常高效。又如for循環開始后,由于單次for循環處理兩個樣本點,所以利用SIMD指令可同時檢查兩個樣本點的全局過濾器條件(第11行)。如果兩個樣本點所在簇的中心點的位置在本次迭代都不會發生變化,則可以跳過這兩個樣本點;反之,如果兩個樣本點中至少有一個樣本點通過全局過濾器,即樣本點所在簇的中心點的位置會在本次迭代中發生改變,則繼續收緊兩個樣本點的距離上限,并使用SIMD指令再次同時檢查它們的全局過濾器條件(第12~13行)。與樣本點初始化階段的距離計算過程相同,第18~23行詳細闡述了利用向量化后的樣本點數據計算兩個樣本點與通過篩選后的中心點的距離,并為兩個樣本點中通過全局過濾器的樣本點找到離它最近的中心點,從而在單次for循環中完成了兩個樣本點的分配,循環次數也因此減半,但是總的迭代次數不變。以上向量化操作由ARMv8架構提供的內部指令(ARM Neon Intrinsic)實現。

算法2 并行加速陰陽K-means算法
陰陽K-means算法在NUMA系統上運行時,所有樣本點數據將被加載到單個NUMA節點的本地內存中,導致其他NUMA節點中的處理器核心對樣本點數據的訪問都是遠程內存訪問,過多的遠程內存訪問會明顯降低算法的性能。陰陽K-means算法的執行特性表明單個NUMA節點中的各個處理器核心對樣本點數據的處理是獨立的,從而每個NUMA節點可以只持有部分樣本點數據。理想的情況是每個NUMA節點需要的樣本點數據都預先放在該節點本地內存中,從而每個NUMA節點都能快速且高效地訪問自己所需的樣本點數據[27]。此外,每個NUMA節點都需要持有所有中心點數據的一個副本。
根據對NUMA親和性優化的分析,可以通過以下三個步驟實現具有NUMA親和性的陰陽K-means算法:
1)線程綁定:該步驟指定線程如何綁定到處理器核心。OpenMP[28]4.0或更高版本提供了兩個環境變量——OMP_PLACES和OMP_PROC_BIND,這兩個環境變量通常結合使用。OMP_PLACES用于指定線程綁定到的機器上的位置;OMP_PROC_BIND指定綁定策略(線程關聯策略),它規定了線程如何分配到位置。在算法2的并行加速陰陽K-means算法實現中,指定線程綁定到處理器核心上,并且將每個線程分配到靠近父線程所綁定處理器核心位置的處理器核心。為了充分利用處理器的并行能力,使用與處理器核心(core)數相同的OpenMP線程數。基于以上策略將每個OpenMP線程逐一綁定到不同的處理器核心(第1行)。
2)數據分配:線程綁定完成后,每個NUMA節點上運行的線程已經確定。因此,每個節點所需的樣本點數據同樣是明確的。OpenMP在NUMA節點之間提供了一種名為First Touch Policy的數據分配方法,當線程第一次訪問某個樣本點數據時,會將該樣本點數據所在的頁分配到距離該線程最近的內存中。所以在陰陽K-means算法開始之前,可以并行地對所有的樣本點數據進行一次初始化(第2行),從而每個線程處理的樣本點數據都存放在該線程所屬節點的本地內存中。
3)工作負載分配:將所有的樣本點處理任務均勻地分成m個部分,然后靜態地分配給m個OpenMP線程(第8行),其中m大小等于眾核處理器中計算核心的數量。
完成以上三個步驟后,每個線程都綁定了一個處理器核心,且每個線程需要的樣本點數據都放在該處理器核心所屬NUMA節點的本地內存中。因此消除了各個NUMA節點對樣本點數據的遠程內存訪問,有效地降低了內存訪問開銷。
假設有n個樣本點,每個樣本點有d個維度,這些樣本點在內存中最常見的數據布局方式如圖3所示,樣本點與樣本點以及單個樣本點的每個維度之間都是順序存儲的。但是,當使用式(1)計算多個相鄰樣本點到一個中心點的歐幾里得距離時,由于基于樣本點數據進行向量化時需要訪問多個樣本點的相同維度來構造一個向量,經典內存數據布局將帶來一個明顯的缺點:相鄰樣本點的相同維度數據不連續,向量化存儲訪問開銷較大。

圖3 經典內存數據布局Fig.3 Classic memory data layout
如圖4所示,為了使兩個樣本點的相同維度的內存訪問連續,可以直觀地將所有樣本點的相同維度聚合在一起。但是,當訪問樣本點的下一個維度時,聚合內存數據布局將面臨兩次內存訪問之間的跨度太大的問題,數據的空間局部性較差。并且,在進行NUMA親和性優化時該內存數據布局不利于將樣本點數據分配到各個NUMA節點。

圖4 聚合內存數據布局Fig.4 Aggregated memory data layout
為了提高數據訪問局部性以及便于在各個NUMA節點間分配樣本點數據,提出了圖5所示的分組聚合內存數據布局。以ARMv8架構處理器的128位向量處理單元為例,將兩個樣本點的數據劃分為一組,并在組內使用聚合內存布局。至此,無論是訪問兩個樣本點的相鄰維度還是相鄰的一組樣本點,內存訪問都是完全連續的。

圖5 分組聚合內存數據布局Fig.5 Group aggregated memory data layout
為了證明向量化和NUMA親和性內存訪問優化對陰陽K-means算法的有效性,實驗使用四個大型的真實世界數據集,其中兩個來自UCI機器學習數據庫[29],另外兩個來自Kaggle[30]。表1列出了四個數據集的參數。

表1 真實數據集參數
實驗使用以下兩種陰陽K-means算法的實現:
1)YYG: 只使用全局過濾器的陰陽K-means算法。
2)YYGG: 使用全局過濾器和組過濾器的陰陽K-means算法。
此外,還使用一些縮寫來表示應用不同優化方法的YYG(G)實現:
1)YYG(G)-SV: 標量陰陽K-means算法在CPU上的多線程開源實現[25]。
2)YYG(G)-VV: 僅采用向量化優化的多線程陰陽K-means算法。
3)YYG(G)-NUMA: 同時采用了向量化優化和NUMA親和性優化的多線程陰陽K-means算法。
所有代碼都是用C/C++編寫,使用GCC編譯器編譯,并使用-O3級別的編譯優化以及Intrinsic指令來實現向量化。算法的收斂條件是每個簇的中心點位置不再變化或迭代次數達到1 000。在中心點數量k∈{64,128,256,512}的條件下使用四個數據集進行了實驗,每種情況下算法的性能測試重復三次,結果取平均值。所有數據都存儲為雙精度浮點值。實驗使用以下三個平臺:
1)Phytium FT-2000+平臺: 由一塊Phytium FT-2000+ARMv8 CPU構成,有8個NUMA節點,64個核心,核心的主頻為2.3 GHz,每個核心配備一個硬件線程、32 KB私有L1指令緩存以及32 KB私有L1數據緩存。每4個核心共享2 MB的L2緩存[22]。
2)Marvell ThunderX平臺: 兩塊Marvell ThunderX ARMv8 CPU構成一個NUMA系統,每塊CPU配備48個2.0 GHz的核心,每個核心配備一個硬件線程、78 KB的L1指令緩存和32 KB的L1數據緩存。48個核心共享16 MB的L2緩存[31]。
3)Intel Xeon Gold 6240平臺: 兩塊CPU構成一個NUMA系統,每塊CPU配備18個2.6 GHz的核心,每個核心配備兩個硬件線程,并分別配備576 KB的L1指令緩存和數據緩存。18個核心共享18 MB的L2緩存、24.75 MB的L3緩存。支持AVX-512向量指令集[32]。
為了驗證向量化的效果,對YYG(G)-SV和YYG(G)-VV這兩個算法進行了多次比較實驗。同時為了在NUMA節點上獲得穩定的性能,執行YYG(G)-SV和YYG(G)-VV算法時使用了Linux命令numactl--interleave=all來使所有數據均勻分布到各個節點。Phytium、Marvell和Intel三個平臺的向量化優化實驗結果所表現出的特征類似,分別最高獲得了約4.6、1.9以及8.5的性能加速比。考慮到篇幅原因,本小節僅以Phytium FT-2000+為例進行詳細分析。
表2展示了FT-2000+上向量化陰陽K-means算法相對于標量陰陽K-means算法的加速比。

表2 Phytium FT-2000+平臺上向量化YYG(G)算法 相對于標量YYG(G)算法的加速比Tab.2 Speedup of vectorized version of YYG(G) over scalar version on Phytium FT-2000+
由表可概括出如下幾個關鍵特征:
1)向量化優化效果明顯:在FT-2000+上,大多數情況下向量化算法的性能都優于標量算法,向量化算法的加速比最高達到約4.6。
2)向量化優化效果YYG算法明顯強于YYGG算法:可以從兩個方面來解釋,一方面,前者只有全局過濾器,因此算法中的邏輯判斷語句較少; 另一方面,后者篩選了部分中心點,進一步減少了樣本點與中心點的距離計算次數,導致計算量進一步降低,所以向量化之后的性能提升相對較小。
3)超線性加速比:當不存在NUMA架構的影響時,FT-2000+上ARMv8的128位向量單元對應的向量化理論加速比是2。實驗結果顯示某些情況下向量化加速比遠超過了2,主要是原始標量算法與優化向量算法受到NUMA架構的影響不一樣導致的。為了屏蔽NUMA架構對向量化優化效果的測量,本文又在FT-2000+的單個NUMA節點上進行了實驗,所有情況下向量化加速比都恢復到了小于2的正常范圍。因此,超線性加速比表明向量化算法比標量算法受NUMA架構的影響更小。
4)低維度數據集實驗效果不佳:在實驗中觀察到向量化算法在數據集YP上的性能比標量算法差。向量化在這些情況下不起作用的主要原因是該數據集的維度相對較低,從而計算量較小,指令流水線加載和排空的時間占總時間的比例較大,導致向量化性能不佳。因此,向量化更適合高維度的數據集。根據本文的測試結果顯示,使用維度超過100的數據集的陰陽K-means算法在向量化后可以獲得良好的性能。
5)加速比隨k值的增大而增大:實驗結果表明,隨著k值的增大,加速比整體呈上升趨勢。對于向量化算法來說,k值增大意味著中心點數量增多,從而樣本點與中心點之間的距離計算顯著增多,程序并行部分運行時間占總時間的比例增大,處理器數量不變,根據Amdahl定律可知加速比增大。
為了評估NUMA親和性內存訪問優化帶來的性能提升,將YYG(G)-NUMA算法與YYG(G)-VV進行性能比較。實驗結果顯示,當NUMA節點數較少時,如Intel平臺、Marvell平臺,以及僅使用兩個NUMA節點的Phytium平臺,NUMA親和性優化效果甚微。但是,隨著NUMA節點數的增加,NUMA親和性優化性能顯著增加。
表3和表4分別展示了同時使用向量化和NUMA親和性優化的YYG(G)算法相對于僅使用向量化優化的YYG(G)算法在8節點和4節點的Phytium平臺上的加速比。使用8節點的FT-2000+處理器時,YYG(G)-NUMA在四個數據集上的性能始終優于YYG(G)-VV,加速比為1.11~1.51。而對于使用4節點的FT-2000+而言,NUMA優化的加速比略低于8節點,最高只達到了1.34。因此,在NUMA節點數較多時,通過NUMA親和性內存訪問優化策略,可以有效地避免大量遠程內存訪問,實現性能提升。此外,隨著k值的增大,加速比整體呈下降趨勢,k值的增大意味著中心點數量增加,由于中心點只保存在主節點內存中,且中心點被所有NUMA節點共享,所以中心點數據增多導致其他NUMA節點訪問主節點的內存訪問開銷占比增大,從而加速比下降。

表3 Phytium FT-2000+平臺上具有NUMA親和性的向量化 YYG(G)相對于向量化YYG(G)的加速比(8節點)Tab.3 Speedup of NUMA-aware vectorized version of YYG(G) over vectorized version on Phytium FT-2000+(8 nodes)
結合向量化優化和NUMA親和性優化,YYG(G)-NUMA算法在Phytium平臺、Marvell平臺以及Intel平臺上相對于標量算法的整體加速比分別如表5、表6和表7中所示。當數據集的維度比較高時,并行加速陰陽K-means算法相比傳統陰陽K-means算法能實現較大的性能提升,這表明本文的優化方法可以有效地提高陰陽K-means算法在眾核處理器上的性能,在Phytium平臺、Marvell平臺以及Intel平臺上最高可分別獲得約5.6、2.0和8.7的加速比。

表5 Phytium FT-2000+平臺上具有NUMA親和性的 向量化YYG(G)相對于標量YYG(G)的加速比Tab.5 Speedup of NUMA-aware vectorized version of YYG(G) over scalar version on Phytium FT-2000+

表6 Marvell ThunderX平臺上具有NUMA親和性的 向量化YYG(G)相對于標量YYG(G)的加速比Tab.6 Speedup of NUMA-aware vectorized version of YYG(G) over scalar version on Marvell ThunderX
對比測試結果可知,本文算法在Intel x86平臺上的收益顯著高于Phytium和Marvell兩個ARMv8平臺,主要原因有兩點:一是Intel Xeon處理器配備了512位向量處理單元,相對于ARMv8架構128位向量處理單元,其向量寬度更寬,因此向量化優化收益更大;二是Intel平臺每個核心配備了更大的L1緩存、L2緩存以及獨有的L3緩存,使得本文數據布局優化的效果更加明顯,從而提升了整體優化方案的效果。
除上述實驗外,在Phytium平臺上以數據集WWP為例設計了一組實驗來研究陰陽K-means算法并行加速實現的強擴展性。分別測試了并行加速陰陽K-means算法和傳統陰陽K-means算法從單個NUMA節點擴展到8個NUMA節點時的強擴展效率,實驗結果如圖6和圖7所示。圖6表明向量化和NUMA親和性優化都提高了YYG的擴展效率,YYG-SV的擴展效率經過向量化后從30%提升到60%以上,經過NUMA親和性優化后進一步提升到85%以上。圖7中傳統的YYGG-SV算法的擴展效率已經高達90%以上,向量化大幅度提升了算法在單NUMA節點上的性能,從而計算出的擴展效率下降,但經過NUMA親和性優化后,擴展效率恢復到僅使用多線程技術實現的水平。因此,并行加速陰陽K-means算法依然具有較好的擴展性。

圖6 YYG算法的擴展性Fig.6 Strong scaling of YYG algorithm

圖7 YYGG算法的擴展性Fig.7 Strong scaling of YYGG algorithm
本文使用向量化和NUMA親和性內存訪問優化對眾核處理器上的陰陽K-means算法進行了優化,同時為了獲得更好的數據局部性,還提出了一種新的內存數據布局來支持上述優化。實驗結果表明,使用高維數據集時,向量化可以有效地提高陰陽K-means算法在眾核處理器上的性能;NUMA親和性內存訪問優化在具有較多NUMA節點的NUMA系統上可以進一步提升陰陽K-means算法的性能。在ARMv8和x86眾核平臺上,本文提出的陰陽K-means算法的并行加速實現最終分別獲得了最高約5.6倍與8.7倍的加速。
未來的研究擬將本文的陰陽K-means并行加速實現擴展到分布式系統中,以實現更大規模的并行。