基于奇异值分解的压缩感知观测矩阵优化算法

李　周，崔　琛

(国防科技大学 电子对抗学院，合肥 230037)(*通信作者电子邮箱17756587331@163.com)

0　引言

压缩感知(Compressive Sensing, CS)[1-3]在采样率远小于奈奎斯特采样率的条件下获取信号的离散样本，然后通过非线性优化算法重构出原始信号。重构信号所需的采样率并不取决于信号的带宽，而与信号的稀疏度密切相关，因此CS有效降低了信号获取、存储及传输的代价，引起了越来越多学者的关注。

信号的稀疏表示、观测矩阵的构造、重构算法是CS理论中三个主要研究方向，其中观测矩阵是影响压缩感知性能的关键因素之一[4]。观测矩阵构造的目的是如何采样得到观测值，并能从观测值中重构出原始信号。文献[5-7]对精确重构所需观测矩阵的约束条件进行了研究：文献[5]提出了零空间性质(Null Space Property, NSP)，但在存在噪声的情况下NSP并不能保证信号的精确重构；文献[6]提出了有限等距性质(Restricted Isometry Property, RIP)，但判断观测矩阵是否符合RIP需要计算观测矩阵n列中任意组合的K列，即需要计算观测矩阵各列的组合，实际用于观测矩阵的分析非常困难；文献[7]提出了相关性的概念，其含义是观测矩阵与稀疏基之间的相关程度，通过降低相关性可以减少CS的重构误差和精确重构原始信号所需的观测数。

文献[7]通过线性收缩Gram矩阵中绝对值大于限定阈值的非对角元来降低相关性，取得了较好的实验效果。但是，该算法在收缩Gram矩阵时会改变矩阵的秩[8-9]；同时在求解观测矩阵时需要求解稀疏基的逆矩阵，但由于稀疏基可能奇异，此时需要利用稀疏基的Moore-Penrose广义伪逆来求解新的观测矩阵，造成了算法中计算量大且精度受限的问题[4,8]。文献[9]延续了文献[7]之前的思路，不同之处在于增加了一个保留前一次矩阵优化信息的操作，其目的在于使每一次矩阵更新的变化量减少，该算法继承了文献[7]的缺点，同时需要根据经验手动设置一个权重系数，用来衡量本次矩阵优化结果与前一次矩阵优化结果之间的比重。

文献[10]提出将Gram矩阵非对角元素的平方和作为整体互相关系数，用来表示观测矩阵的整体性能；在此基础上，文献[11]引入了α紧框架的概念，平均化Gram矩阵的非零特征值，减小了整体相关系数。但在由Gram矩阵求解观测矩阵时，文献[10-11]仍采用了求稀疏基Moore-Penrose广义伪逆矩阵的做法，同样造成了算法中计算量大且精度受限的问题。

针对现有算法从优化后的Gram矩阵求解观测矩阵时出现的相关系数较大与需要求广义伪逆矩阵的问题，本文借鉴了K-SVD(K-Singular Value Decomposition)算法中逐行更新优化目标矩阵的思想，在利用现有算法得到优化后的Gram矩阵的基础上，提出了一种基于奇异值分解(Singular Value Decomposition, SVD)的观测矩阵优化算法：该算法通过求解等价变换后的目标函数对观测矩阵行向量的导数得到目标函数取极值时行向量的值，并通过对误差矩阵进行奇异值分解在上述行向量的值中选出使得目标函数取最值时行向量的解析式，之后在每轮迭代中对观测矩阵的每一行分别使用上述的行向量解析式进行优化。

1　观测矩阵相关性和相应算法

1.1　观测矩阵相关性

假定离散信号r∈Rl×h，若r中最多含有K个非零值且K≤h，那么r就称作K-稀疏的。将K-稀疏的信号集合记为：

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

(1)

若r本身为非稀疏的，但可以经过一个稀疏基Ψ作如下变换：

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ΦΨ称为感知矩阵。

对于任意两个不同的稀疏信号r1,r2∈ΛK，必须满足Xr1≠Xr2，否则仅仅根据观测值无法区分r1,r2，因此感知矩阵需要满足一定的条件：对于任意K-稀疏的信号，只要其稀疏度K满足下式，信号即可准确地恢复出来[12]。

(4)

其中μ(X)指X任意两列xi和xj之间的内积绝对值的最大值[7]，计算公式如下：

(5)

然而μ(X)仅表示观测矩阵与稀疏基之间列向量最大的相关性，是一种局部的相关性，因为在X只有某两列的相关性比较大而其余各列之间相关性比较小的情况下，就会出现相关系数μ(X)很大而感知矩阵X的性能并不差的情况。

因此文献[7]提出了表示总体相关性的t-平均相关性μt-av，定义为X所有列向量两两之间的内积绝对值中大于限定阈值t的部分的平均值，或者是X所有列向量两两之间的内积绝对值中最大的t%部分的均值，其定义式如下：

(6)

其中:gij为Gram矩阵G=XTX中的元素，gij为矩阵X的第i列xi与第j列xj的内积，Nav为Hav中的元素数，即Gram矩阵非对角元|gij|大于t的数目或者所有|gij|中最大的t%部分的元素数目。仅阐述μt-av为X任意两列内积绝对值大于t的平均值时Hav的定义：

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

(7)

1.2　观测矩阵优化的相关算法及问题分析

为降低X任意两列的相关性，即减少Gram矩阵非对角元素的值，文献[7]采用如下阈值函数对Gram矩阵G中绝对值大于限定阈值的非对角元进行线性收缩：

(8)

其中：t为阈值，γ为收缩因子。

(9)

其中：A=ΦΨ。在Ψ逆矩阵存在时，Φ=AΨ-1；但在稀疏基Ψ可能奇异导致其逆矩阵不存在时，需要利用稀疏基的Moore-Penrose广义伪逆来求解观测矩阵Φ=AΨ+，从而带来计算量与精度的问题[4,8]。

文献[6]的阈值函数只能使局部比较大的非对角元素减小。为整体减小Gram矩阵中的非对角元素，文献[11]通过平均化Gram矩阵的非零特征值，使得矩阵非零特征值的平方和减小，降低了整体相关系数。但在由Gram矩阵求解观测矩阵时，仍然需要求解稀疏基的广义逆矩阵，故存在着与文献[7]相同的问题[4,8]。

1.3　问题的解决思路

在现有算法优化Gram矩阵后，本文借鉴了文献[13]提出的K-SVD算法中逐行优化目标矩阵的思想从优化后的Gram矩阵求解观测矩阵。

K-SVD是一种通过逐行优化字典矩阵对信号进行稀疏表示的方法，具体目标为：

(10)

其中:Y为要表示的信号，D为所求的字典，X为稀疏矩阵。X与Y按列对应，表示D中的元素以xi为系数线性组合就可得到Y，而K-SVD的目的是在X和Y已知的情况下更新字典D满足上述条件。

(11)

式(11)可以看作把第k个分量剥离后表达式会产生空洞，目的是找到一个新向量以更好地填补这个洞：假设除第k项外其余项是固定的，之后对Ek进行奇异值分解得到Ek=UΛVT，其中U和V的列矢量均是正交基，Λ是对角矩阵。若Λ的对角元素从大到小排列，则表示Ek的能量分量主轴在相应几个正交方向上由大到小分配，取U的第一个列向量来表示di，取V的第一个列向量与Λ的第一个元素的乘积表示xi，至此完成了字典一个条目的更新。

2　观测矩阵的优化算法

在利用现有算法得到优化后的Gram矩阵的基础上，借鉴K-SVD算法中逐行更新优化目标矩阵的思想，本章利用稀疏基的奇异值分解对目标函数进行等价变换，通过求解等价变换后的目标函数对观测矩阵行向量的导数得到目标函数取极值时行向量的值，并通过对误差矩阵进行奇异值分解在上述行向量的值中选出使得目标函数取最值时行向量的解析式，最后利用所求出的行向量解析式逐行迭代优化观测矩阵。

2.1　观测矩阵行向量的优化

(12)

假设稀疏基Ψ的秩为NΨ，其奇异值分解为：

(13)

将奇异值分解式代入目标函数中，可得：

(14)

根据F范数的酉不变性，在式(14)中左乘VΨT，右乘VΨ，同时设ΦUΨ其中为的前NΨ列。设其中为的前NΨ行与前NΨ列，将和代入式(14)可得：

(15)

由式(15)易得：

(16)

(17)

(18)

式(18)对ωj求导可得：

(19)

导数置0便可得到一系列极值点：

(20)

易得：

(21)

(22)

2.2　观测矩阵的迭代优化算法

在现有算法得到优化后的Gram矩阵的基础上，本小节利用2.1节求出的观测矩阵行向量的解析式，采用逐行更新的方法优化观测矩阵。

输入初始观测矩阵Φ∈Rm×n，离散傅里叶变换基Ψ∈Rn×h，迭代次数Iter，退出阈值μ0;

输出观测矩阵。

fork=1:Iter

1)

2)

forj=1:m

end for

3)

4)

计算，若μt-av与上一轮迭代的μt-av之差小于μ0则退出循环，否则转至1)继续循环

end for

3　实验结果及其分析

本文在Matlab平台上对算法进行仿真实验，仿真选择高斯随机矩阵作为初始观测矩阵，离散傅里叶变换基作为稀疏基，比较文献[7]所提的算法CSElad、文献[11]所提的算法CSTsiligianni和经文献[11]中的方法优化Gram矩阵后使用本文所提的算法CSImproved-Tsiligianni在相关性、重构误差和运行时间三方面的性能。

仿真中相关参数如下：m=30，n=l=100，稀疏度K=10，迭代次数Iter=100。Gram矩阵收缩变换时非对角线的限定阈值t=30%，收缩因子γ=0.5，CSImproved-Tsiligianni的退出阈值μ0=10-2。为减少随机性对实验结果造成的影响，每次仿真均为1 000次的蒙特卡罗仿真。

3.1　相关性

仿真实验一 ：CSElad、CSTsiligianni和CSImproved-Tsiligianni都是迭代进行的算法，去除改进算法第四步迭代退出的步骤，将三种算法每次迭代时Gram矩阵的t-平均相关性μt-av进行对比，可以看出在各个算法迭代时t-平均相关性的变化趋势，进而对比出三种算法在t-平均相关性这一方面的优劣性。t-平均相关性μt-av中的参数t=20%，其含义是Gram矩阵非对角元中最大的20%部分的均值。

图1　不同算法的μt-av随迭代的变化Fig. 1　Changes in μt-av of different algorithms with iterations

由图1可知，随着迭代的进行，CSImproved-Tsiligianni的μt-av逐渐减小而趋于稳定，而CSElad和CSTsiligianni的μt-av变化范围较大;同时，改进算法的μt-av比改进之前的算法小，这是由于改进算法利用优化后的观测矩阵行向量的解析式逐行优化观测矩阵。

3.2　重构误差

(23)

仿真实验二：观测次数m对压缩感知的重构误差的影响。仿真了观测次数m从25增加到45时三种算法重构误差的变化情况，如图2所示。可以看出，在给定信号长度和稀疏度的前提下，m越大，重构误差越小；m越小，重构误差越大。

图2　观测次数变化时三种算法的MSE变化情况Fig. 2　Changes in MSE of three algorithms with observation number

仿真实验三：稀疏度K对算法的重构性能的影响。为了得出稀疏度K对各个算法重构性能的影响，仿真稀疏度K由10变化到20时三种算法的MSE的变化情况，如图3所示。可以看出，在给定观测次数和信号长度的前提下，稀疏度越高，重构误差越大。

图3　稀疏度变化时三种算法的MSE变化情况Fig. 3　Changes in MSE of three algorithms with signal’s sparsity

由仿真实验二与三可知在K或m变化时CSImproved-Tsiligianni的MSE小于CSElad和CSTsiligianni，这是由于CSImproved-Tsiligianni产生的观测矩阵与稀疏基的相关性小于CSElad和CSTsiligianni，可以更好地保持不同K稀疏向量之间的距离。

3.3　运行时间

为了研究本文所提算法的复杂度，该仿真实验统计CSImproved-Tsiligianni优化观测矩阵所需的运行时间，并与CSElad和CSTsiligianni的运行时间进行比较。虽然算法的运行时间不能严格地定义算法复杂度，但仍可以在一定程度上对算法的复杂度作出描述。仿真环境为2.53 GHz英特尔i5四核处理器、4 GB内存Windows 10系统下的Matlab R2014a。图4给出了各算法在信号长度l和观测次数m变化时相应的运行时间。

图4　三种算法的运行时间对比Fig. 4　Running time comparison of three algorithms

由图4可知信号长度和观测次数增加时算法的运行时间随之增加，其中CSElad和CSTsiligianni的运行时间小于CSImproved-Tsiligianni的运行时间，这是因为CSImproved-Tsiligianni在优化观测矩阵时对每一行均采用奇异值分解来求解优化后的行向量，而奇异值分解耗时是较长的。不过相对于在相关性与重构误差方面的优势而言，CSImproved-Tsiligianni增加的运行时间还是可以接受的。

4　结语

本文对压缩感知中观测矩阵的优化问题进行研究，借鉴K-SVD算法中逐行更新目标矩阵的思想，在现有算法得到优化后的Gram矩阵基础上，提出了一种基于奇异值分解的观测矩阵优化算法。仿真结果表明：在可接受的运算量下，本文所提算法在观测矩阵与稀疏基的相关性方面优于改进前的算法，从而提高了重构精度。如何优化阈值函数来构造性能更优的Gram矩阵是下一步的研究方向。

参考文献:

[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]王强,张培林,王怀光,等.压缩感知中测量矩阵构造综述[J].计算机应用,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]郑红,李振.压缩感知理论投影矩阵优化方法综述[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]赵瑞珍,秦周,胡绍海.一种基于特征值分解的测量矩阵优化方法[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.