熊承義，張　靜，高志榮，雷　夢(mèng)

(1中南民族大學(xué) 電子信息工程學(xué)院, 智能無(wú)線通信湖北省重點(diǎn)實(shí)驗(yàn)室，武漢 430074；2中南民族大學(xué) 計(jì)算機(jī)科學(xué)學(xué)院，武漢 430074)

?

壓縮感知A*OMP重構(gòu)算法的并行化與GPU加速實(shí)現(xiàn)

熊承義1，張靜1，高志榮2，雷夢(mèng)1

(1中南民族大學(xué) 電子信息工程學(xué)院, 智能無(wú)線通信湖北省重點(diǎn)實(shí)驗(yàn)室，武漢 430074；2中南民族大學(xué) 計(jì)算機(jī)科學(xué)學(xué)院，武漢 430074)

摘要針對(duì)壓縮感知系統(tǒng)實(shí)時(shí)應(yīng)用的需要，探討了A*OMP算法的并行設(shè)計(jì)及基于GPU的加速方法.將耗時(shí)長(zhǎng)的矩陣逆運(yùn)算轉(zhuǎn)化為可并行的矩陣/向量操作，并結(jié)合算法本身的關(guān)聯(lián)特性，進(jìn)一步采用迭代法實(shí)現(xiàn)以降低其計(jì)算復(fù)雜度.利用GPU高效的并行運(yùn)算能力，將算法中可并行的矩陣/向量計(jì)算映射到GPU上并行執(zhí)行，在面向Matlab的Jacket軟件平臺(tái)上對(duì)整體串行算法進(jìn)行了并行化的設(shè)計(jì)與實(shí)現(xiàn).在NVIDIA Tesla K20Xm GPU和Intel(R) E5-2650 CPU上進(jìn)行了測(cè)試，實(shí)驗(yàn)結(jié)果表明：對(duì)比CPU平臺(tái)的串行實(shí)現(xiàn)，基于GPU的A*OMP算法整體上可獲得約40倍的加速，實(shí)現(xiàn)了在保持系統(tǒng)較高重構(gòu)質(zhì)量的同時(shí)能有效降低計(jì)算時(shí)間，較好地滿足了系統(tǒng)實(shí)時(shí)性的需要.

關(guān)鍵詞A*OMP算法；并行；加速；圖形處理單元

近幾年來(lái)，壓縮感知[1]作為一個(gè)新的采樣理論，在無(wú)線通信、陣列信號(hào)處理、雷達(dá)成像、生物傳感等領(lǐng)域得到廣泛應(yīng)用.壓縮感知理論主要涉及3個(gè)核心問(wèn)題：稀疏表示，觀測(cè)矩陣的設(shè)計(jì)，以及重建優(yōu)化算法.其中，信號(hào)重構(gòu)算法的選擇是壓縮感知的關(guān)鍵部分.目前，兩類主流的重構(gòu)算法有凸優(yōu)化算法(BP算法、迭代硬閾值法)和貪婪算法(正交匹配追蹤算法、子空間追蹤算法)[2].正交匹配追蹤算法(OMP)[3]作為一種經(jīng)典的重構(gòu)算法，具有較低復(fù)雜度的同時(shí)擁有較好的重構(gòu)質(zhì)量，已被大量用于人臉識(shí)別、語(yǔ)音處理、稀疏表示等領(lǐng)域.隨著人們對(duì)信號(hào)重構(gòu)質(zhì)量要求的提高，越來(lái)越多的OMP改進(jìn)算法被提出.A*OMP[4,5]是伊斯坦布爾大學(xué)的Nazim Burak Karahanoglu和Hakan Erdogan于2011年12月提出的一種OMP改進(jìn)算法，它是一種新的采用最佳優(yōu)先搜索方式的半貪婪算法.該算法本質(zhì)上是將A*搜索算法與OMP算法相結(jié)合的算法，算法中包含了A*搜索：一種最佳優(yōu)先搜索技術(shù)，利用最佳優(yōu)先搜索可以評(píng)估多條路徑.A*OMP算法利用A*搜索技術(shù)，通過(guò)解決OMP算法中局部最優(yōu)的問(wèn)題來(lái)提高重構(gòu)質(zhì)量，但它卻是以較高的時(shí)間復(fù)雜度為代價(jià)，導(dǎo)致壓縮感知系統(tǒng)具有較差的實(shí)時(shí)性.因此，用A*OMP算法對(duì)大規(guī)模的信號(hào)進(jìn)行處理時(shí)，在具有較高重構(gòu)質(zhì)量的同時(shí)提高壓縮感知系統(tǒng)的實(shí)時(shí)性成為亟待解決的問(wèn)題.

AccelerEyes公司于2009年7月發(fā)布的Jacket[6]是面向Matlab和其他高級(jí)語(yǔ)言的一個(gè)軟件平臺(tái)，該平臺(tái)能充分利用GPU強(qiáng)大的浮點(diǎn)性能提供可視計(jì)算加速.Jacket軟件平臺(tái)基于CUDA架構(gòu)，它集成了透明的Matlab用戶界面，用戶可以通過(guò)交互的Matlab桌面和命令窗口使用Matlab編輯器和調(diào)試器方式進(jìn)行并行化設(shè)計(jì).另外，Jacket包含了完全屏蔽底層硬件復(fù)雜性的高層接口以及一套運(yùn)行于Matlab環(huán)境中用于優(yōu)化并行計(jì)算的基礎(chǔ)函數(shù)庫(kù).Jacket GPU還提供了Matlab的CPU數(shù)據(jù)類型，可將Matlab程序移植到GPU運(yùn)行，實(shí)現(xiàn)CUDA[7]的加速.

結(jié)合基于CUDA架構(gòu)的GPU的強(qiáng)勁并行加速能力，針對(duì)壓縮感知系統(tǒng)實(shí)時(shí)應(yīng)用的需要，本文探討了A*OMP算法的并行設(shè)計(jì)及基于GPU的加速方法.將耗時(shí)長(zhǎng)的矩陣逆運(yùn)算轉(zhuǎn)化為可并行的矩陣/向量操作，并結(jié)合算法本身的關(guān)聯(lián)特性，進(jìn)一步采用迭代法實(shí)現(xiàn)以降低其計(jì)算復(fù)雜度.在NVIDIA CUDA計(jì)算框架上將算法中所有可并行的矩陣/向量計(jì)算映射到GPU上并行執(zhí)行，實(shí)現(xiàn)整體算法的并行加速設(shè)計(jì).在Matlab平臺(tái)上，結(jié)合Jacket GPU軟件實(shí)現(xiàn)代碼的編寫和算法的設(shè)計(jì).在NVIDIA Tesla K20Xm GPU和Intel(R) E5-2650 CPU上進(jìn)行了測(cè)試，實(shí)驗(yàn)結(jié)果表明：對(duì)于基于CPU平臺(tái)的串行實(shí)現(xiàn)，基于GPU的A*OMP算法整體上可獲得約40倍的加速.該并行算法在擁有較高的重構(gòu)質(zhì)量的同時(shí)降低了重構(gòu)時(shí)間，有效地提高了重構(gòu)效率.

1壓縮感知的A*OMP重構(gòu)

設(shè)X∈Rn是一維離散信號(hào)，Ψ∈Rn×n為正交變換矩陣，其中R代表實(shí)數(shù)集.假設(shè)信號(hào)X滿足X=Ψ×Z并且Z只包含k<

Y=ΦX=AZ∈Rm,

(1)

其中A=Φ×Ψ為感知矩陣.

如果測(cè)量矩陣Φ服從有限等距特性(RIP)并且與變換矩陣Ψ具有較低的相關(guān)性，那么就可以從低維度的測(cè)量向量Y中有效地重構(gòu)出高維度的稀疏向量Z，從而得到原始信號(hào)X.

A*OMP算法是一種新的采用最佳優(yōu)先搜索方式的半貪婪算法.該算法本質(zhì)上是將A*搜索算法[8]與OMP算法相結(jié)合的算法.A*算法是一種迭代樹形搜索算法，在每次迭代中，當(dāng)前最優(yōu)的路徑和它所有的子集將會(huì)被添加到搜索樹中，當(dāng)最優(yōu)路徑被發(fā)現(xiàn)是完整路徑的時(shí)候搜索中止.

(1)初始化搜索樹：A*搜索將搜索樹的所有可能的路徑長(zhǎng)度初始化為1.由于每次迭代會(huì)在一條已選擇的部分路徑中添加多個(gè)子集，因此，開始搜索時(shí)，限制初始化子集為I(I<

(2)擴(kuò)展已選擇的部分路徑：在典型的A*搜索中，每次迭代中所有最有可能的路徑的子集都會(huì)被添加到搜索樹中，這將產(chǎn)生大量的可能的子集，導(dǎo)致數(shù)據(jù)量急劇增加.為了減少搜索路徑，將每個(gè)節(jié)點(diǎn)的擴(kuò)展個(gè)數(shù)設(shè)為B(即測(cè)量矩陣與已選擇路徑的殘差內(nèi)積絕對(duì)值最大的B個(gè)原子).為了進(jìn)一步減少內(nèi)存需要，限制樹中搜索路徑的最大數(shù)量為P，當(dāng)超過(guò)這個(gè)限制，最壞路徑(最大成本路徑)從樹中刪除.

(3)選擇最佳路徑：最小的殘差誤差是衡量最好路徑的標(biāo)準(zhǔn).但是，考慮到對(duì)長(zhǎng)度不等的多條路徑的評(píng)估，利用關(guān)于殘差的不同假設(shè)提出如下3種不同的代價(jià)函數(shù)：

(2)

(3)

(4)

圖1　A*OMP算法流程圖Fig.1　Flow chart of A*OMP algorithm

2基于GPU的A*OMP并行實(shí)現(xiàn)

對(duì)比OMP重構(gòu)算法，壓縮感知的A*OMP重構(gòu)算法雖然較高地提高了重構(gòu)質(zhì)量，但它卻是以較高的時(shí)間復(fù)雜度為代價(jià)，大大降低了重構(gòu)效率.在進(jìn)行大尺度圖像重構(gòu)時(shí)，A*OMP算法涉及大量耗時(shí)長(zhǎng)的數(shù)據(jù)操作，且主要是易并行的矩陣/向量運(yùn)算和矩陣逆操作，而GPU比較適合處理大量數(shù)據(jù)高度并行集中的問(wèn)題，算法適合采用CUDA架構(gòu)進(jìn)行加速實(shí)現(xiàn).針對(duì)A*OMP算法的數(shù)據(jù)結(jié)構(gòu)特點(diǎn)，本文重點(diǎn)從矩陣/向量操作和矩陣求逆兩個(gè)方面進(jìn)行并行設(shè)計(jì)，最后完成算法整體的并行設(shè)計(jì)方案.

2.1矩陣/向量計(jì)算的并行化

考慮到A*OMP算法實(shí)現(xiàn)重構(gòu)時(shí)涉及大量耗時(shí)長(zhǎng)的矩陣/向量操作(矩陣/向量乘積、矩陣/向量求和、矩陣/向量求差等操作)，而CPU上運(yùn)行的矩陣/向量操作是逐列或逐元素串行進(jìn)行的，在大尺度的矩陣/向量運(yùn)算上，這樣必然會(huì)導(dǎo)致較長(zhǎng)的重構(gòu)時(shí)間.在這些矩陣運(yùn)算中，矩陣的列與列(或行與行、元素與元素)之間是相互獨(dú)立的.針對(duì)這種結(jié)構(gòu)，結(jié)合Jacket的并行設(shè)計(jì)模型，本文用Jacket的GFOR/GEND[9]命令完成并行設(shè)計(jì)方案.標(biāo)準(zhǔn)的FOR循環(huán)是逐步執(zhí)行每個(gè)迭代，而GFOR/GEND是所有的迭代同步執(zhí)行.矩陣/向量并行設(shè)計(jì)的部分程序如圖2所示.

圖2　矩陣/向量的部分并行代碼Fig.2　Partial parallel code for matrix / vector

2.2矩陣逆運(yùn)算的并行設(shè)計(jì)

(5)

(6)

由公式(6)可知，矩陣逆的計(jì)算可以簡(jiǎn)化為多個(gè)矩陣/向量乘、矩陣/向量求和差的運(yùn)算，其中矩陣/向量并行設(shè)計(jì)參考上文的設(shè)計(jì)方案實(shí)現(xiàn).

圖3　A*OMP并行實(shí)現(xiàn)流程圖Fig.3　Flow chart of A*OMP parallel implementation

考慮到算法并行的可行性，利用Jacket函數(shù)庫(kù)，將A*OMP算法中的搜索樹初始化和擴(kuò)展部分路徑兩個(gè)步驟移植到GPU上并行執(zhí)行，而選擇最佳路徑步驟在CPU端實(shí)現(xiàn).重建開始前，首先要將CPU主存中的測(cè)量矩陣和測(cè)量向量數(shù)據(jù)傳輸?shù)紾PU中，即用Jacket的GSINGLE或GDOUBLE函數(shù)修改Matlab命令.然后用Jacket的GZEROS和GONES函數(shù)在GPU上分配用于存放中間結(jié)果的臨時(shí)存儲(chǔ)空間.算法中比較復(fù)雜的數(shù)據(jù)操作，包括求最大值，矩陣轉(zhuǎn)置等，使用Jacket支持函數(shù)(MAX()、CTRANSPOSE())在GPU上執(zhí)行.如果代碼的操作數(shù)據(jù)在GPU的內(nèi)存中，Matlab會(huì)自動(dòng)調(diào)用Jacket庫(kù)內(nèi)置的函數(shù)執(zhí)行相同的矩陣操作.在GPU完成最后一次迭代后，解決方案會(huì)通過(guò)Matlab命令SINGLE和DOUBLE，將數(shù)據(jù)返回到CPU內(nèi)存，這時(shí)所有的臨時(shí)存儲(chǔ)空間都被自動(dòng)釋放.圖3給出了基于GPU的A*OMP壓縮感知算法的設(shè)計(jì)流程圖.

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

實(shí)驗(yàn)運(yùn)行在CPU為Intel(R) E5-2650，內(nèi)存為4GB的臺(tái)式計(jì)算機(jī)上，GPU采用NVIDIA Tesla K20Xm，該顯卡提供2688個(gè)并行處理核心，雙精度浮點(diǎn)峰值性能1294GFlop/s，單精度浮點(diǎn)峰值性能2268GFlop/s，搭配4GB顯卡內(nèi)存.在Jacket 2.1 for Matlab R2011b平臺(tái)上實(shí)現(xiàn)整個(gè)算法的并行設(shè)計(jì).為了驗(yàn)證本文并行化設(shè)計(jì)的優(yōu)勢(shì)，在圖像重構(gòu)質(zhì)量和算法重構(gòu)時(shí)間兩個(gè)方面進(jìn)行了測(cè)試比較.

在重構(gòu)質(zhì)量的比較中，選擇主觀評(píng)價(jià)和客觀評(píng)價(jià)作為質(zhì)量評(píng)估的兩種方法.主觀評(píng)價(jià)時(shí)，選取256×256尺寸大小的lenna圖像進(jìn)行壓縮感知重構(gòu)，如圖4所示，在相同條件下，對(duì)3種不同的重構(gòu)算法進(jìn)行試驗(yàn)，即A*OMP算法、OMP算法和BP算法.其中，OMP重構(gòu)的PSNR為24.15dB，BP重構(gòu)的PSNR為27.06dB,而A*OMP重構(gòu)的PSNR為31.05dB.客觀評(píng)價(jià)時(shí)，選取6種不同尺寸大小的lenna標(biāo)準(zhǔn)測(cè)試圖像進(jìn)行PSNR值的測(cè)試比較，質(zhì)量結(jié)果如表1所示.其中采樣率設(shè)為0.5，(稀疏向量中非零元素的個(gè)數(shù))取測(cè)量向量長(zhǎng)度的1/8.對(duì)比前兩種重構(gòu)算法，A*OMP表現(xiàn)出更好的重構(gòu)質(zhì)量.

圖4　三種算法重構(gòu)圖像的對(duì)比Fig.4　Comparison of reconstruction quality of three algorithms

在進(jìn)行重構(gòu)時(shí)間的比較中，選取不同尺寸大小(16×16，32×32，64×64，128×128，256×256，512×512)的lenna圖像進(jìn)行試驗(yàn)測(cè)試.對(duì)于256×256和512×512尺寸的兩幅圖像，由于本實(shí)驗(yàn)硬件GPU

表1　三種算法的重構(gòu)圖像的PSNR值比較

內(nèi)存的限制，結(jié)合分塊的思想，將這兩幅圖像均分成多個(gè)64×64的圖像塊來(lái)完成并行設(shè)計(jì).該方法雖然可以較好地降低重構(gòu)時(shí)間，但重構(gòu)的圖像會(huì)產(chǎn)生塊效應(yīng)，導(dǎo)致重構(gòu)圖像具有較差的視覺(jué)效果.特別地，由于對(duì)256×256和512×512尺寸的兩幅圖像進(jìn)行了64×64塊的處理，相當(dāng)于將多個(gè)64×64的圖像進(jìn)行了串行循環(huán)，所以加速效果和單個(gè)64×64尺寸的圖像很接近.此外，為了避免Matlab和Jacket啟動(dòng)時(shí)間影響實(shí)驗(yàn)結(jié)果，所有圖像的重建實(shí)驗(yàn)計(jì)時(shí)11次，丟棄第一次運(yùn)行時(shí)間，選擇剩余10次的平均值作為運(yùn)行時(shí)間.表2給出了不同尺寸圖像在CPU和GPU兩種平臺(tái)下的執(zhí)行時(shí)間和重構(gòu)圖像的PSNR值的比較結(jié)果.其中采樣率設(shè)為0.25，搜索長(zhǎng)度(稀疏向量中非零元素的個(gè)數(shù))取測(cè)量向量長(zhǎng)度的1/8.試驗(yàn)結(jié)果表明：當(dāng)圖像尺寸為16×16大小時(shí)，GPU重建時(shí)間比CPU重建時(shí)間要長(zhǎng).這是因?yàn)閷?shù)據(jù)從CPU內(nèi)存?zhèn)鬏數(shù)紾PU端時(shí)存在較高的通信耗時(shí)，對(duì)于較小的圖像矩陣，這種開銷會(huì)占據(jù)運(yùn)行時(shí)間的大部分，大大降低了GPU的計(jì)算性能.而當(dāng)圖像為128×128大小時(shí)，算法有著較高的加速比，重構(gòu)算法在CPU上的重建需要56032s，而在GPU上僅需要1407s，加速比高達(dá)40倍.可以得出，圖像尺寸是影響加速效果比較直接的因素，在硬件資源充足的條件下，測(cè)試圖像的尺寸越大，加速效果會(huì)越明顯；基于CPU的串行算法與基于GPU的并行加速算法在圖像重建質(zhì)量方面保持一致.

表2　基于CPU的串行實(shí)現(xiàn)和基于GPU的并行實(shí)現(xiàn)的比較

4結(jié)語(yǔ)

基于Jacket技術(shù)，本文在Matlab平臺(tái)上對(duì)A*OMP算法進(jìn)行了并行化設(shè)計(jì)與實(shí)現(xiàn)，由于GPU和CPU內(nèi)存交互通信的損耗，除了像素較小的圖像，與CPU串行實(shí)現(xiàn)相比，算法的GPU并行實(shí)現(xiàn)有明顯的性能優(yōu)勢(shì).實(shí)驗(yàn)結(jié)果表明，圖像尺寸是影響加速效果比較直接的因素，在硬件資源充足的條件下，測(cè)試圖像的尺寸越大，加速效果會(huì)越明顯.本文實(shí)現(xiàn)算法不僅具有較高的重構(gòu)質(zhì)量而且較大程度地降低了重構(gòu)時(shí)間，有效地提高了重構(gòu)效率.

參考文獻(xiàn)

[1]Donoho D L. Compressed sensing [J]. IEEE Transactions on Information Theory,2006,52(4):1289-1306.

[2]戴瓊海,付長(zhǎng)軍,季向陽(yáng).壓縮感知研究[J].計(jì)算機(jī)學(xué)報(bào),2011,34(3):425-434.

[3]Usman M , Prieto C , Odille F ,et al. A computationally efficient OMP-based compressed sensing reconstruction for dynamic MRI [J]. Physics in Medicine & Biology, 2011, 56(7):99-114.

[4]Karahanoglu N B, Erdogan H. A*Orthogonal Matching Pursuit: Best-First Search for Compressed Sensing Signal Recovery[J]. Digital Signal Processing, 2010, 22(4):555-568.

[5]Karahanoglu N B, Erdogan H. Compressed sensing signal recovery via A*Orthogonal Matching Pursuit[C]// IEEE.Acoustics, Speech and Signal Processing (ICASSP), 2011 IEEE International Conference on. New Jersey: IEEE Press, 2011:3732-3735.

[6]張百達(dá),唐玉華,吳俊杰,等. Speeding up the MATLAB complex networks package using graphic processors [J]. Chinese Physics B, 2011, 20(9):460-467.

[7]Rob Farber. CUDA Application Design and Development [M]. Burlington：Morgan Kaufmann，2011.

[8]Dechter R, Pearl J. The Optimality of A*[J]. Symbolic Computation, 1988, 70(5):1841-1863.

[9]Kuldeep Yadav, Ankush Mittal, M A Ansar,et al. Parallel Implementation of Compressed Sensing Algorithm on CUDA- GPU [J].International Journal of Computer Science and Information Security, 2011, 9(3):112.

[10]FANG Y, CHEN L,Wu J,et a1.GPU implementation of orthogonal matching pursuit for compressive sensing[C]//IEEE. Proc ICPADS 2011. New Jersey: IEEE Press, 2011:1044-1047.

[11]Chun-yuan Deng，Yi-min Wei. Some New Results of the Sherman-Morrison-Woodbury Formula[C]//ICMP. Proceeding of The Sixth International Conference of Matrices and Operators.Chengdu:ICMP, 2011:220-223.

Parallelization and GPU Based Realization of A*OMP Algorithm for Compressed Sensing

XiongChengyi1,ZhangJing1,GaoZhirong2,LeiMeng1

(1 College of Electronic and Information Engineering, Hubei Key Lab of Intelligent Wireless Communication,South-Central University for Nationalities, Wuhan 430074, China；2 College of Computer Science, South-Central University for Nationalities, Wuhan 430074, China)

AbstractAiming at the need of real-time application of compressed sensing system, this paper discusses the design of parallel A*OMP algorithm and the acceleration method based on Graphic Processing Unit (GPU). It adopts the idea of iterative method to avoid complex matrix operations in the matrix inversion part，which combine with the correlation properties of the algorithm itself，and the time-consuming matrix inverse operation is transformed into parallel matrix / vector operations. Using GPU′s efficient parallel computing ability，all parallel matrix / vector computations in the algorithm are mapped onto the GPU parallel execution，and the parallel design and implementation of the whole serial algorithm are finished on the Jacket software platform oriented MATLAB. The test on the NVIDIA Tesla K20Xm GPU and Intel (R) E5-2650 CPU was conducted. The experimental results show that A*OMP implementation based on GPU can get 40 times acceleration when comparing to serial implementation in CPU platform. The proposed implementation method could efficiently reduce computation time with keeping good reconstruction quality simultaneously, which meets the real-time requirement of system well.

KeywordsA*OMP algorithm；parallel；acceleration；GPU

收稿日期2015-11-27

作者簡(jiǎn)介熊承義(1969-)，男，教授，碩士生導(dǎo)師，研究方向：圖像處理與模式識(shí)別，E-mail：xiongcy@mail.scuec.edu.cn

基金項(xiàng)目國(guó)家自然科學(xué)基金資助項(xiàng)目(61471400,61201268)；湖北省自然科學(xué)基金資助項(xiàng)目(2013CFC118)，中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)(CZW14018)

中圖分類號(hào)TP391.41

文獻(xiàn)標(biāo)識(shí)碼A

文章編號(hào)1672-4321(2016)02-0079-06