壓縮感知中提升測量矩陣稀疏性的等效字典優化設計

陳映瞳

（西安工業大學基礎學院，西安 710021）

引 言

壓縮感知是一種新型信號觀測與恢復方法，在機械故障診斷［1］、材料的沖擊損傷檢測［2］等領域得到了廣泛應用。這種方法的主體是由測量矩陣、稀疏字典和稀疏重構算法組成的壓縮感知系統，其中稀疏重構算法的穩定性是衡量壓縮感知系統重構性能的重要指標。因此，如何提升稀疏重構算法的穩定性是壓縮感知方法的一個重要研究課題。研究表明等效字典的互相干性值越?。此幕ゲ幌喔尚栽綇姡?，稀疏重構算法的穩定性就越強。因此，提升壓縮感知系統重構性能的一條重要途徑就是減小等效字典的互相干性值。這方面的研究工作最早開始于2007年［3?4］，近年來關于此的研究成果仍不斷涌現。文獻［5］把測量矩陣和稀疏字典放在一起進行優化設計，使得等效字典具有更小的互相干性值，稀疏重構算法則具有更強的穩定性，最終使得壓縮感知系統具有更高的信號重構率。文獻［6］在交替投影算法框架下提出了一種優化設計方法。在每步迭代中，首先固定稀疏字典，通過令等效字典的Gram矩陣逼近某類特定Gram矩陣來優化設計測量矩陣，然后固定測量矩陣，利用一組訓練數據，通過極小化稀疏表示誤差在測量前后的能量和，優化設計稀疏字典。該方法每步迭代的兩個子步驟分別采用不同的目標函數，為了簡化設計過程，文獻［7］采用統一的目標函數優化設計測量矩陣和稀疏字典。為了改進文獻［6］方法不能處理海量訓練數據的不足，基于固定稀疏字典優化設計測量矩陣的方法［8］，通過借鑒文獻［9］中的在線字典學習方法，文獻［10］提出了一種高效的等效字典優化設計方法。文獻［11］采用交替投影方法極小化等效字典的Gram矩陣與某類特定Gram矩陣的差，其間采用共軛梯度法更新等效字典。盡管這些方法設計的等效字典在圖像恢復任務上都取得了較好的數值實驗效果，但在進行優化設計時都沒有考慮到等效字典對后續信號重構效率的影響。

還有一些工作簡化了聯合優化設計方式，即在稀疏字典固定的情形下，僅對測量矩陣進行優化設計。文獻［12］采用交替投影方法極小化等效字典的Gram矩陣與某類特定Gram矩陣的差，但在設計等效字典時沒有考慮它對信號重構效率的影響。文獻［8］在極小化等效字典的Gram矩陣與某類特定Gram矩陣的差時，通過對目標函數添加正則項，使得能夠利用設計好的等效字典較好地重構了不是嚴格稀疏的信號。文獻［8］在進行優化設計時，沒有考慮如何才能利用設計結果更有效地重構信號。為此，文獻［13］對文獻［8］中優化問題的變量（即測量矩陣）賦予了稀疏結構，要求等效字典是按行稀疏的測量矩陣（每行的稀疏度是固定的）和某個正交矩陣（例如離散余弦變換矩陣）的乘積。這種方法能有效提升測量矩陣的稀疏性，但提出的結構有可能阻礙了測量矩陣稀疏性的進一步提升。同時，這種方法對每步迭代中的梯度下降采用回溯準則選取步長，導致方法的時間復雜度較高。和文獻［8］類似，為了能夠利用優化設計的等效字典較好地重構不是嚴格稀疏的信號，文獻［14］在極小化等效字典的Gram矩陣與某類特定Gram矩陣的差時，對目標函數添加了一個包含訓練數據稀疏表示誤差的正則項，但在進行優化設計時都沒有顧及等效字典對信號重構效率的影響。文獻［15］利用光滑化技巧把極小化等效字典互相干性值的問題轉化為一個有約束光滑優化問題，然后利用交替投影方法進行求解。這種方法能顯著減小等效字典的互相干性值，但和文獻［8］一樣，在提高信號重構效率上還有進一步提升空間。

從上述方法分析可知：當前的等效字典優化設計方法普遍沒有考慮等效字典對后續信號重構效率的影響。為此，本文準備固定稀疏字典，極小化等效字典的Gram矩陣與某類特定Gram矩陣的差，同時要求測量矩陣滿足某種基于L1范數的不等式約束。利用L1范數的促稀疏性，可在算法進行過程中自動讓測量矩陣中的部分元素歸零，從而提高后續利用等效字典進行信號重構的效率。然后，基于交替投影方法求解提出的最優化問題，并保證算法的收斂性。和相關方法相比，本文提出的優化設計方法著重關注測量矩陣的稀疏性對后續信號重構效率的影響。利用這種方法優化設計等效字典，有望在信號重構階段獲得更高的重構精度和重構效率。

1 本文優化設計方法

1.1 極小化問題

壓縮感知的測量模型是y=Φx+e，其中y∈RM為測量值，Φ∈RM×N(M?N)為測量矩陣，e為測量噪聲，x∈RN為待重構信號，在稀疏字典Ψ∈RN×L(N≤L)下有稀疏表示x=Ψα+n，其中α為稀疏向量，稀疏度（向量的零范數即指它的非零分量個數），n為稀疏表示誤差。則乘積A=ΦΨ被稱為等效字典。

在稀疏字典Ψ固定的情形下，本文通過求解極小化問題（1）優化設計測量矩陣。

問題（1）的解之所以能更穩定地重構稀疏信號且問題（1）的解有一定稀疏性有以下原因。首先，周知L1范數具有促稀疏性，所以通過約束條件可以迫使測量矩陣的部分元素歸零。而且和文獻［13］稀疏化測量矩陣的方式相比，這種方式沒有硬性規定測量矩陣每一行的稀疏度，所以有可能使得測量矩陣具有更強的稀疏性。其次，集合X描述的是一類理想等效字典的Gram矩陣。這類等效字典的每列具有單位L2范數，且它的互相干性值不超過ξ，ξ為M×L維矩陣所能達到的最小互相干性值。根據壓縮感知中常用稀疏重構算法（例如正交匹配追蹤算法（Orthogonal matching pursuit，OMP）和基追蹤去噪算法（Basis pursuit denoising，BPDN）等）的穩定性保證［16］可知，相較其他M×L維等效字典，這類等效字典搭配常用稀疏重構算法構成的壓縮感知系統，能在最壞意義下具有一定穩定性。因此，通過讓等效字典的Gram矩陣盡可能接近集合X，可使等效字典在后續重構階段帶來更高重構精度。

1.2 求解算法及其收斂性

本節給出極小化問題（1）的解法。將極小化問題（1）改寫成如式（2）所示的等價形式。

本文根據交替投影算法框架求解最優化問題（2）。固定稀疏字典Ψ，設在交替投影算法第k步迭代開始時測量矩陣是Φk，在第k步迭代中交替進行如下兩個步驟。

（1）求Gram矩陣G=(ΦkΨ)T(ΦkΨ)在集合X上的投影，即對G按元素執行如下操作

式中sign(?)是符號函數（規定sign(0)=0）。將投影記為G k。

（2）固定G k求解優化問題（2），即求G k在集合Y上的投影。這時，優化問題（2）可以等價寫為

采用梯度投影算法［17］求解問題（4）。把目標函數記為f(Φ)。按文獻［18］中的算法求Φk在集合Z上的投影，集合為了簡便，投影仍記為Φk，把Φk作為問題（4）的初值。然后選定常數s>0，執行若干步如下操作：按文獻［18］中的算法求在集合Z上的投影?。其中，梯度然后，從Φ開始沿方向前進步長α∈(0，1)，α按Armijo準則選取。設參數β、σ∈(0，1)，找到使得式（5）成立的最小非負整數m，然后得到步長α=βm。

執行完若干步迭代后，把新的測量矩陣記為Φk+1。這樣就完成了交替投影算法的第k步迭代。如此進行若干步交替投影后，輸出最終的測量矩陣Φ。

綜上，本文提出的等效字典優化設計方法的步驟如下。

輸入：測量矩陣Φ∈RM×N，稀疏字典Ψ∈RN×L，交替投影算法的迭代步數Talter，梯度投影算法的迭代步數T，參數s>0和β、σ∈(0，1)，Φ的行數M；Ψ的列數L。

fork=1：Talter

計算A=ΦΨ和G=ATA；對G按元素執行式（3）操作；

按文獻［18］中的算法求Φ在集合Z上的投影，投影仍記為Φ；

fort=1：T

計算B=Φ-s?f(Φ)；

按文獻［18］中的算法求B在集合Z上的投影，投影記為；

利用Φ，Φ?，Ψ，G，β和σ，按Armijo準則選取步長α；

end

end

輸出：測量矩陣Φ。

本文提出的優化設計方法通過計算B的每列在L1范數單位球上的投影，求B在集合Z上的投影。在求B中任意一列在L1范數單位球上的投影時，本文具體采用的是文獻［18］中的算法。

下面從兩個方面說明本文提出的等效字典優化設計方法是收斂的。首先，單看每步交替投影用到的梯度投影算法。因為集合Z是有界集，所以梯度投影算法求解問題（4）產生的序列至少會有1個聚點。根據文獻［17］可知，這個聚點是優化問題（4）的穩定點。其次，來看整個交替投影算法。根據文獻［19］可知，因為集合X和Y都是有界閉集（關于Y的閉性，詳見附錄），所以本文提出的交替投影算法是收斂的。具體來說，設X k∈X是第k步交替投影產生的Gram矩陣，是第k步交替投影產生的測量矩陣對應的Gram矩陣，則序列至 少 會 有 一 個 聚 點且 有：（1），即?是?在集合Y上的投影；（3），即?是在集合X上的投影。

為了便于理解，對上述收斂性結果進行簡要解釋?！凹蟉是RL×L的子集”可以一般化為“G是任一距離空間S的任意非空子集”，“集合X是RL×L的子集”可以一般化為“H是任一距離空間T的任意非空子 集”。問 題（2）的 目 標 函 數 是RL×L×RL×L上 的 連 續 實 值 函 數?(Y k，X k)∈G×H，k=2，3，…，“用梯度投影法求解問題（4）進而得Y k+1”可以抽象成一個從G×H到G×H的集值映射F1(Y，X)=M G(X)×{X}，其中M G(X)表示給定X∈H后f(Y，X)在G上的全局極小點集，具體就是Y k+1∈M G(X k)，(Y k+1，X k)∈F1(Y k，X k)。?(Y k+1，X k)∈G×H，k=2，3，…，執行式（3）操作得到X k+1也可以抽象成一個從G×H到G×H的集值映射F2(Y，X)={Y}×M H(Y)，其中M H(Y)表示給定Y后f(Y，X)在H上的全局極小點集，具體就是{X k+1}=M H(Y k+1)（這里用等號，因為式（3）操作的結果是單點集），(Y k+1，X k+1)∈F2(Y k+1，X k)。因此，對所提交替投影算法，如果設X k∈H是第k步交替投影產生的Gram矩陣是第k步交替投影產生的測量矩陣Φk+1對應的Gram矩陣，則序列{(Y k，X k)}∞k=1滿足(Y k+1，X k+1)∈F(Y k，X k)，k=2，3，…，其中是F1和F2的復合，它的對應法則為

有了上述轉化后，就能應用如下定理1得到所提交替投影法的收斂性結果，即序列至少有一個聚點且有同時還有Y?是X?在集合Y上的投影，X?是Y?在集合X上的投影。另外，定理1可由文獻［20］中的收斂性結果得出。

定理1設S，T都是距離空間，G?S和H?T都是緊集。f是S×T上的連續實值函數。F是從G×H到G×H的集值映射，它在G×H上是閉的。如果下面兩點成立：（1）對G×H中任意兩點(u，v)和(x，y)，(u，v)∈F(x，y)，有f(u，v)≤f(x，y)；（2）如果(x，y)不是f在G×H上的全局極小點，?(u，v)∈F(x，y)，有f(u，v)

2 數值實驗

本節利用圖像恢復問題檢驗本文提出的優化設計方法是否能提高壓縮感知系統的重構精度，以及是否能提升測量矩陣的稀疏性?？偟膩碚f，按如下方法進行檢驗：用文獻［5，7，10，13，15］中的優化設計方法得到測量矩陣和稀疏字典；然后以其為初值，再用本文方法進行優化設計，得到新的測量矩陣（稀疏字典不變）；最后，用兩套測量矩陣和稀疏字典分別進行重構，比較重構效果。另外，用本文方法來處理標準壓縮感知系統的測量矩陣和稀疏字典，即稀疏字典是用KSVD算法［21］從訓練數據中學到的，而測量矩陣是高斯隨機矩陣（即每個元素獨立同分布地滿足正態分布N(0，M-1)，M為矩陣的行數）。下面分別闡述不同情形下的對比辦法和實驗結果。對所有情形，本文提出的方法采用參數T=1 500、Talter=15、β=0.5和σ=0.1。另外，所有代碼均采用MATLAB（R2019b）編寫，采用的筆記本電腦配置是Intel（R）Core（TM）i7?6500U處理器、主頻2.50 GHz、內存8 GB。

2.1 關于文獻［5］的數值實驗

首先，按如下方法獲取訓練數據：從LabelMe圖像數據集［22］中隨機選取400幅圖像，再從每幅圖像中無重疊地提取15個8×8的子塊，把它們放到一起得到6 000個64維訓練數據。再按如下方法獲取測試數據：取10幅經典的512像素×512像素自然圖像；把每幅無重疊地分割成一系列8×8的子塊，得到4 096個64維測試數據。其次，用文獻［5］中的優化設計方法從訓練數據中得到測量矩陣和稀疏字典，采 用 的 參 數 是稀 疏 度6、迭 代 步 數100、初 始 測 量 矩 陣randn(M，N)和N×L的過完備離散余弦變換（Discrete cosine transform，DCT）矩陣形式的初始稀疏字典。再用本文提出的方法對得到的測量矩陣進行優化設計。最后，用兩個測量矩陣分別對每一個測試數據進行測量，對理想測量值都添加高斯白噪聲形式的測量噪聲，使得信噪比等于10 dB。再用OMP算法和兩套等效字典從測量值中恢復測試數據。對原始圖像和重構圖像，用峰值信噪比（Peak signal to noise ratio，PSNR）和結構相似性指標（Structural similarity index，SSIM）［23］評判重構效果，采用MATLAB軟件內置命令PSNR和SSIM計算2個評判指標。另外，對測量矩陣，把絕對值小于10-4的分量認為是零。實驗結果如表1所示，采用文獻［5］方法得到的測量矩陣只包含0.520 8%的零元素，經本文方法處理后，這個比例升至57.552 1%。

表1 使用本文方法前后關于文獻［5］中等效字典的重構效果Table 1 Reconstruction performance for equivalent dictionary in Ref.［5］with and without the proposed method

2.2 關于文獻［7］的數值實驗

對采用文獻［7］方法生成的測量矩陣和稀疏字典，按照針對文獻［5］的比較辦法，比較了在使用本文方法前后等效字典帶來的重構效果和測量矩陣的稀疏性。文獻［7］中的優化設計方法使用參數M=20、N=64、L=100、α=0.4、β=0.8、稀疏度4、迭代步數100、初始測量矩陣randn(M，N)和N×L的過完備DCT矩陣形式的初始稀疏字典。實驗結果如表2所示。經本文方法處理后，采用文獻［7］方法得到的測量矩陣的零元素占比由0.156 3%升至78.828 1%。

表2 使用本文方法前后關于文獻［7］中等效字典的重構效果Table 2 Reconstr uction per for mance for equivalent dictionar y in Ref.［7］with and without the pr oposed method

2.3 關于文獻［10］的數值實驗

對采用文獻［10］方法生成的測量矩陣和稀疏字典，按照針對文獻［5］的比較辦法，比較在本文方法處理前后等效字典帶來的重構效果和測量矩陣的稀疏性。但采用的訓練數據不同。從LabelMe圖像數據集中選取2 920幅圖像，再從每幅圖像中無重疊地提取400個8×8的子塊，把它們放到一起得到1.168×106個64維訓練數據。另外，文獻［10］中的優化設計方法采用的參數是M=20、N=64、L=256、總的迭代步數10、字典學習的迭代步數、稀疏度4、ρ=15和隨機選取L個訓練數據整合在一起然后逐列單位化生成的初始稀疏字典。實驗結果如表3所示，采用文獻［10］方法得到的測量矩陣不包含零元素，而經本文方法處理后，這個比例升至84.921 9%。

表3 使用本文方法前后關于文獻［10］中等效字典的重構效果Table 3 Reconstr uction per formance for equivalent dictionar y in Ref.［10］with and without the pr oposed method

2.4 關于文獻［13］的數值實驗

因為文獻［13］的優化設計方法和本文一樣，在提高壓縮感知系統重構精度的同時也提升了測量矩陣的稀疏性。所以本文準備對同樣的測量矩陣Φ和稀疏字典Ψ，分別采用兩種方法進行優化設計，然后分別用兩套測量矩陣和稀疏字典進行重構，再比較重構效果和2個測量矩陣的稀疏性。采用如下辦法生成Φ和Ψ：按照針對文獻［5］的比較辦法生成訓練數據，然后用KSVD算法生成Ψ，采用的參數是M=20、N=64、L=100、稀疏度4、迭代步數50和N×L的過完備DCT矩陣形式的初始稀疏字典。再對KSVD算法的結果Ψ利用文獻［13］提供的程序ZRE.m生成Φ（程序ZRE.m詳見文獻［24］）。測試數據的生成以及重構過程則和針對文獻［5］的比較辦法中的一樣。另外，文獻［13］采用的參數是κ=20、λ=1.4和DCT矩陣形式的矩陣A。實驗結果如表4所示。經本文方法處理后，按文獻［13］中方法得到的測量矩陣的零元素占比由0.234 4%升至79.765 6%。

表4 使用本文方法前后關于文獻［13］中等效字典的重構效果Table 4 Reconstruction performance for equivalent dictionary in Ref.［13］with and without the proposed method

2.5 關于文獻［15］的數值實驗

對采用文獻［15］方法生成的測量矩陣和稀疏字典，按照針對文獻［5］的比較辦法進行比較。不同的是訓練數據、測試數據、添加測量噪聲的方式和評判標準不一樣。訓練數據按如下方式生成：取10幅經典的512像素×512像素自然圖像，滑動地把每幅圖分割成一系列10×10的子塊，把它們放在一起取其中的10%，得到25 301個100維訓練數據。測試數據由人工生成。利用KSVD算法從訓練數據中得到一個N×N維的稀疏字典Ψ，采用的參數是稀疏度4、迭代步數50和DCT正交矩陣形式的初始稀疏字典。測試數據x按如下方式生成：生成稀疏向量α（稀疏度是4，支撐集服從等概率分布，非零分量獨立同分布地服從(-1，1)上的均勻分布），計算x=Ψα?？偣采? 000個這樣的向量。對測量值逐分量地加上服從正態分布N(0，0.01)的噪聲。因為測試數據不是自然圖像，所以用全體測試數據及其重構結果的均方誤差來評判重構效果。用MATLAB軟件的內置命令immse計算均方誤差（Mean square error，MSE）。另外，文獻［15］中的優化設計方法采用的參數是：η=1.2、ρ=0.5、β=2、T=15、TAM=1000、M=20、N=L=100、初始測量矩陣和KSVD算法生成的初始稀疏字典Ψ。實驗結果如表5所示。對10幅圖像，按文獻［15］中方法得到的測量矩陣都不包含零元素，而經本文方法處理后，測量矩陣的零元素占比（對應表5中從上到下的順序）依次為76.700 0%、77.450 0%、80.250 0%、76.100 0%、78.250 0%、79.850 0%、78.500 0%、76.900 0%、72.700 0%和80.100 0%。

表5 使用本文方法前后關于文獻［15］中等效字典的重構效果Table 5 Reconstruction performance for equivalent dictionary in Ref.［15］with and without the proposed method

2.6 關于標準壓縮感知系統的數值實驗

按照針對文獻［5］的比較辦法進行比較，只是生成待比較的測量矩陣和稀疏字典的方式不同。用KSVD算法從訓練數據中得到稀疏字典，采用的參數是M=12、N=64、L=256、迭代步數50、稀疏度6和N×L維的過完備DCT矩陣形式的初始稀疏字典；而測量矩陣是高斯隨機矩陣。實驗結果如表6所示，生成的標準壓縮感知系統的測量矩陣不包含零元素，而經本文方法處理后，則包含29.166 7%的零元素。

表6 使用本文方法前后關于標準壓縮感知系統中等效字典的重構效果Table 6 Reconstruction performance for equivalent dictionary in standard compressed sensing system with and without the pr oposed method

2.7 重構效果的進一步展示

從2.4節實驗處理的10幅圖像中選擇4幅原始圖像。把每幅原始圖像、利用文獻［13］方法設計的測量矩陣和稀疏字典得到的重構圖像、利用已有稀疏字典和本文方法處理后的測量矩陣得到的重構圖像，以及上述3幅圖像的某個共同細節部分進行并列展示，具體如圖1～4所示。

圖1 圖像“lena”及其重構結果Fig.1 Image“lena”and its reconstruction

從本節所有數值實驗可以看到，對按照文獻［5，7，10，13，15］生成的測量矩陣和稀疏字典，以及標準壓縮感知系統的測量矩陣和稀疏字典，應用本文提出的方法進行優化得到的新測量矩陣搭配上已有稀疏字典能夠在圖像恢復階段帶來更高的信號重構精度，同時還能顯著提高測量矩陣中的零元素占比。

圖2 圖像“jetplane”及其重構結果Fig.2 Image“jetplane”and its reconstruction

圖3 圖像“boats”及其重構效果Fig.3 Image“boats”and its reconstruction

圖4 圖像“bridge”及其重構效果Fig.4 Image“bridge”and its reconstruction

3 結束語

本文利用最優化方法提出了一種壓縮感知等效字典的優化設計方法。提出的最優化問題在使等效字典的Gram矩陣逼近某類特定Gram矩陣的同時，要求測量矩陣滿足某種基于L1范數的不等式約束。然后，根據交替投影算法框架提出了收斂的求解算法。數值實驗表明，就圖像恢復問題而言，對采用一些經典的以及新提出的優化設計方法得到的測量矩陣和稀疏字典，以及標準壓縮感知系統的測量矩陣和稀疏字典，求解提出的最優化問題，得到的新測量矩陣搭配上已有稀疏字典能有效提高信號重構精度，同時顯著提升測量矩陣的稀疏性。

相較同類方法，本文提出的方法著重關注測量矩陣的稀疏性對后續信號重構效率的影響，并首次利用L1范數對測量矩陣賦予稀疏性。利用本文方法優化設計壓縮感知的等效字典，可在信號重構階段獲得精度與效率的雙提高。為了能夠利用本文方法優化設計大規模等效字典，下一步將會研究求解最優化問題（4）的更加高效的迭代算法，或者是最優化問題（4）的閉合形式解的求法。

附錄

現在取K=max{K1，K2}。則?k>K（k是正整數），有