一種用于壓縮感知理論的投影矩陣優化算法

吳光文 張愛軍 王昌明

1 引言

文獻[1～3]在 2006年提出了壓縮感知理論（Compressed Sensing, CS）。該理論指出，對于稀疏信號，可以通過少量的觀測數據恢復原信號，且觀測數據的量遠小于傳統的香農-奈奎斯特采樣定理規定的最小數據量，該理論使得高分辨率信號的采樣能夠突破當前硬件速度限制，實現高速采樣和數據壓縮。

自然界中的多數信號具有一定的規律性，這種規律性體現在數學形式上就是：用一組特定的基底表示信號時，變換系數是稀疏的（絕大多數的系數為零或者絕對值非常小）。稀疏性是CS的理論基礎，利用信號稀疏性，構造 CS系統包括兩個步驟：設計感知機構（編碼機構）和選取合適的重構方法（解碼機構）[4]。精確信號重構（解碼）所需要的數據量依賴于感知矩陣與稀疏字典之間的相關性和信號自身的稀疏度[5]，而投影矩陣的性能是決定編碼質量的關鍵。作為 CS理論研究的一個重要方向，現有的文獻對投影矩陣的約束條件進行了研究，這些約束條件主要包括約束等距性質，零空間性質和相關系數[611]-。然而，判斷投影矩陣是否具備約束等距性質和零空間性質是組合復雜度問題，實際用于投影矩陣性能分析非常困難[12]。

文獻[13,14]引入了相關系數的概念，在CS理論中，相關系數的含義是指投影矩陣與稀疏字典之間列向量內積的最大值，其物理意義是兩者最差相似性[6]。文獻[6～10]的研究表明，通過減少投影矩陣與稀疏字典的相關系數可以提高 CS的重構性能，即相關系數越小，精確重構信號所需的觀測值數目越少，信號適應的稀疏度范圍越大。

文獻[6]提出了一種閾值平均系數方法表示投影矩陣與稀疏字典的相關性，并通過線性收縮 Gram矩陣中絕對值大于限定閾值的非對角元的方法迭代減小相關系數，取得了較好的實驗效果。但是，該方法的計算過程只有收縮系數非常接近1時才近似凸函數[6]。而當線性收縮系數近似為1時，收縮速度非常慢。另外，該算法在使用處理過的Gram矩陣計算降階投影矩陣時會產生干擾，產生絕對值較大的相關系數[4]，本文4.1節的實驗結果表明該算法產生的相關系數比較大，甚至比優化前更大。文獻[7]用等角緊框架 （Equiangular Tight Frame, ETF）Welch界作為投影矩陣和稀疏字典相關系數的極限最小值，將ETF作為優化目標建立一個變形的凸集合，使用梯度下降法逼近最優凸集合中的矩陣，計算投影矩陣。相對于文獻[6]算法，此算法更穩定，但是在算法中步長因子β要根據經驗確定，β的選擇對算法影響比較大[4]。

本文對投影矩陣與稀疏字典的相關性進行研究，優化產生具有最小相關性的投影矩陣，提高信號重構的精度并降低編碼機構對信號稀疏度的要求。設計連續可導的閾值函數，對Gram矩陣的非對角元進行收縮處理。使用基于沃爾夫條件的梯度下降法迭代逼近最優投影矩陣，提高算法的穩定度。

2 壓縮感知的基本理論

CS理論定義[4]：假設信號 x ∈Rn在一組字典Ψ =（Ψ1,Ψ2, …,Ψn）上具有s（ s ＜ ＜n）稀疏度，通過其在投影矩陣 Φ = （ φ1, φ2, …,φn） 上的 m （ m ≥s） 個線性觀測值y （ i）= x, φi,i ∈ { 1,2,… ,m }能夠獲得精確重構的過程，即

式中，Φ為投影矩陣，Ψ為稀疏字典，y為觀測值。

由于信號是稀疏的且信號恢復為病態求逆過程，信號重建可以理解為尋找最小0l范數解的過程。初始信號的稀疏表示向量0α滿足：

選定投影矩陣Φ，信號的重構問題轉化為求解式（3），使用重構算法（如BP, OMP）能夠計算出重構信號的稀疏表示α，進而重構原始信號x[6]：

由式（2），投影矩陣和稀疏字典的相關性為信號準確重建提供保證，相關系數越小，精確重構信號所需的觀測值數目越少，或者說適應信號的稀疏度范圍越大。投影矩陣和稀疏字典的相關系數定義為[6]

其中D是感知矩陣，理論上說，相關系數{}μ D越小，感知（編碼）原始信號后蘊含的信息量越多，但是相關系數有一個下界，即Welch界，如式（5）所示。

3 投影矩陣的優化算法

3.1 感知矩陣相關系數

式中，th[0,1）∈為閾值，式（6）的含義為絕對值大于等于閾值th的G矩陣的非對角元的絕對平均值，能較好地評價投影矩陣的總體相關性。

3.2 Gram矩陣的閾值函數

應用閾值函數處理投影矩陣與稀疏字典生成的Gram 矩陣非對角元，可以使其逼近最優目標。因此，閾值函數對投影矩陣的優化過程起到關鍵作用，文獻[6]中所用的線性優化方法，只有當線性收縮系數γ非常接近1時，才能保證優化過程是收斂的。本文提出一種閾值函數，此函數在單極性區間為凸函數，在（-1,1）區間連續且可導，能夠在保證收斂速度的情況下盡量保留原始數據特征。

式（7）中有兩個參數，th和m, th選恒定值且須th≥ μwelch, μwelch為Welch界。另外式（7）滿足

式（8）表明，式（7）不僅在閾值±th處連續，而且可導。m為大于1的可變參數，當m由小變大時，函數的壓縮程度變大，不同m值對應的函數如圖 1所示。圖1中，th = 0 .2，當m值越大，原始數據的特征改變越大，整體迭代算法的速度越快；m值越小，原始數據的特征改變越少，整體迭代算法的速度越慢，穩定性越好。

圖2顯示了不同m取值時迭代收縮1000步對應的閾值平均相關系數的曲線，稀疏字典為200×400的隨機矩陣，原始的投影矩陣為30×200的高斯分布的隨機矩陣 Φ ，閾值選擇 t h =0.2 ≥ μwelch=0.1758。m取值分別為1.4, 2.0, 8.0。從圖2中發現，3個m取值對應的曲線都收斂，且隨著m值變大收斂速度變快。大量實驗表明，m在區間[1,+∞）中取其它值時，函數曲線變換趨勢和圖2類似。

3.3 更新投影矩陣

優化投影矩陣需要在每個迭代步驟中從收縮非對角元的Gram矩陣中解算出對應的投影矩陣Φ。文獻[6]先用Cholesky分解Gram矩陣，再用Moore-Penrose逆求投影矩陣，此方法會引入干擾誤差，產生較大的相關系數。文獻[7]提出梯度逼近法求解投影矩陣Φ，運算的精度高，算法更加穩定。但是文獻[7]的方法受迭代步長的影響，本文引入基于沃爾夫條件的梯度下降法，用快速逼近每次迭代過程產生的Gram矩陣的方法計算投影矩陣Φ。

式中，kG 為第k步迭代收縮后的Gram矩陣，用梯度下降法由kΦ計算1k+Φ。定義函數，

圖1 Gram矩陣的閾值函數

根據文獻[15]，

為了區別于梯度下降法中迭代過程的標識，將式（11）簡寫作 ? J = 4 Φ Ψ （ ΨTΦTΦ Ψ -）ΨT，梯度下降法的迭代步驟為

其中迭代運算的初始值 Φ（0）=Φk，當式（9）迭代運算達到精度要求時，令=Φ（i+1）, i為梯度下降法迭代的步數。計算式（12）的一個關鍵步驟是求解步長αi，文獻[7]已經證實J（Φk） 并非簡單的二次函數。因此，一些簡單的確定迭代步長的方法在這里不適用。引入文獻[16]提出的基于沃爾夫條件的梯度下降法來解決這個問題。梯度下降法求極值的迭代過程的實質就是保證式（13）成立，

保證式（13）成立的一個實用的規則是沃爾夫條件：

式 （14）中 ， 0 ＜ σ ＜ δ ＜ 1 為 給 定 的 常 數 ，（d（i））T? J（Φ（i）） 為函數J（ Φ） 在d（i）方向上的方向導數。式（14）的第1個不等式叫做Armijo規則，該規則的數學含義是：步長αi越大，函數J（ Φ）的值改變越大。式（14）的第2個條件的含義是：J（ Φ）在點上的方向導數的值必須大于等于δ倍的初始值 Φ（i）的導數值。

圖2 閾值平均相關系數

在實用中σ取值應非常小，但δ要取大值[16]，本文取值 σ =10-4, δ=0.9，使用回溯算法確定符合沃爾夫條件的迭代步長αi，算法如圖3所示。

根據圖3所示算法求得迭代步長αi，將步長代入式（12）即可用梯度下降法求得滿足最小誤差的Φ（i+1），此即優化后的投影矩陣Φ 。投影矩陣優化

k+1算法分為兩大步，第1步更新Gram矩陣，第2步更新投影矩陣。將這兩步細化后的具體的步驟如圖4所示。迭代次數Iter可以是確定的數值，還有另一種方法控制迭代結束，就是判斷連續兩次相關系數之間的差值，如果差值小于設定的值，結束迭代過程。

圖3 回溯算法求步長流程圖

3.4 算法分析

收斂性分析 沃爾夫條件是對梯度下降法中線性搜索方向和步長的限制條件，在式（14）第 1個不等式中，α 為步長，（d（i））T? J （Φ（i）） 為方向導數，這

圖4 投影矩陣優化流程

i個不等式保證函數J的下降幅度與步長αi和方向導數 （d（i））T? J （Φ（i）） 這兩個量是比例關系，即保證每次迭代下降足夠大的量，此式又稱作阿米賀條件（Armijo condition）。式（14）中的第 2 個不等式是曲率條件，保證目標函數在步長αi處的斜率是（d（i））T? J （Φ（i）） 的 δ倍。只要滿足沃爾夫條件，算法的收斂性就會滿足，即沃爾夫條件是梯度下降算法有效性的充分條件，文獻[17]給出了滿足沃爾夫條件的步長αi的存在性證明。

算法時間復雜度 本文用計算機浮點數運算的次數表示算法的時間復雜度。假設投影矩陣Φ為N×M矩陣，Ψ為M×M矩陣，則Gram矩陣計算中 ， 浮 點 運 算 的 次 數 為 M2（2N - 1 ）+ 2M N（2M- 1 ）；更新Gram矩陣的運算中，使用式（7）對非對角元素壓縮，在選定閾值函數后，,（m - 1 ）th是常數，每個非對角元的壓縮需要執行兩次乘法運算，一次開方運算和一次減法運算，此處Gram矩陣為M×M矩陣，但對角元素不參與運算，所以浮點數運算的最大次數是 4 （M2- M ），此時對應所有非對角元全大于閾值th的情況；更新投影矩陣的迭代運算中，假設梯度下降法的迭代步數為K，浮點數運算的次數為 1 2K M2N - 3 KMN。算法時間復雜度可以表示為： I （12K M2N - 3 K MN + 3M2+ 6 M2N - 4 M - 2 M N）, I = I ter 是算法的總體迭代次數。

算法性能邊界 使用本文提出的Gram矩陣非對角元收縮算法對投影矩陣進行優化時，投影矩陣不會產生大量的零值元素，重構算法的適用種類比較多。相關實驗表明，該算法可以配合常用的重構算法對信號壓縮感知處理，如 BP算法，匹配追蹤（Match Pursuit, MP）算法，OMP算法，壓縮采樣匹配追蹤（Compressive Sampling Matching Pursuit,CoSaMP）算法，正則化正交匹配追蹤（Regularized Orthogonal Matching Puisuit, ROMP）算法等，均可以重構優化投影矩陣觀測的原始信號。

4 仿真實驗

為了驗證本文算法的有效性，選擇高斯分布的隨機矩陣作為原始投影矩陣，分別使用文獻[6]方法、文獻[7]方法和本文方法優化原始投影矩陣，然后用這3種優化的投影矩陣和原始的投影矩陣分別進行實驗。首先考察各種投影矩陣的相關系數，進而用各種投影矩陣壓縮感知稀疏向量、1維和2維信號，考察不同投影矩陣對應的重構信號誤差情況。

4.1 優化相關系數實驗

選擇初始投影矩陣Φ為30×200的高斯分布隨機矩陣，Ψ為200×400的稀疏字典。文獻[6]優化方法選擇：閾值th = 0 .2，收縮因子 γ : γ1=0.55,γ2= 0 .75,γ3= 0 .95，迭代1000次；文獻[7]方法選擇：th = 0 .2，迭代次數100，每次迭代過程中求最佳觀測矩陣的迭代次數 50，迭代步長 β : β1=0.01,β2= 0 .02；本文方法選擇：th = 0 .2, σ = 0 .0001,δ= 0 .9，收縮系數 m = 2 ，迭代次數為Iter = 1 00，求最佳觀測矩陣的迭代次數為50。實驗結果為各種算法迭代結束后相關系數μmax和大于設定閾值的相關系數平均值μav，如表1所示。

表1中的實驗結果表明，各種算法優化原始投影矩陣后，大于設定閾值的相關系數平均值μav都有不同程度的下降，說明這3種優化算法都能夠實現降低平均相關性的目的。其中經文獻[6]方法優化投影矩陣后，μmax相對其它兩種優化方法偏大，在γ1=0.55時甚至比原始投影矩陣的μmax更大。即用文獻[6]方法優化投影矩陣后，代表相關性的參數μmax在某些情況下變大了，而非變小。本文在γ取值在 0.50～0.99之間，以不同的步長改變γ的值，做了大量的實驗，實驗結果表明，當0.50 ≤ γ ≤0.70時，μmax比優化前的值大；當0.75≤ γ ≤ 0 .95時，μmax比優化前的值小，且隨著γ變大而μmax逐漸變大；產生這種現象的原因主要是在處理的Gram矩陣計算降階投影矩陣時產生了干擾，因而產生了絕對值較大的相關系數。如果將文獻[6]方法優化的投影矩陣用在信號處理中，信號處理的總體效果變好，因為μav變小了；但是有少數的局部成分會比優化前更差，因為μmax比較大。文獻[7]方法優化投影矩陣后，μmax和μav都有較大的下降，體現了較好的性能，如果將文獻[7]方法優化的投影矩陣用在信號處理中，處理的總體效果和局部效果都比原始投影矩陣好。但是，當迭代步長β取不同值時，對μmax和μav的影響較大，實際應用中，使用者需要根據經驗和實際信號的特征選擇最佳迭代步長，且計算過程中容易陷入局部最優的情況而找不到全局最優值。和其他方法比較，本文提出的優化算法在迭代步數相同、運算量不增加的情況下，maxμ和avμ的值最小，優化的效果最好。

4.2 隨機稀疏向量實驗

為了驗證投影矩陣優化算法在整個壓縮感知過程中的有效性，使用確定稀疏度的稀疏向量作為原始信號。先使用感知矩陣對原始信號編碼產生觀測信號，再使用BP算法和OMP算法對觀測信號進行重構，根據重構的誤差驗證編碼的效果，進而檢驗投影矩陣的性能。實驗選取100000個長度是120，稀疏度是 4的向量（非零值隨機分布），稀疏字典為80120×的單位隨機矩陣，觀測值n取區間[16,40]內的偶數，原始投影矩陣為 n×80的高斯分布隨機矩陣，另外原始投影矩陣分別經過文中提到的3種方法優化，共4種投影矩陣用于實驗，選擇迭代次數為1000次。縱坐標為重構誤差的對數表示，實驗所用的一部分代碼來自SparseLab工具箱[18]。結果如圖5所示。

圖5中的曲線變化趨勢表明，隨觀測值增加，4種投影矩陣對應信號重構誤差逐漸減小。比較圖5（a）4種投影矩陣對應的曲線，3種優化后的投影矩陣對應的效果均優于原始投影矩陣，說明這3種優化方法對投影矩陣的性能具有改善作用。而本文方法的效果比其它算法效果更好。圖 5（b）中的曲線變化趨勢同圖 5（a），說明了本文方法相對于其它算法的優勢不受重構算法的影響。

表1 不同投影矩陣對應的相關系數

圖5 4種投影矩陣對應重構誤差隨觀測值n變化曲線

4.3 信號仿真實驗

使用不同投影矩陣對1維和2維信號分別進行壓縮感知實驗，對觀測值使用 OMP算法重構，比較各種算法的優化效果。1維信號選擇 blocks,Doppler信號，均使用小波變換稀疏原始信號。其中 blocks信號用 haar小波，Doppler信號用Symmlet8小波。比較重構信號，本文方法優化的投影矩陣對應的重構信號最逼近原始信號。其中，重構的Doppler信號如圖6所示，因blocks信號的重構效果圖可以得出同樣的結論，這里省略。

用信噪比和均方根誤差兩項指標對不同算法進行比較，如表 2，表 3所示。從表中看出，相對于其它投影矩陣，使用本文的方法對投影矩陣優化后，重構信號的信噪比最高，均方根誤差最小。

選擇大小為256×256的Lenna灰度圖像，首先使用離散 Symmlet8小波對原始圖像進行稀疏化處理，再使用不同投影矩陣進行壓縮感知處理，最后進行圖像重建，不同的投影矩陣處理的結果如圖7和表4所示。2維信號處理的結果同樣表明本文算法相對于其它優化算法有較大的優越性。

表2 不同投影矩陣對應的重構blocks信號指標

表3 不同投影矩陣對應的重構doppler信號指標

表4 不同投影矩陣對應的重構Lenna圖像指標

5 結束語

本文對壓縮感知投影矩陣的優化問題進行了研究，提出了連續可導的收縮閾值函數，保證了收縮Gram 矩陣非對角元的迭代過程的收斂性。用梯度下降法逼近投影矩陣的方法提高了優化過程中投影矩陣算法的精度和穩定性，克服了用Moore-Penrose逆來求解投影矩陣產生的干擾誤差問題。使用基于沃爾夫條件的梯度下降法，解決了迭代步長對算法性能影響較大的問題。通過對隨機稀疏向量使用各種投影矩陣的壓縮感知實驗，考查使用BP和OMP算法重構信號時產生的誤差，證明了本文所提算法在提高信號重構性能的優越性。通過對基本的小波測試信號和圖像進行壓縮感知處理，使用 OMP算法對感知后的信號進行重構，實驗結果表明本文的算法在壓縮感知處理實際信號時，優于實驗中的其它算法。本文的研究工作還有許多有待改進的地方，例如如何優化收縮閾值函數的數學表達式，以構造相關性更小的投影矩陣；研究如何實現投影矩陣優化和感知信號重構的協同問題。

圖6 不同投影矩陣對應的重構Doppler信號

圖7 不同投影矩陣對應的重構Lenna圖像

[1] Donoho D L, Elad M, and Temlyakov V N. Stable recovery of sparse overcomplete representations in the presence of noise[J]. IEEE Transactions on Information Theory, 2006,52（1）: 6-18.

[2] Candes E J, Romberg J K, and Tao T. Stable signal recovery from incomplete and inaccurate measurements[J].Communications on Pure and Applied Mathematics, 2006,59（8）: 1207-1223

[3] Candes E J and Tao T. Near-optimal signal recovery from random projections: universal encoding strategies[J]. IEEE Transactions on Information Theory, 2006, 52（12）:5406-5425.

[4] 鄭紅, 李振. 壓縮感知理論投影矩陣優化方法綜述[J]. 數據采集與處理, 2014, 52（1）: 43-53.Zheng Hong and Li Zhen. Survey on optimization methods for projection matrix in compress sensing theory[J]. Journal of Data Acquisition and Processing, 2014, 52（1）: 43-53.

[5] 戴瓊海, 付長軍, 季向陽. 壓縮感知研究[J]. 計算機學報,2011, 34（3）: 425-434.Dai Qiong-hai, Fu Chang-jun, and Ji Xiang-yang. Research on compressed sensing[J]. Chinese Journal of Computers,2011, 34（3）: 425-434.

[6] Elad M. Optimized projections for compressed sensing[J].IEEE Transactions on Signal Processing, 2007, 55（12）:5695-5703.

[7] Abolghasemi V, Ferdowsi S, and Sanei S. A gradient-based alternating minimization approach for optimization of the measurement matrix in compressive sensing[J]. Signal Processing, 2012, 92（3）: 999-1009.

[8] 李佳, 王強, 沈毅, 等. 壓縮感知中測量矩陣與重建算法的協同構造[J]. 電子學報, 2013, 41（1）: 29-34.Li Jia, Wang Qiang, Shen Yi, et al.. Collaborative construction of measurement matrix and reconstruction algorithm in compressive sensing[J]. Acta Electronica Sinica,2013, 41（1）: 29-34.

[9] Zhang Qi-heng, Fu Yu-li, Li Hai-feng, et al.. Optimized projection matrix for compressed sensing[J]. Circuit System Signal Processing, 2014, 33（5）: 1627-1636.

[10] Xu Jian-ping, Pi Yi-ming, and Cao Zong-jie. Optimized projection matrix for compressive sensing[J]. EURASIP Journal on Advances in Signal Processing, 2010,doi:10.1155/2010/560349.

[11] 林波, 張增輝, 朱炬波. 基于壓縮感知的 DOA估計稀疏化模型與性能分析[J]. 電子與信息學報, 2014, 36（3）: 589-594.Lin Bo, Zhang Zeng-hui, and Zhu Ju-bo. Sparsity model and performance analysis of DOA estimation with compressive sensing[J]. Journal of Electronics & Information Technology,2014, 36（3）: 589-594.

[12] Donoho D L. For most large underdetermined systems of linear equations the minimal l1-norm solution is also the sparsest solution[J]. Communications on Pure and Applied Mathematics, 2006, 59（6）: 797-829.

[13] Donoho D L and Stark P B. Uncertainty principles and signal recovery[J]. SIAM Journal on Applied Mathematics, 1989,49（3）: 906-931.

[14] Donoho D L and Elad M. Optimally sparse representation in general （nonorthogonal） dictionaries via minimization[J].Proceedings of the National Academy of Science, 2003, 100（5）:2197-2202.

[15] Petersen K B and Pedersen M S. The matrix cookbook[OL].http://www.matrixcookbook.com, 2013.10.

[16] Barth T J, Griebel M, Keyes D E, et al.. Scientific computing with MATLAB and octave[OL]. http://www.springer.com/series/5151, 2013.12.

[17] Jorge N and Wright S J. Numerical Optimization Theoretical and Practical Aspects[M]. 2nd Edition, New York: Springer-Verlag Berlin and Heidelberg GmbH & Co. K, 2006: 30-60.

[18] Stodden V and Donoho D. SparseLab21-core[OL]. http://sparselab.stanford.edu, 2013.10.