基于CUBLAS和CUDA的MNF并行算法設(shè)計(jì)與優(yōu)化

周海芳+高暢+方民權(quán)

摘 要：為實(shí)現(xiàn)高光譜影像數(shù)據(jù)快速降維，基于nVidia 的圖像處理單元（graphic processing unit， GPU）研究最大噪聲分?jǐn)?shù)變換（Maximum Noise Fraction Rotation，MNF Rotation）降維算法的并行設(shè)計(jì)與優(yōu)化，通過(guò)對(duì)加速熱點(diǎn)并行優(yōu)化，擇優(yōu)整合，設(shè)計(jì)并實(shí)現(xiàn)基于CUBLAS（CUDA Basic Linear Algebra Subprograms）庫(kù)的MNF-L（MNF-on-Library）算法和基于CPU/GPU異構(gòu)系統(tǒng)的MNF-C（MNF-on-CUDA）算法.實(shí)驗(yàn)結(jié)果顯示MNF-L算法加速11.5～60.6倍不等，MNF-C算法加速效果最好，加速46.5～92.9倍不等.研究結(jié)果表明了GPU在高光譜影像線性降維領(lǐng)域的巨大優(yōu)勢(shì).

關(guān)鍵詞：圖像處理單元；GPU性能優(yōu)化；高光譜影像降維；最大噪聲分?jǐn)?shù)變換；協(xié)方差矩陣計(jì)算

中圖分類號(hào)：TP391.4 文獻(xiàn)標(biāo)志碼：A

Parallel Algorithm Design and Performance Optimization of Maximum Noise Fraction Rotation Based on CUBLAS and CUDA

ZHOU Haifang ，GAO Chang， FANG Minquan

（School of Computer， National University of Defense Technology， Changsha 410073， China ）

Abstract：To rapidly reduce the huge dimensions of hyperspectral image， this paper investigated the design and optimization of the parallel Maximum Noise Fraction （MNF） Rotation algorithm based on nVidia graphic processing units （GPUs）. In particular， a MNF-L （MNF-on-Library） algorithm based on the CUBLAS library functions and a MNF-C（MNF-on-CUDA） algorithm on the CPU/GPU heterogeneous system was presented by designing mapping schemes and parallel optimizing strategies. Experiment result shows that the MNF-L algorithm can obtain the speedups between 11.5 and 60.6， and the MNF-C algorithm can get the speedups between 46.5 and 92.9. Therefore， GPUs have a great advantage in reducing dimensions of hyperspectral images.

Key words： graphic processing unit； performance optimization for GPU； dimensionality reduction of hyperspectral image； maximum noise fraction rotation； covariance matrix calculation

高光譜遙感影像具有波段連續(xù)、光譜分辨率高的特點(diǎn)，能從其光譜空間中獲取豐富的地物特征信息，因其“圖譜合一”的優(yōu)勢(shì)，廣泛應(yīng)用于農(nóng)業(yè)、林業(yè)、軍事、環(huán)境科學(xué)、地質(zhì)等各領(lǐng)域，具有很好的實(shí)用性和研究?jī)r(jià)值[1-2].在數(shù)據(jù)處理過(guò)程中，連續(xù)波段成像導(dǎo)致高光譜影像數(shù)據(jù)量龐大，且連續(xù)波段之間數(shù)據(jù)相關(guān)性強(qiáng)，信息冗余大，存儲(chǔ)處理困難，為應(yīng)用和分析帶來(lái)不便，如產(chǎn)生維數(shù)災(zāi)難、Hughes現(xiàn)象等[3]，因此，數(shù)據(jù)降維應(yīng)運(yùn)而生.

怎樣將高維空間數(shù)據(jù)映射到低維子空間，同時(shí)保持重要特征不被丟失是高光譜數(shù)據(jù)降維遵循的基本原則.高光譜影像降維主要涉及矩陣操作，如濾波、矩陣乘、協(xié)方差矩陣計(jì)算等，是典型的計(jì)算密集型和訪存密集型過(guò)程，傳統(tǒng)的降維過(guò)程一般采用串行方式進(jìn)行，計(jì)算復(fù)雜度高，耗時(shí)長(zhǎng)，無(wú)法滿足各領(lǐng)域?qū)Ω吖庾V數(shù)據(jù)及時(shí)處理的需求[4-5].

2007年，nVidia公司發(fā)布統(tǒng)一計(jì)算設(shè)備架構(gòu)（Compute Unified Device Architecture， CUDA），GPU并行計(jì)算開始在科學(xué)計(jì)算領(lǐng)域承擔(dān)重要角色[6].CPU+GPU異構(gòu)系統(tǒng)在高性能計(jì)算領(lǐng)域表現(xiàn)突出，天河1A、泰坦等超級(jí)計(jì)算機(jī)相繼成為TOP500榜首[7].

高光譜影像并行處理已有廣泛研究，方民權(quán)[8-9]等人分別基于GPU和MIC研究了高光譜影像主成分分析降維算法，在2個(gè)GPU上獲得了128倍加速比，在3個(gè)MIC上加速133倍；Sergio[10]等人基于GPU集群實(shí)現(xiàn)了高光譜影像實(shí)時(shí)解混；Elmaghrbay[11]等人提出了高光譜影像端元提取的快速GPU算法；Chang[12]等人實(shí)現(xiàn)了一種并行模擬退火方法，并成功應(yīng)用到高維遙感影像數(shù)據(jù)特征提??；Santos[13]等人基于GPU平臺(tái)實(shí)現(xiàn)了高光譜影像有損壓縮算法的并行加速，表明GPU在高效數(shù)據(jù)處理方面的巨大潛力.羅耀華[14]等人首次基于GPU實(shí)現(xiàn)了MNF并行方法研究，在數(shù)據(jù)規(guī)模為1 036 × 235 × 229時(shí)，最高取得了4倍的加速比，但該算法僅對(duì)加速熱點(diǎn)中的協(xié)方差矩陣運(yùn)算進(jìn)行了并行移植，且涉及的優(yōu)化分析較少，加速效果較差.MNF作為高光譜數(shù)據(jù)特征提取的主流算法，鮮有較全面的并行化研究，本文研究正以此為契機(jī)，對(duì)該算法的并行設(shè)計(jì)進(jìn)行深入分析.

最大噪聲分?jǐn)?shù)變換將原始數(shù)據(jù)中的噪聲分離，提取影像數(shù)據(jù)的主要特征，從而表征主要信息[14-15].該算法由兩層主成分變換構(gòu)成，在主成分分析基礎(chǔ)上考慮了噪聲對(duì)圖像的影響，具有更好的效果，是目前主流的線性降維算法，實(shí)現(xiàn)其并行化在高光譜降維處理領(lǐng)域具有重要意義.

CUBLAS庫(kù)函數(shù)[16]和CUDA均可實(shí)現(xiàn)算法在GPU上并行加速，因其突出的加速效果而被廣泛使用.采用CUDA程序設(shè)計(jì)具有很大的靈活性，實(shí)現(xiàn)較為復(fù)雜，需要程序員熟練掌握GPU體系結(jié)構(gòu)知識(shí)及其編程模型；CUBLAS函數(shù)庫(kù)內(nèi)部整合了具體的并行和優(yōu)化細(xì)節(jié)，為使用者提供了方便的接口，但對(duì)加速算法存在限制.

本文以MNF降維算法為基礎(chǔ)，為實(shí)現(xiàn)較好的加速性能，分別基于CUBLAS庫(kù)和CUDA 進(jìn)行GPU上的并行算法設(shè)計(jì)，并分析比較其性能.本文結(jié)構(gòu)如下：第一部分對(duì)MNF算法進(jìn)行熱點(diǎn)分析；第二部分提出并實(shí)現(xiàn)了基于CUBLAS庫(kù)的MNF GPU并行算法MNF-L（MNF-on-Library）；第三部分基于CUDA提出并實(shí)現(xiàn)了MNF的GPU映射和優(yōu)化算法MNF-C（MNF-on-CUDA）；第四部分通過(guò)實(shí)驗(yàn)測(cè)試并分析比較MNF-L與MNF-C的并行性能；最后是總結(jié)和展望.

1 MNF算法及熱點(diǎn)分析

1.1 MNF算法概述

MNF算法實(shí)質(zhì)是兩次層疊的主成分變換，第一層用于分離并重新調(diào)節(jié)數(shù)據(jù)中的噪聲，消除波段間的相關(guān)；第二層對(duì)噪聲白化數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)主成分變換.

用X={X1，X2，X3，…，Xs}={Y1，Y2，Y3，…，YB}T表示高光譜影像數(shù)據(jù)，用S（W×H，寬W、高H）表示圖像像元數(shù)目，B表示波段數(shù).高光譜影像實(shí)質(zhì)是由B幅W×H圖像組成的三維立體模型.在實(shí)際處理中，X可用B行S列的二維矩陣表示.通過(guò)MNF變換，提取主要的m個(gè)特征波段（m

Step1. 對(duì)高光譜矩陣X濾波；

Step2. 計(jì)算濾波噪聲的協(xié)方差矩陣CN；

Step3. 對(duì)CN特征分解，

DN為降序排列的特征值，U為對(duì)應(yīng)的特征向量，得到變換矩陣

Step4. 求原始高光譜X的協(xié)方差矩陣CD；

Step5. 對(duì)CD進(jìn)行變換：

Step6. 對(duì)CD-adj特征值分解，

其中DD-adj為降序排列的特征值，V為所對(duì)應(yīng)的特征向量；

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

實(shí)現(xiàn)高光譜影像串行MNF降維算法，對(duì)W=614，H=1 087，B=224數(shù)據(jù)降維，測(cè)試并統(tǒng)計(jì)各步驟占總計(jì)算時(shí)間的比例（圖1）.圖1數(shù)據(jù)顯示， Step2和Step4中協(xié)方差矩陣計(jì)算占比最大，超過(guò)80%，Step1的濾波和Step8的MNF變換耗時(shí)較為明顯，上述4個(gè)步驟是MNF算法并行設(shè)計(jì)和優(yōu)化的熱點(diǎn).

2 基于CUBLAS的MNF并行算法

為簡(jiǎn)化GPU移植編碼，nVidia引入經(jīng)典BLAS（Basic Linear Algebra Subprograms）庫(kù)的CUDA實(shí)現(xiàn)版本——CUBLAS庫(kù).CUBLAS庫(kù)充分利用了GPU體系結(jié)構(gòu)，從底層匯編實(shí)現(xiàn)了BLAS庫(kù)函數(shù)的高效運(yùn)算，是在GPU上進(jìn)行矩陣運(yùn)算的最佳選擇之一，憑借其在矩陣乘法中突出的性能表現(xiàn)而被廣泛應(yīng)用[16-17].本節(jié)基于CUBLAS庫(kù)設(shè)計(jì)和實(shí)現(xiàn)協(xié)方差矩陣計(jì)算和MNF變換，并提出MNF-L算法.

2.1 基于 CUBLAS庫(kù)實(shí)現(xiàn)協(xié)方差矩陣計(jì)算

向量協(xié)方差計(jì)算及分解變形如公式（7）：

通過(guò)式（7）變形，將協(xié)方差計(jì)算轉(zhuǎn)換為向量求和∑X與向量?jī)?nèi)積∑XY運(yùn)算，其中向量?jī)?nèi)積可轉(zhuǎn)變?yōu)檩斎敫吖庾V矩陣與其轉(zhuǎn)置矩陣乘積過(guò)程.向量和與向量?jī)?nèi)積運(yùn)算可用CUBLAS庫(kù)實(shí)現(xiàn)：

1）調(diào)用CUBLAS庫(kù)中的cublasSasum函數(shù)求各波段向量和，循環(huán)調(diào)用B次.

2）調(diào)用cublasSsyrk（單精度）求各波段向量?jī)?nèi)積.

cublasSsyrk 接口實(shí)現(xiàn)的功能為C = alpha*A*AT + beta*C，令alpha= 1.0f，beta = 0.0f，輸入矩陣A為高光譜影像矩陣X.

將1）和2）的計(jì)算結(jié)果傳回CPU端，并根據(jù)公式（7），在CPU端串行計(jì)算協(xié)方差矩陣中各元素值.

2.2 采用CUBLAS庫(kù)實(shí)現(xiàn)MNF變換

MNF變換Z=TTX，實(shí)質(zhì)為矩陣乘法運(yùn)算，因此可直接調(diào)用CUBLAS庫(kù)中矩陣乘函數(shù)cublasSgemm（單精度）加速M(fèi)NF變換.

接口cublasSgemm 實(shí)現(xiàn)的功能為C = alpha*A*B + beta*C，為了完成Z=TTX 的功能，令alpha= 1.0f，beta = 0.0f，其中輸入?yún)?shù)A為變換矩陣TT，參數(shù)B為高光譜影像X.

2.3 基于CUBLAS庫(kù)的MNF-L （MNF-on-Library）算法

上述兩節(jié)分別介紹了使用CUBLAS庫(kù)函數(shù)加速協(xié)方差矩陣計(jì)算和MNF變換的過(guò)程，而串行算法的另一個(gè)加速熱點(diǎn)——濾波，在CUBLAS庫(kù)函數(shù)中，沒(méi)有對(duì)應(yīng)的函數(shù)可直接加速.因此在MNF-L算法中，該步驟采用3.2節(jié)所述CUDA并行方案.根據(jù)上述方案，整合各熱點(diǎn)優(yōu)化過(guò)程，提出MNF-L并行算法，如圖2.

在該算法過(guò)程中，多次調(diào)用CUBLAS庫(kù)函數(shù)，由于函數(shù)接口只定義了浮點(diǎn)類型的數(shù)據(jù)運(yùn)算，輸入的高光譜數(shù)據(jù)存儲(chǔ)類型為uchar，因此必須轉(zhuǎn)換為浮點(diǎn)數(shù)據(jù)類型后，才能進(jìn)行計(jì)算，引入數(shù)據(jù)類型轉(zhuǎn)換開銷的同時(shí)，增加了數(shù)據(jù)傳輸時(shí)間.

3 基于CUDA的MNF并行算法

利用CUDA編程模型設(shè)計(jì)GPU并行算法時(shí)，需要設(shè)計(jì)者根據(jù)算法特點(diǎn)對(duì)GPU線程組織層次進(jìn)行設(shè)計(jì)和存儲(chǔ)優(yōu)化等，具有很大的靈活性，其優(yōu)化效果很大程度上取決于設(shè)計(jì)者的映射方法，訪存設(shè)計(jì)等諸多細(xì)節(jié).本節(jié)基于CUDA C，對(duì)串行算法中的3個(gè)熱點(diǎn)分別進(jìn)行了映射方案設(shè)計(jì)和優(yōu)化，實(shí)現(xiàn)MNF并行算法.

3.1 協(xié)方差矩陣計(jì)算并行方案與優(yōu)化

協(xié)方差矩陣中的元素表示各向量之間的協(xié)方差，兩個(gè)向量的協(xié)方差計(jì)算公式見2.1節(jié)公式（7）.

其中參與協(xié)方差計(jì)算的是各波段所有像元組成的向量，數(shù)據(jù)量較大，可通過(guò)拆分實(shí)現(xiàn)并行計(jì)算；另協(xié)方差矩陣中i行j列元素的計(jì)算只需波段i和波段j的數(shù)據(jù)，因此協(xié)方差矩陣中各元素的計(jì)算相互獨(dú)立，可并行計(jì)算；這就使得協(xié)方差矩陣計(jì)算存在兩層并行.

3.1.1 GPU上協(xié)方差矩陣計(jì)算的映射方案

GPU包含grid-block-thread 3個(gè)線程組織層次，本文針對(duì)協(xié)方差矩陣運(yùn)算過(guò)程提出3種GPU映射方案.

方案1 GPU中每個(gè)thread負(fù)責(zé)協(xié)方差矩陣中一個(gè)元素的計(jì)算，如圖3所示，啟動(dòng)B個(gè)block，每個(gè)block啟動(dòng)B個(gè)thread，實(shí)質(zhì)是將grid中所有線程組織為B×B大小的矩陣，以此對(duì)應(yīng)協(xié)方差矩陣中每個(gè)元素的計(jì)算任務(wù).由于協(xié)方差矩陣的對(duì)稱性，可先根據(jù)線程id和線程塊id的大小啟動(dòng)下三角（上三角）位置的線程進(jìn)行計(jì)算，最后將結(jié)果矩陣對(duì)稱位置補(bǔ)齊即可.

方案2 GPU上每個(gè)block完成協(xié)方差矩陣中一個(gè)元素的計(jì)算（圖4）.采用二維結(jié)構(gòu)組織線程塊，即啟動(dòng)B×B大小的block方陣，各block分別對(duì)應(yīng)協(xié)方差矩陣中相同位置元素的計(jì)算任務(wù).同方案1類似，采用對(duì)稱補(bǔ)償?shù)姆椒p少實(shí)際參與計(jì)算的線程塊數(shù)目.

方案3 GPU中每個(gè)grid負(fù)責(zé)協(xié)方差矩陣中一個(gè)元素的計(jì)算，即每啟動(dòng)一個(gè)kernel函數(shù)完成協(xié)方差矩陣中一個(gè)元素的計(jì)算，循環(huán) （1+B）*B/2次完成協(xié)方差矩陣中所有元素的計(jì)算.由于各波段數(shù)據(jù)量（S=W×H）龐大，協(xié)方差矩陣中單個(gè)元素計(jì)算涉及高維向量加與內(nèi)積運(yùn)算，將其映射到grid層可減小并行粒度，且成像波段數(shù)B一般為224，使循環(huán)次數(shù)控制在可接受范圍內(nèi).

上述3種方案的本質(zhì)區(qū)別在于協(xié)方差矩陣中單個(gè)元素計(jì)算的映射層次不同，后續(xù)將針對(duì)方案1展開優(yōu)化研究，根據(jù)文獻(xiàn)[8]中相關(guān)內(nèi)容對(duì)方案2，3進(jìn)行優(yōu)化，后文3.1.3節(jié)將比較3種方案的執(zhí)行效果.

3.1.2 GPU上協(xié)方差矩陣計(jì)算的共享存儲(chǔ)優(yōu)化

CUDA線程可以訪問(wèn)不同層次存儲(chǔ)空間的數(shù)據(jù)，如全局，本地，共享，常量或紋理等，不同層次的存儲(chǔ)器訪問(wèn)速度不同.共享存儲(chǔ)器位于片上存儲(chǔ)，同一個(gè)線程塊內(nèi)的線程可以共享訪問(wèn)，其訪存速度比全局訪存快一個(gè)數(shù)量級(jí).

1）方案1中，同一block中的線程存在數(shù)據(jù)復(fù)用，如線程塊i中的線程重復(fù)讀取第i波段的影像數(shù)據(jù)，圖5所示，將該波段數(shù)據(jù)存儲(chǔ)到共享存儲(chǔ)器中，block內(nèi)所有線程直接訪問(wèn)共享存儲(chǔ)，最多可減少（B-1）*S次全局訪存開銷.

基于2.1中式（7）對(duì)協(xié)方差計(jì)算的改寫，在方案1映射結(jié)構(gòu)上予以實(shí)現(xiàn)，其中矩陣（與其轉(zhuǎn)置矩陣）相乘的過(guò)程可進(jìn)行共享存儲(chǔ)優(yōu)化，如圖6顯示，將數(shù)據(jù)分塊讀取至共享存儲(chǔ)單元，進(jìn)行分塊迭代.建立等同于線程塊大?。╰hx*thx）的共享存儲(chǔ)區(qū)，每次讀取對(duì)應(yīng)塊到共享存儲(chǔ)，完成分塊相乘、相加；共享存儲(chǔ)向右滑動(dòng)，重復(fù)上述過(guò)程，直到數(shù)據(jù)末尾.

2）從全局存儲(chǔ)讀取數(shù)據(jù)至共享存儲(chǔ)時(shí)，將右乘矩陣按照轉(zhuǎn)置后的位置進(jìn)行保存，使矩陣行列相乘轉(zhuǎn)變?yōu)樾行邢喑耍苊庾x取共享存儲(chǔ)時(shí)出現(xiàn)bank conflict，提高矩陣相乘的效率.

根據(jù)上述思想在原始方案1基礎(chǔ)上進(jìn)行優(yōu)化，并采用4組數(shù)據(jù)實(shí)驗(yàn)對(duì)比不同優(yōu)化算法的執(zhí)行時(shí)間，圖7所示，優(yōu)化后的執(zhí)行時(shí)間比原始版本低1～2個(gè)數(shù)量級(jí)，同時(shí)使用兩類優(yōu)化方法的加速效果更為顯著；說(shuō)明本文采取的共享存儲(chǔ)優(yōu)化策略是非常有效的.

3.1.3 協(xié)方差矩陣計(jì)算不同方案的性能比較

測(cè)定4組高光譜數(shù)據(jù)在3.1.1節(jié)中3種映射方案下（優(yōu)化后）的協(xié)方差矩陣計(jì)算時(shí)間（不包括數(shù)據(jù)通信時(shí)間），見圖8，其中方案1采用3.1.2節(jié)兩種優(yōu)化方法（記為G-cov），在3種方案中性能最佳，說(shuō)明本文選擇的設(shè)計(jì)方案和進(jìn)行的優(yōu)化改進(jìn)具有顯著的成效.

3.2 噪聲估計(jì)（濾波）的并行映射與優(yōu)化

濾波是對(duì)像素點(diǎn)進(jìn)行處理以得到相應(yīng)點(diǎn)的噪聲估計(jì)值，該步驟需要目標(biāo)點(diǎn)及周圍8個(gè)點(diǎn)參與運(yùn)算.圖像各點(diǎn)的濾波計(jì)算相互獨(dú)立，不存在相互依賴，因此圖像中所有像元點(diǎn)都能并行計(jì)算，且高光譜圖像不同波段也能并行計(jì)算.在GPU中，每個(gè)線程可用來(lái)計(jì)算一個(gè)像元的噪聲估計(jì)值.

3.2.1 噪聲估計(jì)（濾波）的并行映射方案

GPU上每個(gè)線程完成高光譜矩陣中一個(gè)元素的濾波任務(wù)，如圖9（記為method1），邊界噪聲事先置零，不參與映射，啟動(dòng)H-2個(gè)線程塊，每個(gè)線程塊內(nèi)啟動(dòng)W-2個(gè)thread，將噪聲矩陣非邊界元素的濾波任務(wù)映射到相應(yīng)的線程中進(jìn)行計(jì)算，每次完成一個(gè)波段的映射，每個(gè)線程循環(huán)計(jì)算B次，完成所有波段元素的濾波任務(wù).

改變method1中線程塊及線程的組織方式，均采用二維分布，將噪聲矩陣分塊映射，得到映射圖10（記為method2）.

3.2.2 濾波的GPU細(xì)粒度優(yōu)化

針對(duì)濾波過(guò)程，主要采取以下兩種優(yōu)化方法：

1）利用寄存器替換本地存儲(chǔ)訪問(wèn).

原執(zhí)行函數(shù)采用中間值濾波，是對(duì)包括原濾波點(diǎn)在內(nèi)的周圍9個(gè)點(diǎn)進(jìn)行排序，選取中間的點(diǎn)與原濾波點(diǎn)做差，在這個(gè)過(guò)程中，將使用本地存儲(chǔ)器，而無(wú)法使用寄存器，本地存儲(chǔ)實(shí)際存儲(chǔ)在顯存中，訪存延遲大，因此性能較差.

改進(jìn)濾波函數(shù)，采用均值濾波，即將周圍所有點(diǎn)相加求得的平均值與原濾波點(diǎn)比對(duì)，此時(shí)，無(wú)需本地存儲(chǔ)器，而僅需要寄存器參與運(yùn)算，可大大加快濾波過(guò)程.經(jīng)驗(yàn)證，改進(jìn)的濾波函數(shù)不影響處理精度，因此可采用改進(jìn)方法替代原方法.

2）利用共享存儲(chǔ)減少?gòu)?fù)用數(shù)據(jù)的全局訪存.

在濾波過(guò)程中，非邊界的點(diǎn)會(huì)被相鄰的9個(gè)點(diǎn)濾波使用，因此存在9次復(fù)用.利用共享存儲(chǔ)訪存快、線程塊內(nèi)共享的屬性，可以先將線程塊所需要的數(shù)據(jù)讀入共享存儲(chǔ)，各線程需要時(shí)從共享存儲(chǔ)讀取數(shù)據(jù)，此時(shí)只需訪問(wèn)一次全局存儲(chǔ)，減少了8次全局存儲(chǔ)訪問(wèn).

圖9，圖10的方案中，由于block組織方式不同，采取的共享存儲(chǔ)數(shù)據(jù)也有所差異.

method1中每個(gè)線程塊負(fù)責(zé)一行數(shù)據(jù)的濾波，因此，將該行與前、后兩行讀取到共享存儲(chǔ)區(qū)（圖11）.共享存儲(chǔ)大小為3*S個(gè)字節(jié)，線程塊內(nèi)所有線程并行讀取一次（個(gè)別線程讀取兩次）可完成本地到共享空間的轉(zhuǎn)儲(chǔ).

圖12顯示，將method2中各矩陣塊及其周圍元素?cái)?shù)據(jù)讀取到該數(shù)據(jù)塊濾波映射的線程塊中，使該block內(nèi)所有線程可以直接從共享存儲(chǔ)區(qū)讀取計(jì)算數(shù)據(jù).

3.2.3 不同優(yōu)化下的濾波計(jì)算性能對(duì)比

以W=614，H=1 087，B=224的高光譜數(shù)據(jù)為輸入，對(duì)比不同方法的執(zhí)行時(shí)間，圖13前兩組數(shù)據(jù)為method1分別采用中間值濾波和均值濾波的kernel執(zhí)行時(shí)間（不包括通信時(shí)間），結(jié)果表明改進(jìn)后執(zhí)行時(shí)間減少了一半以上，均值濾波中寄存器訪存加速效果顯著.

method1，method2均同時(shí)采用均值濾波和共享存儲(chǔ)優(yōu)化進(jìn)行實(shí)現(xiàn)，比較兩種方法優(yōu)化后的執(zhí)行時(shí)間，如圖13后兩組數(shù)據(jù)顯示，二者執(zhí)行時(shí)間十分接近，較均值濾波優(yōu)化方法提高3～4 ms，加速比例為11.7%～15.6%，說(shuō)明3.2.2節(jié)討論的兩類優(yōu)化均能使算法取得加速效果.

3.3 MNF變換的GPU映射

2.2節(jié)采用CUBLAS庫(kù)中矩陣乘函數(shù)實(shí)現(xiàn)MNF變換，本節(jié)針對(duì)MNF變換中相乘矩陣規(guī)模的特殊性，在傳統(tǒng)矩陣乘并行方案基礎(chǔ)上進(jìn)行并行設(shè)計(jì)、改進(jìn)存儲(chǔ)策略，以充分發(fā)揮加速潛能.

3.3.1 MNF變換的映射方案

圖14為映射圖，根據(jù)分塊計(jì)算的思想，將結(jié)果矩陣縱向分割，每個(gè)線程塊完成一個(gè)分割塊的計(jì)算，即線程塊采用一維分布，各block負(fù)責(zé)blockDim.x列數(shù)據(jù)的計(jì)算.由于S一般遠(yuǎn)遠(yuǎn)超過(guò)block數(shù)量，因此，每個(gè)block執(zhí)行完一塊計(jì)算后，固定跳步到下一塊計(jì)算，直到所有分塊完成.

每個(gè)線程塊采用二維結(jié)構(gòu)組織線程，每個(gè)線程負(fù)責(zé)相應(yīng)位置元素的計(jì)算，如果m（波段數(shù)）大于blockDim.y，各線程最多需循環(huán)計(jì)算（m+blockDim.y-1）/blockDim.y次，通常數(shù)據(jù)降維后波段數(shù)小于32，設(shè)置blockDim.y為16，使每個(gè)block縱向最多映射2次.

3.3.2 MNF變換的優(yōu)化策略和性能比較

單個(gè)線程中實(shí)質(zhì)進(jìn)行的是向量乘法運(yùn)算，需要不斷的讀取變換矩陣各行數(shù)據(jù)與高光譜影像各列數(shù)據(jù)，直接的global 訪存延遲大，效率低，因此考慮使用共享存儲(chǔ)降低讀取數(shù)據(jù)帶來(lái)的延遲，而考慮到變換矩陣規(guī)模較?。ㄒ话阈∮?2×224），沒(méi)有超出常量存儲(chǔ)器的存儲(chǔ)容量（64 kB），因此可直接將變換矩陣整體保存在常量存儲(chǔ)中，進(jìn)一步提升訪問(wèn)速度的同時(shí)，避免各線程重復(fù)讀取數(shù)據(jù)到共享存儲(chǔ).

因此在該步驟中進(jìn)行了如下優(yōu)化：

1）采用常量存儲(chǔ)器存儲(chǔ)變換矩陣；

2）使用共享存儲(chǔ)器存儲(chǔ)原始圖像B×blockDim.x大小的塊數(shù)據(jù).

根據(jù)不同的優(yōu)化方式實(shí)現(xiàn)了4類優(yōu)化版本：

v0）全局訪存，無(wú)共享存儲(chǔ)優(yōu)化；

v1）僅使用共享存儲(chǔ)優(yōu)化；

v2）僅使用常量存儲(chǔ)優(yōu)化；

v3）同時(shí)使用共享存儲(chǔ)和常量存儲(chǔ).

圖15展示了4類版本的性能對(duì)比，4組數(shù)據(jù)測(cè)試下的性能比較結(jié)果保持一致，即優(yōu)化效果v3>v1>v2>v0，v3性能最好，說(shuō)明采用共享存儲(chǔ)和常量存儲(chǔ)均能取得優(yōu)化效果；隨著X數(shù)據(jù)量增加，v1較v2優(yōu)勢(shì)漸增，由于v2僅對(duì)變換矩陣使用常量存儲(chǔ)優(yōu)化，其全局訪存主要來(lái)自X，隨X數(shù)據(jù)量增大，v2全局訪存次數(shù)將顯著增加.

3.4 基于CUDA C的MNF-C（MNF-on-CUDA）

算法

參照前3節(jié)中關(guān)于協(xié)方差矩陣計(jì)算、濾波和MNF變換的GPU并行映射方案和優(yōu)化效果，以最優(yōu)方案替代串行算法中的相應(yīng)步驟，整合提出MNF-C算法，圖16為其流程圖.

其中，第3個(gè)GPU kernel執(zhí)行期間，CPU端計(jì)算與GPU端協(xié)方差矩陣計(jì)算可同時(shí)進(jìn)行，實(shí)現(xiàn)了Host（CPU）與device（GPU）計(jì)算重疊，充分發(fā)揮了異構(gòu)系統(tǒng)的并行優(yōu)勢(shì).

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

4.1 實(shí)驗(yàn)平臺(tái)和測(cè)試數(shù)據(jù)集

本文使用CUDA C實(shí)現(xiàn)上述算法.實(shí)驗(yàn)平臺(tái)為GPU微型超算平臺(tái)，包含2個(gè)8核的Intel（R） Xeon（R） CPU E5-2670和nVidia Kepler架構(gòu)的Tesla K20c GPU，采用Intel（R） C Compiler XE 13.0.0.079和CUDA5.5工具包進(jìn)行編譯.

實(shí)驗(yàn)采用了4組AVIRIS高光譜數(shù)據(jù)，數(shù)據(jù)已經(jīng)經(jīng)過(guò)預(yù)處理，將光譜數(shù)據(jù)轉(zhuǎn)換為unsigned char類型存儲(chǔ)在二進(jìn)制文件中.表1詳細(xì)列出了測(cè)試集的信息.

4.2 MNF-L與MNF-C性能對(duì)比

4.2.1 協(xié)方差矩陣計(jì)算性能對(duì)比

對(duì)比MNF-L中采用CUBLAS庫(kù)函數(shù)和MNF-C中程序語(yǔ)言實(shí)現(xiàn)協(xié)方差矩陣的執(zhí)行效率，采用4組數(shù)據(jù)計(jì)算總執(zhí)行時(shí)間，并分別標(biāo)識(shí)運(yùn)算、通信、創(chuàng)建銷毀等不同階段.如圖17所示，CUDA 實(shí)現(xiàn)版本（記為G-cov）同CUBLAS版本（記為L(zhǎng)-cov）總執(zhí)行時(shí)間幾乎相同，各有兩次稍占優(yōu)勢(shì)，但整體相差細(xì)微.

高光譜數(shù)據(jù)集&實(shí)現(xiàn)方案

G-cov與L-cov執(zhí)行熱點(diǎn)主要集中在cov運(yùn)算和數(shù)據(jù)傳輸，但最耗時(shí)階段形成鮮明對(duì)比，G-cov主要集中在運(yùn)算階段，通信開銷較小，L-cov則與之相反.CUBLAS庫(kù)接口定義的數(shù)據(jù)類型均為浮點(diǎn)型，而高光譜數(shù)據(jù)初始為unsigned char，運(yùn)算前需要進(jìn)行類型轉(zhuǎn)換，數(shù)據(jù)量增加3倍，導(dǎo)致巨大通信開銷.

4.2.2 MNF變換性能對(duì)比

見表2，v3（3.3節(jié)）在MNF-C算法中采用的CUDA C實(shí)現(xiàn)的MNF變換，對(duì)比MNF-L中CUBLAS庫(kù)函數(shù)實(shí)現(xiàn)的執(zhí)行時(shí)間.

CUBLAS庫(kù)函數(shù)加速效果優(yōu)于v3版本，說(shuō)明CUDA C的設(shè)計(jì)和優(yōu)化仍存在加速空間；但在實(shí)驗(yàn)中發(fā)現(xiàn)CUBLAS首次啟動(dòng)需要約200 ms的啟動(dòng)開銷，大于CUDA C中的kernel的啟動(dòng)開銷，如果只進(jìn)行一次矩陣乘法，將難以發(fā)揮其運(yùn)算優(yōu)勢(shì).

同時(shí)隨數(shù)據(jù)規(guī)模不斷增加，CUBLAS的缺點(diǎn)逐漸顯現(xiàn)，如通信開銷、存儲(chǔ)限制等，4.4節(jié)將進(jìn)行具體分析.

4.2.3 并行算法加速比

實(shí)驗(yàn)分別測(cè)試了表1中4組數(shù)據(jù)的串行MNF降維、MNF-L和MNF-C并行降維時(shí)間（表3），計(jì)算MNF-L與MNF-C相對(duì)于串行MNF的加速比（圖18）.圖中數(shù)據(jù)顯示，基于CUBLAS庫(kù)實(shí)現(xiàn)的MNF-L算法可加速11.5～60.6倍不等；MNF-C算法加速效果最好，可加速46.5～92.9倍不等.

文獻(xiàn)[14]對(duì)MNF算法中協(xié)方差矩陣計(jì)算進(jìn)行了熱點(diǎn)加速，使用3組高光譜數(shù)據(jù)進(jìn)行測(cè)試，實(shí)驗(yàn)顯示，隨高光譜圖像數(shù)據(jù)量增大，并行加速效果增加，當(dāng)最大數(shù)據(jù)集為1 036 × 235 × 229時(shí)，加速比為4.58.而本文設(shè)計(jì)實(shí)現(xiàn)的兩種算法各有特色，較文獻(xiàn)[14]，采用了更全面的優(yōu)化策略和測(cè)試數(shù)據(jù)集，不僅協(xié)方差矩陣計(jì)算實(shí)現(xiàn)了更好的加速比，還發(fā)掘出濾波、MNF變換等步驟的加速潛能，加速效果提升了數(shù)十倍，可應(yīng)用于高光譜等較大規(guī)模數(shù)據(jù)的特征提取，將大大提高執(zhí)行效率.

4.3 實(shí)驗(yàn)討論

輸入不同測(cè)試數(shù)據(jù)，計(jì)算MNF-C和MNF-L算法執(zhí)行時(shí)間比值（圖19），實(shí)驗(yàn)結(jié)果顯示，MNF-C算法加速效果優(yōu)于MNF-L算法，隨數(shù)據(jù)集不斷增大，兩種算法性能差距逐漸縮小.

高光譜數(shù)據(jù)集

根據(jù)算法特點(diǎn)，對(duì)MNF-L與MNF-C進(jìn)行比較分析：

1）MNF-L中CUBLAS庫(kù)函數(shù)計(jì)算前需要先將unsigned char數(shù)據(jù)轉(zhuǎn)換為float型，引入了數(shù)據(jù)轉(zhuǎn)換和通信開銷.

2）MNF-C可顯式調(diào)整算法步驟，實(shí)現(xiàn)host與device的計(jì)算重疊，而MNF-L中的CUBLAS庫(kù)函數(shù)隱藏相關(guān)細(xì)節(jié)，無(wú)法實(shí)現(xiàn)不同設(shè)備間的協(xié)同計(jì)算.

3）在協(xié)方差矩陣計(jì)算部分，兩算法加速效果相近.雖然MNF變換中CUBLAS庫(kù)函數(shù)頗具優(yōu)勢(shì)，但該步驟所占時(shí)間比重較小，使得加速程度受限，加之CUBLAS庫(kù)函數(shù)的首次啟動(dòng)開銷，一定程度上削弱了MNF-L算法的優(yōu)勢(shì).

4）隨數(shù)據(jù)量增大，MNF-L加速效果逐漸提升，運(yùn)算優(yōu)勢(shì)一定程度彌補(bǔ)了通信開銷，但對(duì)于更大規(guī)模的測(cè)試數(shù)據(jù)，所需存儲(chǔ)空間容量至少為高光譜數(shù)據(jù)的 4倍（使用雙精度接口函數(shù)達(dá)8倍），將導(dǎo)致顯存不足.

兩類算法各有特點(diǎn)，可應(yīng)用于不同的加速需求領(lǐng)域.基于CUBLAS庫(kù)設(shè)計(jì)的并行算法，實(shí)現(xiàn)簡(jiǎn)單，代碼簡(jiǎn)潔，隱藏了相關(guān)算法的實(shí)現(xiàn)細(xì)節(jié)，對(duì)程序員是透明的，較容易掌握，可用于復(fù)雜和編碼量較大的算法.除了CUBLAS庫(kù)，CUDA還提供了一系列GPU加速庫(kù)函數(shù)，如針對(duì)計(jì)算機(jī)視覺(jué)和計(jì)算流體力學(xué)的cusolverDN庫(kù)、應(yīng)用于化學(xué)反應(yīng)動(dòng)力學(xué)的cusolverSP、cusolverRF庫(kù)，應(yīng)用于深度學(xué)習(xí)的cuDNN等，此類科學(xué)計(jì)算方法一般算法復(fù)雜，為簡(jiǎn)化實(shí)現(xiàn)，采用GPU庫(kù)函數(shù)加速是最為直接和快速的方法.而對(duì)于GPU加速庫(kù)函數(shù)難以實(shí)現(xiàn)的功能或者對(duì)時(shí)效性要求較高的高端應(yīng)用領(lǐng)域，需要獲得盡可能高的性能加速比， 那么CUDA并行設(shè)計(jì)則成為必要選擇.

總之，并行程序設(shè)計(jì)必須以運(yùn)算結(jié)果正確為前提，根據(jù)算法特點(diǎn)、應(yīng)用場(chǎng)景、性能要求等多方面合理選擇實(shí)現(xiàn)方案.

5 結(jié) 論

本文基于CPU/GPU異構(gòu)系統(tǒng)研究了高光譜線性降維方法MNF的并行實(shí)現(xiàn)技術(shù)，提出了切實(shí)有效的MNF并行算法.分別針對(duì)協(xié)方差矩陣計(jì)算、濾波、MNF變換等熱點(diǎn)進(jìn)行GPU并行加速和優(yōu)化研究，設(shè)計(jì)并實(shí)現(xiàn)了基于CUBLAS庫(kù)的MNF-L算法和基于CUDA的MNF-C并行降維算法.實(shí)驗(yàn)結(jié)果顯示MNF-L算法加速11.5～60.6倍，MNF-C加速效果較好，獲得了46.5～92.9的加速比.最后對(duì)基于CUBLAS和CUDA實(shí)現(xiàn)的方法在性能、算法、應(yīng)用領(lǐng)域等方面進(jìn)行了分析討論.

針對(duì)本文應(yīng)用領(lǐng)域——高光譜遙感數(shù)據(jù)降維，后續(xù)工作將圍繞可較好體現(xiàn)高維空間結(jié)構(gòu)的非線性降維算法展開并行優(yōu)化的研究.

參考文獻(xiàn)

[1] 趙英時(shí).遙感應(yīng)用分析原理與方法[M].北京：科學(xué)出版社，2003：5-30.

ZHAO Yingshi.The principle and method of analysis of remote sensing application[M].Beijing：Science Press，2003：5-30.（In Chinese）

[2] 張良培，張立福.高光譜遙感[M].武漢：武漢大學(xué)出版社，2005：30-88.

ZHANG Liangpei，ZHANG Lifu. Hyperspectral remote sensing[M].Wuhan： Wuhan University Press，2005：30-88. （In Chinese）

[3] HUGHES G. On the mean accuracy of statistical pattern recognizers[J]. IEEE Transactions on Information Theory， 1968， 14（1）：55-63.

[4] UTO K， KOSUGI Y， SAITO G.Semi-supervised hyperspectral subspace learning based on a generalized eigenvalue problem for regression and dimensionality reduction[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2014，7（6）： 2583-2599.

[5] JIA Xiuping， JOHN A Richards. Efficient transmission and classification of hyperspectral image data[J]. IEEE Transactions on Geoscience and Remote Sensing，2003，41（5）：1129-1131.

[6] 張見明，余列祥，劉路平.基于GPU 加速的邊界面法正則積分的研究[J].湖南大學(xué)學(xué)報(bào)：自然科學(xué)版，2013，40（3）：41-45.

ZHANG Jianming，YU Liexiang，LIU Luping. Research on regular integration in boundary face method based on GPU-acceleration[J]. Journal of Hunan University：Natural Sciences，2013，40（3）：41-45. （In Chinese）

[7] TOP500. TOP500 Supercomputer Sites[EB/OL].http：//www.top500.org.

[8] 方民權(quán)， 周海芳， 申小龍. 基于GPU的高光譜PCA降維多級(jí)并行算法[J]. 東北大學(xué)學(xué)報(bào)：自然科學(xué)版， 2014，35（S1）： 238-243.

FANG Minquan， ZHOU Haifang， SHEN Xiaolong. Multilevel parallel algorithm of PCA dimensionality reduction for hyperspectral image on GPU[J]. Journal of Northeastern University：Natural Science， 2014， 35（S1）： 238-243. （In Chinese）

[9] 方民權(quán)， 張衛(wèi)民， 張理論， 等. 面向MIC的高光譜影像降維多級(jí)并行算法及性能優(yōu)化[J]. 東北大學(xué)學(xué)報(bào)：自然科學(xué)版， 2014， 35（S1）： 25-30，36.

FANG Minquan， ZHANG Weimin， ZHANG Lilun， et al. Multilevel parallel algorithms and performance optimization of hyperspectral image dimensionality reduction on MIC[J]. Journal of Northeastern University：Natural Science， 2014， 35（S1）： 25-30，36. （In Chinese）

[10]SERGIO Sánchez，RUI Ramalho，LEONEL Sousa. Antonio Plaza，real-time implementation of remotely sensed hyperspectral image unmixing on GPUs[J].Journal of Real-Time Image Processing，2015，10（3）：469-483.

[11]ELMAGHRBAY M， AMMAR R， RAJASEKARAN S. Fast GPU algorithms for endmember extraction from hyperspectral images[C]// Proceedings of the 2012 IEEE Symposium on Computers and Communications （ISCC 12）.Cappadocia，2012： 000631-000636.

[12]CHANG Y L ， CHEN K S ， HUANG B， et al. A parallel simulated annealing approach to band selection for high-dimensional remote sensing images[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2011，4（3）：579-590.

[13]SANTOS L， MAGLIE， VITULLI R，et al. Highly-parallel GPU architecture for lossy hyperspectral image compression[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing， 2013，6（2）：670-681.

[14]羅耀華，郭科，趙仕波.基于GPU 的高光譜遙感MNF并行方法研究[J]. 四川師范大學(xué)學(xué)報(bào)：自然科學(xué)版， 2013，36（3）：476-479.

LUO Yaohua，GUO Ke，ZHAO Shibo.Minimum noise fraction of hyperspectral remote sensing in parallel computing based on GPU[J]. Journal of Sichuan Normal University：Natural Science，2013， 36（3）：476-479. （In Chinese）

[15]LIU Xiang，ZHANG Bing，GAO Lianru，et al. A maximum noise fraction transform with improved noise estimation for hyperspectral images[J].Science in China， Series F： Information Sciences，2009（9）：1578-1587.

[16]MANUEL Carcenac.From tile algorithm to stripe algorithm： a CUBLAS-based parallel implementation on GPUs of Gauss method for the resolution of extremely large dense linear systems stored on an array of solid state devices[J].The Journal of Supercomputing，2014，68（1）：365-413.

[17]ZHANG Bo，YANG Xiang，YANG Fei，et al. The CUBLAS and CULA based GPU acceleration of adaptive finite element framework for bioluminescence tomography[J].Optics Express，2010， 18（19）：20201-20214.