QR分解與特征值優(yōu)化觀測矩陣的算法研究

鄭曉，薄華，孫強(qiáng)

(1.上海海事大學(xué) 信息工程學(xué)院，上海201306;2.西安理工大學(xué)自動化與信息工程學(xué)院，陜西西安710000)

壓縮感知理論的突出優(yōu)點(diǎn)在于能夠同時完成信號獲取和壓縮，它突破了奈奎斯特采樣定理，能夠以遠(yuǎn)少于傳統(tǒng)方法所需的采樣數(shù)來進(jìn)行圖像的精確重構(gòu)。它在減少冗余數(shù)據(jù)、節(jié)省存儲空間上有較大優(yōu)勢。因而，壓縮感知掀起了信號處理領(lǐng)域的革命，被廣泛應(yīng)用于醫(yī)用電子學(xué)、模式識別、數(shù)據(jù)采樣與挖掘、無線通信、信道編碼、天文學(xué)以及雷達(dá)遙感等領(lǐng)域。觀測矩陣的設(shè)計是壓縮感知研究的關(guān)鍵問題之一，構(gòu)造性能好的觀測矩陣對于觀測值的獲取和圖像的精確重構(gòu)都起到了至關(guān)重要的作用。目前，國內(nèi)外提出的優(yōu)化方法主要有:通過減小格拉姆矩陣(Gram)的非對角線元素增大觀測矩陣與稀疏基之間的非相干性，從而實現(xiàn)觀測矩陣的優(yōu)化。用到的方法主要有等角緊框架法、梯度迭代法以及特征值優(yōu)化方法[1-11]等;利用增強(qiáng)觀測矩陣的列獨(dú)立性實現(xiàn)觀測矩陣優(yōu)化[12-13]等等。上面提到的優(yōu)化方法主要是通過減小觀測矩陣與稀疏基之間的互相干性單方面展開[1-11]，或者通過增大觀測矩陣的列獨(dú)立性展開[12-13]。文中將把增強(qiáng)觀測矩陣的列獨(dú)立性和增大觀測矩陣與稀疏基之間的非相干性結(jié)合起來，尋求一個更加優(yōu)化的觀測矩陣，從而利用相同的觀測個數(shù)，實現(xiàn)重構(gòu)質(zhì)量更好、仿真實驗更穩(wěn)定的效果。

1 壓縮感知

壓縮感知(compressed sensing，CS)突破了香農(nóng)采樣定理的瓶頸，使高分辨率信號的采集成為可能，引起廣大學(xué)者的研究。對壓縮感知進(jìn)行數(shù)學(xué)建模，它主要涉及3方面的內(nèi)容[14]:

1)尋求稀疏基Ψ∈RN×N，并將信號X∈RN進(jìn)行稀疏表示，由式(1)可以得到稀疏向量Θ∈RN:

2)設(shè)計一個平穩(wěn)的、與稀疏基ψ不相關(guān)的觀測矩陣Φ∈RM×N，保證稀疏向量Θ能夠從N維降到M維，同時不丟失重要信息，進(jìn)而得到觀測集合Y∈RM，如式(2)所示:

3)設(shè)計快速重構(gòu)算法，利用觀測集合Y并根據(jù)式(3)恢復(fù)原信號:

稀疏表示是指在某個特定變換域中用盡量少的基函數(shù)來盡可能完整地表示原始信號[15]。信號能夠進(jìn)行有效的稀疏表示是壓縮感知的前提。所以，稀疏基的選擇是壓縮感知的先決條件。常用的稀疏基有DCT基、小波基、Chirplet基、Gabor基等，但是這些固定的正交基有時還不足以表示如聲音或自然圖像這些信號所具有的復(fù)雜未知規(guī)則性，不能保證信號在變換域足夠稀疏。稀疏基的構(gòu)造不是文中的研究重點(diǎn)，根據(jù)文中實驗圖像的特點(diǎn)，文中采用較常用的小波基。

觀測矩陣進(jìn)行的線性觀測是把信號稀疏表示和信號重構(gòu)連接起來的橋梁，觀測矩陣的設(shè)計是壓縮感知的重要部分[15]。提出了觀測矩陣設(shè)計時應(yīng)該滿足的3個原則，根據(jù)這3個原則構(gòu)造了一些常用觀測矩陣，例如高斯矩陣、正交觀測矩陣、多項式確定觀測矩陣、循環(huán)矩陣、循環(huán)直積觀測矩陣等。文中采用高斯矩陣作為基本矩陣。

常用的重構(gòu)算法主要分為3類[16]:1)貪婪追蹤算法:匹配追蹤算法、正交匹配追蹤算法(OMP)、分段正交匹配追蹤算法等;2)凸松弛算法:梯度投影法、基追蹤算法內(nèi)點(diǎn)法等;3)組合算法:傅立葉采樣以及鏈?zhǔn)阶粉櫟取Ｎ闹胁捎肙MP算法。

2 文中優(yōu)化算法

現(xiàn)有的觀測矩陣幾乎都存在著一些不足，目前觀測矩陣的研究熱點(diǎn)在于:尋找需要更少觀測值的新矩陣;對觀測矩陣進(jìn)行優(yōu)化，在相同觀測數(shù)的情況下，使觀測矩陣具有更好的性質(zhì);構(gòu)造或優(yōu)化能夠減小計算復(fù)雜度和提高實驗穩(wěn)定性的觀測矩陣。

為了利用更少的觀測值得到更加精確、更加穩(wěn)定的重構(gòu)結(jié)果，文中結(jié)合觀測矩陣列向量的獨(dú)立性以及觀測矩陣與稀疏基之間的非相干性兩方面對觀測矩陣進(jìn)行優(yōu)化，提出一種利用QR分解來增大觀測矩陣的列獨(dú)立性方法，同時提出通過優(yōu)化Gram矩陣特征值增大觀測矩陣與稀疏基之間的非相干性的優(yōu)化方法(以下簡稱QT算法)。

2.1 QR分解增大列獨(dú)立性

觀測矩陣的列獨(dú)立性越大，重構(gòu)所需的觀測值越少，重構(gòu)質(zhì)量也越高。文獻(xiàn)[19]中指出設(shè)計觀測矩陣時，要保證觀測矩陣的最小奇異值大于某一個大于零的常數(shù)。而文獻(xiàn)[20]中指出矩陣的最小奇異值與矩陣的列向量相關(guān)性密切相關(guān)。最小奇異值越大矩陣的列相關(guān)性越弱，列獨(dú)立性越強(qiáng)。因此，觀測矩陣的最小奇異值是影響圖像重構(gòu)質(zhì)量的重要因素。所以，可以在不違背文獻(xiàn)[15]中提出的觀測矩陣應(yīng)該滿足的性質(zhì)前提下，通過增大觀測矩陣的最小奇異值對觀測矩陣進(jìn)行優(yōu)化，最終會得到更好的重構(gòu)效果。

QR分解優(yōu)化不但能夠增大矩陣的最小奇異值，同時能夠保持矩陣的性質(zhì)基本不變[16]。QR分解優(yōu)化主要原理為:首先對觀測矩陣Φ∈RM×N進(jìn)行QR分解，得到正交矩陣Q∈RM×N和上三角矩陣A∈RN×N。然后對A分析，看出其主對角線元素遠(yuǎn)遠(yuǎn)大于非對角線元素，利用該發(fā)現(xiàn)，把A的非對角線元素全置零得A1，再根據(jù)A1得優(yōu)化后的觀測矩陣Φ1∈RM×N。通過QR分解優(yōu)化，Φ的最小奇異值小于 Φ1的最小奇異值[12]。

2.2 基于特征值的優(yōu)化

文獻(xiàn)[1-10]提出觀測矩陣與稀疏基之間的互相干性越小，重構(gòu)所需的觀測值越少，重構(gòu)精度也越大。觀測矩陣與稀疏基之間的互相干性μ(Φ，Ψ)如式(4)所示:

基于特征的優(yōu)化算法是針對Gram矩陣的特征值進(jìn)行的，算法的本質(zhì)在于通過減小Gram矩陣的非對角線元素減小互相干性。

首先，構(gòu)造 D∈RN×N，令 D=ΦΨ，其中觀測矩陣Φ∈RM×N，并且秩為 M，稀疏基 Ψ∈RN×N。對 D 進(jìn)行歸一化得到 D1，然后構(gòu)造 Gram矩陣 G∈RN×N，令G=DT1D。G的非對角線元素的大小和觀測矩陣與稀疏基之間的非相干性的大小有著密切關(guān)系:非對角線元素越小，互相干性越小。文中通過使G的非對角線元素的平方和最小化實現(xiàn)互相干性最小化。

對上面的G進(jìn)行分析，假設(shè)G有M個大于0的特征值lk(k＜M)，根據(jù)矩陣的跡理論以及特征值的相關(guān)理論，加上G是歸一化矩陣，求解觀測矩陣和稀疏基之間最小相關(guān)性問題可以等價為求解式(5)的最優(yōu)解。

式中:gij表示G中的元素，di、dj均為G的列向量。通過式(5)可以看出，當(dāng)lk都為N/M時可得到理論的最優(yōu)解。此時，式(5)等價為式(6)

文中優(yōu)化算法的主體思路是選擇高斯隨機(jī)矩陣為原始觀測矩陣，利用它構(gòu)造Gram矩陣，之后對Gram矩陣優(yōu)化，讓它的特征值逐漸逼近N/M。然后通過優(yōu)化后的Gram矩陣反向求出此時的觀測矩陣，并且對此時的觀測矩陣進(jìn)行QR分解優(yōu)化，增大觀測矩陣的列獨(dú)立性。繼續(xù)用求得的觀測矩陣求解Gram矩陣，在不斷地進(jìn)行觀測矩陣和Gram矩陣的相互求解過程中，當(dāng)達(dá)到設(shè)置的終止條件時，輸出此時的觀測矩陣。通過上述優(yōu)化后，觀測矩陣不但列獨(dú)立性得到了增強(qiáng)，而且它與稀疏基的非相干性也相應(yīng)增大了。另外，由于每次迭代中QR分解優(yōu)化時A只保留相對較大的主對角線元素，觀測矩陣保存了主要信息，又因為G不斷趨向于最優(yōu)化矩陣，加上觀測矩陣與Gram矩陣的不斷相互求解，因此，每次得到的觀測矩陣很相近，相應(yīng)地，仿真結(jié)果穩(wěn)定性也較好。

2.3 算法流程及具體實現(xiàn)

通過算法分析，設(shè)計出QT算法流程如圖1。

圖1 QT算法流程圖Fig.1 The flow of QT algorithm

對上述流程進(jìn)行細(xì)化，具體實現(xiàn)如下:

輸入 Φ∈RN×NΨ∈RN×N，最大迭代次數(shù)設(shè)為100。

輸出 優(yōu)化后的觀測矩陣Φ2。

具體實現(xiàn)步驟如下:

1)構(gòu)造D∈RM×N，并且令 D=ΦΨ，對其進(jìn)行歸一化處理，得到歸一化矩陣D1;

2)由G=DT1D1構(gòu)造出Gram矩陣G，然后對G進(jìn)行奇異值分解，即G=USVT，其中U，V均為酉矩陣，S 為 對 角 矩 陣，并 且 U ∈ RN×N，V ∈ RN×N，S∈RN×N;

3)將S中對角線上的非零元素置為N/M，得到S1，接著得到G1=US1VT;

4)根據(jù)S1=LTL，對S1進(jìn)行分解，得到矩陣L。其中L∈RM×N，并且它的對角線元素為，其他元素為零;

5)設(shè) D2∈RM×N，并且令 D2=LVT，然后通過Φ1=D2Ψ－1，求得 Φ1;

6)對Φ1進(jìn)行近似QR分解的優(yōu)化:即首先對Φ1進(jìn)行QR分解，得到Φ1=EF。其中，正交陣E∈RM×N，而上三角矩陣F∈RM×N;然后將F中的非對角線元素設(shè)置為零，得到F1;最后根據(jù)Φ2=EF1，求得進(jìn)一步更新的觀測矩陣Φ2;

7)利用D2=Φ2Ψ求出D2，接著對D2進(jìn)行歸一化處理，得到D3，然后繼續(xù)利用G2=DT3D3得到新的Gram矩陣G2，求出G2中除主對角線元素外的所有值的平方和(sum)，當(dāng)(sum)－((N/M)2－N)|＜0.1時，輸出此時的觀測矩陣Φ2，否則跳轉(zhuǎn)到②繼續(xù)新的循環(huán)。

3 實驗結(jié)果與分析

為了驗證文中算法的有效性，文中選擇多幅自然圖像，應(yīng)用文中提出的算法進(jìn)行觀測矩陣的構(gòu)造，然后對觀測數(shù)據(jù)，采用OMP算法進(jìn)行重構(gòu)，將重構(gòu)圖像的信噪比和目前的經(jīng)典算法進(jìn)行比較。

采用256×256的自然圖像，初始的觀測矩陣采用高斯矩陣，稀疏基選取小波基，重構(gòu)算法選取OMP算法。進(jìn)行對比算法的觀測矩陣分別是初始高斯矩陣，經(jīng)QR分解優(yōu)化的高斯矩陣[12]，基于特征值優(yōu)化的高斯矩陣[11]。實驗部分分為2個部分:對比不同的觀測長度下，不同觀測矩陣得到的重構(gòu)精度的比較(多次重構(gòu)的平均結(jié)果);比較在相同觀測長度下，不同觀測矩陣得到的重構(gòu)精度(單次重構(gòu)的結(jié)果)。

3.1 不同觀測長度的比較

為了驗證文中構(gòu)造的觀測矩陣的有效性，對于同樣一幅自然圖像，分別選取不同的觀測矩陣長度M，比較采用不同觀測矩陣得到的觀測向量的重構(gòu)精度，并用信噪比表示。由于觀測矩陣選取的隨機(jī)性，所以每次結(jié)果都不盡相同。

為得到一個統(tǒng)計性的結(jié)果，采用50次平均的結(jié)果為最終結(jié)果。表1、2分別表示4種方法在不同觀測個數(shù)時lena和cameralman重構(gòu)圖像的信噪比。

表1 不同觀測個數(shù)的PSNR(lena)Table 1PNSR ofdifferent measurement numbers(lena)dB

表2 不同觀測個數(shù)的PSNR(cameraman)Table 2 The PNSR ofdifferent measurement dB

由表1和表2中的數(shù)據(jù)可以看出，經(jīng)QR分解優(yōu)化后，PSNR相對于未經(jīng)優(yōu)化時提高的幅度不明顯，而經(jīng)過特征值優(yōu)化后，PSNR相對于未經(jīng)優(yōu)化時提高的幅度較大。將QT算法和前3種算法中效果最好的基于特征值優(yōu)化算法相比較:當(dāng)M=80、100、120、140、160 時，lena 的 PSNR 提高情況分別為7.6、3.8、3.5、2.5 和 1.8 dB，而 cameraman 的 PSNR 提高情況分別為 6.5 、3.1 、2.9 、2.2 與 1.2 dB。通過比較分析，隨著M的增大，文中方法的優(yōu)勢越來越小，也就是說M越小，優(yōu)化算法優(yōu)勢更大?？偟膩碚f，文中優(yōu)化方法在保證近似精確重構(gòu)的條件下，減小觀測次數(shù)上更有空間。

3.2 算法的魯棒性比較

由于觀測矩陣的初始矩陣均為高斯噪聲，因此每次仿真試驗的初始矩陣都是隨機(jī)的，從而每次重構(gòu)的結(jié)果都會有所差異。由于QT算法同時考慮了觀測矩陣列的獨(dú)立性，以及觀測矩陣和基函數(shù)的不相關(guān)性，因此，比其他3種算法有更好的魯棒性。依然選取lena和cameraman 2幅圖像進(jìn)行圖示對比。圖2和圖3分別為當(dāng)M=80時對2幅圖像進(jìn)行觀測重構(gòu)之后的單次重構(gòu)精度的曲線圖。

圖2 M=80時的PNSR圖(lena)Fig.2 The graph of PNSR when the value of M is 80(lena)

圖3 M=80時的PNSR圖(cameraman)Fig.3 The graph of PNSR when the value of M is 80(cameraman)

由圖2、3可以看出:當(dāng)M=80時，其他3種方法的PSNR變化幅度比較大，而且經(jīng)QR分解優(yōu)化和特征值優(yōu)化后的PSNR值在個別點(diǎn)上比未經(jīng)優(yōu)化時小。文中優(yōu)化方法每次的PSNR始終大于其他3種方法，并且每次它的波動在5%以內(nèi)，穩(wěn)定性較好。另外，隨著觀測值的增多，其他方法的實驗仿真穩(wěn)定性也會有所提高，和文中優(yōu)化算法在穩(wěn)定性上的趨勢比較接近。由上述實驗可以得到結(jié)論:當(dāng)觀測值較少時，文中提出的優(yōu)化算法無論在PSNR，還是實驗結(jié)果的穩(wěn)定性上，都比其他3類算法有更好的效果。圖4和圖5給出了lena和cameraman在M=120時的4種算法的重構(gòu)結(jié)果。

圖4 不同方法優(yōu)化的重構(gòu)圖(lena)Fig.4 The reconstructed images by different optimization methods(lena)

圖5 不同方法優(yōu)化的重構(gòu)圖(cameraman)Fig.5 The reconstructed images by different optimization methods(cameraman)

從對圖4和圖5進(jìn)行分析可以看出，M=120時，經(jīng)QR分解優(yōu)化效果不大?；谔卣髦祪?yōu)化的重構(gòu)圖相對前2種方法在帽子邊緣和頭發(fā)、臉部這些細(xì)節(jié)部分明顯變得更清晰，重構(gòu)質(zhì)量得到了較大提高。文中算法優(yōu)化的重構(gòu)圖在鼻子、嘴巴以及眼睛這3個細(xì)節(jié)處的重構(gòu)效果更好。綜合以上分析，文中算法在相同的觀測個數(shù)時，重構(gòu)質(zhì)量更高。

4 結(jié)束語

文中提出的優(yōu)化算法，緊密結(jié)合觀測矩陣的列獨(dú)立性和觀測矩陣與稀疏基的非相干性兩方面對觀測矩陣進(jìn)行優(yōu)化。通過對實驗結(jié)果進(jìn)行分析，文中提出的優(yōu)化算法在重構(gòu)質(zhì)量和穩(wěn)定性這2個方面有較大優(yōu)勢，尤其在觀測個數(shù)較少時，優(yōu)勢更明顯。后期的工作將嘗試研究確定性觀測矩陣，如何采用更少的觀測個數(shù)，得到更精確的重構(gòu)結(jié)果是下一步的研究目標(biāo)。

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

[2]DINH M N.Efficientprojection for compressed sensing[C]//7th ACIS International CONFerence on Computer and Information Science.Oregon，America，2008:322-327.

[3]JULIO M D，GUILLERMO S.Learning to sense sparse signals:Simultaneous sensing matrix and sparsifying dictionary optimization[J].IEEE Transactions on Image Processing，2009，18(7):1395-1408.

[4]CAO Z.Optimized projection matrix for compressive sensing[J].EURAIP Journal Advances in Signal Processing，2010(43):55－60.

[5]YU Lifeng，LI Gang，CHANG Liping.Optimizing projection matrix for compressed sensing systems[C]// 8th International Conference on Communnicationand Signal Processing，2011:1-5.

[6]EVAGGLIAT，LISIMACHOS P K，AGGELOS K K.Useof tight frames for optimized compressed sensing[C]//20th European Signal Processing Conference，Bucharest，2012:1439-1443.

[7]VAHID A，SAIDEH F，SAEID S.A gradient-based alternating minimization approach for optimization of measurement matrix in compressive sensing[J].Signal Processing，2012，92(4):999-1009.

[8]THONG T D，LU Gan，NAM H N，et al.Fast and efficientcompressive sensing using structurally random matrices[J].IEEE Transactions on Signal Processing，2012，60(1):139-154.

[9]WEI Chen，MIGUEL R D.On the use of unit－norm tight framesto improve the use ofMSE performance in compressive sensing application[J].Transactions on Signal Processing Letter，2012，19(1):8-11.

[10]TSILIGIANNI E，KONDI L P，KATSAGGELOS A K.Use of tight frames for optimized compressed sensing[C]//20th European Signal Processing ConferenceBucharest，2012:1439-1443.

[11]趙瑞珍，秦周，胡紹海.一種基于特征值分解的測量矩陣優(yōu)化方法[J].信號處理，2012，28(5):653-658.Zhao Ruizhen，QIN Zhou.An method based on eigenvalue decomposition to optimize measurement matrix[J].Signal Processing，2012，28(5):653-658.

[12]傅迎華.可壓縮傳感重構(gòu)算法與近似QR分解[J].計算機(jī)應(yīng)用，2008，28(9):2300-2302.FU Ying hua.Reconstruction algorithm of compressed sensing and QR decomposition [J].Computer Applications，2008，28(9):2300-2302.

[13]彭玉樓，何怡剛，林斌.基于奇異值分解的壓縮感知噪聲信號重構(gòu)算法[J].儀器儀表學(xué)報，2012，33(12):2655-2660.PENG Yulou，YI Gang，LIN Bin.Reconstruction algorithm of compressed sensing noise signal based on singular value decomposition[J].Journal of Scientific Instrument，2012，33(12):2655-2660.

[14]DONOHO.Compressed sensing[J].IEEE Trans on Information Theory，2006，52(4):1289-1306.

[15]YAAKOV TSAIG，DAVID L.Extensions of Compressed sensing[J].Signal Processing，2006，86(3):533-548.

[16]DAVID L.Compressed sensing[J].Information Theory，2006，52(4):1289-1306.