# Sparse Signal Reconstruction Based on Multiparameter Approximation Function with Smoothedl0 Norm

Authors: Xiao-Feng Fang; Jiang-She Zhang; Ying-Qi Li
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The smoothedl0 norm algorithm is a reconstruction algorithm in compressive sensing based on approximate smoothed l0 norm. It introduces a sequence of smoothed functions to approximate the l0 norm and approaches the solution using the specific iteration process with the steepest method. In order to choose an appropriate sequence of smoothed function and solve the optimization problem effectively, we employ approximate hyperbolic tangent multiparameter function as the approximation to the big “steep nature” in l0 norm. Simultaneously, we propose an algorithm based on minimizing a reweighted approximate l0 norm in the null space of the measurement matrix. The unconstrained optimization involved is performed by using a modified quasi-Newton algorithm. The numerical simulation results show that the proposed algorithms yield improved signal reconstruction quality and performance.


---

## Body

## 1. Introduction

Sparse representation of signals has been extensively studied for decades. The concept of signal sparsity andl1 norm based on recovery techniques can be traced back to the work of Logan [1] in 1965, Santosa and Symes [2] in 1986, and Donoho and Stark in 1989 [3]. It is generally agreed that the foundation of the current state of the art in compressive sensing (CS) theory was laid by Candés et al. [4], Donoho [5], and Candés and Tao [6] in 2006. Mathematically, under sparsity assumptions, one would want to recover a signal x∈Rn, for example, the coefficient sequence of the signal in the appropriate basis, by solving the combinatorial optimization problem:
(1)min⁡x∥x∥0s.t.y=Ψx,
where y is an observed signal, Ψ is the dictionary matrix, x means a coefficient vector for the linear combination, and ∥x∥0 is the l0 norm of the vector x indicating the number of its nonzero elements.Exactly solving (1) is reported to be NP-hard [5, 7]; thus a relaxed approach [8–10] was proposed to address this challenge. By replacing nonconvex l0 norm with convex l1 norm, the method finds optimal solution using convex optimization [8, 11–13]. It was also introduced [6, 14–17] that, for many problems, the minimizing of the l1 norm is equivalent to that of the l0 norm under certain conditions. However, experimental results in [13] posed a strong question about the equivalence of minimizing both norms in practical problems. The sparse signal reconstruction (SSR) algorithm based on the optimization of a smoothed approximate l0 norm is studied in [11] where simulation results are compared with corresponding results obtained with several existing SSR algorithms with respect to reconstruction performance and computational complexity. These results favor the use of the approximate l0 norm.In this paper, we present a new signal reconstruction algorithm for CS, which is based on the minimization of a smoothed approximatel0 norm. But it differs from the previous algorithms in several aspects. First, we use approximation hyperbolic tangent function to approximate the l0 norm, which is found to have better approximation quality of l0 norm than the standard Gauss function. Second, the l0 norm minimization algorithm in this paper is carried out in null space of the measurement matrix, in which the constraints on measurements are eliminated and become unconstrained. Third, we use a reweighting technique in the null space of the measurement matrix to force the algorithm to reach the desired sparse solution faster. The rest of the paper is organized as follows. In Section 2, the smoothed l0 norm algorithm is briefly presented. Section 3 is about the signal reconstruction on the base of the reweighted approximation smoothed l0 norm algorithm. Section 4 provides the experimental results. The paper is finally concluded in Section 5.


## 2. Smoothedl0 Norm Algorithm (SL0)

Smoothedl0 norm attempts to solve the problem in (1) by approximating the l0 norm with a continuous function. Consider the continuous Gaussian function fσ(x) with the parameter σ:
(2)fσ(x)=e-x2/2σ2,x∈R,σ∈R+.The parameterσ determines the quality of the approximation and for any family of functions fσ(x) which approximates the Kronecker delta function.Define the continuous multivariate functiong(x) as
(3)g(x)=∑i=1N‍fσ(xi),x∈RN×1.
It follows from (2) and (3) that σ→0 and
(4)lim⁡σ→0fσ(x)={1x=00x≠0.Since the number of entries inx is N and the function g(x) is an indicator of the number of zero entries in x, the l0 norm of reconstructed vector x is approximated by
(5)∥x∥0≈N-g(x).Substituting this approximation into (1) yields the problem
(6)minimise∑i=1N‍(1-e-xi2/2σ2)s.t.y=Ψx.The approach is then to solve the problem (6) for a decreasing sequence of σ’s. The underlying thought is to select a σ which ensures that the initial solution is in the subset of RN×1 over which the approximation is convex and gradually increases the accuracy of the approximation. Here, the final smoothed l0 algorithm is proved to be faster with the possibility of recovering a sparser solution [11].Algorithm 1.
The final smoothedl0 norm (SL0 algorithm) is as follows.(1)
Initialization: choose an arbitrary solution from the feasible set, for example, the minimuml2 norm solution of x=Ψx; that is,
(7)x=ΨT(ΨΨT)-1y.(2)
Choose a decreasing sequence forσi, i=1,2,…,K.(3)
Fork=1,2,…,K, let σ=σk. Minimize function (6) using L iterations of the steepest ascent algorithm.
(a)
LetΔx=[x1exp⁡(-x12/2σk2),…,xnexp⁡(-xn2/2σk2)]T.(b)
Letx←x-μΔx (where μ is a small position constant); then
(8)x⟵x-ΨT(ΨΨT)-1y.(4)
Final answer isx.


## 3. Signal Reconstruction Based on Approximation Function

### 3.1. Approximatel0 Norm with Approximate Hyperbolic Tangent Sequence

The key to SL0 algorithm is to select high quality smoothed continuous function to approximatel0 norm. The optimal value of l0 norm is achieved through the solution of the minimum value of the approximate smoothed continuous function. In general, the approximation of the l0 norm with Gauss function cannot produce desired result. To effectively approximate l0 norm, we analyzed inverse tangent function [18], hyperbolic tangent function, and approximate hyperbolic tangent sequence.Inverse tangent function is as follows:(9)hσ(xi)=2πarctan(xi22σ2).
Hyperbolic tangent function is as follows:
(10)gσ(xi)=exi2/2σ2-e-xi2/2σ2exi2/2σ2+e-xi2/2σ2.
Approximate hyperbolic tangent function is as follows:
(11)fσ(xi)=18exi2/2σ2-e-xi2/2σ2exi2/2σ2+e-xi2/2σ2+78(1-e-xi2/((1/8)σ2)),
where σ is a parameter and xi is a component of the vector x. The following properties can be obtained.(1) For formula (11), let Fσ(x)=∑i=1N‍fσ(xi); then ∥x∥0=lim⁡σ→0Fσ(x) since
(12)lim⁡σ→0fσ(xi)=lim⁡σ→0(18·exi2/2σ2-e-xi2/2σ2exi2/2σ2+e-xi2/2σ2+78·(1-e-xi2/((1/8)σ2)))={0xi=01xi≠0.
According to the definition of l0 norm, when σ goes infinitely to zero, the value of Fσ(x) goes as well to the nonzero elements of the vector x, that is, the l0 norm of vector x. Therefore, it can be shown with the smoothed function ∥x∥0=lim⁡σ→0Fσ(x). Obviously, formulas (9) and (10) have the same approximation property of l0 norm as that of (11). The approximation performance of l0 norm of each formula, however, is different. For further illustration of the advantages of approximate hyperbolic tangent function as smoothed continuous function, we compared the distribution of each approximation function at interval [-1,1] when σ=0.1. Figure 1 is the approximation quality comparison of the smoothed l0 norm of the standard Gauss function, inverse tangent function, hyperbolic tangent function, and approximate hyperbolic tangent function. The estimation of l0 norm by approximate hyperbolic tangent function outsmarted other functions. It shows that the approximate hyperbolic tangent function has better “steep nature" between the intervals of [-0.5,0.5] and the estimation of the l0 norm is more accurate.Figure 1
Comparison of the four approximatel0 norm functions with σ=0.1.(2) Givenσ≥0, approximate hyperbolic tangent function fσ(xi) has better convergence than other approximate functions; that is, fσ(xi)≥(1-φ(xi)), for arbitrary σ, φ, where φ(xi)=e-xi2/2σ2.Proof.
Letu(x)=f(x)-(1-φ(x)), when x=0 or x=∞, and u(x)=0, when 0<|x|<∞,
(13)u(x)=18·ex2/2σ2-e-x2/2σ2ex2/2σ2+e-x2/2σ2+78·(1-e-x2/((1/8)σ2))-1+e-x2/2σ2=18·(1-2e-x2/2σ2ex2/2σ2+e-x2/2σ2)+78·(1-e-x2/((1/8)σ2))-1+e-x2/2σ2=-18·2e-x2/2σ2ex2/2σ2+e-x2/2σ2-78·e-x2/((1/8)σ2)+18·e-x2/2σ2+78·e-x2/2σ2=18·(-2e-x2/2σ2ex2/2σ2+e-x2/2σ2+e-x2/2σ2)+78·(e-x2/2σ2-e-x2/((1/8)σ2))=18·e-x2/2σ2ex2/2σ2+e-x2/2σ2·(-2+ex2/2σ2+e-x2/2σ2)+78·(e-x2/2σ2-e-x2/((1/8)σ2))≥18·11+ex2/σ2(ex2/4σ2+e-x2/4σ2)2≥0.From (9), (10), and (11), it is clear that the smaller the σ is, the more the local extreme values of objective function would be, and it would be difficult to get the global optimal value of the objective function. The σ determines the smoothness of the objective function. The bigger the σ, the smoother the objective function, and the less accurate it would be to estimate the l0 norm. Conversely, the smaller the σ is, the better the approximation of the l0 norm would be. In the process, we built a decreasing sequence σ1,σ2,…,σn to optimize each corresponding objective function until σi is small enough, so as to eliminate the influence of the local extreme value and get the global optimal value of the smoothed function. We, thus, use approximate hyperbolic tangent sequence for the approximation of l0 norm.Problem (1) can be approximated by the following model:
(14)min⁡Fσ(x)s.t.y=Φx.(3) There are lots of algorithms to the solving of problem (14), among which the most representative one is the steepest descent method. But there still exist two drawbacks.(a)
“Notched effect" in search direction strongly hinders the convergence speed.(b)
Search steps cannot be estimated. Usually, it is estimated with experiences, which lacks theoretical support.We, therefore, propose the modified quasi-Newton algorithm to solve the above problem. First, calculate the Newton direction of the approximate hyperbolic tangent function:(15)d=-∇2Fσ(x)-1∇Fσ(x),
where
(16)∇Fσ(x)=[18·4x1/σ2(ex12/2σ2+e-x12/2σ2)2+78·16x1σ2e-8x12/σ2,…,18·4xN/σ2(exN2/2σ2+e-xN2/2σ2)2+78·16xNσ2e-8xN2/σ24x1/σ2(ex12/2σ2+e-x12/2σ2)2]T,(17)∇2(Fσ(x))=[∂2fσ(x1)∂x120⋯00∂2fσ(x2)∂x22⋯0⋮⋮⋮⋮00⋯∂2fσ(xN)∂xN2],(18)∂2fσ(xk)∂xk2=18·(4/σ2)(1+(xk2/σ2))(exk2/2σ2+e-xk2/2σ2)2+78·16σ2σ-8xk2/σ2(1-16xk2σ2).
We know that the Hesse matrix is nonpositive definite in Newton direction d, which cannot ensure that the Newton direction is descent direction, so the Hesse matrix should be modified. Then, set up a new matrix:
(19)G=∇2(Fσ(x))+εkI,
where εk is a group of proper positive numbers, I is identity matrix, and the diagonal elements in matrix G are all positive. Since
(20)lim⁡xk→0(exk2/2σ2+e-xk2/2σ2)2=4,
choose
(21)εk=78·162xk2σ4e-8xk2/σ2,
as modification coefficients; then the diagonal elements in matrix G can be shown as
(22)18·4(1+(xk2/σ2))σ2(exk2/2σ2+e-xk2/2σ2)2+78·16σ2e-8xk2/σ2.
Also, it can be obtained that
(23)d=-∇2Fσ(x)-1·∇Fσ(x)=[-σ2x1+28·4σ2x1e-8x12/σ2σ2+x12+28·4σ2e-8x12/σ2,…,-σ2xN+28·4σ2xNe-8xN2/σ2σ2+xN2+28·4σ2e-8xN2/σ2]T.(4) Usually, parameterσ is chosen as σj=γ·σj-1, j=1,2,…, and J,γ∈(0.5,1). σ1 and σj are chosen as follows.Letx=max⁡i|xi0|, for faster algorithm convergence; parameter σ should meet the following:
(24)hσ1(x)=2πarctan(x22σ12)≤12⟹σ1≥x2,
and choosing σ1=max⁡i(|xi0|/2), when σj→0, we can obtain Fσj(x)→∥x∥0, so when σj is smaller, Fσj(x) can better reflect the sparsity of vector x. Meanwhile, however, it is more sensitive to noise; the value of σj should not be too small. In summary, the steps are as follows.Algorithm 2.
The approximate smoothedl0 norm (ASL0 algorithm) is as follows. (1)
Initialization:x=ΦTy, initial residual r0=0, and smoothed function’s initial parameter σ=1.(2)
Calculate the modified Newton direction with (23) and use modified Newton algorithm to get x=x+d.(3)
Use gradient projection to havex=x-ΦT(ΦΦT)-1(Φx-y), and calculate the residual r=y-Φx.(4)
If the residual of the adjacent iterations meets∥r-r0∥≤ε, then stop. Otherwise, let r0=r, continue, and return to (2).(5)
Further decreaseσ, let σ=γ·σ(0<γ<1), usually set γ=0.7, and return to (2). Otherwise, stop iterating and get the optimal value x.


### 3.2. Reweighting the Approximatel0 Norm Minimization

It is well known that all solutions ofΨx=y can be parameterized as
(25)x=xs+Vnξ,
where xs is a solution of Ψx=y and Vn is a matrix whose columns constitute an orthogonal basis of the null space of Ψ [19]. Vector xs and matrix Vn can be evaluated by using the singular-value decomposition (SVD) or by using the QR decomposition of matrix Ψ [20]. Using (25) the problem in (14) is reduced to
(26)minimizeξ∑i=1N‍fσ(xsi+ViTξ),
where
(27)fσ(xsi+ViTξ)=18·e[xsi+viTξ]2/2σ2-e-[xsi+viTξ]2/2σ2e[xsi+viTξ]2/2σ2+e-[xsi+viTξ]2/2σ2+78·(1-e-[xsi+viTξ]2/((1/8)σ2)).Signal reconstruction based on the solution of the problem in (28) works well, but the technique can be considerably enhanced by incorporating a reweighting strategy, and the reweighted unconstrained problem can be formulated as follows [21]:
(28)minimizeξ∑i=1N‍wi·fσ(xsi+ViTξ),
where wi are positive scalars that are from a weight vector w=[w1w2⋯wN]. Starting with an initial w0=eN, where eN is the all-one vector of dimension N, in the (k+1)th iteration the weight vector is updated to w(k+1) with its ith component given by
(29)wi(k+1)=1|xi(k)|+τ,
where xi(k) denotes the ith component of vector x(k) obtained in the kth iteration as x(k)=xs+Vnξ(k) and τ is a small positive scalar which is used to prevent numerical instability when |xi(k)| approaches zero. As long as σ>0, the objective function in (28) remains differentiable and its gradient can be obtained as
(30)gi(x)=wi(18·4[xsi+viTξ]σ2+78·16[xsi+viTξ]σ2e-8[xsi+viTξ]2/σ2).
For a fixed value of σ, the problem in (26) is solved by using a quasi-Newton algorithm where an approximation of the inverse of the Hessian is obtained by using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) update formula. It can be shown that function (28) remains convex in the region where the largest magnitude of the components of x=xl2+Vnξ is less than σ, where vector xs is chosen to be the least squares solution xl2 of y=Ψx; namely, xs=xl2=ΨT(ΨΨT)-1y. Based on this, a reasonable initial value of σ can be chosen as σ0=max⁡|xl2|+τ, where τ is a small positive scalar. As the algorithm starts at the origin ξ(0)=0, the above choice of σ0 ensures that the optimization starts in a convex region. This greatly facilitates the convergence of the proposed algorithm.Algorithm 3.
The reweighted approximate smoothedl0 norm (RASL0 algorithm) is as follows.(1)
Initialization: inputΨ, xl2, σT, r, τ, and ε, and set ξ(0)=0, w(0)=eN, σ=max⁡|xl2|+τ, and k=0.(2)
Perform theQR decomposition of ΨT=QR and construct Vn using the last columns of Q.(3)
Withw=w(k) and ξ(0) as an initial point, apply the BFGS algorithm to solve the problem in (28), where reweighting with parameter ε is applied using (29) in each iteration; denote the solution by ξ(k).(4)
Computex(k)=xl2+Vnξ(k) and update weight vector to w(k+1) with (29).(5)
Ifσ≤σT, stop and output x(k) as the solution; otherwise set ξ(0)=ξ(k), σ=r·σ, k=k+1, and repeat from step (3).


## 3.1. Approximatel0 Norm with Approximate Hyperbolic Tangent Sequence

The key to SL0 algorithm is to select high quality smoothed continuous function to approximatel0 norm. The optimal value of l0 norm is achieved through the solution of the minimum value of the approximate smoothed continuous function. In general, the approximation of the l0 norm with Gauss function cannot produce desired result. To effectively approximate l0 norm, we analyzed inverse tangent function [18], hyperbolic tangent function, and approximate hyperbolic tangent sequence.Inverse tangent function is as follows:(9)hσ(xi)=2πarctan(xi22σ2).
Hyperbolic tangent function is as follows:
(10)gσ(xi)=exi2/2σ2-e-xi2/2σ2exi2/2σ2+e-xi2/2σ2.
Approximate hyperbolic tangent function is as follows:
(11)fσ(xi)=18exi2/2σ2-e-xi2/2σ2exi2/2σ2+e-xi2/2σ2+78(1-e-xi2/((1/8)σ2)),
where σ is a parameter and xi is a component of the vector x. The following properties can be obtained.(1) For formula (11), let Fσ(x)=∑i=1N‍fσ(xi); then ∥x∥0=lim⁡σ→0Fσ(x) since
(12)lim⁡σ→0fσ(xi)=lim⁡σ→0(18·exi2/2σ2-e-xi2/2σ2exi2/2σ2+e-xi2/2σ2+78·(1-e-xi2/((1/8)σ2)))={0xi=01xi≠0.
According to the definition of l0 norm, when σ goes infinitely to zero, the value of Fσ(x) goes as well to the nonzero elements of the vector x, that is, the l0 norm of vector x. Therefore, it can be shown with the smoothed function ∥x∥0=lim⁡σ→0Fσ(x). Obviously, formulas (9) and (10) have the same approximation property of l0 norm as that of (11). The approximation performance of l0 norm of each formula, however, is different. For further illustration of the advantages of approximate hyperbolic tangent function as smoothed continuous function, we compared the distribution of each approximation function at interval [-1,1] when σ=0.1. Figure 1 is the approximation quality comparison of the smoothed l0 norm of the standard Gauss function, inverse tangent function, hyperbolic tangent function, and approximate hyperbolic tangent function. The estimation of l0 norm by approximate hyperbolic tangent function outsmarted other functions. It shows that the approximate hyperbolic tangent function has better “steep nature" between the intervals of [-0.5,0.5] and the estimation of the l0 norm is more accurate.Figure 1
Comparison of the four approximatel0 norm functions with σ=0.1.(2) Givenσ≥0, approximate hyperbolic tangent function fσ(xi) has better convergence than other approximate functions; that is, fσ(xi)≥(1-φ(xi)), for arbitrary σ, φ, where φ(xi)=e-xi2/2σ2.Proof.
Letu(x)=f(x)-(1-φ(x)), when x=0 or x=∞, and u(x)=0, when 0<|x|<∞,
(13)u(x)=18·ex2/2σ2-e-x2/2σ2ex2/2σ2+e-x2/2σ2+78·(1-e-x2/((1/8)σ2))-1+e-x2/2σ2=18·(1-2e-x2/2σ2ex2/2σ2+e-x2/2σ2)+78·(1-e-x2/((1/8)σ2))-1+e-x2/2σ2=-18·2e-x2/2σ2ex2/2σ2+e-x2/2σ2-78·e-x2/((1/8)σ2)+18·e-x2/2σ2+78·e-x2/2σ2=18·(-2e-x2/2σ2ex2/2σ2+e-x2/2σ2+e-x2/2σ2)+78·(e-x2/2σ2-e-x2/((1/8)σ2))=18·e-x2/2σ2ex2/2σ2+e-x2/2σ2·(-2+ex2/2σ2+e-x2/2σ2)+78·(e-x2/2σ2-e-x2/((1/8)σ2))≥18·11+ex2/σ2(ex2/4σ2+e-x2/4σ2)2≥0.From (9), (10), and (11), it is clear that the smaller the σ is, the more the local extreme values of objective function would be, and it would be difficult to get the global optimal value of the objective function. The σ determines the smoothness of the objective function. The bigger the σ, the smoother the objective function, and the less accurate it would be to estimate the l0 norm. Conversely, the smaller the σ is, the better the approximation of the l0 norm would be. In the process, we built a decreasing sequence σ1,σ2,…,σn to optimize each corresponding objective function until σi is small enough, so as to eliminate the influence of the local extreme value and get the global optimal value of the smoothed function. We, thus, use approximate hyperbolic tangent sequence for the approximation of l0 norm.Problem (1) can be approximated by the following model:
(14)min⁡Fσ(x)s.t.y=Φx.(3) There are lots of algorithms to the solving of problem (14), among which the most representative one is the steepest descent method. But there still exist two drawbacks.(a)
“Notched effect" in search direction strongly hinders the convergence speed.(b)
Search steps cannot be estimated. Usually, it is estimated with experiences, which lacks theoretical support.We, therefore, propose the modified quasi-Newton algorithm to solve the above problem. First, calculate the Newton direction of the approximate hyperbolic tangent function:(15)d=-∇2Fσ(x)-1∇Fσ(x),
where
(16)∇Fσ(x)=[18·4x1/σ2(ex12/2σ2+e-x12/2σ2)2+78·16x1σ2e-8x12/σ2,…,18·4xN/σ2(exN2/2σ2+e-xN2/2σ2)2+78·16xNσ2e-8xN2/σ24x1/σ2(ex12/2σ2+e-x12/2σ2)2]T,(17)∇2(Fσ(x))=[∂2fσ(x1)∂x120⋯00∂2fσ(x2)∂x22⋯0⋮⋮⋮⋮00⋯∂2fσ(xN)∂xN2],(18)∂2fσ(xk)∂xk2=18·(4/σ2)(1+(xk2/σ2))(exk2/2σ2+e-xk2/2σ2)2+78·16σ2σ-8xk2/σ2(1-16xk2σ2).
We know that the Hesse matrix is nonpositive definite in Newton direction d, which cannot ensure that the Newton direction is descent direction, so the Hesse matrix should be modified. Then, set up a new matrix:
(19)G=∇2(Fσ(x))+εkI,
where εk is a group of proper positive numbers, I is identity matrix, and the diagonal elements in matrix G are all positive. Since
(20)lim⁡xk→0(exk2/2σ2+e-xk2/2σ2)2=4,
choose
(21)εk=78·162xk2σ4e-8xk2/σ2,
as modification coefficients; then the diagonal elements in matrix G can be shown as
(22)18·4(1+(xk2/σ2))σ2(exk2/2σ2+e-xk2/2σ2)2+78·16σ2e-8xk2/σ2.
Also, it can be obtained that
(23)d=-∇2Fσ(x)-1·∇Fσ(x)=[-σ2x1+28·4σ2x1e-8x12/σ2σ2+x12+28·4σ2e-8x12/σ2,…,-σ2xN+28·4σ2xNe-8xN2/σ2σ2+xN2+28·4σ2e-8xN2/σ2]T.(4) Usually, parameterσ is chosen as σj=γ·σj-1, j=1,2,…, and J,γ∈(0.5,1). σ1 and σj are chosen as follows.Letx=max⁡i|xi0|, for faster algorithm convergence; parameter σ should meet the following:
(24)hσ1(x)=2πarctan(x22σ12)≤12⟹σ1≥x2,
and choosing σ1=max⁡i(|xi0|/2), when σj→0, we can obtain Fσj(x)→∥x∥0, so when σj is smaller, Fσj(x) can better reflect the sparsity of vector x. Meanwhile, however, it is more sensitive to noise; the value of σj should not be too small. In summary, the steps are as follows.Algorithm 2.
The approximate smoothedl0 norm (ASL0 algorithm) is as follows. (1)
Initialization:x=ΦTy, initial residual r0=0, and smoothed function’s initial parameter σ=1.(2)
Calculate the modified Newton direction with (23) and use modified Newton algorithm to get x=x+d.(3)
Use gradient projection to havex=x-ΦT(ΦΦT)-1(Φx-y), and calculate the residual r=y-Φx.(4)
If the residual of the adjacent iterations meets∥r-r0∥≤ε, then stop. Otherwise, let r0=r, continue, and return to (2).(5)
Further decreaseσ, let σ=γ·σ(0<γ<1), usually set γ=0.7, and return to (2). Otherwise, stop iterating and get the optimal value x.


## 3.2. Reweighting the Approximatel0 Norm Minimization

It is well known that all solutions ofΨx=y can be parameterized as
(25)x=xs+Vnξ,
where xs is a solution of Ψx=y and Vn is a matrix whose columns constitute an orthogonal basis of the null space of Ψ [19]. Vector xs and matrix Vn can be evaluated by using the singular-value decomposition (SVD) or by using the QR decomposition of matrix Ψ [20]. Using (25) the problem in (14) is reduced to
(26)minimizeξ∑i=1N‍fσ(xsi+ViTξ),
where
(27)fσ(xsi+ViTξ)=18·e[xsi+viTξ]2/2σ2-e-[xsi+viTξ]2/2σ2e[xsi+viTξ]2/2σ2+e-[xsi+viTξ]2/2σ2+78·(1-e-[xsi+viTξ]2/((1/8)σ2)).Signal reconstruction based on the solution of the problem in (28) works well, but the technique can be considerably enhanced by incorporating a reweighting strategy, and the reweighted unconstrained problem can be formulated as follows [21]:
(28)minimizeξ∑i=1N‍wi·fσ(xsi+ViTξ),
where wi are positive scalars that are from a weight vector w=[w1w2⋯wN]. Starting with an initial w0=eN, where eN is the all-one vector of dimension N, in the (k+1)th iteration the weight vector is updated to w(k+1) with its ith component given by
(29)wi(k+1)=1|xi(k)|+τ,
where xi(k) denotes the ith component of vector x(k) obtained in the kth iteration as x(k)=xs+Vnξ(k) and τ is a small positive scalar which is used to prevent numerical instability when |xi(k)| approaches zero. As long as σ>0, the objective function in (28) remains differentiable and its gradient can be obtained as
(30)gi(x)=wi(18·4[xsi+viTξ]σ2+78·16[xsi+viTξ]σ2e-8[xsi+viTξ]2/σ2).
For a fixed value of σ, the problem in (26) is solved by using a quasi-Newton algorithm where an approximation of the inverse of the Hessian is obtained by using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) update formula. It can be shown that function (28) remains convex in the region where the largest magnitude of the components of x=xl2+Vnξ is less than σ, where vector xs is chosen to be the least squares solution xl2 of y=Ψx; namely, xs=xl2=ΨT(ΨΨT)-1y. Based on this, a reasonable initial value of σ can be chosen as σ0=max⁡|xl2|+τ, where τ is a small positive scalar. As the algorithm starts at the origin ξ(0)=0, the above choice of σ0 ensures that the optimization starts in a convex region. This greatly facilitates the convergence of the proposed algorithm.Algorithm 3.
The reweighted approximate smoothedl0 norm (RASL0 algorithm) is as follows.(1)
Initialization: inputΨ, xl2, σT, r, τ, and ε, and set ξ(0)=0, w(0)=eN, σ=max⁡|xl2|+τ, and k=0.(2)
Perform theQR decomposition of ΨT=QR and construct Vn using the last columns of Q.(3)
Withw=w(k) and ξ(0) as an initial point, apply the BFGS algorithm to solve the problem in (28), where reweighting with parameter ε is applied using (29) in each iteration; denote the solution by ξ(k).(4)
Computex(k)=xl2+Vnξ(k) and update weight vector to w(k+1) with (29).(5)
Ifσ≤σT, stop and output x(k) as the solution; otherwise set ξ(0)=ξ(k), σ=r·σ, k=k+1, and repeat from step (3).


## 4. Simulation Result

This section describes some experiments testifying to the performances of the proposed algorithm. All experiments are performed under Windows 7 and MATLAB V8.0 (R2012a) running on a Lenovo workstation with an Intel (R) Core (TM), CPU i7 at 2.40 GHz, and 6 GB of memory.Simulation 1: Algorithm Performances Comparison of One-Dimensional Signal. In the experiment, the signal length and number of measurements were set to N=256 and M=100, respectively. A total of 16 sparse signals with sparsity K=4q-3, q=1,2,…,16, were used. Then a K-sparse signal x was constructed as follows.(1)
Setx to zero vector of length N.(2)
Generate a vectoru of length K assuming that each component ui is a random value drawn from a normal distribution N(0,1).(3)
Randomly selectK indices from the set 1,2,…,N, and set xi1=u1, and xi2=u2,…,xiK=uK. The measurement matrix was assumed to be of size M×N and was generated by drawing its elements from N(0,1), followed by a normalization step to ensure that the l2 norm of each column is unity.The performance of iteratively reweighted (IR) algorithm [22, 23] with p=0.1, the SL0 algorithm, ASL0 algorithm and RASL0 algorithm with σT=10-4, r=1/3, τ=0.01, and ϵ=0.08 was measured in terms of the number of perfect reconstructions over 100 runs. The results obtained were plotted in Figure 2. It can be observed that the RASL0 algorithm outperforms the IR algorithm. On comparing the SL0 algorithm and the ASL0 algorithm, we note that the two algorithms are comparable for K smaller than 45 but the RASL0 algorithm performs better for K larger than 45.Figure 2
Number of perfect reconstructions by each algorithm over 100 runs withN=256 and M=100.The mathematical complexity of the SL0 algorithm, IR algorithm, ASL0 algorithm, and RASL0 algorithm was measured in terms of average CPU time over 100 runs for typical instances withM=N/2 and K=round(M/2.5) where N varies in the range between 128 and 512. In Figure 3, it is noted that the RASL0 algorithm and the SL0 algorithms were more efficient than the IR algorithm, and the complexity of the RASL0 algorithm was slightly higher than that of the SL0 algorithm. The moderate increase in the mathematical complexity of the RASL0 algorithm was primarily due to the fact that the objective function in (26) needed to be modified in each iteration using (29). Typically the RASL0 algorithm converges in a small number of iterations. Figure 4 showed how the objective function in (26) converges in 24 iterations, where the parameters were set to N=256, M=100, K=40, and σ=0.0212.Figure 3
Average CPU time required by each algorithm over 100 runs withM=N/2 and K=round(M/2.5).Figure 4
FunctionFσ(ξ) for N=256, M=100, K=40, and σ=0.0212.Simulation 2: Algorithm Performances Comparison of Two-Dimensional Signal. As to examine the algorithms’ performances for different images, the reconstruction performance of each method is evaluated in terms of two indicators: the reconstruction relative error (RE) and peak signal-to-noise ratio (PSNR). The relative error of reconstructed image is defined as follows:
(31)RE(x)=∥x-x¯∥2∥x¯∥2,
where x¯ and x are the original and the reconstructed image, respectively. Apparently, the lower the value of relative error is, the better the reconstructed performance will be. The PSNR of reconstructed images is defined by
(32)PSNR(x)=10log⁡10(m×n×2552∥x-x¯∥2),
where m, n are the size of the image.In experiment, the DCT basis is selected as the sparse representation dictionary; compression ratio isM/N=0.25. Figure 5 shows that the quality of the reconstruction with RASL0 algorithm has been improved and details of images are better reconstructed.The quality of the reconstruction with different algorithms in different images. Original image.
(a)
Original image
(b)
SL0
(c)
ASL0
(d)
RASL0
(e)
Original image
(f)
SL0
(g)
ASL0
(h)
RASL0Table1 shows three different 256×256 images, with the same compression ratio M/N=0.5. From the comparison of the PSNR, the RE, the signal-to-noise ratio (SNR), and the matching degree (MD) among RASL0 algorithm, SL0 algorithm, and ASL0 algorithm, it is clear that the relative errors with RASL0 algorithm are 0.02% smaller than the SL0 algorithm, the SNR ratio is 3 dB higher, the PSNR is improved by 2 db, and the matching degree is also improved.Table 1
Reconstruction quality comparison with different images and different approximation smoothedl0 algorithms.
Image
Algorithm
PSNR (dB)
RE
SNR (dB)
MD
Lena
SL0
29.6238
0.0947
34.3124
0.9831
ASL0
32.2244
0.0742
36.2702
0.9893
RASL0
36.9566
0.0578
39.9432
0.9926
Camera
SL0
26.4476
0.0921
30.6337
0.9846
ASL0
29.0284
0.0747
32.2134
0.9873
RASL0
33.3116
0.0634
36.2727
0.9912
Boat
SL0
28.1573
0.0852
31.3023
0.9813
ASL0
29.5001
0.0786
33.7382
0.9865
RASL0
30.5814
0.0614
34.3743
0.9908Figure6 is the reconstruction effect comparison of image (Lena) with different algorithms. It shows that RASL0 outperforms OMP, GP, GPSR, and NSL0 algorithms. Table 2 is the quality comparison of different algorithms, indicating that RASL0 algorithm is better in all aspects.Table 2
Reconstruction quality comparison with four other algorithms.
Algorithm
PSNR (dB)
RE
Runtime (s)
OMP
30.5951
0.0812
48.32
GP
30.7121
0.0803
12.23
SP
31.2486
0.0702
8.82
GPSR
34.2154
0.0656
11.38
RASL0
36.9566
0.0578
13.36Reconstruction quality comparison with four other algorithms.
(a)
Original image
(b)
OMP
(c)
GP
(d)
SP
(e)
GPSR
(f)
RASL


## 5. Conclusion and Future Work

In this paper, we use approximate hyperbolic tangent function to estimatel0 norm and design a RSNL0 algorithm based on minimizing an approximate l0 norm of the signal in the null space of the measurement matrix, where a reweighting technique is used to force the solution’s sparsity and a quasi-Newton algorithm is used to accelerate the optimization. Simulation results are presented which demonstrate that the proposed algorithm yields improved signal reconstruction performance and requires a reduced amount of computation relative to iteratively reweighted algorithms based on the lp norm, with p<1. When compared with a known algorithm based on a smoothed l0 norm, improved signal reconstruction is achieved although the amount of computation is increased somewhat. However, how to choose the parameter decreasing sequence so as to eliminate the disturbance of the minimal value to the solution of the minimum value would be our future work.



---

*Source: 416542-2014-06-04.xml*