依據列相關性優(yōu)化高斯測量矩陣

，，

(1.北京科技大學 自動化學院，北京 100083; 2.北京科技大學 計算機與通信工程學院，北京 100083)

依據列相關性優(yōu)化高斯測量矩陣

邊勝琴1，2，徐正光1，張利欣1

(1.北京科技大學自動化學院，北京100083; 2.北京科技大學計算機與通信工程學院，北京100083)

為了提高信號重建的精度以及稀疏度適用范圍，提出了一種新的測量矩陣優(yōu)化方法，減小測量矩陣和稀疏變換矩陣的相關性;首先，由測量矩陣和稀疏變換矩陣的乘積構造Gram矩陣；根據Gram矩陣的維數，計算互相關函數的下確界即Welch界；其次，由Welch界確定閾值，收縮Gram矩陣中大于閾值的非對角元；然后，由新得的Gram矩陣和稀疏變換矩陣反解出測量矩陣，迭代更新，從而達到減小相關性，優(yōu)化測量矩陣的目的;實驗結果表明：依據Welch界優(yōu)化測量矩陣，能快速降低壓縮感知矩陣相關性的最大值，提高OMP算法的性能，例如在誤差率為10-0.9時，原高斯隨機矩陣需要23個觀測值，算法優(yōu)化后只需16個觀測值，相對于Elad、Zhao等觀測矩陣優(yōu)化方法，文中提出的算法具有更小的重構誤差，性能和穩(wěn)定性也略有提升。

壓縮感知；測量矩陣；互相關系數；信號重構

0 引言

2004年由Donoho[1]、Candès[2-3]、及華裔科學家Tao(陶哲軒)等人提出的壓縮感知理論(compressive sensing，CS)，突破了傳統的奈奎斯特采樣定理的限制，實現了數據的獲取和壓縮同時進行，避免了對大量無用信息地采集，降低了數據存儲、傳輸和處理的開銷，在信號處理領域產生了極其重要的影響，引起了眾多科學家的關注。目前，在醫(yī)學成像、生物傳感、遙感成像、圖像處理等領域，應用也越來越廣泛。

信號的稀疏表示是壓縮感知的前提，研究人員發(fā)現根據信號的特征，總可以找到它稀疏表示形式。壓縮感知理論對信號進行稀疏表示之后，再通過一個測量矩陣，將其投影到一個低維空間，通過低維觀測值對原始信號進行重構。

1 測量矩陣

測量矩陣的性能和信號的稀疏度范圍密切相關，如何設計或者優(yōu)化測量矩陣，減少測量次數或者增大信號的稀疏度范圍，成為壓縮感知理論的一個重要方向[4-5]?，F有文獻對此展開研究，主要有約束等距性質(restricted isometry property，RIP)，零空間性質和互相關系數[6-13]等。文獻[6]提出了約束等距性質，并指出約束等距性是信號能夠完全重構的充分條件。零空間性質是指同一個測量值通過l1范數優(yōu)化得到的解具有唯一性。判定觀測矩陣是否具備RIP性質和零空間性質，需要遍歷多個子集，計算復雜度高，用來指導測量矩陣的構造不太現實。很多學者嘗試用相關性理論[9-10]，來衡量測量矩陣的性能，弱化RIP條件。研究表明，減小測量矩陣與稀疏變換矩陣的互相關系數，能夠提高壓縮感知算法的重構性能，互相關系數越小，信號精確重建所需的測量次數就越少，信號適用的稀疏度范圍就越大[11]。

感知矩陣優(yōu)化的基本思想是：保持稀疏變換矩陣不變，只對測量矩陣進行優(yōu)化，先通過測量矩陣和稀疏變換矩陣的乘積構造格拉姆(Gram)矩陣，在迭代過程中減少格拉姆矩陣的非對角元的絕對值，由新得的格拉姆矩陣和原稀疏變換矩陣更新測量矩陣。Elad[7]提出了一種閾值平均系數法，通過線性收縮Gram矩陣中絕對值大于閾值的非對角元，減小平均互相關系數，從而優(yōu)化傳感矩陣，但收縮因子和閾值點配置依靠經驗，參數不當，可能需要較多的迭代次數；Vahid[12]采用梯度下降法使Gram矩陣接近單位矩陣，但是，壓縮感知中的Gram矩陣不是滿秩矩陣，和滿秩的單位矩陣只能近似；文獻[8]定義了一種整體互相關系數，在不改變Gram矩陣的跡也即特征值之和的前提下，平均化所有非零特征值，使得非對角元整體平方和最小的方式優(yōu)化觀測矩陣，實驗表明，采用這種方式優(yōu)化，相關系數下降較快，但是容易達到局部最小值。文獻[13-15]提出了等角緊框架理論(Equiangular Tight Frame，ETF)，并給出矩陣列相關性的邊界，邊界最小值也即下確界稱為Welch界，采用構造方法生成的矩陣才具有邊界最小值。

本文對測量矩陣優(yōu)化方法進行了研究，提出了一種新的優(yōu)化方法，減小Gram矩陣中較大的非對角元，經過較少的迭代次數，就能將Gram矩陣的相關系數降到Welch界附近；在相同迭代次數的情況下，與文獻[7-8]中的方法相比，最大相關性和平均相關性，都較低，從而在信號重構時，能夠提高信號的重構精度，降低壓縮感知系統對信號稀疏度的要求。

2 基本原理

2.1 壓縮感知理論

長度為N的信號x∈N，在給定的一組標準正交基Ψ=[Ψ1,Ψ2,…,ΨN]上表示為：

或x=Ψα

(1)

若N×1維的列向量α中只有K(K?N)個非零系數，則稱信號x在基矩陣Ψ下是K稀疏的。若非零系數排序后，按指數級衰減趨于零，則稱信號是可壓縮的。

壓縮感知的過程是：長度為N的信號x，通過測量矩陣Φ∈M×N進行投影，得到觀測值y∈M(M?N)，再由觀測值y對信號進行精確重構。

觀測過程表示為：

y=Φx

(2)

將稀疏表示(1)代入(2)式，得:

y=Φx=ΦΨα=Αα

(3)

Φ∈M×N稱為測量矩陣或投影矩陣，Ψ∈N×N稱為稀疏矩陣，Α=ΦΨ，Α∈M×N稱為感知矩陣。

當原始信號x本身就是K稀疏的信號時，稀疏變換矩陣Ψ是單位矩陣，測量矩陣Φ滿足一定條件，信號x就可以通過最小l0范數優(yōu)化問題精確重構。

s.t.y=Φx

(4)

求解(4)式，一種方法是采用匹配追蹤類算法，通過每次選取局部最優(yōu)解來逐步逼近原始信號，由于復雜度低，獲得廣泛應用，主要算法包括匹配追蹤(matching pursuit，MP)、正交匹配追蹤(orthogonal matching pursuit，OMP)、分級正交匹配追蹤(stagewise orthogonal matching pursuit，StOMP)、正則化的正交匹配追蹤(regularized orthogonal matching pursuit，ROMP)、子空間追蹤(subspace pursuit，SP)、壓縮采樣匹配追蹤(compressive sampling matching pursuit, CoSaMP)等。另外一種方法是將最小l0范數轉化為最小l1范數問題，采用凸優(yōu)化方法求解，主要包括基追蹤法(basic pursuit，BP)，欠定系統局部解法(focal underdetermined system solver，FOCUSS)，梯度投影重建法(gradient projection for sparse reconstruction，GPSR)等。凸優(yōu)化方法在測量矩陣滿足限制等容條件下，能夠精確重建所有的稀疏信號，而且所需的測量次數也較小，但是計算復雜度高，對于大尺度問題，計算時間長，不太實用。

對于測量矩陣Φ來說，列數始終大于行數，列數越小，說明對采樣率的要求就越低，并且要求Φ和Ψ不相干，Φ的行向量不能由Ψ的列向量進行稀疏表示。高斯隨機矩陣、貝努利隨機矩陣、部分傅里葉隨機矩陣及稀疏投影矩陣等，幾乎與所有的稀疏基都不相關，可以作為普適性的觀測矩陣。但是文獻[14]指出，測量矩陣和稀疏字典對應的Gram矩陣并沒有逼近理想模型，信號的重構精度和稀疏度范圍還有可以提升的空間。

2.2 Gram矩陣和互相關系數

定義1(相關性)：矩陣A的互相關系數，也即A的相干性，是矩陣A的任意兩列的內積最大值。

其中:ai和aj都是矩陣Α的列向量。

由定理2可以看出，互相關系數和信號的稀疏度成反比例關系，互相關系數越小，精確重構信號的稀疏度范圍就越大，反之，互相關系數越大，精確重構信號的稀疏度范圍就越小，對維數確定的觀測矩陣來說，所需的觀測值數目相對較大。定理1表明Welch界是相關系數的下確界，只有任意兩列內積都相等時，才能達到，對于高斯隨機矩陣來說，不可能達到最小值，但是可以將Welch界作為調整Gram矩陣的一個重要參數，對觀測矩陣進行優(yōu)化。

2.3 依據Welch界優(yōu)化觀測矩陣

觀測矩陣優(yōu)化的具體過程是：首先，由測量矩陣和稀疏變換矩陣的乘積，計算Gram矩陣及其相關性的下確界(Welch界)，并將2倍Welch界作為優(yōu)化算法的閾值；其次，縮小Gram矩陣中大于閾值的非對角元，縮小到原來的一半，保留小于閾值的元素；再者，對修改后的Gram矩陣進行奇異值(singular value decomposition，SVD)分解，保證奇異值個數不變，將Gram矩陣的秩收縮，再由Gram矩陣和稀疏字典計算出測量矩陣；反復迭代，當Gram矩陣中不存在大于2倍Welch界的非對角元時，調整閾值，對較大的非對角元進行微調，重復上述過程，直到完成迭代次數。

整體算法步驟如下：

1)初始化高斯隨機矩陣Φ∈M×N，Φ為列滿秩矩陣，選擇稀疏矩陣為Ψ∈N×N，感知壓縮矩陣Α=ΦΨ，Α的維數為M×N，迭代次數t=1；

如果Gram矩陣非對角線上的元素都小于給的的閾值，按下式修改G的非對角元，

3 實驗結果

3.1 優(yōu)化前后Gram矩陣比較

構造50*120的高斯隨機矩陣作為測量矩陣，采用絕對稀疏的隨機信號作為重建對象，稀疏矩陣即為單位矩陣，測量矩陣和稀疏矩陣之間的相關性按定理1計算，則最小互相關系數也就是Welch界為：W(Α)=0.1085，采用上述方法對矩陣進行優(yōu)化，大于閾值的元素縮小到原來的0.5倍，小于閾值的元素保持不變，迭代30次，大于閾值的元素不存在時，將絕對值最大的非對角元的0.9倍作為新的閾值，繼續(xù)迭代，直到達到設定的迭代次數。結果如圖1所示，(a)為原始高斯的矩陣的Gram矩陣，(b)為訓練后測量矩陣的Gram矩陣，(c)為訓練后測量矩陣的Gram矩陣的絕對值；參照色彩條狀軸，比較圖(a)(c)可知，Gram矩陣非對角元素的值在逐漸減小，而且十分接近，說明優(yōu)化后矩陣相干性絕對值縮小了，而且范圍比較優(yōu)化前要集中，這對信號重構十分有用。

圖1 Gram矩陣比較

3.2 觀測矩陣優(yōu)化前后直方圖比較

測量矩陣和稀疏字典的乘積Gram矩陣共有14 400個元素，除去對角線全1的元素還剩14 280個元素，采用上述方法對矩陣進行優(yōu)化，將非對角元素的絕對值分成80組，對Gram矩陣非對角元的分布狀況進行分析，如圖2所示。

圖2 Gram矩陣非對角元絕對值的直方圖

橫坐標表示Gram矩陣的非對角元的絕對值，縱坐標表示相應區(qū)間的元素的個數，可以看出，優(yōu)化前在0附近的元素最多，相關性最大值為0.61；迭代100次后在0.16附近的元素最多，相關性最大值縮小到0.2，元素的分布更加密集，說明本文提出的方法確實能夠降低觀測矩陣的相關性，使得相關性分布在Welch值附近。對Gram矩陣采用不同方法優(yōu)化后，最大相關性和平均相關性也不同，對上述矩陣迭代9次進行比較，最大相關性即乘積矩陣中的非對角線元素的最大值，平均相關性是乘積矩陣中前20個較大的非對角元的均值，由表1可以看出，在前9次的迭代運算中，本文方法得到的最大相關性和平均相關性均小于由文獻[7]和文獻[8]得到的值；比較前30次的迭代，本文方法得到的平均相關性始終最低。

3.3 OMP和BP算法優(yōu)化前后比較

為了驗證觀測矩陣優(yōu)化算法在整個壓縮過程中的有效性，本文進行了兩次實驗。

實驗一：選取50 000個向量，向量長度為120，每個稀疏向量包含4個非零項，非零項的位置隨機分布，進行壓縮感知實驗，觀測次數的取值范圍為[16,32]，分別BP算法和OMP算法對觀測信號進行重構，測試重構誤差隨觀測次數的變化，檢驗投影矩陣的性能。

表1 相關性比較

分別采用4種不同觀測矩陣進行比較，原高斯隨機矩陣，文獻[7]優(yōu)化方法、文獻[8]優(yōu)化方法，以及本文優(yōu)化方法，當重構錯誤次數大于200時，迭代終止，前200次的平均誤差作為以后的重構誤差。實驗結果如下圖所示，橫坐標為測量次數，縱坐標為重構誤差的對數表示，實驗部分代碼來自SparseLab實驗箱。

由圖3(a)可以看出，優(yōu)化投影矩陣能大幅提高OMP算法的性能，例如在10-0.9的重構誤差下，原高斯矩陣需要23個觀測值，而文獻[7]和本文方法優(yōu)化的觀測矩陣只需要16個，說明觀測矩陣優(yōu)化后，每個觀測值包含更多的信息量，投影矩陣的優(yōu)化，能夠提升信號的重構性能。圖3(b)表示BP算法在觀測矩陣優(yōu)化下信號重構性能略有提升，4種觀測矩陣重構誤差相差不是很大。兩個圖的共同點是隨著觀測次數的增加，重構誤差率在逐漸減小，相對于其他3種測量矩陣，文中提出的優(yōu)化算法，具有更小的重構誤差率。

圖3 不同優(yōu)化算法隨觀測次數m的變化

實驗二：同樣選取50 000個向量，向量長度為120，固定觀測次數為25次，向量的稀疏度從3逐漸增加到9，測試重構誤差隨稀疏度的變化，檢驗投影矩陣的性能。

原始投影矩陣為25×80的高斯隨機矩陣，對以上4種優(yōu)化算法進行仿真，當重構錯誤次數大于200時，迭代終止，前200次的平均誤差作為以后的重構誤差。橫坐標為信號的稀疏度，縱坐標為重構誤差的對數表示，結果如圖4所示。

圖4 不同優(yōu)化算法隨信號稀疏度K的變化

圖4表明，隨著信號稀疏度的增加，重構誤差率在逐漸增大，觀測矩陣優(yōu)化后，文獻[7] 和本文方法比原高斯隨機矩陣，采用OMP算法重構時，具有更小的重構誤差，文獻[8]提出的優(yōu)化方法在本文的實驗條件下，性能不及原高斯隨機矩陣。使用BP算法重構時，重構誤差隨著稀疏度的增加，差異性越來越小，本文優(yōu)化方法略優(yōu)于其他算法。

4 結論

本文提出了一種優(yōu)化測量矩陣算法，由測量矩陣和稀疏變換矩陣的乘積得出Gram矩陣，Gram矩陣的非對角元即測量矩陣的行向量與稀疏變換矩陣的列向量的相關性，根據代數中的定理可知，滿秩矩陣列向量之間的相關性存在最小值，稱為Welch界。將Gram矩陣的Welch界作為優(yōu)化測量矩陣的目標，縮小Gram矩陣的非對角元；再對調整后的Gram矩陣進行奇異值分解，保證奇異值的個數不變，反解出新的測量矩陣；重復上述過程，直到達到設定的迭代次數。實驗證明，Gram矩陣的相關性能快速地下降到Welch界附近，由此得到的測量矩陣，在進行圖像重建時，壓縮感知性能和重建質量都有不同程度的提高。

[1] Donoho D L, Elad M, 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 EJ, Romberg J K, Tao T. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information [J]. IEEE Transaction on Information Theory, 2006, 52(2): 489-509.

[3] Candes E J, Romberg J K, Tao T. Stable signal recovery from incomplete and inaccurate measurements [J]. Communications on pure and applied mathematics, 2006, 59(8): 1207-1223.

[4] 鄭 紅, 李 振. 壓縮感知理論投影矩陣優(yōu)化方法綜述 [J]. 數據采集與處理, 2014(1): 43-53.

[5] 戴瓊海，付長軍，季向陽. 壓縮感知研究 [J]. 計算機學報, 2011, 34(3): 425-434.

[6] Candes E J. The restricted isometry property and its implications for compressed sensing[J]. Comptes Rendus Mathematique, 2008, 346(9-10): 589-592.

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

[8] 趙瑞珍，秦 周，胡紹海. 一種基于特征值分解的測量矩陣優(yōu)化方法 [J]. 信號處理, 2012, 28(5): 653-658.

[9] Baraniuk R, Davenport M, DeVore R, et al. A simple proof of the restricted isometry property for random matrices [J]. Constructive Approximation, 2008, 28(3): 253-263.

[10] Xu J, Pi Y, Cao Z. Optimized projection matrix for compressive sensing[J]. EURASIP Journal on Advances in Signal Processing, 2010, 2010(1): 560349.

[11] Donoho D L, Elad M. optimally sparse representation in general(non-orthogonal)dictionaries via minimization [J]. Proceeding of the National Academy of Sciences, 2003, 100(5): 2197-2202.

[12] 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.

[13] Tropp J A, Dhillon I S, Heath R W, et al. Designing structured tight frames via an alternating projection method [J]. IEEE Transactions on information theory, 2005, 51(1): 188-209.

[14] Strohmer T, Heath R W. Grassmannian frames with applications to coding and communication [J]. Applied and computational harmonic analysis, 2003, 14(3): 257-275.

[15] Sustik M A, Tropp J A, Dhillon I S, et al. On the existence of equiangular tight frames [J]. Linear Algebra and its applications, 2007, 426(2-3): 619-635.

[16] Tillmann A M, Pfetsch M E. The Computational Complexity of the Restricted Isometry Property, the Nullspace Property, and Related Concepts in Compressed Sensing [J]. IEEE Transactions on Information Theory, 2014, 60(2): 1248-1259.

[17] 吳光文，張愛軍，王昌明. 一種用于壓縮感知理論的投影矩陣優(yōu)化算法 [J]. 電子與信息學報, 2015(7): 1681-1687.

[18] 李哲濤，潘 田，朱更明,等. 低冪平均列相關性測量矩陣構造算法 [J]. 電子學報, 2014(7): 1360-1364.

[19] 柯家龍，李繼樓. 壓縮感知中的投影矩陣優(yōu)化算法 [J]. 計算機技術與發(fā)展, 2015(3): 95-98.

OptimizeMeasurementofGaussianMatrixBasedonColumnCorrelation

Bian Shengqin1，2, Xu Zhengguang1, Zhang Lixin1

(1.School of Automation and Electrical Engineering, Beijing University of Science and Technology, Beijing 100083, China；2.School of Computer and Communication Engineering, Beijing University of Science and Technology, Beijing 100083, China)

In order to improve the accuracy of signal reconstruction and the application range of sparsity, a new method of measuring matrix optimization is proposed. First, multiply measurement matrix and sparse transformation matrix to construct a Gram matrix, and calculate the minimum value of mutual coherence, that is the Welch bound; Secondly, set a threshold based on Welch bound and reduce the elements of the non-diagonal of the Gram matrix; Third, produce new projection matrix from inverse solution of new Gram matrix and sparse transformation matrix iteratively, so as to achieve the purpose of reduction the mutual coherence and optimizing the measurement matrix. Experiments in the last show: measurement matrix based on the Welch optimization can rapidly reduce the maximum value of the compressed sensing correlation matrix and improve the performance of OMP algorithm, such as when error rate is 10-0.9, the original Gauss random matrix need 23 observations, but our optimized matrix only 16 observations. On the whole, when compared with the optimization method of Elad and Zhao, the reconstruction error of the algorithm in this paper is smaller, the performance and stability are slightly improved.

compressive sensing(CS); projection matrix; mutual coherence; signal reconstruction

2017-04-12；

2017-05-08。

邊勝琴(1977-)，女，河北保定人，博士研究生，主要從事壓縮感知方向的研究。

徐正光(1952-)，男，教授，博士研究生導師，主要從事控制科學與工程、圖像處理方向的研究。

1671-4598(2017)11-0141-05

10.16526/j.cnki.11-4762/tp.2017.11.036

TP911.7.1

A