談繼魁，方 勇，霍迎秋

(西北農(nóng)林科技大學(xué)信息工程學(xué)院，陜西楊陵712100)

隨著信息技術(shù)、計(jì)算機(jī)技術(shù)及多媒體技術(shù)的快速發(fā)展，需要處理的數(shù)據(jù)量越來(lái)越大，而傳統(tǒng)的Nyquist采樣定理要求對(duì)信號(hào)的采樣率不得低于信號(hào)帶寬的2倍，這就需要有效的方法來(lái)提高對(duì)信號(hào)的處理能力。壓縮感知(Compressed Sensing，CS)理論是Donobo與Candes等人在2006年提出的一種新的信號(hào)處理理論框架［1］。該理論表明，利用較少的信號(hào)值就可實(shí)現(xiàn)稀疏的或可壓縮信號(hào)的精確重建，這使得其在信號(hào)處理領(lǐng)域中有著非常廣闊的應(yīng)用前景［2］。重建算法是CS理論的核心。在匹配追蹤(Matching Pursuit，MP)算法的基礎(chǔ)上提出來(lái)的正交匹配追蹤(Orthogonal Matching Pursuit，OMP)算法是經(jīng)典的重建算法［3］，其沿用了MP算法選擇原子的準(zhǔn)則，不同的是OMP算法通過(guò)遞歸地對(duì)已選擇的所有原子進(jìn)行正交化處理，以保證每次迭代的最優(yōu)性，從而減少了迭代次數(shù)［4］。

OMP算法提高了算法的性能，但在分解的每一步都對(duì)所選原子做正交化處理，增加了算法的計(jì)算復(fù)雜度［5］。針對(duì)這一問(wèn)題，本文采用并行計(jì)算在圖形處理單元(Graphics Processing Unit，GPU)上加速該算法，首先分析OMP算法的復(fù)雜度，找出影響其復(fù)雜度較高的部分，在GPU平臺(tái)上設(shè)計(jì)其并行方案對(duì)算法進(jìn)行優(yōu)化，以提高OMP算法的運(yùn)行速度。

1 OMP算法分析

令傳感矩陣 Α=(α1，…，αn)，其中 αi∈Rm(R 表示實(shí)數(shù)集)是矩陣Α的第i列，稱(chēng)矩陣Α和向量αi分別為字典和原子。令y=Az，其中z在Ψ域(可逆變換矩陣)中k階稀疏［6］。當(dāng)?shù)螖?shù)t≤k時(shí)，OMP算法計(jì)算步驟如下:

第一步，輸入傳感矩陣 A∈Rm×n，測(cè)量值 y∈Rm，稀疏度為k;令殘差r=y，被選中原子的索引集i=Φ;

第三步，由i=i∪it更新索引集;

按照上述幾步，對(duì)OMP算法的復(fù)雜度作如下分析:

1)投影:投影部分的實(shí)質(zhì)是矩陣向量的乘法運(yùn)算。殘差r在原子上的投影可以表述為〈r，ai'〉/ai'2，〈·，·〉表示兩個(gè)向量的內(nèi)積，·2表示向量的l2模。為簡(jiǎn)化計(jì)算，令ρ =(a12，…，an2)T，則投影部分可以轉(zhuǎn)化為 ATr./ρ，其中．/表示點(diǎn)除。由于矩陣Α的大小為m×n而殘差r是大小為m的列向量，因此，投影部分的復(fù)雜度為Ο(m×n)。

2)選擇最匹配的原子:從未選中的原子中選取與殘差r投影絕對(duì)值最大原子的索引。字典Α每一行有n個(gè)原子，且k?n，故本部分的復(fù)雜度為 Ο(n)。

4)更新殘差:由于t≤k?n2，因此相比于投影部分的復(fù)雜度，這一步的復(fù)雜度非常小甚至可以忽略。

基于上述分析可知，當(dāng)k值比較小的時(shí)候，OMP算法的計(jì)算復(fù)雜度主要取決于投影部分，且其復(fù)雜度為Ο(m×n);而當(dāng)k取值較大的時(shí)候，OMP算法的計(jì)算復(fù)雜度主要取決于矩陣求逆部分，其復(fù)雜度為 Ο(k3)。因此，影響OMP算法計(jì)算效率的關(guān)鍵步驟為投影部分矩陣向量相乘和更新權(quán)值部分矩陣求逆。

2 GPU簡(jiǎn)述

GPU是由流多處理器(Streaming Multiprocessor，SM)構(gòu)成的，而流多處理器是由流處理器(Streaming Processor，SP)構(gòu)成的。運(yùn)行在GPU上的程序稱(chēng)為核函數(shù)。圖1所示為GPU編程模型，GPU線程被組織成Thread-Block-Grid三層結(jié)構(gòu)，每個(gè) Grid包含若干 Block，每個(gè) Block又包含若干個(gè)Thread。GPU運(yùn)行時(shí)采用單指令多線程(Single Instruction，Multiple Thread，SIMT)執(zhí)行模式，這意味著同一個(gè)SM中的SP執(zhí)行的指令完全相同。GPU的存儲(chǔ)器是由1個(gè)全局存儲(chǔ)器和4個(gè)不同類(lèi)型的片上存儲(chǔ)器構(gòu)成，分別為常量存儲(chǔ)器、紋理存儲(chǔ)器、共享存儲(chǔ)器和寄存器［7］。寄存器是線程特定的，即每個(gè)線程只能訪問(wèn)自己的寄存器，而共享存儲(chǔ)器是SM特定的，每個(gè)SM的所有線程都能訪問(wèn)這個(gè)SM的共享存儲(chǔ)器。全局存儲(chǔ)器、常量存儲(chǔ)器和紋理存儲(chǔ)器都可以被所有GPU的線程訪問(wèn)。常量存儲(chǔ)器和紋理存儲(chǔ)器都是只讀的，使用它們可以減少所需要的內(nèi)存帶寬。全局存儲(chǔ)器、常量存儲(chǔ)器和紋理存儲(chǔ)器的使用可以實(shí)現(xiàn)不同的優(yōu)化，例如紋理存儲(chǔ)器可以為特定的數(shù)據(jù)格式提供不同的尋址方式［7-10］。

圖1 GPU編程模型示意圖

本實(shí)驗(yàn)是在NVIDIA GTX480顯卡上完成的，它由15個(gè)有效的SM構(gòu)成，每個(gè)SM由32個(gè)SP構(gòu)成。因此，該顯卡總共含有480個(gè)處理器內(nèi)核(Processor Core)。NVIDIA GTX480的每個(gè)SM包含64 kbyte的寄存器和48 kbyte的共享存儲(chǔ)器，在運(yùn)行時(shí)SM將寄存器平均分配給所有線程。此外，每個(gè)NVIDIA GTX480還包含1.5 Gbyte的全局存儲(chǔ)器和64 kbyte的紋理存儲(chǔ)器。

3 OMP算法在GPU上的并行實(shí)現(xiàn)

通過(guò)分析可知，影響OMP算法復(fù)雜度的關(guān)鍵步驟是矩陣求逆和矩陣向量相乘兩部分。下面分別介紹基于GPU的矩陣求逆和矩陣向量相乘的并行計(jì)算，并給出OMP算法在GPU上并行實(shí)現(xiàn)的方案。

3．1 矩陣逆更新

由式(1)可知求解t×t大小矩陣逆的復(fù)雜度為 Ο(t3)。然而在OMP算法的每次迭代中，At是在上添加一列構(gòu)成的，這表明可以使用矩陣逆更新算法［11］來(lái)求解矩陣的逆從而降低其計(jì)算復(fù)雜度。

式中:d和P的表述如式(3)和式(4)所示

為了進(jìn)一步簡(jiǎn)化計(jì)算，作如下定義

將d和P簡(jiǎn)化為式(7)和式(8)

由此，Ht的表達(dá)式可以簡(jiǎn)化為式(9)

3．2 矩陣向量相乘塊算法

3．3 OMP算法的并行實(shí)現(xiàn)

要在GPU上并行實(shí)現(xiàn)OMP算法，需要在GPU上為輸入變量分配內(nèi)存。輸入變量為和y，輸出變量為z。由前面的分析可知，OMP算法中會(huì)用到原子的范數(shù)，因此定義ρ(長(zhǎng)度為n的向量)用于存放算法中經(jīng)常使用原子的-范數(shù)，定義向量i用于存儲(chǔ)已選定原子的索引。如圖2所示為OMP算法在GPU上并行實(shí)現(xiàn)的流程圖。

圖2 OMP算法GPU并行化實(shí)現(xiàn)流程圖

在OMP算法運(yùn)行前，首先把矩陣AT和向量y中的數(shù)據(jù)從CPU內(nèi)存?zhèn)鬏斨罣PU的全局存儲(chǔ)器;OMP算法執(zhí)行結(jié)束后，將輸出的z從GPU全局存儲(chǔ)器中拷貝到CPU內(nèi)存中。OMP算法的每一次迭代中都會(huì)用到ρ，為了避免頻繁的重新計(jì)算，提前計(jì)算得到ρ并把它保存到GPU的全局存儲(chǔ)器中。

OMP并行過(guò)程包含k次循環(huán)，每次循環(huán)求得一個(gè)原子的索引。在GPU上并行實(shí)現(xiàn)OMP算法方案如下:首先定義一個(gè)核函數(shù)，其中“·* ”表示兩個(gè)向量的點(diǎn)乘。在這個(gè)核函數(shù)中，z用來(lái)區(qū)分原子是否已經(jīng)被選中:如果z中相應(yīng)原子不為零，則表示該原子已被選中，將在這個(gè)原子上的投影值設(shè)置為零，以避免該原子再次被選中。這種方法的優(yōu)點(diǎn)是不需要額外的內(nèi)存來(lái)區(qū)分原子是否被選中。其次，通過(guò)計(jì)算最大投影值來(lái)尋找最匹配原子的過(guò)程，使用CUBLAS庫(kù)函數(shù)來(lái)實(shí)現(xiàn)。第三，采用矩陣逆更新算法和矩陣向量相乘塊算法來(lái)實(shí)現(xiàn)矩陣求逆和矩陣向量相乘在GPU上的并行計(jì)算。第四，定義兩個(gè)核函數(shù)，一個(gè)用于計(jì)算和的積，另一個(gè)用于更新殘差r。

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

本實(shí)驗(yàn)是在Intel(R)core(TM)i7處理器12 Gbyte內(nèi)存，NVIDIA GTX480顯卡上完成的。實(shí)驗(yàn)的基本設(shè)置如下:傳感矩陣Φ服從標(biāo)準(zhǔn)高斯分布;變換矩陣Ψ不影響算法的運(yùn)行速度，因此令A(yù)=Φ。實(shí)驗(yàn)測(cè)試中，每組參數(shù)設(shè)置下的測(cè)量值都能夠正確地恢復(fù)，選取以下6組實(shí)驗(yàn)數(shù)據(jù)進(jìn)行比較分析，其中n和m分別表示傳感矩陣Α的行和列，而k表示矩陣Α的稀疏度。

4．1 矩陣向量相乘

從表1的實(shí)驗(yàn)結(jié)果可以看出:基于GPU的并行矩陣向量相乘塊算法相比其在CPU的串行算法運(yùn)算速度提高了50多倍;且隨著數(shù)據(jù)量的增加，加速效果更好，因而較好地解決了矩陣向量相乘計(jì)算耗時(shí)長(zhǎng)的難題。

表1 矩陣向量相乘的結(jié)果對(duì)比

4．2 矩陣求逆

從表2的實(shí)驗(yàn)結(jié)果可看出:基于GPU并行矩陣求逆算法相比其在CPU上的串行算法運(yùn)行速度提高了近2倍;其加速比只與稀疏度k相關(guān)，而不依賴(lài)于n和m的大小。即k相同時(shí)，并行求逆算法相比其串行算法加速比相同。

表2 矩陣求逆的結(jié)果對(duì)比

4．3 仿真結(jié)果

從表3的實(shí)驗(yàn)結(jié)果可以看出:基于GPU的并行OMP算法相比其CPU上的串行算法運(yùn)行速度提高了30～44倍:當(dāng)n=16 384時(shí)，加速比可以達(dá)到44倍;n=8 192時(shí)，加速比也可達(dá)到30多倍。實(shí)驗(yàn)結(jié)果充分地證明了基于GPU并行計(jì)算OMP算法的有效性和方案的可行性。圖3是對(duì)表3中6組不同參數(shù)設(shè)置下并行OMP算法的加速效果圖。

表3 仿真實(shí)驗(yàn)的結(jié)果對(duì)比

圖3 并行OMP算法的加速效果

5 小結(jié)

本文在GPU平臺(tái)上設(shè)計(jì)了一個(gè)并行OMP算法。通過(guò)分析OMP算法的復(fù)雜度發(fā)現(xiàn)，該算法的計(jì)算瓶頸在于矩陣求逆和投影部分的矩陣向量相乘。因此在矩陣求逆中使用矩陣逆更新算法從而大幅降低算法的計(jì)算復(fù)雜度，并在GPU平臺(tái)上設(shè)計(jì)了矩陣求逆對(duì)應(yīng)的并行算法;此外，利用分塊思想設(shè)計(jì)了一個(gè)并行的矩陣向量相乘塊算法以加速投影模型。實(shí)驗(yàn)結(jié)果表明，相比于串行算法而言，基于GPU的并行OMP算法可獲得30～44倍的加速比，這充分說(shuō)明了并行算法的有效性和方案的可行性。然而，通過(guò)對(duì)比表1和表3的實(shí)驗(yàn)結(jié)果可以發(fā)現(xiàn)，矩陣向量相乘計(jì)算所耗時(shí)間占總運(yùn)行時(shí)間的2/3以上，也就是說(shuō)，找到更好的矩陣向量相乘算法是提高OMP算法效率的關(guān)鍵。因此，對(duì)OMP算法中矩陣向量相乘的進(jìn)一步優(yōu)化，將是今后工作的重點(diǎn)。

［1］喻玲娟，謝曉春.壓縮感知理論簡(jiǎn)介［J］.電視技術(shù)，2008，32(12):16-18.

［2］高睿.基于壓縮傳感的匹配追蹤重建算法研究［D］.北京:北京交通大學(xué)，2009.

［3］ 方紅，楊海蓉.貪婪算法與壓縮感知理論［J］.自動(dòng)化學(xué)報(bào)，2011，37(12):1413-1421.

［4］方紅，章權(quán)兵，韋穗.改進(jìn)的后退型最優(yōu)正交匹配追蹤圖像重建方法［J］.華南理工大學(xué)學(xué)報(bào):自然科學(xué)版，2008，36(8):23-27.

［5］ 吳凌華，張小川.貪婪算法與壓縮感知理論［J］.電訊技術(shù)，2011，51(1):120-124.

［6］ TROPP J，GILBERT A.Signal recovery from random measurements via orthogonal matching pursuit［J］.IEEE Trans.Inform.Theory，2005，12(53):4655-4666.

［7］劉金碩，曾秋，鄒斌，等.快速魯棒特征算法的CUDA加速優(yōu)化［J］.計(jì)算機(jī)科學(xué)，2014，41(4):24-27.

［8］張舒，褚艷利.GPU高性能運(yùn)算之CUDA［M］.北京:中國(guó)水利水電出版社，2009.

［9］ FARBER R.CUDA application design and development［M］.San Francisco:Morgan Kaufmann，2011.

［10］ JASON S，EDWARD K.CUDA by Example［M］.Boston:Addison-Wesley，2010.

［11］ HAGER W W.Updating the Inverse of a Matrix［J］.SIAM Review，2008，31(2):221-239.

［12］ FUJIMOTO N.Faster matrix-vector multiplication on GeForce 8800GTX［C］//Proc.IEEE International Symposium on Parallel and Distributed Processing.［S.l.］:IEEE Press，2008:14-18.

［13］ FANG Yong，LIANG Chen，WU Jiaji，et al.GPU implementation of orthogonal matching pursuit for compressive sensing［C］//Proc.IEEE 17th International Conference on Parallel and Distributed Systems(ICPADS).［S.l.］:IEEE Press，2008:1044-1047.