Atomic Decomposition by Basis Pursuit.
The time-frequency and time-scale communities have recently developed a large number of overcomplete waveform dictionaries --- stationary wavelets, wavelet packets, cosine packets, chirplets, and warplets, to name a few.  Decomposition into overcomplete systems is not unique, and several methods for decomposition have been proposed, including the method of frames (MOF), Matching pursuit (MP), and, for special dictionaries, the best orthogonal basis (BOB).Basis Pursuit (BP) is a principle for decomposing a signal into an "optimal" superposition of dictionary elements, where optimal means having the smallest l1 norm of coefficients among all such decompositions.  We give examples exhibiting several advantages over MOF, MP, and BOB, including better sparsity and superresolution.  BP has interesting relations to ideas in areas as diverse as ill-posed problems, in abstract harmonic analysis, total variation denoising, and multiscale edge denoising.BP in highly overcomplete dictionaries leads to large-scale optimization problems.  With signals of length 8192 and a wavelet packet dictionary, one gets an equivalent linear program of size 8192 by 212,992. Such problems can be attacked successfully only because of recent advances in linear programming by interior-point methods.  We obtain reasonable success with a primal-dual logarithmic barrier method and conjugate-gradient solver.
Introduction
Over the last several years, there has been an explosion of interest in alternatives to traditional
signal representations. Instead of just representing signals as superpositions of
sinusoids (the traditional Fourier representation) we now have available alternate dictionaries
- collections of parameterized waveforms - of which the Wavelets dictionary is only
the best known. Wavelets, Steerable Wavelets, Segmented Wavelets, Gabor dictionaries,
Multi-scale Gabor Dictionaries, Wavelet Packets, Cosine Packets, Chirplets, Warplets, and
a wide range of other dictionaries are now available. Each such dictionary D is a collection
of waveforms (OE fl ) fl2\Gamma , with fl a parameter, and we envision a decomposition of a signal s
as
or an approximate decomposition
where R (m) is a residual. Depending on the dictionary, such a representation decomposes
the signal into pure tones (Fourier dictionary), bumps (wavelet dictionary), chirps (chirplet
dictionary), etc.
Most of the new dictionaries are overcomplete, either because they start out that way,
or because we merge complete dictionaries, obtaining a new mega-dictionary consisting of
several types of waveform (e.g. Fourier & Wavelets dictionaries). The decomposition (1.1)
is then nonunique, because some elements in the dictionary have representations in terms
of other elements.
1.1 Goals of Adaptive Representation
Nonuniqueness gives us the possibility of adaptation, i.e., of choosing among many representations
one that is most suited to our purposes. We are motivated by the aim of
achieving simultaneously the following
ffl Sparsity. We should obtain the sparsest possible representation of the object - the
one with the fewest significant coefficients.
ffl Superresolution. We should obtain a resolution of sparse objects that is much higher-resolution
than that possible with traditional non-adaptive approaches.
An important constraint, which is perhaps in conflict with both the goals:
ffl Speed. It should be possible to obtain a representation in order O(n) or O(n log(n))
time.
1.2 Finding a Representation
Several methods have been proposed for obtaining signal representations in overcomplete
dictionaries. These range from general approaches, like the Method of Frames [8], and the
method of Matching Pursuit [23], to clever schemes derived for specialized dictionaries, like
the method of Best Orthogonal Basis [6]. These methods are described briefly in Section
2.3.
In our view, these methods have both advantages and shortcomings. The principal
emphasis of the proposers of these methods is in achieving sufficient computational speed.
While the resulting methods are practical to apply to real data, we show below by computational
examples that the methods, either quite generally or in important special cases,
lack qualities of sparsity-preservation and of stable super-resolution.
1.3 Basis Pursuit
Basis Pursuit (BP) finds signal representations in overcomplete dictionaries by convex op-
timization: it obtains the decomposition that minimizes the ' 1 norm of the coefficients
occurring in the representation. Because of the non-differentiability of the ' 1 norm, this
optimization principle leads to decompositions that can have very different properties from
the Method of Frames - in particular they can be much sparser. Because it is based on
global optimization, it can stably super-resolve in ways that Matching Pursuit can not.
BP can be used with noisy data by solving an optimization problem trading off a
quadratic misfit measure with an ' 1 norm of coefficients. Examples show that it can stably
suppress noise while preserving structure that is well-expressed in the dictionary under
consideration.
BP is closely connected with linear programming. Recent advances in large-scale linear
programming - associated with interior-point methods - can be applied to BP, and make it
possible, with certain dictionaries, to nearly-solve the BP optimization problem in nearly-
linear time. We have implemented a primal-dual log barrier interior-point method as part
of a computing environment called Atomizer, which accepts any of a wide range of dictio-
naries. Instructions for Internet access of Atomizer are given in Section 6.6. Experiments
with standard time-frequency dictionaries indicate some of the potential benefits of BP.
Experiments with some nonstandard dictionaries - like the stationary wavelet dictionary
and the Heaviside dictionary - indicate important connections between BP and methods
like Mallat and Hwang's Multi-Scale Edge Representation and Osher, Rudin and Fatemi's
Total Variation-based De-Noising methods.
number of dictionaries and existing methods for overcomplete representation. In Section
3 we discuss the principle of Basis Pursuit and its relations to existing methods and to
ideas in other fields. In Section 4 we discuss methodological issues associated with BP -
in particular some of the interesting nonstandard ways it can be deployed. In Section 5 we
describe Basis Pursuit De-Noising, a method for dealing with problem (1.2). In Section 6
we discuss recent advances in large-scale linear programming, and resulting algorithms for
BP. In Section 7 we discuss a number of connections with other work.
Representations
n) be a discrete-time signal of length n; this may also be viewed as a
vector in R n . We are interested in the reconstruction of this signal using superpositions of
elementary waveforms. Traditional methods of analysis and reconstruction involve the use
of orthogonal bases, such as the Fourier basis, various discrete cosine transform bases, and
orthogonal wavelet bases. Such situations can be viewed as follows: given a list of n wave-
forms, one wishes to represent s as a linear combination of these waveforms. The waveforms
in the list, viewed as vectors in R n , are linearly independent, and so the representation is
unique.
2.1 Dictionaries and Atoms
A considerable focus of activity in the recent signal processing literature has been the
development of signal representations outside the basis setting. We use terminology introduced
by Mallat and Zhang [23]. A dictionary is a collection of parameterized waveforms
\Gamma). The waveforms OE fl are discrete-time signals of length n called atoms.
Depending on the dictionary, the parameter fl can have the interpretation of indexing
frequency, in which case the dictionary is a frequency or Fourier dictionary, of indexing
time/scale jointly, in which case the dictionary is a time-scale dictionary, or of indexing
time/frequency jointly, in which case the dictionary is a time-frequency dictionary. Usually
dictionaries are complete or overcomplete, in which case they contain exactly n atoms, or
more than n atoms, but one could also have continuum dictionaries containing an infinity of
atoms, and undercomplete dictionaries for special purposes, containing fewer than n atoms.
Dozens of interesting dictionaries have been proposed over the last few years; we focus in
this paper on a half dozen or so; much of what we do applies in other cases as well.
2.1.1 Trivial Dictionaries
We begin with some overly simple examples. The Dirac dictionary is simply the collection of
waveforms that are zero except in one point: This
is of course also an orthogonal basis of R n - the standard basis. The Heaviside dictionary is
the collection of waveforms that jump at one particular point:
1 ft-flg . Atoms in this dictionary are not orthogonal, but every signal has a representation
2.1.2 Frequency Dictionaries
A Fourier dictionary is a collection of sinusoidal waveforms OE fl indexed by
is an angular frequency variable and - 2 f0; 1g indicates phase type : sine or
cosine. In detail,
For the standard Fourier dictionary, we let fl run through the set of all cosines with Fourier
sines with Fourier frequencies
dictionary consists of n waveforms; it is in fact a basis, and a very
simple one : the atoms are all mutually orthogonal. An overcomplete Fourier dictionary is
obtained by sampling the frequencies more finely. Let ' be a whole number ? 1 and let
' be the collection of all cosines with sines with
1. This is an '-fold overcomplete system. We also use
below complete and overcomplete dictionaries based on discrete cosine transforms and sine
transforms.
2.1.3 Time-Scale Dictionaries
There are several types of Wavelet dictionary; to fix ideas, we consider the Haar dictionary,
with "Father Wavelet" . The dictionary
is a collection of translations and dilations of the basic mother wavelet, together with
translations of a father wavelet. It is indexed by is a scale
location, and - 2 f0; 1g indicates gender. In detail,
a:
For the standard Haar dictionary, we let fl run through the discrete collection of mother
wavelets with dyadic scales a locations that are integer
multiples of the scale b and the collection of father wavelets
at the coarse scale j 0 . This dictionary consists of n waveforms; it is an orthonormal basis.
An overcomplete wavelet dictionary is obtained by sampling the locations more finely : one
location per sample point. This gives the so-called Stationary Haar dictionary, consisting
of O(n log 2 (n)) waveforms. It is called stationary since the whole dictionary is invariant
under circulant shift.
A variety of other wavelet bases are possible. The most important variations are smooth
wavelet bases, using splines or using wavelets defined recursively from two-scale filtering
relations [9]. Although the rules of construction are more complicated (boundary conditions
[25], orthogonality versus bi-orthogonality [9], etc.), these have the same indexing structure
as the standard Haar dictionary. In this paper, we use Symmlet-8 smooth wavelets, i.e.,
Daubechies Nearly Symmetric wavelets with eight vanishing moments; see [9] for examples.
2.1.4 Time-Frequency Dictionaries
Much recent activity in the wavelet communities has focused on the study of time-frequency
phenomena. The standard example, the Gabor dictionary, is due to Gabor (1946); in our
notation, we take is a frequency, - is a location, ' is a phase,
and ffit is the duration, and consider atoms OE
Such atoms indeed consist of frequencies near ! and essentially vanish far away from - .
For fixed ffi t, discrete dictionaries can be built from time-frequency lattices,
with \Delta- and \Delta! chosen sufficiently fine these are complete.
For further discussions see e.g. [8].
Recently, Coifman and Meyer [5] developed the wavelet packet and cosine packet dictionaries
especially to meet the computational demands of discrete-time signal processing.
For 1-d discrete time signals of length n, these dictionaries each contain about n log 2 (n)
waveforms. A wavelet packet dictionary includes, as special cases, a standard orthogonal
wavelets dictionary, the Dirac dictionary, and a collection of oscillating waveforms spanning
a range of frequencies and durations. A cosine packet dictionary contains, as special cases,
the standard orthogonal Fourier dictionary, and a variety of Gabor-like elements sinusoids
of various frequencies weighted by windows of various widths and locations.
In this paper, we often use wavelet packet and cosine packet dictionaries as examples of
overcomplete systems, and we give a number of examples decomposing signals into these
Time
Frequency
(b) Phase Plane
(a) Frequency Domain
|FFT(WaveletPacket(3,3,7))|
Frequency
-0.50.5(c) Time Domain
Time

Figure

2.1: Time-frequency phase plot of a wavelet packet atom.
time-frequency dictionaries. A simple block-diagram helps us visualize the atoms appearing
in the decomposition. This diagram, adapted from Coifman and Wickerhauser [6], associates
with each cosine packet or wavelet packet a rectangle in the time-frequency phase
plane. The association is illustrated in Figure 2.1 for a certain wavelet packet. When a
signal is a superposition of several such waveforms, we indicate which waveforms appear in
the superposition by shading the corresponding rectangles in the time-frequency plane.
2.1.5 Further Dictionaries
We can always merge dictionaries to create mega-dictionaries; examples used below include
mergers of Wavelets with Heavisides.
2.2 Linear Algebra
Suppose we have a discrete dictionary of p waveforms and we collect all these waveforms
as columns of an n by p matrix \Phi, say. The decomposition problem (1.1) can be written
is the vector of coefficients in (1.1). When the dictionary furnishes a
basis, then \Phi is an n by n nonsingular matrix and we have the unique representation
When the atoms are, in addition, mutually orthonormal, then \Phi
the decomposition formula is very simple.
An important (but trivial) comment. Given a dictionary of waveforms, one can distinguish
analysis from synthesis. Synthesis is the operation of building up a signal by
superposing atoms; it involves a matrix that is n by Analysis involves the
operation of associating with each signal a vector of coefficients attached to atoms; it involves
a matrix that is p by n: a = \Phi T s. Synthesis and analysis are very different linear
operations, and we must take care to distinguish them. One should avoid assuming that
the analysis operator ~
us coefficients that can be used as is to synthesize s.
(a) Signal: Carbon (b) Synthesis Phase Plane
Time
Frequency
(c) Analysis Phase Plane
Time
Frequency
Sorted Coefficients
Order
Amplitude
Synthesis: Solid
Analysis: Dashed

Figure

2.2: Analysis versus synthesis of the signal Carbon.
In the overcomplete case we are interested in, p AE n and \Phi is not invertible. There are
then many solutions to (2.2), and a given approach selects a particular solution. One does
not uniquely and automatically solve the synthesis problem by applying a simple, linear,
analysis operator.
We now illustrate the difference between synthesis
Panel 2.2a shows the signal Carbon. Panel 2.2b shows the time-frequency structure of a
sparse synthesis of Carbon, a vector ff yielding using a wavelet packet dictionary.
To visualize the decomposition, we present a phase-plane display with shaded rectangles,
as described above. Panel 2.2c gives an analysis of Carbon, the coefficients a
displayed in a phase-plane. Once again, between analysis and synthesis there is a large
difference in sparsity. In Panel 2.2d we compare the sorted coefficients of the overcomplete
representation (synthesis) with the analysis coefficients.
2.3 Existing Decomposition Methods
There are several currently popular approaches to obtaining solutions to (2.2).
2.3.1 Frames
The Method of Frames (MOF) [8] picks out, among all solutions of (2.2), one whose coefficients
have minimum l 2 norm:
subject to s: (2.3)
The solution of this problem is unique; label it ff y . Geometrically, the collection of all
solutions to (2.2) is an affine subspace in R n ; MOF selects the element of this subspace
closest to the origin. It is sometimes called a minimum-length solution. There is a matrix
\Phi y , the generalized inverse of \Phi, that calculates the minimum-length solution to a system
of linear equations:
(a) Signal: Hydrogen
(b) Ideal Phase Plane
Time
Frequency
(c) Phase Plane by MOF
Time
Frequency

Figure

2.3: MOF representation is not sparse.
For so-called "Tight Frame" dictionaries MOF is available in closed form. Nice example:
the standard wavelet packet dictionary. One can compute that for all vectors v, k\Phi T vk
In short \Phi
. Notice that \Phi T is simply the analysis
operator.
There are two key problems with the Method of Frames. First, MOF is not sparsity-
preserving. If the underlying object has a very sparse representation in terms of the dictio-
nary, then the coefficients found by MOF are likely to be very much less sparse. Each atom
in the dictionary that has nonzero inner product with the signal is, at least potentially, and
also usually, a member of the solution.

Figure

2.3a shows the signal Hydrogen, made of a single atom in a wavelet packet dic-
tionary. The result of a frame decomposition in that dictionary is depicted in a phase-plane
portrait, Figure 2.3c. While the underlying signal can be synthesized from a single atom,
the frame decomposition involves many atoms, and the phase-plane portrait exaggerates
greatly the intrinsic complexity of the object.
Second, MOF is intrinsically resolution-limited. No object can be reconstructed with
features sharper than those allowed by the underlying operator \Phi y \Phi. Suppose the underlying
object is sharply localized: . The reconstruction will not be ff, but instead
\Phi y \Phiff which, in the overcomplete case, will be spatially spread out. Figure 2.4 presents
a signal TwinSine, consisting of the superposition of two sinusoids that are separated by
less than the so-called Rayleigh Distance 2-=n. We analyze these in a 4-fold overcomplete
discrete cosine dictionary. In this case, reconstruction by MOF, Figure 2.4b, is simply
convolution with the Dirichlet kernel. The result is the synthesis from coefficients with a
broad oscillatory appearance, consisting not of two but of many frequencies, and giving no
visual clue that the object may be synthesized from two frequencies alone.
2.3.2 Matching Pursuit
Mallat and Zhang [23] have discussed a general method for approximate decomposition (1.2)
that addresses the sparsity issue directly. Starting from an initial approximation s
(a) Signal: TwinSine
-0.50.51.5Frequency/Nyquist
(b) MOF coefs
Amplitude
-0.50.51.5Frequency/Nyquist
(c) MP Coefs
Amplitude
-0.50.51.5Frequency/Nyquist
(d) BP Coefs
Amplitude

Figure

2.4: Analyzing TwinSine with a 4-fold overcomplete discrete cosine dictionaries.
and residual R builds up a sequence of sparse approximations stepwise. At stage
k, it identifies the dictionary atom that best correlates with the residual and then adds to
the current approximation a scalar multiple of that atom, so that s
i and R steps, one has a representation
of the form (1.2), with residual R = R (m) . A similar algorithm was proposed for Gabor
dictionaries by S. Qian and D. Chen [30].
An intrinsic feature of the algorithm is that when stopped after a few steps, it yields
an approximation using only a few atoms. When the dictionary is orthogonal, the method
works perfectly. If the object is made up of only m - n atoms and the algorithm is run
for m steps, it recovers the underlying sparse structure exactly.
When the dictionary is not orthogonal, the situation is less clear. Because the algorithm
is myopic, one expects that, in certain cases, it might choose wrongly in the first few
iterations and, in such cases, end up spending most of its time correcting for any mistakes
made in the first few terms. In fact this does seem to happen.
To see this, we consider an attempt at super-resolution. Figure 2.4a portrays again
the signal TwinSine consisting of sinusoids at two closely spaced frequencies. When MP is
applied in this case (Figure 2.4c), using the 4-fold overcomplete discrete cosine dictionary,
the initial frequency selected is in between the two frequencies making up the signal. Because
of this mistake, MP is forced to make a series of alternating corrections that suggest
a highly complex and organized structure. MP misses entirely the doublet structure. One
can certainly say in this case that MP has failed to super-resolve.
Second, one can give examples of dictionaries and signals where MP is arbitrarily sub-optimal
in terms of sparsity. While these are somewhat artificial, they have a character not
so different from the super-resolution example.
DeVore and Temlyakov's Example. Vladimir Temlyakov, in a talk at the IEEE Conference
on Information Theory and Statistics, October 1994, described an example in which
the straightforward greedy algorithm is not sparsity-preserving. In our adaptation of this
example, based on Temlyakov's joint work with R.A. DeVore [10], one constructs a dictionary
having atoms. The first n are the Dirac basis; the final atom involves a linear
combination of the first n with decaying weights. The signal s has an exact decomposition
in terms of A atoms; but the greedy algorithm goes on forever, with an error of size
O(1=
m) after m steps. We illustrate this decay in Figure 2.5a. For this example we set
choose the signal s . The dictionary consists of Dirac
elements
with c chosen to normalize OE n+1 to unit norm.
Shaobing Chen's Example. The DeVore-Temlyakov example applies to the original MP
algorithm as announced by Mallat and Zhang in 1992. A later refinement (see also Pati
[29]) involves an extra step of orthogonalization. One takes all m terms that have entered
at stage m and solves the least squares problem
min
ks
a
for coefficients (a (m)
one forms the residual -
a (m)
, which will be
orthogonal to all terms currently in the model. This method is called Orthogonal Matching
Pursuit (OMP) by Pati [29]. The DeVore-Temlyakov example does not apply to OMP,
but Shaobing Chen found in Summer 1993 an example of similar flavor that does. In this
example, a special signal and dictionary are constructed, with the following flavor. The
dictionary is composed of atoms OE fl with ng. The first A atoms come from the
Dirac dictionary, with . The signal is a simple equiweighted linear
combination of the first A atoms:
Dictionary atoms with fl ? A are a
linear combination of the corresponding Dirac ffi fl and s. OMP chooses all atoms except the
first A before ever choosing one of the first A. As a result, instead of the ideal behavior
one might hope for, terminating after just A steps, one gets n steps before convergence,
and the rate is relatively slow. We illustrate the behavior of the reconstruction error in

Figure

2.5b. We chose 1024. The dictionary was OE
and OE
as
one might have hoped for the ideal behavior
2.3.3 Best Orthogonal Basis
For certain dictionaries, it is possible to develop specific decomposition schemes custom-tailored
to the dictionary.
Wavelet packet and cosine packet dictionaries are examples; they have very special
properties. Certain special subcollections of the elements in these dictionaries amount to
orthogonal bases; one gets in this way a wide range of orthonormal bases (in fact
such orthogonal bases for signals of length n).
Coifman and Wickerhauser [6] have proposed a method of adaptively picking from
among these many bases a single orthogonal basis that is the "best basis". If (s[B] I ) I
denotes the vector of coefficients of s in orthogonal basis B, and if we define the "entropy"
(a) MP on DeVore and Temlyakov's example
m, Number of Terms in Reconstruction
Reconstruction
Greedy: Dashed
(b) OMP on Chen's example
m, Number of Terms in Reconstruction
Reconstruction
Greedy: Dashed

Figure

2.5: Counter examples for MP.
I e(s[B] I ), where e(s) is a scalar function of a scalar argument, they give a fast
algorithm for solving
The algorithm in some cases delivers near-optimal sparsity representations. In par-
ticular, when the object in question has a sparse representation in an orthogonal basis
taken from the library, one expects that BOB will work well. However, when the signal
is composed of a moderate number of highly non-orthogonal components, the method
may not deliver sparse representations - the demand that BOB find an orthogonal basis
prevents it from finding a highly sparse representation. An example comes from the signal
WernerSorrows, which is a superposition of several chirps, sinusoids and Diracs; see

Figure

2.6a. When analyzed with a cosine packet dictionary and the original Coifman-
entropy, BOB finds nothing: it chooses a global sinusoid basis as best; the
lack of time-varying structure in that basis means that all chirp and transient structure in
the signal is missed entirely; see Figure 2.6b.
3 Basis Pursuit
We now discuss our approach to the problem of overcomplete representations. We assume
that the dictionary is overcomplete, so that there are in general many representations
The principle of Basis Pursuit is to find a representation of the signal whose coefficients
have minimal ' 1 norm. Formally, one solves the problem
subject to s: (3.1)
From one point of view, (3.1) is very similar to the Method of Frames (2.3): we are simply
replacing the ' 2 norm in (2.3) with the ' 1 norm. However, this apparently slight change
has major consequences. The Method of Frames leads to a quadratic optimization problem
-226
(a) Signal: Werner Sorrows (b) Phase Plane: BOB by C-W Entropy
Time
Frequency
(c) Phase Plane: BOB by l^1 Entropy
Time
Frequency
(d) Phase Plane: BP
Time
Frequency

Figure

2.6: Analyzing the signal WernerSorrows with a cosine packet dictionary.
with linear equality constraints, and so involves essentially just the solution of a system of
linear equations. In contrast, Basis Pursuit requires the solution of a convex, nonquadratic
optimization problem, which involves considerably more effort and sophistication.
3.1 Linear Programming
To explain the last comment, and the name Basis Pursuit, we develop a connection with
linear programming (LP).
The linear program in so-called standard form [7, 16] is a constrained optimization
problem defined in terms of a variable x
subject to
where c T x is the objective function, is a collection of equality constraints, and
is a set of bounds. The main question is, which variables should be zero.
The Basis Pursuit problem (3.1) can be equivalently reformulated as a linear program
in the standard form (3.2) by making the following translations:
Hence, the solution of (3.1) can be obtained by solving an equivalent linear program. (The
equivalence of minimum ' 1 optimizations with linear programming has been known since
the 1950's; see [2]). The connection between Basis Pursuit and linear programming is useful
in several ways.
3.1.1 Solutions as Bases
In the linear programming problem (3.2), suppose A is an n by m matrix with m ? n, and
suppose an optimal solution exists. It is well know that a solution exists in which at most
n of the entries in the optimal x are nonzero. Moreover, in the generic case, the solution
is so-called nondegenerate, and there are exactly n nonzeros. The nonzero coefficients are
associated with n columns of A, and these columns make up a basis of R n . Once the basis
is identified, the solution is uniquely dictated by the basis. Thus finding a solution to the
LP is identical to finding the optimal basis. In this sense, linear programming is truly a
process of Basis Pursuit.
Translating the LP results into BP terminology, we have the decomposition
The waveforms (OE
are linearly independent but not necessarily orthogonal. The collection
is not, in general, known in advance, but instead depends on the problem data (in this
case s). The selection of waveforms is therefore signal-adaptive.
3.1.2 Algorithms
BP is an optimization principle, not an algorithm. Over the last forty years, a tremendous
amount of work has been done on the solution of linear programs. Until the 1980's, most
work focused on variants of Dantzig's simplex algorithm, which many readers have no doubt
studied. In the last ten years, some spectacular breakthroughs have been made by the use
of so-called "interior-point methods", which use an entirely different principle.
From our point of view, we are free to consider any algorithm from the LP literature as
a candidate for solving the BP optimization problem; both the simplex and interior-point
algorithms offer interesting insights into BP. When it is useful to consider BP in the context
of a particular algorithm, we will indicate this by label: either BP-Simplex or BP-Interior.
BP-Simplex. In standard implementations of the simplex method for LP, one first
finds an initial basis B consisting of n linearly independent columns of A for which the
corresponding solution B \Gamma1 b is feasible (non-negative). Then one iteratively improves
the current basis by, at each step, swapping one term in the basis for one term not in
the basis, using the swap that best improves the objective function. There always exists
a swap that improves or maintains the objective value, except at the optimal solution.
Moreover, LP researchers have shown how one can select terms to swap in such a way as to
guarantee convergence to an optimal solution (anti-cycling rules) [16]. Hence the simplex
algorithm is explicitly a process of "Basis Pursuit": iterative improvement of a basis until
no improvement is possible, at which point the solution is achieved.
Translating this LP algorithm into BP terminology, one starts from any linearly independent
collection of n atoms from the dictionary. One calls this the current decomposition.
Then one iteratively improves the current decomposition by swapping atoms in the current
decomposition for new atoms, with the goal of improving the objective function. By application
of anti-cycling rules, there is a way to select swaps that guarantees convergence to
an optimal solution (assuming exact arithmetic).
BP-Interior. The collection of feasible points fx : is a convex polyhedron
in R m (a "simplex"). The simplex method, viewed geometrically, works by walking
around the boundary of this simplex, jumping from one vertex (extreme point) of the
polyhedron to an adjacent vertex at which the objective is better. Interior-point methods
instead start from a point x (0) well inside the interior of the simplex (x (0) AE 0) and go
"through the interior" of the simplex. Since the solution of a LP is always at an extreme
(a) Signal: Carbon
(b) Phase Plane: MOF
Time
Frequency
Phase Plane: BOB
Time
Frequency
Phase Plane: MP
Time
Frequency
Phase Plane: BP
Time
Frequency
Figure

3.1: Analyzing the signal Carbon with a wavelet packet dictionary.
point of the simplex, as the interior-point method converges, the current iterate x (k) approaches
the boundary. One may abandon the basic interior-point iteration and invoke a
"crossover" procedure that uses simplex iterations to find the optimizing extreme point.
Translating this LP algorithm into BP terminology, one starts from a solution to the
overcomplete representation problem \Phia iteratively modifies the
coefficients, maintaining feasibility \Phia applying a transformation that effectively
sparsifies the vector a (k) . At some iteration, the vector has - n significantly nonzero entries,
and it "becomes clear" that those correspond to the atoms appearing in the final solution.
One forces all the other coefficients to zero and "jumps" to the decomposition in terms of
the - n selected atoms. (More general interior-point algorithms start with a (0) ? 0 but
don't require the feasibility \Phia
3.2 Examples
We now give computational examples of BP in action.
3.2.1 Carbon
The synthetic signal Carbon is a composite of 6 atoms: a Dirac, a sinusoid, and 4 mutually
orthogonal wavelet packet atoms, adjacent in the time-frequency plane. The wavelet packet
dictionary of depth employed, based on filters for Symmlets with 8 vanishing
moments. (Information about problem sizes for all examples is given in Table 1).

Figure

3.1 displays the results in phase-plane form; for comparison, we include the phase
planes obtained using MOF, MP, and BOB. First, note that MOF uses all basis functions
that are not orthogonal to the 6 atoms, i.e. all the atoms at times and frequencies that
overlap with some atom appearing in the signal. The corresponding phase plane is very
diffuse or smeared out. Second, MP is able to do a relatively good job on the sinusoid and
the Dirac, but makes mistakes in handling the 4 close atoms. Third, BOB cannot handle
the nonorthogonality between the Dirac and the cosine; it gives a distortion (a coarsening)
(a) Signal: FM
Time
Frequency
(f) PhasePlane: BP
Time
Frequency
Time
Frequency
Time
Frequency
Time
Frequency
Figure

3.2: Analyzing the signal FM-Cosine with a cosine packet dictionary.
of the underlying phase plane picture. Finally, BP finds the "exact" decomposition in the
sense that the four atoms in the quad, the Dirac and the sinusoid are all correctly identified.
3.2.2 TwinSine
Recall that the signal TwinSine in Figure 2.4a consists of 2 cosines with frequencies closer
together than the Rayleigh distance. In Figure 2.4d, we analyze these in the 4-fold over-complete
discrete cosine dictionary. Recall that in this example, MP began by choosing at
the first step a frequency in between the two ideal ones and then never corrected the error.
In contrast, BP resolves the two frequencies correctly.
3.2.3 FM Signal

Figure

3.2a displays the artificial signal FM-Cosine consisting of a frequency-modulated
sinusoid superposed with a pure sinusoid: Figure
3.2b shows the ideal phase plane.
In

Figure

3.2c-f we analyze it using the cosine packet dictionary based on a primitive bell
of width 16 samples. It is evident that BOB cannot resolve the nonorthogonality between
the sinusoid and the FM signal. Neither can MP. However, BP yields a clean representation
of the two structures.
3.2.4 Gong

Figure

3.3a displays the Gong signal, which vanishes until time t 0 and then follows a decaying
sinusoid for t ? t 0 .
In

Figures

3.3b-3.3d, we analyze it with the cosine packet dictionary based on a primitive
bell of width samples. BP gives the finest representation of the decay structure; visually
somewhat more interpretable than the BOB and MP results.
0(a) Signal: Gong
(c) Phase Plane: MOF
Time
Frequency
Time
Frequency
Time
Frequency
Time
Frequency
Figure

3.3: Analyzing the signal Gong with a cosine packet dictionary.
3.3 Comparisons
We briefly compare BP with the three main methods introduced in Section 2.3.
3.3.1 Matching Pursuit
At first glance MP and BP seem quite different. MP is an iterative algorithm, which does
not explicitly seek any overall goal, but merely applies a simple rule repeatedly. In contrast,
BP is a principle of global optimization without any specified algorithm. The contrast of
Orthogonal MP with a specific algorithm, BP-Simplex, may be instructive. Orthogonal
Matching Pursuit starts from an "empty model" and builds up a signal model an atom at
a time, at each step adding to the model only the most important new atom among all
those not so far in the model. In contrast, BP-Simplex, starts from a "full" model (i.e.
representation of the object in a basis) and then iteratively improves the "full" model, by
taking relatively useless terms out of the model, swapping them for useful new ones. Hence,
MP is a sort of build-up approach, while BP-Simplex is a sort of swap-down approach.
To make BP and BOB most comparable, suppose that they are both working with a
cosine packet dictionary, and note that the ' 1 -norm of coefficients is what Coifman and
[6] call an "additive measure of information". So suppose we apply the
. Then the two methods
compare as follows: in BOB, we are optimizing E only over orthogonal bases taken from the
dictionary, while in BP we are optimizing E over all bases formed from the dictionary.
This last remark suggests that it might be interesting to apply the BOB procedure with
the ' 1 norm as entropy in place of the standard Coifman-Wickerhauser entropy. In Figure
2.6c we try this on the WernerSorrows example of Section 2.3.3. The signal is analyzed in
a cosine packet dictionary, with primitive bell width 16. The ' 1 entropy results in a time-varying
basis that reveals clearly some of the underlying signal structure. The ' 1 entropy
Phase Plane: BP Iteration
Time
Frequency
Time
Frequency
Time
Frequency
Time
Frequency
Time
Frequency
Time
Frequency
Figure

3.4: Phase plane evolution at BP-Interior iteration.
by itself improves the performance of BOB; but BP does better still (Figure 2.6d).
This connection between BP and BOB suggests an interesting algorithmic idea. In the
standard implementation of the simplex method for LP, one starts from an initial basis and
then iteratively improves the basis by swapping one term in the basis for one term not in
the basis, using the swap that best improves the objective function. Which initial basis? It
seems natural in BP-Simplex to use the Coifman-Wickerhauser algorithm and employ as a
start the best orthogonal basis.
With this choice of starting basis, BP can be seen as a method of refining BOB by
swapping in non-orthogonal atoms in place of orthogonal ones whenever this will improve
the objective.
3.3.3 Method of Frames
As already discussed, MOF and BP differ in the replacement of an l 2 objective function
by an l 1 objective. BP-Interior has an interesting relation to the Method of Frames. BP-
Interior initializes with the Method of Frames solution. Hence one can say that BP sequentially
"improves" on the Method of Frames. Figure 3.4 shows a "movie" of BP-Interior
in action on the FM-Cosine example, using a cosine packet dictionary. Six stages in the
evolution of the phase plane are shown, and one can see how the phase plane improves in
clarity, step-by-step.
Variations
The recent development of time-frequency dictionaries motivates most of what we have
done so far. However, the methods we have developed are general and can be applied to
other dictionaries, with interesting results.
4.1 Stationary Smooth Wavelets
The usual (orthonormal) dictionaries of (periodized) smooth wavelets consist of wavelets at
scales indexed by at the j-th scale, there are 2 j wavelets of width
. The wavelets at this scale are all circulant shifts of each other, the shift being n=2 j
samples. Some authors [32] have suggested that this scheme can be less than satisfactory,
essentially because the shift between adjacent wavelets is too large. They would say that
if the important "features" of the signal are (fortuitously) "aligned with" the wavelets
in the dictionary, then the dictionary will provide a sparse representation of the signal;
however, because there are so few wavelets at level j, then most likely, the wavelets in
the dictionary are not "precisely aligned" with features of interest, and the dictionary may
therefore provide a very diffuse representation.
The stationary wavelet dictionary has, at the j-th level, n (not are
all the circulant shifts of the basic wavelet of width - n=2 j . Since this dictionary always
contains wavelets "aligned with" any given feature, the hope is that such a dictionary
provides a superior representation.
Panel 4.1a shows the signal HeaviSine, and 4.1b shows the result of BP with the
Stationary Symmlet-8 dictionary mentioned in Section 2.1; the coefficients are displayed
in a multi-resolution fashion, where at level j all the coefficients of scale 2 j =n are plotted
according to spatial position.
There is a surprisingly close agreement of the BP representation in a stationary wavelet
dictionary with ideas about signal representation associated with the "Multi-Scale Edges"
ideas of Mallat and Hwang [22]. The Multi-Scale Edge method analyzes the continuous
wavelet transform (CWT) at scale 2 \Gammaj and identifies the maxima of this transform. Then
it selects maxima that are "important" by thresholding based on amplitude. These "im-
portant" maxima identify important features of the signal. Mallat and Hwang proposed
an iterative method that reconstructs an object having the same values of the CWT at
"maxima". This is almost (but not quite) the same thing as saying that one is identifying
"important" wavelets located at the corresponding maxima, and reconstructing the object
using just those maxima.
Panel 4.1c shows a CWT of HeaviSine based on the same Symmlet-8 wavelet, again in
multi-resolution fashion; Panel 4.1d shows the maxima of the CWT. At fine scales, there is
virtually a 1-1 relationship between the maxima of the transform and the wavelets selected
by BP; compare panel 4.1b. So in a stationary wavelet dictionary, the global optimization
principle BP yields results that are close to certain heuristic methods.
An important contrast: Meyer has a counterexample to multi-scale edge approaches,
showing that the Mallat-Hwang approach may fail in certain cases [26]; but there can be
no such counterexamples to BP.
4.2 Dictionary Mergers
An important methodological tool is the ability to combine dictionaries to make bigger,
more expressive dictionaries. We mention here two possibilities. Examples of such decompositions
are given in Section 5 below.
Jump+Sine. Merge the Heaviside dictionary with a Fourier dictionary. Either dictionary
can efficiently represent objects that the other cannot; for example, Heavisides have
difficulty representing sinusoids, while sinusoids have difficulty representing jumps. Their
combination might therefore be able to offer the advantages of both.
(a) Signal: HeaviSine
-22
Position
log(resolution)
(b) Coefs from BP on HeaviSine
-22
Position
log(resolution)
(c) Coefs from CWT on HeaviSine
-22
Position
log(resolution)
(d) Mutiscale Edges Representation of HeaviSine

Figure

4.1: Analyzing the signal HeaviSine with a stationary wavelet dictionary.
Jump+Wavelet. For similar reasons, one might want to merge Heavisides with Wavelets.
In fact, we have found it sometimes preferable instead to merge "tapered heavisides" with
wavelets; these are step discontinuities that start at 0, jump at time t 0 to a level one unit
higher, and later decay to the original 0 level.
De-Noising
We now adapt BP to the case of noisy data. We assume data of the form
where (z i ) is a standard white Gaussian noise, oe ? 0 is a noise level, and s is the clean
signal. In this setting, s is unknown, while y is known. We don't want to get an exact
decomposition of y, so we don't apply BP directly. Instead decompositions like (1.2) become
relevant.
5.1 Proposal
Basis Pursuit De-Noising (BPDN) refers to solution of
min
The solution a (-) is a function of the parameter -. It yields a decomposition into
signal-plus-residual:
where s . The size of the residual is controlled by -. As - ! 0, the residual goes
to zero and the solution behaves exactly like BP applied to y. As - !1, the residual gets
large; we have r (-) ! y and s (-) ! 0.
Recently Michael Saunders and Shaobing Chen have shown that (5.1) is equivalent to
the following perturbed linear program:
subject to Ax
Perturbed linear programming is really quadratic pro-
gramming, but retains structure similar to linear programming. Hence we can have a similar
classification of algorithms, into BPDN-Simplex and BPDN-Interior-Point types. (In
quadratic programming, "simplex like" algorithms are usually called Active Set algorithms,
so our label is admittedly nonstandard.)
5.2 Choice of -
Assuming the dictionary is normalized so that kOE fl, we set - to the value
where p is the cardinality of the dictionary.
This can be motivated as follows. In the case of a dictionary that is an orthonormal basis,
a number of papers [11, 14] have carefully studied an approach to de-noising by so-called
"soft-thresholding in an orthonormal basis". In detail, suppose that \Phi is an orthogonal
matrix, and define empirical OE-coefficients by
Define the soft threshold nonlinearity j - and define the thresholded
empirical coefficients by
This is soft thresholding of empirical orthogonal coefficients. The papers just cited show
that thresholding at - n has a number of optimal and near-optimal properties as regards
mean-squared error.
We claim that (again in the case of an ortho-basis) the thresholding estimate -
ff is also
the solution of (5.1). Observe that the soft thresholding nonlinearity solves the scalar
minimum problem:
Note that, because of the orthogonality of \Phi, so we can rewrite
(5.1) in this case as
min
Now applying (5.2) coordinatewise establishes the claim.
The scheme we have suggested here - to be applied in overcomplete as well as orthogonal
settings - therefore includes soft-thresholding in ortho-bases as a special case. Formal
arguments similar to those in [13] can be used to give a proof that mean-squared error
properties of the resulting procedure are near-optimal under certain conditions.
(d) Recovered: BOB
(f) Recovered: BP
(c) Recovered: MOF
(a) Signal: Gong
(b) The Noised:

Figure

5.1: De-Noising noisy Gong with a cosine packet dictionary.
5.3 Examples
We present two examples of BPDN in action with time-frequency dictionaries. We compare
BPDN with three other de-noising methods adapted from MOF, MP and BOB. Method-
of-Frames De-Noising (MOFDN) refers to minimizing the least square fit error plus a l 2
penalizing term:
min a
ks
2where - is a penalizing parameter; we chose - in these examples to be oe
log(p). Matching
Pursuit De-Noising (MPDN) runs Matching Pursuit until the coefficient associated with
the selected atom gets below the threshold
oe. The Best Orthogonal Basis De-Noising
(BOBDN) is a thresholding scheme in the best orthogonal basis chosen by the BOB
algorithm with a special entropy [12].
5.3.1 Gong

Figure

5.1 displays de-noising results on the signal Gong, at signal to noise ratio 1, using a
cosine packet dictionary. Panel a) displays the noiseless signal and panel b) displays a noisy
version. Panels c)-f) display de-noising results for MOF, BOB, MP, and BP, respectively.
BP outperforms the other methods visually.
5.3.2 TwinSine

Figure

5.2 employs the signal TwinSine, described earlier, to investigate super-resolution
in the noisy case. Panels a) and b) give the noiseless and noisy TwinSine, respectively.
Using a 4-fold overcomplete discrete cosine dictionary, reconstructions by the MOF, MP,
and by BPDN are given. MOF gives a reconstruction that is inherently resolution-limited
and oscillatory. As in the noiseless case, MP gives a reconstruction that goes wrong at
selects the average of the two frequencies in the TwinSine signal. BP correctly
resolves the non-negative doublet structure.
(a) TwinSine
(b) Noised TwinSine,
(c) DCT transform
Frequency/Nyquist
(d) MOF Coefs
Frequency/Nyquist
Frequency/Nyquist
(f) BP Coefs
Frequency/Nyquist

Figure

5.2: De-Noising noisy TwinSine-2 with a 4-fold overcomplete discrete cosine dictionary

5.4 Total Variation De-Noising
Recently, Rudin, Osher and Fatemi [28] have called attention to the possibility of de-noising
images using total-variation penalized least-squares. More specifically, they propose the
optimization problem
min g2
where TV (g) is a discrete measure of the total variation of g. A solution of this problem is
the de-noised object. Li and Santosa [20] have developed an alternative algorithm for this
problem based on interior-point methods for convex optimization.
For the 1-dimensional case (signals rather than images) it is possible to implement
what amounts to total variation de-noising by applying BPDN with a Heaviside dictionary.
Indeed, if s is an arbitrary object, it has a unique decomposition in Heavisides (recall
(2.1)). Suppose that the object is 0 at and that the decomposition is
; then the total variation is given by
Moreover to get approximate equality even for objects not obeying zero-boundary condi-
tions, one has only to normalize OE 0 appropriately. Consequently, total variation de-noising
is essentially a special instance of our proposal (5.1).
We have studied BPDN in the Heaviside dictionary, thereby obtaining essentially a
series of tests of TV De-Noising. For comparison, we considered also soft thresholding
in orthogonal wavelet dictionaries based on the S8-Symmlet smooth wavelet. We also
constructed a new dictionary, based on the Jump+Wave merger of S8-Symmlet wavelets
with "Smoothly Tapered Heavisides", which is to say, atoms OE fl that jump at a given point fl
and then decay smoothly away from the discontinuity. For comparability with the Heaviside
dictionary, we normalized the Jump+Wave dictionary so that every kOE fl k TV - 1.
(a) Signal: Blocks
(d) BPDeNoise: Heaviside
(f) BPDeNoise: Jump+Wave
(c) Sorted Coefs
Order
Amplitude
Dotted: Heaviside
Wave
Jump+Wave
bpfig54.m 16-May-95 Figure 5.4: TV DeNoise
Signal

Figure

5.3: De-Noising noisy Blocks.
A typical result, for the object Blocky, is presented in Figure 5.3. From the point
of view of visual appearance, total variation reconstruction (panel d) far outperforms the
other methods.
Of course, the object Blocky has a very sparse representation in terms of Heavisides.
When we consider an object like Cusp, which is piecewise smooth rather than piecewise
constant, the object will no longer have a sparse representation. On the other hand, using
the Jump+Wave dictionary based on a merger of wavelets with tapered Heavisides will lead
to a sparse representation - see Figure 5.4c. One can predict that a Heaviside dictionary
will perform less well than this merged dictionary.
This completely obvious comment, translated into a statement about total variation de-
noising, becomes a surprising prediction. One expects that the lack of sparse representation
of smooth objects in the Heaviside dictionary will translate into worse performance of TV
de-noising than of BPDN in the merged Jump+Wave dictionary.
To test this, we conducted experiments. Figure 5.4 compares
de-noising, and BPDN in the merged Jump+Wave dictionary. TV De-Noising now exhibits
visually distracting stairstep artifacts; the dictionary Jump+Wave seems to us to behave
much better.
6 Solutions of Large-Scale Linear Programs
As indicated in Section 3.1, the optimization problem (3.1) is equivalent to a linear program
(3.2). Also, as in Section 5.1, the optimization problem (5.1) is equivalent to a perturbed
linear program (5.3). The problems in question are large-scale; we have conducted decompositions
of signals of length in a wavelet packet dictionary, leading to a linear
program of size 8192 by 212; 992.
Over the last ten years there has been a rapid expansion in the size of linear programs
that have been successfully solved using digital computers. A good overview of the recent
rapid progress in this field and the current state of the art is afforded by the article of Lustig,
(a) Signal: Cusp
(c) Sorted Coefs
Order
Amplitude
Dotted: Heaviside
Wave
Jump+Wave
bpfig56.m 16-May-95 Figure 5.6: Dictionary Merge
Signal

Figure

5.4: De-Noising noisy Cusp.
Marsten and Shanno [21] and the accompanying discussions by Bixby [1], Saunders [31],
Todd [33], and Vanderbei [34]. Much of the rapid expansion in the size of linear programs
solved is due to the "Interior Point revolution" initiated by Karmarkar's proof that a pseudo-polynomial
time algorithm could be based on an interior-point method [18]. Since then a
very wide array of interior-point algorithms have been proposed and considerable practical
[21] and theoretical [27] understanding is now available. In this section we describe our
algorithm and our experience with it.
6.1 Duality Theory
We consider the linear program in the standard form
subject to
This is often called the primal linear program. The primal linear program is equivalent to
the dual linear program
subject to A T y
x is called the primal variable; y and z are called the dual variables. The term "primal
refers to the quantity kb \Gamma Axk 2 ; the term "dual infeasibility" refers to kc \Gamma
the term "duality gap" refers to the difference between the primal objective and
the dual objective: c T y.
A fundamental theorem of linear programming states that (x; solves the linear
program (6.1) if and only if the primal infeasibility, the dual infeasibility and the duality
gap are all zero. Therefore when (x; are nearly primal feasible and nearly dual feasible,
the duality gap offers a good description about the accuracy of (x; as a solution: the
smaller the duality gap is, the closer (x; are to the optimal solution.
6.2 A Primal-Dual Log-Barrier LP Algorithm
Mathematical work on interior-point methods over the last ten years has led to a large
variety of approaches, with names like projective scaling, (primal/dual) affine scaling, (pri-
mal/dual) logarithmic barrier and predictor-corrector. We cannot summarize all these ideas
here; many of them are mentioned in [21] and others are covered in the references of that
article.
Our approach is based on a primal-dual log-barrier algorithm. In order to regularize
standard LP, Gill et al. [15] proposed solving the following perturbed LP:
subject to Ax
where fl and ffi are normally small (e.g. regularization parameters. (We comment
that such a perturbed LP with solves the BPDN problem (5.1)). The main steps of
the interior point algorithm are as follows:
1. Set parameters: the feasibility tolerance FeaTol, the duality gap tolerance PDGapTol,
the two regularization parameters fl and ffi .
2. Initialize x ? 0;
3. Loop
(a) Set
where X and Z are diagonal matrices composed from x and z.
(b) Solve
for 4y and set
(c) Calculate the primal and dual step sizes ae p ; ae d and update the variables:
(d) Increase k by 1.
4. Until the following three conditions are satisfied:
(a) Primal Infeasibility
FeaTol.
(b) Dual
FeaTol.
(c) Duality
For fuller discussions of this and related algorithms, again see [15] or references there.
While in principle we could have based our approach on other interior-point schemes,
the primal-dual approach naturally incorporates several features we found useful. First, the
iterates z do not have to be feasible. We are only able to choose a starting point that
is nearly feasible and remain nearly feasible throughout the sequence of iterations. Second,
after both primal and dual feasibility have been nearly achieved, it is easy to check for
closeness to the solution value; at the limiting solution c T x  and the duality gap
quantifies the distance from this ideal.
6.3 Implementation Heuristics
The primal-dual log barrier algorithm we just described works in a fashion similar to other
interior-point methods [21]. It starts from an initial feasible (or nearly feasible) solution
located at or near the "center" of the feasible region, and iteratively improves the current
solution until the iterates (x; achieve the desired accuracy. It requires a relatively
small number of iterations: for example, a few dozen iterations would be common. Each
iteration requires the solution of a system of equations involving A, A T , and other problem
data like x; z. In the primal-dual log barrier method, the system is (6.4). Thus the
numerical solution to a linear program by interior-point methods amounts to a sequence
of several dozen solutions of special systems of linear equations. This leads to a slogan: if
those systems can be solved rapidly, then it is possible to solve the LP rapidly.
Of course, in general solving systems of equations is not rapid: a general n by n system
takes order O(n 3 ) time to solve by standard elimination methods or by modern
stable factorization schemes [17, 16]. In order for practical algorithms to be based on the
interior-point heuristic, it is necessary to be able to solve the systems of equations much
more rapidly than one could solve general systems. In the current state of the art of linear
programming [31], one attempts to do this by exploiting sparsity of the underlying matrix
A.
However, the optimization problems we are interested in have a key difference from the
successful large-scale applications outlined in [21]. The matrix A we deal with is not at
all sparse; it is generally completely dense. For example, if A is generated from a Fourier
dictionary, most of the elements of A will be of the same order of magnitude. Because of
this density, it is unlikely that existing large-scale interior-point computer codes could be
easily applied to the problems described in this paper.
In our application we have a substitute for sparsity. We consider only dictionaries that
have fast implicit algorithms for \Phia and \Phi T s, and therefore lead to linear programs where
the A matrix admits fast implicit algorithms for both Au and A T v. Now whenever one has
fast implicit algorithms, it is natural to think of solving equations by conjugate-gradient
methods; such methods allow one to solve equations using only products Bv
with various strategically chosen vectors v. Adapting such ideas, one develops fast implicit
algorithms for (ADA and attempts to solve the central equations (6.4) iteratively,
avoiding the costly step of explicitly forming the matrices (ADA
In our application, we do not really need an exact solution of the optimization problem.
Moreover, we have a natural initial solution - from MOF - that would be viewed by some
researchers as already an acceptable method of atomic decomposition. By starting from
this decomposition and applying a strategy based on a limited number of iterations of our
algorithm, we get what we view as an iterative improvement on MOF. Compare Figure 3.4.
CPU Running Time in Seconds

Figure

Signal Problem Size
MOF BOB MP BP

Figure

2.4 TwinSine 256 .3500 - .6667 7.517

Figure

2.6 WernerSorrows 1024 - .9500 - 158.2

Figure

3.1 Carbon 1024 .2000 2.617 2.650 11.70

Figure

3.2 FM-Cosine 1024 1.050 .9333 182.9 150.2

Figure

3.3 Gong 1024 1.433 5.683 50.63 448.2

Figure

4.1 HeaviSine 256 - 26.92

Figure

5.1 Noisy Gong 1024 2.117 6.767 8.600 142.2

Figure

5.2 Noisy TwinSine 256 .4167 - .6833 5.717

Table

1: CPU Running Times of the Examples
We stress that our strategy is to "pursue an optimal basis"; while we would like to reach
the optimal basis, we make no specific claims that we can always reach it in reasonable
time; perhaps the "pursuit" language will help remind one of this fact. We do believe that
the pursuit process, carried out for whatever length of time we are willing to invest in it,
makes a useful improvement over the Method of Frames.
6.4 Routine Settings For BP
Our strategy for routine signal processing by BP is as follows:
ffl We employ the "primal-dual logarithmic barrier method" for perturbed LP [15].
ffl We suppose fast implicit algorithms for Aa and A T s.
ffl We only aim to reach an approximate optimum.
would usually suffice for this.
ffl Each barrier iteration involves approximate solution of the central equations (6.4)
using the conjugate-gradient method, e. g. at accuracy
We refer the reader to [4] for more detailed discussion of the our implementation.
6.5 Complexity Analysis

Table

1 displays the CPU times in seconds spent in running various atomic decomposition
techniques in our experiments; all computation was done on a Sun Sparc20 workstation.
We employ a conjugate-gradient solver for the generalized inverse in the MOF solution
(2.4); the resulting algorithm for MOF has a complexity order O(n log(n)). We implement
Coifman and Wickerhauser's BOB algorithm [6], which also has a complexity of order
O(n log(n)). We observe that BP is typically slower than MOF and BOB. BP is also slower
than MP (which has a quasi-linear complexity, depending on the number of chosen atoms)
except on the FM-Cosine signal in Figure 3.2.
Several factors influence the running time of Basis Pursuit:
1. Problem Sizes. The complexity goes up "quasi-linearly" as the problem size increases
[4].
2. Parameter Settings. The complexity of our primal-dual logarithmic barrier interior-point
implementation depends on both the the accuracy of the solution and the accuracy
of the conjugate-gradient solver. The accuracy of the solution is determined by the two
parameters FeaTol, PDTol controlling the number of barrier iterations, and the parameter
CGAccuracy, which decides the accuracy of the CG solver and consequently the number of
CG iterations. As the required solution accuracy goes up, the complexity goes up drasti-
cally. We recommend setting FeaTol, PDGapTol and CGAccuracy at 10 \Gamma1 for routine signal
processing; we recommend 10 \Gamma2 or 10 \Gamma3 when one is interested in superresolution. We used
the setting 10 \Gamma1 for the computational experiments presented in Figures 2.6, 3.1, 3.2, 3.3,
5.1 and 5.1. In

Figures

2.5, 3.2 and 5.1, we attempted to super-resolve two cosines with
close frequencies; thus we use the setting 10 \Gamma2 . In Figure 4.1, we used the setting 10 \Gamma3 .
3. Signal Complexity. When the signal has a very sparse representation, the algorithm
converges quickly. The signal Carbon, which contains only 6 atoms from a wavelet packet
dictionary, takes about 10 seconds, whereas it takes about 7 minutes for the signal Gong,
which is much more complex.
4. Basis Pursuit versus Basis Pursuit De-Noising. We employ the same interior-point
implementation for BP and BPDN, except for a difference in the value of the regularization
parameter ffi: ffi is small, e.g. 10 \Gamma4 for BP, while BPDN. The choice
it regularizes the central equations to be solved at each barrier iteration. Thus the BPDN
implementation seems to converge more quickly than the BP implementation. For example,
according to our experiments [4], it takes only 3 minutes to perform BPDN on the noisy
Gong signal of length 1024 with a cosine packet dictionary at the parameter setting
it takes about 8 hours to perform BP on the signal Gong at the same parameter setting.
6.6 Reproducible Research
This paper has been written following the discipline of Reproducible Research described in
[3]. As a complement to this article, we are releasing the underlying software environment
by placing it on internet for access either by anonymous FTP or WWW browsers.
Web Browser: http://playfair.stanford.edu/~schen/Atomizer.html
FTP Client: playfair.stanford.edu file: pub/chen-s/Atomizer0600.tar.Z
For reasons of space we refer the reader to [4] for a discussion of related work in statistics
and elsewhere.



--R

Progress in linear programming.
Least Absolute Deviations: Theory
WaveLab and reproducible research.

Remarques sur l'analyze de Fourier

Linear Programming and Extensions.

Ten Lectures on Wavelets.
Some remarks on greedy algorithms.

de-noising in an orthonormal basis chosen from a library of bases
Empirical atomic decomposition.
Wavelet shrinkage: asymptopia?
Solving reduced KKT systems in barrier methods for linear and quadratic programming.
Numerical Linear Algebra and Optimization.
Matrix Computations
A new polynomial-time algorithm for linear programming
A primal-dual interior point algorithm for linear programming
An affine scaling algorithm for minimizing total variation in image enhancement.
Interior point methods for linear program- ming: computational state of the art
detection and processing with wavelets.
Matching Pursuit in a time-frequency dictionary
On finding primal- and dual- optimal bases
Ondelettes sur l'intervalle.
Wavelets: Algorithms and Applications.

Nonlinear total-variation-based noise removal algorithms
Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition.
Signal representation using adaptive normalized Gaussian functions.
Major Cholesky would feel proud.
Shiftable multiscale transforms.
Theory and practice for interior point methods.
Interior point methods: algorithms and formulations.
--TR

--CTR
Anders la Cour-Harbo, Fast estimation of optimal sparseness of music signals, Proceedings of the 24th IASTED international conference on Signal processing, pattern recognition, and applications, p.205-209, February 15-17, 2006, Innsbruck, Austria
M. E. Davies , L. Daudet, Sparse audio representations using the MCLT, Signal Processing, v.86 n.3, p.457-470, March 2006
Yuanqing Li , Andrzej Cichocki , Shun-ichi Amari, Analysis of sparse representation and blind source separation, Neural Computation, v.16 n.6, p.1193-1234, June 2004
Zoltn Szab , Andrs Lrincz, -Sparse representations: generalized sparse approximation and the equivalent family of SVM tasks, Acta Cybernetica, v.17 n.3, p.605-614, January 2006
Tong Zhang, Approximation bounds for some sparse kernel regression algorithms, Neural Computation, v.14 n.12, p.3013-3042, December 2002
Prasanth B. Nair , Arindam Choudhury , Andy J. Keane, Some greedy learning algorithms for sparse regression and classification with mercer kernels, The Journal of Machine Learning Research, 3, 3/1/2003
Shachar Fleishman , Iddo Drori , Daniel Cohen-Or, Bilateral mesh denoising, ACM Transactions on Graphics (TOG), v.22 n.3, July
Peng Xu , Dezhong Yao, Two dictionaries matching pursuit for sparse decomposition of signals, Signal Processing, v.86 n.11, p.3472-3480, November 2006
Joel A. Tropp, Algorithms for simultaneous sparse approximation: part II: Convex relaxation, Signal Processing, v.86 n.3, p.589-602, March 2006
Yaakov Tsaig , David L. Donoho, Breakdown of equivalence between the minimal l1-norm solution and the sparsest solution, Signal Processing, v.86 n.3, p.533-548, March 2006
Chiranjib Bhattacharyya, Second Order Cone Programming Formulations for Feature Selection, The Journal of Machine Learning Research, 5, p.1417-1433, 12/1/2004
Gradient LASSO for feature selection, Proceedings of the twenty-first international conference on Machine learning, p.60, July 04-08, 2004, Banff, Alberta, Canada
Roger Koenker, Quantile regression for longitudinal data, Journal of Multivariate Analysis, v.91 n.1, p.74-89, October 2004
Balaji Krishnapuram , Lawrence Carin , Mario A. T. Figueiredo , Alexander J. Hartemink, Sparse Multinomial Logistic Regression: Fast Algorithms and Generalization Bounds, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.27 n.6, p.957-968, June 2005
Bob L. Sturm , Laurent Daudet , Curtis Roads, Pitch-shifting audio signals using sparse atomic approximations, Proceedings of the 1st ACM workshop on Audio and music computing multimedia, October 27-27, 2006, Santa Barbara, California, USA
Pando Georgiev , Panos Pardalos , Fabian Theis, A bilinear algorithm for sparse representations, Computational Optimization and Applications, v.38 n.2, p.249-259, November  2007
F. Malgouyres, Image Compression Through a Projection onto a Polyhedral Set, Journal of Mathematical Imaging and Vision, v.27 n.2, p.193-200, February  2007
Tong Zhang, On the Dual Formulation of Regularized Linear Systems with Convex Risks, Machine Learning, v.46 n.1-3, p.91-129, 2002
Sayan Mukherjee , Qiang Wu, Estimation of Gradients and Coordinate Covariation in Classification, The Journal of Machine Learning Research, 7, p.2481-2514, 12/1/2006
Yutaka Ohtake , Alexander Belyaev , Hans-Peter Seidel, Sparse surface reconstruction with adaptive partition of unity and radial basis functions, Graphical Models, v.68 n.1, p.15-24, January 2006
Shihao Ji , Lawrence Carin, Bayesian compressive sensing and projection optimization, Proceedings of the 24th international conference on Machine learning, p.377-384, June 20-24, 2007, Corvalis, Oregon
Anna C. Gilbert , S. Muthukrishnan , Martin J. Strauss, Approximation of functions over redundant dictionaries using coherence, Proceedings of the fourteenth annual ACM-SIAM symposium on Discrete algorithms, January 12-14, 2003, Baltimore, Maryland
Lorenzo Peotta , Lorenzo Granai , Pierre Vandergheynst, Image compression using an edge adapted redundant dictionary and wavelets, Signal Processing, v.86 n.3, p.444-456, March 2006
Cheng-En Guo , Song-Chun Zhu , Ying Nian Wu, Modeling Visual Patterns by Integrating Descriptive and Generative Methods, International Journal of Computer Vision, v.53 n.1, p.5-29, June
Arthur E. C. Pece, On the computational rationale for generative models, Computer Vision and Image Understanding, v.106 n.1, p.130-143, April, 2007
Sayan Mukherjee , Ding-Xuan Zhou, Learning Coordinate Covariances via Gradients, The Journal of Machine Learning Research, 7, p.519-549, 12/1/2006
Yaakov Tsaig , David L. Donoho, Extensions of compressed sensing, Signal Processing, v.86 n.3, p.549-571, March 2006
Fabian J. Theis , Gonzalo A. Garca, On the use of sparse signal decomposition in the analysis of multi-channel surface electromyograms, Signal Processing, v.86 n.3, p.603-623, March 2006
R. Gribonval , R. M. Figueras i Ventura , P. Vandergheynst, A simple test to check the optimality of a sparse signal approximation, Signal Processing, v.86 n.3, p.496-510, March 2006
Alexander J. Smola , Bernhard Schlkopf, Bayesian kernel methods, Advanced lectures on machine learning, Springer-Verlag New York, Inc., New York, NY,
Fabian J. Theis , Pando Georgiev , Andrzej Cichocki, Robust sparse component analysis based on a generalized Hough transform, EURASIP Journal on Applied Signal Processing, v.2007 n.1, p.86-86, 1 January 2007
Charles A. Micchelli , Massimiliano Pontil, Feature space perspectives for learning the kernel, Machine Learning, v.66 n.2-3, p.297-319, March     2007
Sylvain Fischer , Rafael Redondo , Laurent Perrinet , Gabriel Cristbal, Sparse approximation of images inspired from the functional architecture of the primary visual areas, EURASIP Journal on Applied Signal Processing, v.2007 n.1, p.122-122, 1 January 2007
L. Daudet , B. Torrsani, Hybrid representations for audiophonic signal encoding, Signal Processing, v.82 n.11, p.1595-1617, November 2002
F. Malgouyres, Rank related properties for Basis Pursuit and total variation regularization, Signal Processing, v.87 n.11, p.2695-2707, November, 2007
Alexander M. Bronstein , Michael M. Bronstein , Michael Zibulevsky, Blind source separation using block-coordinate relative Newton method, Signal Processing, v.84 n.8, p.1447-1459, August 2004
Yee Whye Teh , Max Welling , Simon Osindero , Geoffrey E. Hinton, Energy-based models for sparse overcomplete representations, The Journal of Machine Learning Research, 4, 12/1/2003
Yee Whye Teh , Max Welling , Simon Osindero , Geoffrey E. Hinton, Energy-based models for sparse overcomplete representations, The Journal of Machine Learning Research, v.4 n.7-8, p.1235-1260, October 1 - November 15, 2004
Cynthia Dwork , Frank McSherry , Kunal Talwar, The price of privacy and the limits of LP decoding, Proceedings of the thirty-ninth annual ACM symposium on Theory of computing, June 11-13, 2007, San Diego, California, USA
Masashi Sugiyama , Klaus-Robert Mller, The subspace information criterion for infinite dimensional hypothesis spaces, The Journal of Machine Learning Research, 3, p.323-359, 3/1/2003
Mark D. Plumbley , Samer A. Abdallah , Thomas Blumensath , Michael E. Davies, Sparse representations of polyphonic music, Signal Processing, v.86 n.3, p.417-431, March 2006
Pavel Kisilev , Michael Zibulevsky , Yehoshua Y. Zeevi, A multiscale framework for blind separation of linearly mixed signals, The Journal of Machine Learning Research, 4, 12/1/2003
Alexander J. Smola , Sebastian Mika , Bernhard Schlkopf , Robert C. Williamson, Regularized principal manifolds, The Journal of Machine Learning Research, 1, p.179-209, 9/1/2001
Pavel Kisilev , Michael Zibulevsky , Yehoshua Y. Zeevi, A multiscale framework for blind separation of linearly mixed signals, The Journal of Machine Learning Research, v.4 n.7-8, p.1339-1364, October 1 - November 15, 2004
Gianluca Monaci , scar Divorra Escoda , Pierre Vandergheynst, Analysis of multimodal sequences using geometric video representations, Signal Processing, v.86 n.12, p.3534-3548, December 2006
Model , Michael Zibulevsky, Signal reconstruction in sensor arrays using sparse representations, Signal Processing, v.86 n.3, p.624-638, March 2006
Mrio A. T. Figueiredo, Adaptive Sparseness for Supervised Learning, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.25 n.9, p.1150-1159, September
Joseph F. Murray , Kenneth Kreutz-Delgado, Learning Sparse Overcomplete Codes for Images, Journal of VLSI Signal Processing Systems, v.46 n.1, p.1-13, January   2007
Joseph F. Murray , Kenneth Kreutz-Delgado, Learning Sparse Overcomplete Codes for Images, Journal of VLSI Signal Processing Systems, v.45 n.1-2, p.97-110, November  2006
Alex J. Smola , Bernhard Schlkopf, A tutorial on support vector regression, Statistics and Computing, v.14 n.3, p.199-222, August 2004
Zhe Chen , Simon Haykin, On different facets of regularization theory, Neural Computation, v.14 n.12, p.2791-2846, December 2002
Ivana Radulovic , Pascal Frossard, Multiple description coding with redundant expansions and application to image communications, Journal on Image and Video Processing, v.2007 n.1, p.8-8, January 2007
Olga Sorkine , Daniel Cohen-Or , Dror Irony , Sivan Toledo, Geometry-Aware Bases for Shape Approximation, IEEE Transactions on Visualization and Computer Graphics, v.11 n.2, p.171-180, March 2005
Kjersti Engan , Karl Skretting , John Hkon Husy, Family of iterative LS-based dictionary learning algorithms, ILS-DLA, for sparse signal representation, Digital Signal Processing, v.17 n.1, p.32-49, January, 2007
Gunnar Rtsch , Sebastian Mika , Bernhard Schlkopf , Klaus-Robert Mller, Constructing Boosting Algorithms from SVMs: An Application to One-Class Classification, IEEE Transactions on Pattern Analysis and Machine Intelligence, v.24 n.9, p.1184-1199, September 2002
Hongxing Zou , Dianjun Wang , Xianda Zhang , Yanda Li, Nonnegative time-frequency distributions for parametric time-frequency representations using semi-affine transformation group, Signal Processing, v.85 n.9, p.1813-1826, September 2005
matching pursuit;time-frequency analysis;overcomplete signal representation;multiscale edges;time-scale analysis;interior-point methods for linear programming;denoising;wavelet packets;wavelets;total variation denoising;ell1 norm optimization;cosine packets