This paper introduces a new, fast and accurate algorithm

for solving problems in the area of compressed sensing,

and more generally, in the area of signal and image reconstruction

from indirect measurements. This algorithm

is inspired by recent progress in the development of novel

first-order methods in convex optimization, most notably

Nesterov’s smoothing technique. In particular, there is a

crucial property thatmakes thesemethods extremely efficient

for solving compressed sensing problems. Numerical

experiments show the promising performance of our

method to solve problems which involve the recovery of

signals spanning a large dynamic range

Bobin, J.

Candès, E. J.

English

Caltech Authors - Main

A FAST AND ACCURATE FIRST-ORDER ALGORITHM
FOR COMPRESSED SENSING
J. Bobin and E. J. Cande`s
Applied and Computational Mathematics (MC 217-50), California Institute of Technology
Pasadena, CA 91125, USA
ABSTRACT
This paper introduces a new, fast and accurate algorithm
for solving problems in the area of compressed sensing,
and more generally, in the area of signal and image re-
construction from indirect measurements. This algorithm
is inspired by recent progress in the development of novel
first-order methods in convex optimization, most notably
Nesterov’s smoothing technique. In particular, there is a
crucial property that makes these methods extremely effi-
cient for solving compressed sensing problems. Numeri-
cal experiments show the promising performance of our
method to solve problems which involve the recovery of
signals spanning a large dynamic range.
Index Terms— Compressed sensing, sparsity, 1 and
total-variation minimization, smoothing technique in op-
timization.
1. INTRODUCTION
Compressed sensing (CS) [3, 2, 4] is a new sampling
theory based on the revelation that one can exploit spar-
sity or compressibility when acquiring signals of gene-
ral interest, and that one can design nonadaptive sam-
pling techniques that condense the information in a com-
pressible signal into a small amount of data. There are
early indications showing that this revelationmay change
the way engineers think about signal acquisition in areas
ranging from analog-to-digital conversion, digital optics,
magnetic resonance imaging, and seismics.
In compressed sensing, one acquires a signal x ∈ Rn
by collecting data of the form b = Ax+n : x is the signal
of interest (more precisely, its coefficient sequence in a
representation where it is assumed to be fairly sparse), A
is an m × n “sampling” matrix, and n is a noise term.
A standard approach in CS attempts to reconstruct x by
solving
min
x
f(x) s. t. ‖b−Ax‖2 < , (1)
where 2 is an estimated upper bound on the noise po-
wer. The choice of the regularizing function f depends
on prior assumptions about the signal x of interest : if
x is (approximately) sparse, an appropriate convex func-
tion is the 1 norm (as advocated by the CS theory), if
x is a piecewise constant object, the total-variation (TV)
norm provides accurate recovery results, and so on.
Solving large-scale problems such as (1) (think of
the unknown as a mega-pixel image) is challenging. Al-
though, one cannot review the vast literature on this sub-
ject, the majority of the algorithms that have been pro-
posed are unable to solve these problems accuratelywith
low computational complexity. On the one hand, interior-
point methods are accurate but problematic because they
need to solve large systems of linear equations along the
way (the Newton step). On the other hand, newly intro-
duced techniques based upon iterative thresholding, see
[5, 1, 6] and the many earlier references therein converge
slowly in the sense that they require a very large num-
ber of iterations when high accuracy is required. This is
important as in many applications of interest, the signals
exhibit a wide dynamic range.
In this paper, we build upon the recent work of Nes-
terov [8], which develops a series of first-order methods
with improved convergence rates, by adapting these ideas
to develop algorithms for solving signal recovery pro-
blems. These novel algorithms show a lot of promise ;
they are fast, accurate and seem robust in the sense that
their performance does not depend on the fine tuning of
various controlling parameters.
2. FAST FIRST-ORDER ALGORITHMS
2.1. General formulation
Consider the saddle point problem
min
x∈Qp
max
u∈Qd
〈u,Wx〉, (2)
where x ∈ Rn, u ∈ Rp and W ∈ Rp×n. We will re-
fer to Qp and Qd as the primal and dual feasible sets.
Put f(x) = maxu∈Qd〈u,Wx〉 ; the function f is convex
but generally nonsmooth. If W = I and Qd = {u :
‖u‖∞ ≤ 1}, then f(x) = ‖x‖1 and (2) is the CS re-
covery problem if we set Qp = {x : ‖b − Ax‖2 < }.
1457978-1-4244-5654-3/09/$26.00 ©2009 IEEE ICIP 2009
TV minimization can also be cast as (2). In [8], Nesterov
proposed to substitute f by the smooth approximation
fμ(x) = max
u∈Qd
〈u,Wx〉 − μpd(u), (3)
where pd(u) is a prox-function for Qd ; that is, pd(u) is
continuous and strongly convex on Qd, with convexity
parameter σd (we will assume that pd vanishes on Qd).
Nesterov proved that fμ is continuously differentiable,
that ∇fμ(x) = W ∗uμ(x) where uμ(x) is the optimal
solution of (3), and that ∇fμ is Lipschitz with constant
Lμ =
1
μσd
‖W‖2 (‖W‖ is the operator norm ofW ). Nes-
terov’s algorithm minimizes fμ over Qp by iteratively
estimating three sequences {xk}, {yk} and {zk} while
smoothing the feasible set Qp. The algorithm depends
on two scalar sequences {αk} and {τk} discussed below,
and takes the following form :
Initialize x0. For k ≥ 0,
1. Compute ∇fμ(xk).
2. Compute yk :
yk =Argminx∈Qp
Lμ
2
‖x − xk‖
2
2
+ 〈∇fμ(xk), x− xk〉.
3. Compute zk :
zk =Argminx∈Qp
Lμ
2σp
pp(x)+
Pk
i=0 αi〈∇fμ(xi), x−xi〉.
4. Update xk :
xk =τkzk + (1 − τk)yk .
Stop when a given criterion is valid.
In this algorithm the function pp(x) is a prox-function
for the primal feasible set Qp with strong convexity pa-
rameter σp. At step k, yk is the current guess at the op-
timal solution. If we only performed the second step of
the algorithm with yk−1 instead of xk, we would obtain
a standard first-order technique. The novelty is that the
sequence zk “keeps in mind” the previous iterations, and
the point xk at which the gradient of fμ is evaluated hap-
pens to be a subtle average between zk and yk. It has been
shown in [8] that if αk = 1/2(k+1) and τk = 2/(k+3),
then the algorithm converges to x = Argminx∈Qp fμ(x)
with the convergence rate
fμ(yk)− fμ(x) ≤ 4Lμpp(x
)
(k + 1)2σp
. (4)
The decay is far better than what is achievable via
standard gradient-based optimization techniques (k−2
vs. k−1).
2.2. The choice of the primal set prox-function
A good primal prox-function is a smooth function
that is likely to have some positive effect near the so-
lution. In the setting of (1), a suitable smoothing prox-
function may be
pp(x) =
1
2
||b −Ax||22 +
ρ
2
‖x− x0‖22 (5)
for some x0 ∈ Rn, e.g. an initial guess of the solution.
Notice that the bound on the error at iteration k in (4) is
proportional to pp(x) ; choosing x0 wisely (a good first
guess) can make pp(x) small ; when nothing is known
about the solution as in our later experiments, a natural
choice may be x0 = A∗b.
2.3. Projecting onto Qp
Whenm < n, the feasible set is unbounded which is
a departure fromNesterov’s algorithm. Now applying the
algorithm above requires computing the projection onto
Qp to estimate yk and zk. After an appropriate change of
variables, the minimization problems in steps 2 and 3 can
be recast as
min
x∈Qp
γ
2
‖b−Ax‖22 +
1
2
‖x‖22 + 〈c, x〉. (6)
The constraint x ∈ Qp can be relaxed as follows with
an appropriate choice of λ : minx γ+λ2 ‖b − Ax‖22 +
1
2
‖x‖22 + 〈c, x〉. The solution is then given by
x = (I + βA∗A)−1(βA∗b − c), (7)
where β = γ + λ. Hence, we need to solve a positive
definite system of the form (I + βA∗A)Δx = Δb. Note
that the eigenvalues of this system are of the form 1 +
βλ(A∗A), where λ(A∗A) are the eigenvalues of A∗A.
Random matrices. When A ∈ Rm×n, m < n,
is a random matrix with i.i.d. entries with mean zero
and variance 1, it is well known that the nonzero ei-
genvalues of A∗A/n are all very near the interval
[(1 −
√
m/n)2, (1 +
√
m/n)2] [7]. As a consequence,
apart from many eigenvalues at 1, the eigenvalues of
I + βA∗A are highly clustered. In this case, choosing
β = γ + max {0, ‖A∗A(b −Ac)‖2/− γ − 1} is close
to the optimal value and, hence, steps 2 and 3 can be
performed with only a few Conjugate Gradients (CG)
iterations. To drive this point home, Table 1 reports the
results of a simple experiment in which A is a random
Gaussian matrix with n = 256 and a varying value ofm.
The value of β is set to 0.01/‖A∗A‖, the vectors b and c
are random, and CG iterations are applied to compute an
approximate solution xCG to (7). Each xCG is compared
to the true solution x obtained by solving the system with
a direct solver. The results demonstrate that for random
measurement matrices, steps 2 and 3 can be computed
with a handful of CG iterations with excellent accuracy.
Projections. In a wide range of CS applications, the
matrix A is a projection onto a subspace of Rn as when
the Fourier, the Hadamard or the noiselet transforms are
1458
m/n # iterations ‖(xCG − x)/|x|‖∞
0.01 3 2 10−4
0.04 4 5 10−5
0.08 4 7 10−5
0.16 4 4 10−4
Table 1. Conjugate gradients and projection onto Qp.
Averaged numbers of CG steps and averaged maximum
(entry-wise) relative errors. These averages are over 100
trials for each value ofm/n.
subsampled. Steps 2 and 3 become trivial since one has
an explicit formula
x =
(
I − β
1 + β
A∗A
)
(βA∗b− c).
where β = γ + max {0, ‖A∗A(b−Ac)‖2/− γ − 1}.
In practice, the computational cost of this step is just the
cost of applyingA∗A ; e.g. two FFTs when A is a partial
DFT.
2.4. Nesterov and Continuation
In the spirit of the fixed continuation technique intro-
duced in [6], the nature of the algorithm makes it ame-
nable to some sort of continuation. We discuss two pos-
sibilities.
In Step 2 and 3 of the proposed algorithm, updating
yk and zk is similar to a projected gradient step. For ins-
tance, yk is of the form
yk = PQp
(
xk − 1
Lμ
∇fμ(xk)
)
, (8)
where PQp is the projector onto Qp. As Lμ ∝ μ−1,
the step size scales like μ. Thus the algorithm converges
faster when μ is larger. For 1 minimization, ∇fμ =
1/μW ∗P{x|‖x‖∞≤μ}(Wxk) . The parameter μ is thus
equivalent to a threshold, and the accuracy on each en-
try of the final estimate will also scale like μ. High accu-
racy requires a small value of μ. To overcome this trade-
off between speed and accuracy, we propose a kind of
continuation scheme. The idea is to apply Nesterov’s al-
gorithm sequentially with a decreasing sequence {μt}.
Initialize μ0 and x0 = xμ0 . For t ≥ 1,
1. Apply Nesterov’s algorithm with μ = μt and x0 = xμt−1
2. Decrease the value of μ : μt+1 = ρμt with ρ < 1
Stop when the desired value of μ is reached.
Above, x0 is the point defining the prox-function, see
(5), and xμt is the computed solution for μ = μt. As
μt decreases, this algorithm sequentially computes tigh-
ter approximations xμt to x. Numerical results are given
in Section 3.2.
3. NUMERICAL EXPERIMENTS
We apply the algorithm discussed above to the clas-
sical 1 norm and 2D total-variation norm minimization
problems. Both these examples assume incomplete data.
3.1. Recovery of signals with wide dynamic range
We assume that x is sparse and that A is a partial
Fourier transform (this models the problem of recove-
ring a signal with a sparse spectrum from data sampled
at a rate much lower than the Nyquist rate) and use 1
minimization to recover x from b. The number of en-
tries of x is n = 65, 536, the number of measurements
is m = n/4 and only n/100 	 655 entries of x are
nonzero. The amplitudes of the nonzero components are
selected at random ; the lowest value is equal to 1 and the
highest is equal to 10d where d = 1, · · · , 4. Note that
when d = 3, 4, the signals exhibit a large dynamic range
(DR), which is challenging for most numerical methods.
The stopping criterion is 1− fμ(xk+1)/fμ(xk) < 10−6.
As discussed earlier, the parameter μ is of paramount
importance as it fixes the accuracy of the smooth approxi-
mation. The lower this parameter, the closer xμ to x.
When minimizing the 1 norm, one can check that the
entries of xμ will differ from those of x by an amount
which is on the order of μ. In our experiment, the std. of
the additive Gaussian noise is σ = 0.01 and we set μ =
30σ as this gives a good balance between speed and ac-
curacy (the lowest nonzero entry of x has an amplitude
equal to 1). The 2 norm error on the residual has been
set to  = σ
√
m2 + 2
√
2m.
Figure 1 presents the results of an experiment when
the DR is 104 (40dB) and Table 2 reports on experiments
corresponding to various values of the DR. Each value in
Table 2 is an average computed from 25 random trials.
The results highlight that whatever the DR, the error
on the nonzero entries is at most on the order of μ. Hence,
the algorithm achieves very high accuracy whenever nee-
ded, as in the case of high DR. Since μ is smaller than
the smallest entry, the consequence is that all the spikes
(all the nonzero components) are correctly detected. As
expected, the higher the DR, the higher the number of
iterations. However, the number of iterations is still re-
markably small. For a DR of 40dB, the algorithm returns
a highly accurate solution in just about 1,200 FFTs.
3.2. TV Minimization
The total-variation norm of a 2D digital object xij is
given by
‖x‖TV =
∑
i,j
√
[D1x]2ij + [D2x]
2
ij ,
1459
0 2000 4000 6000 8000 10000 12000 14000 16000
100
101
102
103
104
Fig. 1. Solid line : original signal. Dots - ◦ : recovered
signal.
d # iterations ∞ error off Support(x) Detection rate
1 47 0.11 100%
2 64 0.15 100%
3 109 0.16 100%
4 305 0.16 100%
Table 2. 1-recovery results as a function of dynamic
range. Average number of iterations required to reach
convergence, average maximum size off the support of
x, and average detection rate in%.
where D1 and D2 are the horizontal and vertical diffe-
rences : e. g. [D1x]ij = xi+1,j −xi,j . SetW to be the li-
near map defined via [Wx]ij = ([D1x]ij , [D2x]ij). Then
minimizing the TV norm can be cast as a saddle point
problem since ‖x‖TV = maxu∈Qd〈Wx, u〉, where Qd
is the family u = {uij} with uij ∈ R2 and ‖ui,j‖2 ≤ 1
for each pair (i, j).
We apply the algorithm to solve a TV minimization
problem from incomplete frequency data. Here, we col-
lect noisy Fourier coefficients about the 256 × 256 clas-
sical Logan-Shepp phantom ; these m = 5357 samples
lie on radial lines just as in [3]. The level of the additive
noise is σ = 0.01 ;  =
√
m2 + 2
√
2mσ.
With μ = 10−7, it takes 1092 iterations to reach
convergence1. Applying a continuation technique with
a sequence of μ’s equal to {μ0, μ0/2, · · · , 10−7}, with
μ0 = ‖A∗b‖TV /n 	 10−4, lowers the iteration count ;
convergence is reached in 512 steps. The SNR is equal to
58.2 dB without continuation and equal to 66.4 dB with
continuation. These preliminary results show that conti-
nuation with Nesterov’s algorithm can definitely improve
the overall speed of convergence.
1. convergence is reached when 1− fμ(xk+1)/fμ(xk) < 10−6
Fig. 2. Top : Original 256 × 256 image. Bottom-left :
Recovery without continuation. Bottom-right : Recovery
with continuation.
4. REFERENCES
[1] J.-F. Cai, S. Osher, and Z. Shen. Linearized Bregman ite-
rations for compressed sensing. Math. Comp. (to appear),
2008.
[2] E. Cande`s and T. Tao. Near optimal signal recovery from
random projections : Universal encoding strategies ? IEEE
Trans. on Information Theory, 52(12) :5406–5425, 2006.
[3] E. J. Cande`s, J. Romberg, and T. Tao. Robust uncertainty
principles : Exact signal reconstruction from highly incom-
plete frequency information. IEEE Trans. on Information
Theory, 52(2) :489–509, 2006.
[4] D. L. Donoho. Compressed sensing. IEEE Trans. on Infor-
mation Theory, 52(4) :1289–1306, April 2006.
[5] M.A. Figueiredo, R. Nowak, and S.J.Wright. Gradient pro-
jection for sparse reconstruction : Application to compres-
sed sensing and other inverse problems. IEEE Journal of
Selected Topics in Signal Processing, 1(4) :586 – 597, Dec.
2007.
[6] E. T. Hale, W. Yin, and Y. Zhang. A fixed-point conti-
nuation method for l1-regularized minimization with ap-
plications to compressed sensing. Technical Report - Rice
University, 2007.
[7] Mehta M.L. Random Matrices. Academic Press, 3rd edi-
tion, 1991.
[8] Y. Nesterov. Smooth minimization of non-smooth func-
tions. Math. Program., Serie A, 103 :127–152, 2005.
1460


A fast and accurate first-order algorithm for compressed sensing

A ﬁxed-point continuation method for l1-regularized minimization with applications to compressed sensing.

C a n d `

Compressed sensing.

Gradient projection for sparse reconstruction : Application to compressed sensing and other inverse problems.

Linearized Bregman iterations for compressed sensing.

Near optimal signal recovery from random projections : Universal encoding strategies ?

Near optimal signal recovery from random projections : Universal encoding strategies?

Random Matrices.

Robust uncertainty principles : Exact signal reconstruction from highly incomplete frequency information.

Smooth minimization of non-smooth functions.

https://authors.library.caltech.edu/19577/1/Bobin2009p111492009_16Th_Ieee_International_Conference_On_Image_Processing_Vols_1-6.pdf

A fast and accurate first-order algorithm for compressed sensing

Abstract

Similar works

Full text

Available Versions

Caltech Authors - Main