remarkRemark hypothesisHypothesis claimClaim Fast OMP for Exact Recovery and
ApproximationH. Yu, J. He, and M. Cheng

# Fast Orthogonal Matching Pursuit for Exact Recovery and Sparse Approximation

Huiyuan Yu  Jia He  Maggie X. Cheng Department of Applied Mathematics, College
of Computing, Illinois Institute of Technology, Chicago, IL, 60616. ().
<EMAIL_ADDRESS>
###### Abstract

Orthogonal Matching Pursuit (OMP) has been a powerful method in sparse signal
recovery and approximation. However OMP suffers computational issue when the
signal has large number of non-zeros. This paper advances OMP in two fronts:
it offers a fast algorithm for the orthogonal projection of the input signal
at each iteration, and a new selection criterion for making the greedy choice,
which reduces the number of iterations it takes to recover the signal. The
proposed modifications to OMP directly reduce the computational complexity.
Experiment results show significant improvement over the classical OMP in
computation time. The paper also provided a sufficient condition for exact
recovery under the new greedy choice criterion. For general signals that may
not have sparse representations, the paper provides a bound for the
approximation error. The approximation error is at the same order as OMP but
is obtained within fewer iterations and less time.

###### keywords:

Greedy Algorithm, Compressive Sensing, Sparse Signal Recovery, Approximation,
Orthogonal Matching Pursuit, Restricted Isometry Property

65F10, 65F20, 68Q25, 68R10, 68W25

##  1 Introduction

Let $\bm{x}$ be a $d$-dimensional real signal. Suppose there is a real
measurement matrix $\Phi\in\mathbb{R}^{N\times d}$, through which we can
obtain an $N$-dimensional measurement $\bm{y}=\Phi\bm{x}$. Usually $N<d$,
which presents an underdetermined system. How to reconstruct the original
signal $\bm{x}$ from an underdetermined system? If $\bm{x}$ is sparse, then by
exploiting sparsity, we may be able to find a unique solution. $\bm{x}$ is
called a $k$-sparse signal if $\bm{x}$ has at most $k$ non-zero components.

The measurement matrix $\Phi$ is also called a dictionary, and each column
$\bm{\varphi}$ of the dictionary called an atom. Let
$\mathcal{J}=\\{1,\ldots,d\\}$ represent the index set of all atoms in the
dictionary. If the dictionary is overcomplete, there are many representations
of $\bm{y}=\sum\limits_{\gamma\in\mathcal{J}}a_{\gamma}\bm{\varphi}_{\gamma}$.
Intuitively, we would like to find the sparsest solution:
$\min\limits_{\bm{x}}\|\bm{x}\|_{0}\textrm{ subject to }\bm{y}=\Phi\bm{x}$,
but it is an NP-hard problem. Different optimization principles lead to
different sparse representations of $\bm{y}$, for example, basis pursuit (BP)
[12, 7, 6] and the method of frames (MOF) [9] among many others [14, 23]:

  * •

Find a representation of the input signal whose coefficients have the minimal
$\ell_{1}$ norm.

(BP) $\min_{\bm{x}}\|\bm{x}\|_{1}\;\;\textrm{ subject to }\;\bm{y}=\Phi\bm{x}$

  * •

Find a representation of the input signal whose coefficients have the minimal
$\ell_{2}$ norm.

(MOF) $\min_{\bm{x}}\|\bm{x}\|_{2}\;\;\textrm{ subject to
}\;\bm{y}=\Phi\bm{x}$

BP and MOF both provide convex relaxation to the $\ell_{0}$ norm minimization
problem, however, neither of them provides the sparsest solution, except those
satisfying the sparsity condition specified in [11].

Matching Pursuit (MP) [16] uses an iterative procedure that directly addresses
the sparsity issue. Orthogonal Matching Pursuit (OMP) [19, 21, 22] inherits
the greedy approach from MP that selects an atom with the maximal correlation
with the residual at present, but improves over the standard MP by adding
least square minimization at each iteration. Let $\Gamma$ be the index set of
atoms found so far, the least square estimation is used for computing the
orthogonal projection of the input signal $\bm{y}$ onto the subspace spanned
by the atoms indexed by $\Gamma$:

(OMP)
$\min_{\bm{x}_{\Gamma}}\|\bm{y}-\Phi_{\Gamma}\bm{x}_{\Gamma}\|_{2}^{2}\;\;\text{
with }|\Gamma|\leq k$

The complete OMP algorithm is included in the Appendix for comparison purpose.

OMP has been shown to have better results than MP. Many variations of OMP have
been developed [13, 15, 8, 18, 17]. Under certain conditions OMP provides
recovery guarantee [4, 22, 21, 10, 3, 5]. The excellent performance of OMP
results from the orthogonal projection of $\bm{y}$ onto the subspace spanned
by the atoms selected so far. The least square solution is obtained by
$\bm{x}_{\Gamma}=\Phi_{\Gamma}^{+}\bm{y}$. As $\Gamma$ increases, solving the
least square problem significantly increases the computational load. In this
paper, we propose two fundamental improvements over classical OMP to avoid the
high complexity of computing pseudo inverse over an increasing-sized matrix:

  * •

When solving the least square problem at each iteration, instead of computing
$\Phi_{\Gamma}^{+}\bm{y}$ over the entire support $\Gamma$, it uses successive
regression over a single atom. It makes the same greedy choice as OMP does at
each iteration, but is much faster due to reduced computation load. The
proposed algorithm is called OMP-SR.

  * •

Selecting a block of atoms instead of one atom at each iteration. An important
invariant held by the classical OMP is that only one optimal atom can be
selected in each iteration so that under certain condition OMP finds $k$ non-
zeros in $k$ iterations. Without this invariant, the sufficient condition may
not lead to exact recovery. The blocked version will extend the greedy choice
to multiple atoms at each iteration but still preserve the convergence
property of OMP. The proposed algorithm is called blocked successive
regression (BSR).

In general, the measurement $\bm{y}$ is often with noise. A general signal may
be represented as the linear combination of atoms from the dictionary with
additive noise,

$\bm{y}=\Phi{\bm{x}}+\bm{\varepsilon}.$

We are interested in the best approximation of $\bm{y}$ using a linear
combination of atoms. Best approximation is the one with the smallest
approximation error measured by the $\ell_{2}$ norm of the residual, and
hence, the optimization principle is,

(Sparse Approximation) $\min_{\bm{x}}\|\bm{y}-\Phi\bm{x}\|_{2}\;\;\textrm{
subject to }\|\bm{x}\|_{0}\leq\\!k$

###  1.1 Major Results

  * •

Improvement on OMP, which preserved the excellent stability and performance
guarantee of OMP but are more computationally efficient and hence much faster.

  * •

Theories of sufficient conditions for both exact recovery and sparse
approximation. These sufficient conditions are specified in terms of the
properties of the dictionary. Compared to the popular Restricted Isometry
Property (RIP), these conditions are more computationally feasible to check
[20, 2, 1].

  * •

Approximation error of BSR for general signals. The error bound extends the
result of [21] to the blocked version and is shown to be at the same order as
OMP, but is obtained with fewer number of iterations. Table 3 shows the
approximation error, iterations, and running time.

###  1.2 Notation

  * •

$\Phi^{\top}$: transpose of matrix $\Phi$

  * •

$\Phi^{+}$: pseudo inverse of matrix $\Phi$

  * •

$(\Phi^{\top}\Phi)^{-1}$: inverse of matrix $(\Phi^{\top}\Phi)$

  * •

$\|A\|_{p\rightarrow q}=\sup\limits_{{\bm{x}}\neq
0}\frac{\|A{\bm{x}}\|_{q}}{\|{\bm{x}}\|_{p}}$: operator norm of matrix $A$.

  * •

$\|A\|_{p\rightarrow p}$ is abbreviated to $\|A\|_{p}$ .

##  2 Orthogonal Matching Pursuit through Successive Regression (OMP-SR)

OMP-SR is a fast implementation of OMP. When solving the least square problem
at each iteration of OMP, it avoids the expensive computation for the pseudo
inverse of $\Phi_{J^{t-1}}$; instead, it only projects onto the atom selected
in the current iteration via univariate regression, and then updates the
coefficients of atoms selected in previous iterations through a backtracking
procedure: $b_{l}=\beta_{l}-\sum\limits_{k=l+1}^{t}b_{k}\gamma_{l,k}$(see
Algorithm 1), where
$\beta_{t}=\frac{\langle\bm{z}_{t},\bm{y}\rangle}{\langle\bm{z}_{t},\bm{z}_{t}\rangle}$
is the coefficient newly obtained in the current iteration, $b_{l}$ is the
updated coefficient for the atoms selected in previous iterations. Note that
the inner product $\langle\bm{z}_{l},\bm{z}_{l}\rangle$ does not need to be
recomputed. It only needs to be computed once, that is when we compute
$\beta_{l}$ in the $l$-th iteration.

Algorithm 1 OMP-SR

Initialization: $\bm{a}_{0}=\bm{z}_{0}=\bm{1}$, $\bm{r}^{0}=\bm{y}$,
$J^{0}=\phi$

for $t=1$ to $\kappa$ do

Choose $j^{t}=\mathrm{arg}\max\limits_{j\in\mathcal{J}\setminus
J^{t-1}}\left|\sum\limits_{i=1}^{N}\Phi_{i,j}\bm{r}_{i}^{t-1}\right|$

Let $\bm{a}_{t}$ be the $j^{t}$-th column of matrix $\Phi$.

Regress $\bm{a}_{t}$ on $\bm{z}_{l}$ and get coefficients
$\gamma_{l,t}=\frac{\langle\bm{z}_{l},\bm{a}_{t}\rangle}{\langle\bm{z}_{l},\bm{z}_{l}\rangle},\;\text{
for }l=0,\dots,t-1$

Compute $\bm{z}_{t}=\bm{a}_{t}-\sum\limits_{l=0}^{t-1}\gamma_{l,t}\bm{z}_{l}$

Regress $\bm{y}$ on $\bm{z}_{t}$ to get
$\beta_{t}=\frac{\langle\bm{z}_{t},\bm{y}\rangle}{\langle\bm{z}_{t},\bm{z}_{t}\rangle}$

Let $b_{t}=\beta_{t}$

if $t>1$ then

for $l=t-1$ to $1$ do

$b_{l}=\beta_{l}-\sum\limits_{k=l+1}^{t}b_{k}\gamma_{l,k}$

end for

end if

Update index set $J^{t}=J^{t-1}\bigcup\\{j^{t}\\}$

Update residual $\bm{r}^{t}=\bm{y}-\sum\limits_{l=1}^{t}b_{l}\Phi_{j^{l}}$

end for

Let $\bm{x}_{j^{t}}=b_{t}$ for $t=1,\ldots,\kappa$, and let $\bm{x}_{j}=0$ for
$j\notin J^{\kappa}$

Return $\bm{x}$

OMP-SR selects the same atom and generates the same residual as OMP does at
each iteration, and therefore returns the same result as OMP. OMP-SR starts to
show performance gain over OMP when the number of non-zeros in $\bm{x}$
increases due to not having to compute the pseudo inverse of a growing matrix.

##  3 Blocked Successive Regression (BSR)

BSR builds on the idea of successive regression in OMP-SR but selects a block
of atoms at each iteration. The block size $c$ is a hyper parameter, usually
decided by a grid search. The algorithm is still greedy in nature: in each
iteration it selects the atoms that have the largest correlations with the
residual measured by the $\ell_{2}$ norm. Each iteration of BSR performs an
orthogonal projection of $\bm{y}$ over $c$ newly selected atoms, instead of
over $|\Gamma|$ atoms, which could be costly as $|\Gamma|$ increases with
iterations. Subsequently the coefficients for atoms selected in previous
iterations are updated through
${b}_{i}={\beta}_{i}-\sum\limits_{k=l+1}^{t}\sum\limits_{j\in\Gamma_{k}}{b}_{j}\gamma_{i,j}$
(see Algorithm 2).

Algorithm 2 BSR

Initialization: $\bm{r}^{0}=\bm{y}$, $\Gamma=\phi$, $\bm{z}_{0}=\bm{1}$

for $t=1$ to $\kappa$ do

$\Gamma_{t}=\mathrm{arg}\max\limits_{\begin{subarray}{c}{|\Omega|=c}\\\
{\Omega\subset\mathcal{J}\setminus\Gamma}\end{subarray}}\left\|\Phi^{\top}_{\Omega}\bm{r}^{t-1}\right\|_{2}$

for each $j\in\Gamma_{t}$ do

Let $\bm{a}_{j}$ be the $j$-th column of matrix $\Phi$.

Regress $\bm{a}_{j}$ on $\bm{z}_{0}$ to get coefficient
$\gamma_{0,j}=\frac{\langle\bm{z}_{0},\bm{a}_{j}\rangle}{\langle\bm{z}_{0},\bm{z}_{0}\rangle}$

Compute $\bm{z}_{j}=\bm{a}_{j}-\gamma_{0,j}\bm{z}_{0}$

if $t>1$ then

Regress $\bm{a}_{j}$ on $Z_{\Gamma_{l}}$ to get coefficients
$\bm{\gamma}_{\Gamma_{l},j}=Z_{\Gamma_{l}}^{+}\bm{a}_{j},\;\text{for
}l=1,\ldots,t-1$

Compute
$\bm{z}_{j}=\bm{z}_{j}-\sum\limits_{l=1}^{t-1}\sum\limits_{i\in\Gamma_{l}}\gamma_{i,j}\bm{z}_{i}$

end if

end for

Regress $\bm{y}$ on $Z_{\Gamma_{t}}$ to get coefficients
$\bm{\beta}_{\Gamma_{t}}=Z_{\Gamma_{t}}^{+}\bm{y}$

Let $\bm{b}_{\Gamma_{t}}=\bm{\beta}_{\Gamma_{t}}$

if $t>1$ then

for $l=t-1$ to $1$ do

for $i\in\Gamma_{l}$ do

${b}_{i}={\beta}_{i}-\sum\limits_{k=l+1}^{t}\sum\limits_{j\in\Gamma_{k}}{b}_{j}\gamma_{i,j}$

end for

end for

end if

Update index set $\Gamma=\Gamma\bigcup\Gamma_{t}$

Update residual $\bm{r}^{t}=\bm{y}-\Phi_{\Gamma}\bm{b}_{\Gamma}$

Break if $\|\bm{r}^{t}\|_{2}\leq\delta$

end for

Let $\bm{x}_{\Gamma}=\bm{b}_{\Gamma}$, and let
$\bm{x}_{\mathcal{J}\setminus\Gamma}=0$

Return $\bm{x}$

Making the greedy choice
$\Gamma_{t}=\mathrm{arg}\max\limits_{\begin{subarray}{c}{|\Omega|=c}\\\
{\Omega\subset\mathcal{J}\setminus\Gamma}\end{subarray}}\left\|\Phi^{\top}_{\Omega}\bm{r}^{t-1}\right\|_{2}$
does not incur combinatorial complexity. In fact, it only needs to compute the
absolute inner product values
$|\langle{\bm{\varphi}}_{j},{\bm{r}}^{t-1}\rangle|$ for
${\bm{\varphi}}_{j}\in\Phi\setminus\Phi_{\Gamma}$ and take the largest $c$
columns.

The BSR algorithm halts if the residual becomes too small or it has exhausted
$\kappa$ iterations, which amounts to two of the three halting rules listed in
[17] for matching pursuit type of algorithms.

##  4 Conditions for Exact Recovery

###  4.1 Background

Assume there are $k$ non-zero entries in a $d$-dimensional signal $\bm{x}$,
and $k\ll d$. Let $\Lambda_{\text{opt}}=\\{i_{1},\ldots,i_{k}\\}$ be the set
of indices for the non-zero entries of $\bm{x}$. Without loss of generality,
we can partition the measurement matrix as $\Phi=[\Phi_{\text{opt}}|\Psi]$ so
that $\Phi_{\text{opt}}$ has $k$ columns,
$\Phi_{\text{opt}}=[\varphi_{i_{1}},\ldots,\varphi_{i_{k}}]$, and $\Psi$ has
the remaining $d-k$ columns.

In the absence of noise, the measured signal $\bm{y}$ has a sparse
representation:
$\bm{y}=\Phi\bm{x}=\sum\limits_{j\in\Lambda_{\mathrm{opt}}}a_{j}\bm{\varphi}_{j}$.
Exact recovery aims to recover the coefficients $a_{j}$ for all atoms indexed
by $\Lambda_{\mathrm{opt}}$, which are the non-zero entries in $\bm{x}$.

####  4.1.1 The Exact Recovery Condition of OMP

In [21], the exact recovery condition for OMP is established. Under this
condition, OMP can find all optimal atoms. Let $\bm{r}$ be the residual at the
current iteration. Initially before the first iteration, $\bm{r}$ is set to be
$\bm{r}=\bm{y}$.

To see how the OMP algorithm locates the columns of $\Phi_{\mathrm{opt}}$, we
define the greedy choice ratio as follows:

(1) $\rho(\bm{r})\stackrel{{\scriptstyle\text{ def
}}}{{\scalebox{2.0}[1.0]{=}}}\frac{\left\|{\Psi^{\top}}\bm{r}\right\|_{\infty}}{\left\|\Phi_{\mathrm{opt}}^{\top}\bm{r}\right\|_{\infty}}=\frac{\max_{\bm{\psi}}|\langle\bm{\psi},\bm{r}\rangle|}{\left\|\Phi_{\mathrm{opt}}^{\top}\bm{r}\right\|_{\infty}},$

where $\bm{\psi}$ ranges over columns of $\Psi$.

A sufficient condition for OMP algorithm to succeed is that $\rho(\bm{r})<1$
holds for all iterations. When $\rho(\bm{r})<1$ holds, the greedy algorithm
will pick a column from $\Phi_{\mathrm{opt}}$ that has not been selected by
the algorithm. The columns that have been selected so far will not be the
column that gives the maximum absolute inner product
$\left\|\Phi_{\mathrm{opt}}^{\top}\bm{r}\right\|_{\infty}$, because the
residual $\bm{r}$ is orthogonal to them. Each iteration of OMP will select one
new optimal column, so the OMP algorithm is guaranteed to locate $k$ optimal
columns in $k$ iterations.

If the ratio fails to hold, e.g., if $\rho(\bm{r})=1$, the algorithm is not
guaranteed to make the correct choice, since the absolute inner product of the
residual with a non-optimal column is the same as the absolute inner product
with an optimal column and there is no reason to prefer one over another. A
non-optimal column could be selected in this case. If $\rho(\bm{r})>1$, the
algorithm will certainly select a non-optimal column. In summary, if
$\rho(\bm{r})<1$ fails to hold at any iteration, the algorithm could fail to
locate $k$ optimal columns in $k$ iterations.

####  4.1.2 The Exact Recovery Condition of OMP-SR

Algorithm OMP-SR essentially is an OMP algorithm with fast implementation: it
starts with the same initial residual $\bm{r}^{0}$ and selects the same atom
in the next iteration, so the residual $\bm{r}^{t}$ after the $t$-th iteration
is the same. Since $\bm{r}^{t}$ is used as input to the next iteration when
choosing a column, the next iteration will results in the same residual
$\bm{r}^{t+1}$. By induction, after $k$ iterations the algorithm returns the
same result as OMP does. The exact recovery condition for OMP-SR is the same
as for OMP.

###  4.2 The Exact Recovery Conditions of BSR

BSR is essentially a greedy algorithm, which makes a greedy choice at each
iteration, except that BSR selects a block of columns at each iteration with a
fixed block size $c$ ($c\geq 1$). If $c=1$, BSR reduces to OMP-SR. We have
learned that under the condition of $\rho(\bm{r})<1$, OMP and OMP-SR can find
one optimal column in each iterations. Then for BSR, under what condition will
each iteration of BSR only select the optimal columns from
$\Phi_{\mathrm{opt}}$ except the last iteration? This is the best case, in
which BSR can locate all optimal columns within $\lceil k/c\rceil$ iterations.
We call the condition for the best case as the strong exact recovery condition
for BSR.

####  4.2.1 A Strong Exact Recovery Condition for BSR

Recall that $\Phi=[\Phi_{\text{opt}}|\Psi]$ so that $\Phi_{\text{opt}}$ has
the $k$ optimal columns, and $\Psi$ has the remaining $d-k$ columns. Let
$\bm{r}$ denote the residual at the current iteration before the greedy choice
is made.

For a fixed block size $c$, the greedy choice ratio is defined as follows:

(2) $\rho_{c}(\bm{r})\stackrel{{\scriptstyle\text{ def
}}}{{\scalebox{2.0}[1.0]{=}}}\frac{\max\limits_{\Omega_{1}}\left\|\Phi^{\top}_{\Omega_{1}}\bm{r}\right\|_{2}}{\max\limits_{\Omega_{2}}\left\|\Phi^{\top}_{\Omega_{2}}\bm{r}\right\|_{2}},$

such that $|\Omega_{1}|=|\Omega_{2}|=c$ and
$|\Omega_{2}\cap\Lambda_{\mathrm{opt}}|>|\Omega_{1}\cap\Lambda_{\mathrm{opt}}|$,
i.e., $\Omega_{2}$ has at least one more optimal column than $\Omega_{1}$.
Given a $k$-sparse signal, BSR can recover the signal within $\lceil
k/c\rceil$ iterations if the following condition holds.

######  Theorem 4.1 (The strong exact recovery condition for BSR).

A sufficient condition for BSR to recover a $k$-sparse signal within $\lceil
k/c\rceil$ iterations is that

(3) $\rho_{c}(\bm{r})<1$

holds for all iterations.

Proof: In all but the last iteration, BSR algorithm will select $c$ columns
from $\Phi_{\mathrm{opt}}$ that have not been selected by the algorithm. This
can be observed from two aspects: 1) the columns that have been selected so
far will not contribute to the maximum $\ell_{2}$-norm
$\max\limits_{\Omega_{2}}\left\|\Phi^{\top}_{\Omega_{2}}\bm{r}\right\|_{2}$
since the residual $\bm{r}$ is orthogonal to them, therefore the maximum value
have to be achieved by the new columns the algorithm has not seen so far; and
2) with $\rho_{c}(\bm{r})<1$, the BSR algorithm will select as many optimal
columns as possible in each iteration. Therefore the $c$ columns selected by
the algorithm in each iteration (except the last) will not include any column
from $\Psi$, and they all have to come from $\Phi_{\mathrm{opt}}$. In the last
iteration it is possible to select non-optimal columns from $\Psi$ because all
optimal columns in $\Phi_{\mathrm{opt}}$ have been exhausted. Therefore the
BSR algorithm is guaranteed to locate $k$ optimal columns in $\lceil
k/c\rceil$ iterations. $\square$

####  4.2.2 A Weak Exact Recovery Condition for BSR

What is the condition for $\rho_{c}(\bm{r})<1$ to hold in Theorem 4.1? In the
absence of a straightforward answer, we first discuss the condition for BSR to
recover a $k$-sparse signal within $k$ iterations, then revisit the condition
(3).

We call the condition for BSR to recover a $k$-sparse signal within $k$
iterations the weak exact recovery condition for BSR. For the weak condition,
we use the same greedy choice ratio as defined in section 4.1.1:
$\rho(\bm{r})\stackrel{{\scriptstyle}}{{=}}\frac{\left\|{\Psi}^{\top}\bm{r}\right\|_{\infty}}{\left\|\Phi_{\mathrm{opt}}^{\top}\bm{r}\right\|_{\infty}}$.

######  Theorem 4.2 (The Weak Exact Recovery Condition for BSR).

A sufficient condition for BSR to recover a $k$-sparse signal within $k$
iterations is that

(4) $\rho(\bm{r})<1$

holds for all iterations.

Proof: If $\rho(\bm{r})<1$ holds, then in each iteration BSR can choose at
least one optimal column from $\Phi_{\text{opt}}$. This is because the BSR
algorithm chooses the $c$ columns that give the maximum $\ell_{2}$-norm
$\max\limits_{|\Omega|=c}\left\|{{(\Phi_{\mathrm{opt}}})^{\top}_{\Omega}}\bm{r}\right\|_{2}$.
The $c$ columns that achieved the maximum $\ell_{2}$-norm have the largest $c$
absolute inner products with residual $\bm{r}$, which certainly include the
largest value. If the condition $\rho(\bm{r})<1$ holds in all iterations, the
BSR algorithm is guaranteed to select at least one new column from
$\Phi_{\mathrm{opt}}$. Therefore, for a $k$-sparse signal, BSR can find all
optimal columns in at most $k$ iterations. $\square$

Although intuitive, the condition in (4) expressed in terms of the greedy
choice ratio cannot be checked before we know the residuals in all iterations.
We need to establish a sufficient condition for the exact recovery by BSR in
terms of the property of the dictionary $\Phi$.

BSR will select at least one optimal column at each iteration but can also
possibly select some non-optimal columns. We can split the columns of $\Psi$
into two parts: $\Psi=[\Psi_{J}|\Psi_{\scriptsize\hbox{\vbox{\hrule
height=0.6pt\kern 0.60275pt\hbox{\kern-1.39998pt${J}$\kern 0.70004pt}}}}]$,
where $\Psi_{J}$ are the non-optimal columns that have been selected by BSR
algorithm so far, and $\Psi_{\scriptsize\hbox{\vbox{\hrule height=0.6pt\kern
0.60275pt\hbox{\kern-1.39998pt${J}$\kern 0.70004pt}}}}$ include the remaining
columns.

Let matrix $X$ be the submatrix of $\Phi$ that includes all columns of
$\Phi_{\text{opt}}$ and the columns in $\Psi$ that have been selected by BSR
at the previous iterations, i.e., $X=[\Phi_{\text{opt}}|\Psi_{J}]$.

We first prove a lemma that can lead to the sufficient condition
$\rho(\bm{r})<1$.

Let $\Pi$ denote the index set for the columns in $\Phi_{\mathrm{opt}}$ that
have not been selected by the algorithm so far, so $|\Pi|\leq k$. Let
$(\cdot)_{\Pi}$ denote the columns in the matrix indexed by $\Pi$, and
$(\cdot)_{\Pi,:}$ denotes the rows of the matrix indexed by $\Pi$.

######  Lemma 4.3.

If
$\max\limits_{\bm{\psi}}\left\|\left(X^{+}\right)_{\Pi,:}\bm{\psi}\right\|_{1}<1$,
where vector $\bm{\psi}$ ranges over columns of
$\Psi_{\scriptsize\hbox{\vbox{\hrule height=0.6pt\kern
0.60275pt\hbox{\kern-1.39998pt${J}$\kern 0.70004pt}}}}$, then the residual
$\bm{r}$ satisfies $\rho(\bm{r})<1$.

Proof:

$\begin{aligned}
\rho(\bm{r})&=\frac{\left\|\Psi^{\top}\bm{r}\right\|_{\infty}}{\left\|\Phi_{\mathrm{opt}}^{\top}\bm{r}\right\|_{\infty}}\\\
&\overset{(a)}{=}\frac{\left\|\Psi_{\scriptsize\hbox{\vbox{\hrule
height=0.6pt\kern 0.60275pt\hbox{\kern-1.39998pt${J}$\kern
0.70004pt}}}}^{\top}\bm{r}\right\|_{\infty}}{\left\|\Phi_{\mathrm{opt}}^{\top}\bm{r}\right\|_{\infty}}\\\
&\overset{(b)}{=}\frac{\left\|\Psi_{\scriptsize\hbox{\vbox{\hrule
height=0.6pt\kern 0.60275pt\hbox{\kern-1.39998pt${J}$\kern
0.70004pt}}}}^{\top}X(X^{\top}X)^{-1}X^{\top}\bm{r}\right\|_{\infty}}{\left\|\Phi_{\mathrm{opt}}^{\top}\bm{r}\right\|_{\infty}}\\\
&\overset{(c)}{=}\frac{\left\|\Psi_{\scriptsize\hbox{\vbox{\hrule
height=0.6pt\kern 0.60275pt\hbox{\kern-1.39998pt${J}$\kern
0.70004pt}}}}^{\top}\left(X(X^{\top}X)^{-1}\right)_{\Pi}X_{\Pi}^{\top}\bm{r}\right\|_{\infty}}{\left\|\Phi_{\Pi}^{\top}\bm{r}\right\|_{\infty}}\\\
&\overset{(d)}{\leq}\left\|\Psi_{\scriptsize\hbox{\vbox{\hrule
height=0.6pt\kern 0.60275pt\hbox{\kern-1.39998pt${J}$\kern
0.70004pt}}}}^{\top}\left(\left(X^{+}\right)^{\top}\right)_{\Pi}\right\|_{\infty}\\\
&\overset{(e)}{=}\left\|\left(X^{+}\right)_{\Pi,:}\Psi_{\scriptsize\hbox{\vbox{\hrule
height=0.6pt\kern 0.60275pt\hbox{\kern-1.39998pt${J}$\kern
0.70004pt}}}}\right\|_{1}\\\
&\overset{(f)}{=}\max_{\bm{\psi}}\left\|\left(X^{+}\right)_{\Pi,:}\bm{\psi}\right\|_{1},\\\
\end{aligned}$

where $\bm{\psi}$ ranges over columns of $\Psi_{\scriptsize\hbox{\vbox{\hrule
height=0.6pt\kern 0.60275pt\hbox{\kern-1.39998pt${J}$\kern 0.70004pt}}}}$.

$(a)$ is because
$\left\|\Psi^{\top}\bm{r}\right\|_{\infty}=\left\|\Psi_{\scriptsize\hbox{\vbox{\hrule
height=0.6pt\kern 0.60275pt\hbox{\kern-1.39998pt${J}$\kern
0.70004pt}}}}^{\top}\bm{r}\right\|_{\infty}$, since columns in $\Psi_{J}$ are
orthogonal to $\bm{r}$.

$(b)$ is due to that residual $\bm{r}$ lies in the column span of $X$. Sparse
representation of
$\bm{y}=\sum\limits_{j\in\Lambda_{\mathrm{opt}}}a_{j}\bm{\varphi}_{j}$, so
$\bm{y}$ lies in the column span of $\Phi_{\mathrm{opt}}$. Since $X$ includes
all columns of $\Phi_{\mathrm{opt}}$, $\bm{y}$ also lies in the column span of
$X$. On the other hand, $X$ includes all the columns found by the algorithm so
far, so $\Phi_{\Gamma}\bm{s}_{\Gamma}$ also lies in the column span of $X$.
Therefore the residual $\bm{r}=\bm{y}-\Phi_{\Gamma}\bm{s}_{\Gamma}$ also lies
in the column span of $X$.

$X(X^{\top}X)^{-1}X^{\top}$ is a projector onto the column span of $X$. Since
$\bm{r}$ lies in the column span of $X$, so the projection of $\bm{r}$ onto
the column span of $X$ is $\bm{r}$ itself:

$X(X^{\top}X)^{-1}X^{\top}\bm{r}=\bm{r}.$

In the vectors $X^{\top}\bm{r}$ and $\Phi_{\mathrm{opt}}^{\top}\bm{r}$, all
the nonzero elements are from $\Phi^{\top}_{\Pi}\bm{r}$. Thus,

$\left\|X^{\top}\bm{r}\right\|_{\infty}=\left\|\Phi^{\top}_{\mathrm{opt}}\bm{r}\right\|_{\infty}=\left\|X^{\top}_{\Pi}\bm{r}\right\|_{\infty}=\left\|\Phi^{\top}_{\Pi}\bm{r}\right\|_{\infty}.$

This leads to $(c)$.

$(d)$ is direct from the definition of matrix norm. Suppose that $A$ is a
matrix, and $\bm{v}$ is a vector with nonzero components. For $1\leq
p\leq\infty$, we have

$\|A\|_{p}=\sup\limits_{\bm{v}\neq 0}\frac{\|A\bm{v}\|_{p}}{\|\bm{v}\|_{p}}.$

Therefore,

(5) $\|A\bm{v}\|_{p}\leq\|A\|_{p}\|\bm{v}\|_{p}.$

$(e)$ is due to that the matrix norm $\|\cdot\|_{\infty}$ takes the maximum
absolute row sum, and $\|\cdot\|_{1}$ takes the maximum absolute column sum,
so $\|A^{\top}\|_{\infty}=\|A\|_{1}$ for a given matrix $A$.

$(f)$ follows the definition of matrix norm $\|\cdot\|_{1}$ to get the maximum
absolute column sum of the matrix.

Therefore, if
$\max\limits_{\bm{\psi}}\left\|\left(X^{+}\right)_{\Pi,:}\bm{\psi}\right\|_{1}<1$,
where $\bm{\psi}$ ranges over columns of $\Psi_{\scriptsize\hbox{\vbox{\hrule
height=0.6pt\kern 0.60275pt\hbox{\kern-1.39998pt${J}$\kern 0.70004pt}}}}$,
then $\rho(\bm{r})<1$. $\square$

Although condition
$\max\limits_{\bm{\psi}\in\Psi_{\scriptsize\hbox{\vbox{\hrule
height=0.6pt\kern 0.60275pt\hbox{\kern-1.39998pt${J}$\kern
0.70004pt}}}}}\left\|\left(X^{+}\right)_{\Pi,:}\bm{\psi}\right\|_{1}<1$ is
expressed in terms of the property of the dictionary, this condition still
cannot be checked without executing the algorithm. In practice it is unlikely
that the optimal columns are know a priori, so the submatrices
$X,\Psi_{\scriptsize\hbox{\vbox{\hrule height=0.6pt\kern
0.60275pt\hbox{\kern-1.39998pt${J}$\kern 0.70004pt}}}}$ cannot be located
before the execution of the algorithm. More practical methods are needed to
check the sufficient condition without the execution of the algorithm.

In [21], a fundamental property of the dictionary $\Phi$, called coherence is
defined as:

(6) $\mu\overset{\text{def}}{=}\max\limits_{j\neq
k}|\langle\bm{\varphi}_{j},\bm{\varphi}_{k}\rangle|$

Coherence $\mu$ is the maximum absolute value of pairwise inner product
between the columns of the dictionary.

For a positive integer $m$, the cumulative coherence function, $\mu_{1}(m)$ of
the dictionary, is defined as

(7)
$\mu_{1}(m)\overset{\text{def}}{=}\max\limits_{|\Lambda|=m}\max\limits_{\bm{\psi}}\sum\limits_{j\in\Lambda}|\langle\bm{\varphi}_{j},\bm{\psi}\rangle|$

where $\Lambda$ is the set of indices for any $m$ columns of $\Phi$, and
$\bm{\psi}$ ranges over the columns of $\Phi$ not indexed by $\Lambda$.
$\mu_{1}(m)$ is the maximum cumulative coherence from any $m$ columns of
$\Phi$.

Next we use the cumulative coherence property of the dictionary to derive a
sufficient condition.

######  Lemma 4.4.

$\max\limits_{\bm{\psi}\in\Psi_{\scriptsize\hbox{\vbox{\hrule
height=0.6pt\kern 0.60275pt\hbox{\kern-1.39998pt${J}$\kern
0.70004pt}}}}}\left\|\left(X^{+}\right)_{\Pi,:}\bm{\psi}\right\|_{1}<1$
whenever $\mu_{1}(l)+\mu_{1}(n)<1$ holds, where $n$ is the number of columns
in $X$, and $l=\min(|\Pi|,k-1)$.

Proof:

$\begin{aligned} &\max_{\bm{\psi}\in\Psi_{\scriptsize\hbox{\vbox{\hrule
height=0.6pt\kern 0.60275pt\hbox{\kern-1.39998pt${J}$\kern
0.70004pt}}}}}\left\|\left(X^{+}\right)_{\Pi,:}\bm{\psi}\right\|_{1}\\\
=&\max_{\bm{\psi}\in\Psi_{\scriptsize\hbox{\vbox{\hrule height=0.6pt\kern
0.60275pt\hbox{\kern-1.39998pt${J}$\kern
0.70004pt}}}}}\left\|\left(\left(X^{\top}X\right)^{-1}X^{\top}\right)_{\Pi,:}\bm{\psi}\right\|_{1}\\\
=&\max_{\bm{\psi}\in\Psi_{\scriptsize\hbox{\vbox{\hrule height=0.6pt\kern
0.60275pt\hbox{\kern-1.39998pt${J}$\kern
0.70004pt}}}}}\left\|\left(\left(X^{\top}X\right)^{-1}\right)_{\Pi,:}X^{\top}\bm{\psi}\right\|_{1}\\\
\leq&\left\|\left(\left(X^{\top}X\right)^{-1}\right)_{\Pi,:}\right\|_{1}\cdot\left(\max_{\bm{\psi}\in\Psi_{\scriptsize\hbox{\vbox{\hrule
height=0.6pt\kern 0.60275pt\hbox{\kern-1.39998pt${J}$\kern
0.70004pt}}}}}\left\|X^{\top}\bm{\psi}\right\|_{1}\right)\\\
\leq&\frac{\mu_{1}(n)}{1-\mu_{1}(l)},\\\ \end{aligned}$

where $n$ is the number of columns in $X$, and $l=\min(|\Pi|,k-1)$.

It can be shown that

  1. 1).

$\max\limits_{\bm{\psi}\in\Psi_{\scriptsize\hbox{\vbox{\hrule
height=0.6pt\kern 0.60275pt\hbox{\kern-1.39998pt${J}$\kern
0.70004pt}}}}}\left\|X^{\top}\bm{\psi}\right\|_{1}\leq\mu_{1}(n)$, and

  2. 2).

$\left\|\left(\left(X^{\top}X\right)^{-1}\right)_{\Pi,:}\right\|_{1}\leq\frac{1}{1-\mu_{1}(l)}$.

1) can be proved following the definition of $\mu_{1}(\cdot)$,

$\displaystyle\max_{\bm{\psi}\in\Psi_{\scriptsize\hbox{\vbox{\hrule
height=0.6pt\kern 0.60275pt\hbox{\kern-1.39998pt${J}$\kern
0.70004pt}}}}}\left\|X^{\top}\bm{\psi}\right\|_{1}$
$\displaystyle=\max_{\bm{\psi}\in\Psi_{\scriptsize\hbox{\vbox{\hrule
height=0.6pt\kern 0.60275pt\hbox{\kern-1.39998pt${J}$\kern
0.70004pt}}}}}\sum_{\bm{\varphi_{j}}\in
X}|\langle\bm{\varphi}_{j},\bm{\psi}\rangle|$ $\displaystyle\leq\mu_{1}(n).$

To prove 2), recall that the columns in the dictionary $\Phi$ are normalized
unit vectors (i.e., $\|\bm{\varphi}\|_{2}=1$, for $\bm{\varphi}\in\Phi$). We
can represent

$X^{\top}X=I_{n}+A,$

where $I_{n}$ is an identity matrix of size $n\times n$, and $A$ is a matrix
of size $n\times n$ collecting the pairwise inner products, so
$A_{jk}=\langle\bm{\varphi}_{j},\bm{\varphi}_{k}\rangle,\textrm{ for }j\neq
k$, and the diagonal $A_{jj}=0$, and hence $\|A\|_{1}\leq\mu_{1}(n-1)$
following the definition of the cumulative coherence function
$\mu_{1}(\cdot)$.

When $\|A\|_{1}<1$, the Neumann series
$\sum\limits_{k=0}^{\infty}(-A)^{k}=(I_{n}+A)^{-1}$ (add reference here).
Therefore,
$(X^{\top}X)^{-1}=(I_{n}+A)^{-1}=\sum\limits_{k=0}^{\infty}(-A)^{k}$. Note
that having $\mu_{1}(n-1)<1$ is sufficient to make $\|A\|_{1}<1$.

$\displaystyle\left\|\left(\left(X^{\top}X\right)^{-1}\right)_{\Pi,:}\right\|_{1}$
$\displaystyle=\left\|\left(\sum\limits_{k=0}^{\infty}(-A)^{k}\right)_{\Pi,:}\right\|_{1}$
$\displaystyle\overset{(a)}{\leq}\sum\limits_{k=0}^{\infty}\|A_{\Pi,:}\|^{k}_{1}$
$\displaystyle\overset{(b)}{=}\frac{1}{1-\|A_{\Pi,:}\|_{1}}$
$\displaystyle\overset{(c)}{\leq}\frac{1}{1-\mu_{1}(l)}.$

$(a)$ is due to $\|(-A)^{k}\|_{1}\leq\|A\|_{1}^{k}$; $(b)$ is due to
$\|A_{\Pi,:}\|_{1}\leq\|A\|_{1}<1$; $(c)$ is due to
$\|A_{\Pi,:}\|_{1}\leq\mu_{1}(l)$. At the beginning before the first
iteration, $X=\Phi_{\mathrm{opt}}$, $\Pi=\Lambda_{\mathrm{opt}}$, and $A$ is a
$k\times k$ matrix with zeros at the diagonal positions, so
$\|A_{\Pi,:}\|_{1}\leq\mu_{1}(k-1)$. After the first iteration, $|\Pi|\leq
k-1$. Therefore, $\|A_{\Pi,:}\|_{1}\leq\mu_{1}(l)$ for all iterations, with
$l=\min(|\Pi|,k-1)$.

Therefore $\max\limits_{\bm{\psi}\in\Psi_{\scriptsize\hbox{\vbox{\hrule
height=0.6pt\kern 0.60275pt\hbox{\kern-1.39998pt${J}$\kern
0.70004pt}}}}}\left\|\left(X^{+}\right)_{\Pi,:}\bm{\psi}\right\|_{1}<1$
whenever $\mu_{1}(l)+\mu_{1}(n)<1$. $\square$

Lemma 4.4 and Lemma 4.3 together lead to the following conclusion: the
residual $\bm{r}$ satisfies $\rho(\bm{r})<1$ whenever

(8) $\mu_{1}(l)+\mu_{1}(n)<1.$

####  4.2.3 Revisit Theorem 4.1: the Sufficient Condition for $\rho_{c}({\bf
r})<1$

Now we are ready to show that $\mu_{1}(l)+\mu_{1}(n)<1$ is also sufficient for
$\rho_{c}({\bm{r}})<1$ to hold in Theorem 4.1.

Proof:

From the definition of $\rho_{c}({\bm{r}})$, it is suggested that $\Omega_{1}$
and $\Omega_{2}$ both can have columns from $\Psi$ and $\Phi_{\mathrm{opt}}$.
If there is a column $\bm{\varphi}_{j}$ from $\Phi_{\Omega_{1}}$ and a column
$\bm{\varphi}_{l}$ from $\Phi_{\Omega_{2}}$ such that
$|\langle\bm{\varphi}_{j},\bm{r}\rangle|>|\langle\bm{\varphi}_{l},\bm{r}\rangle|$,
then the greedy algorithm would have selected $\bm{\varphi}_{j}$. In this
case, we would switch the memberships of two columns in $\Omega_{1}$ and
$\Omega_{2}$ as long as such switching does not violate the condition
$|\Omega_{2}\cap\Lambda_{\mathrm{opt}}|>|\Omega_{1}\cap\Lambda_{\mathrm{opt}}|$.
The only case that the switching is forbidden is that $\bm{\varphi}_{j}$ is a
column in $\Psi$ and $\bm{\varphi}_{l}$ is a column in $\Phi_{\mathrm{opt}}$
so that switching would have violated the above condition. There will be at
least one such column in $\Psi$ and one column in $\Phi_{\mathrm{opt}}$ that
cannot be switched.

Let $\Theta_{1}\subseteq\Omega_{1}$ and $\Theta_{2}\subseteq\Omega_{2}$ be the
indices for such non-switchable columns in $\Psi$ and $\Phi_{\mathrm{opt}}$,
respectively. We have
$|\Theta_{1}|=|\Theta_{2}|=|\Omega_{2}\cap\Lambda_{\mathrm{opt}}|-|\Omega_{1}\cap\Lambda_{\mathrm{opt}}|\geq
1$. We define

(9) $\rho^{\prime}_{c}(\bm{r})\stackrel{{\scriptstyle\text{ def
}}}{{\scalebox{2.0}[1.0]{=}}}\frac{\max\limits_{\Theta_{1}}\left\|\Psi^{\top}_{\Theta_{1}}{\bm{r}}\right\|_{2}}{\max\limits_{\Theta_{2}}\left\|(\Phi_{\mathrm{opt}})^{\top}_{\Theta_{2}}\bm{r}\right\|_{2}},$

Note that $\rho^{\prime}_{c}(\bm{r})<1$ is sufficient for
$\rho_{c}(\bm{r})<1$, since for the other columns in
$\Phi_{\Omega_{1}\setminus\Theta_{1}}$, if it happens to have a larger
absolute inner product with $\bm{r}$, we can switch it with a column in
$\Phi_{\Omega_{2}}$ to improve the maximum value obtained by $\Omega_{2}$ in
(2) since switching does not violate the condition
$|\Omega_{2}\cap\Lambda_{\mathrm{opt}}|>|\Omega_{1}\cap\Lambda_{\mathrm{opt}}|$.

Next we show that the sufficient condition for $\rho(\bm{r})<1$ is also
sufficient for $\rho^{\prime}_{c}(\bm{r})<1$.

(10)
$\displaystyle\frac{\max\limits_{\Theta_{1}}\left\|\Psi^{\top}_{\Theta_{1}}{\bm{r}}\right\|_{2}}{\max\limits_{\Theta_{2}}\left\|(\Phi_{\mathrm{opt}})^{\top}_{\Theta_{2}}\bm{r}\right\|_{2}}<1$
$\displaystyle\Longleftrightarrow$
$\displaystyle\frac{\max\limits_{\Theta_{1}}\sum\limits_{{\bm{\psi}}_{j}\in\Psi_{\Theta_{1}}}\left\|{\bm{\psi}}_{j}^{\top}{\bm{r}}\right\|_{2}^{2}}{\max\limits_{\Theta_{2}}\sum\limits_{{\bm{\varphi}}_{l}\in(\Phi_{\mathrm{opt}})_{\Theta_{2}}}\left\|{\bm{\varphi}}_{l}^{\top}\bm{r}\right\|_{2}^{2}}<1\;.$

Now we sort in descending order the absolute inner products
$|\langle\bm{\psi}_{j},\bm{r}\rangle|$ and
$|\langle\bm{\varphi}_{l},\bm{r}\rangle|$, respectively. Let the index set
$\Theta_{1}=\\{j_{(1)},j_{(2)},j_{(3)},\ldots\\}$, and
$\Theta_{2}=\\{l_{(1)},l_{(2)},l_{(3)},\ldots\\}$ with subscript $(i)$
indicating the $i$-th largest value. It is sufficient to show that
$\rho^{\prime}_{c}(\bm{r})<1$ as long as each of the following
$\rho^{\prime}_{(i)}(\bm{r})<1$, which is defined as

(11) $\rho^{\prime}_{(i)}(\bm{r})\stackrel{{\scriptstyle\text{ def
}}}{{\scalebox{2.0}[1.0]{=}}}\frac{|\langle{\bm{\psi}}_{j_{(i)}},{\bm{r}}\rangle|}{|\langle{\bm{\varphi}}_{l_{(i)}},{\bm{r}}\rangle|},$

where $i=1,\ldots,|\Theta_{1}|$.

If $\rho(\bm{r})<1$ holds for the current iteration, then
$\rho^{\prime}_{(1)}(\bm{r})<1$ also holds, since
$\rho(\bm{r})=\rho^{\prime}_{(1)}(\bm{r})$. After removing $j_{(1)}$ from
$\Theta_{1}$ and $l_{(1)}$ from $\Theta_{2}$, the the second largest becomes
the largest, then we are facing the same problem again:

$\rho^{\prime}_{(2)}(\bm{r})=\frac{\max_{\bm{\psi}}|\langle\bm{\psi},\bm{r}\rangle|}{\max_{\bm{\varphi}}|\langle\bm{\varphi},\bm{r}\rangle|},$

where $\bm{\psi}$ ranges over columns in $\Psi$ excluding
$\bm{\psi}_{j_{(1)}}$, and $\bm{\varphi}$ ranges over columns in
$\Phi_{\mathrm{opt}}$ excluding $\bm{\varphi}_{j_{(1)}}$.

If condition $\mu_{1}(l)+\mu_{1}(n)<1$ is satisfied, then
$\rho^{\prime}_{(2)}(\bm{r})<1$ also holds, and so on. Since
$\rho^{\prime}_{(i)}(\bm{r})<1$ holds for all $i$ until exhausting all indices
in $\Theta_{1}$ and $\Theta_{2}$, and each term in the summation in (10) is
nonnegative, then $\rho^{\prime}_{c}(\bm{r})<1$ holds. Therefore
$\mu_{1}(l)+\mu_{1}(n)<1$ is also sufficient for $\rho_{c}(\bm{r})<1$.

#### Upper bound of $n+l$

It is trivial to show that $\mu_{1}(l)+\mu_{1}(n)\leq\mu(n+l)$. Therefore, it
is necessary to find the upper bound of $n+l$.

Recall that $n$ is the number of columns in $X$, and $l=\min(|\Pi|,k-1)$,
where $\Pi$ denotes the index set for the columns in $\Phi_{\mathrm{opt}}$
that have not been selected by the algorithm so far.

Before the first iteration, $n=k$, and $l=k-1$, $n+l=2k-1$. If
$\mu_{1}(l)+\mu_{1}(n)<1$ holds, then $\rho_{c}(\bm{r})<1$ holds, so the
greedy algorithm would select as many optimal columns as possible. As long as
$|\Pi|\geq c$, it would select $c$ optimal columns. Assume $k\geq c\geq 1$.
After the first iteration, $n=k$ and $l=k-c$, therefore $n+l\leq 2k-1$. In
subsequent iterations, $l$ decreases faster than $n$ increases, thus $n+l$ is
not increasing. Therefore the upper bound for $n+l$ is

(12) $n+l\leq 2k-1.$

Combining Theorem 4.1 and upper bound (12), we have the following theorem:

######  Theorem 4.5 (The strong exact recovery condition for BSR).

Suppose that $\mu$ is the coherence of the dictionary as defined in (6). A
sufficient condition for BSR to recover a $k$-sparse signal within $\lceil
k/c\rceil$ iterations is that

(13) $\mu(2k-1)<1.$

##  5 Approximation to General Signals

General signals may not have an exact sparse representation. In this case, the
best we can obtain is its approximation. Suppose that for an arbitrary input
signal $\bm{y}$, there exists a $k$-term optimal sparse approximation
$\bm{s}_{\mathrm{opt}}$, i.e.,
$\bm{y}=\bm{s}_{\mathrm{opt}}+\bm{\varepsilon},\text{ and
}\bm{s}_{\mathrm{opt}}=\sum\limits_{j\in\Lambda_{\mathrm{opt}}}a_{j}\bm{\varphi}_{j},$
where $|\Lambda_{\mathrm{opt}}|=k$. The difference between $\bm{y}$ and its
optimal approximation $\bm{s}_{\mathrm{opt}}$ is the irreducible error
$\bm{\varepsilon}$. We can analyze the convergence property of the BSR
algorithm from the perspectives of a sufficient condition under which the
algorithm can find all optimal atoms within $k$ iterations and the
approximation error when the algorithm terminates. The approximation error
will be measured by $\ell_{2}$ norm.

###  5.1 Recovery Condition for General Signals

Suppose at the current iteration, the algorithm found a sparse representation
$\bm{s}=\sum\limits_{\bm{\varphi}_{j}\in\Phi}b_{j}\bm{\varphi}_{j}$, then the
residual is $\bm{r}=\bm{y}-\bm{s}$. The algorithm can find at least one
optimal column in $\Lambda_{\mathrm{opt}}$ in the next iteration if the greedy
choice ratio $\rho(\bm{r})<1$.

$\begin{aligned}
\rho(\bm{r})&=\frac{\left\|\Psi^{\top}\bm{r}\right\|_{\infty}}{\left\|\Phi_{\mathrm{opt}}^{\top}\bm{r}\right\|_{\infty}}\\\
&=\frac{\left\|\Psi^{\top}(\bm{y}-\bm{s})\right\|_{\infty}}{\left\|\Phi_{\mathrm{opt}}^{\top}(\bm{y}-\bm{s})\right\|_{\infty}}\\\
&=\frac{\left\|\Psi^{\top}(\bm{y}-\bm{s}_{\mathrm{opt}}+\bm{s}_{\mathrm{opt}}-\bm{s})\right\|_{\infty}}{\left\|\Phi_{\mathrm{opt}}^{\top}(\bm{y}-\bm{s}_{\mathrm{opt}}+\bm{s}_{\mathrm{opt}}-\bm{s})\right\|_{\infty}}\\\
&\overset{(a)}{=}\frac{\left\|\Psi^{\top}(\bm{y}-\bm{s}_{\mathrm{opt}}+\bm{s}_{\mathrm{opt}}-\bm{s})\right\|_{\infty}}{\left\|\Phi_{\mathrm{opt}}^{\top}(\bm{s}_{\mathrm{opt}}-\bm{s})\right\|_{\infty}}\\\
&\leq\frac{\left\|\Psi^{\top}(\bm{y}-\bm{s}_{\mathrm{opt}})\right\|_{\infty}}{\left\|\Phi_{\mathrm{opt}}^{\top}(\bm{s}_{\mathrm{opt}}-\bm{s})\right\|_{\infty}}+\frac{\left\|\Psi^{\top}(\bm{s}_{\mathrm{opt}}-\bm{s})\right\|_{\infty}}{\left\|\Phi_{\mathrm{opt}}^{\top}(\bm{s}_{\mathrm{opt}}-\bm{s})\right\|_{\infty}}\\\
&=\rho_{\bm{\varepsilon}}+\rho_{\mathrm{t}},\\\ \end{aligned}$

where (a) is due to that the columns in $\Phi_{\mathrm{opt}}$ are orthogonal
to the vector $\bm{y}-\bm{s}_{\mathrm{opt}}$. Therefore
$\Phi_{\mathrm{opt}}^{\top}(\bm{y}-\bm{s}_{\mathrm{opt}})$ yields an all-zero
vector.

The second term $\rho_{\mathrm{t}}$ is caused by the gap between what the
algorithm has found so far and the optimal approximation. Obviously
$\rho_{\mathrm{t}}$ is the same $\rho$ we would get for a signal that has an
exact sparse representation, i.e., when $\bm{y}=\bm{s}_{\mathrm{opt}}$,
therefore $\rho_{\mathrm{t}}\leq\frac{\mu_{1}(n)}{1-\mu_{1}(l)}$ from (4.3)
and (4.4).

The first term $\rho_{\bm{\varepsilon}}$ is caused by the irreducible error
existing in the signal. $\rho_{\bm{\varepsilon}}$ can also be bounded.

$\begin{aligned}
\rho_{\bm{\varepsilon}}&=\frac{\left\|\Psi^{\top}(\bm{y}-\bm{s}_{\mathrm{opt}})\right\|_{\infty}}{\left\|\Phi_{\mathrm{opt}}^{\top}(\bm{s}_{\mathrm{opt}}-\bm{s})\right\|_{\infty}}\\\
&\overset{(a)}{\leq}\frac{\left\|\bm{y}-\bm{s}_{\mathrm{opt}}\right\|_{2}}{\left\|\Phi_{\mathrm{opt}}^{\top}(\bm{s}_{\mathrm{opt}}-\bm{s})\right\|_{\infty}}\\\
&\overset{(b)}{\leq}\frac{\sqrt{k}\left\|\bm{y}-\bm{s}_{\mathrm{opt}}\right\|_{2}}{\left\|\Phi_{\mathrm{opt}}^{\top}(\bm{s}_{\mathrm{opt}}-\bm{s})\right\|_{2}}\\\
&\overset{(c)}{\leq}\frac{\sqrt{k}\left\|\bm{y}-\bm{s}_{\mathrm{opt}}\right\|_{2}}{\sigma_{k}(\Phi_{\mathrm{opt}})\left\|\bm{s}_{\mathrm{opt}}-\bm{s}\right\|_{2}}\\\
&\overset{(d)}{\leq}\frac{\sqrt{k}\left\|\bm{y}-\bm{s}_{\mathrm{opt}}\right\|_{2}}{\sqrt{1-\mu_{1}(k-1)}\left\|\bm{s}_{\mathrm{opt}}-\bm{s}\right\|_{2}}\\\
\end{aligned}$

(a) is due to that
$\left\|\Psi^{\top}(\bm{y}-\bm{s}_{\mathrm{opt}})\right\|_{\infty}\leq\left\|\bm{y}-\bm{s}_{\mathrm{opt}}\right\|_{2}$:

$\displaystyle\left\|\Psi^{\top}(\bm{y}-\bm{s}_{\mathrm{opt}})\right\|_{\infty}$
$\displaystyle=\max\limits_{\bm{\psi}\in\Psi}|\bm{\psi}^{\top}(\bm{y}-\bm{s}_{\mathrm{opt}})|$
$\displaystyle\leq\max\limits_{\bm{\psi}\in\Psi}\|\bm{\psi}\|_{2}\left\|\bm{y}-\bm{s}_{\mathrm{opt}}\right\|_{2}.$
$\displaystyle=\left\|\bm{y}-\bm{s}_{\mathrm{opt}}\right\|_{2}$

$\|\bm{\psi}\|_{2}=1$ since columns in $\Phi$ are normalized.

(b) is because for a $k\times 1$ vector $\bm{v}$,
$\|\bm{v}\|_{2}\leq\sqrt{k}\|\bm{v}\|_{\infty}$. Therefore
$\left\|\Phi_{\mathrm{opt}}^{\top}(\bm{s}_{\mathrm{opt}}-\bm{s})\right\|_{2}\leq\sqrt{k}\left\|\Phi_{\mathrm{opt}}^{\top}(\bm{s}_{\mathrm{opt}}-\bm{s})\right\|_{\infty}$.

(c) is due to the singular value definition. $\Phi_{\mathrm{opt}}$ has $k$
columns. The $k$-th singular value of $\Phi_{\mathrm{opt}}$ is defined as

$\sigma_{k}(\Phi_{\mathrm{opt}})=\inf\limits_{\bm{u}\neq
0}\frac{\|\Phi_{\mathrm{opt}}^{\top}\bm{u}\|_{2}}{\|\bm{u}\|_{2}},$

and therefore,

(14)
$\sigma_{k}(\Phi_{\mathrm{opt}})\|\bm{s}_{\mathrm{opt}}-\bm{s}\|_{2}\leq\|\Phi_{\mathrm{opt}}^{\top}(\bm{s}_{\mathrm{opt}}-\bm{s})\|_{2}.$

(d) is due to the relation of singular values to eigenvalues. Let matrix
$Z=\Phi_{\mathrm{opt}}^{\top}\Phi_{\mathrm{opt}}$, then
$\lambda_{k}(Z)=\sigma_{k}^{2}(\Phi_{\mathrm{opt}})$. By the Gershgorin circle
theorem, every eigenvalue of the square matrix $Z$ lies within at least one of
the Gershgorin discs, i.e.,

$\begin{aligned} |\lambda-Z_{ii}|\leq\sum\limits_{j\neq
i}|Z_{ij}|=\sum\limits_{j\neq
i}|\langle\bm{\varphi}_{i},\bm{\varphi}_{j}\rangle|\leq\mu_{1}(k-1),\end{aligned}$

and $Z_{ii}=1$ since columns in $\Phi$ are normalized. Thus,
$-\mu_{1}(k-1)\leq 1-\lambda_{k}\leq\mu_{1}(k-1)$.

Therefore, $\sigma_{k}(\Phi_{\mathrm{opt}})\geq\sqrt{1-\mu_{1}(k-1)}$.
$\square$

Combining $\rho_{\bm{\varepsilon}}$ and $\rho_{\mathrm{t}}$, we have

(15)
$\rho(\bm{r})\leq\frac{\sqrt{k}\left\|\bm{y}-\bm{s}_{\mathrm{opt}}\right\|_{2}}{\sqrt{1-\mu_{1}(k-1)}\left\|\bm{s}_{\mathrm{opt}}-\bm{s}\right\|_{2}}+\frac{\mu_{1}(n)}{1-\mu_{1}(l)}.$

######  Theorem 5.1 (General Recovery Condition of BSR).

The BSR algorithm can find all optimal atoms from $\Phi_{\mathrm{opt}}$ within
$k$ iterations provided that the following condition holds for all iterations:

(16)
$\frac{\sqrt{k}\left\|\bm{y}-\bm{s}_{\mathrm{opt}}\right\|_{2}}{\sqrt{1-\mu_{1}(k-1)}\left\|\bm{s}_{\mathrm{opt}}-\bm{s}\right\|_{2}}+\frac{\mu_{1}(n)}{1-\mu_{1}(l)}<1.$

If the irreducible error $\|\bm{y}-\bm{s}_{\mathrm{opt}}\|_{2}$ is nonzero, it
imposes a more stringent condition on the dictionary $\Phi$ regarding the
orthogonality between atoms, as it requires the coherence $\mu$ to be very
small to satisfy condition (16). Compared to signals that have exact sparse
representations, the general signals are more difficult to recover. Noise
works adversely for BSR the same way it does for OMP and OMP-SR.

###  5.2 Approximation Error of BSR

Let $\bm{s}_{\kappa}$ represent the approximation solution found by the BSR
algorithm, either after the fixed number of iterations, or after it breaks out
from the for loop due to the residual $\|\bm{r}\|_{2}$ being too small. In the
latter case, the residual is already bounded by $\delta$ given in the stopping
criterion. We focus on the first case when the BSR algorithm has finished
$\kappa$ iterations, and ask what the upper bound for
$\|\bm{y}-\bm{s}_{\kappa}\|_{2}$ will be. If the algorithm has found all $k$
optimal atoms when it terminates, the approximation error is bounded by
$\|\bm{y}-\bm{s}_{\mathrm{opt}}\|_{2}$. It is only worth discussing when the
algorithm has not found all optimal atoms after $\kappa$ iterations.

We first consider a baseline solution that has only used the optimal atoms
from $\Phi_{\mathrm{opt}}$ so far, so that
$\bm{s}_{\kappa}=\sum\limits_{j\in\Lambda_{\mathrm{opt}}}b_{j}\bm{\varphi}_{j}$.
From (16) we know that as long as the following condition holds, the algorithm
will find at least one optimal atom in the next iteration:

$\frac{\sqrt{k}\;\|\bm{y}-\bm{s}_{\mathrm{opt}}\|_{2}}{\sqrt{1-\mu_{1}(k-1)}\;\|\bm{s}_{\mathrm{opt}}-\bm{s}\|_{2}}<\frac{1-\mu_{1}(l)-\mu_{1}(n)}{1-\mu_{1}(l)}$

We need $\mu_{1}(l)+\mu_{1}(n)<1$ to ensure that the right hand side of the
inequality is positive.

The above inequality is equivalent to

$\begin{aligned}
\|\bm{s}_{\mathrm{opt}}-\bm{s}\|_{2}&>\frac{\sqrt{k}\left(1-\mu_{1}(l)\right)\;\|\bm{y}-\bm{s}_{\mathrm{opt}}\|_{2}}{\sqrt{1-\mu_{1}(k-1)}(1-\mu_{1}(l)-\mu_{1}(n))}\\\
&\overset{(a)}{\geq}\frac{\sqrt{k(1-\mu_{1}(k-1))}}{1-\mu_{1}(l)-\mu_{1}(n)}\|\bm{y}-\bm{s}_{\mathrm{opt}}\|_{2}\\\
&=h\|\bm{y}-\bm{s}_{\mathrm{opt}}\|_{2}\\\ \end{aligned}$

where $h=\frac{\sqrt{k(1-\mu_{1}(k-1))}}{1-\mu_{1}(l)-\mu_{1}(n)}$.

(a) is due to that $l\leq k-1$, therefore $\mu_{1}(l)\leq\mu_{1}(k-1)$ based
on the definition in (7), therefore $1-\mu_{1}(l)\geq 1-\mu_{1}(k-1)$.

If the algorithm fails to find an optimal atom in the next iteration, it is
because the above condition is broken, then we have

(17) $\|\bm{s}_{\mathrm{opt}}-\bm{s}\|_{2}\leq
h\|\bm{y}-\bm{s}_{\mathrm{opt}}\|_{2}.\\\ $

Let $\bm{u}=\bm{y}-\bm{s}_{\mathrm{opt}}$ and
$\bm{v}=\bm{s}_{\mathrm{opt}}-\bm{s}$. Then we have $\|\bm{v}\|_{2}\leq
h\|\bm{u}\|_{2}$ and $\bm{y}-\bm{s}=\bm{u}+\bm{v}$. In addition, $\bm{u}$ and
$\bm{v}$ are orthogonal to each other since $\bm{u}$ is the residual of
projecting $\bm{y}$ onto the the column span of $\Phi_{\mathrm{opt}}$ and
$\bm{v}$ is a vector in the column span of $\Phi_{\mathrm{opt}}$. From the
Pythagorean Theorem that says if $\bm{u}\perp\bm{v}$, then
$\|\bm{u}+\bm{v}\|_{2}=\sqrt{\|\bm{u}\|_{2}^{2}+\|\bm{v}\|_{2}^{2}}$. Thus,

(18)
$\|\bm{y}-\bm{s}\|_{2}\leq\sqrt{1+h^{2}}\|\bm{y}-\bm{s}_{\mathrm{opt}}\|_{2}.\\\
$

This means if the algorithm fails to find an optimal atom in the next
iteration, the current residual is upper bounded by
$\sqrt{1+h^{2}}\|\bm{y}-\bm{s}_{\mathrm{opt}}\|_{2}$. The future iterations
will not increase the residual due to the least square principle.

So far we have established the error bound for a baseline solution that only
uses atoms from $\Phi_{\mathrm{opt}}$ up to the current iteration. If BSR
algorithm has used non-optimal columns in addition to those optimal atoms, the
least square estimate will only produce smaller residual due to the least
square principle, i.e.,

$\|\bm{y}-\bm{s}_{\kappa}\|_{2}\leq\|\bm{y}-\bm{s}\|_{2}\leq\sqrt{1+h^{2}}\|\bm{y}-\bm{s}_{\mathrm{opt}}\|_{2}.$

Therefore we can conclude that by the time the BSR algorithm has finished, the
error bound is upper bounded by
$\sqrt{1+h^{2}}\|\bm{y}-\bm{s}_{\mathrm{opt}}\|_{2}$.

######  Theorem 5.2 (The Approximation Error of BSR).

Suppose that $\mu_{1}(l)+\mu_{1}(n)<1$. Let
$h=\frac{\sqrt{k(1-\mu_{1}(k-1))}}{1-\mu_{1}(l)-\mu_{1}(n)}$. The residual
from the approximation solution $\bm{s}_{\kappa}$ found by the BSR algorithm
satisfies

(19)
$\|\bm{y}-\bm{s}_{\kappa}\|_{2}\leq\sqrt{1+h^{2}}\|\bm{y}-\bm{s}_{\mathrm{opt}}\|_{2}\\\
$

when the algorithm fails to find all optimal atoms, and

(20)
$\|\bm{y}-\bm{s}_{\kappa}\|_{2}\leq\|\bm{y}-\bm{s}_{\mathrm{opt}}\|_{2}\\\ $

when the algorithm succeeds in finding all optimal atoms.

##  6 Experiments

Data used in the experiments are posted at github.

|  |  |
---|---|---|---
|  |  |

Figure 1: The images used for experiment, (a) phantom, (b) transaxial CT, (c)
trees, (d) letters, (e)–(h) MNIST dataset handwritten digits.

###  6.1 Sparse Signal Recovery

We first show sparse signal recovery performance when the signal has a sparse
representation. The first experiment is on image data, where the non-zero
elements constitute the content of an image, and exhibit continuity in the
true signal $\bm{x}$. The images are shown in Figure 1. The second experiment
is on signals defined on graph structures, where the non-zero elements are
distributed among the nodes of a graph. We used synthetic data defined on a
binary tree, and data that are collected from IEEE 118-bus power system and
IEEE 1354-bus power system, where the true signal $\bm{x}$ consists of the
values of the state variables of a power system. Since the algorithms do not
depend on the signal structure to find the non-zeros, they worked well with
both types of data.

Table 1 and Table 2 show the performance of the proposed OMP-SR and BSR, and
we report the number of iterations, the recovered optimal atoms, normalized
MSE (NMSE), and running time in seconds. Image data are reported in Table 1,
and graph data are reported in Table 2.

If a $k$-sparse signal can be recovered by OMP within $k$ iterations, it can
certainly be recovered by BSR within $k$ iterations. Those that cannot be
recovered by OMP within $k$ iterations are shown to take far less than $k$
iterations and far less time by the BSR algorithm to fully recover. Since OMP-
SR picks the same atoms as OMP does, we reported the result of OMP-SR in the
same row as OMP and only reported its time (in blue text). It is observed that
OMP-SR is faster than OMP, and BSR is faster than OMP and OMP-SR in every
single test. The block size $c$ in BSR is a hyper parameter searched from {2,
3, 4, 8}.

The third experiment is to show the relation between $k$ and running time.
Image data for the phantom and the MNIST handwritten digit ’7’ were used. We
created different versions from the original image to have different image
sizes $d$ and different $k/d$ ratios. Figure 2 shows how the running times and
iterations of OMP, OMP-SR, and BSR increase as the number of non-zeros $k$
increases. As $k$ increases, BSR (blue curve) shows significant improvement in
running time over OMP (red curve), while OMP-SR is in-between.

|
---|---
|

Figure 2: Running time (top) and iterations (bottom) used by the algorithms to
recover $k$ non-zeros in the signal. OMP and OMP-SR use the same number of
iterations. The datasets used: (a),(c) phantom; (b),(d) MNIST dataset
handwritten digit ’7’.

###  6.2 Sparse Approximation

The fourth experiment is for the general signal’s sparse approximation. We
added noise $\bm{\varepsilon}$ to the Phantom image, and reported the
approximation error when the measurements are subject to increasing levels of
noise. Table 3 shows that at each noise level, BSR found the $k$ non-zeros
with fewer iterations than OMP and significantly less running time.

Table 1: Image datasets. Reported NMSE and time in seconds. OMP-SR running
time is highlighted in blue.

Data | k | method  | ite | found  | NMSE | time
---|---|---|---|---|---|---
MNIST(3) $392\\!\\!\times\\!\\!784\\!\\!$ | 126 | OMP | 126 | 110 | $0.0656$ | 0.4077 (0.2977)
OMP | 142 | 126 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 0.5799  (0.3280)
BSR | 51 | 126 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 0.1298
MNIST(5) $392\\!\\!\times\\!\\!784\\!\\!$ | 162 | OMP | 162 | 99 | $0.7058$ | 0.5902 (0.4204)
OMP | 784 | 162 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 18.0441  (12.993832)
BSR | 85 | 162 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 0.3926
MNIST(8) $392\\!\\!\times\\!\\!784\\!\\!$ | 174 | OMP | 174 | 92 | $0.8275$ | 0.7703 (0.4934)
OMP | 546 | 174 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 11.0087  (6.5816)
BSR | 84 | 174 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 0.3855
MNIST(9) ${392\\!\\!\times\\!\\!784}$ | 130 | OMP | 130 | 125 | 0.0727 | 0.4453  (0.3190)
OMP | 135 | 130 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 0.5620 (0.3917)
BSR | 35 | 130 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 0.1111
Phantom ${4512\\!\\!\times\\!\\!9024}$ | 641 | OMP | 641 | 638 | 0.0278 | 89.0610  (70.8633)
OMP | 644 | 641 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 97.9924 (75.9530)
BSR | 81 | 641 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 9.1253
Transaxial CT   ${4225\\!\\!\times\\!\\!8450}$ |  1089  | OMP | 1089 | 1064 | $0.0675$ | 282.6071 (250.1557)
OMP | 1115 | 1089 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 302.0011 (263.0615)
BSR | 57 | 1089 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 14.0619
Trees ${19200\\!\\!\times\\!\\!38400\\!\\!}$ |  4670  | OMP | 4670 | 4652 |  $\\!\\!0.00114$  | 2391.0661  (702.4279)
OMP | 4688 | 4670 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 2571.1333  (754.8234)
BSR | 117 | 4670 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 20.0294
Letters ${5712\\!\\!\times\\!\\!11424}$ | 851 | OMP | 851 | 851 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 191.7733 (129.9197)
BSR | 107 | 851 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 16.4433

Table 2: Synthetic data for signals defined on graph structures. Reported NMSE
and time in seconds. OMP-SR running time is highlighted in blue.

Data | k | method | ite | found | NMSE | time
---|---|---|---|---|---|---
Binary Tree ${256\\!\\!\times\\!\\!512}$ | 70 | OMP | 70 | 70 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 0.1105 (0.0798)
BSR | 25 | 70 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 0.024294
118 Bus ${59\\!\\!\times\\!\\!118}$ | 14 | OMP | 14 | 14 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 0.0284 (0.0172)
BSR | 5 | 14 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 0.003876
118 Bus ${59\\!\\!\times\\!\\!118}$ | 100 | OMP | 100 | 86 | 0.8406 | 0.0706 (0.0582)
OMP | 118 | 100 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 0.0922 (0.0659)
BSR | 40 | 100 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 0.025842
1354 Bus ${677\\!\\!\times\\!\\!1354}$ | 270 | OMP | 270 | 215 | $0.2029$ | 3.7723 (2.1686)
OMP | 336 | 270 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 5.4101 (3.4590)
BSR | 96 | 270 | $<\\!\\!1e\\!\\!-\\!\\!11$ | 0.982119

Table 3: Results for the phantom image with increasing noise level
$\|\bm{\varepsilon}\|_{2}$. Reported normalized approximation error
$\frac{\|\bm{y}-\Phi\bm{x}\|_{2}}{\|\bm{y}\|_{2}}$, and running time in
seconds. OMP-SR running time is highlighted in blue.

Noise | k | method  | ite | found |  $\frac{\|\bm{y}-\Phi\bm{x}\|_{2}}{\|\bm{y}\|_{2}}$  | time
---|---|---|---|---|---|---
$\|\bm{\varepsilon}\|_{2}=0.1$ | 641 | OMP | 641 | 638 | 0.0231 | 88.7496  (70.3730)
OMP | 644 | 641 | $0.0001$ | 94.6122 (72.1749)
BSR | 81 | 641 | $0.0001$ | 9.1572
$\|\bm{\varepsilon}\|_{2}\\!=\\!50$ | 641 | OMP | 641 | 638 | 0.0777 | 87.1276  (70.6804)
OMP | 644 | 641 | $0.0727$ | 95.4393 (71.5545)
BSR | 81 | 641 | $0.0720$ | 9.7602
$\|\bm{\varepsilon}\|_{2}\\!=\\!100$ | 641 | OMP | 641 | 631 | 0.1470 | 88.9356  (72.0881)
OMP | 684 | 641 | $0.1337$ | 165.5339 (138.7913)
BSR | 101 | 641 | $0.1287$ | 13.8060
$\|\bm{\varepsilon}\|_{2}\\!=\\!150$ | 641 | OMP | 641 | 621 | 0.2687 | 89.2221  (70.3730)
OMP | 789 | 641 | $0.1671$ | 134.5601  (101.7830)
BSR | 197 | 641 | $0.1671$ | 25.4323

##  7 Discussion and Future Work

OMP has the advantage of simplicity. A greedy algorithm such as OMP is easy to
implement but difficult to analyze. This work offered significant performance
improvement over the classical OMP with theoretical analysis for convergence
and approximation error bound. In addition, the proposed extensions for OMP
come from a principled approach. They work well when combined with other
heuristic or ensemble approaches. One possible future work direction is to
improve the greedy choice by leveraging the structure in the signal model.

The iterative approach of OMP (including ours) takes more iterations when $k$
increases, and the cost of orthogonal projection also increases. Our methods
reduced the cost of orthogonal projection in each iteration, but the total
cost still increases with iterations. The minimal $\ell_{1}$ norm solution is
the sparsest only when the signal is sparse enough [11]. Therefore another
future work direction is to identify the specific measurement matrix property
that drives sparsity during $\ell_{1}$ norm minimization and use that to
improve the greedy choice in an iterative procedure.

## References

  * [1] A. S. Bandeira, E. Dobriban, D. G. Mixon, and W. F. Sawin, Certifying the restricted isometry property is hard, IEEE Transactions on Information Theory, 59 (2013), pp. 3448–3450.
  * [2] R. Baraniuk, M. Davenport, R. DeVore, and M. Wakin, A simple proof of the restricted isometry property for random matrices, Constructive Approximation, 28 (2008), pp. 253–263, https://doi.org/10.1007/s00365-007-9003-x, https://doi.org/10.1007/s00365-007-9003-x.
  * [3] E. J. Candes, The restricted isometry property and its implications for compressed sensing, Comptes rendus mathematique, 346 (2008), pp. 589–592.
  * [4] E. J. Candes and T. Tao, Near-optimal signal recovery from random projections: Universal encoding strategies?, IEEE Transactions on Information Theory, 52 (2006), pp. 5406–5425.
  * [5] L.-H. Chang and J.-Y. Wu, An improved rip-based performance guarantee for sparse signal recovery via orthogonal matching pursuit, IEEE Transactions on Information Theory, 60 (2014), pp. 5702–5715.
  * [6] S. Chen and D. Donoho, Basis pursuit, in Proceedings of 1994 28th Asilomar Conference on Signals, Systems and Computers, vol. 1, 1994, pp. 41–44.
  * [7] S. S. Chen, D. L. Donoho, and M. A. Saunders, Atomic decomposition by basis pursuit, SIAM Journal on Scientific Computing, 20 (1998), pp. 33–61.
  * [8] W. Dai and O. Milenkovic, Subspace pursuit for compressive sensing signal reconstruction, IEEE Transactions on Information Theory, 55 (2009), pp. 2230–2249, https://doi.org/10.1109/TIT.2009.2016006.
  * [9] I. Daubechies, Time-frequency localization operators: a geometric phase space approach, IEEE Transactions on Information Theory, 34 (1988), pp. 605–612.
  * [10] S. Dirksen, G. LecuÃ©, and H. Rauhut, On the gap between restricted isometry properties and sparse recovery conditions, IEEE Transactions on Information Theory, 64 (2018), pp. 5478–5487.
  * [11] D. L. Donoho, For most large underdetermined systems of linear equations the minimal $\ell_{1}$-norm solution is also the sparsest solution, Communications on Pure and Applied Mathematics, 59 (2006), pp. 797–829.
  * [12] D. L. Donoho and M. Elad, Optimally sparse representation in general (nonorthogonal) dictionaries via $\ell^{1}$ minimization, Proceedings of the National Academy of Sciences, 100 (2003), pp. 2197–2202.
  * [13] D. L. Donoho, Y. Tsaig, I. Drori, and J. Starck, Sparse solution of underdetermined systems of linear equations by stagewise orthogonal matching pursuit, IEEE Transactions on Information Theory, 58 (2012), pp. 1094–1121, https://doi.org/10.1109/TIT.2011.2173241.
  * [14] H. Ge and W. Chen, Recovery of signals by a weighted $\ell_{2}\/\ell_{1}$ minimization under arbitrary prior support information, Signal Processing, 148 (2018), pp. 288–302.
  * [15] S. Kwon, J. Wang, and B. Shim, Multipath matching pursuit, IEEE Transactions on Information Theory, 60 (2014), pp. 2986–3001, https://doi.org/10.1109/TIT.2014.2310482.
  * [16] S. Mallat and Z. Zhang, Matching pursuits with time-frequency dictionaries, IEEE Transactions on Signal Processing, 41 (1993), pp. 3397–3415.
  * [17] D. Needell and J. A. Tropp, Cosamp: Iterative signal recovery from incomplete and inaccurate samples, Applied and computational harmonic analysis, 26 (2009), pp. 301–321.
  * [18] D. Needell and R. Vershynin, Uniform uncertainty principle and signal recovery via regularized orthogonal matching pursuit, Foundations of Computational Mathematics, 9 (2009), pp. 317–334.
  * [19] Y. Pati, R. Rezaiifar, and P. Krishnaprasad, Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition, in Proceedings of 27th Asilomar Conference on Signals, Systems and Computers, vol. 1, 1993, pp. 40–44.
  * [20] A. M. Tillmann and M. E. Pfetsch, The computational complexity of the restricted isometry property, the nullspace property, and related concepts in compressed sensing, IEEE Transactions on Information Theory, 60 (2013), pp. 1248–1259.
  * [21] J. A. Tropp, Greed is good: Algorithmic results for sparse approximation, IEEE Transactions on Information theory, 50 (2004), pp. 2231–2242.
  * [22] J. A. Tropp and A. C. Gilbert, Signal recovery from random measurements via orthogonal matching pursuit, IEEE Transactions on information theory, 53 (2007), pp. 4655–4666.
  * [23] P. Yin, Y. Lou, Q. He, and J. Xin, Minimization of $\ell_{1-2}$ for compressed sensing, SIAM Journal on Scientific Computing, 37 (2015), pp. A536–A563.

## Appendix: OMP Algorithm

The OMP algorithm from [19] is listed here for reference

Algorithm 3 OMP [19]

Input: $\Phi\in\mathbb{R}^{N\times d}$, $\bm{y}\in\mathbb{R}^{N}$, $\kappa$

Initialization:  $\bm{r}^{0}=\bm{y}$, $\Gamma=J^{0}=\phi$,
$\mathcal{J}=\\{1,\ldots,d\\}$

for $t=1$ to $\kappa$ do

Choose
$\lambda^{t}=arg\max\limits_{j\in\mathcal{J}\setminus\Gamma}\left|<\varphi_{j},\bm{r}^{t-1}>\right|$

Update index set $\Gamma=\Gamma\bigcup\\{\lambda^{t}\\}$

Estimate $\hat{\bm{s}}^{t}=arg\min_{\bm{s}}\|\bm{y}-\Phi\bm{s}\|_{2}^{2}$ by
setting the $l$-th entry $s_{l}$=0 for $l\notin\Gamma$

Update residual $\bm{r}^{t}=\bm{y}-\Phi\hat{\bm{s}}^{t}$

end for

Return $\left\\{\hat{\bm{s}}^{t}\right\\}_{t=1}^{\kappa}$