New Nonsmooth Equations-Based Algorithms for ℓ₁-Norm Minimization and Applications

Abstract

Recently, Xiao et al. proposed a nonsmooth equations-based method to solve the ℓ₁-norm minimization problem (2011). The advantage of this method is its simplicity and lower storage. In this paper, based on new nonsmooth equations reformulation, we investigate new nonsmooth equations-based algorithms for solving ℓ₁-norm minimization problems. Under mild conditions, we show that the proposed algorithms are globally convergent. The preliminary numerical results demonstrate the effectiveness of the proposed algorithms.

1. Introduction

We consider the ℓ₁-norm minimization problem min_𝑥𝑓(𝑥)≜12‖𝐴𝑥−𝑏‖²+𝜌‖𝑥‖₁, (1.1) where 𝑥∈ℝ^𝑛, 𝑏∈ℝ^𝑚, 𝐴∈ℝ^𝑚×𝑛, and 𝜌 is a nonnegative parameter. Throughout the paper, we use ‖𝑣‖=∑^𝑛_𝑖=1|𝑣_𝑖|² and ‖𝑣‖₁=∑_𝑖|𝑣_𝑖| to denote the Euclidean norm and the ℓ₁-norm of vector 𝑣∈𝑅^𝑛, respectively. Problem (1.1) has many important practical applications, particularly in compressed sensing (abbreviated as CS) [1] and image restoration [2]. It can also be viewed as a regularization technique to overcome the ill-conditioned, or even singular, nature of matrix 𝐴, when trying to infer 𝑥 from noiseless observations 𝑏=𝐴𝑥 or from noisy observations 𝑏=𝐴𝑥+𝜉, where 𝜉 is the white Gaussian noise of variance 𝜎² [3–5].

The convex optimization problem (1.1) can be cast as a second-order cone programming problem and thus could be solved via interior point methods. However, in many applications, the problem is not only large scale but also involves dense matrix data, which often precludes the use and potential advantage of sophisticated interior point methods. This motivated the search of simpler first-order algorithms for solving (1.1), where the dominant computational effort is a relatively cheap matrix-vector multiplication involving 𝐴 and 𝐴^𝑇. In the past few years, several first-order algorithms have been proposed. One of the most popular algorithms falls into the iterative shrinkage/thresholding (IST) class [6, 7]. It was first designed for wavelet-based image deconvolution problems [8] and analyzed subsequently by many authors, see, for example, [9–11]. Figueiredo et al. [12] studied the gradient projection and Barzilai-Borwein method [13] (denoted by GPSR-BB) for solving (1.1). They reformulated problem (1.1) as a box-constrained quadratic program and solved it by a gradient projection and Barzilai-Borwein method. Wright et al. [14] presented sparse reconstruction algorithm (denoted by SPARSA) to solve (1.1). Yun and Toh [15] proposed a block coordinate gradient descent algorithm for solving (1.1). Yang and Zhang [16] investigated alternating direction algorithms for solving (1.1).

Quite recently, Xiao et al. [17] developed a nonsmooth equations-based algorithm (called SGCS) for solving ℓ₁-norm minimization problems in CS. They reformulated the box-constrained quadratic program obtained by Figueiredo et al. [12] into a system of nonsmooth equations and then applied the spectral gradient projection method [18] to solving the nonsmooth equation. The main advantage of the SGCS is its simplicity and lower storage. The difference between the above algorithms and SGCS is that SGCS did not use line search to decrease the value of objective function at each iteration and instead used a projection step to accelerate the iterative process. However, each projection step in SGCS requires two matrix-vector multiplication involving 𝐴 or 𝐴^𝑇, which means that each iteration requires matrix-vector multiplication involving 𝐴 or 𝐴^𝑇 four times, while each iteration in GPSR-BB and IST is only two times. This may bring in more computational complexity. In addition, the dimension of the system of nonsmooth equations is 2𝑛, which is twice of the original problems. These drawbacks motivate us to study new nonsmooth equations-based algorithms for the ℓ₁-norm minimization problem.

In this paper, we first reformulate problem (1.1) into a system of nonsmooth equations. This system is Lipschitz continuous and monotone and many effective algorithms (see, e.g., [18–22]) can be used to solve it. We then apply spectral gradient projection (denoted by SGP) method [18] to solve the resulting system. Similar to SGCS, each iteration in SGP requires matrix-vector multiplication involving 𝐴 or 𝐴^𝑇 four times. In order to reduce the computational complexity, we also propose a modified SGP (denoted by MSGP) method to solve the resulting system. Under mild conditions, the global convergence of the proposed algorithms will be ensured.

The remainder of the paper is organized as follows. In Section 2, we first review some existing results of nonsmooth analysis and then derive an equivalent system of nonsmooth equations to problem (1.1). We verify some nice properties of the resulting system in this section. In Section 3, we propose the algorithms and establish their global convergence. In Section 4, we apply the proposed algorithms to some practical problems arising from compressed sensing and image restoration and compare their performance with that of SGCS, SPARSA, and GPSR-BB.

Throughout the paper, we use ⟨⋅,⋅⟩ to denote the inner product of two vectors in 𝑅^𝑛.

2. Preliminaries

By nonsmooth analysis, a necessary condition for a vector 𝑥∈ℝ^𝑛 to be a local minima of nonsmooth function 𝑓∶ℝ^𝑛→ℝ is that (0,…,0)^𝑇∈𝜕𝑓(𝑥), (2.1) where 𝜕𝑓(𝑥) denotes the subdifferential of 𝑓 at 𝑥 [23]. If 𝑓 is convex, then (2.1) is also sufficient for 𝑥 to be a solution of (1.1). The subdifferential of the absolute value function |𝑡| is given by the signum function sign(𝑡), that is 𝜕|𝑡|(𝑡){1}[]{}=sign∶=,𝑡>0,−1,1,𝑡=0,−1,𝑡<0. (2.2) For problem (1.1), the optimality conditions therefore translate to ∇_𝑖𝑓(𝑥)+𝜌sign𝑥_𝑖=0,𝑥_𝑖∇>0,_𝑖𝑓(𝑥)≤𝜌,𝑥_𝑖=0, (2.3) where ∇_𝑖𝑓(𝑥)=𝜕𝑓(𝑥)/𝜕𝑥_𝑖, 𝑖=1,…,𝑛. It is clear that the function defined by (1.1) is convex. Therefore a point 𝑥^∗∈𝑅^𝑛 is a solution of problem (1.1) if and only if it satisfies ∇_𝑖𝑓𝑥^∗+𝜌=0,if𝑥^∗_𝑖∇>0,_𝑖𝑓𝑥^∗−𝜌=0,if𝑥^∗_𝑖<0,−𝜌≤∇_𝑖𝑓𝑥^∗≤𝜌,if𝑥^∗_𝑖=0. (2.4) Formally, we call the above conditions the optimality conditions for problem (1.1).

For any given 𝜏>0, we define a mapping 𝐻^𝜏=(𝐻^𝜏₁,𝐻^𝜏₂,…,𝐻^𝜏_𝑛)^𝑇∶ℝ^𝑛→ℝ^𝑛 by 𝐻^𝜏_𝑖(𝑥)≜max𝜏∇_𝑖𝑓(𝑥)−𝜌,min𝑥_𝑖,𝜏∇_𝑖𝑓(𝑥)+𝜌. (2.5) Then 𝐻^𝜏 is a continuous mapping and is closely related to problem (1.1). It is generally not differentiable in the sense of Fréchet derivative but semismooth in the sense of Qi and Sun [24]. The following proposition shows that the ℓ₁-norm minimization problem (1.1) is equivalent to a nonsmooth equation. It can be easily obtained by the use of the optimality conditions and the convexity of the function 𝑓 defined by (1.1).

Proposition 2.1. Let 𝜏>0 be any given constant. A point 𝑥^∗∈ℝ^𝑛 is a solution of problem (1.1) if and only if it satisfies 𝐻^𝜏𝑥^∗=0. (2.6)

The above proposition has reformulated problem (1.1) as a system of nonsmooth equations. Compared with the nonsmooth equation reformulation in [17], the dimension of (2.6) is only half of the dimension of the equation in [17].

Given 𝑎,𝑏,𝑐,𝑑∈ℝ. It is easy to verify that (see, e.g. [25]) {}{}=()(),{}{}=()()min𝑎,𝑏−min𝑐,𝑑1−𝑠)(𝑎−𝑐+𝑠𝑏−𝑑max𝑎,𝑏−max𝑐,𝑑1−𝑡)(𝑎−𝑐+𝑡𝑏−𝑑 (2.7) with {}{}𝑠=0,𝑎≤𝑏,𝑐≤𝑑;1,𝑎>𝑏,𝑐>𝑑;min𝑎,𝑏−min𝑐,𝑑+𝑐−𝑎,𝑏−𝑑+𝑐−𝑎otherwise,{}{}𝑡=0,𝑎≥𝑏,𝑐≥𝑑;1,𝑎<𝑏,𝑐<𝑑;max𝑎,𝑏−max𝑐,𝑑+𝑐−𝑎,𝑏−𝑑+𝑐−𝑎otherwise. (2.8) It is clear that 0≤𝑠, 𝑡≤1. By (2.5), we have for any 𝑥,𝑦∈ℝ^𝑛, it holds that 𝐻^𝜏_𝑖(𝑥)−𝐻^𝜏_𝑖𝑦=max𝜏∇_𝑖𝑓(𝑥)−𝜌,min𝑥_𝑖,𝜏∇_𝑖𝑓(𝑥)𝜏∇+𝜌−max_𝑖𝑓𝑦−𝜌,min𝑦_𝑖,𝜏∇_𝑖𝑓𝑦+𝜌=𝜏1−𝑡_𝑖∇_𝑖𝑓(𝑥)−∇_𝑖𝑓𝑦+𝑡_𝑖min𝑥_𝑖,𝜏∇_𝑖𝑓(𝑥)+𝜌−min𝑦_𝑖,𝜏∇_𝑖𝑓𝑦+𝜌=𝜏1−𝑡_𝑖∇_𝑖𝑓(𝑥)−∇_𝑖𝑓𝑦+𝑡_𝑖1−𝑠_𝑖𝑥_𝑖−𝑦_𝑖+𝜏𝑠_𝑖∇_𝑖𝑓(𝑥)−∇_𝑖𝑓𝑦=𝑡_𝑖1−𝑠_𝑖𝑥_𝑖−𝑦_𝑖+𝜏1−𝑡_𝑖+𝑡_𝑖𝑠_𝑖∇_𝑖𝑓(𝑥)−∇_𝑖𝑓𝑦, (2.9) where 0≤𝑠_𝑖, 𝑡_𝑖≤1. Define two diagonal matrixes 𝑆 and 𝑇 by 𝑆=diag𝑠₁,𝑠₂,…,𝑠_𝑛,𝑇=diag𝑡₁,𝑡₂,…,𝑡_𝑛. (2.10) Then we obtain 𝐻^𝜏(𝑥)−𝐻^𝜏𝑦=𝑇(𝐼−𝑆)𝑥−𝑦+𝜏(𝐼−𝑇+𝑇𝑆)∇𝑓(𝑥)−∇𝑓𝑦. (2.11) Since ∇𝑓(𝑥)=𝐴^𝑇(𝐴𝑥−𝑏), we get 𝐻^𝜏(𝑥)−𝐻^𝜏𝑦=𝑇(𝐼−𝑆)+𝜏(𝐼−𝑇+𝑇𝑆)𝐴^𝑇𝐴𝑥−𝑦. (2.12) The next proposition shows the Lipschitz continuity of 𝐻^𝜏 defined by (2.5).

Proposition 2.2. For each 𝜏>0, there exists a positive constant 𝐿(𝜏) such that 𝐻^𝜏(𝑥)−𝐻^𝜏𝑦≤𝐿(𝜏)𝑥−𝑦,∀𝑥,𝑦∈ℝ^𝑛. (2.13)

Proof. By (2.10) and (2.12), we have 𝐻^𝜏(𝑥)−𝐻^𝜏𝑦≤𝑇(𝐼−𝑆)+𝜏(𝐼−𝑇+𝑇𝑆)𝐴^𝑇𝐴𝑥−𝑦≤1+𝜏𝐴^𝑇𝐴𝑥−𝑦. (2.14) Let 𝐿(𝜏)≜1+𝜏‖𝐴^𝑇𝐴‖. Then (2.13) holds. The proof is complete.

The following proposition shows another good property of the system of nonsmooth equations (2.6).

Proposition 2.3. There exists a constant 𝜏^∗>0 such that for any 0<𝜏≤𝜏^∗, the mapping 𝐻^𝜏∶ℝ^𝑛→ℝ^𝑛 is monotone, that is 𝐻^𝜏(𝑥)−𝐻^𝜏𝑦,𝑥−𝑦≥0,∀𝑥,𝑦∈ℝ^𝑛. (2.15)

Proof. Let 𝐷_𝑖𝑖 be the 𝑖th diagonal element of 𝐴^𝑇𝐴. It is clear that 𝐷_𝑖𝑖>0, 𝑖=1,…,𝑛. Set 𝜏^∗≜min_𝑖{1/𝐷_𝑖𝑖}. Note that 𝐴^𝑇𝐴 is symmetric and positive semidefinite. Consequently, for any 𝜏∈(0,𝜏^∗], matrix 𝑇(𝐼−𝑆)+𝜏(𝐼−𝑇+𝑇𝑆)𝐴^𝑇𝐴 is also positive semidefinite. Therefore, it follows from (2.12) that 𝐻^𝜏(𝑥)−𝐻^𝜏𝑦,𝑥−𝑦≥0. (2.16) This completes the proof.

3. Algorithms and Their Convergence

In this section, we describe the proposed algorithms in detail and establish their convergence. Let 𝜏>0 be given. For simplicity, we omit 𝜏 and abbreviate 𝐻^𝜏(⋅) as 𝐻(⋅).

Algorithm 3.1 . (spectral gradient projection method (abbreviated as SGP)). Given initial point 𝑥₀∈ℝ^𝑛 and constants 𝜃₀=1, 𝑟>0, 𝜈≥0, 𝜎>0, 𝛾∈(0,1). Set 𝑘∶=0.

Step 1. Compute 𝑑_𝑘 by 𝑑_𝑘=−𝜃_𝑘𝐻𝑥_𝑘, (3.1) where for each 𝑘≥1, 𝜃_𝑘 is defined by 𝜃_𝑘=𝑠^𝑇_𝑘−1𝑠_𝑘−1𝑦^𝑇_𝑘−1𝑠_𝑘−1 (3.2) with 𝑠_𝑘−1=𝑥_𝑘−𝑥_𝑘−1 and 𝑦_𝑘−1=𝐻(𝑥_𝑘)−𝐻(𝑥_𝑘−1)+𝑟‖𝐻(𝑥_𝑘)‖^𝜈𝑠_𝑘−1. Stop if 𝑑_𝑘=0.

Step 2. Determine steplength 𝛼_𝑘=𝛾^𝑚𝑘 with 𝑚_𝑘 being the smallest nonnegative integer 𝑚 such that −𝐻𝑥_𝑘+𝛾^𝑚𝑑_𝑘,𝑑_𝑘≥𝜎𝛾^𝑚𝐻𝑥_𝑘+𝛾^𝑚𝑑_𝑘𝑑_𝑘. (3.3) Set 𝑧_𝑘∶=𝑥_𝑘+𝛼_𝑘𝑑_𝑘. Stop if ‖𝐻(𝑧_𝑘)‖=0.

Step 3. Compute 𝑥_𝑘+1=𝑥_𝑘−𝐻𝑧_𝑘,𝑥_𝑘−𝑧_𝑘𝐻𝑧_𝑘²𝐻𝑧_𝑘. (3.4) Set 𝑘∶=𝑘+1, and go to Step 1.

Remark 3.2. (i) The idea of the above algorithm comes from [18]. The major difference between Algorithm 3.1 and the method in [18] lies in the definition of 𝑦_𝑘−1. The choice of 𝑦_𝑘−1 in Step 1 follows from the modified BFGS method [26]. The purpose of the term 𝑟‖𝐻(𝑥_𝑘)‖^𝜈𝑠_𝑘−1 is to make 𝑦_𝑘−1 be closer to 𝐻(𝑥_𝑘)−𝐻(𝑥_𝑘−1) as 𝑥_𝑘 tends to a solution of (2.6).
(ii) Step 3 is called the projection step. It is originated in [20]. The advantage of the projection step is to make 𝑥_𝑘+1 closer to the solution set of (2.6) than 𝑥_𝑘. We refer to [20] for details.
(iii) Since −⟨𝐻(𝑥_𝑘),𝑑_𝑘⟩=‖𝐻(𝑥_𝑘)‖‖𝑑_𝑘‖, by the continuity of 𝐻, it is easy to see that inequality (3.3) holds for all 𝑚 sufficiently large. Therefore Step 2 is well defined and so is Algorithm 3.1.

The following lemma comes from [20].

Lemma 3.3. Let 𝐻∶ℝ^𝑛→ℝ^𝑛 be monotone and 𝑥,𝑦∈ℝ^𝑛 satisfy ⟨𝐻(𝑦),𝑥−𝑦⟩>0. Let 𝑥⁺𝐻𝑦=𝑥−,𝑥−𝑦𝐻𝑦²𝐻𝑦. (3.5) Then for any 𝑥^∗∈ℝⁿ satisfying H(𝑥^∗)=0, it holds that 𝑥⁺−𝑥^∗2≤𝑥−𝑥^∗2−𝑥⁺−𝑥². (3.6)

The following theorem establishes the global convergence for Algorithm 3.1.

Theorem 3.4. Let {𝑥_𝑘} be generated by Algorithm 3.1 and 𝑥^∗ a solution of (2.6). Then one has 𝑥_𝑘+1−𝑥^∗2≤𝑥_𝑘−𝑥^∗2−𝑥_𝑘+1−𝑥_𝑘². (3.7) In particular, {𝑥_𝑘} is bounded. Furthermore, it holds that either {𝑥_𝑘} is finite and the last iterate is a solution of the system of nonsmooth equations (2.6), or the sequence is infinite and lim_k→∞‖𝑥_𝑘+1−𝑥_𝑘‖=0. Moreover, {𝑥_𝑘} converges to some solution of (2.6).

Proof. The proof is similar to that in [18]. We omit it here.

Remark 3.5. The computational complexity of each of SGP’s steps is clear. In large-scale problems, most of the work is matrix-vector multiplication involving 𝐴 and 𝐴^𝑇. Steps 1 and 2 of SGP require matrix-vector multiplication involving 𝐴 or 𝐴^𝑇 two times each, while each iteration in GPSR-BB involves matrix-vector multiplication only two times. This may bring in more computational complexity. Therefore, we give a modification of SGP. The modified algorithm, which will be called MSGP in the rest of the paper, coincides with SGP except at Step 3, whose description is given below.

Algorithm 3.6 (modified spectral gradient projection Method (abbreviated as MSGP)). Given initial point 𝑥₀∈ℝ^𝑛 and constants 𝜃₀=1, 𝑟>0, 𝜎>0, 𝛾∈(0,1) a positive integer 𝑀. Set 𝑘∶=0.

Step 3. Let 𝑚=𝑘/𝑀. If 𝑚 is a positive integer, compute 𝑥_𝑘+1=𝑥_𝑘−𝐻𝑧_𝑘,𝑥_𝑘−𝑧_𝑘𝐻𝑧_𝑘²𝐻𝑧_𝑘; (3.8) otherwise, let 𝑥_𝑘+1=𝑧_𝑘. Set 𝑘∶=𝑘+1, and go to Step 1.

Lemma 3.7. Assume that {𝑥_𝑘} is a sequence generated by Algorithm 3.6 and 𝑥^∗∈ℝ^𝑛 satisfies 𝐻(𝑥^∗)=0. Let 𝜆_max(𝐴^𝑇𝐴) be the maximum eigenvalue of 𝐴^𝑇𝐴 and 𝜏∈(0,1/𝜆_max(𝐴^𝑇𝐴)]. Then it holds that 𝑥_𝑘+1−𝑥^∗≤𝑥_𝑘−𝑥^∗,𝑘=0,1,2,…. (3.9)

Proof. Let 𝑥_𝑘+1 be generated by (3.8). It follows from Lemma 3.3 that (3.9) holds. In the following, we assume that 𝑥_𝑘+1=𝑧_𝑘. Then, we obtain 𝑥_𝑘+1−𝑥^∗=𝑥_𝑘+𝛼_𝑘𝑑_𝑘−𝑥^∗=𝑥_𝑘−𝛼_𝑘𝜃_𝑘𝐻𝑥_𝑘−𝑥^∗+𝛼_𝑘𝜃_𝑘𝐻𝑥^∗=𝑥_𝑘−𝑥^∗−𝛼_𝑘𝜃_𝑘𝐻𝑥_𝑘−𝐻𝑥^∗. (3.10) This together with (2.12) implies that 𝑥_𝑘+1−𝑥^∗=𝑥_𝑘−𝑥^∗−𝛼_𝑘𝜃_𝑘𝑇_𝑘𝐼−𝑆_𝑘+𝜏𝐼−𝑇_𝑘+𝑇_𝑘𝑆_𝑘𝐴^𝑇𝐴𝑥_𝑘−𝑥^∗=1−𝛼_𝑘𝜃_𝑘𝑥_𝑘−𝑥^∗+𝛼_𝑘𝜃_𝑘𝐼−𝑇_𝑘+𝑇_𝑘𝑆_𝑘𝐼−𝜏𝐴^𝑇𝐴𝑥_𝑘−𝑥^∗≤1−𝛼_𝑘𝜃_𝑘𝑥_𝑘−𝑥^∗+𝛼_𝑘𝜃_𝑘𝐼−𝜏𝐴^𝑇𝐴𝑥_𝑘−𝑥^∗. (3.11) Let 𝜏∈(0,1/𝜆_max(𝐴^𝑇𝐴)]. Then we get 𝑥_𝑘+1−𝑥^∗≤𝑥_𝑘−𝑥^∗. (3.12) This completes the proof.

Now we establish a global convergence theorem for Algorithm 3.6.

Theorem 3.8. Let 𝜆_max(𝐴^𝑇𝐴) be the maximum eigenvalue of 𝐴^𝑇𝐴 and 𝜏∈(0,1/𝜆_max(𝐴^𝑇𝐴)]. Assume that {𝑥_𝑘} is generated by Algorithm 3.6 and 𝑥^∗ is a solution of (2.6). Then one has 𝑥_𝑘+1−𝑥^∗≤𝑥_𝑘−𝑥^∗,𝑘=0,1,2,…. (3.13) In particular, {𝑥_𝑘} is bounded. Furthermore, it holds that either {𝑥_𝑘} is finite and the last iterate is a solution of the system of nonsmooth equations (2.6), or the sequence is infinite and lim_𝑘→∞‖𝑥_𝑘+1−𝑥_𝑘‖=0. Moreover, {𝑥_𝑘} converges to some solution of (2.6).

Proof. We first note that if the algorithm terminates at some iteration 𝑘, then we have 𝑑_𝑘=0 or ‖𝐻(𝑧_𝑘)‖=0. By the definition of 𝜃_𝑘, we have 𝐻(𝑥_𝑘)=0 if 𝑑_𝑘=0. This shows that either 𝑥_𝑘 or 𝑧_𝑘 is a solution of (2.6).
Suppose that 𝑑_𝑘≠0 and ‖𝐻(𝑧_𝑘)‖≠0 for all 𝑘. Then an infinite sequence {𝑥_𝑘} is generated. It follows from (3.3) that 𝐻𝑧_𝑘,𝑥_𝑘−𝑧_𝑘=−𝛼_𝑘𝐻𝑧_𝑘,𝑑_𝑘≥𝜎𝛼²_𝑘𝐻𝑧_𝑘𝑑_𝑘>0. (3.14) Let 𝑥^∗ be an arbitrary solution of (2.6). By Lemmas 3.7 and 3.3, we obtain 𝑥_𝑘+1−𝑥^∗≤𝑥_𝑘−𝑥^∗,𝑥_𝑚𝑀+1−𝑥^∗2≤𝑥_𝑚𝑀−𝑥^∗2−𝑥_𝑚𝑀+1−𝑥_𝑚𝑀², (3.15) where 𝑚 is a nonnegative integer. In particular, the sequence {‖𝑥_𝑘−𝑥^∗‖} is nonincreasing and hence convergent. Moreover, the sequence {𝑥_𝑘} is bounded, and lim_𝑚→∞𝑥_𝑚𝑀+1−𝑥_𝑚𝑀=0. (3.16) Following from (3.8) and (3.14), we have 𝑥_𝑚𝑀+1−𝑥_𝑚𝑀=𝐻𝑧_𝑚𝑀,𝑥_𝑚𝑀−𝑧_𝑚𝑀𝐻𝑧_𝑚𝑀≥𝜎𝛼²_𝑚𝑀𝑑_𝑚𝑀. (3.17) This together with (3.16) yields lim_𝑚→∞𝛼_𝑚𝑀𝑑_𝑚𝑀=0. (3.18) Now we consider the following two possible cases:(i)liminf_𝑚→∞‖𝐻(𝑥_𝑚𝑀)‖=0;(ii)liminf_𝑚→∞‖𝐻(𝑥_𝑚𝑀)‖=𝜖>0.
If (i) holds, then by the continuity of 𝐻 and the boundedness of {𝑥_𝑚𝑀}, it is clear that the sequence {𝑥_𝑚𝑀} has some accumulation point 𝑥^∗ such that 𝐻(𝑥^∗)=0. Since the sequence {‖𝑥_𝑘−𝑥^∗‖} converges, it must hold that {𝑥_𝑘} converges to 𝑥^∗.
If (ii) holds, then by the boundedness of {𝑥_𝑚𝑀} and the continuity of 𝐻, there exist a positive constant 𝐶 and a positive integer 𝑚₀ such that 12𝜖≤𝐻𝑥_𝑚𝑀≤𝐶,∀𝑚≥𝑚₀. (3.19) On the other hand, from (3.2) and the definitions of 𝑠_𝑘−1 and 𝑦_𝑘−1, we have 𝜃_𝑚𝑀=𝑠^𝑇_𝑚𝑀−1𝑠_𝑚𝑀−1𝑦^𝑇_𝑚𝑀−1𝑠_𝑚𝑀−1=𝑠^𝑇_𝑚𝑀−1𝑠_𝑚𝑀−1𝐻𝑥_𝑚𝑀−𝐻𝑥_𝑚𝑀−1,𝑥_𝑚𝑀−𝑥_𝑚𝑀−1+𝑟𝐻(𝑥_𝑚𝑀)^𝜈𝑠^𝑇_𝑚𝑀−1𝑠_𝑚𝑀−1, (3.20) which together with (3.19) and Propositions 2.2 and 2.3 implies 1𝐿+𝑟𝐶^𝜈≤𝜃_𝑚𝑀≤2^𝜈𝑟𝜖^𝜈,∀𝑚≥𝑚₀. (3.21) Consequently, we obtain from (3.1), (3.19), and (3.21) 𝑑_𝑚𝑀=𝜃_𝑚𝑀𝐻𝑥_𝑚𝑀≥𝜖2(𝐿+𝑟𝐶^𝜈),𝑑_𝑚𝑀≤2^𝜈𝐶𝑟𝜀^𝜈,∀𝑚≥𝑚₀. (3.22) Therefore, it follows from (3.18) that lim_𝑚→∞𝛼_𝑚𝑀=0. By the line search rule, we have for all 𝑚 sufficiently large, 𝑚_𝑘−1 will not satisfy (3.3). This means −𝐻𝑥_𝑚𝑀+𝛾^𝑚𝑘−1𝑑_𝑚𝑀,𝑑_𝑚𝑀<𝜎𝛾^𝑚𝑘−1𝐻𝑥_𝑚𝑀+𝛾^𝑚𝑘−1𝑑_𝑚𝑀𝑑_𝑚𝑀. (3.23) Since {𝑥_𝑚𝑀} and {𝑑_𝑚𝑀} are bounded, we can choose subsequences of {𝑥_𝑚𝑀} and {𝑑_𝑚𝑀} converging to 𝑥^∗∗ and 𝑑^∗∗, respectively. Taking the limit in (3.23) for the subsequence, we obtain −𝐻𝑥^∗∗,𝑑^∗∗≤0. (3.24) However, it is not difficult to deduce from (3.1), (3.19), and (3.21) that −𝐻𝑥^∗∗,𝑑^∗∗>0. (3.25) This yields a contradiction. Consequently, liminf_𝑚→∞‖𝐻(𝑥_𝑚𝑀)‖=𝜖>0 is not possible. The proof is then complete.

4. Applications to Compressed Sensing and Image Restoration

In this section, we apply the proposed algorithms, that is, SGP and MSGP, to solve some practical problems arising from the compressed sensing and image restoration. We will compare the proposed algorithms with SGCS, SPARSA, and GPSR-BB. The system of nonsmooth equations in SGCS is 𝐹(𝑧)≜min{𝑧,𝜏(𝐻𝑧+𝑐)}=0, (4.1) where 𝑧, 𝑐, 𝐻 are defined as those in [17]. The test problems are associated with applications in the areas of compressed sensing and image restoration. All experiments were carried out on a Lenovo PC (2.53 GHz, 2.00 GB of RAM) using Matlab 7.8. The parameters in SGCS are specified as follows: 𝜏=7,𝛽=1,𝛾=1.2,𝜉=10⁻⁴,𝜌=0.1. (4.2) The parameters in SGP and MSGP are specified as follows: 𝜏=7,𝜎=1,𝑟=0.8,𝛾=0.5,𝑀=10. (4.3) Throughout the experiments, we choose the initial iterate to be 𝑥⁰=0.

In our first experiment, we consider a typical CS scenario, where the goal is to reconstruct a length-𝑛 sparse signal (in the canonical basis) from 𝑚 observations, where 𝑚<𝑛. The 𝑚×𝑛 matrix 𝐴 is obtained by first filling it with independent samples of the standard Gaussian distribution and then orthonormalizing the rows. Due to the storage limitations of PC, we test a small size signal with 𝑚=1024, 𝑛=4096. The observed vector is 𝑏=𝐴𝑥_orig+𝜉, where 𝜉 is Gaussian white noise with variance 𝜎²=10⁻⁴ and 𝑥_orig is the original signal with 50 randomly placed ±1 spikes and with zeros in the remaining elements. The regularization parameter is chosen as 𝜌=0.05‖𝐴^𝑇𝑏‖_∞. We compare the performance of SGP and MSGP with that of SGCS, SPARSA, and GPSR-BB by solving the problem and choose 𝜈=1 in SGP and MSGP algorithms. We measure the quality of restoration by means of mean squared error (MSE) to the original signal 𝑥_orig defined by MSE=1𝑛𝑥−𝑥_orig², (4.4) where 𝑥 is the restored signal. To perform this comparison, we first run the SGCS algorithm and stop the algorithm if the following inequality is satisfied: 𝑥^𝑘+1−𝑥^𝑘𝑥^𝑘<10⁻⁵, (4.5) and then run each of the other algorithms until each reaches the same value of the objective function reached by SGCS.

The original signal and the estimation obtained by solving (1.1) using the MSGP method are shown in Figure 1. We can see from Figure 1 that MSGP does an excellent job at locating the spikes with respect to the original signal. In Figure 2, we plot the evolution of the objective function versus iteration number and CPU time, for these algorithms. It is readily to see that MSGP worked faster than other algorithms.

(a)

(b)

(c)

Figure 1 

From top to bottom: original signal, observation, and reconstruction obtained by MSGP.

(a)

(b)

Figure 2 

The objective function plotted against iteration number and CPU time for SGCS, SPARSA, GPSR-BB, SGP, and MSGP.

In the second experiment, we test MSGP for three image restoration problems based on the images as House, Cameraman, and Barbara. House and Cameraman images are of size 256×256 and the other is of size 512×512. All the pixels are contaminated by Gaussian noise with the standard deviation of 0.05 with blurring. The blurring function is chosen to be a two-dimensional Gaussian, ℎ=1𝑖,𝑗1+𝑖²+𝑗², (4.6) truncated such that the function has a support of 9×9. The image restoration problem has the form (1.1), where 𝜌=0.0005 and 𝐴=HW are the composition of the 9×9 uniform blur matrix and the Haar discrete wavelet transform (DWT) operator. We compare the performance of MSGP with that of SGCS, SPARSA, and GPSR-BB by solving the problem and choose 𝜈=0 in the MSGP method. As usual, we measure the quality of restoration by signal-to-noise ratio (SNR) defined as SNR=10×log₁₀𝑥_orig²𝑥_orig−𝑥², (4.7) where 𝑥_orig and 𝑥 are the original and restored images, respectively. We first run SGCS and stop the process if the following inequality is satisfied: 𝑥^𝑘+1−𝑥^𝑘𝑥^𝑘<10⁻⁵, (4.8) and then run the other algorithms until their objective function value reach SGCS's value. Table 1 reports the number of iterations (Iter), the CPU time in seconds (Time), and the SNR to the original images (SNR).

Table 1 

Test results for SGCS, SPARSA, GPSR-BB, SGP, and MSGP in image restoration.

It is easy to see from Table 1 that the MSGP is competitive with the well-known algorithms: SPARSA and GPSR-BB, in computing time and number of iterations and improves the SGCS greatly. Therefore we conclude that the MSGP provides a valid approach for solving ℓ₁-norm minimization problems arising from image restoration problems.

Preliminary numerical experiments show that SGP and MSGP algorithms have improved SGCS algorithm greatly. This may be because the system of nonsmooth equations solved here has lower dimension than that in [17] and the modification to projection steps that we made reduces the computational complexity.

Acknowledgments

The authors would like to thank the anonymous referee for the valuable suggestions and comments. L. Wu was supported by the NNSF of China under Grant 11071087; Z. Sun was supported by the NNSF of China under Grants 11126147 and 11201197.