Subspace Pursuit for Compressive Sensing Signal Reconstruction

0.43 MB·English

by  Wei Dai

1
Subspace Pursuit for Compressive Sensing Signal
Reconstruction
Wei Dai and Olgica Milenkovic
Department of Electrical and Computer Engineering
University of Illinois at Urbana-Champaign
9 Abstract—We propose a new method for reconstruction of strongpriorknowledgeofsignalsparsityallowsforrecovering
0 sparse signals with and without noisy perturbations, termed the x using m N projections only. One of the outstanding
0 subspace pursuit algorithm. The algorithm has two important results in CS≪theory is that the signal x can be reconstructed
2 characteristics: low computational complexity, comparable to using optimization strategies aimed at finding the sparsest
that of orthogonal matching pursuit techniques when applied

n to very sparse signals, and reconstruction accuracy of the same signal that matches with the m projections. In other words,
a orderasthatofLPoptimizationmethods.Thepresentedanalysis the reconstruction problem can be cast as an l minimization
0
J shows that in the noiseless setting, the proposed algorithm can problem [2]. It can be shown that to reconstruct a K-sparse

8 exactly reconstruct arbitrary sparse signals provided that the
signalx, l minimizationrequiresonly m=2K randompro-
sensing matrix satisfies the restricted isometry property with a 0
jectionswhenthe signalandthe measurementsare noise-free.
] constant parameter. In the noisy setting and in the case that
A the signal is not exactly sparse, it can be shown that the mean Unfortunately, the l0 optimization problem is NP-hard. This
N squarederrorofthereconstructionisupperboundedbyconstant issuehasledtoalargebodyofworkinCStheoryandpractice
multiplesof themeasurement and signal perturbation energies. centeredaroundthedesignofmeasurementandreconstruction
.
cs IndexTerms—Compressivesensing,orthogonal matchingpur- algorithms with tractable reconstruction complexity.
[ suit, reconstruction algorithms, restricted isometry property, The work by Donoho and Candès et. al. [1], [3], [4], [5]
sparse signal reconstruction. demonstrated that CS reconstruction is, indeed, a polynomial
3
time problem – albeit under the constraint that more than
v
1 I. INTRODUCTION 2K measurementsareused.Thekeyobservationbehindthese
1 Compressive sensing (CS) is a sampling method closely findings is that it is not necessary to resort to l0 optimization
8 to recover x from the under-determined inverse problem; a
connected to transform coding which has been widely used
0
much easier l optimization, based on Linear Programming
. in modern communication systems involving large scale data 1
3 (LP) techniques, yields an equivalent solution, as long as the
samples. A transform code converts input signals, embedded
0 sampling matrix Φ satisfies the so called restricted isometry
in a high dimensional space, into signals that lie in a space
8
property (RIP) with a constant parameter.
0 of significantly smaller dimensions. Examples of transform
While LP techniques play an important role in designing
: coders include the well known wavelet transforms and the
v computationally tractable CS decoders, their complexity is
ubiquitous Fourier transform.
i
X Compressive sensing techniques perform transform cod- still highly impractical for many applications. In such cases,
the need for faster decodingalgorithms- preferablyoperating
r ing successfully whenever applied to so-called compressible
a in linear time - is of critical importance, even if one has
and/orK-sparsesignals,i.e.,signalsthatcanberepresentedby
to increase the number of measurements. Several classes of
K N significant coefficients over an N-dimensionalbasis.
≪ low-complexity reconstruction techniques were recently put
Encoding of a K-sparse, discrete-time signal x of dimension
forward as alternatives to linear programming (LP) based
N isaccomplishedbycomputingameasurementvectorythat
recovery, which include group testing methods [6], and al-
consists of m N linear projections of the vector x. This
≪ gorithms based on belief propagation [7].
can be compactly described via
Recently, a family of iterative greedy algorithms received
y=Φx. significant attention due to their low complexity and simple
geometricinterpretation.They include the OrthogonalMatch-
Here, Φ represents an m N matrix, usually over the field
× ing Pursuit (OMP), the Regularized OMP (ROMP) and the
of real numbers. Within this framework, the projection basis
Stagewise OMP (StOMP) algorithms. The basic idea behind
isassumedtobeincoherentwith thebasisinwhichthesignal
these methods is to find the support of the unknown signal
has a sparse representation [1].
sequentially.Ateachiterationofthealgorithms,oneorseveral
Althoughthe reconstructionof the signal x RN from the
∈ coordinates of the vector x are selected for testing based
possiblynoisyrandomprojectionsisanill-posedproblem,the
on the correlation values between the columns of Φ and
This work is supported by NSF Grants CCF 0644427, 0729216 and the the regularized measurement vector. If deemed sufficiently
DARPAYoungFaculty Awardofthesecondauthor. reliable, the candidate columnindices are subsequentlyadded
WeiDaiandOlgicaMilenkovic arewiththeDepartmentofElectrical and
to the current estimate of the support set of x. The pursuit
ComputerEngineering, University ofIllinois atUrbana-Champaign, Urbana,
IL61801-2918USA(e-mail: weidai07@ uiuc.edu;[email protected]). algorithms iterate this procedure until all the coordinates in
2
the correct support set are included in the estimated support candidate word is identified. SP employs a search strategy in
set.ThecomputationalcomplexityofOMPstrategiesdepends which a constant number of vectors is expurgated from the
on the number of iterations needed for exact reconstruction: candidatelist. Thisfeatureis mainlyintroducedforsimplicity
standard OMP always runs through K iterations, and there- of analysis: one can easily extend the algorithm to include
fore its reconstruction complexity is roughly O(KmN) (see adaptiveexpurgationstrategiesthat do notnecessarily operate
Section IV-C for details). This complexity is significantly on fixed-sized lists.
smaller than that of LP methods, especially when the signal Incompressivesensing,themajorchallengeassociatedwith
sparsity level K is small. However, the pursuit algorithms do sparse signal reconstruction is to identify in which subspace,
not have provable reconstruction quality at the level of LP generated by not more than K columns of the matrix Φ,
methods. For OMP techniques to operate successfully, one the measured signal y lies. Once the correct subspace is
requires that the correlation between all pairs of columns determined, the non-zero signal coefficients are calculated by
of Φ is at most 1/2K [8], which by the Gershgorin Circle applying the pseudoinversion process. The defining character
Theorem [9] represents a more restrictive constraint than of the SP algorithm is the method used for finding the K
the RIP. The ROMP algorithm [10] can reconstruct all K- columns that span the correct subspace: SP tests subsets of
sparse signals provided that the RIP holds with parameter K columns in a group, for the purpose of refining at each
δ 0.06/√logK,whichstrengthenstheRIP requirements stage an initially chosen estimate for the subspace. More
2K
≤
for l -linear programming by a factor of √logK. specifically,thealgorithmmaintainsalistofK columnsofΦ,
1
The main contribution of this paper is a new algorithm, performs a simple test in the spanned space, and then refines
termed the subspace pursuit (SP) algorithm. It has provable thelist. Ify doesnotlieinthecurrentestimateforthecorrect
reconstruction capability comparable to that of LP methods, spanning space, one refines the estimate by retaining reliable
and exhibits the low reconstruction complexity of matching candidates, discarding the unreliable ones while adding the
pursuit techniques for very sparse signals. The algorithm can same number of new candidates. The “reliability property” is
operate both in the noiseless and noisy regime, allowing captured in terms of the order statistics of the inner products
for exact and approximate signal recovery, respectively. For ofthereceivedsignalwiththecolumnsofΦ,andthesubspace
any sampling matrix Φ satisfying the RIP with a constant projection coefficients.
parameter independent of K, the SP algorithm can recover As a consequence,the maindifferencebetweenROMP and
arbitrary K-sparse signals exactly from its noiseless mea- the SP reconstruction strategy is that the former algorithm
surements. When the measurements are inaccurate and/or the generates a list of candidates sequentially, without back-
signal is not exactly sparse, the reconstruction distortion is tracing: it starts with an empty list, identifies one or several
upper bounded by a constant multiple of the measurement reliable candidates during each iteration, and adds them to
and/or signal perturbation energy. For very sparse signals the already existing list. Once a coordinate is deemed to be
with K const √N, which, for example, arise in certain reliable and is added to the list, it is not removed from it
≤ ·
communicationscenarios,thecomputationalcomplexityofthe until the algorithm terminates. This search strategy is overly
SP algorithm is upper bounded by O(mNK), but can be restrictive, since candidates have to be selected with extreme
further reduced to O(mNlogK) when the nonzero entries caution. In contrast, the SP algorithm incorporates a simple
of the sparse signal decay slowly. method for re-evaluating the reliability of all candidates at
The basic idea behind the SP algorithm is borrowed from each iteration of the process.
coding theory, more precisely, the A∗ order-statistic algo- At the time of writing this manuscript, the authors became
rithm [11] for additive white Gaussian noise channels. In awareofthe relatedworkbyJ.Tropp,D.NeedellandR. Ver-
this decoding framework, one starts by selecting the set of shynin[12],describingasimilarreconstructionalgorithm.The
K most reliable information symbols. This highest reliability main difference between the SP algorithm and the CoSAMP
information set is subsequently hard-decision decoded, and algorithmof[12]isinthemannerinwhichnewcandidatesare
the metric of the parity checks corresponding to the given addedtothelist.Ineachiteration,intheSPalgorithm,onlyK
information set is evaluated. Based on the value of this newcandidatesareadded,while the CoSAMPalgorithmadds
metric,someofthelow-reliabilitysymbolsinthemostreliable 2K vectors. This makes the SP algorithm computationally
information set are changed in a sequential manner. The more efficient, but the underlying analysis more complex. In
algorithm can thereforebe seen as operatingon an adaptively addition, the restricted isometry constant for which the SP
modified coding tree. If the notion of “most reliable symbol” algorithm is guaranteed to converge is larger than the one
is replaced by “column of sensing matrix exhibiting highest presentedin [12]. Finally, this paper also contains an analysis
correlation with the vector y”, the notion of “parity-check of the number of iterations needed for reconstruction of a
metric” by “residual metric”, then the above method can be sparse signal (see Theorem 6 for details), for which there is
easily changed for use in CS reconstruction. Consequently, no counterpart in the CoSAMP study.
one can perform CS reconstruction by selecting a set of K The remainder of the paper is organized as follows. Sec-
columns of the sensing matrix with highest correlation that tion II introduces relevant concepts and terminology for de-
spanacandidatesubspaceforthesensedvector.Ifthedistance scribingtheproposedCSreconstructiontechnique.SectionIII
of the received vector to this space is deemed large, the containsthealgorithmicdescriptionoftheSPalgorithm,along
algorithm incrementally removes and adds new basis vectors with a simulation-based study of its performance when com-
according to their reliability values, until a sufficiently close paredwithOMP,ROMP,andLPmethods.SectionIVcontains
3
themainresultofthepaperpertainingtothenoiselesssetting: for K m, 0 δ 1, if for all index sets I 1, ,N
≤ ≤ ≤ ⊂{ ··· }
a formal proof for the guaranteed reconstruction performance such that I K and for all q R|I|, one has
| |≤ ∈
and the reconstruction complexity of the SP algorithm. Sec-
(1 δ) q 2 Φ q 2 (1+δ) q 2.
tion V contains the main result of the paper pertaining to the − k k2 ≤k I k2 ≤ k k2
noisy setting. Concluding remarks are given in Section VI, We define δ , the RIP constant, as the infimum of all
K
while proofs of most of the theorems are presented in the parameters δ for which the RIP holds, i.e.
Appendix of the paper.
δ :=inf δ : (1 δ) q 2 Φ q 2 (1+δ) q 2,
K − k k2 ≤k I k2 ≤ k k2
II. PRELIMINARIES n I K, q R|I| .
A. Compressive SensingandtheRestricted Isometry Property ∀| |≤ ∀ ∈
o
Let supp(x) denote the set of indices of the non-zero Remark 1 (RIP and eigenvalues): If a sampling matrix
coordinates of an arbitrary vector x = (x1,...,xN), and let Φ Rm×N satisfies the RIP with parameters (K,δK), then
∈
supp(x) = denotethesupportsizeofx,orequivalently, for all I 1, ,N such that I K, it holds that
0
|itsl nor|m1.kA·kssumenextthatx RN is anunknownsignal ⊂{ ··· } | |≤
0 ∈ 1 δ λ (Φ∗Φ ) λ (Φ∗Φ ) 1+δ ,
with supp(x) K, and let y Rm be an observation of x − K ≤ min I I ≤ max I I ≤ K
| | ≤ ∈
via M linear measurements, i.e., where λ (Φ∗Φ ) and λ (Φ∗Φ ) denote the minimal
min I I max I I
and maximal eigenvalues of Φ∗Φ , respectively.
y=Φx, I I
Remark 2 (Matrices satisfying the RIP): Mostknownfam-
where Φ Rm×N is henceforth referred to as the sampling
∈ ilies of matrices satisfying the RIP property with optimal or
matrix.
near-optimal performance guarantees are random. Examples
We are concerned with the problem of low-complexity
include:
recovery of the unknown signal x from the measurement y.
1) Random matrices with i.i.d. entries that follow either
A naturalformulationof the recoveryproblemis within an l
0
the Gaussian distribution, Bernoulli distribution with
norm minimization framework, which seeks a solution to the
zero mean and variance 1/n, or any other distribution
problem
that satisfies certain tail decay laws. It was shown in
min x subjecttoy=Φx.
k k0 [13] that the RIP for a randomly chosen matrix from
Unfortunately,the abovel0 minimizationproblemisNP-hard, such ensembles holds with overwhelming probability
and hence cannot be used for practical applications [3], [4]. whenever
m
Onewaytoavoidusingthiscomputationallyintractablefor- K C ,
mulationis to considera l -regularizedoptimizationproblem, ≤ log(N/m)
1
where C is a function of the RIP constant.
min x subjecttoy=Φx,
k k1 2) Random matrices from the Fourier ensemble. Here,
where one selects m rows from the N N discrete Fourier
N transform matrix uniformly at ran×dom. Upon selection,
x = x
k k1 | i| the columns of the matrix are scaled to unit norm. The
Xi=1 resulting matrix satisfies the RIP with overwhelming
denotes the l1 norm of the vector x. probability, provided that
Themainadvantageofthel minimizationapproachisthat
1 m
it is a convex optimization problem that can be solved effi- K C ,
≤ (logN)6
ciently by linear programming (LP) techniques. This method
is therefore frequently referred to as l -LP reconstruction [3], where C depends only on the RIP constant.
1
[13], and its reconstruction complexity equals O m2N3/2 There exists an intimate connection between the LP recon-
wheninteriorpointmethodsareemployed[14].See[15],[16], struction accuracy and the RIP property, first described by
(cid:0) (cid:1)
[17] for other methodsto furtherreducethe complexityof l - Candés and Tao in [3]. If the sampling matrix Φ satisfies the
1
LP. RIP with constants δ , δ , and δ , such that
K 2K 3K
The reconstruction accuracy of the l -LP method is de-
1 δ +δ +δ <1, (1)
scribed in terms of the restricted isometry property (RIP), K 2K 3K
formally defined below. then the l -LP algorithm will reconstructall K-sparse signals
1
Definition 1 (Truncation): Let Φ Rm×N, x RN and exactly. This sufficient condition (1) can be improved to
∈ ∈
I 1, ,N . The matrix Φ consists of the columns of
⊂ { ··· } I δ <√2 1, (2)
Φ with indices i I, and xI is composed of the entries of x 2K −
∈
indexed by i I. The space spanned by the columns of ΦI as demonstrated in [18].
∈
is denoted by span(ΦI). Forsubsequentderivations,weneedtworesultssummarized
Definition 2 (RIP): A matrix Φ Rm×N is said to satisfy in the lemmas below. The first part of the claim, as well as a
∈ related modification of the second claim also appeared in [3],
theRestrictedIsometryProperty(RIP)withparameters(K,δ)
[10]. For completeness, we include the proof of the lemma in
1Weinterchangeably usebothnotations inthepaper. Appendix A.
4
Lemma 1 (Consequences of the RIP): III. THE SPALGORITHM
1) (Monotonicity of δK) For any two integers K K′, ThemainstepsoftheSPalgorithmaresummarizedbelow.2
≤
δK δK′.
≤ Algorithm 1 Subspace Pursuit Algorithm
2) (Near-orthogonalityofcolumns)LetI,J 1, ,N Input: K, Φ, y
⊂{ ··· }
be two disjointsets, I J =φ. Supposethatδ < Initialization:
|I|+|J|
1. For arbitrary vectors a R|I| and b R|J|, 1) T0 = K indices corresponding to the largest magni-
T ∈ ∈ {
tude entries in the vector Φ∗y .
|hΦIa,ΦJbi|≤δ|I|+|J|kak2kbk2, 2) yr0 =resid y,ΦTˆ0 . }
Iteration: At the ℓth iteration, go through the following steps
and (cid:0) (cid:1)
kΦ∗IΦJbk2 ≤δ|I|+|J|kbk2. 1) Tm˜ℓag=nitTudℓ−e1entr{iKes iinndthiceevseccotorrresΦp∗oyndℓ−in1g.to the largest
2) Set x =Φ†Sy. r
The lemma implies that δ δ δ , which conse- p T˜ℓ (cid:9)
K ≤ 2K ≤ 3K 3) Tℓ = K indicescorrespondingto the largestelements
quentlysimplifies(1)toδ <1/3.Both(1)and(2)represent {
3K of x .
p
sufficient conditions for exact reconstruction. }
4) yℓ =resid(y,Φ ).
In orderto describe the main steps of the SP algorithm,we 5) Ifr yℓ > yTℓ−ℓ1 , let Tℓ = Tℓ−1 and quit the
introduce next the notion of the projection of a vector and its r 2 r 2
iteration.
residue. (cid:13) (cid:13) (cid:13) (cid:13)
Output: (cid:13) (cid:13) (cid:13) (cid:13)
Definition 3 (Projection and Residue): Let y Rm and
∈ 1) The estimated signal xˆ, satisfying xˆ = 0
Φ Rm×|I|.SupposethatΦ∗Φ isinvertible.Theprojection {1,···,N}−Tℓ
ofIy∈onto span(ΦI) is defineId aIs and xˆTℓ =Φ†Tℓy.
y =proj(y,Φ ):=Φ Φ†y,
p I I I
A schematic diagram of the SP algorithm is depicted in
where Fig. 1(b). For comparison, a diagram of OMP-type methods
Φ† :=(Φ∗Φ )−1Φ∗ is provided in Fig. 1(a). The subtle, but important, differ-
I I I I
ence between the two schemes lies in the approach used to
denotesthe pseudo-inverseof the matrix Φ , and ∗ standsfor generate Tℓ, the estimate of the correct support set T. In
I
matrix transposition. OMP strategies, during each iteration the algorithm selects
The residue vector of the projection equals one or several indices that represent good partial support
set estimates and then adds them to Tℓ. Once an index is
yr =resid(y,ΦI):=y yp. includedinTℓ,itremainsinthissetthroughouttheremainder
−
of the reconstruction process. As a result, strict inclusion
Wefindthefollowingpropertiesofprojectionsandresidues rules are needed to ensure that a significant fraction of the
of vectors useful for our subsequent derivations. newly added indices belongs to the correct support T. On
Lemma 2 (Projection and Residue): the other hand, in the SP algorithm, an estimate Tℓ of size
K is maintained and refined during each iteration. An index,
1) (Orthogonality of the residue) For an arbitrary vector
whichisconsideredreliablein someiterationbutshownto be
y Rm, and a sampling matrix Φ Rm×K of full
∈ I ∈ wrongatalateriteration,canbeaddedtoorremovedfromthe
column rank, let y =resid(y,Φ ). Then
r I estimatedsupportsetatanystageoftherecoveryprocess.The
Φ∗y =0. expectationis thattherecursiverefinementsoftheestimate of
I r
the support set will lead to subspaces with strictly decreasing
2) (Approximation of the projection residue) Consider a distance from the measurement vector y.
matrix Φ Rm×N. Let I,J 1, N be two We performed extensive computer simulations in order to
∈ ⊂ { ··· }
disjoint sets, I J =φ, and suppose that δ < 1. compare the accuracy of different reconstruction algorithms
|I|+|J|
Furthermore,lety span(Φ ),y =proj(y,Φ )and empirically. In the compressive sensing framework, all sparse
I p J
T ∈
y =resid(y,Φ ). Then signalsare expectedtobe exactlyreconstructedas longasthe
r J
level of the sparsity is below a certain threshold. However,
δ
y |I|+|J| y , (3) the computational complexity to test this uniform reconstruc-
k pk2 ≤ 1 δmax(|I|,|J|) k k2 tion ability is O NK , which grows exponentially with K.
−
Instead,forempiricaltesting,weadoptthesimulationstrategy
and (cid:0) (cid:1)
described in [5] which calculates the empirical frequency of
δ
1 |I|+|J| y y y . (4)
− 1 δ k k2 ≤k rk2 ≤k k2 2In Step 3) of the SP algorithm, K indices with the largest correlation
(cid:18) − max(|I|,|J|)(cid:19) magnitudesareusedtoformTℓ.InCoSaMP[12],2Ksuchindicesareused.
This small difference results in different proofs associated with Step 3) and
The proof of Lemma 2 can be found in Appendix B. different RIPconstants thatguarantee successful signalreconstruction.
5
Reconstruction Rate (500 Realizations): m=128, N=256
1
0.9
n0.8
o
nstructi0.7
o
ec0.6
R
act 0.5
Ex
(a) Iterations in OMP, Stagewise OMP, and Regularized OMP: in each of 0.4
iteration, one decides on a reliable set of candidate indices to be added ncy
intothelistTℓ−1;onceacandidateisadded,itremainsinthelistuntilthe ue0.3
q
algorithm terminates. Fre0.2 Linear Programming (LP)
Subspace Pursuit (SP)
0.1 Regularized OMP
Standard OMP
0
0 10 20 30 40 50 60 70
Signal Sparsity: K
(a)SimulationsforGaussiansparsesignals:OMPandROMPstarttofailwhen
K≥19andwhenK≥22respectively, ℓ1-LPbeginstofailwhenK≥35,
andtheSPalgorithm fails onlywhenK≥45.
Reconstruction Rate (500 Realizations): m=128, N=256
1
(b)IterationsintheproposedSubspacePursuitAlgorithm:alistofK can-
0.9
didates,whichisallowedtobeupdatedduringtheiterations,ismaintained.
Figure 1: Description of reconstruction algorithms for K- on0.8
sparse signals: thoughboth approacheslook similar, the basic nstructi0.7
ideas behind them are quite different. eco0.6
R
act 0.5
Ex
of 0.4
exactreconstructionfortheGaussianrandommatrixensemble. ncy
The steps of the testing strategy are listed below. ue0.3
q
e
1) For given values of the parameters m and N, choose a Fr0.2 Linear Programming (LP)
Subspace Pursuit (SP)
signal sparsity level K such that K m/2; 0.1 Regularized OMP
2) Randomly generate a m N sampl≤ing matrix Φ from Standard OMP
× 0
the standard i.i.d. Gaussian ensemble; 0 10 20 30 40 50 60 70
Signal Sparsity: K
3) Select a support set T of size T = K uniformly at
random,andgeneratethesparses|ign|alvectorxbyeither (b) Simulations for zero-one sparse signals: both OMP and ROMP starts to
failwhenK≥10,ℓ1-LPbeginstofailwhenK≥35,andtheSPalgorithm
one of the following two methods: fails whenK≥29.
a) Draw the elements of the vector x restricted to T
Figure 2: Simulations of the exact recovery rate: compared
from the standard Gaussian distribution; we refer
with OMPs, the SP algorithm has significantly larger critical
to this type of signal as a Gaussian signal. Or,
sparsity.
b) set all entries of x supported on T to ones; we
refer to this type of signal as a zero-one signal.
Note that zero-one sparse signals are of special interest technique is plotted as well.
for the comparative study, since they represent a partic- Figure2depictstheempiricalfrequencyofexactreconstruc-
ularly challenging case for OMP-type of reconstruction tion. The numerical values on the x-axis denote the sparsity
strategies. level K, while the numerical values on the y-axis represent
4) Compute the measurement y = Φx, apply a recon- the fraction of exactly recovered test signals. Of particular
struction algorithm to obtain xˆ, the estimate of x, and interest is the sparsity level at which the recovery rate drops
compare xˆ to x; below100%-i.e.thecriticalsparsity-which,whenexceeded,
5) Repeat the process 500 times for each K, and then leadstoerrorsinthereconstructionalgorithmappliedtosome
simulate the same algorithm for different values of m of the signals from the given class.
and N.
The simulation results reveal that the critical sparsity of
The improved reconstruction capability of the SP method, the SP algorithm by far exceeds that of the OMP and ROMP
compared with that of the OMP and ROMP algorithms, is techniques, for both Gaussian and zero-one inputs. The re-
illustratedby two examplesshownin Fig. 2. Here, the signals construction capability of the SP algorithm is comparable to
aredrawnbothaccordingtotheGaussianandzero-onemodel, thatoftheLPbasedapproach:theSPalgorithmhasa slightly
and the benchmark performance of the LP reconstruction highercriticalsparsityforGaussiansignals,butalsoaslightly
6
lower critical sparsity for zero-one signals. However, the SP
algorithms significantly outperforms the LP method when
it comes to reconstruction complexity. As we analytically
demonstrate in the exposition to follow, the reconstruction
complexityoftheSPalgorithmforbothGaussianandzero-one
sparse signals is O(mNlogK), whenever K O √N ,
≤
while thecomplexityof LPalgorithmsbasedoninteri(cid:16)orpoi(cid:17)nt
methods is O m2N3/2 [14] in the same asymptotic regime.
IV(cid:0). RECOV(cid:1)ERY OF SPARSE SIGNALS
For simplicity, we start by analyzing the reconstruction Figure 3: After each iteration, a K-dimensional hyper-plane
performance of SP algorithms applied to sparse signals in closer to y is obtained.
the noiseless setting. The techniquesused in this context, and
the insights obtained are also applicable to the analysis of
SP reconstruction schemes with signal or/and measurement
whichimpliesthatateachiteration,theSPalgorithmidentifies
perturbations.Notethatthroughouttheremainderofthepaper,
a K-dimensional space that reduces the reconstruction error
we use the notation Si (S D,L , i Z+) stacked over
∈ { } ∈ ofthevectorx.SeeFig.3foranillustration.Thisobservation
an inequality sign to indicate that the inequality follows from
is formally stated as follows.
Definition(D) or Lemma (L) i in the paper.
Theorem 2: AssumethattheconditionsofTheorem1hold.
A sufficient condition for exact reconstruction of arbitrary
For each iteration of the SP algorithm, one has
sparse signals is stated in the following theorem.
Theorem 1: Let x RN be a K-sparse signal, and let
its corresponding mea∈surement be y = Φx Rm. If the kxT−Tℓk2 ≤cKkxT−Tℓ−1k2, (6)
∈
sampling matrix Φ satisfies the RIP with constant
and
δ3K <0.165, (5) yℓ cK yℓ−1 < yℓ−1 , (7)
r 2 ≤ 1 2δ r 2 r 2
thenthe SPalgorithmis guaranteedto exactlyrecoverx from − 3K
(cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13)
y via a finite number of iterations. where (cid:13) (cid:13) (cid:13) (cid:13) (cid:13) (cid:13)
Remark 3: TherequirementonRIPconstantcanberelaxed c = 2δ3K(1+δ3K). (8)
to K (1 δ )3
3K
−
δ <0.205,
3K
if we replace the stopping criterion yℓ yℓ−1 To prove Theorem 2, we need to take a closer look at the
r 2 ≤ r 2
with yℓ = 0. This claim is supported by substituting operationsexecutedduringeachiterationofthe SPalgorithm.
r 2 (cid:13) (cid:13) (cid:13) (cid:13)
δ < 0.205 into Equation (6). Howev(cid:13)er, (cid:13)for simp(cid:13)licity(cid:13)of During one iteration, two basic sets of computations and
a3nKalys(cid:13)(cid:13)is,w(cid:13)(cid:13)eadopt yℓ yℓ−1 fortheiterationstopping comparisons are performed: first, given Tℓ−1, K additional
criterion. r 2 ≤ r 2 candidateindicesforinclusionintotheestimateofthesupport
(cid:13)(cid:13) (cid:13)(cid:13) (cid:13)(cid:13) (cid:13)(cid:13) setareidentified;andsecond,givenT˜ℓ,K reliableindicesout
Remark 4: In the original version of this manuscript, we ofthetotal2K indicesareselectedtoformTℓ.InSubsections
proved the weaker result δ3K 0.06. At the time of revision IV-A and IV-B, we provide the intuition for choosing the se-
≤
of the paper, we were given access to the manuscript [19] by lectionrules. Now,letx be the residuesignalcoefficient
T−T˜ℓ
NeedelandTropp.Usingsomeoftheprooftechniquesintheir vector correspondingto the support set estimate T˜ℓ. We have
work, we managed to improve the results in Theorem 3 and
the following two theorems.
thereforethe RIP constant of the originalsubmission. The in-
Theorem 3: It holds that
terested reader is referred to http://arxiv.org/abs/0803.0811v2
for the first version of the theorem. This paper contains only 2δ
the proof of the stronger result. xT−T˜ℓ 2 ≤ (1 δ3K )2 kxT−Tℓ−1k2.
3K
−
(cid:13) (cid:13)
This sufficientconditionis provedby applyingTheorems2 (cid:13) (cid:13)
and 6. The computationalcomplexityis related to the number The proof of the theorem is postponed to Appendix D.
ofiterationsrequiredforexactreconstruction,andisdiscussed Theorem 4: The following inequality is valid
attheendofSectionIV-C.Beforeprovidingadetailedanalysis
of the results, let us sketch the main ideas behind the proof. 1+δ
x 3K x .
We denote by xT−Tℓ−1 and xT−Tℓ the residual signals k T−Tℓk2 ≤ 1 δ3K T−T˜ℓ 2
−
based upon the estimates of supp(x) before and after the (cid:13) (cid:13)
(cid:13) (cid:13)
ℓth iteration of the SP algorithm. Provided that the sampling
The proof of the result is deferred to Appendix E.
matrix Φ satisfies the RIP with constant (5), it holds that
BasedonTheorems3and4,onearrivesattheresultclaimed
kxT−Tℓk2 <kxT−Tℓ−1k2, in Equation (6).
7
Furthermore, according to Lemmas 1 and 2, one has to analytically justify due to the following fact. Although it
is clear that for all indices j / T, the values of v ,y are
yℓ = resid(y,Φ ) ∈ |h j i|
r 2 k Tℓ k2 upper bounded by δK+1 x , it may also happen that for all
k k
(cid:13)(cid:13) (cid:13)(cid:13) =kresid(ΦT−TℓxT−Tℓ,ΦTℓ) i∈T, the values of |hvi,yi| are small as well. Dealing with
+resid(Φ x ,Φ ) maximumcorrelationsin this scenario cannotbe immediately
Tℓ Tℓ Tℓ k2 proved to be a good reconstruction strategy. The following
D3
= kresid(ΦT−TℓxT−Tℓ,ΦTℓ)+0k2 example illustrates this point.
(4) Example 1: Without loss of generality, let
≤ kΦT−TℓxT−Tℓk2 T ={1,··· ,K}. Let the vectors vi (i∈ T) be orthonormal,
(6) andlettheremainingcolumnsvj,j / T,ofΦ beconstructed
∈
≤ 1+δK ·cKkxT−Tℓ−1k2, (9) randomly, using i.i.d. Gaussian samples. Consider the
following normalized zero-one sparse signal
where the secondpequality holds by the definition of the
1
residue,while(4)and(6)refertothelabelsoftheinequalities y= v .
i
used in the bounds. In addition, √K
i∈T
X
yrℓ−1 2 =kresid(y,ΦTℓ−1)k2 Then, for K sufficiently large,
(cid:13)(cid:13) (cid:13)(cid:13) =(4)kresid(ΦT−Tℓ−1xT−Tℓ−1,ΦTℓ−1)k2 |hvi,yi|= √1K ≪1, forall1≤i≤K.
1 δ δ
≥ −1K δ−K 2K kΦT−Tℓ−1xT−Tℓ−1k2 It is straightforwardto envision the existence of an index j ∈/
− T, such that
1 2δ
≥ 1−−δK2K 1−δKkxT−Tℓ−1k2 |hvj,yi|≈δK+1 > √1K.
1 2δ p
≥ √−1 δ2K kxT−Tℓ−1k2. (10) The latter inequality is critical, because achieving very small
− K values for the RIP constant is a challenging task.
Uponcombining(9)and(10),oneobtainsthefollowingupper
This example represents a particularly challenging case for
bound
the OMP algorithm. Therefore, one of the major constraints
1 δ2
yℓ − Kc yℓ−1 imposed on the OMP algorithm is the requirement that
r 2 ≤ 1 2δ K r 2
p− 2K 1
(cid:13)(cid:13) (cid:13)(cid:13) L1 1 c (cid:13)(cid:13)yℓ−1(cid:13)(cid:13) . mi∈aTx|hvi,yi|= √K >mj∈/aTx|hvj,yi|≈δK+1.
≤ 1 2δ K r 2
3K
Finally,elementarycalcula−tionsshow(cid:13)(cid:13)thatw(cid:13)(cid:13)henδ3K ≤0.165, TwohimcheedtecthaiyssrfeaqsutiaresmKenitn,cδrKea+se1s.has to be less than 1/√K,
c
K In contrast, the SP algorithm allows for the existence of
<1,
1 2δ3K some index j / T with
− ∈
which completes the proof of Theorem 2.
max v ,y < v ,y .
i j
i∈T |h i| |h i|
A. Why Does Correlation Maximization Work for the SP As long as the RIP constant δ3K is upper bounded by the
Algorithm? constant given in (5), the indices in the correct support of
x, that account for the most significant part of the energy
Both in the initialization step and during each iteration
of the signal, are captured by the CM procedure. Detailed
of the SP algorithm, we select K indices that maximize
descriptions of how this can be achieved are provided in the
the correlations between the column vectors and the residual
proofs of the previously stated Theorems 3 and 5.
measurement.Henceforth,thisstepisreferredtoascorrelation
Letusfirstfocusontheinitializationstep.Bythedefinition
maximization(CM).Considertheidealcasewhereallcolumns
of the set T0 in the initialization stage of the algorithm, the
of Φ are orthogonal3. In this scenario, the signal coefficients
set of the K selected columns ensures that
can be easily recoveredby calculating the correlations v ,y
i
- i.e., all indices with non-zero magnitude are in the hcorrecit Φ∗ y Φ∗y D2(1 δ ) x . (11)
support of the sensed vector. Now assume that the sampling k T0 k2 ≥k T k2 ≥ − K k k2
matrix Φ satisfies the RIP. Recall that the RIP (see Lemma Now, if we assume that the estimate T0 is disjoint from
1) implies that the columns are locally near-orthogonal.Con- the correct support, i.e., that T0 T = φ, then by the near
sequently, for any j not in the correct support, the magnitude orthogonality property of Lemma 1, one has
T
of the correlation v ,y is expected to be small, and more
precisely,upperbouhndjedibyδ x .Thisseemstoprovide kΦ∗T0yk2 =kΦ∗T0ΦTxTk2 ≤δ2Kkxk2.
K+1k k2
a very simple intuition why correlation maximization allows The last inequality clearly contradicts (11) whenever δ
K
≤
for exact reconstruction. However, this intuition is not easy δ <1/2. Consequently, if δ <1/2, then
2K 2K
3Ofcourse,inthiscasenocompressionispossible. T0 T =φ,
6
\
8
and at least one correct element of the support of x is in T0.
This phenomenon is quantitatively described in Theorem 5.
Theorem 5: After the initialization step, one has
1 δ 2δ
xT0TT 2 ≥ −1K+−δ 2K kxk2,
K
(cid:13) (cid:13)
and (cid:13) (cid:13)
8δ 8δ2
kxT−T0k2 ≤ 12+Kδ− 2K kxk2.
p 2K
The proof of the theorem is postponed to Appendix C.
To studytheeffectof correlationmaximizationduringeach
iteration, one has to observe that correlation calculations are
performed with respect to the vector
yrℓ−1 =resid(y,ΦTℓ−1)
Figure 4: The projection coefficient vector x′ is a smeared
instead of being performed with respect to the vector y. p
version of the vector xTTT′.
As a consequence, to show that the CM process captures a
significant part of residual signal energy requires an analysis
including a number of technical details. These can be found
changes in their magnitudes. Fortunately, the energy of this
in the Proof of Theorem 3. smear, i.e., ǫ , is proportional to the norm of the residual
k k2
signal x , which can be proved to be small according to
T−T˜ℓ
B. Identifying Indices Outside of the Correct Support Set the analysis accompanying Theorem 3. As long as the smear
atNleoastet Kthaotfththereemadreo2nKotibnedloicnegstiontthheecsoertrTe˜cℓt,saumppoonrgtwsehtiTch. igsoondotesstiemveartee,oxfpT ≈T˜xℓT˜vℓi,aothneelasrhgoeusltdprboejecatbiolenctooeoffibctaieinntsa.
In order to expurgate those indices from T˜ℓ, or equivalently, Thisintuitiveexplanationisformalizedinthepreviouslystated
T
Theorem 4.
in order to find a K-dimensional subspace of the space
span(Φ )closesttoy,weneedtoestimatetheseK incorrect
T˜ℓ
indices. C. Convergence of the SP Algorithm
Define ∆T :=T˜ℓ Tℓ−1. This set contains the K indices
Inthissubsection,we upperboundthenumberofiterations
−
which are deemed incorrect. If ∆T T =φ, our estimate of
needed to reconstruct an arbitrary K-sparse signal using the
incorrectindicesisperfect.However,sometimes∆T T =φ.
SP algorithm.
T 6
This means that among the estimated incorrect indices, there
Given an arbitrary K-sparse signal x, we first arrange its
T
aresomeindicesthatactuallybelongtothecorrectsupportset
elements in decreasing order of magnitude. Without loss of
T. The question of interest is how often these correct indices
generality, assume that
are erroneously removed from the support estimate, and how
quickly the algorithm manages to restore them back. x x x >0,
1 2 K
| |≥| |≥···≥| |
We claim thatthe reductioninthe normintroducedby
sucherroneousexpurgationissmall. Tk·hke2intuitiveexplanation and that xj =0, ∀j >K. Define
for this claim is as follows. Let us assume that all the min x
indicesinthesupportofxhavebeensuccessfullycaptured,or ρ := |xK| = 1≤i≤K | i|. (12)
equivalently, that T T˜ℓ. When we project y onto the space min kxk2 K x2
⊂ i=1 i
span(Φ ), it can be shown thatits correspondingprojection
T˜ℓ Letn denotethe numberof iteqraPtionsof the SP algorithm
coefficient vector x satisfies it
p
needed for exact reconstruction of x. Then the following
x =x , theorem upper bounds n in terms of c and ρ . It can be
p T˜ℓ it K min
viewed as a bound on the complexity/performance trade-off
and that it contains at least K zeros. Consequently, the K
for the SP algorithm.
indices with smallest magnitude - equal to zero - are clearly
Theorem 6: The number of iterations of the SP algorithm
not in the correct support set.
is upper bounded by
However, the situation changes when T * T˜ℓ, or equiva-
lently, when T T˜ℓ =φ. After the projection, one has logρmin 1.5 K
− 6 nit min − +1, · .
≤ logc logc
x =x +ǫ (cid:18) − K − K(cid:19)
p T˜ℓ
for some nonzero ǫ RT˜ℓ . View the projection coefficient This result is a combination of Theorems 7 and (12),4
| |
∈ described below.
vector x as a smeared version of x (see Fig. 4 for
p T˜ℓ
illustration): the coefficients indexed by i ∈/ T may become 4TheupperboundinTheorem7isalsoobtained in[12]while theonein
non-zero; the coefficients indexed by i T may experience Theorem8isnot.
∈
9
Theorem 7: One has the idea behind this example, another upper bound on n is
it
logρ described in Theorem 8 and proved in Appendix F.
min
nit − +1. Itisclearthatthenumberofiterationsrequiredforexactre-
≤ logc
K
− constructiondependsonthevaluesoftheentriesofthesparse
signal.We thereforefocusourattentiononthefollowingthree
Theorem 8: It can be shown that
particular classes of sparse signals.
1.5 K
nit · . 1) Zero-one sparse signals. As explained before, zero-one
≤ logc
− K signals represent the most challenging reconstruction
category for OMP algorithms. However, this class of
The proof of Theorem 7 is intuitively clear and presented
signals has the best upper bound on the convergence
below, while the proof of Theorem 8 is more technical and
rate of the SP algorithm.Elementarycalculationsreveal
postponed to Appendix F.
that ρ =1/√K and that
Proof of Theorem 7: The theorem is proved by contra- min
diction. Consider Tℓ, the estimate of T, with logK
n .
it
logρmin ≤ 2log(1/cK)
l= − +1 .
logc 2) Sparse signals with power-law decaying entries (also
(cid:24) − K (cid:25)
known as compressible sparse signals). Signals in this
Suppose that T *Tℓ, or equivalently, T Tℓ =φ. Then
− 6 category are defined via the following constraint
x = x2 x c i−p,
k T−Tℓk2 i | i|≤ x·
si∈XT−Tℓ for some constants c > 0 and p > 1. Compressible
x
(12)
sparse signals have been widely considered in the CS
min x = ρ x .
≥ i∈T | i| mink k2 literature, since most practical and naturally occurring
However, according to Theorem 2, signals belong to this class [13]. It follows from Theo-
rem 7 that in this case
x (c )ℓ x
k T−Tℓk2 ≤ K k k2 plogK
<ρminkxk2, nit ≤ log(1/cK)(1+o(1)),
where the last inequality follows from our choice of ℓ such where o(1) 0 when K .
that (cK)ℓ < ρmin. This contradicts the assumption T * Tℓ 3) Sparse signa→lswith expone→ntia∞llydecayingentries. Sig-
and therefore proves Theorem 7. nals in this class satisfy
AdrawbackofTheorem7isthatitsometimesoverestimates x c e−pi, (13)
i x
the number of iterations, especially when ρ 1. The | |≤ ·
min
example to follow illustrates this point. ≪ forsomeconstantscx >0andp>0.Theorem6implies
Example 2: Let K = 2, x = 210, x = 1, x = = that
1 2 3
···
xN =0.SupposethatthesamplingmatrixΦsatisfiestheRIP pK (1+o(1)) if 0<p 1.5
with cK = 12. Noting that ρmin . 2−10, Theorem 6 implies nit ≤(log1(.15/KcK) if p>1.5≤ ,
that log(1/cK)
nit 11. where again o(1) 0 as K .
≤ → →∞
Simulation results, shown in Fig. 5, indicate that the above
Indeed,ifwetakeacloselookatthestepsoftheSPalgorithm,
analysis gives the right order of growth in complexity with
we can verify that
respect to the parameter K. To generate the plots of Fig.
n 1.
it
≤ 5, we set m = 128, N = 256, and run simulations for
After the initialization step, by Theorem 5, it can be shown different classes of sparse signals. For each type of sparse
that signal, we selected different values for the parameter K, and
for each K, we selected 200 different randomly generated
Gaussian sampling matrices Φ and as many different support
8δ 8δ2
kxT−T0k2 ≤ 12+Kδ− 2K kxk2 <0.95kxk2. sets T. The plots depict the average number of iterations
p 2K versus the signal sparsity level K, and they clearly show that
As a result, the estimate T0 must contain the index one and n = O(log(K)) for zero-one signals and sparse signals
it
kxT−T0k2 ≤1. After the first iteration, since with coefficients decaying according to a power law, while
n = O(K) for sparse signals with exponentially decaying
kxT−T1k2 ≤cKkxT−T0k<0.95<mi∈iTn|xi|, coitefficients.
With the bound on the number of iterations required for
we have T T1.
⊂ exactreconstructionat hand,the computationalcomplexityof
This example suggests that the upper bound (7) can be the SP algorithm can be easily estimated: it equals the com-
tightened when the signal components decay fast. Based on plexityofoneiterationmultipliedbythenumberofiterations.
10
m=128, N=256, # of realizations=200 The case that requires special attention during analysis
Zero−one sparse signal is K2 > O(N). Again, if compressible sparse signals are
6
Power law decaying sparse signal: p=2 considered, the complexity of projections can be significantly
Exponentially decaying sparse signal: p=log(2)/2
reducedifonereusestheresultsfrompreviousiterationsatthe
nit5
ns : O(K) current iteration. If exponentially decaying sparse signals are
ber of Iteratio 34 O(log(K)) cmthoeonsstsidisgeinrgeandli,fiacosannentompisaaeryt—owfarntehtdetuocsiegonnthalyleareenfcdfoevctreteirvaetthteshiegennraeelrsgisdeputaiacrlasilotlyyf
m to K′ K. In both cases, the complexity depends on the
u
N ≪
e specific implementation of the CM and projection operations
ag 2
er and is beyond the scope of analysis of this paper.
Av
One advantage of the SP algorithm is that the number of
1
iterations required for recovery is significantly smaller than
that of the standard OMP algorithm for compressible sparse
0
0 5 10 15 20 25 30 signals. To the best of the authors’ knowledge, there are no
K
known results on the number of iterations of the ROMP and
Figure 5: Convergence of the subspace pursuit algorithm for StOMPalgorithmsneededforrecoveryofcompressiblesparse
different signals. signals.
V. RECOVERY OF APPROXIMATELY SPARSE SIGNALS
In each iteration, CM requires mN computations in general. FROM INACCURATE MEASUREMENTS
Forsomemeasurementmatriceswithspecialstructures,forex-
We first consider a sampling scenario in which the signal
ample,sparsematrices,thecomputationalcostcanbereduced
x is K-sparse, but the measurement vector y is subjected to
significantly. The cost of computing the projections is of the
an additive noise component, e. The following theorem gives
order of O K2m , if one uses the Modified Gram-Schmidt
a sufficient condition for convergence of the SP algorithm in
(MGS) algorithm [20, pg. 61]. This cost can be reduced
furtherby“r(cid:0)eusing(cid:1)”thecomputationalresultsofpastiterations termsofthe RIP constantδ3K, aswell asan upperboundson
the recovery distortion that depends on the energy (l -norm)
withinfutureiterations.Thisispossiblebecausemostpractical 2
of the error vector e.
sparse signals are compressible, and the signal support set
Theorem 9 (Stability under measurement perturbations):
estimates in different iterations usually intersect in a large
Let x RN be such that supp(x) K, and let its
number of indices. Though there are many ways to reduce ∈ | | ≤
correspondingmeasurementbe y=Φx+e, where e denotes
the complexity of both the CM and projection computation
the noise vector. Suppose that the sampling matrix satisfies
steps,weonlyfocusonthemostgeneralframeworkoftheSP
the RIP with parameter
algorithm, and assume that the complexity of each iteration
equals O mN +mK2 . As a result, the total complexity δ <0.083. (14)
3K
of the SP algorithm is given by O m N +K2 logK for
(cid:0) (cid:1) ThenthereconstructiondistortionoftheSPalgorithmsatisfies
compressible sparse signals, and it is upper bounded by
(cid:0) (cid:0) (cid:1) (cid:1)
O m N +K2 K for arbitrary sparse signals. When the x xˆ c′ e ,
signal is very sparse, in particular, when K2 O(N), the k − k2 ≤ Kk k2
tot(cid:0)al c(cid:0)omplexity(cid:1) of(cid:1)SP reconstruction is uppe≤r bounded by where
1+δ +δ2
O(mNK) for arbitrary sparse signals and by O(mNlogK) c′ = 3K 3K.
K δ (1 δ )
for compressible sparse signals (we once again point out that 3K − 3K
most practical sparse signals belong to this signal category
The proof of the theorem is given in Section V-A.
[13]).
We also study the case where the signal x is only approxi-
ThecomplexityoftheSPalgorithmiscomparabletoOMP-
mately K-sparse, and the measurement y is contaminated by
type algorithms for very sparse signals where K2 O(N).
≤ a noise vector e. To simplify the notation, we henceforth use
For the standard OMP algorithm, exact reconstructionalways
x to denote the vector obtained from x by maintaining the
requiresK iterations.Ineachiteration,theCMoperationcosts K
K entries with largest magnitude and setting all other entries
O(mN) computations and the complexity of the projection
in the vector to zero. In this setting, a signal x is said to be
is marginal compared with the CM. The corresponding total
approximately K-sparse if x x = 0. Based on Theorem
complexity is therefore always O(mNK). For the ROMP − K 6
9, we can upper boundthe recoverydistortion in terms of the
and StOMP algorithms, the challenging signals in terms of
l and l norms of x x and e, respectively, as follows.
convergenceratearealsothesparsesignalswithexponentially 1 2 − K
Corollary 1: (Stability under signal and measurement per-
decayingentries.Whenthepin(13)issufficientlylarge,itcan
turbations) Let x RN be approximately K-sparse, and let
be shown that both ROMP and StOMP also need O(K) iter- ∈
y = Φx+e. Suppose that the sampling matrix satisfies the
ations for reconstruction. Note that CM operation is required
RIP with parameter
inbothalgorithms.Thetotalcomputationalcomplexityisthen
O(mNK). δ <0.083.
6K