聂栋栋 弓耀玲
(燕山大学理学院 河北秦皇岛 066004) (niedd@ysu.edu.cn)
自压缩感知(compressed sensing, CS)信号处理理论 [1] 提出以来，信号重构作为其关键的环节，一直是该理论的研究重点.信号重构的目标是根据非自适应线性测量的低维数据，准确而有效地恢复原始采样获取的高维数据.在压缩感知理论框架下，信号在特定基下或一定的变换下可稀疏表示，则信号重构问题可转化为稀疏信号的重构.如何求欠定方程组最稀疏的解是稀疏重构问题的本质，即 l 0 范数的最小化问题.
然而直接求解欠定方程组的最稀疏解是NP难问题，处理起来非常棘手.在实际中，贪婪算法是最常见的近似求解方法 [2-5] ，如正交匹配追踪法(orthogonal matching pursuit, OMP) [2] 、梯度追踪法(gradient pursuit, GP) [3] 、子空间追踪法(subspace pursuit, SP) [4] 等，其基本思想是通过逐步迭代找到未知信号的支撑集，原信号的估计值即可利用最小二乘法得到.贪婪类算法需要较少的计算量，但是通常估计精度不高且需要较多的观测值.
另一种思路是寻找新的模型替代 l 0 范数的求解.这种方法分为2步：
1) 需要找到合适的函数近似估计 l 0 范数，由于 l 0 范数是 x 的不连续函数，实际计算中不好处理，很多文献常用 l 1 范数来近似 l 0 范数，如基追踪法(basis pursuit, BP) [6] ，其依据是 l 1 范数最小化与 l 0 范数最小化在满足一定条件时具有相同的稀疏解，可以通过线性规划方法求解.然而， l 1 范数不能充分反映向量 x 的稀疏特性，重构时需要较多的测量值，计算量较大.高斯函数类也是近似 l 0 范数经典的函数之一，基于光滑 l 0 范数算法(smoothed l 0 norm, SL0) [7] 是用一个连续平滑函数的极限来近似 l 0 范数，如文献[7-8]采用
(1)
对 l 0 范数进行近似，其中 σ 为趋于0的正数.文献[9-10]则在式(1)的基础上进一步改进，采用式双曲正切函数
(2)
近似 l 0 范数，其中 σ 为很小的正数.在文献[11-13]中，用一组反正切函数
(3)
近似 l 0 范数，其中 σ 为接近于0的正常数.文献[14]是采用复合三角函数近似 l 0 范数.但以上函数表达式均稍复杂，使得迭代计算复杂度较高.文献[15]的快速稀疏信号恢复算法(AL0算法)提出用函数
(4)
近似 l 0 范数，其中 
受此启发，本文在保证精确度的条件下，对近似 l 0 范数的函数进一步简化，从而简化计算.
2) 对新的函数模型选择不同的线性规划或优化算法求解.对于算法实现，NSL0算法选用修正牛顿法求解最优化问题，需要对海森矩阵进行修正，构造一个可近似海森矩阵逆的正定对称阵作为迭代方向，一定程度上提高了收敛速度.AL0算法将问题分解为2个最小化问题分别求解，然后交替迭代，将解投影到可行域.本文选取高精度的牛顿迭代法，将问题转为无约束优化问题后，对无约束问题直接应用牛顿算法，迭代得到最优解，不需要进一步对解进行投影.而且得到的牛顿方向为下降方向，不需要进一步修正，从而减少了误差影响，保证了精度.
信号的恢复质量在算法的不断改进下有所提升，但仍存在存储量大、计算复杂度高、精确度不够高等问题，需要进一步研究和改进.本文在之前 l 0 范数工作的基础上，结合AL0算法快速收敛和牛顿迭代法高精度的优点，提出新的基于近似 l 0 范数的稀疏信号重构算法，用一个简单的分式函数的和来近似估计 l 0 范数，通过牛顿迭代算法求得稀疏解.最后，通过数值仿真和已有的相似算法及经典的OMP算法进行比较，并分析了本文算法恢复的信号在重构误差、信噪比方面的优势.
令 z ∈ 
n 是传统采样得到的信号，经过正交基或紧框架 ψ 变换后的系数是稀疏的，即 z = ψx ，其中 x ∈ n 且 
为稀疏度.在无噪声干扰的情形下，测量向量 y 可以表示为 y = Φψx ，其中 y ∈ m ， Φ 为 m × n 维测量矩阵， ψ 为 n × n 维基矩阵.令 A = Φψ ，则重构问题简化为在矩阵 A 和测量向量 y 已知的条件下，求方程组 y = Ax 最稀疏解的问题，即 l 0 范数的最小化问题：
(5)
在实际应用中往往会不可避免地受噪声的影响，测量值是不精确的，令 γ 表示噪声，则测量值进一步表示为 y = Ax + γ .
本文采用一个形式较为简单的分式函数
(6)
来近似 l 0 范数，其中 δ 是很小的正数，实验中参数 δ 可取为一组下降且趋于0的序列.显然有:
(7)
令:
(8)
根据 l 0 范数的定义可知，当 δ 无限趋近于0时， F ( x )的值就无限接近于向量 x 中非0元的个数，即 x 的 l 0 范数，故有:
(9)
因此,我们可以用式(8)来近似 l 0 范数进行迭代计算.式(6)的优势在于形式简单，尤其在算法迭代时简化了计算，同时又不降低精度.
本文将最优化问题式(5)等价变换为拉格朗日形式最小化问题：
(10)
其中, λ 是拉格朗日乘子.然后对式(10)进行求解.由2.1节可知，式(10)可转化为
(11)
式(11)是无约束的最优化问题，最优解 x * 满足优化条件:
x L ( x * ; λ )=0.
(12)
进一步有:
(13)
其中, f ′( x )为 f ( x )的导数， A T 为 A 的转置.
当 δ →0，
x L ( x , λ )= A T Ax - A T y + 
(14)
为方便表示，令:
X f ≜ 
(15)
H f = A T A + λ X f ,
(16)
其中, 
为对角矩阵，对角线上元素为 
则有:
x L ( x ; λ )= H f x - A T y .
(17)
由式(17)可以得到Hessian矩阵可表示为
(18)
显然，Hessian矩阵 H f 为正定矩阵，故得到的牛顿方向必然为下降方向，可以利用牛顿迭代 [16] 求解式(11)得：
(19)
式(19)即算法的迭代式， h x 为迭代步长，初始解 x 0 取最小二乘解.
鉴于求逆运算会随着矩阵维数的增高计算量迅速加大，为尽量降低算法的计算复杂度，考虑在算法迭代中设法降低维数.
由于稀疏化后的信号向量的大部分值为0或接近于0，故可以只迭代计算对应变量的关键元素，如变量中绝对值最大的 l 个分量，忽略0元素及接近于0的元素.
令：
(20)
其中，参数 α ∈[0,1)，通常可取 α =0.01.
x S =( x i ), i ∈ S ,
(21)
A S =( a i ), i ∈ S ,
(22)
其中， a i 为 A 的第 i 列，则算法迭代的 H f 可替换为
(23)
相应地:
(24)
其中， 
为 S 的补集，从而使得逆矩阵的求解维数下降，较好地减少了计算复杂度.
在算法实现过程中，直接计算式(19)时，涉及到矩阵和向量的乘积及矩阵的求逆，假设 A 是 n × n 矩阵， y 为 n ×1向量，则 A T y 的计算复杂度为 O ( n 2 )，对 H f 求逆的计算复杂度为 O ( n 3 ).
经过降维后，令 l 为式(20)中 S 的元素个数，则如式(23)，对 H f 求逆的计算复杂度变为 O ( l 3 )，而通常情况下， l ≪ n ，从而使得本文算法每次迭代的计算复杂度维持在 O ( n 2 )，这与SL0，NSL0，OMP，AL0算法的计算复杂度均是相当的.
本文算法的伪代码描述如下：
算法1 . 基于近似 l 0 范数的稀疏信号重构算法.
输入：矩阵 A 、测量向量 y 、正则化参数 λ 、参数 δ 、终止条件 δ min ；
输出：重构信号 x .
初始化：初始解 x 0 = A T ( AA T ) -1 y .
迭代：
步骤1. 令:
S ={ t :| x t |>0},
x S =( x i ), i ∈ S ,
A S =( a i ), i ∈ S ;
步骤2. 更新矩阵：
步骤3. 迭代 x ，
步骤4. 更新 δ ， δ = δ ×0.5；
步骤5. 判断算法是否终止，若 δ > δ min 不满足输出迭代的结果 x ，否则返回步骤1继续迭代.
在不考虑噪声存在的情况下，算法的理论证明如下.
定理1 . 给定矩阵 A =( A S A N )，向量 
其中 x S 指 x 的非0元素部分， A S 指 A 中与 x S 对应相乘的列， y = Ax ,
R ≜ 
如果 Q ≜ 
且 
是可逆的，则 H -1 A T y = x .
证明. 由矩阵逆的分块形式知：
(25)
又 
因此有:
(26)
将式(25)带入式(26)并经过推导可以得到:
证毕.
定理2 . 给定矩阵 A ，向量 x 及其支撑集 S ，对角阵 X f ≜ 
满足 
且 
则对于任何正则化常数 λ ，有：
( A T A + λ X f ) -1 A T y = x .
证明. 对于任意支撑集 S ,存在初等矩阵 P S 使得矩阵 A 和向量 x 满足定理1(记为 
和 
即 
则:
( A T A + λ X f ) -1 A T y = 
(27)
又 
有:
(28)
其中, μ 是一个对角矩阵，对角元素和 
中的0元素相对应， ≜ 
显然， 
假设 Q 和 可逆，则由定理1有:
( A T A + λ X f ) -1 A T y = 
证毕.
为验证本文算法的有效性，对一维随机稀疏信号进行实验，并与经典算法OMP及相似的AL0算法等进行了比较.下面所有数值仿真中，矩阵 A 取高斯随机矩阵, 测量值向量通过 y = Ax + γ 生成, 其中 γ 表示高斯白噪声.实验采用信噪比(signal noise radio, SNR )作为评估算法精度的性能指标，定义为
(29)
其中, x opt 表示对源信号 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 不同噪声均方差下各算法重构信号的信噪比
本文在之前 l 0 范数工作的基础上，结合AL0算法快速收敛和牛顿迭代法高精度的优点，用一个简单的分式函数的和来近似估计 l 0 范数，通过牛顿迭代算法求得稀疏解，从而对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 smoothed l 0 norm[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 smoothed l 0 norm and revised Newton method[J]. Journal of Computer-Aided Design & Computer Graphics, 2012, 24(4): 478-484 (in Chinese)
(赵瑞珍, 林婉娟, 李浩, 等. 基于光滑 l 0 范数和修正牛顿法的压缩感知重建算法[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 approximate l 0 norm[J]. Acta Electronica Sinica, 2012, 46(6): 1185-1189 (in Chinese)
(王军华, 黄知涛, 周一宇, 等. 基于近似 l 0 范数的稳健稀疏重构算法[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 on l 0 norm approximation minimization[J]. Computer Engineering and Application, 2015, 51(10): 200-204 (in Chinese)
(李颖, 王泽, 王军华, 等. 基于 l 0 范数近似最小化的稀疏信号重构方法[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 approximate l 0 norm 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
Nie Dongdong and Gong Yaoling
( College of Science , Yanshan University , Qinhuangdao , Hebei 066004)
Abstract The signal reconstruction algorithm is the key to compressed sensing. Signal reconstruction based on approximate l 0 norm chooses a continuous function to estimate l 0 norm, thus the minimization problem of l 0 norm is transformed into an optimization problem of a smooth function. It is critical for the signal reconstruction algorithm to select the appropriate smooth function and optimization algorithm. To improve the accuracy of the sparse signal recovered in the compression sense, the sum of a simple fractional function is proposed to approximate l 0 norm on the basis of previous work in the paper. Then the sparse solution of an unconstrained optimization problem of the function is solved by Newton iterative algorithm, which effectively integrated the advantages of the fast convergence of approximate l 0 norm algorithm and the high precision of Newton iteration algorithm. Thus, the minimization of l 0 norm is approximated smoothly and efficiently within less time. The performance of the proposed algorithm is tested and compared with some existing similar algorithms in the case of different compression ratio, sparseness and noise levels in the simulation experiments. Simulation results show that the performance of the proposed algorithm is better than the existing similar algorithms, and the precision of reconstructed signal is greatly improved, which improves the signal recovery quality in compressed sensing effectively under the same conditions.
Key words compressed sensing; signal reconstruction; l 0 norm; continuous function; Newton iteration
摘 要 信号重构算法是压缩感知的关键.基于近似 l 0 范数的信号重构选取一个连续函数近似估计 l 0 范数，从而将 l 0 范数最小化问题转化为平滑函数的优化问题.该算法的关键在于选择合适的平滑函数和优化算法.为了提高压缩感知中稀疏信号恢复的精度，在之前工作的基础上，提出用一个简单的分式函数的和来近似估计 l 0 范数.然后通过牛顿迭代算法求解该函数的无约束优化问题的稀疏解，整合了似零范数算法快速收敛和牛顿迭代法精度高的优点.这样就可以在较少的时间内平滑且有效地近似 l 0 范数的最小化问题.仿真实验测试了所提算法在不同的压缩比、稀疏度及噪声水平情况下的性能，并与现有的同类算法进行了比较.结果表明：所提算法比现有的同类算法性能更好，重建信号的精度有了较大的提升，这有效地提高了在同等条件下压缩感知信号的恢复质量.
关键词 压缩感知；信号恢复； l 0 范数；连续函数；牛顿迭代
中图法分类号 TP391; TN911.73
收稿日期 ：2016-11-15；
修回日期： 2017-09-25
基金项目 ：燕山大学基础研究专项课题(理工A类)(15LGA016)
This work was supported by the Special Project for Basic Research of Yanshan University (Science & Engineering A) (15LGA016).
DOI: 10.7544/issn1000-1239.2018.20160829
Nie Dongdong , born in 1977. PhD, associate professor. Her main research interests include digital image processing and machine learning.
Gong Yaoling , born in 1990. Master. Her main research interests include image processing and compressed sensing.