基于奇異值分解的壓縮感知觀測矩陣優(yōu)化算法

李　周，崔　琛

(國防科技大學(xué) 電子對抗學(xué)院，合肥 230037)(*通信作者電子郵箱17756587331@163.com)

0　引言

壓縮感知(Compressive Sensing, CS)[1-3]在采樣率遠(yuǎn)小于奈奎斯特采樣率的條件下獲取信號的離散樣本，然后通過非線性優(yōu)化算法重構(gòu)出原始信號。重構(gòu)信號所需的采樣率并不取決于信號的帶寬，而與信號的稀疏度密切相關(guān)，因此CS有效降低了信號獲取、存儲及傳輸?shù)拇鷥r(jià)，引起了越來越多學(xué)者的關(guān)注。

信號的稀疏表示、觀測矩陣的構(gòu)造、重構(gòu)算法是CS理論中三個(gè)主要研究方向，其中觀測矩陣是影響壓縮感知性能的關(guān)鍵因素之一[4]。觀測矩陣構(gòu)造的目的是如何采樣得到觀測值，并能從觀測值中重構(gòu)出原始信號。文獻(xiàn)[5-7]對精確重構(gòu)所需觀測矩陣的約束條件進(jìn)行了研究：文獻(xiàn)[5]提出了零空間性質(zhì)(Null Space Property, NSP)，但在存在噪聲的情況下NSP并不能保證信號的精確重構(gòu)；文獻(xiàn)[6]提出了有限等距性質(zhì)(Restricted Isometry Property, RIP)，但判斷觀測矩陣是否符合RIP需要計(jì)算觀測矩陣n列中任意組合的K列，即需要計(jì)算觀測矩陣各列的組合，實(shí)際用于觀測矩陣的分析非常困難；文獻(xiàn)[7]提出了相關(guān)性的概念，其含義是觀測矩陣與稀疏基之間的相關(guān)程度，通過降低相關(guān)性可以減少CS的重構(gòu)誤差和精確重構(gòu)原始信號所需的觀測數(shù)。

文獻(xiàn)[7]通過線性收縮Gram矩陣中絕對值大于限定閾值的非對角元來降低相關(guān)性，取得了較好的實(shí)驗(yàn)效果。但是，該算法在收縮Gram矩陣時(shí)會改變矩陣的秩[8-9]；同時(shí)在求解觀測矩陣時(shí)需要求解稀疏基的逆矩陣，但由于稀疏基可能奇異，此時(shí)需要利用稀疏基的Moore-Penrose廣義偽逆來求解新的觀測矩陣，造成了算法中計(jì)算量大且精度受限的問題[4,8]。文獻(xiàn)[9]延續(xù)了文獻(xiàn)[7]之前的思路，不同之處在于增加了一個(gè)保留前一次矩陣優(yōu)化信息的操作，其目的在于使每一次矩陣更新的變化量減少，該算法繼承了文獻(xiàn)[7]的缺點(diǎn)，同時(shí)需要根據(jù)經(jīng)驗(yàn)手動(dòng)設(shè)置一個(gè)權(quán)重系數(shù)，用來衡量本次矩陣優(yōu)化結(jié)果與前一次矩陣優(yōu)化結(jié)果之間的比重。

文獻(xiàn)[10]提出將Gram矩陣非對角元素的平方和作為整體互相關(guān)系數(shù)，用來表示觀測矩陣的整體性能；在此基礎(chǔ)上，文獻(xiàn)[11]引入了α緊框架的概念，平均化Gram矩陣的非零特征值，減小了整體相關(guān)系數(shù)。但在由Gram矩陣求解觀測矩陣時(shí)，文獻(xiàn)[10-11]仍采用了求稀疏基Moore-Penrose廣義偽逆矩陣的做法，同樣造成了算法中計(jì)算量大且精度受限的問題。

針對現(xiàn)有算法從優(yōu)化后的Gram矩陣求解觀測矩陣時(shí)出現(xiàn)的相關(guān)系數(shù)較大與需要求廣義偽逆矩陣的問題，本文借鑒了K-SVD(K-Singular Value Decomposition)算法中逐行更新優(yōu)化目標(biāo)矩陣的思想，在利用現(xiàn)有算法得到優(yōu)化后的Gram矩陣的基礎(chǔ)上，提出了一種基于奇異值分解(Singular Value Decomposition, SVD)的觀測矩陣優(yōu)化算法：該算法通過求解等價(jià)變換后的目標(biāo)函數(shù)對觀測矩陣行向量的導(dǎo)數(shù)得到目標(biāo)函數(shù)取極值時(shí)行向量的值，并通過對誤差矩陣進(jìn)行奇異值分解在上述行向量的值中選出使得目標(biāo)函數(shù)取最值時(shí)行向量的解析式，之后在每輪迭代中對觀測矩陣的每一行分別使用上述的行向量解析式進(jìn)行優(yōu)化。

1　觀測矩陣相關(guān)性和相應(yīng)算法

1.1　觀測矩陣相關(guān)性

假定離散信號r∈Rl×h，若r中最多含有K個(gè)非零值且K≤h，那么r就稱作K-稀疏的。將K-稀疏的信號集合記為：

ΛK={r:‖r‖0≤K}

(1)

若r本身為非稀疏的，但可以經(jīng)過一個(gè)稀疏基Ψ作如下變換：

r=Ψs

(2)

如果變換后的s符合‖s‖0≤K≤h，即r可以用已知的稀疏基中少量的元素線性表示，那么r也被稱作K-稀疏的。

CS理論的測量過程[1]可以表示為：

y=Φr=ΦΨs=Xs

(3)

其中：r為原始信號，Ψ∈Rn×h為稀疏基，s為變換后的稀疏信號，Φ∈Rm×n為觀測矩陣，XΦΨ稱為感知矩陣。

對于任意兩個(gè)不同的稀疏信號r1,r2∈ΛK，必須滿足Xr1≠Xr2，否則僅僅根據(jù)觀測值無法區(qū)分r1,r2，因此感知矩陣需要滿足一定的條件：對于任意K-稀疏的信號，只要其稀疏度K滿足下式，信號即可準(zhǔn)確地恢復(fù)出來[12]。

(4)

其中μ(X)指X任意兩列xi和xj之間的內(nèi)積絕對值的最大值[7]，計(jì)算公式如下：

(5)

然而μ(X)僅表示觀測矩陣與稀疏基之間列向量最大的相關(guān)性，是一種局部的相關(guān)性，因?yàn)樵赬只有某兩列的相關(guān)性比較大而其余各列之間相關(guān)性比較小的情況下，就會出現(xiàn)相關(guān)系數(shù)μ(X)很大而感知矩陣X的性能并不差的情況。

因此文獻(xiàn)[7]提出了表示總體相關(guān)性的t-平均相關(guān)性μt-av，定義為X所有列向量兩兩之間的內(nèi)積絕對值中大于限定閾值t的部分的平均值，或者是X所有列向量兩兩之間的內(nèi)積絕對值中最大的t%部分的均值，其定義式如下：

(6)

其中:gij為Gram矩陣G=XTX中的元素，gij為矩陣X的第i列xi與第j列xj的內(nèi)積，Nav為Hav中的元素?cái)?shù)，即Gram矩陣非對角元|gij|大于t的數(shù)目或者所有|gij|中最大的t%部分的元素?cái)?shù)目。僅闡述μt-av為X任意兩列內(nèi)積絕對值大于t的平均值時(shí)Hav的定義：

Hav={(i,j):gij>t,i≠j}

(7)

1.2　觀測矩陣優(yōu)化的相關(guān)算法及問題分析

為降低X任意兩列的相關(guān)性，即減少Gram矩陣非對角元素的值，文獻(xiàn)[7]采用如下閾值函數(shù)對Gram矩陣G中絕對值大于限定閾值的非對角元進(jìn)行線性收縮：

(8)

其中：t為閾值，γ為收縮因子。

(9)

其中：A=ΦΨ。在Ψ逆矩陣存在時(shí)，Φ=AΨ-1；但在稀疏基Ψ可能奇異導(dǎo)致其逆矩陣不存在時(shí)，需要利用稀疏基的Moore-Penrose廣義偽逆來求解觀測矩陣Φ=AΨ+，從而帶來計(jì)算量與精度的問題[4,8]。

文獻(xiàn)[6]的閾值函數(shù)只能使局部比較大的非對角元素減小。為整體減小Gram矩陣中的非對角元素，文獻(xiàn)[11]通過平均化Gram矩陣的非零特征值，使得矩陣非零特征值的平方和減小，降低了整體相關(guān)系數(shù)。但在由Gram矩陣求解觀測矩陣時(shí)，仍然需要求解稀疏基的廣義逆矩陣，故存在著與文獻(xiàn)[7]相同的問題[4,8]。

1.3　問題的解決思路

在現(xiàn)有算法優(yōu)化Gram矩陣后，本文借鑒了文獻(xiàn)[13]提出的K-SVD算法中逐行優(yōu)化目標(biāo)矩陣的思想從優(yōu)化后的Gram矩陣求解觀測矩陣。

K-SVD是一種通過逐行優(yōu)化字典矩陣對信號進(jìn)行稀疏表示的方法，具體目標(biāo)為：

(10)

其中:Y為要表示的信號，D為所求的字典，X為稀疏矩陣。X與Y按列對應(yīng)，表示D中的元素以xi為系數(shù)線性組合就可得到Y(jié)，而K-SVD的目的是在X和Y已知的情況下更新字典D滿足上述條件。

(11)

式(11)可以看作把第k個(gè)分量剝離后表達(dá)式會產(chǎn)生空洞，目的是找到一個(gè)新向量以更好地填補(bǔ)這個(gè)洞：假設(shè)除第k項(xiàng)外其余項(xiàng)是固定的，之后對Ek進(jìn)行奇異值分解得到Ek=UΛVT，其中U和V的列矢量均是正交基，Λ是對角矩陣。若Λ的對角元素從大到小排列，則表示Ek的能量分量主軸在相應(yīng)幾個(gè)正交方向上由大到小分配，取U的第一個(gè)列向量來表示di，取V的第一個(gè)列向量與Λ的第一個(gè)元素的乘積表示xi，至此完成了字典一個(gè)條目的更新。

2　觀測矩陣的優(yōu)化算法

在利用現(xiàn)有算法得到優(yōu)化后的Gram矩陣的基礎(chǔ)上，借鑒K-SVD算法中逐行更新優(yōu)化目標(biāo)矩陣的思想，本章利用稀疏基的奇異值分解對目標(biāo)函數(shù)進(jìn)行等價(jià)變換，通過求解等價(jià)變換后的目標(biāo)函數(shù)對觀測矩陣行向量的導(dǎo)數(shù)得到目標(biāo)函數(shù)取極值時(shí)行向量的值，并通過對誤差矩陣進(jìn)行奇異值分解在上述行向量的值中選出使得目標(biāo)函數(shù)取最值時(shí)行向量的解析式，最后利用所求出的行向量解析式逐行迭代優(yōu)化觀測矩陣。

2.1　觀測矩陣行向量的優(yōu)化

(12)

假設(shè)稀疏基Ψ的秩為NΨ，其奇異值分解為：

(13)

將奇異值分解式代入目標(biāo)函數(shù)中，可得：

(14)

根據(jù)F范數(shù)的酉不變性，在式(14)中左乘VΨT，右乘VΨ，同時(shí)設(shè)ΦUΨ其中為的前NΨ列。設(shè)其中為的前NΨ行與前NΨ列，將和代入式(14)可得：

(15)

由式(15)易得：

(16)

(17)

(18)

式(18)對ωj求導(dǎo)可得：

(19)

導(dǎo)數(shù)置0便可得到一系列極值點(diǎn)：

(20)

易得：

(21)

(22)

2.2　觀測矩陣的迭代優(yōu)化算法

在現(xiàn)有算法得到優(yōu)化后的Gram矩陣的基礎(chǔ)上，本小節(jié)利用2.1節(jié)求出的觀測矩陣行向量的解析式，采用逐行更新的方法優(yōu)化觀測矩陣。

輸入初始觀測矩陣Φ∈Rm×n，離散傅里葉變換基Ψ∈Rn×h，迭代次數(shù)Iter，退出閾值μ0;

輸出觀測矩陣。

fork=1:Iter

1)

2)

forj=1:m

end for

3)

4)

計(jì)算，若μt-av與上一輪迭代的μt-av之差小于μ0則退出循環(huán)，否則轉(zhuǎn)至1)繼續(xù)循環(huán)

end for

3　實(shí)驗(yàn)結(jié)果及其分析

本文在Matlab平臺上對算法進(jìn)行仿真實(shí)驗(yàn)，仿真選擇高斯隨機(jī)矩陣作為初始觀測矩陣，離散傅里葉變換基作為稀疏基，比較文獻(xiàn)[7]所提的算法CSElad、文獻(xiàn)[11]所提的算法CSTsiligianni和經(jīng)文獻(xiàn)[11]中的方法優(yōu)化Gram矩陣后使用本文所提的算法CSImproved-Tsiligianni在相關(guān)性、重構(gòu)誤差和運(yùn)行時(shí)間三方面的性能。

仿真中相關(guān)參數(shù)如下：m=30，n=l=100，稀疏度K=10，迭代次數(shù)Iter=100。Gram矩陣收縮變換時(shí)非對角線的限定閾值t=30%，收縮因子γ=0.5，CSImproved-Tsiligianni的退出閾值μ0=10-2。為減少隨機(jī)性對實(shí)驗(yàn)結(jié)果造成的影響，每次仿真均為1 000次的蒙特卡羅仿真。

3.1　相關(guān)性

仿真實(shí)驗(yàn)一 ：CSElad、CSTsiligianni和CSImproved-Tsiligianni都是迭代進(jìn)行的算法，去除改進(jìn)算法第四步迭代退出的步驟，將三種算法每次迭代時(shí)Gram矩陣的t-平均相關(guān)性μt-av進(jìn)行對比，可以看出在各個(gè)算法迭代時(shí)t-平均相關(guān)性的變化趨勢，進(jìn)而對比出三種算法在t-平均相關(guān)性這一方面的優(yōu)劣性。t-平均相關(guān)性μt-av中的參數(shù)t=20%，其含義是Gram矩陣非對角元中最大的20%部分的均值。

圖1　不同算法的μt-av隨迭代的變化Fig. 1　Changes in μt-av of different algorithms with iterations

由圖1可知，隨著迭代的進(jìn)行，CSImproved-Tsiligianni的μt-av逐漸減小而趨于穩(wěn)定，而CSElad和CSTsiligianni的μt-av變化范圍較大;同時(shí)，改進(jìn)算法的μt-av比改進(jìn)之前的算法小，這是由于改進(jìn)算法利用優(yōu)化后的觀測矩陣行向量的解析式逐行優(yōu)化觀測矩陣。

3.2　重構(gòu)誤差

(23)

仿真實(shí)驗(yàn)二：觀測次數(shù)m對壓縮感知的重構(gòu)誤差的影響。仿真了觀測次數(shù)m從25增加到45時(shí)三種算法重構(gòu)誤差的變化情況，如圖2所示。可以看出，在給定信號長度和稀疏度的前提下，m越大，重構(gòu)誤差越?。籱越小，重構(gòu)誤差越大。

圖2　觀測次數(shù)變化時(shí)三種算法的MSE變化情況Fig. 2　Changes in MSE of three algorithms with observation number

仿真實(shí)驗(yàn)三：稀疏度K對算法的重構(gòu)性能的影響。為了得出稀疏度K對各個(gè)算法重構(gòu)性能的影響，仿真稀疏度K由10變化到20時(shí)三種算法的MSE的變化情況，如圖3所示?？梢钥闯觯诮o定觀測次數(shù)和信號長度的前提下，稀疏度越高，重構(gòu)誤差越大。

圖3　稀疏度變化時(shí)三種算法的MSE變化情況Fig. 3　Changes in MSE of three algorithms with signal’s sparsity

由仿真實(shí)驗(yàn)二與三可知在K或m變化時(shí)CSImproved-Tsiligianni的MSE小于CSElad和CSTsiligianni，這是由于CSImproved-Tsiligianni產(chǎn)生的觀測矩陣與稀疏基的相關(guān)性小于CSElad和CSTsiligianni，可以更好地保持不同K稀疏向量之間的距離。

3.3　運(yùn)行時(shí)間

為了研究本文所提算法的復(fù)雜度，該仿真實(shí)驗(yàn)統(tǒng)計(jì)CSImproved-Tsiligianni優(yōu)化觀測矩陣所需的運(yùn)行時(shí)間，并與CSElad和CSTsiligianni的運(yùn)行時(shí)間進(jìn)行比較。雖然算法的運(yùn)行時(shí)間不能嚴(yán)格地定義算法復(fù)雜度，但仍可以在一定程度上對算法的復(fù)雜度作出描述。仿真環(huán)境為2.53 GHz英特爾i5四核處理器、4 GB內(nèi)存Windows 10系統(tǒng)下的Matlab R2014a。圖4給出了各算法在信號長度l和觀測次數(shù)m變化時(shí)相應(yīng)的運(yùn)行時(shí)間。

圖4　三種算法的運(yùn)行時(shí)間對比Fig. 4　Running time comparison of three algorithms

由圖4可知信號長度和觀測次數(shù)增加時(shí)算法的運(yùn)行時(shí)間隨之增加，其中CSElad和CSTsiligianni的運(yùn)行時(shí)間小于CSImproved-Tsiligianni的運(yùn)行時(shí)間，這是因?yàn)镃SImproved-Tsiligianni在優(yōu)化觀測矩陣時(shí)對每一行均采用奇異值分解來求解優(yōu)化后的行向量，而奇異值分解耗時(shí)是較長的。不過相對于在相關(guān)性與重構(gòu)誤差方面的優(yōu)勢而言，CSImproved-Tsiligianni增加的運(yùn)行時(shí)間還是可以接受的。

4　結(jié)語

本文對壓縮感知中觀測矩陣的優(yōu)化問題進(jìn)行研究，借鑒K-SVD算法中逐行更新目標(biāo)矩陣的思想，在現(xiàn)有算法得到優(yōu)化后的Gram矩陣基礎(chǔ)上，提出了一種基于奇異值分解的觀測矩陣優(yōu)化算法。仿真結(jié)果表明：在可接受的運(yùn)算量下，本文所提算法在觀測矩陣與稀疏基的相關(guān)性方面優(yōu)于改進(jìn)前的算法，從而提高了重構(gòu)精度。如何優(yōu)化閾值函數(shù)來構(gòu)造性能更優(yōu)的Gram矩陣是下一步的研究方向。

參考文獻(xiàn):

[1]CANDES E J, ROMBERG J, TAO T. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information [J]. IEEE Transactions on Information Theory, 2006, 52(2): 489-509.

[2]DONOHO D L. Compressed sensing [J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289-1306.

[3]CANDES E J, TAO T. Near-optimal signal recovery from random projections: universal encoding strategies? [J]. IEEE Transactions on Information Theory, 2006, 52(12): 5406-5425.

[4]王強(qiáng),張培林,王懷光,等.壓縮感知中測量矩陣構(gòu)造綜述[J].計(jì)算機(jī)應(yīng)用,2017,37(1):188-196. (WANG Q, ZHANG P L, WANG H G, et al. A survey on the construction of measurement matrices in compressive sensing [J]. Journal of Computer Applications, 2017, 37(1): 188-196.)

[5]COHEN A, DAHMEN W, DEVORE R. Compressed sensing and best k-term approximation [J]. Journal of the American Mathematical Society, 2009, 22(1): 211-231.

[6]CANDES E J, TAO T. Decoding by linear programming [J]. IEEE Transactions on Information Theory, 2005, 51(12): 4203-4215.

[7]ELAD M. Optimized projections for compressed sensing [J]. IEEE Transactions on Signal Processing, 2007, 55(12): 5695-5702.

[8]鄭紅,李振.壓縮感知理論投影矩陣優(yōu)化方法綜述[J].數(shù)據(jù)采集與處理,2014,29(1):43-53. (ZHENG H, LI Z. Survey on optimization methods for projection matrix in compress sensing theory [J]. Journa1 of Data Acquisition and Processing, 2014, 29(1): 43-53.)

[9]XU J, PI Y, CAO Z. Optimized projection matrix for compressive sensing [J]. EURASIP Journal on Advances in Signal Processing, 2012, 2010: Article No. 43.

[10]趙瑞珍,秦周,胡紹海.一種基于特征值分解的測量矩陣優(yōu)化方法[J].信號處理,2012,28(5):653-658. (ZHAO R Z, QIN Z, HU S H. An optimization method for measurement matrix based on eigenvalue decomposition [J]. Signal Processing, 2012, 28(5): 653-658.)

[11]TSILIGIANNI E, KONDI L P, KATSAGGELOS A K. Use of tight frames for optimized compressed sensing [C]// Proceedings of the 20th European Signal Processing Conference (EUSIPCO). Piscataway, NJ: IEEE, 2012: 1439-1443.

[12]DONOHO D L, ELAD M. Optimally sparse representation in general (nonorthogonal) dictionaries via1minimization [J]. Proceedings of the National Academy of Sciences of the United States of America, 2003, 100(5): 2197-2202.

[13]AHARON M, ELAD M, BRUCKSTEIN A.rmK-SVD: an algorithm for designing overcomplete dictionaries for sparse representation [J]. IEEE Transactions on Signal Processing, 2006, 54(11): 4311-4322.

[14]ABOLGHASEMI V, FERDOWSI S, SANEI S. A gradient-based alternating minimization approach for optimization of the measurement matrix in compressive sensing [J]. Signal Processing, 2012, 92(4): 999-1009.

[15]KWON S, WANG J, SHIM B. Multipath matching pursuit [J]. IEEE Transactions on Information Theory, 2014, 60(5): 2986-3001.