A simple sparse approximation problem requests an approximation of a given input signal as a linear combination of T elementary signals drawn from a large, linearly dependent collection. An important generalization is simultaneous sparse approximation. Now one must approximate several input signals at once using different linear combinations of the same T elementary signals. This formulation appears, for example, when analyzing multiple observations of a sparse signal that have been contaminated with noise. A new approach to this problem is presented here: a greedy pursuit algorithm called simultaneous orthogonal matching pursuit. The paper proves that the algorithm calculates simultaneous approximations whose error is within a constant factor of the optimal simultaneous approximation error. This result requires that the collection of elementary signals be weakly correlated, a property that is also known as incoherence. Numerical experiments demonstrate that the algorithm often succeeds, even when the inputs do not meet the hypotheses of the proof

Gilbert, A. C.

Strauss, M. J.

Tropp, J. A.

English

Caltech Authors - Main

SIMULTANEOUS SPARSE APPROXIMATION VIA GREEDY PURSUIT
J. A. Tropp, A. C. Gilbert
Department of Mathematics
The University of Michigan
Ann Arbor, MI 48109
M. J. Strauss
Departments of Mathematics and EECS
The University of Michigan
Ann Arbor, MI 48109
ABSTRACT
A simple sparse approximation problem requests an ap-
proximation of a given input signal as a linear combination
of T elementary signals drawn from a large, linearly de-
pendent collection. An important generalization is simul-
taneous sparse approximation. Now one must approximate
several input signals at once using different linear combi-
nations of the same T elementary signals. This formulation
appears, for example, when analyzing multiple observations
of a sparse signal that have been contaminated with noise.
A new approach to this problem is presented here: a
greedy pursuit algorithm called Simultaneous Orthogonal
Matching Pursuit. The paper proves that the algorithm cal-
culates simultaneous approximations whose error is within
a constant factor of the optimal simultaneous approxima-
tion error. This result requires that the collection of ele-
mentary signals be weakly correlated, a property that is also
known as incoherence. Numerical experiments demonstrate
that the algorithm often succeeds, even when the inputs do
not meet the hypotheses of the proof.
1. INTRODUCTION
We work in the complex inner-product space Cd, which is
called the signal space. We write 〈·, ·〉 for the usual inner
product and ‖·‖2 for the associated norm. The symbol ‖·‖p
indicates the p vector norm, while ‖·‖p,q is the norm on
linear operators mapping p to q. We use ∗ for the complex
conjugate transpose of vectors and matrices.
A dictionary D is a finite collection of unit-norm ele-
mentary signals, called atoms, that spans the signal space.
Each atom is denoted ϕω, where ω is drawn from an in-
dex set Ω. The number of atoms N is typically much larger
than the dimension d of the signal space. We also define the
d×N dictionary matrix Φ whose columns are atoms.
Suppose that S is a d × K matrix whose columns are
input signals. We wish to approximate all K input signals
using different linear combinations of the same T atoms.
Typically, T is much smaller than the dimension of the sig-
nal space, so the approximation is sparse. More precisely,
the simultaneous sparse approximation problem (SSA) elic-
its an N×K coefficient matrix C that solves the mathemat-
ical program
min
C
‖S − Φ C‖2F subject to
the matrix C has at most T nonzero rows. (SSA)
The squared Frobenius matrix norm ‖·‖2F returns the sum of
the squares of the entries in a matrix.
The (SSA) problem arises if we are given multiple ob-
servations of a sparse input signal that are contaminated
with noise. For example, the k-th input signal might have
the form
sk = x + νk
where νk is a realization of some random process and where
x can be expressed using a linear combination of T atoms.
The goal is to identify the atoms that comprise x.
To solve (SSA), we propose a greedy pursuit method,
Simultaneous Orthogonal Matching Pursuit (S-OMP). For
general dictionaries, (SSA) cannot be solved without check-
ing every combination of T nonzero rows. This follows
from results in [1]. Nevertheless, we have been able to
prove that S-OMP correctly solves the simultaneous sparse
approximation problem, provided that the atoms are weakly
correlated. To quantify this property, we define the coher-
ence parameter of the dictionary,
µ
def= max
λ=µ
|〈ϕλ,ϕω〉| .
When the coherence parameter is small, each pair of atoms
is nearly orthogonal.
This paper provides the first proof that any algorithm
can obtain provably good solutions to (SSA). A simple ver-
sion of our result follows1. Suppose that the set Λopt in-
dexes the T atoms that appear in some solution to (SSA).
Then we may define the d×K matrix Aopt whose k-th col-
umn is the best approximation of the k-th input signal using
the T atoms listed in Λopt.
1Note that the present result does not yield an optimal bound for the
constant, even when the dictionary is orthonormal. A more subtle analysis
is necessary to achieve the improvement.
V - 7210-7803-8874-7/05/$20.00 ©2005 IEEE ICASSP 2005
Theorem 1 Assume that µT < 12 . After T steps, suppose
that Simultaneous Orthogonal Matching Pursuit returns the
approximation AT . Then the output error is bounded by
‖S − AT ‖F ≤
√
1 + C(D ,K, T ) ‖S − Aopt‖F ,
and the constant is no worse than
C(D ,K, T ) ≤ K T (1− µT )
(1− 2µT )2 .
In particular, if the optimal approximation error is zero, S-
OMP returns an approximation that achieves zero error.
The algorithm S-OMP performs much better in practice
than our theory predicts. Not only can we recover the input
signals when T is large, the error does not grow as quickly
as our bounds would suggest. We have discovered that mod-
erate noise levels create surprising difficulties for our algo-
rithm. To our knowledge, these are the first numerical ex-
periments performed on (SSA).
The rest of the paper expands on the claims of the in-
troduction. Section 2 provides a rigorous statement of our
greedy pursuit algorithm. In Section 3, we sketch the proof
that the algorithm constructs approximate solutions to (SSA),
and we discusses several other factors that affect its perfor-
mance. The paper concludes with Section 4, which summa-
rizes our numerical experiments.
2. THE ALGORITHM
Let us continue with a formal description of the algorithm.
Algorithm 2 (S-OMP)
INPUT:
• A d×K matrix S of input signals
• The number T of atoms in the approximation
OUTPUT:
• A set ΛT containing T indices
• A d×K approximation matrix AT
• A d×K residual matrix RT
PROCEDURE:
1. Initialize the residual matrix R0 = S , the index set
Λ0 = ∅, and the iteration counter t = 1.
2. Find an index λt that solves the easy optimization
problem
max
ω∈Ω
K∑
k=1
|〈Rt−1 ek,ϕω〉| .
We use ek to denote the k-th canonical basis vector.
3. Set Λt = Λt−1 ∪ {λt}.
4. Determine the orthogonal projector Pt onto the span
of the atoms indexed in Λt.
5. Calculate the new approximation and residual:
At = Pt S
Rt = S − At.
6. Increment t, and return to Step 2 if t ≤ T .
This procedure reduces to standard Orthogonal Matching
Pursuit [1] when K = 1.
Step 2 of the algorithm is referred to as the greedy se-
lection. The intuition behind maximizing the sum of abso-
lute correlations is that we wish to find an atom that con-
tributes the most energy to as many of the input signals as
possible. Note that this absolute sum can also be written as
‖Rt∗ϕω‖1. In contrast, Leviatan and Temlyakov [2] have
studied a greedy algorithm for (SSA) that picks an atom by
maximizing ‖Rt∗ϕω‖∞.
Steps 4 and 5 have been written to emphasize the con-
ceptual structure of the algorithm. It is possible to imple-
ment them much more efficiently using standard techniques
for least-squares problems. See [3, Ch. 5] for extensive de-
tails. It is important to note that each column of the residual
Rt is orthogonal to the atoms indexed in Λt. Therefore, no
atom is ever chosen twice.
3. PROOF OF CORRECTNESS
We will develop a condition which guarantees that S-OMP
selects an optimal atom at iteration t. From this condition,
it is easy to prove that Simultaneous Orthogonal Matching
Pursuit can compute approximate solutions to (SSA).
Theorem 3 Assume that µT < 12 , and fix a signal matrix
S . At iteration t, suppose that S-OMP has chosen t optimal
atoms, and let At be the current approximation of the signal
matrix. At iteration (t + 1), greedy selection will identify
another optimal atom provided that
‖S − At‖2F > ‖S − Aopt‖2F +
T (1− µT )
(1− 2µT )2 ‖Φ
∗(S − Aopt)‖2∞,∞ (1)
where Aopt denotes an optimal approximation of the signal
matrix using T atoms.
In words, the algorithm selects another optimal atom when-
ever the current approximation is somewhat worse than an
optimal approximation. We interpret ‖Φ∗(S − Aopt)‖∞,∞
as the maximum total correlation between a fixed atom and
the residuals left over from the optimal approximation.
V - 722
Proof. Suppose that some solution of (SSA) involves the
T atoms indexed in Λopt. Define the d × T matrix Φopt
whose columns are the atoms listed by Λopt. Let the d ×
(N − T ) matrix Ψopt contain the remaining atoms. Recall
the definition Rt = S − At.
First, observe that each row of the matrix (Φopt∗Rt) lists
the inner products between a fixed atom in Λopt and the
columns of Rt. The rows of the matrix (Ψopt∗Rt) have an
analogous interpretation. The (∞,∞) matrix norm returns
the maximum absolute row sum of its argument, and so the
algorithm chooses another optimal atom if and only if the
ratio
ρ
def=
‖Ψopt∗Rt‖∞,∞
‖Φopt∗Rt‖∞,∞
(2)
is strictly less than one. We must ensure that ρ < 1.
Rewrite Rt = (S − Aopt) + (Aopt − At). Substitute
this expression into (2). The term Φopt
∗(S − Aopt) van-
ishes from the denominator because of orthogonality. Ap-
ply the triangle inequality to the numerator to see that ρ is
no greater than
‖Ψopt∗(Aopt − At)‖∞,∞
‖Φopt∗(Aopt − At)‖∞,∞
+
‖Ψopt∗(S − Aopt)‖∞,∞
‖Φopt∗(Aopt − At)‖∞,∞
.
(3)
Now we bound the first fraction in (3). Let Φopt
+ de-
note the generalized inverse of Φopt. Using an argument
analogous to that in [4, Thm. 3.1], we discover that
‖Ψopt∗(Aopt − At)‖∞,∞
‖Φopt∗(Aopt − At)‖∞,∞
≤ ∥∥Φopt+Ψopt
∥∥
1,1
. (4)
It can be shown that the second fraction in (3) is no greater
than
‖Ψopt∗(S − Aopt)‖∞,∞∥∥Φopt+
∥∥−1
2,1
‖Aopt − At‖F
. (5)
In the numerator, we replace Ψopt with Φ = [Φopt Ψopt],
using the fact that the columns of (S−Aopt) are orthogonal
to the columns of Φopt.
We substitute the bounds (4) and (5) into (3) and per-
form some algebraic manipulations to find a condition that
ρ < 1. To complete the argument, we introduce the coher-
ence estimates developed in [5, Sec. 3].

Corollary 4 Assume that µT < 12 . Given any input matrix
S , Simultaneous Orthogonal Matching Pursuit will always
construct a T -term approximation AT that satisfies the er-
ror bound
‖S − AT ‖2F ≤ ‖S − Aopt‖2F +
T (1− µT )
(1− 2µT )2 ‖Φ
∗(S − Aopt)‖2∞,∞
where Aopt is an optimal T -term approximation of S .
From here, one reaches Theorem 1 by noting that
‖Φ∗(S − Aopt)‖2∞,∞ ≤ K ‖S − Aopt‖2F
because the columns of Φ all have unit norm.
4. NUMERICAL EXPERIMENTS
To evaluate the performance of our algorithm, we have tested
it with three types of input signal. Each type is a variant on
the form sk = xk + νk, where νk is random noise and
where xk can be expressed using a linear combination of
T atoms (possibly the same for each k as given in the intro-
ductory example). In each dimension d, our dictionary is the
collection of d complex exponentials and d impulses. That
is, ϕω[t] = e2πi tω/d for ω = 1, . . . , d and ϕλ[t] = δλ[t]
for λ = 1, . . . , d. Note that this dictionary has coherence
µ = 1/
√
d.
We begin with input signals of the form
sk =
T∑
j=1
αjk ϕωjk .
For each signal sk, we select T atoms independently and
uniformly from the dictionary. The coefficients αjk are cho-
sen from iid normal distributions. Our algorithm is, there-
fore, searching for the best T atoms with which to repre-
sent K signals, each of which is a linear combination of
T atoms. We have observed that our algorithm always re-
covers T atoms from the collection of approximately KT
distinct atoms that participate in the K input signals. All
of the error in the residual is due to the fact that the input
signals involve more atoms that we are allowed to use.
The second type of input signal has the form
sk =
T∑
j=1
αjk ϕωj
For all K signals, we use the same core of T atoms, but the
coefficients αjk are chosen from iid normal distributions.
For these experiments, we fixed the dimension of the sig-
nal space at d = 128 and the number of signals at K = 2.
We vary the value of T to explore how many core atoms we
can successfully recover with our algorithm. For each set
of parameters, we performed 1000 independent trials. We
computed the Hamming distance between the set of recov-
ered atoms and the core set. Hamming distance zero means
that we recover the entire core set, while distance one means
that we fail to recover any of the core atoms. In Figure 1,
we plot the average Hamming distance as a function of T .
The error bars mark one standard deviation from the mean.
We can see from this figure that our theoretical bounds are
far too pessimistic. Even for T = 90 (out of a possible 128
atoms), we typically recover most of the core set.
V - 723
60 70 80 90 100 110 120
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
number of vectors T
av
era
ge
 Ha
mm
ing
 di
sta
nc
e
Fig. 1. The average Hamming distance between the core
set of the vectors and the recovered set as a function of the
number of vectors T in the core set. (Input type two)
The third input type has the form
sk =
T∑
j=1
αj ϕωj + νk.
That is, we choose T atoms at random and form a linear
combination with random coefficients αj ∈ {±1}. Then
we construct K input signals by corrupting the original sig-
nal with iid additive white Gaussian noise νk. For these
experiments, we fix the dimension d = 256; we vary T
from 2 to 4; we vary K from 2 to 6; and we examine SNR
values of 10, 13, 16, and 20 dB. For each parameter set,
we perform 1000 trials. Figure 2 displays the average Ham-
ming distance as a function of the number of signals. For
each value of T , we use a distinct line type (e.g., dashed)
so the four dashed lines correspond to the four SNR val-
ues. Naturally, the Hamming distance increases as SNR de-
creases. Observe that, independent of the number of core
atoms T and the SNR, we recover the core signal better
when we have more observations. Furthermore, the pres-
ence of noise has a significant effect on the performance of
the algorithm. The previous example showed that we can
often recover core sets of atoms that are almost as large as
the dimension of the signal space. Yet for moderate SNR
(e.g., 13dB), we cannot reliably recover three atoms in a
256-dimensional signal space. With the parameter settings
we have chosen, our theoretical results predict that
‖S − At‖2F∑K
k=1 ‖νk‖22
≤ 1 + 3KT.
To see if this bound accurately predicts the dependence on
K and T , we plot in Figure 3 the total relative error as a
function of the number of signals K. For each T , we use
a different line type. The two groups of lines represent the
extreme SNR values (10 and 20 dB). The plot shows that
the size of T has a negligible effect on the error. That is, the
theoretical bounds reflect a dependence on T that is absent
in the empirical evidence. Again, our algorithm performs
better than the theoretical results might lead us to believe.
2 2.5 3 3.5 4 4.5 5 5.5 6
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
number of signals
ha
m
m
in
g 
di
st
an
ce
T=2
T=3
T=4
Fig. 2. The average Hamming distance between the core set
of vectors and the recovered set as a function of the number
of signals and the SNR. (Input type three)
2 2.5 3 3.5 4 4.5 5 5.5 6
1.5
2
2.5
3
3.5
4
4.5
5
5.5
6
number of signals
to
ta
l r
el
at
ive
 e
rro
r
T=1
T=2
T=3
T=4
Fig. 3. The total relative error as a function of the number
of signals and the number of core vectors for two values of
SNR. (Input type three)
5. REFERENCES
[1] G. Davis, S. Mallat, and M. Avellaneda, “Greedy adap-
tive approximation,” J. Constr. Approx., vol. 13, pp.
57–98, 1997.
[2] D. Leviatan and V. N. Temlyakov, “Simultaneous
approximation by greedy algorithms,” IMI Report
2003:02, Univ. of South Carolina at Columbia, 2003.
[3] G. H. Golub and C. F. Van Loan, Matrix Computations,
Johns Hopkins University Press, 3rd edition, 1996.
[4] J. A. Tropp, “Greed is good: Algorithmic results for
sparse approximation,” IEEE Trans. Inform. Theory,
vol. 50, no. 10, October 2004, To appear.
[5] J. A. Tropp, “Just relax: Convex programming methods
for subset selection and sparse approximation,” ICES
Report 04-04, The University of Texas at Austin, 2004.
V - 724


Simultaneous sparse approximation via greedy pursuit

Greed is good: Algorithmic results for sparse approximation,”

Greedy adaptive approximation,”

Just relax: Convex programming methods for subset selection and sparse approximation,”

Simultaneous approximation by greedy algorithms,”

http://authors.library.caltech.edu/9043/1/TROicassp05.pdf

Simultaneous sparse approximation via greedy pursuit

Abstract

Similar works

Full text

Available Versions

Caltech Authors - Main