原文链接：Thresholded Smoothed l0 Norm for Accelerated Sparse Recovery http://ieeexplore.ieee.org/document/7069222/

带阈值的平滑l0范数加速稀疏恢复

Han Wang, Qing Guo, Member, IEEE, Gengxin Zhang, Guangxia Li, and Wei Xiang, Senior Member, IEEE

译者：柳如风

 

摘要：平滑l0范数(Smoothed l0 norm,SL0)是一种快速并可扩展到复频域的稀疏恢复算法，它适用于许多实际的实时应用。在这封信中，我们提出了一种增强的算法名为“带阈值的平滑l0范数(Thresholded Smoothed l0 Norm,T-SL0)”来加速SL0的迭代过程。T-SL0引入一个迭代效率指示量并和预置的阈值进行比较来实时决定是否要继续迭代。通过确定和避开低效率的迭代，我们的方法比原始SL0算法收敛更快。展示的实验结果表明我们的方法在不损失精度的情况下能有效地加速SL0算法。

关键词：压缩感知，稀疏恢复，平滑l0范数

I. 介绍

　　稀疏性是许多自然信号固有的属性，从医学图像到天文学数据。得益于压缩感知理论[1]-[3]，这些稀疏信号能够高效地压缩然后通过许多种稀疏恢复算法来恢复。

　　平滑l0范数(SL0)[4]是一种对稀疏恢复有效的启发式算法。不像大多数使用凸松弛或者贪婪基追踪算法，SL0利用了分层次的非凸性(Graduated Non-Convexity,GNC)[5]方法，它是著名的Simulated Annealing[6]方法的确定性版本，为了寻找最稀疏的解去优化非凸的l0范数。与基于l1范数的比较(FISTA[7]或SPGL1[8]等)，由于l0范数的目标函数，SL0能利用弱稀疏性约束在噪声环境中稳定地恢复[9]。得益于它已被证实的快速和复频域可扩展性[10]，SL0能被应用到实时信号处理中，例如视频重构核信道估计。然而，由于不可弯曲的策略（inflexible,(fixed-count loop)），可以看出这算法在一些迭代中收敛缓慢，而在其他的迭代很快速，它表示SL0中的一些迭代是低效的。观测到的这个促使我们通过避开无效的迭代来加速算法。

　　为了加速SL0，我们首先分析了迭代中的变量依赖关系，并利用一个变量来表示实时的收敛效率。然后我们比较指示量和阈值，如果指示量的绝对值大于阈值，当前的迭代可以认为是高效的然后执行它。否则，它将被忽略。由于引入的阈值，增强的SL0算法相比原始的被极大地加速了。

　　基于以上想法，这封信提出了一种增强的算法名为“带阈值的平滑l0范数(T-SL0)”。实验结果显示我们的方法能够在不损失精度的前提下两倍加速算法。算法剩下的部分在下面组织。在第二段，原始SL0算法被简单地回顾。然后，在第三段，我们分析了迭代过程中的变量依赖关系，介绍了有条件迭代的基本概念并提出了T-SL0算法。最后，实验结果被展示在第四段中来评价我们的方法的性能。

II. 平滑l0范数算法

　　基本的想法在CS理论里是从压缩版本的y等中提取出一个系数的向量x，为了寻找问题${P_0}:{\min _x}{\left\| x \right\|_0},s.t.y = Ax$的一个解，其中${l_0}$范数${\left\| x \right\|_0}$表示x中非零元素的数目。

　　因为目标函数是非凸的（从而棘手），SL0算法旨在解决平滑版本${P_0}$的问题。SL0使用凸函数$N - {F_\sigma }(x)$近似${\left\| x \right\|_0}$，其中N是向量x的维度，${F_\sigma }(x) = \sum\nolimits_{n = 1}^N {{f_\sigma }\left( {{x_n}} \right)} $，${f_\sigma }$属于一簇当形状参数$\sigma  \to 0$时近似于克罗内克函数(${\delta _{ij}} = \left\{ {\begin{array}{*{20}{c}}
{0{\rm{, }}if{\rm{ }}i \ne j{\rm{ }}}\\
{1{\rm{,  }}if{\rm{ }}i = j}
\end{array}} \right.$)的凸函数。从而，原始问题${P_0}$等价于问题$Q:{\lim _{\sigma  \to 0}}{\max _x}{F_\sigma }(x),s.t.y = Ax$。

　　然而，对于小的$\sigma $，函数${F_\sigma }$包含了许多局部最大值且它的最大化不容易，相反，当$\sigma $足够大时它没有局部最大值。从而，SL0利用GNC，一种启发式方法能够通过逐渐收缩凸平滑函数的形状参数来最优化非凸函数，来避免${F_\sigma }$陷入局部最大值。特别地，不是问题Q，SL0解一系列的问题${Q_\sigma }:{\max _x}{F_\sigma }(x),s.t.y = Ax$，通过每次循环逐渐减少$\sigma $。每个子问题${Q_\sigma }$使用梯度上升法迭代几次就能解出。从而，SL0算法包含了两步循环：外循环逐渐减少$\sigma $使${f_\sigma }$趋向于克罗内克函数，同时内循环对给定的$\sigma $使用梯度上升法通过很少的迭代次数解${Q_\sigma }$，避免陷入局部最大值。SL0的实现算法通过使用高斯簇函数${f_\sigma }(x) = {e^{ - {{\left| x \right|}^2}/2{\sigma ^2}}}$在图1中给出。

图1

III. 带阈值的平滑l0范数算法

A. 分析变量依赖关系

　　为了揭露SL0迭代过程中的变量依赖关系，首先定义$g(x) \buildrel \Delta \over =  - x{e^{ - {{(\left| x \right|/\sqrt 2 {\sigma ^{(j)}})}^2}}}$，和它的偏导$\partial g(x)/\partial \bar x = \left( {{{\left( {\left| x \right|/{\sigma ^{(j)}}} \right)}^2} - 1} \right){e^{ - {{\left( {\left| x \right|/\sqrt 2 {\sigma ^{\left( j \right)}}} \right)}^2}}}$。然后在第i次迭代中，步进向量的第n个元素$\Delta {x^{(i)}}$能被展开为$\Delta x_n^{(i)} = g\left( {x_n^{(i - 1)}} \right)$。此外，在一下分析中，我们假设源的稀疏度不超过理论边界[9]保证稳定地恢复。

　　让$I = LJ$是总的迭代次数，并使用拉格朗日均值定理，有

$\Delta x_n^{(i)} - \Delta x_n^{(I)} = \frac{{\partial g}}{{\partial \bar x}}\left( {\xi _n^{(i)}} \right)\left( {x_n^{(i - 1)} - x_n^{(I - 1)}} \right)$ ,                                                            (1)

其中$\xi _n^{(i)}$落在$x_n^{(i - 1)}$和$x_n^{(I - 1)}$之间。当迭代终止时，有$\Delta x_n^{(I)} \to 0$和$\Delta x_n^{(I - 1)} \to {x_n}$，其中${x_n}$是x的第n个元素。从而，(1)能被如下展开

$\Delta {x^{(i)}} = {D^{(i)}}\left( {{x^{(i - 1)}} - x} \right)$,                                                         (2)

其中${D^{(i)}} = diag\left( {\frac{{\partial g}}{{\partial \bar x}}\left( {\xi _1^{(i)}} \right),...,\frac{{\partial g}}{{\partial \bar x}}\left( {\xi _N^{(i)}} \right)} \right)$是对角矩阵，它的逆是${\left( {{D^{(i)}}} \right)^{ - 1}}$。从而：

${x^{(i - 1)}} - x = {\left( {{D^{(i)}}} \right)^{ - 1}}\Delta {x^{(i)}}$,                                                            (3)

对(2)和(3)的两边取l2范数得到

$\left\{ {\begin{array}{*{20}{c}}
{{{\left\| {\Delta {x^{(i)}}} \right\|}_2} \le \mathop {\max }\limits_{1 \le n \le N} \left\{ {\left| {\frac{{\partial g}}{{\partial \bar x}}\left( {\xi _n^{(i)}} \right)} \right|} \right\} \cdot {{\left\| {{x^{(i - 1)}} - x} \right\|}_2}}\\
{{{\left\| {{x^{(i - 1)}} - x} \right\|}_2} \le \mathop {\max }\limits_{1 \le n \le N} \left\{ {\left| {\frac{{\partial g}}{{\partial \bar x}}{{\left( {\xi _n^{(i)}} \right)}^{ - 1}}} \right|} \right\} \cdot {{\left\| {\Delta {x^{(i)}}} \right\|}_2}}
\end{array}} \right.$,                                                             (4)

很容易看出

$0 \le \frac{{{{\left\| {\Delta {x^{(i)}}} \right\|}_2}}}{{{{\left\| {{x^{(i - 1)}} - x} \right\|}_2}}} \le 1$,                                                            (5)

给定$RMSE = {\left\| {\hat x - x} \right\|_2}/\sqrt N $，(5)式写为

${\left\| {\Delta {x^{(i)}}} \right\|_2} = {c^{(i)}}RMS{E^{(i - 1)}}$,                                                            (6)

其中尺度因子${c^{(i)}} \in \left[ {0,\sqrt N } \right],i = 1,...,I$。因为GNC基于迭代是准稳态的，被恢复的向量 ${x^{(i - 1)}}$，连同对应的$RMS{E^{(i - 1)}}$和步长向量$\Delta {x^{(i)}}$哪一个是被${x^{(i - 1)}}$确定的，在连续迭代中微小地改变[4]。因此尺度因子${c^{(i)}}$也在当前迭代和前一次迭代中微小地改变，即$r_c^{(i)} = {c^{(i)}}/{c^{(i - 1)}} \approx 1$。

　　${c^{(i)}}$的稳定属性暗示我们能够越过(6)式宽松的边界并预测${\left\| {\Delta {x^{(i)}}} \right\|_2}$，第i次迭代的步长向量的欧几里得长度，是与第i-1次的RMSE成比例。由于x未知从而RMSE不能再迭代过程中得到，步长大小${\left\| {\Delta x} \right\|_2}$能够用作指示量来反映解的实时RMSE等级。

B. 我们方法的主要思想

　　在SL0算法中，恢复向量随着每次迭代逐渐减少RMSE逐渐接近它的真实值x。然而，需要注意的是在图1的内循环中无论效率总是执行L次。例如不可调节策略可能导致一些低效的等级。我们大量仿真结果显示了RMSE在一些迭代中缓慢减少而另一些快速减少，因此迭代效率在每次迭代改变。鉴于这种观察背后的基本思想提出T-SL0算法引入迭代阈值以便只执行高效的迭代同时避开低效的。

　　特别地，我们考虑RMSE的相对改变，这反映了逐次迭代的效率，如下：

$\begin{array}{l}
R_{RMSE}^{(i - 1)} = \frac{{RMS{E^{(i - 1)}} - RMS{E^{(i - 2)}}}}{{RMS{E^{(i - 2)}}}}\\
{\rm{         }} = \frac{{\frac{1}{{r_c^{(i)}}}{{\left\| {\Delta {x^{(i)}}} \right\|}_2} - {{\left\| {\Delta {x^{(i - 1)}}} \right\|}_2}}}{{{{\left\| {\Delta {x^{(i - 1)}}} \right\|}_2}}}\\
{\rm{         }} \approx \frac{{{{\left\| {\Delta {x^{(i)}}} \right\|}_2} - {{\left\| {\Delta {x^{(i - 1)}}} \right\|}_2}}}{{{{\left\| {\Delta {x^{(i - 1)}}} \right\|}_2}}} = R_{step}^{(i)}
\end{array}$,                                                            (7)

其中$R_{step}^{(i)}$是步长的相对改变。等式(7)意味着${R_{step}}$能被视为实时收敛效率的指示。

图2

　　和图1的原始SL0算法比较，提出的算法通过(7)式计算$R_{step}^{(i)}$并将它和预置的阈值$\lambda $比较。如果$\left| {R_{step}^{(i)}} \right| \ge \lambda $，由于GNC的准稳态性质，这意味着RMSE在第(i-1)次迭代快速下降并且可能在第i次迭代也一样。否则，如果$\left| {R_{step}^{(i)}} \right| < \lambda $，在当前循环的第i次和下一次迭代可以看作低效的，因此避开它进入下一次循环。

　　提出的算法T-SL0的细节在图2中。我们的方法分享和图1中SL0相同的框架，除了在内循环中使用了附加的评价步骤。可以看到计算$R_{step}^{(i)}$引入了很少的计算复杂性，同时新的算法通过避开无意义的迭代获得了非凡的速度增加。

　　以上的迭代效率阈值$\lambda $决定了多少迭代将被执行。如果$\lambda $选的太小或太大，迭代数也会剩余或者不足。因此，阈值应该合适选择。在大多数应用中，例如视频重构或者信道估计中，恢复算法主要是一帧一帧执行的。源的稀疏性（视频采样或者信道脉冲响应）关于帧长度缓慢改变，从而迭代处理在连续帧缓慢改变。所以，一个敏感的阈值确定方式是从一些帧选择“平均迭代效率”作为阈值。

　　特别地，我们能够从恢复第一帧过程中训练SL0算法，然后获取$R_{step}^{(i)},i = 1,...,I$的序列。使

$\lambda  = \left| {\frac{1}{I}\sum\limits_{i = 1}^I {R_{step}^{(i)}} } \right|$,                                                            (8)

作为第一帧的平均迭代效率。取$\lambda $作为T-SL0的阈值，大约一半的迭代次数可视为无效的并能被忽略。从而，稀疏信号恢复的整个迭代过程能被加速约两倍。

图3

IV. 实验结果

　　在这段中，我们将验证以上的分析并通过实验比较T-SL0和SL0的性能和恢复能力。在我们的实验中，源稀疏向量x是使用复伯努利-高斯模型产生，即${x_n}\~pCN(0,{\sigma _s}) + (1 - p)CN(0,{\sigma _n})$，其中$CN(u,\sigma )$表示均值为$\mu $、标准差为$\sigma $的复高斯分布，$p$表示${x_n}$是标准差为${\sigma _s}$的复信号的概率，同时$(1 - p)$表示${x_n}$是标准差为${\sigma _n}$的复高斯噪声采样的概率。稀疏性意味着$p \ll 1$。令$N = 1000,p = 0.1,{\sigma _s} = 1,{\sigma _n} = {10^{ - 3}}$，从而稀疏源x是K稀疏、N维的带复高斯噪声的复列向量，其中K=kN=100。测量能够实现为${y^{M \times 1}} = {A^{M \times N}}{x^{N \times 1}}$，其中A是实际的傅里叶转换矩阵并且M=400是测量数量。SL0和T-SL0的参数选择为：$L = 3,q = 0.95,{\mu _0} = 2$和${\sigma ^{(J)}} = 4{\sigma _n} = 4 \times {10^{ - 3}}$。我们的实验运行在MATLAB R2014b，使用Intel Core i7 4770K 3.90 GHz处理器和8G内存。

　　图3画出了SL0迭代过程中多个变量的值。很明显步长尺寸${\left\| {\Delta {x^{(i)}}} \right\|_2}$和对应的$RMS{E^{(i - 1)}}$同时下降导致尺度因子${c^{(i)}} = {\left\| {\Delta {x^{(i)}}} \right\|_2}/RMS{E^{(i - 1)}}$维持稳定，从而使$r_c^{(i)} = {c^{(i)}}/{c^{(i - 1)}} \approx 1$ 。概括地，${\left\| {\Delta x} \right\|_2}$和$RMSE$相对改变，即$R_{step}^{(i)}$和$R_{step}^{(i - 1)}$，互相配合良好。这些结果确认了第III段的(6)和(7)式，并且透露了收敛效率指示${R_{step}}$确实反映了实时迭代效率。

　　在下面的实验，我们首先用SL0恢复了稀疏向量x，确定迭代效率阈值$\lambda $，然后再使用T-SL0恢复x。SL0和T-SL0的迭代恢复处理画在图4。其中显示两个算法都随迭代次数收敛。更大的$\left| {{R_{step}}} \right|$有更快的RMSE下降。然而，由于使用阈值$\lambda $，T-SL0只在高的$\left| {{R_{step}}} \right|$时执行。特别在实验中，T-SL0仅在168次迭代(14.9ms)时终止而不是SL0需要的354次迭代(29.1ms)，同时产生了与SL0同等级的RMSE为。

　　图5比较了T-SL0算法在不同阈值的RMSE下降，其中标准阈值从(8)式得来，并且对原始SL0算法“没有阈值”参考。图5中的实验结果显示了更低的阈值可能导致过多的迭代次数（能够从末端看到），同时高的阈值导致迭代数目不足以接近真实的RMSE。因此，(8)式中给出的“平均迭代效率”是阈值的合理选择。

　　为了更进一步评价恢复能力，我们在图6中比较T-SLO，SL0，FISTA和SPGL1的相位转变。相空间$\left( {\delta ,\rho } \right) \in {\left[ {0,1} \right]^2}$被分为均匀分布的网格，其中$\delta  = M/N$、$\rho  = K/M$是归一化的不确定测量值和稀疏度，分别在每个网格点，求解实现了${10^4}$次独立的稀疏恢复，且平均的信号运行时间是17.8ms(T-SL0),34.9ms(SL0),87.2ms(FISTA),63.5ms(SPGL1)。基于以上重复试验，一条相转变曲线分割空间为两部分。也就是，右下部分是统计上可恢复的，左上部分是不可恢复的。试验结果显示了T-SL0的恢复能力与SL0有细微的不同，且在同阶相同大小计算负载下要优于基于FISTA和SPGL1的l1范数求解法。

　　从而，以上的实验结果清晰地表明了这些，由于引入了阈值$\lambda $，T-SL0能够增加SL0的执行速度两倍且不损失精度。

图4

图5

图6

IV. 结论

　　在这封信中，我们提出了一种对实时稀疏信号恢复的带阈值的平滑l0范数算法。基于变量依赖关系的分析结果，我们引入了一个指示量和一个阈值来评价实时迭代效率。通过比较指示量和阈值，T-SL0只执行高效的迭代，因此收敛速度比原始的SL0算法更快。我们的方法共享了与SL0相同的框架，除了增加在迭代过程中的评价处理。通过引入很少的复杂性同时通过避开无意义的迭代大大地加速了算法。实验结果也显示了我们的方法能够在不损失精度的情况下显著地加速SL0算法。

参考文献

[1] E. Candes, J. Romberg, and T. Tao, “Robust uncertainty principles: Exact

signal reconstruction from highly incomplete frequency information,”

IEEE Trans. Inf. Theory, vol. 52, no. 2, pp. 489–509, Feb. 2006.

[2] E. Candes and T. Tao, “Near-optimal signal recovery from random

projections: Universal encoding strategies?” IEEE Trans. Inf. Theory,

vol. 52, no. 12, pp. 5406–5425, Dec. 2006.

[3] D. Donoho, “Compressed sensing,” IEEE Trans. Inf. Theory, vol. 52,

no. 4, pp. 1289–1306, Apr. 2006.

[4] H. Mohimani, M. Babaie-Zadeh, and C. Jutten, “A fast approach for

overcomplete sparse decomposition based on smoothed _0 norm,” IEEE

Trans. Signal Process., vol. 57, no. 1, pp. 289–301, Jan. 2009.

[5] A. Blake and A. Zisserman, Visual Reconstruction. Cambridge, MA,

USA: MIT Press, 1987, vol. 2.

[6] S. Kirkpatrick, C. D. Gelatt, andM. P. Vecchi, “Optimization by simulated

annealing,” Science, vol. 220, no. 4598, pp. 671–680, 1983.

[7] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding algorithm

for linear inverse problems,” SIAM J. Imaging Sci., vol. 2, no. 1,

pp. 183–202, 2009.

[8] E. van den Berg and M. P. Friedlander, “Probing the Pareto frontier for

basis pursuit solutions,” SIAM J. Sci. Comput., vol. 31, no. 2, pp. 890–912,

2008.

[9] M. Babaie-Zadeh and C. Jutten, “On the stable recovery of the sparsest

overcomplete representations in presence of noise,” IEEE Trans. Signal

Process., vol. 58, no. 10, pp. 5396–5400, Oct. 2010.

[10] G. H. Mohimani, M. Babaie-Zadeh, and C. Jutten, “Complex-valued

sparse representation based on smoothed _0 norm,” in Proc. IEEE

ICASSP, Las Vegas, NV, USA, Mar. 2008, pp. 3881–3884.