arXiv:math/0503066v2 [math.NA] 7 Dec 2005
Stable Signal Recovery from
Incomplete and Inaccurate Measurements
Emmanuel Candes
†
, Justin Romberg
†
, and Terence Tao
♯
† Applied and Computational Mathematics, Caltech, Pasadena, CA 91125
♯ Department of Mathematics, University of California, Los Angeles, CA 90095
February, 2005; Revised June 2005
Abstract
Suppose we wish to recover a vector x
0
∈ R
m
(e.g. a digital signal or image) from
incomplete and contaminated observations y = Ax
0
+ e; A is a n by m matrix with far
fewer rows than columns (n ≪ m) and e is an error term. Is it possible to recover x
0
accurately based on the data y?
To recover x
0
, we consider the solution x
♯
to the ℓ
1
-regularization problem
min kxk
ℓ
1
subject to kAx − yk
ℓ
2
≤ ǫ,
where ǫ is the size of the error term e. We show that if A obeys a uniform uncertainty
principle (with unit-normed columns) and if the vector x
0
is sufficiently sparse, then
the solution is within the noise level
kx
♯
− x
0
k
ℓ
2
≤ C · ǫ.
As a first example, suppose that A is a Gaussian random matrix, then stable recovery
occurs for a lmost all such A’s provided that the number of nonzeros of x
0
is of about the
same order as the number of o bservations. As a s e c ond instance, suppose o ne observes
few Fourier samples of x
0
, then stable recovery occurs for almost a ny set of n coefficients
provided that the number of nonzeros is of the order of n/[log m]
6
.
In the case where the error term vanishes, the r e c overy is of course exact, and this
work actually provides novel insights on the exact recovery phenomenon discussed in
earlier pa pers. The methodology also explains why one can also very nearly recover
approximately sparse signals.
Keywords. ℓ
1
-minimization, basis pursuit, restricted orthonormality, sparsity, singular
values of random matrices.
Acknowledgments. E. C. is partially supported by a National Science Foundation grant
DMS 01-40698 (FRG) and by an Alfred P. Sloan Fellowship. J. R. is supported by National
Science Foundation grants DMS 01-40698 and IT R ACI-0204932. T. T. is supported in part
by grants from the Packard Foundation.
1
1 Introduction
1.1 Exact recovery of sparse signals
Recent papers [2–5,10] have developed a series of powerful results about the exact recovery
of a finite signal x
0
∈ R
m
from a very limited number of observations. As a representative
result from this literature, consider the problem of recovering an unknown sparse signal
x
0
(t) ∈ R
m
; that is, a signal x
0
whose support T
0
= {t : x
0
(t) 6= 0} is assumed to have
small card inality. All we know about x
0
are n linear measurements of the form
y
k
= hx
0
, a
k
i k = 1, . . . , n or y = Ax
0
,
where the a
k
∈ R
m
are known test signals. Of special interest is the vastly underdetermined
case, n ≪ m, wh ere th ere are many more unknowns than observations. At first glance, this
may seem impossible. However, it turns out that one can actually recover x
0
exactly by
solving the convex p rogram
1
(P
1
) min kxk
ℓ
1
subject to Ax = y, (1)
provided that the matrix A ∈ R
n×m
obeys a uniform uncertainty principle.
The uniform uncertainty principle, introduced in [5] and refined in [4], essentially states that
the n × m measurement matrix A obeys a “restricted isometry hypothesis.” To introduce
this notion, let A
T
, T ⊂ {1, . . . , m} be the n × |T | submatrix obtained by extracting the
columns of A corresponding to the indices in T . Then [4] defines th e S-restricted isometry
constant δ
S
of A which is the smallest quantity such that
(1 −δ
S
) kck
2
ℓ
2
≤ kA
T
ck
2
ℓ
2
≤ (1 + δ
S
) kck
2
ℓ
2
(2)
for all subsets T with |T | ≤ S and coefficient sequences (c
j
)
j∈T
. This property essentially
requires that every set of columns with cardinality less than S approximately behaves like
an orthonormal system. It was shown (also in [4]) that if S verifies
δ
S
+ δ
2S
+ δ
3S
< 1, (3)
then solving (P
1
) recovers any sparse signal x
0
with support size obeying |T
0
| ≤ S.
1.2 Stable recovery from imperfect measurements
This paper d evelops results for the “imperfect” (and far m ore realistic) scenarios where the
measurements are noisy and the signal is not exactly sparse. Everyone would agree that
in most practical situations, we cannot assume that Ax
0
is known with arbitrary precision.
More appr op riately, we will assume instead that one is given “noisy” data y = Ax
0
+e, where
e is some unknown perturbation bounded by a known amount kek
ℓ
2
≤ ǫ. To be broadly
applicable, our recovery procedure must be stable: small changes in the observations should
result in small changes in the recovery. This wish, however, may be quite hopeless. How can
we possibly hope to recover our signal when not only the available information is severely
incomplete but in addition, the few available observations are also inaccurate?
1
(P
1
) can even be recast as a linear program [6].
2
Consider nevertheless (as in [12] for example) the convex program searching, among all
signals consistent with the data y, for that with minimum ℓ
1
-norm
(P
2
) min kxk
ℓ
1
subject to kAx − yk
ℓ
2
≤ ǫ. (4)
The first result of this paper shows that contrary to the belief expressed above, the solution
to (P
2
) recovers an unknown sparse object with an error at most proportional to the noise
level. Our condition for stable recovery again involves the restricted isometry constants.
Theorem 1 Let S be such that δ
3S
+ 3δ
4S
< 2. Then for any signal x
0
supported on T
0
with |T
0
| ≤ S and any perturbation e with kek
ℓ
2
≤ ǫ, the solution x
♯
to (P
2
) obe ys
kx
♯
− x
0
k
ℓ
2
≤ C
S
· ǫ, (5)
where the constant C
S
may only depend on δ
4S
. For reasonable values of δ
4S
, C
S
is well
behaved; e.g. C
S
≈ 8.82 for δ
4S
= 1/5 and C
S
≈ 10.47 for δ
4S
= 1/4.
It is interesting to note that for S obeying the condition of the theorem, the reconstruction
from noiseless data is exact. It is quite possible that for some matrices A, this condition
tolerates larger values of S than (3).
We would like to offer two comments. First, the matrix A is rectangular with many more
columns than rows. As such, most of its singular values are zero. As emphasized earlier,
the fact that the severely ill-posed matrix inversion keeps the perturbation from “blowing
up” is rather remarkable and perhaps unexpected.
Second, no recovery method can perform fundamentally better for arbitrary perturbations
of s ize ǫ. To see why this is true, suppose one had available an oracle letting us know, in
advance, the support T
0
of x
0
. With this additional information, the problem is well-posed
and one could reconstruct x
0
by the method of Least-Squares for examp le,
ˆx =
(
(A
∗
T
0
A
T
0
)
−1
A
∗
T
0
y on T
0
0 elsewhere.
In the abs en ce of any other information, one could easily argue that no method would
exhibit a fundamentally better performance. Now of course, ˆx −x
0
= 0 on th e complement
of T
0
while on T
0
ˆx − x
0
= (A
∗
T
0
A
T
0
)
−1
A
∗
T
0
e,
and since by hypothesis, the eigenvalues of A
∗
T
0
A
T
0
are well-behaved
2
kˆx − x
0
k
ℓ
2
≈ kA
∗
T
0
ek
ℓ
2
≈ ǫ,
at least for pertur bations concentrated in the row s pace of A
T
0
. In short, obtaining a
reconstruction with an error term whose size is guaranteed to be proportional to th e noise
level is the best one can hope for.
Remarkably, not only can we recover sparse input vectors but one can also stably recover
approximately sparse vectors, as we have the following companion theorem.
2
Observe the role played by the singular values of A
T
0
in the analysis of the oracle error.
3
Theorem 2 Suppose that x
0
is an arbitrary vec tor in R
m
and let x
0,S
be the truncated
vector corre sponding to the S largest values of x
0
(in absolute value). U nder the hypothesis
of Theorem 1, the solution x
♯
to (P
2
) obe ys
kx
♯
− x
0
k
ℓ
2
≤ C
1,S
· ǫ + C
2,S
·
kx
0
− x
0,S
k
ℓ
1
√
S
. (6)
For reasonable values of δ
4S
the constants in (5) are well behaved; e.g. C
1,S
≈ 12.04 and
C
1,S
≈ 8.77 for δ
4S
= 1/5.
Roughly speaking, the theorem says that minimizing ℓ
1
stably recovers the S-largest entries
of an m-dimensional unknown vector x from n measurements only.
We now specialize this result to a commonly discussed mo del in mathematical signal pro-
cessing, namely, the class of compressible signals. We say that x
0
is compressible if its
entries obey a power law
|x
0
|
(k)
≤ C
r
· k
−r
, (7)
where |x
0
|
(k)
is the kth largest value of x
0
(|x
0
|
(1)
≥ |x
0
|
(2)
≥ . . . ≥ |x
0
|
(m)
), r > 1, and
C
r
is a constant which depends only on r. Such a model is appropriate for the wavelet
co efficients of a piecewise smooth signal, f or example. If x
0
obeys (7), then
kx
0
− x
0,S
k
ℓ
1
√
S
≤ C
′
r
· S
−r+1/2
.
Observe now that in this case
kx
0
− x
0,S
k
ℓ
2
≤ C
′′
r
· S
−r+1/2
,
and for generic elements obeying (7), there are no fundamentally better estimates available.
Hence, we see that with n measurements only, we achieve an approximation error which is
almost as good as that one would obtain by knowing everything about the signal x
0
and
selecting its S-largest entries.
As a last remark, we would like to point out that in the noiseless case, Theorem 2 improves
upon an earlier result from Cand`es and Tao, see also [8]; it is sharper in the sense that 1)
this is a deterministic statement and there is no probability of failure, 2) it is universal in
that it holds for all signals, 3) it gives upper estimates with better bounds and constants,
and 4) it holds for a wider range of values of S.
1.3 Examples
It is of course of interest to know which matrices obey the uniform u ncertainty principle
with good isometry constants. Using tools from random matrix theory, [3,5,10] give several
examples of matrices such that (3) holds for S on the order of n to within log factors.
Examples include (proofs and additional discussion can be foun d in [5]):
• Random matrices with i.i.d. entries. Suppose the entries of A are i.i.d. Gaussian with
mean zero and variance 1/n, then [5, 10, 17] show that the condition for Theorem 1
holds with overwhelming pr obab ility when
S ≤ C · n/ log(m/n).
4
In fact, [4] gives numerical values for the constant C as a function of the ratio n/m.
The same conclusion applies to binary matrices with independent entries taking values
±1/
√
n with equal probability.
• Fourier ensemble. Suppose now that A is obtained by selecting n rows from the
m × m discrete Fourier transform and r en ormalizing the columns so that they are
unit-normed. If the rows are selected at random, the condition for Theorem 1 holds
with overwhelming probability f or S ≤ C · n/(log m)
6
[5]. (For simplicity, we have
assumed that A takes on real-valued entries although our theory clearly accommodates
complex-valued matrices so that our discussion holds for both complex and real-valued
Fourier transforms.)
This case is of special interest as reconstructing a digital signal or image from incom-
plete Fourier data is an important inverse problem with applications in biomedical
imaging (MRI and tomography), Astrophysics (interferometric imaging), and geo-
physical exploration.
• General orthogonal measurement ensembles. Suppose A is obtained by selecting n
rows from an m by m orthonormal matrix U and renormalizing the columns so that
they are unit-normed. Then [5] shows that if the rows are selected at random, the
condition for Theorem 1 holds with overwhelming probability provided
S ≤ C ·
1
µ
2
·
n
(log m)
6
, (8)
where µ :=
√
m max
i,j
|U
i,j
|. Observe that for the Four ier matrix, µ = 1, and thus
(8) is an extension of the Fourier ensemble.
This fact is of significant p ractical relevance because in many situations, signals of
interest may not be sparse in the time domain but rather may be (approximately)
decomposed as a sparse superposition of waveforms in a fixed orthonormal basis Ψ;
e.g. in a nice wavelet basis. S uppose that we use as test signals a set of n vectors taken
from a second orthonormal basis Φ. We then solve (P
1
) in the coefficient domain
(P
′
1
) min kαk
ℓ
1
subject to Aα = y,
where A is obtained by extracting n rows from the orthonormal matrix U = ΦΨ
∗
. The
recovery condition then depends on the mutual coherence µ between the measurement
basis Φ and the sparsity basis Ψ which measures the similarity between Φ and Ψ;
µ(Φ, Ψ) =
√
m max |hφ
k
, ψ
j
i|, φ
k
∈ Φ, ψ
j
∈ Ψ.
1.4 Prior work and innovations
The problem of recovering a sparse vector by minimizing ℓ
1
under linear equality constraints
has recently received much attention, mostly in the context of Basis Pursuit, where the goal
is to uncover sparse signal decompositions in overcomplete dictionaries. We refer the reader
to [11, 13] and the references therein for a full discussion.
We would especially like to note two works by Donoho, Elad, and Temlyakov [12], and Tropp
[18] that also study the recovery of sparse signals from noisy observations by solving (P
2
)
(and other closely related optimization programs), and give conditions for stable recovery.
In [12], the sparsity constraint on the underlying signal x
0
depends on the magnitude of
5
