康 凱，王粒賓，鐘子發(fā)

(1.電子工程學(xué)院 通信對抗系,安徽 合肥 230037；2.第二炮兵工程大學(xué)信息工程系,陜西 西安 710025；3.解放軍61922部隊(duì),北京 100094)

信號稀疏重構(gòu)是壓縮感知(compressive sensing,簡稱CS)[1]的核心問題之一.假設(shè)稀疏線性模型為b=Ax，其中，b為線性測量值，A為欠定的感知矩陣，x為待重構(gòu)信號.在標(biāo)準(zhǔn)的CS框架下，重構(gòu)x時(shí)僅僅限制x中非零元素的個(gè)數(shù)(即稀疏度)不能超過某一閾值k，并沒有對非零元素的結(jié)構(gòu)加以約束，也就是說，非零元素可以出現(xiàn)在x中的任何位置.但是在某些應(yīng)用中，如多帶信號重構(gòu)[2]及子空間信號重構(gòu)[3]等，x的非零元素呈聚集分布，此時(shí)稱x滿足塊稀疏(block sparsity)模型.顯然，塊稀疏的x可以在標(biāo)準(zhǔn)CS框架下進(jìn)行重構(gòu)，但是將這種塊稀疏結(jié)構(gòu)加以利用，可以獲得更好的重構(gòu)性能，因此，塊稀疏信號重構(gòu)已成為最近的研究熱點(diǎn).

塊稀疏信號重構(gòu)算法可分為2類：一類為貪婪追蹤類算法，如塊正交匹配追蹤[4-5](block orthogonal matching pursuit,簡稱Block OMP)和塊壓縮采樣匹配追蹤[6](block compressive sampling matching pursuit,簡稱Block CoSaMP)；另一類為凸松弛類算法，如基于二階錐規(guī)劃(second-order cone programming,簡稱SOCP)和內(nèi)點(diǎn)法(interior point,簡稱IP)的SOCP+IP算法[3]，基于半定規(guī)劃(semi-definite programming,簡稱SDP)和內(nèi)點(diǎn)法的SDP+IP算法[7].貪婪追蹤類算法的特點(diǎn)是計(jì)算量小，但重構(gòu)誤差大，且其計(jì)算量隨著塊稀疏度的增大而急劇增大；凸松弛類算法的重構(gòu)精度高，但是現(xiàn)有基于SOCP和SDP的算法計(jì)量大，特別在大尺度問題下，計(jì)算量將是無法克服的困難.

應(yīng)用凸松弛模型，作者提出一種基于交替方向法(alternating direction method,簡稱ADM)[8-9]的快速塊稀疏重構(gòu)算法(block-sparse recovery based on ADM,簡稱ADM-BSR)，通過仿真說明該算法的有效性.

1 問題說明

設(shè)線性測量值b∈RM，感知矩陣A∈RM×N，塊稀疏信號x∈RN，加性噪聲n∈RM，塊稀疏模型如圖1所示[7].

圖1 塊稀疏模型示意圖

圖1中，d為塊長度，m為塊數(shù)，滿足N=dm.x[i]=(xd(i-1)+1,xd(i-1)+2,…,xdi)T表示第i個(gè)塊，A[i]=(ad(i-1)+1,ad(i-1)+2,…,adi)表示對應(yīng)于x[i]的子矩陣.塊稀疏的數(shù)學(xué)模型為

(1)

塊稀疏重構(gòu)的數(shù)學(xué)模型為[3,7]

(2)

其中:ε為噪聲n的均方差.直接求解上式是NP-難的[7]，一類近似求解方法是對式(2)進(jìn)行凸松弛，轉(zhuǎn)化為凸優(yōu)化問題[3,7]，即

(3)

式(3)稱為塊稀疏重構(gòu)的l2/l1模型.l2/l1模型進(jìn)行塊稀疏重構(gòu)的優(yōu)勢在于高精度，但是現(xiàn)有的SOCP+IP算法和SDP+IP算法計(jì)算量大，例如SOCP+IP的計(jì)算復(fù)雜度高達(dá)O(N3)[10]，由此可見在大尺度問題下上述2種方法不再適用.

針對此，該文目標(biāo)是基于l2/l1模型提出一種快速塊稀疏重構(gòu)算法，在保持較高重構(gòu)精度的前提下，提升計(jì)算效率.在此需要說明的是，在統(tǒng)計(jì)學(xué)領(lǐng)域Yuan等利用變量的組聚集特性，提出一種稱之為group lasso[11]的問題模型.在group lasso中，目標(biāo)函數(shù)為

(4)

其中:τ為正則化參數(shù).

式(3)、(4)二者的不同點(diǎn)在于：式(3)中的參數(shù)ε具有明確的物理意義，其值大小與噪聲的均方差成比例，設(shè)置容易；而式(4)中的τ值則很難選擇.因此，在信號處理領(lǐng)域研究式(3)的求解更具有實(shí)際意義.

2 ADM框架

設(shè)f(x):RM→R和g(y):RN→R均為凸函數(shù)，B∈RP×M，C∈RP×N，w∈RP.考慮下述凸優(yōu)化問題

(5)

其中:變量x和y在目標(biāo)函數(shù)中是可分離且可單獨(dú)求解的，而在等式約束下則是相互關(guān)聯(lián)的.

(5)式的增廣拉格朗日函數(shù)為

(6)

其中:λ∈RP為乘子，β為懲罰因子.

利用ADM得到(6)式的迭代形式為

(7)

已經(jīng)證明，在每次迭代中，若xk+1和yk+1能夠精確求解，則ADM是收斂于式(5)的全局最優(yōu)解[8].

3 基于ADM的塊稀疏信號重構(gòu)算法

在塊稀疏信號重構(gòu)l2/l1模型的基礎(chǔ)上，利用ADM提出一種快速算法，稱之為ADM-BSR算法，下面詳細(xì)介紹ADM-BSR算法的推導(dǎo)過程.對式(3)進(jìn)行變量分裂，可得

(8)

式(8)的增廣拉格朗日函數(shù)為

(9)

利用ADM求解式(9)的迭代形式為

(10)

其中:PBε為到集合Bε?{a∈RM：‖a‖2≤ε}上的正交投影算子.

式(10)的4個(gè)子問題中，nk+1、λk+1和βk+1均可得到閉式解，而xk+1只能通過迭代進(jìn)行求解.考慮到x[i](i=1,…,n)是x中互不重疊的子塊，因此選擇BCD法[12]求解xk+1.為方便描述，定義dk=nk+1-b-λk/βk，目標(biāo)函數(shù)轉(zhuǎn)化為

(11)

可見F(x)是一個(gè)非光滑的凸函數(shù).分別針對x[i](i=1,…,m)求次梯度，可得[13]

(12)

其中

定義

(13)

因此，由式(12)可得求解x[i]的KKT條件為

(14)

當(dāng)x[i]≠0時(shí)，定義

(15)

對于隨機(jī)測量矩陣，對其列進(jìn)行標(biāo)準(zhǔn)正交化處理后，可假設(shè)AT[i]A[i]=I.將式(14)代入式(12)，則有

(16)

由式(15)、(16)可知，si與x[i]共線，因此有

(17)

4 ADM-BSR算法分析

4.1 計(jì)算復(fù)雜度

對于某一具體問題，無法估計(jì)ADM-BSR算法所需要的迭代次數(shù).但是對于每次迭代， ADM-BSR算法僅需要向量與標(biāo)量乘積運(yùn)算、向量之間相加運(yùn)算、子矩陣A[l]和AT[l]與向量的乘積運(yùn)算、矩陣A與向量的乘積運(yùn)算，不存在矩陣與矩陣乘積運(yùn)算以及矩陣求逆運(yùn)算，因此，相比于SOCP+IP算法和SDP+IP算法，ADM-BSR算法的計(jì)算量大大減少.

4.2 收斂性

首先，分析采用BCD方法求解xk+1的收斂性.Tseng在文獻(xiàn)[12]中指出，若目標(biāo)函數(shù)由一個(gè)光滑的凸函數(shù)和一個(gè)非光滑的塊可分凸函數(shù)組成，在Gauss-Seidel循環(huán)規(guī)則下BCD方法可以保證得到全局最優(yōu)解.顯然，解xk+1的目標(biāo)函數(shù)F(x)滿足這個(gè)要求，并且在ADM-BSR中也可采取Gauss-Seidel循環(huán)規(guī)則，因此從理論上可以保證隨著迭代次數(shù)的增加，xk+1可無限接近于F(x)的全局最優(yōu)解.但是，在算法的實(shí)際實(shí)施過程中迭代次數(shù)往往有限，因此xk+1將不再是F(x)的最優(yōu)解，而是次優(yōu)解.這種情況下，子問題不能精確求解時(shí)，ADM的收斂性將難以證明.

然而，ADM-BSR算法利用BCD方法迭代求解xk+1時(shí)，以每個(gè)子塊x[i]的KKT條件作為終止條件，可以保證xk+1非常接近于F(x)的最優(yōu)解，且在迭代過程中不會出現(xiàn)發(fā)散現(xiàn)象，因此可以預(yù)測ADM-BSR算法是收斂的，大量仿真實(shí)驗(yàn)也證明了這點(diǎn).

5 仿真實(shí)驗(yàn)

與文獻(xiàn)[3-7]相同，進(jìn)行綜合實(shí)驗(yàn).5.1節(jié)的仿真用于驗(yàn)證利用信號塊稀疏特性可以提高重構(gòu)精度；5.2節(jié)通過比較ADM-BSR算法、Block OMP算法[4]和Block CoSaMP算法[6]的仿真實(shí)驗(yàn)，說明ADM-BSR在計(jì)算效率上的優(yōu)勢.

5.1 塊稀疏重構(gòu)和標(biāo)準(zhǔn)稀疏重構(gòu)的精度比較

在綜合實(shí)驗(yàn)中，采取如下的仿真場景：感知矩陣A為一個(gè)M×N維隨機(jī)矩陣，其元素獨(dú)立且服從均值為0、方差為1/(2N)的高斯分布，其列進(jìn)行標(biāo)準(zhǔn)正交化處理.線性測量值b=Axtrue+n，其中n為均值為0、方差為σ2的高斯白噪聲.塊長度設(shè)為d，則xtrue共有m=N/d個(gè)塊，xtrue的K個(gè)非零塊在[1,m]中隨機(jī)選取，非零塊中元素值服從標(biāo)準(zhǔn)高斯分布.與文獻(xiàn)[14]相同，定義重構(gòu)誤差為均方誤差(mean-squared error,簡稱MSE)，即

(18)

問題參數(shù)設(shè)置如下：測量維數(shù)M=2 048，信號維數(shù)N=4 096，噪聲均方差σ=0.1，塊稀疏度K=16，塊長度d=64，標(biāo)準(zhǔn)稀疏度k=Kd=1 024.

圖2 標(biāo)準(zhǔn)稀疏重構(gòu)與塊稀疏重構(gòu)性能比較

從仿真結(jié)果可知，標(biāo)準(zhǔn)l1模型下重構(gòu)誤差為MSEl1=2.55×10-2，而l2/l1模型下重構(gòu)誤差僅為MSEl2/l1=4.028×10-4.因此，利用信號的塊稀疏特性，可以大大提高重構(gòu)精度.

5.2 ADM-BSR與Block OMP、Block CoSaMP計(jì)算效率的比較

仿真場景設(shè)置與5.1節(jié)相同.問題參數(shù)為M=2 048，N=4 096，σ=0.1，d=64，在塊稀疏度K從6變化到30過程中計(jì)算3種算法的計(jì)算耗時(shí)和均方誤差.設(shè)置Block CoSaMP和ADM-BSR的迭代次數(shù)為20，Block OMP的迭代次數(shù)與稀疏度K相同.仿真環(huán)境為Matlab R2008a，P4 2.6 G(雙核)，1 G內(nèi)存，Windows XP操作系統(tǒng).

圖3為3種算法的計(jì)算耗時(shí)隨塊稀疏度變化的曲線，圖4為3種算法的均方誤差隨塊稀疏度變化的曲線.仿真曲線為50次平均的結(jié)果.

圖3 3種算法的計(jì)算耗時(shí)隨塊稀疏度變化的曲線

圖4 3種算法的均方誤差隨塊稀疏度變化的曲線

從圖3、4可以看出：1)在計(jì)算耗時(shí)方面，ADM-BSR算法的計(jì)算耗時(shí)隨著塊稀疏度增加而基本保持不變，而Block CoSaMP和Block OMP則隨著塊稀疏度的增加急劇增加，自K=12(塊稀疏比δ=K/m=0.18)起，ADM-BSR算法的計(jì)算耗時(shí)將小于貪婪追蹤類算法.這直觀說明了在高塊稀疏度下，ADM-BSR的計(jì)算優(yōu)勢十分明顯；2)在重構(gòu)誤差方面，總的來看，Block OMP的MSE最大，ADM-BSR和Block CoSaMP的相對較小.對于ADM-BSR和Block CoSaMP，在低塊稀疏度時(shí)，兩者的MSE非常接近，但隨著塊稀疏度增加，ADM-BSR的MSE小于Block CoSaMP的.

6 結(jié)束語

作者研究了塊稀疏信號的重構(gòu)問題.基于塊稀疏重構(gòu)的l2/l1模型，利用ADM提出一種ADM-BSR算法，分析了算法的計(jì)算復(fù)雜度和收斂性.通過仿真實(shí)驗(yàn)可以看出ADM-BSR算法的重構(gòu)精度高于基于貪婪追蹤的Block CoSaMP和Block OMP，且在塊稀疏度較高時(shí)，計(jì)算耗時(shí)也遠(yuǎn)小于這2種算法.

參考文獻(xiàn):

[1] Candes E, Wakin B M. An introduction to compressive sampling[J].IEEE Signal Processing Magazine，2008,25(2):21-30.

[2] Mishali M, Eldar C Y. Blind multi-band signal reconstruction: compressed sensing for analog signals[J].IEEE Transactions on Signal Processing,2009,57(3):993-1009.

[3] Eldar C Y, Mishali M. Robust recovery of signals from a structured union of subspaces[J].IEEE Transactions on Information Theory,2009,55(11):5302-5316.

[4] Eldar C Y, Kuppinger P, Bolcsker H. Block-sparse signals: uncertainty relations and efficient recovery[J].IEEE Transactions Signal Processing,2010,57(6):3042-3054.

[5] Lv X, Wan C, Bi G. Block orthogonal greedy algorithm for stable recovery of block-sparse signal representations[J].Signal Processing,2010,90(12):3265-3277.

[6] Baraniuk G R, Cevher V, Duarte F M, et al. Model-based compressive sensing[J].IEEE Transactions on Information Theory,2010,56(4):1982-2001.

[7] Stojnic M, Paryaresh F, Hassibi B. On the reconstruction of block-sparse signals with an optimal number of measurements[J].IEEE Transactions Signal Processing,2009,57(8):3075-3085.

[8] Yang J, Zhang Y. Alternating direction algorithms for l1-problems in compressive sensing[R].Houston: Rice University,2009.

[9] Afonso V M, Bioucas-Dias M J, Figueiredo A M. Fast image reconvery using variable splitting and constrained optimization[J].IEEE Transactions on Image Processing,2010,19(9):2345-2356.

[10] Rakotomamonjy A. Surverying and comparing simultaneous sparse approximation (or group lasso)algorithm[J].Signal Processing,2011,91:1505-1526.

[11] Yuan M, Lin Y. Model selection and estimation in regression with grouped variables[J].Journal of Royal Statistical Society, Series B,2006,68(1):49-67.

[12] Tseng P. Convergence of block coordinate descent method for non-differentiable minimization[J].Journal of Optimization Theory and Application,2001,109:475-494.

[13] Bertsekas D, Nedic A, Ozdaglar A. Convex analysis and optimization[M].Nashua: Athena Scientific,2003.

[14] Wright J S, Nowak D R, Figueiredo A M. Sparse reconstruction by separable approximation[J].IEEE Transactions on Image Processing,2009,57(7):2479-2493.