康 凱,王粒賓,鐘子發(fā)
(1.電子工程學(xué)院 通信對(duì)抗系,安徽 合肥 230037;2.第二炮兵工程大學(xué)信息工程系,陜西 西安 710025;3.解放軍61922部隊(duì),北京 100094)
信號(hào)稀疏重構(gòu)是壓縮感知(compressive sensing,簡(jiǎn)稱(chēng)CS)[1]的核心問(wèn)題之一.假設(shè)稀疏線性模型為b=Ax,其中,b為線性測(cè)量值,A為欠定的感知矩陣,x為待重構(gòu)信號(hào).在標(biāo)準(zhǔn)的CS框架下,重構(gòu)x時(shí)僅僅限制x中非零元素的個(gè)數(shù)(即稀疏度)不能超過(guò)某一閾值k,并沒(méi)有對(duì)非零元素的結(jié)構(gòu)加以約束,也就是說(shuō),非零元素可以出現(xiàn)在x中的任何位置.但是在某些應(yīng)用中,如多帶信號(hào)重構(gòu)[2]及子空間信號(hào)重構(gòu)[3]等,x的非零元素呈聚集分布,此時(shí)稱(chēng)x滿(mǎn)足塊稀疏(block sparsity)模型.顯然,塊稀疏的x可以在標(biāo)準(zhǔn)CS框架下進(jìn)行重構(gòu),但是將這種塊稀疏結(jié)構(gòu)加以利用,可以獲得更好的重構(gòu)性能,因此,塊稀疏信號(hào)重構(gòu)已成為最近的研究熱點(diǎn).
塊稀疏信號(hào)重構(gòu)算法可分為2類(lèi):一類(lèi)為貪婪追蹤類(lèi)算法,如塊正交匹配追蹤[4-5](block orthogonal matching pursuit,簡(jiǎn)稱(chēng)Block OMP)和塊壓縮采樣匹配追蹤[6](block compressive sampling matching pursuit,簡(jiǎn)稱(chēng)Block CoSaMP);另一類(lèi)為凸松弛類(lèi)算法,如基于二階錐規(guī)劃(second-order cone programming,簡(jiǎn)稱(chēng)SOCP)和內(nèi)點(diǎn)法(interior point,簡(jiǎn)稱(chēng)IP)的SOCP+IP算法[3],基于半定規(guī)劃(semi-definite programming,簡(jiǎn)稱(chēng)SDP)和內(nèi)點(diǎn)法的SDP+IP算法[7].貪婪追蹤類(lèi)算法的特點(diǎn)是計(jì)算量小,但重構(gòu)誤差大,且其計(jì)算量隨著塊稀疏度的增大而急劇增大;凸松弛類(lèi)算法的重構(gòu)精度高,但是現(xiàn)有基于SOCP和SDP的算法計(jì)量大,特別在大尺度問(wèn)題下,計(jì)算量將是無(wú)法克服的困難.
應(yīng)用凸松弛模型,作者提出一種基于交替方向法(alternating direction method,簡(jiǎn)稱(chēng)ADM)[8-9]的快速塊稀疏重構(gòu)算法(block-sparse recovery based on ADM,簡(jiǎn)稱(chēng)ADM-BSR),通過(guò)仿真說(shuō)明該算法的有效性.
設(shè)線性測(cè)量值b∈RM,感知矩陣A∈RM×N,塊稀疏信號(hào)x∈RN,加性噪聲n∈RM,塊稀疏模型如圖1所示[7].
圖1 塊稀疏模型示意圖
圖1中,d為塊長(zhǎng)度,m為塊數(shù),滿(mǎn)足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)表示對(duì)應(yīng)于x[i]的子矩陣.塊稀疏的數(shù)學(xué)模型為
(1)
塊稀疏重構(gòu)的數(shù)學(xué)模型為[3,7]
(2)
其中:ε為噪聲n的均方差.直接求解上式是NP-難的[7],一類(lèi)近似求解方法是對(duì)式(2)進(jìn)行凸松弛,轉(zhuǎn)化為凸優(yōu)化問(wèn)題[3,7],即
(3)
式(3)稱(chēng)為塊稀疏重構(gòu)的l2/l1模型.l2/l1模型進(jìn)行塊稀疏重構(gòu)的優(yōu)勢(shì)在于高精度,但是現(xiàn)有的SOCP+IP算法和SDP+IP算法計(jì)算量大,例如SOCP+IP的計(jì)算復(fù)雜度高達(dá)O(N3)[10],由此可見(jiàn)在大尺度問(wèn)題下上述2種方法不再適用.
針對(duì)此,該文目標(biāo)是基于l2/l1模型提出一種快速塊稀疏重構(gòu)算法,在保持較高重構(gòu)精度的前提下,提升計(jì)算效率.在此需要說(shuō)明的是,在統(tǒng)計(jì)學(xué)領(lǐng)域Yuan等利用變量的組聚集特性,提出一種稱(chēng)之為group lasso[11]的問(wèn)題模型.在group lasso中,目標(biāo)函數(shù)為
(4)
其中:τ為正則化參數(shù).
式(3)、(4)二者的不同點(diǎn)在于:式(3)中的參數(shù)ε具有明確的物理意義,其值大小與噪聲的均方差成比例,設(shè)置容易;而式(4)中的τ值則很難選擇.因此,在信號(hào)處理領(lǐng)域研究式(3)的求解更具有實(shí)際意義.
設(shè)f(x):RM→R和g(y):RN→R均為凸函數(shù),B∈RP×M,C∈RP×N,w∈RP.考慮下述凸優(yōu)化問(wèn)題
(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].
在塊稀疏信號(hào)重構(gòu)l2/l1模型的基礎(chǔ)上,利用ADM提出一種快速算法,稱(chēng)之為ADM-BSR算法,下面詳細(xì)介紹ADM-BSR算法的推導(dǎo)過(guò)程.對(duì)式(3)進(jìn)行變量分裂,可得
(8)
式(8)的增廣拉格朗日函數(shù)為
(9)
利用ADM求解式(9)的迭代形式為
(10)
其中:PBε為到集合Bε?{a∈RM:‖a‖2≤ε}上的正交投影算子.
式(10)的4個(gè)子問(wèn)題中,nk+1、λk+1和βk+1均可得到閉式解,而xk+1只能通過(guò)迭代進(jìn)行求解.考慮到x[i](i=1,…,n)是x中互不重疊的子塊,因此選擇BCD法[12]求解xk+1.為方便描述,定義dk=nk+1-b-λk/βk,目標(biāo)函數(shù)轉(zhuǎn)化為
(11)
可見(jiàn)F(x)是一個(gè)非光滑的凸函數(shù).分別針對(duì)x[i](i=1,…,m)求次梯度,可得[13]
(12)
其中
定義
(13)
因此,由式(12)可得求解x[i]的KKT條件為
(14)
當(dāng)x[i]≠0時(shí),定義
(15)
對(duì)于隨機(jī)測(cè)量矩陣,對(duì)其列進(jìn)行標(biāo)準(zhǔn)正交化處理后,可假設(shè)AT[i]A[i]=I.將式(14)代入式(12),則有
(16)
由式(15)、(16)可知,si與x[i]共線,因此有
(17)
對(duì)于某一具體問(wèn)題,無(wú)法估計(jì)ADM-BSR算法所需要的迭代次數(shù).但是對(duì)于每次迭代, 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ì)算量大大減少.
首先,分析采用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)滿(mǎn)足這個(gè)要求,并且在ADM-BSR中也可采取Gauss-Seidel循環(huán)規(guī)則,因此從理論上可以保證隨著迭代次數(shù)的增加,xk+1可無(wú)限接近于F(x)的全局最優(yōu)解.但是,在算法的實(shí)際實(shí)施過(guò)程中迭代次數(shù)往往有限,因此xk+1將不再是F(x)的最優(yōu)解,而是次優(yōu)解.這種情況下,子問(wèn)題不能精確求解時(shí),ADM的收斂性將難以證明.
然而,ADM-BSR算法利用BCD方法迭代求解xk+1時(shí),以每個(gè)子塊x[i]的KKT條件作為終止條件,可以保證xk+1非常接近于F(x)的最優(yōu)解,且在迭代過(guò)程中不會(huì)出現(xiàn)發(fā)散現(xiàn)象,因此可以預(yù)測(cè)ADM-BSR算法是收斂的,大量仿真實(shí)驗(yàn)也證明了這點(diǎn).
與文獻(xiàn)[3-7]相同,進(jìn)行綜合實(shí)驗(yàn).5.1節(jié)的仿真用于驗(yàn)證利用信號(hào)塊稀疏特性可以提高重構(gòu)精度;5.2節(jié)通過(guò)比較ADM-BSR算法、Block OMP算法[4]和Block CoSaMP算法[6]的仿真實(shí)驗(yàn),說(shuō)明ADM-BSR在計(jì)算效率上的優(yōu)勢(shì).
在綜合實(shí)驗(yàn)中,采取如下的仿真場(chǎng)景:感知矩陣A為一個(gè)M×N維隨機(jī)矩陣,其元素獨(dú)立且服從均值為0、方差為1/(2N)的高斯分布,其列進(jìn)行標(biāo)準(zhǔn)正交化處理.線性測(cè)量值b=Axtrue+n,其中n為均值為0、方差為σ2的高斯白噪聲.塊長(zhǎng)度設(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,簡(jiǎn)稱(chēng)MSE),即
(18)
問(wèn)題參數(shù)設(shè)置如下:測(cè)量維數(shù)M=2 048,信號(hào)維數(shù)N=4 096,噪聲均方差σ=0.1,塊稀疏度K=16,塊長(zhǎng)度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.因此,利用信號(hào)的塊稀疏特性,可以大大提高重構(gòu)精度.
仿真場(chǎng)景設(shè)置與5.1節(jié)相同.問(wèn)題參數(shù)為M=2 048,N=4 096,σ=0.1,d=64,在塊稀疏度K從6變化到30過(guò)程中計(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í)隨著塊稀疏度增加而基本保持不變,而B(niǎo)lock CoSaMP和Block OMP則隨著塊稀疏度的增加急劇增加,自K=12(塊稀疏比δ=K/m=0.18)起,ADM-BSR算法的計(jì)算耗時(shí)將小于貪婪追蹤類(lèi)算法.這直觀說(shuō)明了在高塊稀疏度下,ADM-BSR的計(jì)算優(yōu)勢(shì)十分明顯;2)在重構(gòu)誤差方面,總的來(lái)看,Block OMP的MSE最大,ADM-BSR和Block CoSaMP的相對(duì)較小.對(duì)于ADM-BSR和Block CoSaMP,在低塊稀疏度時(shí),兩者的MSE非常接近,但隨著塊稀疏度增加,ADM-BSR的MSE小于Block CoSaMP的.
作者研究了塊稀疏信號(hào)的重構(gòu)問(wèn)題.基于塊稀疏重構(gòu)的l2/l1模型,利用ADM提出一種ADM-BSR算法,分析了算法的計(jì)算復(fù)雜度和收斂性.通過(guò)仿真實(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.