# Nonmonotone Adaptive Barzilai-Borwein Gradient Algorithm for Compressed Sensing

Authors: Yuanying Qiu; Jianlei Yan; Fanyong Xu
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

We study a nonmonotone adaptive Barzilai-Borwein gradient algorithm forl1-norm minimization problems arising from compressed sensing. At each iteration, the generated search direction enjoys descent property and can be easily derived by minimizing a local approximal quadratic model and simultaneously taking the favorable structure of the l1-norm. Under some suitable conditions, its global convergence result could be established. Numerical results illustrate that the proposed method is promising and competitive with the existing algorithms NBBL1 and TwIST.


---

## Body

## 1. Introduction

In recent years, algorithms for finding sparse solutions to underdetermined linear systems of equations have been intensively investigated in signal processing and compressed sensing. The fundamental principle of compressed sensing (CS) is that a sparse signalx-∈Rn can be recovered from the underdetermined linear system b=Ax-, where A∈Rm×n (often m≪n). By defining l0 norm (∥x∥0) of a vector as the number of nonzero components in x, one natural way to reconstruct x- from the system is to solve the following problem:
(1)min⁡x∈Rn∥x∥0s.t.b=Ax
via certain reconstruction technique. However, the l0 norm problem is combinatorial and generally computationally intractable. A fundamental decoding model in CS is to replace l0 norm by l1 norm, which is defined as ∥x∥1=∑i=1n‍|x(i)|. The resulting adaptation of (1) is the so-called basis pursuit (BP) problem [1]
(2)min⁡x∈Rn∥x∥1s.t.b=Ax.
It is shown that, under some reasonable conditions, problem (2) can produce the desired solutions with high probability [2]. When b contains some noise in most practical applications, the constraint in (2) should be relaxed to the penalized least squares problem:
(3)min⁡x∈Rn12∥b-Ax∥22+μ∥x∥1.
Here, μ>0 is related to the Lagrange multiplier of the constraint in (2).It follows from some existing results that if a signal is sparse or approximately sparse in some orthogonal basis, then an accurate recovery can be obtained whenA is a random matrix projection [3]. Quite a number of algorithms have been proposed and studied for solving the aforementioned l1-problems arising in CS. Recently, some first-order methods are popular for solving (3), such as the projection steepest descent method [4] and the gradient projection algorithm (GPSR) proposed by Figueiredo et al. [5]. Moreover, based on a smoothing technique studied in [6], a fast and accurate first-order algorithm called NESTA was proposed in [7], and so on. By an operator splitting technique, Hale et al. derive the iterative shrinkage/thresholding fixed-point continuation algorithm (FPC) [8]. One most widely studied first-order method is the iterative shrinkage/thresholding (IST) method [9–11], which is designed for wavelet-based image deconvolution. TwIST [12, 13] and FISTA [14] speed up the performance of IST and have virtually the same complexity but with better convergence properties. Another closely related method is the sparse reconstruction algorithm SpaRSA [15], which is to minimize nonsmooth convex problem with separable structures. In [16], the authors proposed nonsmooth equations-based method for l1-norm problems. SPGL1 [17] solves a least squares problem with l1-norm constraint by the spectral gradient projection method with an efficient Euclidean projection on l1-norm ball. In [18], Yun and Toh studied a block coordinate gradient descent (CGD) method for solving (3). Recently, the alternating direction method (ADM) has received much attention for solving total variation regularization problems for image restoration and is also capable of solving the l1-norm regularization problems in CS [19, 20].Very recently, Xiao et al. propose a Barzilai-Borwein gradient algorithm [21] for solving l1 regularized nonsmooth minimization problems (NBBL1) [22], in which they approximate f locally by a convex quadratic model at each iteration, where the Hessian is replaced by the multiples of a spectral coefficient with an identity matrix. Motivated by them, we propose a nonmonotone adaptive Barzilai-Borwein gradient algorithm for l1-norm minimization in compressed sensing, which is based on a new quasi-Newton equation [23] and a new adaptive spectral coefficient. Under reasonable assumptions, its convergence result could be established. Numerical experiments illustrate that the proposed method is efficient to recover a sparse signal arising in compressive sensing and outperforms NBBL1.A full description of the proposed algorithm is presented in the next section. Meanwhile, we establish its global convergence under some suitable conditions. In Section3, some numerical results were reported to illustrate the efficiency of the proposed method. Finally, we have a conclusion section.


## 2. Proposed Algorithm and Convergence Result

In this section, we construct an iterative algorithm to solve thel1-norm regularization problems arising from the spare solution recovery in compressed sensing. Before stating the steps of our method, we first give a brief description of preliminary results for the following unconstrained optimization:
(4)min⁡x∈ℝnf(x),
where f:ℝn→ℝ is a continuously differentiable function.In [23], Wei et al. proposed a new quasi-Newton equation and then derived a new conjugacy condition by using this new quasi-Newton equation. Using the Taylor formula for the objective function f(x),
(5)f(x)≈fk+gkT(x-xk)+12(x-xk)TGk(x-xk),
where fk (resp., Gk) denotes the function value (resp., Hessian matrix) and gk denotes ∇f(xk) at xk. Hence, substituting x=xk-1,
(6)fk-1≈fk-gkTsk-1+12sk-1TGksk-1.
Therefore,
(7)sk-1TGksk-1≈2(fk-1-fk)+2gkTsk-1=2(fk-1-fk)+(gk-1+gk)Tsk-1+sk-1Tyk-1.
Consider Bk as a new approximation of Gk such that
(8)sk-1TBksk-1=sk-1Tyk-1+ϑk-1,
where sk-1=xk-xk-1, yk-1=gk-gk-1, y-k-1=yk-1+(ϑk-1/∥sk-1∥2)sk-1, and ϑk-1=2(fk-1-fk)+(gk+gk-1)Tsk-1. This suggests the following new quasi-Newton equation:
(9)Bksk-1=y¯k-1.
In [24], Li et al. make a modification of the y¯k in (9) as follows:
(10)yk-1*=yk-1+max⁡{ϑk-1,0}∥sk-1∥2sk-1,
and Yuan and Wei [25] make some further studies on it. Observe that this new quasi-Newton equation contains not only gradient value information but also function value information at the present and the previous step. In general, such Bk+1 will be produced by updating Bk with some typical and popular formulae such as BFGS, DFP, and SR1. Furthermore, let the approximation Hessian Bk be a diagonal matrix with positive components; that is, Bk=λkI with an identity matrix I and λk>0. Then, the quasi-Newton condition (9) possesses the following form:
(11)λkIsk-1=yk-1*.
Multiplying both sides by sk-1T, it follows that
(12)λkBB1*=sk-1Tyk-1*∥sk-1∥22,
and, multiplying both sides by yk-1T, it gives
(13)λkBB2*=∥yk-1*∥22sk-1Tyk-1*,
where yk-1* is defined as (10). If sk-1Tyk-1*>0 holds, then the matrix λkI is positive definite which ensures that the search direction -λk-1gk is descent at current point.Now, we focus our attention on thel1-norm minimization problem (3). The algorithm can be described as the iterative form
(14)xk+1=xk+αkdk,
where αk is a stepsize and dk is a search direction defined by minimizing a quadratic approximated model of F(x)=(1/2)∥Ax-b∥22+μ∥x∥1. Since l1-term is not differentiable, hence, at current xk, objective function F is approximated by the quadratic approximation Qk,
(15)F(xk+d)=f(xk+d)+μ∥xk+d∥1≈f(xk)+∇f(xk)Td+λk2∥d∥22+μ[∥xk∥1+∥xk+hd∥1-∥xk∥1h]Δ__Qk(d),
where f(x)=∥Ax-b∥22 and h is a small positive number. The term in [·] can be considered as an approximate Taylor expansion of ∥xk+d∥1 with a small h, and the case h=1 reduces the equivalent form ∥xk+d∥1. Minimizing (15) yields
(16)min⁡d∈RnQk(d)⟺min⁡d∈Rn∇f(xk)Td+λk2∥d∥22+μ∥xk+hd∥1h⟺min⁡d∈Rnh2λk(∇f(xk)Td+λk2∥d∥22+μ∥xk+hd∥1h)⟺min⁡d∈Rn12∥xk+hd-(xk-hλk∇f(xk))∥22+μhλk∥xk+hd∥1⟺min⁡d∈Rn∑i=1n‍{12(xki+hdi-(xki-hλk∇fi(xk)))2dddddddddd+μhλk∥xki+hdi∥1},
where xki, di, and ∇fi(xk) denote the ith component of xk, d, and ∇f(xk), respectively. The favorable structure of (16) admits the explicit solution
(17)xki+hdki=max⁡{|xki-hλk∇fi(xki)|-μhλk,0}×xki-(h/λk)∇fi(xk)|xki-(h/λk)∇fi(xk)|.
Hence, the search direction at current point is
(18)dk=-1h[xki-hλk∇fi(xk)|xki-hλk∇fi(xk)|xk-max⁡{|xki-hλk∇fi(xki)|-μhλk,0}sssssssaa×xki-(h/λk)∇fi(xk)|xki-(h/λk)∇fi(xk)|[xki-hλk∇fi(xk)|xki-hλk∇fi(xk)|xk-max⁡{|xki-hλk∇fi(xki)|-μhλk,0}].In this paper, we adopt the following adaptive Barzilai-Borwein step in (18):
(19)λkABB={λkBB1*,ifλkBB2*/λkBB1*>0.9λkBB2*,otherwise,
where λkBB1* and λkBB2* are defined in (12) and (13), respectively.In the light of all derivations above, we now describe the nonmonotone adaptive Barzilai-Borwein gradient algorithm (abbreviated as NABBL1) (see Algorithm1).Algorithm 1:NABBL1 for (3).
Initialization: Give initial point x0∈Rn, set parameters μ>0, α¯>0 and ρ∈(0,1),
δ
∈
(
0,1
), h∈(0,1] and positive integer m¯. Set k=0.
Step  1. If ∥dk∥2=0, then stop. Otherwise, continue.
Step  2. Compute dk via (18).
Step  3. Compute αk via following nonmonotone line search
F(xk+αkdk)≤max⁡0≤j≤m(k)F(xk-j)+δαk(∇f(xk)Tdk+μ∥xk+hdk∥1-μ∥xk∥1h),
where the smallest nonnegative integerjk such as the stepsize αk=α¯ρjk satisfing above,
m
(
0
)
=
0 and 0≤m(k)≤min⁡{m(k-1)+1,m¯}.
Step  4. Let xk+1=xk+αkdk.
Step  5. Let k=k+1. Go to Step  1.Remark 1.
We have shown that ifλk>0, then the generated direction is descent. However, in this case, the condition λk>0 may fail to be fulfilled and the hereditary descent property is not guaranteed any more. To cope with this defect, we should keep the sequence {λk} uniformly bounded; that is, for sufficiently small λmin⁡>0 and sufficiently large λmax⁡>0, the λk is forced as
(20)λk=min⁡{λmax⁡,max⁡{λk,λmin⁡}}.
This approach ensures that λk is bounded from zero and subsequently ensures that dk is descent at per-iteration.
We prepare to show our main global convergence result of algorithm NABBL1. The desirable convergence is directly from Theorem 3.3 of [22]; we state it for completeness here.Theorem 2.
Let the sequences{xk} and {dk} be generated by Algorithm 1. Then, there exists a subsequence K such that
(21)lim⁡k→∞,k∈K∥dk∥2=0.


## 3. Experimental Results

In this section, we describe some experiments to illustrate the good performance of the algorithm NABBL1 for reconstructing sparse signals. These experiments are all tested in Matlab R2010a. The relative error is used to measure the quality of the reconstructive signals which is defined as(22)RelErr=∥x^-x-∥2∥x-∥2,
where x^ denotes the reconstructive signal and x- denotes the original signal.In our experiments, we consider a typical compressive sensing scenario, where the goal is to reconstruct ann length sparse signal from m observations. The random A is the Gaussian matrix whose elements are generated from shape i.i.d. normal distributions 𝒩(0,1) (randn (m,n) in Matlab). In real applications, the measurement b is usually contaminated by noise; that is, b=Ax+η, where η is the Gaussian noise distributed as 𝒩(0,σ2I).We test a small size signal withn=211, m=29; the original contains randomly k=26 nonzero elements. The proposed algorithm starts at a zero point and terminates when the relative change of two successive points is sufficiently small; that is,
(23)∥xk-xk-1∥2∥xk-1∥2<τ.
In this experiment, we take τ=10-4, h=0.8, λmin⁡=10-30, and λmax⁡=1030. In the line search, we choose α¯=10-2, ρ=0.35, δ=10-4, and m¯=5. The original signal, the limited measurement, and the reconstructed signal when the noise level σ2=10-3 are given in Figure 1.(a) Original signal with length of 2048 and 64 nonzero elements. (b) Noisy measurement with length of 512. (c) Recovered signal by NABBL1 (red circle) versus original signal (blue peaks).
(a)
(b)
(c)Comparing (a) to (c) in Figure1, we clearly see that the original sparse signal is restored almost exactly. We see that all the blue peaks are circled by the red circles, which illustrates that the original signal has been found almost exactly. Altogether, this simple experiment shows that our algorithm performs quite well and provides an efficient approach to recover large sparse nonnegative signal.We choose four different signals with noise level ofσ2=10-3 compared with algorithms NBBL1 [22] and TwIST [13] in our next experiment. In order to test the speed of the algorithms more fairly, we list the average of the five results in Table 1. Numerical results are listed in Table 1, in which we report the CPU time in seconds (time) required for the whole reconstructing process and the relative error (RelErr). From Table 1, we can see that algorithm NABBL1 is faster than algorithms NBBL1 and TwIST, and the number of iterations of algorithm NABBL1 is less than that of the algorithms NBBL1 and TwIST with different signals.Table 1
Test results of NABBL1 and NBBL1 with different combinations of(n,m,k).
n
m
k
NABBL1
NBBL1
TwIST
Time
RelErr
Time
RelErr
Time
RelErr
8192
2048
256
2.9844
1.73
e
-
02
3.5313
1.97
e
-
02
4.0156
4.83
e
-
02
4096
1024
128
0.6875
1.93
e
-
02
0.9063
2.09
e
-
02
1.7500
4.96
e
-
02
2048
512
64
0.2344
1.66
e
-
02
0.3125
1.90
e
-
02
0.5625
4.50
e
-
02
1024
256
32
0.0625
1.55
e
-
02
0.0938
2.23
e
-
02
0.1718
4.54
e
-
02From Figure2, NABBL1 usually decreases relative errors faster than NBBL1 and TwIST throughout the entire iteration process. We conclude that NABBL1 provides an efficient approach for solving l1 regularized nonsmooth problem from compressed sensing and is competitive with or performs better than NBBL1 and TwIST.Comparison result of NABBL1, NBBL1, and TwIST withσ2=10-3. (a) Original signal with length of 2048 and 64 nonzero elements. (b) Original signal with length of 4096 and 128 nonzero elements. The x-axes represent the CPU time in seconds. The y-axes represent the relative error.
(a)
(b)


## 4. Conclusion

In this paper, we proposed a nonmonotone adaptive Barzilai-Borwein algorithm (NABBL1) for solving al1 regularized least squares problem arising from spare solution recovery in compressed sensing. At each iteration, the generated search direction enjoys descent property and can be easily derived by minimizing a local approximal quadratic model and simultaneously taking the favorable structure of the l1-norm. Numerical results illustrate that the proposed method is promising and competitive with the existing algorithms NBBL1 and the two-step IST (TwIST). Our future topic is to extend NABBL1 method for solving matrix trace norm minimization problems in compressed sensing or some minimization problems in computed tomography reconstruction [26, 27].



---

*Source: 410104-2014-04-07.xml*