基于近似l0范数的稀疏信号重构

聂栋栋 弓耀玲

(燕山大学理学院 河北秦皇岛 066004) (niedd@ysu.edu.cn)

自压缩感知(compressed sensing, CS)信号处理理论[1]提出以来，信号重构作为其关键的环节，一直是该理论的研究重点.信号重构的目标是根据非自适应线性测量的低维数据，准确而有效地恢复原始采样获取的高维数据.在压缩感知理论框架下，信号在特定基下或一定的变换下可稀疏表示，则信号重构问题可转化为稀疏信号的重构.如何求欠定方程组最稀疏的解是稀疏重构问题的本质，即l0范数的最小化问题.

然而直接求解欠定方程组的最稀疏解是NP难问题，处理起来非常棘手.在实际中，贪婪算法是最常见的近似求解方法[2-5]，如正交匹配追踪法(orthogonal matching pursuit, OMP)[2]、梯度追踪法(gradient pursuit, GP)[3]、子空间追踪法(subspace pursuit, SP)[4]等，其基本思想是通过逐步迭代找到未知信号的支撑集，原信号的估计值即可利用最小二乘法得到.贪婪类算法需要较少的计算量，但是通常估计精度不高且需要较多的观测值.

另一种思路是寻找新的模型替代l0范数的求解.这种方法分为2步：

1) 需要找到合适的函数近似估计l0范数，由于l0范数是x的不连续函数，实际计算中不好处理，很多文献常用l1范数来近似l0范数，如基追踪法(basis pursuit, BP)[6]，其依据是l1范数最小化与l0范数最小化在满足一定条件时具有相同的稀疏解，可以通过线性规划方法求解.然而，l1范数不能充分反映向量x的稀疏特性，重构时需要较多的测量值，计算量较大.高斯函数类也是近似l0范数经典的函数之一，基于光滑l0范数算法(smoothedl0norm, SL0)[7]是用一个连续平滑函数的极限来近似l0范数，如文献[7-8]采用

(1)

对l0范数进行近似，其中σ为趋于0的正数.文献[9-10]则在式(1)的基础上进一步改进，采用式双曲正切函数

(2)

近似l0范数，其中σ为很小的正数.在文献[11-13]中，用一组反正切函数

(3)

近似l0范数，其中σ为接近于0的正常数.文献[14]是采用复合三角函数近似l0范数.但以上函数表达式均稍复杂，使得迭代计算复杂度较高.文献[15]的快速稀疏信号恢复算法(AL0算法)提出用函数

(4)

2) 对新的函数模型选择不同的线性规划或优化算法求解.对于算法实现，NSL0算法选用修正牛顿法求解最优化问题，需要对海森矩阵进行修正，构造一个可近似海森矩阵逆的正定对称阵作为迭代方向，一定程度上提高了收敛速度.AL0算法将问题分解为2个最小化问题分别求解，然后交替迭代，将解投影到可行域.本文选取高精度的牛顿迭代法，将问题转为无约束优化问题后，对无约束问题直接应用牛顿算法，迭代得到最优解，不需要进一步对解进行投影.而且得到的牛顿方向为下降方向，不需要进一步修正，从而减少了误差影响，保证了精度.

信号的恢复质量在算法的不断改进下有所提升，但仍存在存储量大、计算复杂度高、精确度不够高等问题，需要进一步研究和改进.本文在之前l0范数工作的基础上，结合AL0算法快速收敛和牛顿迭代法高精度的优点，提出新的基于近似l0范数的稀疏信号重构算法，用一个简单的分式函数的和来近似估计l0范数，通过牛顿迭代算法求得稀疏解.最后，通过数值仿真和已有的相似算法及经典的OMP算法进行比较，并分析了本文算法恢复的信号在重构误差、信噪比方面的优势.

1 压缩感知理论

令z∈n是传统采样得到的信号，经过正交基或紧框架ψ变换后的系数是稀疏的，即z=ψx，其中x∈n且为稀疏度.在无噪声干扰的情形下，测量向量y可以表示为y=Φψx，其中y∈m，Φ为m×n维测量矩阵，ψ为n×n维基矩阵.令A=Φψ，则重构问题简化为在矩阵A和测量向量y已知的条件下，求方程组y=Ax最稀疏解的问题，即l0范数的最小化问题：

(5)

在实际应用中往往会不可避免地受噪声的影响，测量值是不精确的，令γ表示噪声，则测量值进一步表示为y=Ax+γ.

2 基于近似l0范数的稀疏信号重构

2.1 l0范数近似

本文采用一个形式较为简单的分式函数

(6)

来近似l0范数，其中δ是很小的正数，实验中参数δ可取为一组下降且趋于0的序列.显然有:

(7)

令:

(8)

根据l0范数的定义可知，当δ无限趋近于0时，F(x)的值就无限接近于向量x中非0元的个数，即x的l0范数，故有:

(9)

因此,我们可以用式(8)来近似l0范数进行迭代计算.式(6)的优势在于形式简单，尤其在算法迭代时简化了计算，同时又不降低精度.

2.2 算法实现

本文将最优化问题式(5)等价变换为拉格朗日形式最小化问题：

(10)

其中,λ是拉格朗日乘子.然后对式(10)进行求解.由2.1节可知，式(10)可转化为

(11)

式(11)是无约束的最优化问题，最优解x*满足优化条件:

(12)

进一步有:

(13)

其中,f′(x)为f(x)的导数，AT为A的转置.

当δ→0，

(14)

为方便表示，令:

(15)

Hf=ATA+λXf,

(16)

(17)

由式(17)可以得到Hessian矩阵可表示为

(18)

显然，Hessian矩阵Hf为正定矩阵，故得到的牛顿方向必然为下降方向，可以利用牛顿迭代[16]求解式(11)得：

(19)

式(19)即算法的迭代式，hx为迭代步长，初始解x0取最小二乘解.

鉴于求逆运算会随着矩阵维数的增高计算量迅速加大，为尽量降低算法的计算复杂度，考虑在算法迭代中设法降低维数.

由于稀疏化后的信号向量的大部分值为0或接近于0，故可以只迭代计算对应变量的关键元素，如变量中绝对值最大的l个分量，忽略0元素及接近于0的元素.

令：

(20)

其中，参数α∈[0,1)，通常可取α=0.01.

xS=(xi),i∈S,

(21)

AS=(ai),i∈S,

(22)

其中，ai为A的第i列，则算法迭代的Hf可替换为

(23)

相应地:

(24)

2.3 算法的计算复杂度分析

在算法实现过程中，直接计算式(19)时，涉及到矩阵和向量的乘积及矩阵的求逆，假设A是n×n矩阵，y为n×1向量，则ATy的计算复杂度为O(n2)，对Hf求逆的计算复杂度为O(n3).

经过降维后，令l为式(20)中S的元素个数，则如式(23)，对Hf求逆的计算复杂度变为O(l3)，而通常情况下，l≪n，从而使得本文算法每次迭代的计算复杂度维持在O(n2)，这与SL0，NSL0，OMP，AL0算法的计算复杂度均是相当的.

2.4 算法描述

本文算法的伪代码描述如下：

算法1. 基于近似l0范数的稀疏信号重构算法.

输入：矩阵A、测量向量y、正则化参数λ、参数δ、终止条件δmin；

输出：重构信号x.

初始化：初始解x0=AT(AAT)-1y.

迭代：

步骤1. 令:

S={t:|xt|>0},
xS=(xi),i∈S,
AS=(ai),i∈S;

步骤2. 更新矩阵：

步骤3. 迭代x，

步骤4. 更新δ，δ=δ×0.5；

步骤5. 判断算法是否终止，若δ>δmin不满足输出迭代的结果x，否则返回步骤1继续迭代.

2.5 算法的理论证明

在不考虑噪声存在的情况下，算法的理论证明如下.

证明. 由矩阵逆的分块形式知：

(25)

(26)

将式(25)带入式(26)并经过推导可以得到:

证毕.

(ATA+λXf)-1ATy=x.

(27)

(28)

证毕.

3 仿真实验

为验证本文算法的有效性，对一维随机稀疏信号进行实验，并与经典算法OMP及相似的AL0算法等进行了比较.下面所有数值仿真中，矩阵A取高斯随机矩阵, 测量值向量通过y=Ax+γ生成, 其中γ表示高斯白噪声.实验采用信噪比(signal noise radio,SNR)作为评估算法精度的性能指标，定义为

(29)

其中,xopt表示对源信号x的稀疏估计值.相对误差定义为

(30)

实验中随机产生稀疏信号x和矩阵A，实验结果在参数设定的情况下进行100次求平均.实验条件为ThinkpadT410i笔记本电脑Microsoft Windows 7操作系统MatlabR2010a 处理平台的运行环境.

实验1. 测试不同算法的整体恢复性能.在相同取值下，信号长度n=256，稀疏度k=20，测量值向量长度m=128，噪声均方差取为0.01，分别测试本文算法(Proposed)、AL0算法、SL0算法、NSL0算法及OMP算法所重建的信号在信噪比、重构误差和重构时间方面的对比.结果如表1所示，可以看出，本文算法重建的信号信噪比提升幅度较大，相对误差也较低，较优于其他算法.

Table 1 Performance Comparison under Different Algorithms表1 各算法恢复性能比较

实验2. 测试在不同数量测量值的情况下，各算法的性能.信号长度n=256，稀疏度k=20，测量值数量分别是100，128，150，192时，测试各算法的恢复效果.结果如图1～3所示.

图1是信噪比随测量值数量的变化曲线图.测试结果表明：在不同压缩比下，本文算法所重建信号的信噪比均高于对比算法，信号的精确度有较大的提升.

图2是相对误差随测量值数量变化的曲线图.由图2可以看出，在压缩比较低时，各算法的重建误差都比较大，但是随着测量值数量的增加，本文算法重建信号的相对误差比其他算法都有所降低，优势进一步体现.

Fig. 1 SNR comparison under different number of measured value图1 不同测量数值各算法重构信号的信噪比

Fig. 2 Relative error comparison under different number of measured value图2 不同测量数值各算法重构信号的相对误差

Fig. 3 Time comparison under different number of measured value图3 不同测量数值各算法重构时间

图3则展示了各算法的重构时间，在各压缩比下，本文算法的运算时间要少于NSL0算法和SL0算法，和AL0算法相差不大.实验3. 测试算法在不同稀疏度情况下的恢复精度.信号长度n=256，压缩比m/n=0.5，稀疏度由15～35变化，分别测试各算法的恢复效果.结果如图4所示，仿真结果表明：在不同稀疏度下，由本文算法所重建信号的信噪比均高于对比算法，尤其在稀疏度较低时，本文算法恢复的信号精度提高较大.

Fig. 4 SNR comparison under different sparseness图4 不同稀疏度下各算法重构信号的信噪比

实验4. 测试算法在不同噪声水平下的恢复精度.信号长度n=256，压缩比m/n=0.5，噪声均方差由0.004～0.02变化，分别测试各算法的恢复效果,结果如图5所示.可以看出，在不同噪声水平下，与对比算法相比，本文算法依然保持优越性，重建信号的信噪比仍提高较大.进一步说明在噪声干扰较小时，本文算法的性能要较好.

Fig. 5 SNR comparison under different noise levers图5 不同噪声均方差下各算法重构信号的信噪比

4 总 结

本文在之前l0范数工作的基础上，结合AL0算法快速收敛和牛顿迭代法高精度的优点，用一个简单的分式函数的和来近似估计l0范数，通过牛顿迭代算法求得稀疏解，从而对AL0算法恢复信号的精度不够高的缺陷得以改进.最后通过数值仿真，在恢复精度上和经典的OMP算法、相似的现有算法进行了比较.数值仿真结果表明：在不同压缩比、稀疏度及噪声干扰下，本文算法恢复精度均有较大提升，有效地改善了重建效果，提高了信号的重建质量.

[1]Candes E J, Wakin M B. An introduction to compressive sampling[J]. IEEE Signal Processing Magazine, 2008, 25(2): 21-30

[2]Tropp J A, Gilbert A C. Signal recover from random measurements via orthogonal matching pursuit[J]. IEEE Trans Information Theory, 2007, 53(12): 4655-4666

[3]Blumensath T, Davies M E. Gradient pursuits[J]. IEEE Trans on Signal Processing, 2008, 56(6): 2370-2382

[4]Dai Wei, Milenkovic O. Subspace pursuit for compressive sensing signal reconstruction[J]. IEEE Trans Information Theory, 2009, 55(5): 2230-2249

[5]Pei Tingrui, Yang Shu, Li Zhetao, et al. Detouring matching pursuit algorithm in compressed sensing[J]. Journal of Computer Research and Development, 2014, 51(9): 2101-2107 (in Chinese)

(裴廷睿, 杨术, 李哲涛, 等. 压缩感知中迂回式匹配追踪算法[J]. 计算机研究与发展, 2014, 51(9): 2101-2107)

[6]Chen S S B, Donoho D L, Saunders M A. Atomic decomposition by basis pursuit[J]. SIAM Review, 2001, 43(1): 129-159

[7]Mohimani H, Babaie-Zadeh M, Tutten C. A fast approach for overcomplete sparse decomposition based on smoothedl0norm[J]. IEEE Trans Information Theory, 2009, 57(1): 289-301

[8]Liu Wanjun, Wang Hongzhi, Wang Xianlong. Iterative shrinkage SL0 compressed sensing recovery algorithm[J]. Journal of Changchun University of Technology, 2014, 35(4): 434-437 (in Chinese)

(刘婉军, 王宏志, 王贤龙. 迭代收缩SL0压缩感知恢复算法[J]. 长春工业大学学报, 2014, 35(4): 434-437)

[9]Zhao Ruizhen, Lin Wanjuan, Li Hao, et al. Reconstruction algorithm for compressive sensing based on smoothedl0norm and revised Newton method[J]. Journal of Computer-Aided Design & Computer Graphics, 2012, 24(4): 478-484 (in Chinese)

(赵瑞珍, 林婉娟, 李浩, 等. 基于光滑l0范数和修正牛顿法的压缩感知重建算法[J]. 计算机辅助设计与图像学学报, 2012, 24(4): 478-484)

[10]Yang Lianglong, Zhao Shengmei, Zheng Baoyu, et al. The improved reconstruction algorithm for compressive sensing on SL0[J]. Signal Processing, 2012, 28(6): 834-841 (in Chinese)

(杨良龙, 赵生妹, 郑宝玉, 等. 基于SL0压缩感知信号重建的改进算法[J]. 信号处理, 2012, 28(6): 834-841)

[11]Wang Junhua, Huang Zhitao, Zhou Yiyu, et al. Robust sparse recovery based on approximatel0norm[J]. Acta Electronica Sinica, 2012, 46(6): 1185-1189 (in Chinese)

(王军华, 黄知涛, 周一宇, 等. 基于近似l0范数的稳健稀疏重构算法[J]. 电子学报, 2012, 40(6): 1185-1189)

[12]Xiang Gao, Zhang Xiaoling, Shi Jun. Complex-valued sparse reconstruction via arctangent regularization[J]. Signal Processing, 2014, 104(6): 450-463

[13]Li Ying, Wang Ze, Wang Junhua, et al. Sparse signal reconstruction based onl0norm approximation minimization[J]. Computer Engineering and Application, 2015, 51(10): 200-204 (in Chinese)

(李颖, 王泽, 王军华, 等. 基于l0范数近似最小化的稀疏信号重构方法[J]. 计算机工程与应用, 2015, 51(10): 200-204)

[14]Qi Huanfang, Xu Yuanhao. Improved SL0 algorithm for compressive sensing signal reconstruction[J]. Electronic Science & Technology, 2015, 28(4): 27-30 (in Chinese)

(齐焕芳, 徐源浩. 用于压缩感知信号重建的SL0改进算法[J]. 电子科技学报, 2015, 28(4): 27-30)

[15]Wu Feiyun, Zhou Yuehai, Tong Feng. A fast sparse signal recovery algorithm based on approximatel0norm and hybrid optimization[J]. Acta Automatica Sinica, 2014, 40(10): 2145-2150 (in Chinese)

(伍飞云, 周跃海, 童峰. 基于似零范数和混合优化的压缩感知信号快速重构算法[J]. 自动化学报, 2014, 40(10): 2145-2150)

[16]Shi Jun, Ding Renhuan, Xiang Gao, et al. Complex-valued sparse recovery via double-threshold sigmoid penalty[J]. Signal Processing, 2015, 114: 231-244