This paper discusses a new greedy algorithm for solving the sparse approximation problem over quasi-incoherent  dictionaries. These dictionaries consist of waveforms that are uncorrelated "on average," and they provide a natural generalization of incoherent dictionaries. The algorithm provides strong guarantees on the quality of the approximations it produces, unlike most other methods for sparse approximation. Moreover, very efficient implementations are possible via approximate nearest-neighbor data structure

Gilbert, A. C.

Muthukrishnan, S.

Strauss, M. J.

Tropp, J. A.

English

Caltech Authors - Main

IMPROVED SPARSE APPROXIMATION OVER QUASI-INCOHEmNT DICTIONARIES 
J A .  Ropp' A .  C. Gilbed S. Mtrthiilcrishnnd M. J Str.auss§ 
ABSTRACT 
This paper discusses a new greedy algorithm for solving the 
sparse approximation problem over quasi-incoherent dictio- 
naries. These dictionaries consist of waveforms that are un- 
correlated "on average," and they provide a natural general- 
ization of incoherent dictionaries. The algorithm provides 
strong guarantees on the quality of the approximations it 
produces, unlike most other methods for sparse approxima- 
tion. Moreover, very efficient implementations are possible 
via approximate nearest-neighbor data structures. 
1. INTRODUCTION 
Sparse approximation is the problem of finding a concise 
representation of a given signal as a linear comhination of a 
few elementary signals chosen from a rich collection. It has 
shown empirical promise in image processing tasks such as 
feature extraction, because thc approximation cannot suc- 
ceed unless it discovers structure latent in the image. For 
example, Starck, Donoho and Candts have used sparse ap- 
proximation to extract features from noisy astronomical pho- 
tograph and volumetric data [I ] .  Nevertheless, it has been 
difficult to estahlish that proposed algorithms actually solve 
the sparse approximation problem. This paper makes an- 
other step in that direction by describing a greedy algorithm 
that computes solutions with provable quality guarantees. 
A dicrionay 9 for the signal space Rd is a collection of 
vectnrs that spans the entire space. The vectors are called 
atonrs, and we write them as 'PA. The index X may parame- 
terize the timekcale or time/frequency localization of each 
atom, or it may he a label without any additional meaning. 
The number of atoms is often much larger than the signal 
dimension. 
The sparse approxinialion problem with respect to 59 is 
to compute a good representation of each input signal as 
a short linear comhination of atoms. Specifically, for an 
'Inst. for Comp. Eng. and Sci. (ICES), The University of Texas at 
Austin, Austin, TX 78712, j t r o p p a i c e s . u t e x a s  .edu. 
tAT&T Labs-Research. 180 Park Avenue, Flurham Park, NJ 07932. 
a9ilbertaresearch.att.com 
iAT&l Labs-Research & Rutgers University, 180 Park Avenue, 
Plorllam Park, NJ 07932. muthuoresearch.att .con. Supported in 
part by NSF CCR 00-87022 and NSF 1 l R  0220280. 
gAT&T Labs-Rzsearch, 180 Park Avenue, Florham Park. NI 07932. 
mstrauss~research.att.com. 
arbitrary signal z, we search for an na-term superposition 
aopt = h 'PA 
A, 
which minimizes 111 - aupfIJp. We must determine both the 
optimal vectors, 711 atoms whose indices are listed by 
as well as the optimal coefficients bA. 
If 9 is an orthonormal basis, it is computationally easy 
to find aopt. For the indices AOpt, simply take ni atoms with 
the largest inner products /(I, pA)l and form 
A, 
Unfortunately, it  can he difficult or impossible to choose 
an appropriate orthonormal basis for a given situation. For 
example, if the signals contain both harmonic and impulsive 
components, a single onhonormal basis will not represent 
them both efficiently. We have much more freedom with a 
redundant dictionary, since it may include a rich collection 
of waveforms which can provide concise representations of 
many different structures. 
The price that we pay for additional flexibility is an 
increased cost to determine these concise representations. 
For general redundant dictionaries, it is computationally in- 
feasible to search all possible m-term representations. In  
fact, if 9 is an arbitrary dictionary, finding the hest m- 
term representation of an arbitrary signal is NP-hard [2]. 
There are algorithms with provable approximation guaran- 
tees for specific dictionaries, e.g. Villemoes' algorithm for 
Haar wavelet packets [3]. There are also some well-known 
heuristics, such as Matching Pursuit (MP) [4], Onhogo- 
nal Matching Pursuit (OMP) [SI and m-fold Matching Pur- 
suit [6]. Several other methods rely on the Basis Pursuit 
paradigm, which advocates minimizing the C, norm of the 
coefficients in the representation instead of minimizing the 
sparsity directly [7]. 
Some theoretical progress has already been made for 
dictionaries with low coherence. The coherence parame- 
ter ) I  equals the maximal inner product hetween two dis- 
tinct atoms. For example, the union of spikes and sines is 
a dictionary with p = m. The authors in [6] have pre- 
sented an efficient two-stage algorithm for the approximate 
representation of any signal over a sufficiently incoherent 
0-7803-7750-8/03/$17.00 02003 IEEE 1 - 3 7  
dictionary. This is the first known algorithm which prov- 
ably approximates the solution to the sparse problem for any 
class of general dictionaneS. In addition, this algorithm is 
highly efficient. For a suitably incohcrent dictionaty, it is 
also known that Basis Pursuit can resolve the subclass of 
signals which have an a a c t  sparse representation [SI. 
This article offers a number of improvements to [h]. 
Specifically, we present a modified version of the algorithm 
in [h], which calculates significantly more accurate sparse 
representations for incoherent dictionaries and also applies 
to a much larger class of redundant dictionaries. Unlike an 
incoherent dictionary where all the inner products are small, 
the dictionaries we consider only need to have small inner 
products "on average." In addition, our analysis is simpler, 
Of course, the new algorithm can he implemented just as 
efficiently as the ones in [GI. 
2. ALGORITHM AND ANALYSIS 
2.1. Overall results 
For an incoherent dictionary, we have 
Theorem 1 Fi.r a dictionarv 9 M~illi coherence /i, 14'2 seek 
U N  m-term rrpirsentution of  an arbitran signal x, Mherr 
ni < +I$,- ' .  There is air ulgorirlrai rharpiudiices an I I I - ~ ~ ~ W  
rrpresentalion a,,, for x with enur 
In  comparison, the algorithm of [6] requires that 111 < 
- :;2 /L-' and produces approximations with error 
When rii 5 ~ i - " ~ ,  the resulting constant of approximation 
is ahout 46. Meanwhile, the algorithm descrihed here pro- 
duces approximatitins with error 
IIz - a,, 1 1 2  i 3 11% - aol,t.l12 
so long as 4 5 'in 6 / A - ' / * .  
Theorem 1 is a special case of a result for general dic- 
tionaries. To state the full theorem, we need to  horrow a 
definition from [9]. Let @he a matrix whose columns are 
the atoms of 9. Then the Gram matrix G %' Q* Q contains 
all the inner products hetween atoms!. Let I GI denote its en- 
ttywise ahsolute vdluc. Now, we define the Buhel,fitncrion 
~ ~ ( T J I )  of the dictionaty to he the maximum sum of any 111. 
(,nondiagonal) elements from a single row of IGI. In other 
words, the Babel function quantifies thc maximum total co- 
herence between a fixed atom and a collection of I I I  other 
'Here and elsewhere, dcnolesthe coujugatc transpose. 
atoms'. The B a k l  function is a more subtle way of de- 
scribing the dictionary than the coherence, since coherence 
only reflects the largest inner product. Clearly, 
p I ( i 7 " )  5 p111. (1  1 
Tlut is, the cumulative coherence always dominates the Ba- 
he1 function. When the Babel function grows slowly, we say 
informally that the dictionary is qirasi-incoliermf. 
Theorem 2 So long CIS p1 ( i n )  < f, oiir algorithm pie-  
drrces an ni-terrri app~nrireation a,,, which sati;fies 
where aopt is the opti~~ial  in-remi appm.ririration 
Obviously, Theorem 1 follows directly from Theorem 2 
by application of the hound ( : I ) .  
We can easily construct a dictionary for which we need 
the more general theorem. Let each atom he a linear com- 
bination of two impulses: 
(PI; = 9 6 k  + +6at-1 fo rk  = 1, .  . . , d .  
Then thc coherence / L  = e, which means that Theorem 1 
applies only when 111 5 3. Meanwhile, the Babel function 
p t ( n i )  = e < foreveryo,  2 2. Therefore,thegeneral 
theorem shows that approximation succeeds for any 111, and 
the error hound is 
111 - am/I2 5 112 - aoptl12 
Another consequence of Theorem 2 is that the algorithm 
can recover any signal which has an exact ?ii-tenn rcpresen- 
tation, so long as p1(111) < f. In fact, the analysis of191 
shows that Orthogonal Matching Pursuit alone can accom- 
plish this task. Donoho and Elad have proven that Basis 
Pursuit can recover exactly sparse signals under an identical 
condition [ IO] .  But they have not offered an approximation 
guarantee for general signals, and the algorithms associated 
with the Basis Pursuit paradigm are typically very slow. The 
algorithm descrihed here is significantly faster. 
2.2. A Structural Lemma 
An important ingedient in the analysis is our_eeneralization 
of Parseval's Theorem to an arbitrary dictionary. A similar 
lemma appears implicitly in the analysis of [6]. 
'The suhscnpt in the notalion serves to dstinpulsh tho Babel function 
fmm the coherence parameter and to remind us that i t  reprcsciits the CI 
n o m .  i.e. an ahsolute sum. 
1 - 3 8  
Lemma 3 (Approximate Parseval) Let A be a collection 
q f n t  aloiiis. Let P.h be the ortlrogonal pmjector o m  the 
span nf tlre aloira listed [JJJ A. For ei'ew signal x, we have 
Pmof Define the matrix @A whose columns are the m 
atoms listed in A.  Then the orthogoiial projector onto the 
span of the atoms in A can he written a. 
P* = (@;@A)-l'*@; 
Denote the smallest and largest singular values of @,\ hy 
umin and umilx. It follows that the smallest and largest 
eigenvalues of (@I @A)-'/' are respectively equal to 1 /U,,, 
and l / ~ ~ j , ~ .  Therefore, 
Since the squared singular values of @A are equal to the 
eigenvalucs of GA %' @;@A, we can estimate umin and 
cmax by examining this Gram matrix. The GerSgorin Disc 
Theorem [ I  11 states that every eigenvalue of a square matrix 
A lies in one of the ' i n  discs 
Since GA is a principal submdlrix of the dictionary Gram 
matrix G, we can use the Bahel function to hound the sum: 
C I ( G A ) , +  5 ~ 1 ( 7 1 ! )  
i#k 
All atoms have unit norm, so the diagonal of GI\ is identi- 
cally one. We may conclude that 
u l i ,  2 1 - p~(irr)  and 
5 1 + PI (m). 
Finishbywriting II@;zII; =EA I(x,ph)Iz. 0 
The Babel function of an orthonormal hasis is zero. In 
this case. Lemma 3 reduces to Parseval's Theorem. 
2.3. Analysis of two-phase greedy pursuit 
The overall algorithm is a two-phase greedy pursuit. First, 
wc initialize ao = 0 and perform Orthogonal Matching Pur- 
suit until we reach a K-term representation aK with a red- 
sonahle error guarantee. (The optimal error is necessary to 
determine K. For now, assume an oracle provides it. Sec- 
tion 2.4 discusses how to avoid the trip to Delphi.) In the 
second stage, we use (n! - K)-fold Matching Pursuit to ac- 
quire the remaining atoms. l h e  algorithm Eturns the best 
approximation to the signal over the ni chosen atoms. This 
procedure is different from the one given in [6], where the 
algorithm returns the sum of the hest approximation of the 
signal over the first K atoms and the best approximation of 
the li-term residual over the last (na - IC) atoms. 
Now, we sketch the analysis. We require the following. 
Theorem 4 191 For anv signal I, Ortliogonal Matching 
Pursuit can calculate an uppimiinant a K  that consists of 
0 5 h- 5 711 atoinspom the optinral ii!-tenii approxiinant 
aoFt and that satis/ies the ermr botrnd 
After the first phase is complete, we have selected K op- 
timal atoms, A,<, and the I<-term approximant satisfies the 
error bound (2). The second phase chooses E (111 .  - K )  
atoms that have the large? inner products with the residual. 
More precisely, we find K indices, Ac, to maximize 
A F  
We may assume that the atoms chosen in the second phase 
are distinct from those chosen in the first phase because the 
residual (z - a ~ )  is orthogonal to each atom that partie- 
ipates in the approximant a ~ .  The ni atoms we have se- 
lected are indexed hy A,,, d" AK U A+ and 
AS# 
Finally, the algorithm returns a,, which is the best approx- 
imation lo the signal using the atoms in Ant. 
Theorem 5 The vL-terni appmxinration pmducerl hv this 
two-phase alporitlinr satisje~v 
Pmqf The difference between the actual error and the opti- 
mal error is 
1 - 3 9  
since the atoms in A,,,, cany less energy than those in A,n. 
Applying I.emma 3 again, 
The first theorem provides a hound on (Iz - a ~ \ l i  in terms 
of llz - aoptll;. Comhining the two estimates and rearrang- 
ing, we reach the stated result. 0 
2.4. Implementation 
Implementing the algorithm which we have described re- 
quires foreknowledge of the optimal error. l he re  arc two 
ways to escape the need for omniscience. For the first op- 
tion, simply execute the algorithm (m -t 1) times, switching 
from the first phase to the second at each I< = 0,. .  in. 
Then select the hest of the representations. A second option 
is to guess the optimal error. This can he accomplished by 
running the algorithm with guesses taken from a geometric 
progression ranging from E ,  the machine precision, to I ~ Z ( / ~ .  
This requires only log, 1 1 ~ 1 1 ~  / E  attempts, and we may use 
the best of the representations. Both techniques are emhar- 
rassingly parallel, although efficient serial versions are also 
possible. 
Both phases of the algorithm require the determination 
of maximum inner products. At each step, the Orthogonal 
Matchina Pursuit phase chooses an atom with maximal in- 
ner product against the residual. The’second phase can be 
implemented hy selecting a maximal inner product, remov- 
ing that atom from the dictionary and iterating. ‘Therefore, 
in both stages, we use a data structure that preprocesses 
the dictionary and supports two queries: retum an atom 
whose ahsolute inner product with the’residual is maximal 
and delete an atom. An identity for unit vectors states that 
2 
(U, U) = 1 - 4 11% - V1l2 
To find the maximum ahsolute inner product between a sig- 
nal z and the dictionary, we can normalize the signal as 2 
and solve 
iniri [\(PA i ~ 1 1 ;  
This minimization can he performed approximately using 
a nearest-neighbor data stmcture for vectors under the Eu- 
clidean metric 1121. Building the data structure requires 
time and space poly(l9l/~1), where 7 is the precision re- 
quired in the approximation. But each query costs only 
d+polylog( 153 / I ] )  units oftime3. This query does not nec- 
essarily return a vector with the largest inner product, hul i t  
always returns a vector that is nearly as good. The anal- 
ysis of our algorithm changes slightly if we use a nearest- 
neighbor data structure to estimate maximum inner prod- 
ucts. The details are somewhat technical, so we will relegate 
them to a longer version of this article. 
3. REFERENCES 
[ I ]  J. L. Starck, D. I>. Donoho, and E. J. Candts, “Astro- 
nomical image representation by the Curvelet Trans- 
form,” Astronorii.v and Aslrvpli.vsics, 2002, Accepted. 
121 G. Davis, S. Mallat, and M. Avellaneda, “Greedy 
adaptive approximation,” Consti: Appmx., vol. 13, pp. 
57-98, 1997. 
[ 3 ]  L. E Villemoes, “Best approximation with walsh 
atoms,“ Consti: .4pplBX., vol. 13, pp. 329-355, 1997. 
[4] S. Mallat and Z.  Zhang, “Matching pursuits with time 
frequency dictionaries,” IEEE Trans. Signal Process- 
ingvol .41 ,  no. 12,pp. 3397-3415, 1993. 
[5]  G. Davis, S. Mallat, and Z. Zhang, “Adaptive time- 
[6] A. C. Gilbert, S. Muthukrishnan, and M. J. Strauss, 
“Approximation of functions over redundant dictio- 
naries using coherence,” in The 14/11 Annual ACM- 
SIAM $viiipofiimi on Disaule AlgoriAnis, Jan. 2003 
frequency decompositions,” Op/ica/ Eng., July 1 994. 
[71 S. S. Chen, D. L. Donoho, and M. A. Saunders, 
“Atomic decomposition by Basis Pursuit,“ SIAM .I 
Sci. Coirip., vol. 20, no. I ,  pp. 33-61, 1999, electronic. 
[El D. L.. Donoho and X. Hno, “Uncertainty principles 
and ideal atomic decomposition,” 
Tli., vol. 47, pp. 2845-2862, Nov. 2I)OI. 
[9] J. A. Tropp, “Greed is good: Algorithmic results for 
sparse approximation;. ICES Report 0304, The Uni- 
versity of Tcxas at Austin, 2003. 
1101 D. I>. Donoho and M. Elad, “Maximal sparsity rcpre- 
seutation hy Cl minimization,” I’roc. Nu//. Acud. Sf:;., 
vol. inn, pp. 2197-2202, MX. 2003. 
[I I1  R. A. Horn and C. R. Johnson, Ma/rix Analysis, Cam- 
bridge University Press, 1985. 
117, Ph.D. thesis, Stanford, 2000. 
[ 121 P. Indyk, High-Diniensinnal Corniptirurinnal Gennie- 
lSaying that a function IS poly(t) means that it lies in O ( l n )  for some 
c1 > 0. Similarly, each fiincliun in polylop(l.) lies in O(10g“ t) for some 
a > I). 
1 -40  


Improved sparse approximation over quasi-incoherent dictionaries

Approximation of functions over redundant dictionaries using coherence,”

Astronomical image representation by the Curvelet Transform,” Astronorii.v and Aslrvpli.vsics,

Greed is good: Algorithmic results for sparse approximation;. ICES Report 0304, The University of Tcxas at Austin,

Matching pursuits with time frequency dictionaries,”

Villemoes, “Best approximation with walsh atoms,“

http://authors.library.caltech.edu/9047/2/TROicip03.pdf

Improved sparse approximation over quasi-incoherent dictionaries

Abstract

Similar works

Full text

Available Versions

Caltech Authors - Main