求解基于lp+l2范数问题的共轭梯度法
Conjugate Gradient Method for lp+l2 Norm Problems

作者: 詹佳明 , 乌彩英 :内蒙古大学数学科学学院,内蒙古 呼和浩特;

关键词: 稀疏信号恢复三项共轭梯度法lp+l2范数全局收敛性Sparse Signal Recovery Conjugate Gradient Method lp+l2 Norm Global Convergence

摘要:
本文针对稀疏信号的重构问题提出新的函数模型,通过光滑化绝对值函数光滑我们的模型,基于三项共轭梯度法对稀疏信号进行恢复,并试验了不同参数值对稀疏信号恢复效果的影响,数值实验表明本文算法的数值有效性。

Abstract: A new model is proposed for sparse signal reconstruction. We smooth our model by smoothing absolute value function. Then we use a new tri-term conjugate gradient method to restore sparse signal. The effect of different parameters on sparse signal recovery is tested. The numerical results show that our algorithm is efficient.

1. 引言

压缩感知技术在信号重构问题中有着广泛的应用,而稀疏性是压缩感知中的重要问题。2006年,E. Candes,J. Romberg和T. Tao在文献 [1] 中提出,若矩阵A满足RIP条件,则可以通过下述模型(1)精确恢复信号:

min x R n x 0 , s .t . A x = b . (1)

其中 A R m × n ( m < n ) 为观测矩阵, b R m 为原始信号, x 0 表示0范数,即稀疏信号x中非零元素的个

数。但该问题是NP-hard问题 [2] ,因此学者们考虑了如下的凸优化问题:

min x R n x 1 , s .t . A x = b . (2)

(2)称作1范数模型,在适当的假设下,由文献 [3] 的定理1.3可知,模型(2)可以较精确地恢复原始信号。

由于问题(2)的凸性,有许多有效算法可求解之,如基追踪法 [4] ,迭代阈值方法 [5] 。近年来,p ( 0 < p < 1 )拟范数模型受到学者们的青睐,因p范数较1范数更能得到稀疏解,尤其是徐宗本提出p = 0.5时具有较好的计算效果 [6] 。p范数模型为:

min x R n x p , s .t . A x = b . (3)

由于 0 < p < 1 ,故问题(3)为非凸问题,且与0范数同样是NP-hard问题。文献 [1] 中作者利用光滑逼近 | | 函数进而光滑逼近p范数。受到该文的启发,本文利用文献 [7] [8] 中的光滑逼近 | | 函数技术光滑逼近p范数,且结合p范数与2范数各自的优点,利用p范数与2范数加权的方法弥补由于p过小而引起的数值不稳定,从而应用共轭梯度法进行信号恢复,并在适当的假设下,证明了算法的全局收敛性,同时进行了数值测试,测试结果证明我们的方法是有效的。

2. 模型及其性质

我们考虑如下模型:

min x R n λ 1 x p + λ 2 x 2 , s .t . A x = b .

其中 λ 1 , λ 2 R 。该问题的正则化问题为:

min x R n λ 1 x p + λ 2 x 2 + 1 2 A x b 2 2 . (4)

其模型较p范数模型(3),加入了一项2范数项进行调节p过小时目标函数的非凸程度,以p = 0.3时为例,依次做 R = λ 2 / λ 1 = 0 , 1 , 10 , 100 时的范数图像,如图1图4所示:

Figure 1. R = 0 norm image

图1. R = 0范数图像

Figure 2. R = 1 norm image

图2. R = 1范数图像

Figure 3. R = 10 norm image

图3. R = 10范数图像

Figure 4. R = 100 norm image

图4. R = 100范数图像

可见通过适当调整 λ 1 , λ 2 取值在保证产生稀疏解的同时可以增强问题的凸性。由p范数的定义 x p = ( i = 1 n | x i | p ) 1 p 可知,由于 | x | 的非光滑性,导致 x p 是非光滑的,因而我们通过光滑绝对值函数 | x | x p 进行光滑逼近。文献 [8] 中对绝对值函数提出如下两个光滑函数:

φ 1 ( ε , t ) = { t , t ε 2 , t 2 ε + ε 4 , ε 2 < t < t , t ε 2 . ε 2 φ 2 ( ε , t ) = { t 2 2 ε , | t | ε , | t | ε 2 , | t | > ε . (5)

利用 φ i ( i = 1 , 2 ) ,我们得到的模型为:

min λ 1 ϕ i ( ε , x ) + λ 2 x 2 , μ + 1 2 A x b 2 2 . (6)

其中,

ϕ i ( x ) = [ j = 1 n φ i ( ε , x j ) p ] 1 p , i = 1 , 2 , (7)

x 2 , μ = x 1 2 + x 2 2 + + x n 2 + μ 2 , μ > 0.

为叙述方便,令:

G i ( x ) = λ 1 ϕ i ( ε , x ) + λ 2 x 2 , μ + 1 2 A x b 2 2 , ( i = 1 , 2 ) . (8)

计算函数 G i ( x ) 的梯度:

G i ( x ) = λ 1 ϕ i ( ε , x ) + λ 2 x 2 , μ + A T ( A x b ) . (9)

引理1 设D为 R n 中的紧集,若 h : D R 是D上的光滑函数,则h在D上满足Lipschitz条件,即:

h ( x ) h ( y ) L x y , x , y D .

其中 L > 0 称为Lipschitz常数。

证明:对 x , y D ,由 h ( x ) 的光滑性及Lagrange中值定理可知:

| h ( x ) h ( y ) | = | h T ( ξ ) ( x y ) | h ( ξ ) x y ,其中 ξ 位于 x , y 之间。

因D为紧集,则 L > 0 h ( ξ ) L ,从而 h ( x ) h ( y ) L x y

引理2设 D R n 为紧集, f , g : D R 在D上Lipschitz连续,则 f g 在D上Lipschitz连续。

证明:因D紧且f和g均在D上Lipschitz连续,则 M > 0 , L 1 > 0 , L 2 > 0 ,满足:

| f ( x ) | M , | g ( x ) | M , | f ( x ) f ( y ) | L 1 x y , | g ( x ) g ( y ) | L 2 x y , x , y D .

x , y D ,有:

| f ( x ) g ( x ) f ( y ) g ( y ) | = | f ( x ) g ( x ) f ( y ) g ( x ) + f ( y ) g ( x ) f ( y ) g ( y ) | | g ( x ) | | f ( x ) f ( y ) | + | f ( y ) | | g ( x ) g ( y ) | M L 1 x y + M L 2 x y = M ( L 1 + L 2 ) x y .

引理3 若 φ i ( i = 1 , 2 ) 由(5)定义,则 φ i ( ε , t ) 在紧集D上Lipschitz连续。

证明:对于 i = 1 的情况,我们有:

φ i ( ε , t ) = { sgn ( t ) , | t | ε 2 , 2 t ε , | t | < ε 2 .

ε 2 > t 1 > t 2 > ε 2 时,有:

| φ 1 ( ε , t 1 ) φ 1 ( ε , t 2 ) | | 2 t 1 ε 2 t 2 ε | = 2 ε | t 1 t 2 | .

t 1 > ε 2 > t 2 > ε 2 时,有:

| φ 1 ( ε , t 1 ) φ 1 ( ε , t 2 ) | = | 1 2 t 2 ε | < | 2 t 1 ε 2 t 2 ε | = 2 ε | t 1 t 2 | .

ε 2 > t 1 > ε 2 > t 2 时,有:

| φ 1 ( ε , t 1 ) φ 1 ( ε , t 2 ) | = | 2 t 1 ε + 1 | < | 2 t 1 ε 2 t 2 ε | = 2 ε | t 1 t 2 | .

t 1 > ε 2 > ε 2 > t 2 时,有:

| φ 1 ( ε , t 1 ) φ 1 ( ε , t 2 ) | = 2 = 2 ε ε 2 ε | t 1 t 2 | .

t 1 > t 2 > ε 2 ε 2 > t 1 > t 2 时,有:

| φ 1 ( ε , t 1 ) φ 1 ( ε , t 2 ) | = 0 2 ε | t 1 t 2 | .

L 1 = 2 ε ,同理可证存在常数 L 2 使得 φ 2 ( t ) 在紧集D上Lipschitz连续。

最终令 L m = max { L 1 , L 2 } ,我们得到:

| φ i ( ε , t 1 ) φ i ( ε , t 2 ) | L m | t 1 t 2 | , t 1 , t 2 D , i = 1 , 2.

引理4 函数 ϕ i ( ε , x ) 由(7)所定义,则 ϕ i ( ε , x ) 在紧集D上Lipschitz连续。

证明:计算得:

[ ϕ i ( ε , x ) ] j = ( j = 1 n φ i p ( ε , x j ) ) 1 p p φ i p 1 ( ε , x j ) φ ( ε , x j ) .

φ i 连续可微,且 ε 0 , φ i 0 ,我们有 ( j = 1 n φ i p ( ε , x j ) ) 1 p p φ i p 1 均连续可微。假设集合D有界,则D为紧集,由引理1知 ( j = 1 n φ i p ( ε , x j ) ) 1 p p φ i p 1 在紧集D上Lipschitz连续。

同时由引理3可知 φ ( ε , x j ) 在紧集D上Lipschitz连续,则由引理2可知:

( j = 1 n φ i p ( ε , x j ) ) 1 p p φ i p 1 ( ε , x j ) φ ( ε , x j ) 在紧集D上Lipschitz连续,记Lipschitz常数为 L D ,那么有:

ϕ i ( ε , x ) ϕ i ( ε , y ) j = 1 n [ ϕ i ( ε , x ) ] j [ ϕ i ( ε , y ) ] j n L D x y .

ϕ i ( ε , x ) 在紧集D上Lipschitz连续。

引理5 x 2 在紧集D上Lipschitz连续。

证明:计算得 x 2 = [ x 1 x 2 , x 2 x 2 , , x n x 2 ] ,由引理1可知 x j x 在水平集上光滑,则 x j x 在紧集D上Lipschitz连续,记Lipschitz常数为 L c ,则:

x 2 y 2 n L c x y .

引理6 目标函数梯度 G i ( x ) 由(9)所定义,则 G i ( x ) 在紧集D上Lipschitz连续。

证明:计算得 G i ( x ) = λ 1 φ i ( ε , x ) + λ 2 x 2 , μ + A T ( A x b ) ,则:

G i ( x ) G i ( y ) = λ 1 ϕ i ( ε , x ) + λ 2 x 2 , μ + A T ( A x b ) λ 1 ϕ i ( ε , y ) λ 2 y 2 , μ A T ( A y b ) λ 1 n L D x y + λ 2 n L C x y + A 2 2 x y = ( λ 1 n L D + λ 2 n L C + A 2 2 ) x y .

3. 基于非精确线搜索的三项共轭梯度法

求问题 min x R n f ( x ) 的共轭梯度法的具体迭代过程如下:

x k + 1 = x k + α k d k , d k = { g k , k = 0 , g k + β k 1 d k 1 , k > 0.

其中 x k 代表当前迭代点, d k 代表当前搜索方向, α k 代表搜索步长, g k = G ( x k ) ,不同的 β k 代表不同

的共轭梯度算法。Shouqiang Du和Miao Chen在文献 [9] 中提出了一种新型三项共轭梯度法。具体迭代过程如下:

d k = { g k , k = 0 , g k + β k 1 d k 1 θ k 1 y k 1 , k > 0.

其中新加入项 θ k 1 y k 1 可以每步迭代确保目标函数值严格下降,使计算效率更高.本文将利用该共轭梯度法解决信号恢复问题。

算法1 (基于非精确线搜索的共轭梯度法)

步0给出初始参数 μ , ε , η , δ , λ 1 , λ 2 , ρ , σ , x 0 , μ 0 , γ ,置 k : = 0

步1若满足终止准则 G i ( x ) < ε ,停止计算,输出结果 x k x * ;否则,转步2。

步2计算搜索方向:

d k = { g k , k = 0 , g k + β k 1 d k 1 θ k 1 y k 1 , k > 0.

其中,

y k = g k g k 1 , β k = g k T y k 1 η g k 1 T d k 1 + μ | g k T y k 1 | , θ k = g k T d k 1 η g k 1 T d k 1 + μ | g k T y k 1 | .

步3令 α k = max { ρ m : m = 0 , 1 , 2 , } ,同时满足:

G ( x k + α k d k ) > G ( x k ) 2 σ ( 1 γ μ 0 ) α k G ( x k ) .

步4令 x k + 1 = x k + α k d k ,计算 g k + 1 = G ( x k + 1 )

步5令 k : = k + 1 ,转步1。

假设1令 c > 0 为常数,假设水平集 L c = { x k | G ( x k ) c } 是紧集。

定理1 若 { x k } 是由算法1产生的序列,则:

lim k 0 G ( x k ) = 0.

证明:由文献 [9] 中定理6,本文引理6及假设1可得。

4. 数值实验

本节实验在Intel Core i7-6500U 2.50GHz CPU,8G RAM,Windows10 64位操作系统,MATLAB R2016a环境下进行,所有数值实验结果为测试十次取平均值。

我们对两个不同的绝对值近似函数 φ 1 φ 2 分别在信号维数n取2048,4096,8192三种情况下,测试随机信号在无噪声和存在噪声干扰时的恢复效果.通过对比信号恢复时间t,相对误差e和迭代步数k来比较范数p和R对算法对信号恢复效果的影响。参数取值如下:A为 m × n 的随机高斯矩阵,x为n维随机初始信号, x 0 = z e r o s ( n , 1 ) 为初始点, g s = n o r m n d ( 0 , 0.01 , n / 2 , 1 ) 是一个平均值是0,标准差是0.01的对观测向量b的噪声干扰,p取0.3,0.4,0.5,0.6,R取0,1,10,100 (当R = 0时代表2范数项为0,即目标函数为p范数模型),其他参数取值如下:

μ = 1 e 5 , ε = 1 e 5 , η = 1 , δ = 2 , λ = 1 e 6 A T b , ρ = 0.5 , σ = 0.002 , γ = 0.2.

1) 我们先对绝对值近似函数 φ 1 的信号恢复问题进行数值实验,结果如表1表6所示:

Table 1. Signal dimension n = 2048, no noise interference

表1. 信号维数n = 2048,无噪声干扰

Table 2. Signal dimension n = 2048, Gauss noise interference

表2. 信号维数n = 2048,高斯噪声干扰

Table 3. Signal dimension n = 4096, no noise interference

表3. 信号维数n = 4096,无噪声干扰

Table 4. Signal dimension n = 4096, Gauss noise interference

表4. 信号维数n = 4096,高斯噪声干扰

Table 5. Signal dimension n = 8192, no noise interference

表5. 信号维数n = 8192,无噪声干扰

Table 6. Signal dimension n = 8192, Gauss noise interference

表6. 信号维数n = 8192,高斯噪声干扰

2) 我们再对绝对值近似函数 φ 2 的信号恢复问题进行数值实验,结果如表7表12所示:

Table 7. Signal dimension n = 2048, no noise interference

表7. 信号维数n = 2048,无噪声干扰

Table 8. Signal dimension n = 2048, Gauss noise interference

表8. 信号维数n = 2048,高斯噪声干扰

Table 9. Signal dimension n = 4096, no noise interference

表9. 信号维数n = 4096,外界无噪声干扰

Table 10. Signal dimension n = 4096, Gauss noise interference

表10. 信号维数n = 4096,高斯噪声干扰

Table 11. Signal dimension n = 8192, no noise interference

表11. 信号维数n = 8192,无噪声干扰

Table 12. Signal dimension n = 8192, Gauss noise interference

表12. 信号维数n = 8192,高斯噪声干扰

5. 主要结论

通过在不同维度n = 2048,4096,8192下做信号恢复效果分析,我们主要得出以下结论:

1) 相对于p范数模型, l p + l 2 范数模型具有更好的恢复效果,二者在时间和迭代步数相近的情况下, l p + l 2 范数模型对初始信号的精度更高。

2) 做 λ 2 / λ 1 的值对信号恢复效果影响的横向对比, λ 2 / λ 1 的取值对信号恢复的时间和迭代步数影响并不明显;但是 λ 2 / λ 1 = 1 λ 2 / λ 1 = 10 时信号的恢复具有更高精度,这是由于当 λ 2 / λ 1 取值为1~10时2范

数项对p范数的调节使得函数模型的非凸性程度降低,因而提高信号恢复的精度。

3) 做p的值对信号恢复效果影响的纵向对比,显然范数p = 0.6时信号恢复的精度最高,p = 0.3时精度最低。

注:1) 文献 [10] 中提出的模型为

min x R n λ x 1 τ x 2 , s .t . A x = b .

即利用 l 1 l 2 范数加权从而产生稀疏阶,而本文考虑p范数比1范数更易得到稀疏解,从而利用 l p + l 2 范数加权进行信号恢复。

2) 文献 [10] 中的光滑化绝对值函数是文献 [8] 中 φ 3

3) 文献 [11] 中的模型为:

min x R n λ 1 x p p + λ 2 x 2 2 + 1 2 A x b 2 2 .

λ 1 x p 2 + λ 2 x 2 2 λ 1 x p + λ 2 x 2 不易产生稀疏解。取p = 0.3,R = 100对比图形如图5图6所示:

Figure 5. λ 1 x p p + λ 2 x 2 2 norm image

图5. λ 1 x p p + λ 2 x 2 2 范数图像

Figure 6. λ 1 x p + λ 2 x 2 norm image

图6. λ 1 x p + λ 2 x 2 范数图像

基金项目

内蒙古自然科学基金(2018MS01016)。

NOTES

*通讯作者。

文章引用: 詹佳明 , 乌彩英 (2019) 求解基于lp+l2范数问题的共轭梯度法。 应用数学进展, 8, 512-523. doi: 10.12677/AAM.2019.83057

参考文献

[1] Candes, E.J., Romberg, J. and Tao, T. (2006) Robust Uncertainty Principles: Exact Signal Reconstruction from Highly Incomplete Frequency Information. IEEE Transactions on Information Theory, 52, 489-509.
https://doi.org/10.1109/TIT.2005.862083

[2] Ge, D., Jiang, X. and Ye, Y. (2011) A Note on the Complexity of Lp Minimization. Mathematical Programming, 129, 285-299.
https://doi.org/10.1007/s10107-011-0470-2

[3] Candes, E.J. and Tao, T. (2004) Near Optimal Signal Recovery from Random Projections: Universal Encoding Strategies. IEEE Transactions on Information Theory, 52, 5406-5425.
https://doi.org/10.1109/TIT.2006.885507

[4] Chen, S.S., Donoho, D.L. and Saunders, M.A. (1998) Atomic De-composition by Basis Pursuit. SIAM Journal on Scientific Computing, 20, 33-61.
https://doi.org/10.1137/S1064827596304010

[5] Daubechies, I., Defries, M. and De Mol, C. (2004) An Iterative Thresholding Algorithm for Linear Inverse Problems with a Sparsity Constraint. Communications on Pure and Applied Mathematics, 57, 1413-1457.
https://doi.org/10.1002/cpa.20042

[6] Xu, Z.B., Chang, X., Xu, F., et al. (2012) Regularization: A Thresholding Representation Theory and a Fast Solver. IEEE Transactions on Neural Networks and Learning Systems, 23, 1013-1027.
https://doi.org/10.1109/TNNLS.2012.2197412

[7] Jeevan, K., Pant and Lu, W. (2014) New Improved Algorithms for Compressive Sensing Based on l(p) Norm. Circuits and Systems II: Express Briefs. IEEE Transactions, 61, 198-202.

[8] Saheya, B., Yu, C.-H. and Chen, J.-S. (2018) Numerical Comparisons Based on Four Smoothing Func-tions for Absolute Value Equation. Journal of Applied Mathematics and Computing, 56, 131-149.
https://doi.org/10.1007/s12190-016-1065-0

[9] Bakhtawar, B., Zabidin, S., Ahmad, A., et al. (2017) A New Modified Three-Term Conjugate Gradient Method with Sufficient Descent Property and Its Global Convergence. Journal of Mathematics, 2017, 1-12.

[10] Liu, C., Chen, Q., Zhou, B., et al. (2016) L1- and L2-Norm Joint Regularization Based Sparse Signal Reconstruction Scheme. Applied Physics A, 45, 313-323.

[11] Zhang, Y. and Ye, W.Z. (2017) Sparse Recovery by the Iteratively Reweighted, Algorithm for Elastic, Minimization. Optimization, 66, 1-11.
https://doi.org/10.1080/02331934.2017.1359590

分享
Top