# A New Smoothed L0 Regularization Approach for Sparse Signal Recovery

Authors: Jianhong Xiang; Huihui Yue; Xiangjun Yin; Linyu Wang
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Sparse signal reconstruction, as the main link of compressive sensing (CS) theory, has attracted extensive attention in recent years. The essence of sparse signal reconstruction is how to recover the original signal accurately and effectively from an underdetermined linear system equation (ULSE). For this problem, we propose a new algorithm called regularization reweighted smoothedL0 norm minimization algorithm, which is simply called RRSL0 algorithm. Three innovations are made under the framework of this method: (1) a new smoothed function called compound inverse proportional function (CIPF) is proposed; (2) a new reweighted function is proposed; and (3) a mixed conjugate gradient (MCG) method is proposed. In this algorithm, the reweighted function and the new smoothed function are combined as the sparsity promoting objective, and the constraint condition y-Φx22 is taken as a deviation term. Both of them constitute an unconstrained optimization problem under the Tikhonov regularization criterion and the MCG method constructed is used to optimize the problem and realize high-precision reconstruction of sparse signals under noise conditions. Sparse signal recovery experiments on both the simulated and real data show the proposed RRSL0 algorithm performs better than other popular approaches and achieves state-of-the-art performances in signal and image processing.


---

## Body

## 1. Introduction

CS [1, 2] has been successfully applied in a multitude of scientific fields, ranging from image processing tasks to radar to coding theory, making the potential impact of advancements in theory and practice rather large. In fact, various CS tasks eventually boil down to the sparse signal recovery problem in the following underdetermined linear system: (1)y=Φx+w,where Φ∈Rm×n is a sensing matrix (also known as measurement matrix). x∈Rn is the sparse vector (or signal) to be solved. y∈Rm is the vector (or signal) of measurements. Moreover, the w∈Rm denotes the additive noise.For solving the ULSE in (1), we try to recover the sparse signal x from given {y,Φ}. In this case, the sensing matrix Φ contains more columns than rows, which means there would be more than one solution that satisfies the constraint Φx-y22≤ϵ (ϵ is a very small constant). This makes the recovery of sparse signal x an ill-posed problem. Luckily, since the target signal itself is sparse, the most straightforward method is to use its sparsity to improve the problem. That is, the problem of recovering the target signal x can be converted into solving the L0-norm minimization problem as follows:(2)P0argminx∈Rnx0,s.t.Φx-y22≤ϵ.This rather wonderful transformation is actually supported by outstanding theory [3]. The theory points out that, under the condition of no noise, the sparsest solution is indeed a real signal if x is sufficiently sparse and Φ satisfies the restricted isometry property (RIP) [4]:(3)1-δKx22≤Φx22≤1+δKx22,where K is the sparsity of signal x. In (2), L0-norm is nonconvex, which leads to a NP-hard problem. In fact, there are two alternative approaches to solve this problem [5, 6]: (i)
Greedy search(ii)
Relaxation method for theP0Greedy search requires known sparsity as a constraint and the main methods are the approximate algorithms based on greedy matching pursuit (GMP) algorithms, such as OMP [7], StOMP [8], ROMP [9], CoSaMP [10], GOMP [11], and SP [12] algorithms. The objective function of these algorithms is given by [6](4)Pfargminx∈Rn12Φx-y22,s.t.x0≤K.where K is the sparsity of signal x as shown in (3).Based on this, the characteristics of GMP algorithms can be summarized as follows [6]:(i)
Sparsity is used as a prior information(ii)
Least square error is employed as an iterative criterionThe main advantage of GMP algorithms is that it is simple to calculate, but its application range is limited due to its low reconstruction accuracy in the case of noise.At present, relaxation method for theP0 is a main method. Relaxation method is mainly divided into two categories: L1-norm minimization methods and smoothed L0-norm minimization methods. The representative algorithm of the former is BP algorithm [13], and the latter is SL0 algorithm. Both are optimized using an approximate or equivalent L0-norm objective function. Therefore, the effect of sparse signal recovery is similar to that of GMP algorithms.The above sparse signal recovery methods focus more on signal recovery in the absence of noise; they work not well in noise case. However, sparse signal recovery under noise is a very realistic and inevitable problem. Fortunately, the regularization mechanism makes it possible to solve this problem. The regularization mechanism relaxesP0 into the following unconstrained recovery problem [6]: (5)Pgargminx∈Rn12Φx-y22+λgx,where λ>0 is the parameter that balances the trade-off between the deviation term Φx-y22 and the sparsity regularizer g(x). The sparse prior information is enforced by the regularizer g(x). Proper g(x) is crucial to the successful recovery of sparse signal: it should support sparse solutions and ensure that problem Pg can be effectively solved simultaneously.For regularization, there are many sparsity regularizers for relaxing theL0-norm, among which the convex L1-norm [14, 15] and the nonconvex Lp-norm to the p-th power [16–18] are the main methods:(i)
L1-norm:
(1)
x1=∑i=1nxi.(ii)
Lp-norm to the p-th power:
(1)
xpp=∑i=1n|xi|p;(2)
xpp≈∑i=1nxi2+ϵ2p/2;(3)
xpp≈∑i=1nxi+ϵp-1xi.In the case of no noise,L1-norm is equivalent to L0-norm, and since L1-norm is the only norm with sparsity and convexity, the sparse solution can be approximated by convex optimization method. However, the equivalence between L1-norm and L0-norm is no longer valid in the case of noise, so the effect of using L1-norm to promote sparsity is not obvious in this case. Compared with L1-norm, the nonconvex Lp- norm to the p-th power is closer to the L0-norm, so Lp-norm minimization performs better than L1-norm minimization in recovering sparse signal. Based on the idea of approximation, a Gauss function is employed to approximate the L0-norm in [19]. Moreover, [20] uses the Gauss function as a sparsity regularizer. The approximation of L0-norm by Gauss function can be expressed as follows [6]:(6)x0≈limσ→0Fσx=limσ→0∑i=1nfσxi,fσxi=1-exp⁡-xi22σ2.According to this equation, it is obvious that(7)fσxi≈1ifxi≫σ0ifxi≪σ.When σ is a small enough positive value, the sparse regularizer almost equivalent to L0-norm, so it can promote sparsity. Furthermore, the Gauss function can be optimized by optimization methods because it is differentiable and smoothed. Based on this idea, a hyperbolic tangent (tanh) function is proposed in [21]:(8)x0≈limσ→0∑i=1nexpxi2/2σ2-exp-xi2/2σ2expxi2/2σ2+exp-xi2/2σ2.The tanh function proposed in [21] is a smoothed function, and it performs better in approximating L0-norm than the Gauss function in [20]. Based on this, the tanh function has better performance in sparse signal recovery. Indeed, we confirmed this view by a large number of simulation experiments.In this paper, we propose a CIPF function as new sparsity regularizer and show their effectiveness and advantages over other popular regularizers in promoting sparse solutions with both theoretical analyses and experimental evaluations.This paper is organized as follows. Section2 introduces the main contributions of the proposed RRSL0 algorithm. In Section 3, the steps of the RRSL0 algorithm and the selection of related parameters are described in detail. Section 4 provides the experimental results to evaluate the performance of the proposed method. The paper is finally concluded in Section 5.


## 2. Main Contributions of the Proposed RRSL0 Algorithm

In this paper, based on thePg in (5), we propose a new objective function, which is given by(9)argminx∈Rn12Φx-y22+λWHσx.From the equation, we can see that there is an obvious difference compared with other normalized objective functions. The obvious difference is objective regularizer. In this paper, as an objective regularizer, we not only propose a smoothed function approximating L0 norm but also propose a reweighted function to promote sparsity. This section focuses on the relevant contents of reweighted function W=[w1,w2,…,wn]T and Hσ(x).


### 2.1. New Smoothed Function: CIPF

According to [22], some properties of the smoothed functions are summarized in the following.Property. Let f:R→[-∞,+∞] and define fσ(x)≈fσ(x/σ) for any σ>0. The function f has the Property, if(a)
f is real analytic on (x0,∞) for some x0;(b)
∀x≥0, f′′(x)≥-ϵ0, where ϵ0>0 is some constant;(c)
f is convex on R;(d)
f(x)=0↔x=0;(e)
limx→+∞f(x)=1.It follows immediately fromProperty that {fσ(x)} converges to L0-norm as σ→0+; i.e.,(10)limσ→+∞fσx=0ifx=01otherwise.Based on theProperty, this paper proposes a new smoothed function model called CIPF, which satisfies the Property and better approximates L0-norm. The smoothed function model is given as(11)fσx=cx2cx2+σ2.In (11), regulatory factor c is a large constant (because it satisfies σ2/c→0 for any σ). By experiments, the factor c is determined to be 10, which is the good result of the simulation. σ represents a smoothed factor, and when it is smaller, it will make the proposed model closer to the L0-norm. Obviously,(12)limσ→0⁡fσx=0,x=01,x≠0or approximately (13)fσx≈0,x≪σ1,x≫σis satisfied. Let(14)Hσx=∑i=1nfσxi=∑i=1ncxi2cxi2+σ2where Hσ(x)≈x0 for small values of σ, and the approximation tends to equality when σ→0.Figure1 shows the situation of CIPF model approximates the L0-norm. Obviously, the CIPF model makes a better approximation.Figure 1
Different function used in the literature to approximate theL0-norm and some of them are plotted in this figure and display L0.5-norm for comparison.In conclusion, the merits of CIPF model can be summarized as follows:(i)
It closely approximatesL0-norm;(ii)
It is simpler than Gauss and tanh function model.These merits make it possible to reduce the computational complexity on the premise of ensuring the sparse signal reconstruction accuracy, which is of practical significance for sparse signal reconstruction.


### 2.2. New Reweighted Function

Candès et al. [23] proposed the reweighted L1-norm minimization method, which employs the reweighted norm to enhance the sparsity of the solution. Moreover, they provided an analytical result of the improvement in the sparsity recovery by incorporating reweighted function to the objective function. Pant et al. [24] applied another reweighted smoothed L0-norm minimization method, which uses a similar reweighted function to improve sparsity. The reweighted function can be summarized as follows:(i)
Candès et al.:(15)wi=1xixi≠0∞xi=0;(ii)
Pant et al.:wi=1/|xi|+ζ, ζ is a small enough positive constant.From the two reweighted functions, we can find a phenomenon: a large signal entryxi is reweighted with a small value wi; on the contrary, a small signal entry xi is reweighted with a large value wi. By analysis, the large wi force the solution x to concentrate on the indices where wi is small, and by construction these correspond precisely to the indices where x is nonzero.Combined with the above idea, we propose a new reweighted function, which is given by(16)wi=1exi2,s.t.i=1,2,…,n.As for Candès et al., when signal entry is zero or close to zero, the value ofwi will be very large, which is not suitable for computation by computer. Although Pant et al. noticed the problem and improved the reweighted function to avoid this problem, the constant ζ depends on experience. Actually, the proposed reweighted function can avoid the two problems; therefore, the proposed reweighted function can play a better effect.


### 2.3. Representation of Solution Set in Nullspace Measurement Matrix

As fory=Φx, it is well known that all of its solutions can be parameterized as(17)x=xs+Vξwhere xs denotes particular solution of system equation, which can be described as xs=ΦT(ΦΦT)-1y. V∈Rn×(n-m) is a matrix, and its column vectors vi, i=1,2,…,n-m, are made by orthonormal basis of nullspace measurement matrix Φ. V can be evaluated by singular-value decomposition or, more efficiently, QR decomposition of Φ and ξ∈Rn-m. Adopting the solution can diminish the constraint of the optimization and the dimensions of space of solution are reduced from n to n-m, thereby reducing the complexity of computation.


## 2.1. New Smoothed Function: CIPF

According to [22], some properties of the smoothed functions are summarized in the following.Property. Let f:R→[-∞,+∞] and define fσ(x)≈fσ(x/σ) for any σ>0. The function f has the Property, if(a)
f is real analytic on (x0,∞) for some x0;(b)
∀x≥0, f′′(x)≥-ϵ0, where ϵ0>0 is some constant;(c)
f is convex on R;(d)
f(x)=0↔x=0;(e)
limx→+∞f(x)=1.It follows immediately fromProperty that {fσ(x)} converges to L0-norm as σ→0+; i.e.,(10)limσ→+∞fσx=0ifx=01otherwise.Based on theProperty, this paper proposes a new smoothed function model called CIPF, which satisfies the Property and better approximates L0-norm. The smoothed function model is given as(11)fσx=cx2cx2+σ2.In (11), regulatory factor c is a large constant (because it satisfies σ2/c→0 for any σ). By experiments, the factor c is determined to be 10, which is the good result of the simulation. σ represents a smoothed factor, and when it is smaller, it will make the proposed model closer to the L0-norm. Obviously,(12)limσ→0⁡fσx=0,x=01,x≠0or approximately (13)fσx≈0,x≪σ1,x≫σis satisfied. Let(14)Hσx=∑i=1nfσxi=∑i=1ncxi2cxi2+σ2where Hσ(x)≈x0 for small values of σ, and the approximation tends to equality when σ→0.Figure1 shows the situation of CIPF model approximates the L0-norm. Obviously, the CIPF model makes a better approximation.Figure 1
Different function used in the literature to approximate theL0-norm and some of them are plotted in this figure and display L0.5-norm for comparison.In conclusion, the merits of CIPF model can be summarized as follows:(i)
It closely approximatesL0-norm;(ii)
It is simpler than Gauss and tanh function model.These merits make it possible to reduce the computational complexity on the premise of ensuring the sparse signal reconstruction accuracy, which is of practical significance for sparse signal reconstruction.


## 2.2. New Reweighted Function

Candès et al. [23] proposed the reweighted L1-norm minimization method, which employs the reweighted norm to enhance the sparsity of the solution. Moreover, they provided an analytical result of the improvement in the sparsity recovery by incorporating reweighted function to the objective function. Pant et al. [24] applied another reweighted smoothed L0-norm minimization method, which uses a similar reweighted function to improve sparsity. The reweighted function can be summarized as follows:(i)
Candès et al.:(15)wi=1xixi≠0∞xi=0;(ii)
Pant et al.:wi=1/|xi|+ζ, ζ is a small enough positive constant.From the two reweighted functions, we can find a phenomenon: a large signal entryxi is reweighted with a small value wi; on the contrary, a small signal entry xi is reweighted with a large value wi. By analysis, the large wi force the solution x to concentrate on the indices where wi is small, and by construction these correspond precisely to the indices where x is nonzero.Combined with the above idea, we propose a new reweighted function, which is given by(16)wi=1exi2,s.t.i=1,2,…,n.As for Candès et al., when signal entry is zero or close to zero, the value ofwi will be very large, which is not suitable for computation by computer. Although Pant et al. noticed the problem and improved the reweighted function to avoid this problem, the constant ζ depends on experience. Actually, the proposed reweighted function can avoid the two problems; therefore, the proposed reweighted function can play a better effect.


## 2.3. Representation of Solution Set in Nullspace Measurement Matrix

As fory=Φx, it is well known that all of its solutions can be parameterized as(17)x=xs+Vξwhere xs denotes particular solution of system equation, which can be described as xs=ΦT(ΦΦT)-1y. V∈Rn×(n-m) is a matrix, and its column vectors vi, i=1,2,…,n-m, are made by orthonormal basis of nullspace measurement matrix Φ. V can be evaluated by singular-value decomposition or, more efficiently, QR decomposition of Φ and ξ∈Rn-m. Adopting the solution can diminish the constraint of the optimization and the dimensions of space of solution are reduced from n to n-m, thereby reducing the complexity of computation.


## 3. New Algorithm for CS: RRSL0

### 3.1. RRSL0 Algorithm and Its Steps

As explained above, the proposed objective function can be eventually described as(18)argminξ∈Rn-m12y-Φxs+Vξ22+λWTHσxs+Vξ,where Hσ(xs+Vξ)=c·(xs+Vξ)2/c·(xs+Vξ)2+σ2 is a differentiable smoothed accumulated function. Reweighted function W=1/e(xs+Vξ)2. Let xs=[xs(1),xs(2),…,xs(n)]T, V=[v1,v2,…,vn-m]T; therefore, we have (19)G=∂Hσxs+Vξ∂ξ=VT2cσ2xs+Vξcxs+Vξ2+σ22(20)U=∂G∂ξ=VT2cσ4-6c2σ2xs+Vξcxs+Vξ2+σ23V.According to (19), the gradient of the proposed objective function is written as(21)g=VTΦTy-Φxs+Vξ+λ1exs+Vξ2T2cσ2xs+Vξcxs+Vξ2+σ22.Combining (20) and (21), the second derivative is expressed as(22)H=∂g∂ξ=VTΦTΦ+λWTUV.Solving problem of ULSE is to solve the optimization problem in (18). As for this problem, there are many methods, such split Bregman methods [25–27], FISTA [28], alternating direction methods [29], and gradient descent [30]. In order to reduce the computational complexity, this paper adopts conjugate gradient (CG) method to optimize the proposed objective function. There are two classical methods in CG: FR conjugate gradient (FR-CG) method [31] and PRP conjugate gradient (PRP-CG) method [32]. This paper combines FR-CG with PRP-CG, forming a mixed CG algorithm; here we called it MCG method. Based on this, the MCG method is employed to solve the optimization problem in this paper.The problem firstly can be solved by using a sequentialσ-continuation strategy as detailed in the next paragraph.Givenσ a small target value σmin, and a sufficiently large initial value σmax, after referring to the annealing mechanism in simulated annealing [33], this paper proposes a monotonically decreasing sequence σt∣t=2,3,…,T, which is generated as(23)σt=σmaxa-γt-1,s.t.t=1,2,3,…,T.where γ=loga(σmax/σmin)/T-1, a can be a constant that is larger than 1, and T is the maximum of iterations. Using such a monotonically decreasing sequence can avoid the case of too small σ leading to local optimum.According to CG algorithm, the solutionx is updated as(24)xt+1=xt+αtdt,where the parameter d(t) can be given by(25)dt=-gt,t=1-gt+βt-1dt-1,t≥2the parameter β(t-1) is given as(26)βt-1=βPRPt-1,0<βPRPt-1<βFRt-1,βFRt-1,else.where βFR(t-1)=gt22/gt-122 and βPRP(t-1)=gtTgt-gt-1/gt-122, and the parameter α(t) is updated as(27)αt=gt22dtTHtdtFrom above equations, we can see thatα(t) is positive if H(t) is positive definite (PD). As shown in (22), ΦTΦ is PD, and W is PD, so H(t) is PD if U(t) is PD. To get the PD of U(t), this paper makes the following corrections to the diagonal elements of matrix U(t):(28)ui=uiui>η,ηui⩽η.where η is a small positive constant with a value of 10-5. By this processing in (28), the matrix U(t) is PD, H(t) is PD as well, and hence the optimization direction of each component of x keeps consistency.The above is the optimization process analysis of MCG method. For MCG method, the direction of each iteration of this method is the combination of the last iteration direction and the negative gradient direction, which can overcome the sawtooth phenomenon caused by the steepest descent method.According to the explanation above, we can conclude the steps of proposed RRSL0 algorithm, which is given in Algorithm1. As for σ, it can be shown that function Hσ(x) remains convex in the region where the largest magnitude of the component of x=xs+Vξ is less than σ. As the algorithm starts at the original value ξ(0)=0, the above choice of σ1 ensures that the optimization starts in a convex region. This greatly facilitates the convergence of the RRSL0 algorithm.Algorithm 1:Regularized reweighted smoothed L0-norm minimization (RRSL0) algorithm using MCG method.
Initialization:Φ,x,β,σmax,σmin,T and λ Step1: Set ξ(0)=0,W(0)=In,σ1=σmax,t=0; Step2: QR decomposition ΦT=QR, the V equals to the n-m columns of Q, so the initial value x ( 0 ) = x s + V ξ ( 0 ); Step3: For t=1,2,…,T (1) Set σ=σt, compute σt by equation (23); (2) Set s=0, and iterative termination threshold e; (3) While Res>e (a) Computeξ(s+1) by equations (24), (25), (26), (27), (28); (b) ComputeRes=α(s)d(s)22; (4) Let ξ(t)=ξ(s+1), and compute x(t)=xs+Vξ(t); Step4: Output x=x(t).


### 3.2. Selection of Parameters

The selection of parametersλ and σ will affect the performance of the RRSL0 algorithm; thus this paper discusses the selection of the two above parameters in this section.


#### 3.2.1. Selection of Parameterλ

The choice of parameterλ is closely related to the noise level. If the noise level is known, the discrepancyprinciple (DP) can be used to determine the parameter λ. However, in practice, the noise level cannot be measured accurately. In this case, it is necessary to use data-drivenprinciple to select parameter λ, such as BICprinciple. Ref [34], etc. prove the rationality of BICprinciple under appropriate conditions, so this paper uses this principle to determine λ. The BICprinciple is given by (29)BICλk≔12Φx-y22+ln⁡nmWFσx.Then we determine λ by argminλkBICλk.


#### 3.2.2. Selection of Parameterσ

According to (23), descending sequence of σ is generated by σt=σmaxσmin/σmaxt-1/T-1 (it is obtained through simplification of (23)). Parameter σmin and parameter σmax should be appropriately selected. The selection of σmin and σmax is discussed below.For initialσ value σmax, here let x~=maxx; in order to make the algorithm converge quickly, let parameter σmax satisfies(30)Hσx~=cx~2cx~2+σmax2≤b⇒σmax≥1-bbcx~.From the equation, we can see constantb satisfies 1-b/b≥0; thus 0≤b≤1; here we define constant b as 0.5. Hence σmax=cmax⁡x.As for final valueσmin, when σmin→0, Hσmin(x)→x0. That is, the smaller σmin, the more Hσmin(x) can reflect the sparsity of signal x, but at the same time it is also more sensitive to noise; therefore, the value σmin should not be too small. Combining [19], we choose σmin=0.01.


## 3.1. RRSL0 Algorithm and Its Steps

As explained above, the proposed objective function can be eventually described as(18)argminξ∈Rn-m12y-Φxs+Vξ22+λWTHσxs+Vξ,where Hσ(xs+Vξ)=c·(xs+Vξ)2/c·(xs+Vξ)2+σ2 is a differentiable smoothed accumulated function. Reweighted function W=1/e(xs+Vξ)2. Let xs=[xs(1),xs(2),…,xs(n)]T, V=[v1,v2,…,vn-m]T; therefore, we have (19)G=∂Hσxs+Vξ∂ξ=VT2cσ2xs+Vξcxs+Vξ2+σ22(20)U=∂G∂ξ=VT2cσ4-6c2σ2xs+Vξcxs+Vξ2+σ23V.According to (19), the gradient of the proposed objective function is written as(21)g=VTΦTy-Φxs+Vξ+λ1exs+Vξ2T2cσ2xs+Vξcxs+Vξ2+σ22.Combining (20) and (21), the second derivative is expressed as(22)H=∂g∂ξ=VTΦTΦ+λWTUV.Solving problem of ULSE is to solve the optimization problem in (18). As for this problem, there are many methods, such split Bregman methods [25–27], FISTA [28], alternating direction methods [29], and gradient descent [30]. In order to reduce the computational complexity, this paper adopts conjugate gradient (CG) method to optimize the proposed objective function. There are two classical methods in CG: FR conjugate gradient (FR-CG) method [31] and PRP conjugate gradient (PRP-CG) method [32]. This paper combines FR-CG with PRP-CG, forming a mixed CG algorithm; here we called it MCG method. Based on this, the MCG method is employed to solve the optimization problem in this paper.The problem firstly can be solved by using a sequentialσ-continuation strategy as detailed in the next paragraph.Givenσ a small target value σmin, and a sufficiently large initial value σmax, after referring to the annealing mechanism in simulated annealing [33], this paper proposes a monotonically decreasing sequence σt∣t=2,3,…,T, which is generated as(23)σt=σmaxa-γt-1,s.t.t=1,2,3,…,T.where γ=loga(σmax/σmin)/T-1, a can be a constant that is larger than 1, and T is the maximum of iterations. Using such a monotonically decreasing sequence can avoid the case of too small σ leading to local optimum.According to CG algorithm, the solutionx is updated as(24)xt+1=xt+αtdt,where the parameter d(t) can be given by(25)dt=-gt,t=1-gt+βt-1dt-1,t≥2the parameter β(t-1) is given as(26)βt-1=βPRPt-1,0<βPRPt-1<βFRt-1,βFRt-1,else.where βFR(t-1)=gt22/gt-122 and βPRP(t-1)=gtTgt-gt-1/gt-122, and the parameter α(t) is updated as(27)αt=gt22dtTHtdtFrom above equations, we can see thatα(t) is positive if H(t) is positive definite (PD). As shown in (22), ΦTΦ is PD, and W is PD, so H(t) is PD if U(t) is PD. To get the PD of U(t), this paper makes the following corrections to the diagonal elements of matrix U(t):(28)ui=uiui>η,ηui⩽η.where η is a small positive constant with a value of 10-5. By this processing in (28), the matrix U(t) is PD, H(t) is PD as well, and hence the optimization direction of each component of x keeps consistency.The above is the optimization process analysis of MCG method. For MCG method, the direction of each iteration of this method is the combination of the last iteration direction and the negative gradient direction, which can overcome the sawtooth phenomenon caused by the steepest descent method.According to the explanation above, we can conclude the steps of proposed RRSL0 algorithm, which is given in Algorithm1. As for σ, it can be shown that function Hσ(x) remains convex in the region where the largest magnitude of the component of x=xs+Vξ is less than σ. As the algorithm starts at the original value ξ(0)=0, the above choice of σ1 ensures that the optimization starts in a convex region. This greatly facilitates the convergence of the RRSL0 algorithm.Algorithm 1:Regularized reweighted smoothed L0-norm minimization (RRSL0) algorithm using MCG method.
Initialization:Φ,x,β,σmax,σmin,T and λ Step1: Set ξ(0)=0,W(0)=In,σ1=σmax,t=0; Step2: QR decomposition ΦT=QR, the V equals to the n-m columns of Q, so the initial value x ( 0 ) = x s + V ξ ( 0 ); Step3: For t=1,2,…,T (1) Set σ=σt, compute σt by equation (23); (2) Set s=0, and iterative termination threshold e; (3) While Res>e (a) Computeξ(s+1) by equations (24), (25), (26), (27), (28); (b) ComputeRes=α(s)d(s)22; (4) Let ξ(t)=ξ(s+1), and compute x(t)=xs+Vξ(t); Step4: Output x=x(t).


## 3.2. Selection of Parameters

The selection of parametersλ and σ will affect the performance of the RRSL0 algorithm; thus this paper discusses the selection of the two above parameters in this section.


### 3.2.1. Selection of Parameterλ

The choice of parameterλ is closely related to the noise level. If the noise level is known, the discrepancyprinciple (DP) can be used to determine the parameter λ. However, in practice, the noise level cannot be measured accurately. In this case, it is necessary to use data-drivenprinciple to select parameter λ, such as BICprinciple. Ref [34], etc. prove the rationality of BICprinciple under appropriate conditions, so this paper uses this principle to determine λ. The BICprinciple is given by (29)BICλk≔12Φx-y22+ln⁡nmWFσx.Then we determine λ by argminλkBICλk.


### 3.2.2. Selection of Parameterσ

According to (23), descending sequence of σ is generated by σt=σmaxσmin/σmaxt-1/T-1 (it is obtained through simplification of (23)). Parameter σmin and parameter σmax should be appropriately selected. The selection of σmin and σmax is discussed below.For initialσ value σmax, here let x~=maxx; in order to make the algorithm converge quickly, let parameter σmax satisfies(30)Hσx~=cx~2cx~2+σmax2≤b⇒σmax≥1-bbcx~.From the equation, we can see constantb satisfies 1-b/b≥0; thus 0≤b≤1; here we define constant b as 0.5. Hence σmax=cmax⁡x.As for final valueσmin, when σmin→0, Hσmin(x)→x0. That is, the smaller σmin, the more Hσmin(x) can reflect the sparsity of signal x, but at the same time it is also more sensitive to noise; therefore, the value σmin should not be too small. Combining [19], we choose σmin=0.01.


## 3.2.1. Selection of Parameterλ

The choice of parameterλ is closely related to the noise level. If the noise level is known, the discrepancyprinciple (DP) can be used to determine the parameter λ. However, in practice, the noise level cannot be measured accurately. In this case, it is necessary to use data-drivenprinciple to select parameter λ, such as BICprinciple. Ref [34], etc. prove the rationality of BICprinciple under appropriate conditions, so this paper uses this principle to determine λ. The BICprinciple is given by (29)BICλk≔12Φx-y22+ln⁡nmWFσx.Then we determine λ by argminλkBICλk.


## 3.2.2. Selection of Parameterσ

According to (23), descending sequence of σ is generated by σt=σmaxσmin/σmaxt-1/T-1 (it is obtained through simplification of (23)). Parameter σmin and parameter σmax should be appropriately selected. The selection of σmin and σmax is discussed below.For initialσ value σmax, here let x~=maxx; in order to make the algorithm converge quickly, let parameter σmax satisfies(30)Hσx~=cx~2cx~2+σmax2≤b⇒σmax≥1-bbcx~.From the equation, we can see constantb satisfies 1-b/b≥0; thus 0≤b≤1; here we define constant b as 0.5. Hence σmax=cmax⁡x.As for final valueσmin, when σmin→0, Hσmin(x)→x0. That is, the smaller σmin, the more Hσmin(x) can reflect the sparsity of signal x, but at the same time it is also more sensitive to noise; therefore, the value σmin should not be too small. Combining [19], we choose σmin=0.01.


## 4. Numerical Simulation and Analysis

The numerical simulation platform is MATLAB 2017b, which is installed on the WINDOWS 10, 64-bit operating system. The CPU of simulation computer is Intel (R) Core (TM) i5-3230M, and the frequency is 2.6 GHz. In this section, the performance of RRSL0 algorithm is verified by signal and image recovery in the noiseless case and noise case.Here some state-of-the-art algorithms are selected for comparison. The parameters are selected to obtain best performance for each algorithm: for BPDN algorithm [35], λ=σN2log⁡(n); for SL0 algorithm, δmax=2max⁡x, δmin=0.01, and scale factors are set as L=5, ρ=0.8; for L2-SL0 algorithm, δmax=4max⁡x, δmin=0.01, L=10, and ρ=0.8; for Lp-RLS algorithm, T=80, p1=1, pT=0.1, ϵ1=1, ϵT=0.01, and Et=10-25; and for RRSL0 algorithm, smoothed factor σmax=cmax⁡x, σmin=0.01, and iterations T=80. All experiments are based on the 100 trials.


### 4.1. Comparison of Reconstruction Performance

For signal recovery under no noise conditions, we evaluate performance of algorithms byreconstructionsuccessrate (RSR), normlizedmeansquareerror (NMSE), and cpurunningtime (CRT). RSR is defined as the ratio of the number of successful experiments to the total number of experiments, and if x-x^2/x2<10-3, the recovery is considered to be a success. NMSE is defined as x-x^2/x2. CRT is measured with tic and toc. The case of signal recovery under noise is considered by NMSE and signaltonoiseratio (SNR). SNR is defined as 20log⁡(x-x^2/x2).


#### 4.1.1. Signal Recovery in Noiseless Case

In this section, the SL0, BPDN, L2-SL0, and Lp-RLS are used for comparison. Here we fixn=256 and m=100 and the sparsity K=4q+1, q=1,2,…,15 for Figure 2 and let n=[170,220,270,320,370,420,470,520], m=n/2, K=n/5 for Table 1. For every experiment, we generate a pair of x,Φ: Φ is a m×n random Gauss matrix with normalized rows; the nonzero entries of the sparse signal x∈Rn are i.i.d. generated according to the Gaussian distribution N(0,1).Table 1
Signal CRT analysis for BPDN, SL0, L2-SL0,Lp-RLS, and proposed RRSL0 with signal length changes according to sequence [170,220,270,320,370,420,470,520] while 100 runs.
Signal length (n) CPU running time (seconds) BPDN SL0 L2-SL0 Lp-RLS RRSL0 170 0.1958 0.0572 0.0909 0.0636 0.0628 220 0.2899 0.1385 0.2298 0.1503 0.1425 270 0.4953 0.2285 0.4257 0.3050 0.2908 320 0.7675 0.3196 0.6385 0.5121 0.5088 370 1.0586 0.4559 0.9260 0.9023 0.8924 420 1.4767 0.6133 1.1327 1.0906 1.0170 470 1.9414 0.7956 1.4782 1.4178 1.3438 520 2.6194 1.0383 2.0894 1.9097 1.8818Figure 2
Signal recovery in noiseless case. RSR and NMSE of recovery signal are displayed at intervals of 5 with sparsity from 1 to 61. The comparison algorithms are BPDN, SL0, L2-SL0, Lp-RLS, and RRSL0. The experimental results were obtained by 100 independent repeated experiments.
(a)
RSR of recovery signal (b)
NMSE of recovery signalFigure2(a) shows the RSR of all algorithms. When K is large enough and continues to increase, the RSR of all algorithms tends to decrease until they drop to zero. In fact, the speed of changes in the RSR reflects the robustness of the algorithm. Specifically, the slower the RSR changes, the better the robustness of the algorithm. Obviously, the RRSL0 algorithm has a more slowly changing trend than other comparison algorithms. By comparison, SL0 is superior to BPDN and L2-SL0, but inferior to Lp-RLS and RRSL0. BPDN, L2-SL0, Lp-RLS, and RRSL0 are the unconstrained optimization algorithms. When recovering signal in noiseless case, the deviation term Φx-y22 is equivalent to the error between original signal and recovered signal, which interferes with the robustness of the algorithm. Luckily, the Lp-RLS adopts FR-CG algorithm to optimize, which improves the stability of the optimized solution. The proposed RRSL0 algorithm applies MCG algorithm to enhance the reconstruction performance.The NMSE of all algorithms is shown in Figure2(b). NMSE reflects the signal reconstruction accuracy, and the lower the NMSE, the higher the accuracy of the reconstructed signal. As shown in this figure, all algorithms perform well when the sparsity is lower than 29. NMSE of BDPN increases fastest, followed by L2-SL0, SL0, Lp-RLS, and proposed RRSL0. When sparsity is increasing, L2-SL0 does not perform well, while Lp-RLS and RRSL0 perform well. This result indirectly verified that MCG and FR-CG method are superior to steepest descent method in SL0 and L2-SL0; furthermore, the term Φx-y22 has a relatively large influence on the recovery results. In addition, the experiment also verified that RRSL0 performs better than Lp-RLS and other algorithms.Table1 shows the CRT of all algorithms. The n changes according to a given sequence [170,220,270,320,370,420,470,520]. From the table, for any n, SL0 has the shortest computation time, followed by RRSL0, Lp-RLS, and L2-SL0, and BPDN takes the longest computation time. BPDN algorithm is generally implemented by quadratic programming method, and the computational complexity of this method is very high, thus resulting in a large increase in the overall computation time of the algorithm. Furthermore, In L2-SL0, the step factor needs to be obtained by a one-dimensional linear search method, while Lp-RLS and RRSL0 do not. Compared with Lp-RLS, RRSL0 is more prominent in the decrease of computation time. In order to verify the difference in computation time between the two algorithms, the convergence comparison of the two algorithms is shown in Figure 3. From this figure, we can draw that the convergence rate of RRSL0 is much faster than Lp-RLS.Figure 3
Convergence verification forLp-RLS and proposed RRSL0 with signal length n=256, m=100, and K=20 with 100 runs.


#### 4.1.2. Signal Recovery in Noise Case

In this section, we discuss signal recovery performance in noise case. We add noisew to the measurement vector y; moreover, w is randomly formed and follows the Gauss distribution of δN×N(0,1). In order to analyze the antinoise performance of the RRSL0 algorithm more closely to the real situation, we constructed a certain signal as an experimental object in the experiments in this section. The certain signal is given by(31)X=α1sin⁡2πf1Tst+β1cos⁡2πf2Tst+α2sin⁡2πf3Tst+β2cos⁡2πf4Tstwhere α1=0.2, α2=0.1, β1=0.3, and β2=0.4. f1=50Hz, f2=100Hz, f3=200Hz, and f4=300Hz. Here t is a sequence with t=[1,2,3,…,n], Ts is sampling interval valued 1/fs. fs is sampling frequency valued 800Hz. The object that needs to be reconstructed can be expressed as(32)Y=ψφ+δN×N0,1m×1=Φ+δN×N0,1m×1.where X=φ, φT is the Fouriertransform basis (here called sparse basis matrix) given by FourierIn×n, In×n is a unit matrix, and ∈Rn is a sparse signal in frequency domain obtained by the Fouriertransform of X. Moreover, ψ∈Rm×n is a random matrix generated by Gauss distribution, sensing matrix Φ=ψφ. Y∈Rm. Here, let n=128 and m=64. Therefore, the sparse signal  can be recovered from given {Y,Φ} by CS recovery methods, and then the original signal X can be obtained by X=φ.Figure4 shows the signal recovery effect by different algorithms under noise intensity δN = 0.2. Meanwhile, time-frequency characteristics are shown in Figure 5. Obviously, BPDN and SL0 perform not well while L2-SL0, Lp-RLS, and proposed RRSL0 perform quite well. This verifies that the regularization mechanism has good antinoise effect. The performance of each algorithm under different noise intensities is shown in Table 2. When δN=0, the SL0 outperforms other algorithms, but with the increase of δN, the effect of SL0 is getting worse and worse. This result further illustrates that the traditional constrained sparse recovery algorithm does not have the performance of antinoise. For BPDN, L2-SL0, Lp-RLS, and RRSL0, they all applied regularization mechanism, and they are indeed superior to SL0 in noise case. Therefore, the proposed RRSL0 in this paper has the best antinoise performance.Table 2
NMSE analysis for BPDN, SL0, L2-SL0,Lp-RLS, and proposed RRSL0 with noise intensity changes according to sequence δN=[0,0.01,0.05,0.1,0.2,0.5] while 100 runs.
Noise intensity (δN) NMSE of signal recovery BPDN L0 L2-SL0 Lp-RLS RRSL0 0 2.553e-05 1.311e-07 4.563e-07 1.563e-07 2.338e-07 0.01 4.659e-04 3.054e-04 2.352e-04 7.930e-05 6.032e-05 0.05 1.607e-03 4.079e-03 9.012e-04 2.384e-04 1.150e-04 0.1 7.633e-03 8.424e-03 3.309e-03 7.793e-04 5.522e-04 0.2 1.381e-02 2.090e-02 6.091e-03 1.130e-03 8.892e-04 0.5 3.063e-02 6.908e-02 1.083e-02 5.743e-03 3.041e-03Figure 4
Signal recovery effect by different algorithms when noise intensityδN = 0.2. (a) signal recovery by the BPDN algorithm; (b) signal recovery by the SL0 algorithm; (c) signal recovery by the L2-SL0 algorithm; (d) signal recovery by the Lp-RLS algorithm; (e) signal recovery by the RRSL0 algorithm.
(a) (b) (c) (d) (e)Figure 5
Reconstruction time-frequency characteristics of different algorithms when noise intensityδN = 0.2. (a) time-frequency characteristics of the original signal; (b) BPDN; (c) SL0; (d) L2-SL0; (e) Lp-RLS; (f) RRSL0.
(a) (b) (c) (d) (e) (f)


### 4.2. Comparison of Image Recovery Performance

Real images are considered to be approximately sparse under some proper basis, such as the Discrete Cosine Transform (DCT) basis and Discrete Wavelet Transform (DWT) basis. Here we choose DWT basis to recover these images. We compare the recovery performances based on the 2 real images in Figure6: Boat and Barbara, and the size of these images is 256×256. The compression ratio (CR, it is defined as m/n) is 0.5, and noise δN equals 0.01. We still choose SL0, BPDN, L2-SL0, and Lp-RLS to make comparisons. For image recovery, the object of image processing is given by(33)Y=ΦX+bHere Y, X, and b are matrices; among these, Y,b∈Rm×n, X∈Rn×n. In order to meet the basic requirements of CS, we make the following processing:(34)Y=Y1,Y2,…,Ynwhere Yi, i=1,2,…,n, is the column vector of Y. Therefore we can process images by CS.Figure 6
Original images.
(a)
Original Boat (b)
Original BarbaraFor performance of image recovery, we valuate it byPeakSignaltoNoiseRatio (PSNR) and StructuralSimilarityIndex (SSIM). PSNR is defined as(35)PSNR=10log⁡2552MSEwhere MSE=x-x^22, and SSIM is defined as(36)SSIMp,q=2μp+μq+c12σpq+c2μp2+μq2+c1σp2+σq2+c2.Among this,μp is the mean of image p, μq is the mean of image q, σp is the variance of image p, σq is the variance of image q, and σpq is covariance between image p and image q. Parameters c1=z1L and c2=z2L, where z1=0.01,z2=0.03, and L is dynamic range of pixel values. Range of SSIM is [-1,1]; when these two images are the same, the SSIM equals to 1.Figure6 shows the original images. Figure 7 shows the recovery effect of Boat and Barbara with noise intensity δN=0.01 and δN=0.1. In the Boat, the recovered images by SL0 and BPDN have the obvious water ripples while recovered images by other algorithms have no such water ripples. Similarly, for Barbara, recovered images by SL0 and BPDN are blurred compared with recovered images by other algorithms. L2-SL0, Lp-RLS, and RRSL0 algorithms are also effective for noisy image recovery. For L2-SL0, Lp-RLS, and RRSL0 algorithms, their recovery effects are very similar. In order to further analyze the advantages and disadvantages of these algorithms, we analyze PSNR and SSIM of the images recovered by these algorithms and results are shown in Tables 3–6. By observation and analysis, Lp-RLS performs better than L2-SL0 on the whole (except SSIM when δN=0.05, because the robustness of Lp-RLS is not strong); at the same time, RRSL0 outperforms Lp-RLS. Hence, the RRSL0 proposed by this paper is superior to other selected algorithms in image processing.Table 3
PSNR analysis of recovered Boat image by BPDN, SL0, L2-SL0,Lp-RLS, and RRSL0 with noise intensity changes according to sequence δN=[0,0.01,0.05,0.1,0.2,0.5].
Noise intensity (δN) PSNR of Boat image recovery BPDN SL0 L2-SL0 L p-RLS RRSL0 0 28.4928 31.5901 31.7175 32.2332 33.3495 0.01 27.9775 29.8589 30.2423 30.9742 31.5447 0.05 24.0531 27.1226 27.7101 28.2422 29.2505 0.1 20.3769 21.6485 23.9331 25.8639 26.4745 0.2 16.3044 17.9303 20.6744 23.0823 23.7415 0.5 9.2283 10.7863 13.4863 14.4974 16.3646Table 4
SSIM analysis of recovered Boat image by BPDN, SL0, L2-SL0,Lp-RLS, and RRSL0 with noise intensity changes according to sequence δN=[0,0.01,0.05,0.1,0.2,0.5].
Noise intensity (δN) SSIM of Boat image recovery BPDN SL0 L2-SL0 L p-RLS RRSL0 0 0.9865 0.9901 0.9902 0.9914 0.9934 0.01 0.9764 0.9852 0.9865 0.9887 0.9898 0.05 0.9418 0.9714 0.9756 0.9789 0.9827 0.1 0.8827 0.9083 0.9403 0.9631 0.9669 0.2 0.7452 0.8077 0.8868 0.9284 0.9396 0.5 0.3571 0.4427 0.5966 0.6492 0.7351Table 5
PSNR analysis of recovered Barbara image by BPDN, SL0, L2-SL0,Lp-RLS, and RRSL0 with noise intensity changes according to sequence δN=[0,0.01,0.05,0.1,0.2,0.5].
Noise intensity (δN) PSNR of Barbara image recovery BPDN SL0 L2-SL0 L p-RLS RRSL0 0 28.6613 29.5761 31.7921 32.3566 33.4074 0.01 27.7887 30.4746 30.9875 31.7365 32.6767 0.05 26.0316 26.1561 26.3147 27.3533 29.8613 0.1 20.7076 21.6914 23.8343 26.9408 27.5167 0.2 16.6564 16.6721 19.1297 19.2624 19.7801 0.5 9.1316 9.3525 13.1838 15.7972 16.9070Table 6
SSIM analysis of recovered Barbara image by BPDN, SL0, L2-SL0,Lp-RLS, and RRSL0 with noise intensity changes according to sequence δN=[0,0.01,0.05,0.1,0.2,0.5].
Noise intensity (δN) SSIM of Barbara image recovery BPDN SL0 L2-SL0 L p-RLS RRSL0 0 0.9834 0.9920 0.9926 0.9936 0.9984 0.01 0.9817 0.9904 0.9914 0.9927 0.9942 0.05 0.9797 0.9744 0.9753 0.9739 0.9888 0.1 0.9161 0.9314 0.9544 0.9728 0.9810 0.2 0.8165 0.8170 0.9051 0.9142 0.9278 0.5 0.3863 0.3992 0.6309 0.7712 0.8190Figure 7
Images recovery effect by BPDN, SL0, L2-SL0,Lp-RLS, and RRSL0 algorithms under different noise intensities. In subfigures (a), (b), (c), and (d), from left to right are images recovered by BPDN, SL0, L2-SL0, Lp-RLS, and RRSL0 algorithms.
(a)
Recovered Boat with noise intensityδN=0.01 (b)
Recovered Barbara with noise intensityδN=0.01 (c)
Recovered Boat with noise intensityδN=0.1 (d)
Recovered Barbara with noise intensityδN=0.1


## 4.1. Comparison of Reconstruction Performance

For signal recovery under no noise conditions, we evaluate performance of algorithms byreconstructionsuccessrate (RSR), normlizedmeansquareerror (NMSE), and cpurunningtime (CRT). RSR is defined as the ratio of the number of successful experiments to the total number of experiments, and if x-x^2/x2<10-3, the recovery is considered to be a success. NMSE is defined as x-x^2/x2. CRT is measured with tic and toc. The case of signal recovery under noise is considered by NMSE and signaltonoiseratio (SNR). SNR is defined as 20log⁡(x-x^2/x2).


### 4.1.1. Signal Recovery in Noiseless Case

In this section, the SL0, BPDN, L2-SL0, and Lp-RLS are used for comparison. Here we fixn=256 and m=100 and the sparsity K=4q+1, q=1,2,…,15 for Figure 2 and let n=[170,220,270,320,370,420,470,520], m=n/2, K=n/5 for Table 1. For every experiment, we generate a pair of x,Φ: Φ is a m×n random Gauss matrix with normalized rows; the nonzero entries of the sparse signal x∈Rn are i.i.d. generated according to the Gaussian distribution N(0,1).Table 1
Signal CRT analysis for BPDN, SL0, L2-SL0,Lp-RLS, and proposed RRSL0 with signal length changes according to sequence [170,220,270,320,370,420,470,520] while 100 runs.
Signal length (n) CPU running time (seconds) BPDN SL0 L2-SL0 Lp-RLS RRSL0 170 0.1958 0.0572 0.0909 0.0636 0.0628 220 0.2899 0.1385 0.2298 0.1503 0.1425 270 0.4953 0.2285 0.4257 0.3050 0.2908 320 0.7675 0.3196 0.6385 0.5121 0.5088 370 1.0586 0.4559 0.9260 0.9023 0.8924 420 1.4767 0.6133 1.1327 1.0906 1.0170 470 1.9414 0.7956 1.4782 1.4178 1.3438 520 2.6194 1.0383 2.0894 1.9097 1.8818Figure 2
Signal recovery in noiseless case. RSR and NMSE of recovery signal are displayed at intervals of 5 with sparsity from 1 to 61. The comparison algorithms are BPDN, SL0, L2-SL0, Lp-RLS, and RRSL0. The experimental results were obtained by 100 independent repeated experiments.
(a)
RSR of recovery signal (b)
NMSE of recovery signalFigure2(a) shows the RSR of all algorithms. When K is large enough and continues to increase, the RSR of all algorithms tends to decrease until they drop to zero. In fact, the speed of changes in the RSR reflects the robustness of the algorithm. Specifically, the slower the RSR changes, the better the robustness of the algorithm. Obviously, the RRSL0 algorithm has a more slowly changing trend than other comparison algorithms. By comparison, SL0 is superior to BPDN and L2-SL0, but inferior to Lp-RLS and RRSL0. BPDN, L2-SL0, Lp-RLS, and RRSL0 are the unconstrained optimization algorithms. When recovering signal in noiseless case, the deviation term Φx-y22 is equivalent to the error between original signal and recovered signal, which interferes with the robustness of the algorithm. Luckily, the Lp-RLS adopts FR-CG algorithm to optimize, which improves the stability of the optimized solution. The proposed RRSL0 algorithm applies MCG algorithm to enhance the reconstruction performance.The NMSE of all algorithms is shown in Figure2(b). NMSE reflects the signal reconstruction accuracy, and the lower the NMSE, the higher the accuracy of the reconstructed signal. As shown in this figure, all algorithms perform well when the sparsity is lower than 29. NMSE of BDPN increases fastest, followed by L2-SL0, SL0, Lp-RLS, and proposed RRSL0. When sparsity is increasing, L2-SL0 does not perform well, while Lp-RLS and RRSL0 perform well. This result indirectly verified that MCG and FR-CG method are superior to steepest descent method in SL0 and L2-SL0; furthermore, the term Φx-y22 has a relatively large influence on the recovery results. In addition, the experiment also verified that RRSL0 performs better than Lp-RLS and other algorithms.Table1 shows the CRT of all algorithms. The n changes according to a given sequence [170,220,270,320,370,420,470,520]. From the table, for any n, SL0 has the shortest computation time, followed by RRSL0, Lp-RLS, and L2-SL0, and BPDN takes the longest computation time. BPDN algorithm is generally implemented by quadratic programming method, and the computational complexity of this method is very high, thus resulting in a large increase in the overall computation time of the algorithm. Furthermore, In L2-SL0, the step factor needs to be obtained by a one-dimensional linear search method, while Lp-RLS and RRSL0 do not. Compared with Lp-RLS, RRSL0 is more prominent in the decrease of computation time. In order to verify the difference in computation time between the two algorithms, the convergence comparison of the two algorithms is shown in Figure 3. From this figure, we can draw that the convergence rate of RRSL0 is much faster than Lp-RLS.Figure 3
Convergence verification forLp-RLS and proposed RRSL0 with signal length n=256, m=100, and K=20 with 100 runs.


### 4.1.2. Signal Recovery in Noise Case

In this section, we discuss signal recovery performance in noise case. We add noisew to the measurement vector y; moreover, w is randomly formed and follows the Gauss distribution of δN×N(0,1). In order to analyze the antinoise performance of the RRSL0 algorithm more closely to the real situation, we constructed a certain signal as an experimental object in the experiments in this section. The certain signal is given by(31)X=α1sin⁡2πf1Tst+β1cos⁡2πf2Tst+α2sin⁡2πf3Tst+β2cos⁡2πf4Tstwhere α1=0.2, α2=0.1, β1=0.3, and β2=0.4. f1=50Hz, f2=100Hz, f3=200Hz, and f4=300Hz. Here t is a sequence with t=[1,2,3,…,n], Ts is sampling interval valued 1/fs. fs is sampling frequency valued 800Hz. The object that needs to be reconstructed can be expressed as(32)Y=ψφ+δN×N0,1m×1=Φ+δN×N0,1m×1.where X=φ, φT is the Fouriertransform basis (here called sparse basis matrix) given by FourierIn×n, In×n is a unit matrix, and ∈Rn is a sparse signal in frequency domain obtained by the Fouriertransform of X. Moreover, ψ∈Rm×n is a random matrix generated by Gauss distribution, sensing matrix Φ=ψφ. Y∈Rm. Here, let n=128 and m=64. Therefore, the sparse signal  can be recovered from given {Y,Φ} by CS recovery methods, and then the original signal X can be obtained by X=φ.Figure4 shows the signal recovery effect by different algorithms under noise intensity δN = 0.2. Meanwhile, time-frequency characteristics are shown in Figure 5. Obviously, BPDN and SL0 perform not well while L2-SL0, Lp-RLS, and proposed RRSL0 perform quite well. This verifies that the regularization mechanism has good antinoise effect. The performance of each algorithm under different noise intensities is shown in Table 2. When δN=0, the SL0 outperforms other algorithms, but with the increase of δN, the effect of SL0 is getting worse and worse. This result further illustrates that the traditional constrained sparse recovery algorithm does not have the performance of antinoise. For BPDN, L2-SL0, Lp-RLS, and RRSL0, they all applied regularization mechanism, and they are indeed superior to SL0 in noise case. Therefore, the proposed RRSL0 in this paper has the best antinoise performance.Table 2
NMSE analysis for BPDN, SL0, L2-SL0,Lp-RLS, and proposed RRSL0 with noise intensity changes according to sequence δN=[0,0.01,0.05,0.1,0.2,0.5] while 100 runs.
Noise intensity (δN) NMSE of signal recovery BPDN L0 L2-SL0 Lp-RLS RRSL0 0 2.553e-05 1.311e-07 4.563e-07 1.563e-07 2.338e-07 0.01 4.659e-04 3.054e-04 2.352e-04 7.930e-05 6.032e-05 0.05 1.607e-03 4.079e-03 9.012e-04 2.384e-04 1.150e-04 0.1 7.633e-03 8.424e-03 3.309e-03 7.793e-04 5.522e-04 0.2 1.381e-02 2.090e-02 6.091e-03 1.130e-03 8.892e-04 0.5 3.063e-02 6.908e-02 1.083e-02 5.743e-03 3.041e-03Figure 4
Signal recovery effect by different algorithms when noise intensityδN = 0.2. (a) signal recovery by the BPDN algorithm; (b) signal recovery by the SL0 algorithm; (c) signal recovery by the L2-SL0 algorithm; (d) signal recovery by the Lp-RLS algorithm; (e) signal recovery by the RRSL0 algorithm.
(a) (b) (c) (d) (e)Figure 5
Reconstruction time-frequency characteristics of different algorithms when noise intensityδN = 0.2. (a) time-frequency characteristics of the original signal; (b) BPDN; (c) SL0; (d) L2-SL0; (e) Lp-RLS; (f) RRSL0.
(a) (b) (c) (d) (e) (f)


## 4.1.1. Signal Recovery in Noiseless Case

In this section, the SL0, BPDN, L2-SL0, and Lp-RLS are used for comparison. Here we fixn=256 and m=100 and the sparsity K=4q+1, q=1,2,…,15 for Figure 2 and let n=[170,220,270,320,370,420,470,520], m=n/2, K=n/5 for Table 1. For every experiment, we generate a pair of x,Φ: Φ is a m×n random Gauss matrix with normalized rows; the nonzero entries of the sparse signal x∈Rn are i.i.d. generated according to the Gaussian distribution N(0,1).Table 1
Signal CRT analysis for BPDN, SL0, L2-SL0,Lp-RLS, and proposed RRSL0 with signal length changes according to sequence [170,220,270,320,370,420,470,520] while 100 runs.
Signal length (n) CPU running time (seconds) BPDN SL0 L2-SL0 Lp-RLS RRSL0 170 0.1958 0.0572 0.0909 0.0636 0.0628 220 0.2899 0.1385 0.2298 0.1503 0.1425 270 0.4953 0.2285 0.4257 0.3050 0.2908 320 0.7675 0.3196 0.6385 0.5121 0.5088 370 1.0586 0.4559 0.9260 0.9023 0.8924 420 1.4767 0.6133 1.1327 1.0906 1.0170 470 1.9414 0.7956 1.4782 1.4178 1.3438 520 2.6194 1.0383 2.0894 1.9097 1.8818Figure 2
Signal recovery in noiseless case. RSR and NMSE of recovery signal are displayed at intervals of 5 with sparsity from 1 to 61. The comparison algorithms are BPDN, SL0, L2-SL0, Lp-RLS, and RRSL0. The experimental results were obtained by 100 independent repeated experiments.
(a)
RSR of recovery signal (b)
NMSE of recovery signalFigure2(a) shows the RSR of all algorithms. When K is large enough and continues to increase, the RSR of all algorithms tends to decrease until they drop to zero. In fact, the speed of changes in the RSR reflects the robustness of the algorithm. Specifically, the slower the RSR changes, the better the robustness of the algorithm. Obviously, the RRSL0 algorithm has a more slowly changing trend than other comparison algorithms. By comparison, SL0 is superior to BPDN and L2-SL0, but inferior to Lp-RLS and RRSL0. BPDN, L2-SL0, Lp-RLS, and RRSL0 are the unconstrained optimization algorithms. When recovering signal in noiseless case, the deviation term Φx-y22 is equivalent to the error between original signal and recovered signal, which interferes with the robustness of the algorithm. Luckily, the Lp-RLS adopts FR-CG algorithm to optimize, which improves the stability of the optimized solution. The proposed RRSL0 algorithm applies MCG algorithm to enhance the reconstruction performance.The NMSE of all algorithms is shown in Figure2(b). NMSE reflects the signal reconstruction accuracy, and the lower the NMSE, the higher the accuracy of the reconstructed signal. As shown in this figure, all algorithms perform well when the sparsity is lower than 29. NMSE of BDPN increases fastest, followed by L2-SL0, SL0, Lp-RLS, and proposed RRSL0. When sparsity is increasing, L2-SL0 does not perform well, while Lp-RLS and RRSL0 perform well. This result indirectly verified that MCG and FR-CG method are superior to steepest descent method in SL0 and L2-SL0; furthermore, the term Φx-y22 has a relatively large influence on the recovery results. In addition, the experiment also verified that RRSL0 performs better than Lp-RLS and other algorithms.Table1 shows the CRT of all algorithms. The n changes according to a given sequence [170,220,270,320,370,420,470,520]. From the table, for any n, SL0 has the shortest computation time, followed by RRSL0, Lp-RLS, and L2-SL0, and BPDN takes the longest computation time. BPDN algorithm is generally implemented by quadratic programming method, and the computational complexity of this method is very high, thus resulting in a large increase in the overall computation time of the algorithm. Furthermore, In L2-SL0, the step factor needs to be obtained by a one-dimensional linear search method, while Lp-RLS and RRSL0 do not. Compared with Lp-RLS, RRSL0 is more prominent in the decrease of computation time. In order to verify the difference in computation time between the two algorithms, the convergence comparison of the two algorithms is shown in Figure 3. From this figure, we can draw that the convergence rate of RRSL0 is much faster than Lp-RLS.Figure 3
Convergence verification forLp-RLS and proposed RRSL0 with signal length n=256, m=100, and K=20 with 100 runs.


## 4.1.2. Signal Recovery in Noise Case

In this section, we discuss signal recovery performance in noise case. We add noisew to the measurement vector y; moreover, w is randomly formed and follows the Gauss distribution of δN×N(0,1). In order to analyze the antinoise performance of the RRSL0 algorithm more closely to the real situation, we constructed a certain signal as an experimental object in the experiments in this section. The certain signal is given by(31)X=α1sin⁡2πf1Tst+β1cos⁡2πf2Tst+α2sin⁡2πf3Tst+β2cos⁡2πf4Tstwhere α1=0.2, α2=0.1, β1=0.3, and β2=0.4. f1=50Hz, f2=100Hz, f3=200Hz, and f4=300Hz. Here t is a sequence with t=[1,2,3,…,n], Ts is sampling interval valued 1/fs. fs is sampling frequency valued 800Hz. The object that needs to be reconstructed can be expressed as(32)Y=ψφ+δN×N0,1m×1=Φ+δN×N0,1m×1.where X=φ, φT is the Fouriertransform basis (here called sparse basis matrix) given by FourierIn×n, In×n is a unit matrix, and ∈Rn is a sparse signal in frequency domain obtained by the Fouriertransform of X. Moreover, ψ∈Rm×n is a random matrix generated by Gauss distribution, sensing matrix Φ=ψφ. Y∈Rm. Here, let n=128 and m=64. Therefore, the sparse signal  can be recovered from given {Y,Φ} by CS recovery methods, and then the original signal X can be obtained by X=φ.Figure4 shows the signal recovery effect by different algorithms under noise intensity δN = 0.2. Meanwhile, time-frequency characteristics are shown in Figure 5. Obviously, BPDN and SL0 perform not well while L2-SL0, Lp-RLS, and proposed RRSL0 perform quite well. This verifies that the regularization mechanism has good antinoise effect. The performance of each algorithm under different noise intensities is shown in Table 2. When δN=0, the SL0 outperforms other algorithms, but with the increase of δN, the effect of SL0 is getting worse and worse. This result further illustrates that the traditional constrained sparse recovery algorithm does not have the performance of antinoise. For BPDN, L2-SL0, Lp-RLS, and RRSL0, they all applied regularization mechanism, and they are indeed superior to SL0 in noise case. Therefore, the proposed RRSL0 in this paper has the best antinoise performance.Table 2
NMSE analysis for BPDN, SL0, L2-SL0,Lp-RLS, and proposed RRSL0 with noise intensity changes according to sequence δN=[0,0.01,0.05,0.1,0.2,0.5] while 100 runs.
Noise intensity (δN) NMSE of signal recovery BPDN L0 L2-SL0 Lp-RLS RRSL0 0 2.553e-05 1.311e-07 4.563e-07 1.563e-07 2.338e-07 0.01 4.659e-04 3.054e-04 2.352e-04 7.930e-05 6.032e-05 0.05 1.607e-03 4.079e-03 9.012e-04 2.384e-04 1.150e-04 0.1 7.633e-03 8.424e-03 3.309e-03 7.793e-04 5.522e-04 0.2 1.381e-02 2.090e-02 6.091e-03 1.130e-03 8.892e-04 0.5 3.063e-02 6.908e-02 1.083e-02 5.743e-03 3.041e-03Figure 4
Signal recovery effect by different algorithms when noise intensityδN = 0.2. (a) signal recovery by the BPDN algorithm; (b) signal recovery by the SL0 algorithm; (c) signal recovery by the L2-SL0 algorithm; (d) signal recovery by the Lp-RLS algorithm; (e) signal recovery by the RRSL0 algorithm.
(a) (b) (c) (d) (e)Figure 5
Reconstruction time-frequency characteristics of different algorithms when noise intensityδN = 0.2. (a) time-frequency characteristics of the original signal; (b) BPDN; (c) SL0; (d) L2-SL0; (e) Lp-RLS; (f) RRSL0.
(a) (b) (c) (d) (e) (f)


## 4.2. Comparison of Image Recovery Performance

Real images are considered to be approximately sparse under some proper basis, such as the Discrete Cosine Transform (DCT) basis and Discrete Wavelet Transform (DWT) basis. Here we choose DWT basis to recover these images. We compare the recovery performances based on the 2 real images in Figure6: Boat and Barbara, and the size of these images is 256×256. The compression ratio (CR, it is defined as m/n) is 0.5, and noise δN equals 0.01. We still choose SL0, BPDN, L2-SL0, and Lp-RLS to make comparisons. For image recovery, the object of image processing is given by(33)Y=ΦX+bHere Y, X, and b are matrices; among these, Y,b∈Rm×n, X∈Rn×n. In order to meet the basic requirements of CS, we make the following processing:(34)Y=Y1,Y2,…,Ynwhere Yi, i=1,2,…,n, is the column vector of Y. Therefore we can process images by CS.Figure 6
Original images.
(a)
Original Boat (b)
Original BarbaraFor performance of image recovery, we valuate it byPeakSignaltoNoiseRatio (PSNR) and StructuralSimilarityIndex (SSIM). PSNR is defined as(35)PSNR=10log⁡2552MSEwhere MSE=x-x^22, and SSIM is defined as(36)SSIMp,q=2μp+μq+c12σpq+c2μp2+μq2+c1σp2+σq2+c2.Among this,μp is the mean of image p, μq is the mean of image q, σp is the variance of image p, σq is the variance of image q, and σpq is covariance between image p and image q. Parameters c1=z1L and c2=z2L, where z1=0.01,z2=0.03, and L is dynamic range of pixel values. Range of SSIM is [-1,1]; when these two images are the same, the SSIM equals to 1.Figure6 shows the original images. Figure 7 shows the recovery effect of Boat and Barbara with noise intensity δN=0.01 and δN=0.1. In the Boat, the recovered images by SL0 and BPDN have the obvious water ripples while recovered images by other algorithms have no such water ripples. Similarly, for Barbara, recovered images by SL0 and BPDN are blurred compared with recovered images by other algorithms. L2-SL0, Lp-RLS, and RRSL0 algorithms are also effective for noisy image recovery. For L2-SL0, Lp-RLS, and RRSL0 algorithms, their recovery effects are very similar. In order to further analyze the advantages and disadvantages of these algorithms, we analyze PSNR and SSIM of the images recovered by these algorithms and results are shown in Tables 3–6. By observation and analysis, Lp-RLS performs better than L2-SL0 on the whole (except SSIM when δN=0.05, because the robustness of Lp-RLS is not strong); at the same time, RRSL0 outperforms Lp-RLS. Hence, the RRSL0 proposed by this paper is superior to other selected algorithms in image processing.Table 3
PSNR analysis of recovered Boat image by BPDN, SL0, L2-SL0,Lp-RLS, and RRSL0 with noise intensity changes according to sequence δN=[0,0.01,0.05,0.1,0.2,0.5].
Noise intensity (δN) PSNR of Boat image recovery BPDN SL0 L2-SL0 L p-RLS RRSL0 0 28.4928 31.5901 31.7175 32.2332 33.3495 0.01 27.9775 29.8589 30.2423 30.9742 31.5447 0.05 24.0531 27.1226 27.7101 28.2422 29.2505 0.1 20.3769 21.6485 23.9331 25.8639 26.4745 0.2 16.3044 17.9303 20.6744 23.0823 23.7415 0.5 9.2283 10.7863 13.4863 14.4974 16.3646Table 4
SSIM analysis of recovered Boat image by BPDN, SL0, L2-SL0,Lp-RLS, and RRSL0 with noise intensity changes according to sequence δN=[0,0.01,0.05,0.1,0.2,0.5].
Noise intensity (δN) SSIM of Boat image recovery BPDN SL0 L2-SL0 L p-RLS RRSL0 0 0.9865 0.9901 0.9902 0.9914 0.9934 0.01 0.9764 0.9852 0.9865 0.9887 0.9898 0.05 0.9418 0.9714 0.9756 0.9789 0.9827 0.1 0.8827 0.9083 0.9403 0.9631 0.9669 0.2 0.7452 0.8077 0.8868 0.9284 0.9396 0.5 0.3571 0.4427 0.5966 0.6492 0.7351Table 5
PSNR analysis of recovered Barbara image by BPDN, SL0, L2-SL0,Lp-RLS, and RRSL0 with noise intensity changes according to sequence δN=[0,0.01,0.05,0.1,0.2,0.5].
Noise intensity (δN) PSNR of Barbara image recovery BPDN SL0 L2-SL0 L p-RLS RRSL0 0 28.6613 29.5761 31.7921 32.3566 33.4074 0.01 27.7887 30.4746 30.9875 31.7365 32.6767 0.05 26.0316 26.1561 26.3147 27.3533 29.8613 0.1 20.7076 21.6914 23.8343 26.9408 27.5167 0.2 16.6564 16.6721 19.1297 19.2624 19.7801 0.5 9.1316 9.3525 13.1838 15.7972 16.9070Table 6
SSIM analysis of recovered Barbara image by BPDN, SL0, L2-SL0,Lp-RLS, and RRSL0 with noise intensity changes according to sequence δN=[0,0.01,0.05,0.1,0.2,0.5].
Noise intensity (δN) SSIM of Barbara image recovery BPDN SL0 L2-SL0 L p-RLS RRSL0 0 0.9834 0.9920 0.9926 0.9936 0.9984 0.01 0.9817 0.9904 0.9914 0.9927 0.9942 0.05 0.9797 0.9744 0.9753 0.9739 0.9888 0.1 0.9161 0.9314 0.9544 0.9728 0.9810 0.2 0.8165 0.8170 0.9051 0.9142 0.9278 0.5 0.3863 0.3992 0.6309 0.7712 0.8190Figure 7
Images recovery effect by BPDN, SL0, L2-SL0,Lp-RLS, and RRSL0 algorithms under different noise intensities. In subfigures (a), (b), (c), and (d), from left to right are images recovered by BPDN, SL0, L2-SL0, Lp-RLS, and RRSL0 algorithms.
(a)
Recovered Boat with noise intensityδN=0.01 (b)
Recovered Barbara with noise intensityδN=0.01 (c)
Recovered Boat with noise intensityδN=0.1 (d)
Recovered Barbara with noise intensityδN=0.1


## 5. Conclusions

In this paper, we propose the RRSL0 algorithm to recover sparse signal from given{y,Φ} in noise case. For the RRSL0 algorithm, it is constructed under the framework of Tikhonov regularization, in which the constraint condition y-Φx22≤ϵ is the deviation term, and the regularization term is replaced by a reweighted smoothed function WTHσ(x). As a key part of the RRSL0 algorithm, the reweighted smoothed function is to promote sparsity and provide guarantee for robust and accurate signal recovery. Furthermore, We consider the value of λ and the initial value σmax and final value σmin. Combined with previous literatures, we used α-method to determine the value of λ and deduced the method of initial value and final value of σ. Sparse signal recovery experiments on both the simulated signal and real images show the proposed RRSL0 algorithm performs better than the L1 or Lp regularization methods and classical L0 regularization methods. In addition, we would also like to apply the proposed algorithm to other CS applications such as the RPCA [36, 37] and SAR imaging [38].



---

*Source: 1978154-2019-07-17.xml*