集成眾核上快速獨(dú)立成分分析降維并行算法

方民權(quán)　張衛(wèi)民　周海芳

(國(guó)防科學(xué)技術(shù)大學(xué)計(jì)算機(jī)學(xué)院　長(zhǎng)沙　410073)(fmq@hpc6.com)

集成眾核上快速獨(dú)立成分分析降維并行算法

方民權(quán)張衛(wèi)民周海芳

(國(guó)防科學(xué)技術(shù)大學(xué)計(jì)算機(jī)學(xué)院長(zhǎng)沙410073)(fmq@hpc6.com)

摘要高光譜遙感影像快速獨(dú)立成分分析(fast independent component analysis, FastICA)降維過(guò)程包含大規(guī)模矩陣計(jì)算及大量迭代計(jì)算.通過(guò)熱點(diǎn)分析，面向集成眾核(many integrated core, MIC)架構(gòu)設(shè)計(jì)了協(xié)方差矩陣計(jì)算、白化處理和ICA迭代等熱點(diǎn)并行方案，提出和實(shí)現(xiàn)一種M-FastICA并行降維算法，并構(gòu)建算法性能模型；基于集成眾核研究并行程序優(yōu)化策略，針對(duì)各熱點(diǎn)并行方案提出一系列優(yōu)化策略，特別是創(chuàng)新性地提出一種下三角陣負(fù)載均衡方法，并量化測(cè)試其優(yōu)化效果.實(shí)驗(yàn)結(jié)果顯示M-FastICA算法最高可加速42倍，比24核CPU多線(xiàn)程并行快2.2倍；探討了波段數(shù)與并行程序性能的關(guān)系；實(shí)驗(yàn)數(shù)據(jù)驗(yàn)證了算法性能模型的準(zhǔn)確性.

關(guān)鍵詞集成眾核；獨(dú)立成分分析；高光譜影像降維；性能模型；下三角陣負(fù)載均衡

高光譜遙感技術(shù)廣泛應(yīng)用于軍事、農(nóng)業(yè)、環(huán)境科學(xué)、地質(zhì)、海洋學(xué)等領(lǐng)域，這些領(lǐng)域大都要求及時(shí)處理[1].但高光譜影像數(shù)據(jù)有波段多、數(shù)據(jù)量大、相關(guān)性強(qiáng)、冗余多等特點(diǎn)，直接處理將導(dǎo)致樣本類(lèi)別訓(xùn)練困難、維數(shù)災(zāi)難、空空間現(xiàn)象等嚴(yán)重的計(jì)算問(wèn)題[2-4].因此，國(guó)內(nèi)外專(zhuān)家學(xué)者處理高光譜影像數(shù)據(jù)過(guò)程常常包括降維步驟.降維是通過(guò)特定的映射，將高維影像數(shù)據(jù)轉(zhuǎn)換成低維影像數(shù)據(jù)，并保持信息量基本不變.

高光譜影像降維方法可分為線(xiàn)性和非線(xiàn)性2類(lèi)，線(xiàn)性降維有主成分分析(principal component analysis, PCA)[5]、獨(dú)立成分分析(independent component analysis, ICA)[6]等方法;非線(xiàn)性降維包括基于核的方法和基于流行學(xué)習(xí)的方法，如等距映射法(isometric feature mapping, Isomap)[7]、局部線(xiàn)性嵌入(locally linear embedding, LLE)[8]等.由于高光譜影像波段多、數(shù)據(jù)量大等特征，降維處理過(guò)程復(fù)雜且耗時(shí).提高降維速度是重要研究熱點(diǎn)，王和勇等人[9]研究基于聚類(lèi)和改進(jìn)距離的LLE方法以提高5.5倍降維速度；而并行降維是加快降維過(guò)程的另一個(gè)重要研究方向.

自Intel公司2012年推出集成眾核(many integ-rated core, MIC)產(chǎn)品Phi[10]，集成眾核已成為目前高性能計(jì)算領(lǐng)域的研究熱點(diǎn)，搭載MIC的天河2號(hào)超級(jí)計(jì)算機(jī)已連續(xù)6屆登臨TOP500榜首[11].

高光譜影像并行處理在傳統(tǒng)并行系統(tǒng)已有成熟的應(yīng)用：Valencia等人[12]基于異構(gòu)MPI在網(wǎng)絡(luò)機(jī)群系統(tǒng)研究了高光譜影像處理的技術(shù)；Plaza等人[13]提出基于神經(jīng)網(wǎng)絡(luò)的高光譜影像并行分類(lèi)算法.基于CPUGPU異構(gòu)系統(tǒng)，高光譜影像并行處理也有部分研究：Sánchez等人[14]對(duì)高光譜解混進(jìn)行了GPU移植；Plato?等人[15]和Ramalho等人[16]分別用GPU加速了非負(fù)矩陣分解和獨(dú)立成分分析等降維方法.但基于CPU+MIC這一新型高性能體系結(jié)構(gòu)還沒(méi)有相關(guān)研究.本文將面向集成眾核架構(gòu)探討加速高光譜影像降維的可行性.

本文以基于負(fù)熵最大的快速獨(dú)立成分分析(fast independent component analysis, FastICA)算法[17-18]為對(duì)象，面向集成眾核架構(gòu)展開(kāi)并行算法研究，本文貢獻(xiàn)包括：1)分析算法加速熱點(diǎn)并設(shè)計(jì)相應(yīng)的并行方案；2)針對(duì)各熱點(diǎn)并行方案提出并驗(yàn)證了多種建設(shè)性的優(yōu)化策略，特別是創(chuàng)新性地提出一種下三角陣負(fù)載均衡策略；3)提出并實(shí)現(xiàn)了一種并行降維算法M-FastICA；4)設(shè)計(jì)了該并行算法的性能預(yù)測(cè)模型，并驗(yàn)證了其有效性.

1集成眾核體系結(jié)構(gòu)與優(yōu)化策略

1.1集成眾核體系結(jié)構(gòu)

英特爾集成眾核MIC架構(gòu)集成了超過(guò)50個(gè)x86核心，由片上高速環(huán)形互聯(lián)總線(xiàn)聯(lián)接如圖1所示.MIC共有8個(gè)基于GDDR5的存儲(chǔ)控制器、16個(gè)通道,峰值帶寬可達(dá)350 GBps[10].

Fig. 1　Many integrated core architecture.圖1　MIC架構(gòu)

MIC的每個(gè)計(jì)算核心擁有1個(gè)512 b的向量處理單元(vector processing unit, VPU)、1個(gè)x86架構(gòu)標(biāo)量處理單元、32 KB的L1 cache(指令和數(shù)據(jù)各1個(gè))和512 KB的L2 cache.MIC支持硬件多線(xiàn)程技術(shù)，每個(gè)核心有4個(gè)硬件線(xiàn)程，但發(fā)射1條指令需要2個(gè)時(shí)鐘[10].

MIC搭載的嵌入式Linux uOS使其具有5種不同的工作模式：CPU原生模式、CPU為主MIC為輔模式、對(duì)等模式、MIC為主CPU為輔模式和MIC原生模式[10].本文使用CPU為主MIC為輔模式.

1.2集成眾核優(yōu)化策略

MIC是眾核協(xié)處理器，要獲得最優(yōu)性能，可以采取的一般優(yōu)化策略包括并行度、內(nèi)存空間、數(shù)據(jù)通信、Cache訪(fǎng)問(wèn)、向量化、負(fù)載均衡、線(xiàn)程擴(kuò)展性等方面優(yōu)化.對(duì)于CPUMIC異構(gòu)系統(tǒng)，如何合理分配CPU和MIC的計(jì)算任務(wù)、如何設(shè)計(jì)host端和device端的通信這些都是異構(gòu)系統(tǒng)優(yōu)化過(guò)程中所必須考慮的重要因素.

2獨(dú)立成分分析算法與加速熱點(diǎn)

2.1基于負(fù)熵最大的快速獨(dú)立成分分析

FastICA算法有基于負(fù)熵最大、基于似然最大、基于峭度3種形式，本文研究基于負(fù)熵最大的FastICA算法.該算法以負(fù)熵最大為搜尋方向，實(shí)現(xiàn)獨(dú)立源的順序提取，充分體現(xiàn)傳統(tǒng)線(xiàn)性變換中的投影追蹤思想，該算法還采用定點(diǎn)迭代的優(yōu)化算法，使收斂更加快速、穩(wěn)健[17-18].

對(duì)于高光譜影像數(shù)據(jù)X(W×H×B，寬W、高H、波段(band,B))，X可看成是B行S列(S=W×H)的2維矩陣，矩陣的1行表示1個(gè)波段的高光譜影像數(shù)據(jù).通過(guò)獨(dú)立成分分析方法降維，可以獲得m個(gè)IC圖像，其中m

算法1. FastICA高光譜影像降維方法[17-18].

輸入：高光譜影像數(shù)據(jù)X;

輸出：獨(dú)立成分Y.

Step1. 計(jì)算X的協(xié)方差矩陣；

Step2. 計(jì)算協(xié)方差矩陣的特征值及特征向量，并根據(jù)閾值T選取最終IC圖像數(shù)量m；

Step3. 計(jì)算白化矩陣M=DVT，其中D為X的特征值矩陣，V為對(duì)應(yīng)的特征向量矩陣；

Step4. 高光譜影像數(shù)據(jù)白化處理Z=MX，這里的X為減去均值之后的X(即零均值處理)；

Step5. 迭代開(kāi)始：

Step5.1. 初始化Wi；

Step5.4. 歸一化Wi=Wi；

Step5.5. 判斷是否收斂，不收斂則轉(zhuǎn)Step5.2，若收斂，如果所有的Wi已經(jīng)計(jì)算結(jié)束(即i=m)，則退出，否則令i=i+1，轉(zhuǎn)到Step5.1計(jì)算新的Wi；

Step6. 計(jì)算獨(dú)立成分Y=WZ.

2.2加速熱點(diǎn)分析

實(shí)現(xiàn)串行高光譜影像FastICA降維算法，對(duì)W=614,H=1087,B=224的高光譜影像降維，測(cè)試并統(tǒng)計(jì)計(jì)算部分(不含IO)各步驟占總計(jì)算時(shí)間的比率如表1所示.表1數(shù)據(jù)顯示，Step1,Step3,Step4,Step5比重較大，共占總計(jì)算時(shí)間的99%，是本文并行研究的關(guān)鍵熱點(diǎn).分析算法過(guò)程，計(jì)算量主要集中在Step1中協(xié)方差矩陣計(jì)算、Step4中白化處理時(shí)的矩陣運(yùn)算和Step5.2中Wi計(jì)算中的大量向量運(yùn)算.此外Step2中特征值計(jì)算取決于B，變化不大可不考慮并行；Step3白化矩陣很小，計(jì)算量有限，不必做并行考慮；而Step6中的IC變換與W,H,m正相關(guān)，作為次要熱點(diǎn)考慮其并行化.

Table 1　Time Distribution of FastICA

3FastICA集成眾核并行算法

3.1協(xié)方差矩陣任務(wù)劃分

協(xié)方差矩陣是對(duì)稱(chēng)矩陣，僅需計(jì)算下三角陣的協(xié)方差，然后在對(duì)稱(chēng)的上三角陣填入對(duì)應(yīng)的值.傳統(tǒng)的并行方案是對(duì)最外層循環(huán)進(jìn)行任務(wù)劃分，由于是三角矩陣，將任務(wù)調(diào)度過(guò)程設(shè)定為動(dòng)態(tài)調(diào)度，當(dāng)線(xiàn)程數(shù)遠(yuǎn)小于B時(shí)，能獲得較好的并行度和負(fù)載均衡；當(dāng)線(xiàn)程數(shù)與B相當(dāng)或比B大時(shí)，這種方案就會(huì)導(dǎo)致負(fù)載不均衡問(wèn)題.本文采用的高光譜影像數(shù)據(jù)的B=224，與集成眾核的硬件并發(fā)線(xiàn)程總數(shù)228相當(dāng)，即采用傳統(tǒng)方案并行計(jì)算協(xié)方差矩陣可能存在負(fù)載不均衡問(wèn)題.為解決這一問(wèn)題，本文提出一種下三角陣負(fù)載均衡方法，在本文第5節(jié)性能優(yōu)化部分詳細(xì)描述.

3.2白化處理與IC變換任務(wù)分配

白化處理包括白化矩陣計(jì)算和白化處理過(guò)程，其中白化矩陣維度較小、計(jì)算量也較小、無(wú)需并行，主要考慮白化處理Z=MX.針對(duì)白化處理過(guò)程，首先抽象并構(gòu)建白化處理宏觀(guān)模型(對(duì)算法核心計(jì)算的宏觀(guān)描述),如圖2所示，其中X需要進(jìn)行零均值處理，構(gòu)建宏觀(guān)模型時(shí)暫不考慮.

Fig. 2　Macroscopic model for whitening processing.圖2　白化處理宏觀(guān)模型

上述宏觀(guān)模型可視為矩陣乘法，傳統(tǒng)的并行矩陣乘法策略是在最外層循環(huán)做并行.但針對(duì)高光譜影像降維過(guò)程，降維后取得的波段數(shù)目m較小，而使用了硬件多線(xiàn)程并發(fā)技術(shù)的集成眾核并行度遠(yuǎn)大于m，將導(dǎo)致嚴(yán)重的負(fù)載不均衡.因此，本文對(duì)像元素S進(jìn)行任務(wù)劃分，由于S=W×H足夠大，可滿(mǎn)足大量線(xiàn)程的并行均衡負(fù)載要求；也可采用循環(huán)合并策略進(jìn)一步擴(kuò)大并行度，具體優(yōu)化策略參見(jiàn)本文第5節(jié)性能優(yōu)化部分.

IC變換(Y=WZ)的模型與白化處理(Z=MX)模型基本一致，因此采用與白化處理相同的并行任務(wù)劃分方案和優(yōu)化策略.

3.3集成眾核上的ICA迭代

ICA迭代的5個(gè)步驟中，計(jì)算量主要集中在Step5.2中.該運(yùn)算比較復(fù)雜，首先對(duì)其抽象，得到如圖3所示的宏觀(guān)模型，然后針對(duì)該模型進(jìn)行并行任務(wù)分割研究.

Fig. 3　Model of key step for ICA iteration.圖3　ICA迭代關(guān)鍵步驟模型

Fig. 4　Parallel model of ICA iteration.圖4　ICA迭代并行模型

3.4M-FastICA算法

基于上述熱點(diǎn)并行任務(wù)劃分策略，結(jié)合串行算法，對(duì)白化處理和ICA迭代并行宏觀(guān)模型進(jìn)行演繹，獲得熱點(diǎn)詳細(xì)并行方案.通過(guò)第5節(jié)對(duì)比各熱點(diǎn)并行方法在CPU與MIC上執(zhí)行性能和CPU+MIC協(xié)同工作策略，提出M-FastICA并行降維算法，如圖5所示，其中由于Wi正交化、歸一化及判斷等運(yùn)算量較小，ICA迭代過(guò)程可完全Offload到MIC.算法將協(xié)方差矩陣計(jì)算、ICA迭代Offload到MIC上進(jìn)行并行計(jì)算；白化處理、IC變換、轉(zhuǎn)置等CPU計(jì)算占優(yōu)勢(shì)的步驟仍在CPU上并行處理；另外充分利用了MIC與CPU協(xié)同計(jì)算、傳輸計(jì)算重疊等關(guān)鍵異構(gòu)優(yōu)化策略，圖5中縱向可視為時(shí)間軸，橫向重疊部分可同時(shí)執(zhí)行.其原理將在第5節(jié)詳細(xì)討論.

Fig. 5　Flow chart of M-FastICA.圖5　M-FastICA流程圖

4算法性能模型

Fig. 6　Performance model of M-FastICA.圖6　M-FastICA算法性能模型

本節(jié)旨在構(gòu)建M-FastICA算法性能模型，通過(guò)輸入高光譜影像數(shù)據(jù)參數(shù)(W,H,B,m等)，經(jīng)性能模型分析，預(yù)測(cè)出M-FastICA并行算法的降維時(shí)間，如圖6所示：

研究M-FastICA降維算法，建立性能模型：

T=TCOV+Teigen+Twhiteprocess+TICAiter+TICtrans；

(1)

TCOV=TCOVc+TCOVt=(BWH+

B(B+1)WH)vCOVc+(BWH+4BB)vPCI-e；

(2)

Twhiteprocess=(2mBWH)vwhiteprocess；

(3)

TICAiter=TZT×2+TeveryICA×niter≈

2(4mWH)vPCI-e+(4mWH)vICAiter×niter；

(4)

TICtrans=(2mmWH)vICtrans.

(5)

TCOV: 協(xié)方差矩陣計(jì)算時(shí)間；

TCOVc: MIC并行計(jì)算協(xié)方差矩陣時(shí)間；

TCOVt: 高光譜影像數(shù)據(jù)和協(xié)方差矩陣通信時(shí)間；

vCOVc: MIC計(jì)算協(xié)方差矩陣速率(圖7)；

vPCI-e: PCI-e總線(xiàn)傳輸速率，v1版本取3 GBps；

Teigen: 特征值和波段選擇時(shí)間，為0.28±0.05 s；

Twhiteprocess: 白化處理時(shí)間；

vwhiteprocess: CPU白化處理速率，如圖7所示；

TICAiter: ICA迭代時(shí)間；

TZT: 白化結(jié)果矩陣Z及其轉(zhuǎn)置矩陣的傳輸時(shí)間；

TeveryICA: 每次ICA迭代時(shí)間；

niter: 迭代次數(shù)；

vICAiter: MIC上并行ICA迭代性能，如圖7所示；

TICtrans: IC變換時(shí)間；

vICtrans: CPU處理IC變換的速率，如圖7所示:

Fig. 7　Performance of different steps.圖7　各步驟計(jì)算性能

MIC是高度并行計(jì)算協(xié)處理器，對(duì)并行度非常敏感，且在運(yùn)算規(guī)模較大時(shí)能獲得較好的吞吐率；問(wèn)題規(guī)模影響并行度，對(duì)MIC性能發(fā)揮有著較大的影響.本文通過(guò)高光譜影像數(shù)據(jù)縮放，獲得成比例的數(shù)據(jù)規(guī)模，測(cè)試并統(tǒng)計(jì)計(jì)算單元(CPU,MIC)性能.在利用該模型預(yù)測(cè)性能時(shí)，根據(jù)高光譜影像數(shù)據(jù)規(guī)模，從統(tǒng)計(jì)曲線(xiàn)中獲得估計(jì)速率，再根據(jù)算法模型公式預(yù)測(cè)并行降維時(shí)間.

5面向集成眾核的性能優(yōu)化

本節(jié)首先對(duì)協(xié)方差矩陣計(jì)算、白化處理、ICA迭代等展開(kāi)研究集成眾核性能優(yōu)化，然后研究了并行算法在CPU+MIC異構(gòu)系統(tǒng)中的協(xié)同計(jì)算.

5.1協(xié)方差矩陣計(jì)算優(yōu)化

針對(duì)協(xié)方差矩陣計(jì)算過(guò)程，面向集成眾核架構(gòu)進(jìn)行優(yōu)化研究，下列優(yōu)化方法可加速協(xié)方差矩陣計(jì)算：

1) 數(shù)據(jù)類(lèi)型選擇. 高光譜影像屬于圖像數(shù)據(jù)，由0～255等像素信息組成，可以采用unsigned char類(lèi)型(替代原始的float類(lèi)型)數(shù)據(jù)存儲(chǔ)信息，由于單位數(shù)據(jù)的存儲(chǔ)減少，可降低協(xié)處理器訪(fǎng)存開(kāi)銷(xiāo)和主機(jī)端設(shè)備端間通信開(kāi)銷(xiāo).

2) 計(jì)算分解. 相對(duì)于CPU強(qiáng)大的事務(wù)處理能力，作為協(xié)處理器的MIC更擅長(zhǎng)于較簡(jiǎn)單的計(jì)算.將復(fù)雜的計(jì)算任務(wù)分解為數(shù)個(gè)簡(jiǎn)單的計(jì)算，使其適用于MIC輕量級(jí)協(xié)處理器.協(xié)方差矩陣中每個(gè)元素的計(jì)算：

(6)

通過(guò)式(6)的協(xié)方差計(jì)算公式變形，將協(xié)方差的計(jì)算變形為累加和點(diǎn)乘運(yùn)算；在此基礎(chǔ)上還能將所有協(xié)方差矩陣元素的累加提取出來(lái)單獨(dú)計(jì)算，即將累加和矩陣乘運(yùn)算分開(kāi)進(jìn)行.該計(jì)算分解過(guò)程還可減少計(jì)算量，使計(jì)算量從2B(B+1)WH降為BWH+B(B+1)WH，減少了近半.算法2和算法3分別描述了計(jì)算分解前后協(xié)方差矩陣計(jì)算的詳細(xì)過(guò)程.

算法2. 計(jì)算分解前的協(xié)方差矩陣計(jì)算.

輸入：高光譜影像數(shù)據(jù)X；

輸出：協(xié)方差矩陣cov.

parallel fori=0 to (B-1)*B為波段數(shù)*

forj=0 toi*下三角*

fork=0 toS*S=W×H像元數(shù)*

x=X[i][k]；

y=X[j][k]；

m+=x;n+=y；

r+=x×y；

end for

cov[i][j]=(r-m×nS)(S-1)；*計(jì)算協(xié)方差*

cov[j][i]=cov[i][j]；*上三角陣賦值*

end for

end parallel for

算法3. 計(jì)算分解后的協(xié)方差矩陣計(jì)算.

輸入：高光譜影像數(shù)據(jù)X；

輸出：協(xié)方差矩陣cov.

parallel fori=0 to (B-1)

fork=0 toS

m[i]+=X[i][k]；*計(jì)算累加*

end for

end parallel for

parallel fori=0 to (B-1)

forj=0 toi

fork=0 toS

r+=X[i][k]×X[j][k]；*計(jì)算矩陣乘*

end for

cov[i][j]=(r-m[i]×m[j]S)(S-1)；*計(jì)算協(xié)方差*

cov[j][i]=cov[i][j]；*上三角陣賦值*

end for

end parallel for

Fig. 8　Load balancing for lower triangular matrix.圖8　下三角陣負(fù)載均衡

3) 下三角陣負(fù)載均衡. 協(xié)方差矩陣是對(duì)稱(chēng)陣，僅需計(jì)算下三角陣，由于MIC核心較多，2.2節(jié)的并行方案無(wú)法合理地做到下三角陣的負(fù)載均衡，如線(xiàn)程0執(zhí)行1個(gè)協(xié)方差計(jì)算，而線(xiàn)程223需要執(zhí)行224個(gè)協(xié)方差計(jì)算.為實(shí)現(xiàn)下三角陣的負(fù)載均衡，本文設(shè)計(jì)了一種新的負(fù)載均衡方法，如圖8所示.其基本思想是將三角形狀的2維計(jì)算任務(wù)映射為1維計(jì)算任務(wù)，然后對(duì)得到的1維計(jì)算任務(wù)進(jìn)行并行劃分，最終實(shí)現(xiàn)均衡負(fù)載的目的.本文通過(guò)圖8中行③和行④的公式，將1維索引映射到1維數(shù)組中，把2層循環(huán)合并為1層循環(huán)，不僅可以實(shí)現(xiàn)負(fù)載均衡，還可增大并行度.算法4采用了下三角陣負(fù)載均衡策略.

算法4. 下三角負(fù)載均衡后的協(xié)方差矩陣計(jì)算.

輸入：高光譜影像數(shù)據(jù)X；

輸出：協(xié)方差矩陣cov.

parallel fori=0 to (B-1)

fork=0 toS

m[i]+=X[i][k]；*計(jì)算累加*

end for

end parallel for

parallel forii=0 to (n-1)

i=ceil((sqrt(9.0+8×ii)-3)2)；

j=ii-(1+i)×i2；

fork=0 toS

r+=X[i][k]×X[j][k]；*計(jì)算矩陣乘*

end for

cov[i][j]=(r-m[i]×m[j]S)(S-1)；*計(jì)算協(xié)方差*

cov[j][i]=cov[i][j]；*上三角陣賦值*

end parallel for

程序經(jīng)上述優(yōu)化方法優(yōu)化后進(jìn)行實(shí)驗(yàn)測(cè)試，其效果如圖9所示，圖9中0表示未采取優(yōu)化手段，1～3分別表示相繼采用優(yōu)化手段1～3后的時(shí)間消耗.

Fig. 9　Optimization effect of COV calculation.圖9　COV計(jì)算優(yōu)化效果

5.2白化處理集成眾核優(yōu)化

針對(duì)白化處理過(guò)程，本文采取的優(yōu)化方法包括：

1) 并行循環(huán)選擇.高光譜影像白化處理過(guò)程有其獨(dú)特的特點(diǎn)，最外層循環(huán)量是根據(jù)貢獻(xiàn)率所取得的波段數(shù)目m，該值往往為十幾到幾十不等，該層循環(huán)并行度較小，不適合在眾核處理器(MIC上有50多個(gè)核，每個(gè)核可并發(fā)4個(gè)線(xiàn)程)上并行.選擇適合眾核架構(gòu)核數(shù)的并行度的循環(huán)開(kāi)發(fā)并行至關(guān)重要.本文選擇第2層循環(huán)開(kāi)發(fā)并行，其并行度(S=W×H)巨大.

2) 數(shù)據(jù)類(lèi)型選擇.同協(xié)方差矩陣計(jì)算過(guò)程.

3) 循環(huán)交換.當(dāng)存儲(chǔ)訪(fǎng)問(wèn)不連續(xù)時(shí)，cache命中率低，將導(dǎo)致大量的內(nèi)存訪(fǎng)問(wèn)開(kāi)銷(xiāo)；最內(nèi)層循環(huán)無(wú)法向量化時(shí)，MIC的重要部件向量處理單元(vector process unit, VPU)得不到有效利用，將大大影響性能.若程序存在上述2方面情況，需要循環(huán)交換.白化處理由于存儲(chǔ)訪(fǎng)問(wèn)不連續(xù)而采取循環(huán)交換措施，通過(guò)交換內(nèi)2層循環(huán)使數(shù)據(jù)連續(xù)訪(fǎng)問(wèn).

4) 除最內(nèi)層的嵌套循環(huán)展開(kāi).白化處理共3層循環(huán)，最外層取值m(降維后取得的波段數(shù))較小而不適合MIC并行，通過(guò)將外2層嵌套循環(huán)展開(kāi)，增加最外層循環(huán)次數(shù)，使其適用于MIC等眾核架構(gòu)，還能減少并行區(qū)域fork-join開(kāi)銷(xiāo).在嵌套循環(huán)展開(kāi)的同時(shí)保持最內(nèi)層循環(huán)不變，以此來(lái)保證向量化和連續(xù)存儲(chǔ)訪(fǎng)問(wèn).

5) 高光譜影像矩陣轉(zhuǎn)置.白化處理對(duì)高光譜影像矩陣的訪(fǎng)問(wèn)是按列的，為了使訪(fǎng)問(wèn)存儲(chǔ)連續(xù)，可以采用矩陣轉(zhuǎn)置的策略.通過(guò)轉(zhuǎn)置，原始矩陣行列交換，實(shí)現(xiàn)存儲(chǔ)連續(xù)訪(fǎng)問(wèn).該方案將引入轉(zhuǎn)置開(kāi)銷(xiāo)，若轉(zhuǎn)置開(kāi)銷(xiāo)大于性能收益，則不可取(效果統(tǒng)計(jì)時(shí)已包含轉(zhuǎn)置開(kāi)銷(xiāo)).

優(yōu)化方法3和優(yōu)化方法5無(wú)法重復(fù)使用，因其目的是一致的，即連續(xù)訪(fǎng)存.圖10所示6種不同優(yōu)化方法組合效果，其中組0沒(méi)有優(yōu)化，組1用了方法1，組2用了方法1和方法2，組3采取方法1至方法3，方法1～4組成組4，組5用了方法1、方法2、方法5優(yōu)化方法，組6采用方法1、方法2、方法4、方法5.通過(guò)一系列的優(yōu)化，組6獲得最佳性能.

Fig. 10　Optimization effect of whitening process.圖10　白化處理優(yōu)化效果

5.3ICA迭代優(yōu)化

ICA迭代部分包含大量的小規(guī)模迭代計(jì)算，本文對(duì)單次迭代展開(kāi)優(yōu)化研究，主要采取3種優(yōu)化策略：

2) ICA迭代過(guò)程中，參與計(jì)算的Wi是W矩陣的1列，迭代過(guò)程中Wi參與大量計(jì)算，因此本文采用中間數(shù)組存儲(chǔ)Wi中元素參與計(jì)算.C語(yǔ)言中矩陣是按行存儲(chǔ)的，列訪(fǎng)問(wèn)導(dǎo)致Cache命中率下降.利用中間數(shù)組tmp[]存儲(chǔ)Wi的元素，那么運(yùn)算時(shí)訪(fǎng)存是順序的，Cache命中率高，從而提升程序性能.

3) 使用較底層的intrinsics函數(shù)，顯式地指定數(shù)據(jù)向量化處理過(guò)程，以充分發(fā)揮MIC中向量處理單元的性能.intrinsics函數(shù)是類(lèi)匯編函數(shù)，指定程序中的向量運(yùn)算，比如向量?jī)?nèi)積(c=A·B)運(yùn)算時(shí)，普通C程序代碼和intrinsics代碼如下所示：

C版向量?jī)?nèi)積:

for(i=0;i

c+=a[i]+b[i]；

endfor

intrinsics版向量?jī)?nèi)積:

_m512_a,_b,_c；

for(i=0;i

_a=_mm512_load_ps(&a[i])；

_b=_mm512_load_ps(&b[i])；

_c=_mm512_add_ps(_c,_mm512_mul_ps(_a,_b))；

endfor

c=_mm512_reduce_add_ps(_c)；

其中intrinsics代碼中，1條指令完成16個(gè)浮點(diǎn)數(shù)的運(yùn)算.

在程序中分別實(shí)現(xiàn)上述優(yōu)化方法，測(cè)試出各優(yōu)化方法使用前后ICA迭代步驟的時(shí)間消耗，獲得圖11所示的優(yōu)化效果，其中橫坐標(biāo)0表示未采取優(yōu)化措施，1～3分別表示相繼采用優(yōu)化方法1～3后的時(shí)間.

Fig. 11　Optimization effect of ICA iteration.圖11　ICA迭代優(yōu)化效果

1) 異構(gòu)協(xié)同計(jì)算優(yōu)化.CPU+MIC架構(gòu)是典型的異構(gòu)系統(tǒng)，CPU和MIC可同時(shí)進(jìn)行計(jì)算或通信，通過(guò)協(xié)同計(jì)算可掩蓋部分步驟處理時(shí)間.本文算法中，協(xié)方差矩陣計(jì)算(COV)和高光譜影像數(shù)據(jù)X轉(zhuǎn)置的過(guò)程可同時(shí)執(zhí)行，如圖12所示；白化結(jié)果矩陣Z的傳輸可與矩陣Z轉(zhuǎn)置同時(shí)執(zhí)行，如圖13所示.通過(guò)上述CPU與MIC協(xié)同工作，可以掩蓋X和Z的轉(zhuǎn)置開(kāi)銷(xiāo).

Fig. 12　Cooperative computing between COV and X transposition.圖12　COV與X轉(zhuǎn)置協(xié)同計(jì)算

Fig. 13　Transpose Z and transfer Z at the same time.圖13　Z轉(zhuǎn)置與通信重疊

Fig. 14　Performance ratios of MICCPU for different steps.圖14　各步驟MICCPU性能對(duì)比

2) CPU與MIC計(jì)算的權(quán)衡.CPU和MIC都是異構(gòu)系統(tǒng)中重要的計(jì)算資源，由于體系結(jié)構(gòu)的不同，其擅長(zhǎng)領(lǐng)域也不同.本文通過(guò)實(shí)驗(yàn)統(tǒng)計(jì)，獲得各步驟MIC與CPU上執(zhí)行時(shí)間比值(圖14中橫軸1～6分別表示COV、X轉(zhuǎn)置、白化處理、Z轉(zhuǎn)置、ICA迭代和IC變換).比值大于1時(shí)CPU計(jì)算有優(yōu)勢(shì)，小于1則Offload到MIC能獲得更好的性能.觀(guān)察圖14，當(dāng)數(shù)據(jù)較大時(shí)，COV過(guò)程O(píng)ffload到MIC合適，ICA迭代過(guò)程O(píng)ffload到MIC可獲得更好的性能，其他過(guò)程在CPU端處理更劃算.

6實(shí)驗(yàn)結(jié)果與分析

6.1實(shí)驗(yàn)準(zhǔn)備

本文測(cè)試平臺(tái)是天河2號(hào)超級(jí)計(jì)算機(jī).天河2號(hào)單個(gè)節(jié)點(diǎn)搭載2個(gè)12核Xeon E5-2692v2 CPU和57核Intel Phi 31S1P協(xié)處理器.根據(jù)熱點(diǎn)并行方案和優(yōu)化策略，采用OpenMP和LEO(language extensions for offload)的C語(yǔ)言擴(kuò)展來(lái)實(shí)現(xiàn)M-FastICA并行算法，采用Intel C Compiler 2013 sp1.1.106編譯器、O3編譯選項(xiàng)，最后在天河2號(hào)單節(jié)點(diǎn)上測(cè)試性能.

表2羅列了本文采用的AVIRIS高光譜影像數(shù)據(jù)的詳細(xì)信息，共4組.

Table 2　Information of Hyperspectral Image Data

由于高光譜影像數(shù)據(jù)和ICA迭代過(guò)程中取得的隨機(jī)數(shù)不同，導(dǎo)致ICA迭代次數(shù)不同，不同的迭代次數(shù)將導(dǎo)致消耗的時(shí)間也不同.為了獲得準(zhǔn)確的性能對(duì)比數(shù)據(jù)，本文約定：同一數(shù)據(jù)取固定的迭代次數(shù).通過(guò)多次實(shí)驗(yàn)求平均值的方法，確定各組數(shù)據(jù)的迭代次數(shù)，如表3所示：

Table 3　Iteration Number of Hyperspectral Image Data

6.2并行算法加速比

本實(shí)驗(yàn)分別在MIC和CPU上測(cè)試并行程序計(jì)算時(shí)間，與打開(kāi)最佳優(yōu)化開(kāi)關(guān)(O3)的串行程序?qū)Ρ龋@得加速比繪制成圖15.其中MIC啟動(dòng)時(shí)間和MIC上存儲(chǔ)分配時(shí)間是一次性開(kāi)銷(xiāo)，不計(jì)入本文時(shí)間統(tǒng)計(jì)范疇.圖15中處理高光譜影像數(shù)據(jù)1時(shí)，CPU上的加速比高于MIC，可知CPU比MIC更擅長(zhǎng)于處理小規(guī)模運(yùn)算；高光譜影像數(shù)據(jù)2～4中，MIC獲得了比CPU并行更高的性能，說(shuō)明MIC適合規(guī)模較大的計(jì)算；在處理高光譜影像數(shù)據(jù)4時(shí)，M-FastICA算法可獲得42倍加速比，比24核CPU并行快2.2倍.

Fig. 15　Speed-up of MIC and CPU.圖15　MIC與CPU加速比

6.3波段數(shù)與性能的探討

本文選用的AVIRIS高光譜影像數(shù)據(jù)中B=224，剛好與Intel Phi 31S1P協(xié)處理器的硬件線(xiàn)程數(shù)量228接近，因此會(huì)提出疑問(wèn)，如果B變化時(shí)是否會(huì)影響并行程序性能？

當(dāng)B發(fā)生變化時(shí)，各熱點(diǎn)的計(jì)算量也會(huì)產(chǎn)生變化，因此并行程序性能會(huì)有所變化，但總體穩(wěn)定.究其原因，本文對(duì)各熱點(diǎn)的并行程序進(jìn)行了大量深入優(yōu)化，比如下三角矩陣均衡負(fù)載、循環(huán)合并等，已經(jīng)將受B限制的并行度進(jìn)行了擴(kuò)展，并行粒度不會(huì)受B變化而發(fā)生較大變化.

Fig. 16　Speed-ups with different bands.圖16　波段數(shù)與加速比的關(guān)系

本文通過(guò)控制變量法驗(yàn)證B與性能的關(guān)系.選取輸入高光譜影像數(shù)據(jù)3，取其中前112,140,168,196,224波段的影像數(shù)據(jù)作為輸入數(shù)據(jù)，通過(guò)實(shí)驗(yàn)統(tǒng)計(jì)并行算法對(duì)最優(yōu)串行的加速比，如圖16所示.圖16中B=112和B=140時(shí)，其加速比較小，原因是這2組輸入降維時(shí)最終IC圖像數(shù)量為5,6，而其他幾組輸入降維時(shí)最終IC圖像數(shù)量為13,14，這會(huì)影響白化處理、ICA迭代和IC變化等過(guò)程的計(jì)算量及性能.而當(dāng)IC圖像數(shù)量不變時(shí)，程序性能基本穩(wěn)定，不會(huì)因B的變化而影像程序性能.

6.4算法性能模型驗(yàn)證

將高光譜影像數(shù)據(jù)的參數(shù)(B,W,H,m)輸入M-FastICA算法性能模型，通過(guò)式(1)～(5)計(jì)算可得到預(yù)測(cè)時(shí)間；執(zhí)行優(yōu)化后的并行程序，獲得實(shí)際運(yùn)行時(shí)間，如圖17所示.對(duì)比預(yù)測(cè)時(shí)間和執(zhí)行時(shí)間的誤差，圖18所示為誤差百分比.圖17顯示預(yù)測(cè)時(shí)間與真實(shí)運(yùn)行時(shí)間基本一致，圖18顯示預(yù)測(cè)性能和實(shí)測(cè)性能最大誤差小于5%，說(shuō)明本文提出的M-FastICA算法性能模型能較準(zhǔn)確地預(yù)測(cè)算法性能.

Fig. 17　Time of real and perdiction for M-FastICA model.圖17　M-FastICA模型預(yù)測(cè)與實(shí)測(cè)時(shí)間

Fig. 18　Bias ratio of the M-FastICA model.圖18　M-FastICA性能模型預(yù)測(cè)誤差

6.5并行算法瓶頸分析

Fig. 19　Time distribution of M-FastICA.圖19　M-FastICA時(shí)間分布

統(tǒng)計(jì)獨(dú)立成分分析降維算法各個(gè)步驟消耗時(shí)間占總時(shí)間的百分比，繪制成圖19(橫軸1～5分別表示獨(dú)立成分分析的5個(gè)關(guān)鍵步驟：協(xié)方差矩陣計(jì)算、特征值特征向量計(jì)算、白化處理、ICA迭代、IC變換).對(duì)比圖19和表1，我們發(fā)現(xiàn)并行后各步驟耗時(shí)比重有所改變，但協(xié)方差矩陣計(jì)算和ICA變換依舊是算法性能瓶頸.要解決這一瓶頸，可以通過(guò)擴(kuò)展計(jì)算資源來(lái)解決，比如使用MIC集群，這也是下一步工作的重點(diǎn).

7結(jié)論

本文基于集成眾核MIC研究高光譜影像FastICA降維方法，設(shè)計(jì)了協(xié)方差矩陣計(jì)算、白化處理、ICA迭代、IC變換等熱點(diǎn)并行方案和優(yōu)化策略；提出和實(shí)現(xiàn)了基于MIC的M-FastICA算法，給出算法性能模型并驗(yàn)證其準(zhǔn)確性.實(shí)驗(yàn)結(jié)果顯示M-FastICA算法獲得了良好的性能提升，對(duì)比最優(yōu)串行程序最高可加速42倍，比24核CPU多線(xiàn)程快2.2倍.

本文研究顯示，高光譜影像FastICA降維算法可在集成眾核架構(gòu)上獲得良好的性能，但仍存在性能瓶頸，可考慮進(jìn)一步將算法擴(kuò)展到MIC集群以突破性能瓶頸.

參考文獻(xiàn)

[1]Green R O, Eastwood M L, Sarture C M, et al. Imaging spectroscopy and the airborne visibleinfrared imaging spectrometer(AVIRIS)[J]. Remote Sensing of Environment, 1998, 65(3): 227-248

[2]Green A A, Berman M, Switzer P, et al. A transformation for ordering multispectral data in terms of image quality with implications for noise removal[J]. IEEE Trans on Geoscience and Remote Sensing, 2000, 26(1): 65-74

[3]Kaarna A, Zemcik P, Kalviainen H, et al. Compression of multispectral remote sensing images using clustering and spectral reduction[J]. IEEE Trans on Geoscience and Remote Sensing, 2000, 38(2): 1073-1082

[4]Scott D W, Thompson J R. Probability density estimation in higher dimensions[C]Proc of the 15th Symp on the Interface of Computer Science and Statistics. Amsterdam: North-Holland, 1983: 173-179

[5]Jolliffe I T. Principal Component Analysis[M]. Berlin: Springer, 2002

[6]Hyv?rinen A, Oja E. Independent component analysis: Algorithm and applications[J]. Neural Networks, 2000, 13(45): 411-430

[7]Tenenbaum J, De Silva V, Langford J C. A global geometric framework for nonlinear dimensionality reduction[J]. Science, 2000, 290(5500): 2319-2323

[8]Roweis S T, Saul L K. Nonlinear dimensionality reduction by locally linear embedding[J]. Science, 2000, 290(5500): 2323-2326

[9]Wang Heyong, Zheng Jie, Yao Zheng’an, et al. Application of dimension reduction on using improved LLE based on clustering[J]. Journal of Computer Research and Development, 2006, 43(8): 1485-1490 (in Chinese)(王和勇, 鄭杰, 姚正安, 等. 基于聚類(lèi)和改進(jìn)距離的LLE方法在數(shù)據(jù)降維中的應(yīng)用[J]. 計(jì)算機(jī)研究與發(fā)展, 2006, 43(8): 1485-1490)

[10]Jeffers J, Reinders J. Intel Xeon Phi Coprocessor High Performance Programming[M]. Amsterdam: Elsevier, 2013(in Chinese)(James J, James R. Intel Xeon Phi協(xié)處理器高性能編程指南[M]. 陳健, 等譯. 北京: 人民郵電出版社, 2013)

[11]Meuer H, Strohmaier E, Dongarra J, et al. TOP500 supercomputer sites[EBOL]. 2014 [2014-09-26]. http:www.top500.orglists201406

[12]Valencia D, Lastovetsky A, O’Flynn M, et al. Parallel processing of remotely sensed hyperspectral images on heterogeneous networks of workstations using HeteroMPI[J]. International Journal of High Performance Computing Applications, 2008, 22(4): 386-407

[13]Plaza J, Plaza A, Pérez R, et al. Parallel classification of hyperspectral images using neural networks[J]. Computational Intelligence for Remote Sensing Studies in Computational Intelligence, 2008, 133: 193-216

[14]Sánchez S, Ramalho R, Sousa L, et al. Real-time implementation of remotely sensed hyperspectral image unmixing on GPUs[JOL]. Journal of Real-Time Image Processing, 2012 [2014-09-26]. http:link.springer.comarticle10.1007s11554-012-0269-2

[15]Plato? J, Gajdo? P, Kr?mer P, et al. Non-negative matrix factorization on GPU[C]Proc of the 2nd Int Conf on

Networked Digital Technologies. Berlin: Springer, 2010: 21-30

[16]Ramalho R, Tomas P, Sousa L. Efficient independent component analysis on a GPU[C]Proc of the 10th IEEE Int Conf on Computer and Information Technology. Piscataway, NJ: IEEE, 2010: 1128-1133

[17]Hyvrinen A. Fast and robust fixed-point algorithms for independent component analysis[J]. IEEE Trans on Neural Networks, 1999, 10(3): 626-634

[18]Hyvrinen A. The fixed-point algorithm and maximum likelihood estimation for independent component analysis[J]. Neural Processing Letters, 1999, 10(1): 1-5

Fang Minquan, born in 1989. PhD candidate. Student member of China Computer Federation. His main research interests include parallel computing and data assimilation.

Zhang Weimin, born in 1966. Professor and PhD supervisor. His main research interests include parallel computing and data assimilation.

Zhou Haifang, born in 1975. Professor. Member of China Computer Federation. Her main research interests include parallel computing and image processing.

Parallel Algorithm of Fast Independent Component Analysis for Dimensionality Reduction on Many Integrated Core

Fang Minquan, Zhang Weimin, and Zhou Haifang

(CollegeofComputer,NationalUniversityofDefenseTechnology,Changsha410073)

AbstractThere are massive matrix and iterative calculations in fast independent component analysis (FastICA) for hyperspectral image dimensionality reduction. By analyzing hotspots of FastICA algorithm, we design the parallel schemes of covariance matrix calculating, whitening processing and ICA iteration on many integrated core (MIC), implement and validate an M-FastICA algorithm. Further, we present a performance model for M-FastICA. We present a series of optimization methods for the parallel schemes of different hotspots: reforming the arithmetic operations, interchanging and unrolling loops, transposing matrix, using intrinsics and so on. In particular, we propose a novel method to balance the loads when dealing with the lower triangular matrix. Then we measure the performance effects of such optimization methods. Our experiments show that the M-FastICA algorithm can reach a maximum speed-up of 42X times in our test, and it runs 2.2X times faster than the CPU parallel version on 24 cores. We also investigate how the speed-ups change with the bands. The experiment results validate our performance model with an acceptable accuracy and thus can provide a roofline for our optimization effort.

Key wordsmany integrated core (MIC); independent component analysis (ICA); dimensionality reduction of hyperspectral image; performance model; load balancing of lower triangular matrix

收稿日期：2014-09-29；修回日期：2015-04-16

基金項(xiàng)目：國(guó)家自然科學(xué)基金項(xiàng)目(61272146，41375113); 湖南省研究生創(chuàng)新資助項(xiàng)目(CX2015B030)

中圖法分類(lèi)號(hào)TP391.4

This work was supported by the National Natural Science Foundation of China (61272146，41375113) and the Graduate Innovative Foundation of Hunan Province of China (CX2015B030).