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

熊承義，張　靜，高志榮，雷　夢

(1中南民族大學(xué) 電子信息工程學(xué)院, 智能無線通信湖北省重點實驗室，武漢 430074；2中南民族大學(xué) 計算機科學(xué)學(xué)院，武漢 430074)

?

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

熊承義1，張靜1，高志榮2，雷夢1

(1中南民族大學(xué) 電子信息工程學(xué)院, 智能無線通信湖北省重點實驗室，武漢 430074；2中南民族大學(xué) 計算機科學(xué)學(xué)院，武漢 430074)

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

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

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

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

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

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

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

Y=ΦX=AZ∈Rm,

(1)

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

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

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

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

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

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

(2)

(3)

(4)

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

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

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

2.1矩陣/向量計算的并行化

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

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

2.2矩陣逆運算的并行設(shè)計

(5)

(6)

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

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

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

3實驗結(jié)果及分析

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

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

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

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

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

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

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

4結(jié)語

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

參考文獻

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

[2]戴瓊海,付長軍,季向陽.壓縮感知研究[J].計算機學(xué)報,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]張百達,唐玉華,吳俊杰,等. 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

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

基金項目國家自然科學(xué)基金資助項目(61471400,61201268)；湖北省自然科學(xué)基金資助項目(2013CFC118)，中央高?；究蒲袠I(yè)務(wù)費專項(CZW14018)

中圖分類號TP391.41

文獻標(biāo)識碼A

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