基于非凸復合函數的稀疏信號恢復算法

周潔容 李海洋 凌 軍 陳 浩 彭濟根

根據奈奎斯特(Nyquist)采樣定律,想要實現信號的無失真輸出,采樣頻率必須在信號帶寬的兩倍以上.然而,在大多數情形下,按此定律采樣所獲得的信息是冗余的,這不僅造成采樣的浪費,而且處理較大帶寬的信號時會給硬件系統帶來巨大壓力.壓縮感知(Compressed sensing,CS)[1]的出現使該問題的解決成為可能.壓縮感知是一種新型的信號采樣及重構理論,利用少量的測量值就可以實現稀疏或可壓縮信號的精確重構.而稀疏信號的重構在眾多科學研究和工程應用中十分重要,如物理學中的量子態層析成像[2]、天體物理學成像[3]、磁共振成像[4]、信號處理[5]、雷達成像[6]等,CS 理論的引入加快了上述應用的研究和發展.

數學上,稀疏信號重構是指從欠定線性系統bAx中恢復原始信號x,其中b∈Rn是測量向量,A∈Rn×m(n

其中,‖x‖0表示向量x中非零元素的個數[7].

已經證明,(P0)模型的直接求解是NP(Nondeterministic polynomial)難的[8],其計算量會隨著稀疏向量維數的增加而增大,模型的抗噪能力也很差.為此人們提出了多種方法對 (P0) 模型進行求解,主要方法可歸類于啟發式方法、凸松弛方法、非凸松弛方法三種類型.典型的啟發式算法有正交匹配追蹤[9]、閾值算法[10]、子空間追蹤算法[11]、分級正交匹配追蹤[12]等,這些算法重構理論簡單、速度較快,但往往需要更多觀測,算法的重構精度較低、收斂速度較慢,并且在有噪聲情況下信號的恢復精度較低,因此其應用范圍有限.典型的凸松弛算法有梯度投影算法[13]、BP (Basis pursuit)算法[14]、BPDN(Basis pursuit denoising)算法[15]、IRLS (Iterative reweighed least squares)算法[16]、Bregma 迭代法[17]等,其中最經典的是BP 算法,文獻[18]證明了在測量矩陣滿足有限等距性質[19]的條件下BP 模型與(P0)模型等價,但在多種情形下BP 模型中的L1范數不能充分反映信號的稀疏性特征,往往難以獲得稀疏解.為此,人們提出了用非凸泛函替代L0范數的方法,稱之為非凸松弛方法.

相較于L1范數,許多非凸泛函能更好地近似L0范數,從而更好地反映信號的稀疏性特征.典型的非凸松弛算法有NSL0 (Newt-on smoothL0norm)算法[20]、SL0 (SmoothedL0)算法[21]、CTNRAL0 (Composite trigono-metric function nullspace reweighted appro-ximateL0norm)算法[22]、SCSA (Successive concave sparsity approximation)算法[23]等.易見,松弛方法的核心是尋找L0范數的逼近函數,通過極小化該逼近函數尋得最稀疏解.顯然這種近似模型對 (P0) 模型的逼近性能取決于所選取的近似函數對L0范數的逼近程度.因此,如何構造具有更優逼近性能的近似函數已成為稀疏信號重構問題研究中的重要問題.

最近,文獻[23?24]分別構造了兩種近似函數gσ(x)1?e?|x|/σ、hσ(x)|x|/(|x|+σ)以實現對L0范數的逼近,取得了較好的重構效果.一個有趣的想法是,如果將上述每個泛函對信號的作用看成是一次前向(Forward)處理過程,那么我們是否可以通過泛函的復合實現對信號的深度作用,從而提高信號重構的性能? 基于這種想法,本文將上述兩個泛函進行復合,構建一種新的非凸松弛模型,給出該模型的理論分析,對該模型提出一種新的近似算法,并通過數值實驗驗證算法的有效性以及相較于SL0 算法、IRLS 算法、BP 算法、SCSA 算法等經典算法的優越性.

1 一種用于逼近 L0 范數的非凸復合指數泛函

設δ為Kronecker 函數,即

這樣,下面的優化問題可看成是 (P0) 問題的近似模型

文獻[21,23?24]分別采用以下三種DA 函數來實現對L0范數的逼近

并通過對相應近似模型的求解實現了信號的精確重構.

顯而易見,近似式(3)對 (P0) 模型的逼近性能取決于所選取的DA 函數對L0范數的逼近程度.因此,如何構造具有更優逼近性能的DA 函數已成為稀疏信號重構問題研究中的重要課題.注意到,上述DA 函數hσ(x) 和gσ(x)都能很好地實現對L0范數的逼近.一個自然的問題是,是否可以將這兩個函數進行復合以實現對L0范數的深度逼近? 針對這個問題,本文對上述兩個函數hσ(x) 和gσ(x) 的復合進行考察

其中,“?”表示函數復合運算,顯然這是一個DA 函數.不僅如此,可以從幾何圖像和理論分析上說明,在固定參數σ后該DA 函數相對于hσ(x)、gσ(x)、pσ(x)具有更好的對L0范數的逼近性能.

1)幾何圖像分析

圖1 展示了pσ(x)、hσ(x)、gσ(x) 以及復合函數fσ(x)這4 種DA 函數的幾何圖像,其中參數σ0.1.可以看出,在x0附近,pσ(x) 函數曲線最平坦,而本文所提出的復合函數fσ(x) 的圖像曲線相對其他三種函數最為陡峭,這表明函數fσ(x)具有對L0范數最好的逼近效果.

圖1 4 種函數在 σ0.1 時的一元函數分布Fig.1 The unary distribution of the four functions at σ 0.1

圖2 展示了pσ(x)、hσ(x)、gσ(x)和fσ(x) 4 種DA 函數在二維情形的等高幾何圖像,即pσ(x)hσ(x)gσ(x)fσ(x)1的圖像,其中參數σ0.1.可以看出,在4 種函數與直線x1x2的交點中,函數fσ(x) 與x1x2所表示的直線得到的交點的分量絕對值最小,這表明復合指數函數fσ(x) 能更好地逼近L0范數.

圖2 pσ(x)、hσ(x)、gσ(x)和函數 fσ(x)在 σ0.1 時的二元函數分布Fig.2 The bivariate distribution ofpσ(x),hσ(x),gσ(x)and the function fσ(x)atσ 0.1

2)理論分析

上述從幾何直觀上展示了復合指數函數相對于其他DA 函數的更好逼近性能.不僅如此,也可以從理論上證明,對區間 (0,0.9)內任意非零參數σ,復合指數函數相對于pσ(x)、hσ(x)、gσ(x) 都具有更好的逼近性能.不失一般性,僅在第一象限內進行討論,即假定x ≥0.對區間內任意σ(0),因為

首先,固定函數值比較對應變量的大小.

假設hσ(x)gσ(x)r,其中 0

因為

所以有R(r)為增函數,即對任意r有R(r)>R(0)0,得xh >xg.這說明相較于hσ(x),函數gσ(x) 更貼近坐標軸.

下面固定變量x,比較函數值pσ(x)、gσ(x)、fσ(x)的大小.因為

所以,為討論gσ(x)、pσ(x) 的大小關系,只需要比較指數的大小.易得x趨于 0(至少在x ≤2σ時),對任意非零σ恒有

即gσ(x)≥pσ(x).

同理,對

進行討論,由于

所以,在x0附近(至少在|x|≤0.1 范圍內)有fσ(x)≥gσ(x).因此有

說明fσ(x) 的函數曲線更陡峭.

利用函數的對稱性,其他象限也可類似討論.所以由上述函數的幾何圖像和理論分析也能得出fσ(x)相對于hσ(x)、gσ(x)、pσ(x) 具有更強的稀疏性.

2 基于復合指數函數的稀疏模型

以上分析表明,由式(7)定義的復合函數fσ(x)具有對L0范數更優的逼近性能.為此,考慮由該函數所誘導的優化模型

下面對優化模型 (CE) 與 (P0) 之間的關系進行研究,為此首先證明以下引理.

引理 1.設模型 (CE)的最優解為x?,其稀疏度為k,則矩陣A中對應解x?支撐集的子矩陣Ak是列滿秩的,其中x?的支撐集是指集合i1,···,m}.

證明.利用反證法,假設矩陣Ak不是列滿秩的,即矩陣的列向量線性相關,則存在非零向量v,其支撐集包含在向量x?的支撐集內,并且使得Akv0.顯然,向量x?±v滿足約束條件Axb.不失一般性,假設向量v中最大的分量絕對值不超過向量x?中最小的分量絕對值.

根據函數fσ(x) 的嚴格凹性,有

這表明Fσ(x?+v)

上述引理1 表明,模型 (CE)的最優解x?滿足即模型的最優解包含在有限集Γ是列滿秩矩陣}中.

定理1.對任意感知矩陣A和測量向量b,存在一個依賴于A和b的常數σ(A,b),當0<σ <σ(A,b)時,模型 (CE)與模型 (P0) 等價.

證明.要證明模型 (CE)與模型 (P0) 等價,則只需證明模型 (CE) 的最優解為模型(P0) 的最優解即可.

根據引理1,可以得到以下等價關系

說明對任意σ,(CE)模型的最優解x?∈Γ,其中 Γ 為有限集.

下面證明存在一個常值σ(A,b),當0<σ <σ(A,b)時,存在唯一一個點∈Γ均是 (CE) 最小化模型的解,下面利用反證法進行求證.

3 模型轉換

為克服函數Fσ(x) 在原點不可微帶來的不便,可以將稀疏優化模型 (CE) 轉換為以下優化問題.令未知量xu ?v,其中u,v∈Rm,u擁有x中所有正元,其余元素為零;v擁有x中所有負元的絕對值,其余元素也為零.用z[uT,vT]T∈R2m表示拼接向量.經過替換,易得

此時,約束條件Axb轉換為[A,?A]zb.那么無噪模型 (CE) 轉換為

為了便于求解式(14),利用文獻[25]中介紹的MM優化方法. MM 優化方法指當目標函數較難實現優化時,通常可以選擇更容易優化的替代目標函數,當替代函數滿足一定的條件時,其最優解能夠無限逼近原目標函數的最優解.通過利用凹函數Fσ(z)的一階判別條件以及 MM 優化方法對目標函數Fσ(z) 進行放縮,可得

記

即

其次,因為

有

再利用假設式(20)有

由此可知,對任意ε2>0,存在K2>0,對任意k >K2,有

故結合式(23)、式(25)可知,對任意z ∈?,對于式(22)可得

即對任意z ∈?,有任意ε3>0,存在K3>0,對任意k >K3,有

令εmin(ε1,ε2,ε3),Kmax(K1,K2,K3),當k >K時,對不等式(29)有

由式(30)可知

對左端加減同一項,有

對式(31)移項,有

再對式(32)不等號左端進行放縮,得

最后結合式(32)、(33)得

優化模型(21)簡記為

4 算法設計

為給出模型(14)的有效求解算法,需要對算法初始值的選擇進行分析,首先給出復合函數的如下性質.

其中,關于 e?|x|/σ(|x|+σ)/(|x|+σ)3、e?|x|/σ(|x|+σ)/(|x|+σ)4的任意階導數滿足有關 (|x|+σ)?1的次方數不低于 3次,則(n ≥2)關于參數σ?1的次方數不低于 3 次.由此可得

利用式(35)、(36),有

引理2 表明,在時,最小化Fσ(x) 近似于最小化‖x‖1. ‖x‖1的凸性使得算法陷入局部極值的可能性減少,從而提高算法重構精度.鑒于此,本文在設計算法時取初始值為以下問題

的解.

在此基礎上,提出如下算法.

步驟1.輸入:A,b,ε1,ε2,ε3.

步驟2.初始化:

步驟3.算法迭代:

步驟4.輸出稀疏向量:x=xk.

其中,步驟 3 算法迭代環節由內、外兩部分循環構成.外循環為步驟 a)～b) 以及 f)～h),循環利用參數σ實現復合函數對L0范數的逐次逼近.內循環即步驟c)～e),結合外點罰函數法[26]和共軛梯度法迭代求解模型(14),具體為利用外點罰函數法引入正則參數實現模型(14)的無約束轉換,再通過共軛梯度法求解無約束優化問題.最后在步驟j)～k),針對共軛梯度法求得的非稀疏解,利用加權L1范數最小化[27]進行稀疏化處理.其中,d1和d2分別表示外循環和內循環連續迭代的解之間的相對誤差,并用于判斷循環是否停止.值得關注的是,整個算法的核心環節即步驟d),是有關本文模型(14)的凸優化求解過程,具體步驟如下.

模型利用外點罰函數法引入正數M,本文取M1,其放大系數取為 5,將式(15)轉化為以下無約束問題

其中

式(38)中,I是大小為 2m×1 的單位向量.

其中,V(z)(min(0,z1),···,min(0,z2m))T.

其中,步長因子

下降方向

綜上,本文提出的基于非凸復合函數的稀疏信號恢復算法(Non-convex composite sparse,NCCS)參考了文獻[23]中的SCSA 算法,主要包括三個部分.一是算法的初始值選擇,取值為L1范數最小化問題的解,參考文獻[29],其計算量為O(m3).二是無約束問題的求解,通過外點罰函數法和共軛梯度法迭代求解最優值.由式(39)～式(44)可以看出,第二部分每次迭代的主要操作是內循環的兩部分矩陣乘法,即

計算復雜度分別為 O(4m+2mn)和O(4m2+4m2n),則第二部分的計算量為O(l(4m+2mn+4m2+4m2n)),其中l表示內循環的迭代次數.三是迭代結果的稀疏化,利用加權L1范數最小化[27]對凸優化結果進行稀疏化處理,其計算量為O(m3).

5 實驗仿真與結果分析

為驗證本文提出的算法在重構性能上的優越性,本節設計了幾組有關SL0 算法[21]、BP 算法[15]、IRLS算法[16]、SCSA 算法[23]和本文介紹的NCCS 這5 種算法在重構性能上的對照實驗.在仿真實驗中,感知矩陣A和稀疏原信號x的具體取值為:矩陣A的大小取為 250×500,其中矩陣元素服從零均值、單位方差的高斯分布,矩陣的列具有單位L2范數;稀疏原信號x的維數取為 500,非零分項服從正態分布.實驗將討論信號稀疏度在區間 [20,110] 范圍內各算法的重構情況.對于等式約束Axb,在已知感知矩陣A和測量向量b的情況下,分別利用上述5 種算法對稀疏原信號進行 100 次仿真實驗,討論算法的平均重構性能.為展示本文所提出的模型相對于最新文獻[23]所提算法的先進性,本文對5 種算法采用了與文獻[23]相同的參數選擇,即

1)對SL0算法:令mu2,σmin10?4,c0.8,L8.

2)對IRLS 算法:令p0.5.

3)對BP 算法:令l11.0.1.

4)對SCSA 算法:令ε110?3,ε210?2,c

5)本文NCCS 算法:令ε110?3,ε210?2,ε310?1,c0.1,σ0min{2‖x0‖0,α},其中α是待定數.仿真實驗一將給出模擬實驗,以獲得最佳實驗數值.

實驗中的重構性能包括:

1)重構信噪比(Signal noise ratio):

2)重構誤差(Mean square error):

3)歸一化均方差(Normalized mean square error):

4)支撐集恢復成功率(Recovery success rate of the support set):

仿真實驗一.待定數α的最佳選值.

首先對NCCS 算法中待定數α進行模擬實驗,觀察不同大小的數值α對算法運行時間的影響,以獲得更優的實驗參數.為節省計算時長,考慮α大小為 0.1 到 0.9,間隔為 0.1,分別進行 100 次仿真實驗.從圖3 可以看出,各α值對算法運行時間的影響差距不大,總體趨勢是算法運行時間隨α的增加而減少;在α0.8時,算法在稀疏度k ≤100 時皆能保證運行時間最短.信號的稀疏度k是指信號向量中非零元素的個數,為便于編程,在后續實驗中待定數α全部選定為0.8.

圖3 待定數 α 對NCCS 算法運行時間的影響Fig.3 The influence of undetermined number α on the running time of NCCS algorithm

仿真實驗二.NCCS 算法對稀疏原信號重構的實驗.

圖4 是稀疏原信號以及由NCCS 算法獲得的恢復信號的仿真結果.圖4 的稀疏度取值為k65.由圖4 可見,恢復信號和稀疏原信號基本吻合,實驗結果顯示算法的重構誤差為 2.6933× 10?17. 由此可以看出,本文提出的算法恢復出的稀疏信號很接近稀疏原信號.

圖4 NCCS 算法的一維信號重構仿真圖,信號大小為500×1,稀疏度為65Fig.4 One-dimensional signal reconstruction simulation diagram of NCCS algorithm,the signal size is 500 × 1,the sparsity is 65

仿真實驗三.各算法的重構性能和稀疏度的變化關系.

圖5 為5 組實驗算法的重構誤差 (MSE) 和稀疏度k的變化關系.其中IRLS 算法的重構誤差最高;由于最速下降法在迭代過程中存在“鋸齒效應”[30],實驗中SL0 算法的重構誤差也很大,且SL0 與BP、SCSA 這3 種算法都存在重構誤差隨著稀疏度的增加而增大的趨勢;而NCCS 算法的重構誤差相對最小,整體變化較穩定.

圖5 SL0、IRLS、BP、SCSA、NCCS 5 種算法的重構誤差和稀疏度的變化關系Fig.5 The relationship between the reconstruction error and sparsity of the five algorithms of SL0,IRLS,BP,SCSA,and NCCS

圖6 表示隨著稀疏度在區間 [20,110] 的等量增加,5 種算法的重構信噪比 (SNR) 的變化情況.其中SL0、BP 和IRLS 3 種算法的稀疏度越高,重構信噪比越低;NCCS 算法與SCSA 算法的重構信噪比在稀疏度區間 [20,110] 上表現較穩定.在實驗的各組稀疏度下,NCCS 算法的重構信噪比遠高于SL0、BP 和IRLS 這3 種算法,略高于SCSA 算法.所以使用NCCS 算法對信號的恢復精確度有一定的提高.

圖6 SL0、IRLS、BP、SCSA、NCCS 5 種算法的重構信噪比和稀疏度的變化關系Fig.6 The relationship between the reconstructed signalto-noise ratio and sparsity of the five algorithms of SL0,IRLS,BP,SCSA,and NCCS

圖7 展示5 種算法的重構運行時間和稀疏度k的變化關系.實驗結果顯示各算法的運行時間隨稀疏度的增大而增加,稀疏度區間在 [20,90] 時,各算法的實驗平均運行時間能保持在 2 s 以內.在稀疏度k >90 后,SCSA 算法重構運行時間迅速增加;當稀疏度k >120 后,該算法的平均重構運行時間不低于 20s.本文提出的NCCS 算法由于參數σ值的多次衰減導致算法進行多次迭代,因此也造成了運行時間相較于剩下的3 種算法要高,但對比最新文獻[23]的SCSA 算法,NCCS 算法有更好的實驗效果.

圖7 SL0、IRLS、BP、SCSA、NCCS 5 種算法的運行時間和稀疏度的變化關系Fig.7 The relationship between the running time and sparsity of the five algorithms of SL0,IRLS,BP,SCSA,and NCCS

圖8 是各算法的支撐集恢復成功率 (RSS) 和稀疏度k的變化關系.實驗結果顯示SL0、BP 和IRLS這3 種算法的稀疏向量支撐集恢復成功率隨稀疏度的增加而下降,SCSA 算法在稀疏度k100、k110時,其支撐集恢復成功率分別為RSS0.9980、RSS0.9886,考慮為實驗次數的不足而出現的偏差.相比之下,在NCCS 算法的重構仿真實驗中,RSS值保持穩定,信號的支撐集在整個實驗區間內實現完全恢復.

圖8 SL0、IRLS、BP、SCSA、NCCS 5 種算法的支撐集恢復成功率和稀疏度的變化關系Fig.8 The relationship between the recovery success rate of the support set and sparsity of the five algorithms of SL0,IRLS,BP,SCSA,and NCCS

圖9 以及表1 分別為實驗中5 組算法恢復向量的歸一化均方差 (NMSE) 和稀疏度k的變化關系圖像曲線和具體數值記錄.由圖9 和表1 可見,NCCS算法和SCSA 算法的實驗結果遠小于SL0算法、BP算法和IRLS 算法的實驗值.

表1 5 種算法的歸一化均方差的數值記錄Table 1 Numerical records of the normalized mean square error of five algorithms

圖9 SL0、IRLS、BP、SCSA、NCCS 5 種算法的歸一化均方差和稀疏度的變化關系Fig.9 The relationship between the normalized mean square error and sparsity of the five algorithms of SL0,IRLS,BP,SCSA,and NCCS

綜上所述,本文提出的NCCS 算法的綜合重構性能要比實驗的對照算法好,所以利用NCCS 算法能更好地恢復稀疏信號向量.

6 結論

本文運用MM 優化、外點罰函數法、共軛梯度法等方法,借鑒文獻[23]提出的SCSA 算法思想,提出了一種新的稀疏信號重構算法,稱為NCCS 算法.該算法利用逼近性能更優的非凸復合指數函數實現對L0范數的逼近.仿真實驗表明,本文所提出的NCCS 算法不僅有效可行,而且在無噪聲情況下,對比SL0、IRLS、BP、SCSA 這4 種算法,NCCS 算法在稀疏信號恢復實驗中表現出更優越的重構性能.