1

Fast Greedy Approaches for Compressive Sensing
of Large-Scale Signals
Sung-Hsien Hsieh∗,∗∗, Chun-Shien Lu∗∗ , and Soo-Chang Pei∗
Inst. Comm. Eng., National Taiwan University, Taipei, Taiwan
∗∗ Institute of Information Science, Academia Sinica, Taipei, Taiwan

arXiv:1509.03979v2 [] 17 Mar 2016

∗ Graduate

Abstract—Cost-efficient compressive sensing is challenging
when facing large-scale data, i.e., data with large sizes. Conventional compressive sensing methods for large-scale data will
suffer from low computational efficiency and massive memory
storage. In this paper, we revisit well-known solvers called greedy
algorithms, including Orthogonal Matching Pursuit (OMP), Subspace Pursuit (SP), Orthogonal Matching Pursuit with Replacement (OMPR). Generally, these approaches are conducted by
iteratively executing two main steps: 1) support detection and 2)
solving least square problem.
To reduce the cost of Step 1, it is not hard to employ the sensing
matrix that can be implemented by operator-based strategy
instead of matrix-based one and can be speeded by fast Fourier
Transform (FFT). Step 2, however, requires maintaining and
calculating a pseudo-inverse of a sub-matrix, which is random
and not structural, and, thus, operator-based matrix does not
work. To overcome this difficulty, instead of solving Step 2
by a closed-form solution, we propose a fast and cost-effective
least square solver, which combines a Conjugate Gradient (CG)
method with our proposed weighted least square problem to
iteratively approximate the ground truth yielded by a greedy
algorithm. Extensive simulations and theoretical analysis validate that the proposed method is cost-efficient and is readily
incorporated with the existing greedy algorithms to remarkably
improve the performance for large-scale problems.
Index Terms—Compressed/Compressive sensing, Greedy algorithm, Large-scale Data, Least square, Sparsity

s

In this section, we first briefly introduce the background
of compressive sensing (CS) in Sec. I-A. Then, the existing
CS recovery algorithms (for large-scale signals) are discussed
in Sec. I-B. Finally, the overview and contributions of our
proposed method are described in Sec. I-C, followed by the
organization of the remainder of this paper.
A. Background
Compressive sensing (CS) [1][2][3] for sparse signals in
achieving simultaneous data acquisition and compression has
been extensively studied in the literature. CS is recognized to
be composed of fast encoder and slow decoder.
Let x ∈ RN denote a K-sparse 1-D signal to be sensed, let
Φ ∈ RM×N (M < N ) represent a sampling matrix, and let
y ∈ RM be the measurement vector. At the encoder, a signal
x is simultaneously sensed and compressed via Φ to obtain a
so-called measurement vector y as:

Corresponding author: Chun-Shien Lu; e-mail: lcs@iis.sinica.edu.tw

where A = ΦΨ ∈ RM×N . We say that x is K-sparse if s
contains only K non-zero entries (exactly K-sparse) or K
significant components (approximately K-sparse).
At the decoder, the original signal x can be perfectly
recovered by an intuitive solution to CS recovery, called ℓ0 minimization, which is defined as:
min ksk0

I. I NTRODUCTION

y = Φx,

which is usually called a procedure of ranfom projection. The
measurement rate, defined as 0 < M
N < 1, indicates the compression ratio (without quantization) and is a major concern
in many applications. [1][3] show that it is a good choice to
design Φ as a Gaussian random matrix to satisfy either mutual
incoherence property (MIP) or restricted isometry property
(RIP). Moreover, sparsity is an inherent assumption made in
compressed sensing to solve the underdetermined system in
Eq. (1) due to M < N . Nevertheless, for real applications,
natural signals are often not sparse in either the time or space
domain but can be sparsely represented in a transform (e.g.,
discrete cosine transform (DCT) or wavelet) domain. Namely,
x = Ψs, where Ψ is a transform basis (or dictionary) and s is
a sparse representation with respect to Ψ. So, Eq. (1) is also
rewritten as:
y = ΦΨs = As,
(2)

(1)

s.t. ky − Ask2 ≤ ǫ,

(3)

where ǫ is a tolerable error term. Due to M < N , this
system is underdetermined and there exists infinite solutions.
Thus, solving ℓ0 -minimization problem requires combinatorial
search and is NP-hard.
Alternative solutions to Eq. (3) usually are based on two
strategies: convex programming and greedy algorithms. For
convex programming, researchers [1][3] have shown that when
N
) holds, solving ℓ0 -minimization is equivaM ≥ O(K log K
lent to solving ℓ1 -minimization, defined as:
min ksk1
s

s.t. ky − Ask2 ≤ ǫ.

(4)

Typical ℓ1 -minimization models include Basis Pursuit (BP) [4]
and Basis Pursuit De-Noising (BPDN) [5] with the computational complexity of recovery being polynomial. For greedy
approaches, including Orthogonal Matching Pursuit (OMP)
[6], CoSaMP [7], and Subspace pursuit (SP) [8], utilize a
greedy strategy for support detection first and then solve a
least square problem to recover the original signal. The main
difference among these greedy algorithms is how the support
detection step is condcuted.

2

Nevertheless, for large-scale signals (e.g., N > 220 ), both
ℓ1 -minimization and greedy algorithms suffer from high computational complexity and massive memory usage. Ideally,
the memory costs at the encoder and decoder are expected
to approximate O(M ) and O(N ), respectively, which are
the minimum costs to store the original measurement vector
y and signal x. If A is required to be stored completely,
however, it will cost O(M N ) bytes (e.g., when N = 220
and M = 218 , the sensing matrix needs several terabytes). It
often overwhelms the capability of existing hardware devices.
In view of the incoming big data era, such a troublesome
problem needs to pay immediate attention. In this paper,
we say that X ∈ RN1 ×N2 ×···×ND is a kind of big data
if N1 × N2 × · · · × ND is lare enough or more specifically
its size approaches the storage limit of hardwares like PC,
notebook, and so on.
B. Related Work
The existing methods that can deal with compressive sensing of large-scale signals are discussed in this section. As
mentioned above, we focus on computational efficiency and
memory usage. Basically, our survey is conducted from the
aspects of encoder and decoder in compressive sensing. We
mainly discuss block-based, tensor-based, and operator-based
compressive sensing algorithms here.
1) Strategies at Encoder: From Eq. (1), we can see that
both the storage (for Ψ) and computation (for Ψx) costs
require O(M N ) bytes and O(M N ) operations, respectively.
When the signal length becomes large enough, storing Φ and
computing Φx become an obstacle.
In the literature, Gan [11] and Mun and Fowler [12] propose
block-based compressive sensing techniques, wherein a largescale signal is separated into several small block signals,
which are individually sensed via the same but smaller sensing
matrix. The structure of block sensing reduces both storage
and computation costs to O( MN
B ), where B is the number
of blocks. Although block-based compressive sensing can
deal with small blocks quickly and easily, it actually cannot
work for the scenario of medical imaging in that an image
generated from the fast Fourier Transform (FFT) coefficients
of an entire sectional view [18] violates the structure of blockbased sensing.
Shi et al. [19] and Caiafa and Cichocki [9] consider the
problem of large-scale compressive sensing based on tensors.
In other words, the signal is directly sensed and reconstructed
in the original (high) dimensional space instead
√of reshaping
√
N× N
is sensed
to 1-D. For example, a 2-D image X ∈ R
via
Y = Φ1 XΦT2 ,
(5)
√

√

√

√

where Φ1 and Φ2 ∈ R M × N , and Y ∈ R M × M . This
strategy is often called separable sensing [20], [21]. In this
case,
√ both the storage and computation costs are reduced to
O( M N ). [10] further presents a close-from solution for reconstruction from compressive sensing based on assuming the
low-rank structure. It should be noted that since tensor-based
approaches, in fact, change the classical sensing structure (i.e.,
y = Φx) of CS, the decoder no longer follows the conventional

solvers like Eq. (4). Specifically, the measurements in tensorbased approaches form a tensor but conventional solvers only
accept one-dimensional measurement vector.
In addition to block-based and tensor-based approaches,
operator-based approaches are to design Φ as a deterministic
matrix or structurally random matrix, implemented by certain
fast operators. For example, Candes et al. [22] propose the
use of a randomly-partial Fourier matrix as Φ. In this case,
we can implement Φx by D (F F T (x)), where F F T (·) is the
function of fast Fourier transform (FFT) and D (·) denotes
a downsampling operator that outputs an M × 1 vector.
Thus, Φ is not necessarily stored in advance. In addition, the
computation cost also becomes O(N log N ), which especially
outperforms O(M N ) for large-scale signals because M is
positively proportional to N . Do et al. [23] further propose
a kind of random Gaussian-like matrices, called Structurally
Random Matrix (SRM), which benefits from operator-based
strategy and achieves reconstruction performance as good
as random Gaussian matrix. In sum, since operator-based
approaches follow the original CS structure, the decoder is
not necessary to be modified.
2) Strategies at Decoder: For block-based approaches
[11][12], each block can be individually recovered with low
computation cost and memory usage but incurs blocky effects
between boundaries of blocks. In [12], Mun and Fowler propose a method, called BCS-SPL, which further removes blocky
effects by Wiener filtering. In addition to the incapabliity
of sensing medical images like MRI, BCS-SPL is also not
adaptive in that the measurement rates are fixed for different
blocks by ignoring the potential differences in smoothing
blocks that need less measurement rates and complex blocks
that require more measurement rates.
For tensor-based compressive sensing, [9] develops a new
solver called N-way block OMP (N-BOMP). Though NBOMP is indeed faster than conventional CS solvers, its
performance closely depends on the unique sparsity pattern,
i.e., block sparsity, of an image. Specifically, block sparsity
states that the importnat components of an images are clustered
together in blocks. This characteristic seems to only naturally
appear in hyperspectral imaging. In [24], a multiway compressive sensing (MWCS) method for sparse and low-rank tensors
is proposed. MWCS achieves more efficient reconstruction,
but its performance relies heavily on tensor rank estimation,
which is NP-hard. A generalized tensor compressive sensing
(GTCS) method [25], which combines ℓ1 -minimization with
high-order tensors, is beneficial for parallel computation.
For operator-based compressive sensing algorithms, since
the conventional solvers, mentioned in the previous subsection,
still can be used, here we mainly review state-of-the-art convex
optimization algorithms focusing on the large-scale problem,
where only simple operations such as A and AT conducted
by operator are required.
Cevher et al. [26] point out that an optimization algorithm
based on the first-order method such as gradient descent
features nearly dimension-independent convergence rate and
is theoretically robust to the approximations of their oracles.
Moreover, the first-order method such as NESTA [27] often
involves the transpose of sensing matrix, which is easily

3

TABLE I
C OMPARISON BETWEEN EACH ALGORITHM .
algorithms
N-BOMP [9]
[10]
BCS [11]
BCS-SPL [12]
[13]
GPSR [14], SpaRSA [15]
[16][17]

sensing strategy
tensor-based (2D)
tensor-based (2D)
block-based (1D)
block-based (1D)
conventional (1D)
conventional (1D)
conventional (1D)

assumptions
block sparsity
low multilinear-rank
-

implemented by operator. Both GPSR [14] and SpaRSA
[15] are closely related to iterative shrinkage/threshold (IST)
methods and support the operator-based strategy. In addition,
[16][17][28] have shown that algorithms based on solving
fixed-point equation have fast convergence rate, which can
be combined with operator-based strategy too. For example,
Milzareket al. [13] further propose a globalized semismooth
Newton method, where partial DCT matrix is adopted as
the sensing matrix for fast sensing. But, it requires that
signals are sparse in the time/spatial domain leading to limited
applications.
3) Brief Summary of Related Works: Table I depicts the
comparisons among the aforementioned algorithms, where
storage is estimated based on non-operator version. If operator
can be used, the storage of storing a sensing matrix is not
required and, thus, is bounded by O(N ) for each row vector,
which is only related to the minimum requirement for storing
the reconstructed signal. Since the characteristic of compressive sensing states that CS encoder spends lower memory
and computation cost than CS decoder, when taking hardware
implementation in real world into consideration, tensor-based
methods are more complicated than others. For example, the
single-pixel camera designed in [29] uses a DMD array as
a row of A to sense x. By changing the pattern of DMD
array M times, the measurements are collected. This structure,
however, cannot support separable sensing that is commonly
used in tensor-based methods. Block-based CS methods do
not intrinsically overcome large-scale problems and lack convincing theoretical proof about complexity, performance, and
convergence analysis. Operator-based CS methods maintain
the original structures of CS encoder and decoder. Thus, most
of the existing fast algorithms for ℓ1 -minimization can be used
only if all of matrix operations can be executed in an operator
manner. Furthermore, they have strong theoretical validation
since ℓ1 -minimization is a well-known model and has been
developed for years. In fact, the operator-based strategy can
also be employed in tensor-based and block-based compressive
sensing methods to partially reduce their computation cost and
storage usage.
C. Contributions and Overview of Our Method
Up to now, it is still unclear how greedy algorithms can deal
with large-scale problems by utilizing operator-based strategy.
Although SparseLab releases the OMP code combined with
operator, the program still cannot deal with large-scale signals.
This challenge is the objective of this paper and, to our

algorithm type
greedy
closed-form
Landweber-based
Landweber-based
Armijo-based
IST-based
FPC AS-based

storage
√
O(√M N )
O( M N )
O( MBN )
O( MBN )
O(M N )
O(M N )
O(M N )

knowledge, we are the first to explore this issue. In fact, our
idea can help all greedy algorithms to deal with large-scale
signals. We will discuss the problem in detail in Sec. II.
Generally, greedy algorithms are conducted by iteratively
executing two main stages: (a) support detection and (b)
solving least square problem with the known support. To
reduce the cost of support detection, we follow the common
strategy of adopting operator-based, instead of matrix-based,
design of a sensing matrix (e.g., [23]). Therefore, we no longer
discuss this step in this paper as it is not the focus of our
method.
For solving the least square problem, we propose a fast
and cost-effective solver by combining a Conjugate Gradient
(CG) method with a weighted least square model to iteratively
approximate the ground truth. In our method, the memory cost
of solving least square problem is reduced to O(N ), and the
computation cost of CG method is approximately O(N log N )
for finite floating point precision and O(KN log N ) for exact
precision.
In should be noted that although using CG to solve the least
square problem is not new, our extended use of CG brings
additional advantages. For example, Blumensath et al. [30]
proposed “Gradient Pursuits (GP)” in which the memory cost
is dominated by O(M N ) to save Φ (see Table 1 in [30]),
which cannot be stored explicitly for large-scale problem. Our
method extends GP to reduce the memory cost by using SRM
to avoid saving Φ. In addition, we reformulate a least square
problem used in GP into a weighted one and show both models
are equivalent. More specifically, solving weighted least square
problem only requires Φ and can benefit from fast computation
of operator-based approaches (as in SRM). Traditional least
square problem, however, involves sub-matrices of Φ and
cannot directly be conducted by fast operator.
On the other hand, “iteratively reweighted least-square
(IRLS)” was proposed in [31]. Though both IRLS and our
method involve weighting, they are totally different. First,
IRLS uses weighting to approximate ℓ1 -norm solution instead
of ℓ2 -norm solution in original least square problem while our
greedy method still solves ℓ2 -norm solution in the weighted
least square problem. Second, IRLS is not a greedy method.
Moreover, we conduct extensive simulations to demonstrate
that our method can greatly improve OMP [6], Subspace
pursuit (SP) [8], and OMPR [32] in terms of memory usage
and computation cost for large-scale problems.

4

D. Outline of This Paper
The rest of this paper is organized as follows. In Sec. II,
we describe the bottleneck of current greedy algorithms that
is to solve the least square problem. The proposed idea of
fast and cost-effective least square solver for speeding greedy
approaches along with theoretical analysis is discussed in Sec.
III. In Sec. IV, extensive simulations are conducted to show
that our method indeed can be readily incorporated with stateof-the-art greedy algorithms, including OMP, SP, and OMPR,
to improve their performance in terms of the memory and
computation costs. Finally, conclusions are drawn in Sec. V.
II. P ROBLEM S TATEMENT
In this paper, without loss of generality, we focus on a
signal ∈ RN1 ×N2 ×···×ND , where N = N1 and N = N1 × N2
are large enough with respect to D = 1 and D = 2. When
D = 2, the signal usually is reshaped to 1-D form in the
context of compressive sensing. Ideally, the memory cost of a
compressive sensing algorithm should be O(N ), which is the
minimum requirement for saving the original signal, x. The
computational cost, however, depends on an algorithm itself.
Since greedy algorithms share the same framework composed
of support detection and solving least square problem, we shall
focus on reducing the costs of these two procedures.
In this section, we discuss the core of proposed fast and
cost-effective greedy approach and take OMP as an example
for subsequent explanations. We shall point out the dilemma
in terms of memory cost and computation cost when handling
large-scale signals. In fact, both costs suffer from solving the
least square problem, which cannot be conducted by operator
directly.
First, we follow the notations mentioned in the previous
section and briefly introduce OMP [6] in a step-by-step manner
as follows.
1) Initialize the residual measurement r0 = y and initialize
the set of selected supports S0 = {}. Let the initial
iteration counter be i = 1. Let AS be the sub-matrix of
A, where AS consists of the column of A with indices
belonging to the support set S. AT is the transpose of A.
2) Detect supports (or positions of significant components)
by seeking maximum correlation from

(6)
t = argmax | AT ri−1 t |,

calculated, where ŷ = [y T , 0, ..., 0]T and IF F T (·) denotes
| {z }
N −M

inverse FFT function. In this paper, we call it operator-based
OMP (O-OMP).
Unfortunately, the key is that (ATSi ASi )−1 still cannot be quickly computed in terms of operator. Hence,
it requires O(KM ) to store ASi . Instead of calculating
(ATSi ASi )−1 directly, by preserving the Cholesky factorization
of (ATSi−1 ASi−1 )−1 at the (i − 1)th iteration for subsequent
use, Step 3 is accelerated and the memory cost is reduced
to O(K 2 ). Moreover, it is worth mentioning that the sparsity
K of natural signals is often linear to signal length N . For
example, the number of significant DCT coefficients for an
image usually ranges from 0.01N to 0.1N . Also CS has shown
that M must be linear to K log N for successful recovery with
high probability. Under the circumstance, O(N ) is equivalent
to O(K) and O(M ) in the sense of big-O notation. We can
see that when N is large, O(K 2 ) dominates the memory cost
because K 2 ≫ N . Thus, Step 3 makes OMP infeasible for
recovering large-scale signals.
In fact, greedy algorithms share the same operations, i.e.,
AT ri−1 in Step 2 and (ATSi ASi )−1 ATSi y in Step 3, where
the main difference is that the support set Si is found
by different ways, and face the same dilemma. A simple
experiment is conducted and results are shown in Fig. 1
to illustrate the comparison of memory usage among MOMP, O-OMP, and ideal cost (which is defined as N × 8
bytes required for Double data type in Matlab). The OMP
code running in Matlab was downloaded from SparseLab
(http://sparselab.stanford.edu). Obviously, though the memory
cost of O-OMP is reduced without storing A, it still far higher
than that of ideal cost. It is also observed that both M-OMP
and O-OMP exhibit the same slope. Specifically, M-OMP and
O-OMP cost O(M N ) and O(K 2 ), respectively. As mentioned
before, since M, K are linear to N , it means O(K 2 ) = O(N 2 )
and O(M N ) = O(N 2 ) such that both orders of memory cost
of M-OMP and O-OMP are the same and are larger than that
of ideal case. Consequently, solving the least square problem
becomes a bottleneck in greedy approaches. This challenging
issue will be solved in this paper.

3

10

and update the support set Si = Si−1 ∪ {t}.
3) Solve a least square problem si = (ATSi ASi )−1 ATSi y, and
update residual measurement ri = y − ASi si .
4) If i = K, stop; otherwise, i = i + 1 and return to Step 2.
Tropp and Gilbert [6] derive that the computational complexity of OMP is bounded by Step 2 (support detection)
with O(M N ) and Step 3 (solving least square problem)
with O(M K), and the memory cost is O(M N ) when A is
executed in a matrix form. In this paper, we call it matrixbased OMP (M-OMP). As mentioned in Sec. I-B, A can be
designed to be an SRM conducted by operator. Nevertheless,
operator is only helpful for certain operations such as A and
AT . For example, if A is a partial random Fourier matrix,
Ax = D (F F T (x)) and AT y = IF F T (ŷ) can be quickly

Memory Cost (MegaBytes)

t

2

10

1

10

0

10

−1

10

Ideal Cost
O−OMP
M−OMP

−2

10

−3

10

10

11

12

13

14

15

16

17

log2 N
Fig. 1. Memory cost comparison among the matrix-based OMP, operatorbased OMP, and ideal cost under M = N
and K = M
.
4
4

5

III. P ROPOSED M ETHOD FOR S PEEDING R ECOVERY
G REEDY A LGORITHMS

OF

In this section, we first introduce how to determine a
sensing matrix, which can be easily implemented by operator.
Then, we reformulate the least square problem as a weighted
one, which is solved by conjugate gradient (CG) method
to avoid involving the sub-matrices, ASi or A†Si . We first
prove that the solutions to the least square problem and its
weighted counterpart are the same, and then prove that the
solutions to the weighted least square problem and its CGbased counterpart are the same.
A. Sensing Matrix
A random Gaussian matrix is commonly used as the sensing
matrix as it and any orthonomal basis can pair together to
satisfy RIP and MIP in the context of compressive sensing.
The use of random Gaussian matrix as the sensing matrix,
however, leads to the overhead of storage and computation
costs. Although storage consumption can be overcome by
using a seed to generate a random Gaussian matrix, it still
encounters high computational cost.
In [23], Do et al. propose a framework, called Structurally
Random Matrix (SRM), defined as:
Φ = DF R,

(7)

where D ∈ RM×N is a sampling matrix, F ∈ RN ×N is an
orthonormal matrix, and R ∈ RN ×N is a uniform random permutation matrix (randomizer). Since the distributions between
a random Gaussian matrix and SRM’s Φ are verified to be
similar, we choose Eq. (7) as the sensing matrix for our use.
It should be noted that D and R can inherently be replaced by
operators but F depends on what kind of orthonormal basis is
used. It is obvious that any fast transform can be adopted as F .
In our paper, we set F to the Discrete Cosine Transform (DCT)
due to its fast computation and cost-effectiveness. There are
literatures discussing the design of sensing matrix but it is not
the focus of our study here.
B. Reformulating Least Square Problem: Weighted Least
Square Problem
In Sec. II, we describe that the bottleneck of greedy algorithms is the least square problem. To solve the problem,
it is reformulated as a weighted least square problem in our
method. To do that, we first introduce a weighted matrix
W ∈ RN ×N defined as:
(
1, j ∈ Si
,
(8)
WSi [j, j] =
0, j 6∈ Si
where, without loss of generality, Si = {1, 2, ..., i} denotes a
support set at the i-th iteration and WSi [j, j] is the (j, j)th
entry of WSi . As can be seen in Eq. (8), the weighte matrix
W is designed to make the supports unchanged.
Then, we prove that both solutions of the least square
problem and its weighted counterpart are the same.

Theorem 1. Suppose the sub-matrix ASi ∈ RM×K of A has
full column rank with support set Si . Let si ∈ RK be the
solution to the least square problem:
si = argmin ky − ASi ŝk2 ,
ŝ

(9)

and let θi ∈ RN be the solution to the weighted least square
problem:
(10)
θi = argmin ky − AWSi θ̂k2 .
θ̂

We have
si [j] = θi [j],

for 1 ≤ j ≤ i.

−1 T 
ATSi ASi
ASi
T
T
Proof. Let θ∗ = [si 0] =
y be a
0
solution minimizing ky − AWSi θ∗ k2 . Then {θ∗ + v|v ∈
N ull(AWSi )} is the solution set of Eq. (10). Since ASi has
full column rank with rank(ASi ) = K and N ull(AWSi ) =
span (ei+1 , ei+2 , ..., eN ), where ei is a standard basis. Thus,
no matter what v is, the first i entries of θ∗ + v are invariant.
We complete the proof.
C. Reformulating Weighted Least Square as CG-based
Weighted Least Square Problem
We can see from Eq. (10) that the introduce of weighted matrix involves AWSi instead of submatrix ASi so that A can be
calculated by fast operator.
Nonetheless,

−1 T  the closed-form soluT
ASi ASi
ASi
tion of Eq. (10) is
y and still faces the dif0

−1 T
ASi cannot be easily implemented
ficulty in that ATSi ASi
by operator. Instead of seeking closed-form solutions, we aim
to explore the first-order methods (e.g., gradient descent),
which have the following advantages: (1) the operations only
involve A and AT instead of the pseudo-inverse of A, (2) the
convergence rate is nearly dimension-independent [26], and (3)
if A is a sparse matrix, the computation cost can be further
reduced. Conjugate gradient (CG) [33][34] is a well-known
first order method to numerically approximate the solution of
symmetric, positive-definite or positive-semidefinite system of
linear equations. Thus, CG benefits from the advantages of
the first-order method. However, the matrix AWSi in Eq. (10)
is not symmetric. Thus, we reformulate Eq. (10) in terms of
CG as follows. We will prove that both the solutions to the
weighted least square problem and CG-based weighted least
square problem are the same.
Theorem 2. Suppose the sub-matrix ASi ∈ RM×K of A has
full column rank with support set Si . Let θi ∈ RN be the
solution to weighted least square problem defined in Eq. (10).
Let θe ∈ RN be the solution to the CG-based weighted least
square problem reformulated from Eq. (10) as:
θei = argmin kWSTi AT (y − AWSi θ̂)k2 .

(11)

θ̂

Then

θi [j] = θei [j],

for 1 ≤ j ≤ i.

Proof. We first note that now WSTi AT AWSi is symmetirc to

−1 T 
ATSi ASi
ASi
meet the requirement of CG. Since
y is an
0

6

optimal solution for both Eq. (10) and Eq. (11), the solution
set of Eq. (11) can be expressed as:

−1 T 
ATSi ASi
ASi
{
y + v|AWSi v ∈ N ull(WSTi AT )}.
0
It should be noted that N ull(WSTi AT ) = N ull(ATSi ).
In addition, AWSi v ∈ C(AWSi ), where C(AWSi ) denotes the column space
T of AWSi and C(AWSi ) =
C(ASi ). Since C(ASi ) N ull(ATSi ) = {0}, it implies
AWSi v = 0. In other words, v ∈ N ull(AWSi ). Due to
N
i ) = span
 ull(AWS−1
 (ei+1 , ei+2 , ..., eN ), the first i entries of
ATSi ASi
ATSi
y + v are invariant. Similarly, the solution
0
set of Eq. (10) is

−1 T 
ATSi ASi
ASi
{
y + v|v ∈ N ull(AWSi )}.
0

−1 T 
ATSi ASi
ASi
The first i entries of
y + v are also invari0
ant. We complete the proof.
So far, we prove that, with correct support detection, the optimal solution to the CG-based weighted least square problem
in Eq. (11) is equivalent to that to the original least square
problem in Eq. (9). In addition, the matrix WSTi AT AWSi
in Eq. (11) is symmetric and can quickly be solved by CG
method. Nevertheless, [34] points out that the system of linear
equations with a positive-semidefinite matrix diverges unless
some conditions are satisfied. Thus, Theorem 3 further shows
the condition of convergence.
Theorem 3. [34] If WSTi AT y ∈ C(WSTi AT AWSi ), CG
method converges but the solution is not unique.
Now, we check whether the CG-based weighted least square
problem converges. Again, let Si = {1, 2, ..., i} denote a
support set. Then, we have WSTi AT y = [y T ASi 0]T and
 T

A A
C(WSTi AT AWSi ) = C( Si Si ).
0

Because ATSi ASi ∈ Ri×i is full rank, the first i entries of
WSTi AT y must be spanned by the basis of ATSi ASi . The
remaining N − i entries are 0 and trivial. Thus, it implies
the CG part of our CG-based weighted least square solver
satisfies Theorem 3 and converges.
D. Speeding Orthogonal Matching Pursuit and Complexity
Analysis
In the previous section, we descirbe the proposed CGbased weighted least square solver. In this section, we first
show that the CG-based solver can be implemented easily by
operator. Then, we combine it with OMP as a new paradigm
to achieve fast OMP. Moreover, we discuss the convergence
rate of CG and derive the computation complexity of proposed
operator-based OMP via CG (dubbed as CG-OMP). Finally,
we conclude that the memory cost of CG-OMP achieves ideal
O(N ) if Ψ is also conducted by operator.
Algorithm 1 describes the proposed CG-OMP method
(Lines 01 - 10), which employs a CG technique (Lines 11

Algorithm 1 Proposed Orthogonal Matching Pursuit
Input: y, A, K; Output: sK ;
Initialization: i = 1, r0 = y,s0 = AT y, S0 = {};
01. function Proposed Operator-based OMP()
02. for i = 1 to K
03.
t = argmaxt̂ (AT ri−1 )t̂ ;
04.
Si = Si−1 + t;
05.
Assign WSi according to Eq. (8);
06.
si =CG(A, WSi , y, i);
07.
ri = y − Asi ;
07.
i = i + 1;
08. end for
09. Return:sK ;
10. end function
11. function [ ŝ ]=CG(A, WSi , y, i)
12. b = WSTi AT y, H = WSTi AT AWSi ;
13. d0 = r0 = b, ŝ = 0;
14. j = 0;
15. while(krj k2 ≤ ξ)
rjT rj

16.

αj =

17.
18.

ŝ = ŝ + αj dj ;
rj+1 = rj − αj Hdj ;

19.

βj+1 =

dT
j Hdj

;

T
rj+1
rj+1

rjT rj

;

20.
dj+1 = rj+1 + βj+1 dj ;
21.
j = j + 1;
22. end while
23. Return:ŝ;
24. end function

- 23). It is worth mentioning that ξ in Line 15 controls the
precision of CG method. If ξ = 0, it means exact precision
such that the output of CG method is equal to least square
solution. If ξ > 0, CG method attains finite precision. However, the result with finite precision is not certainly worse than
that with exact precision especially under noisy interference,
as later discussed in the 4-th paragraph of Sec. IV-C.
Now we check whether all matrix operations can be implemented by operators in Algorithm 1 in the following.
•

•

•

WSi ∈ RN ×N :
WSi is a diagonal matrix. Thus, WSi x is equal to assign
x[j] = 0 for j 6∈ Si . The memory cost is O(K) to store
the indices of support set and the computation cost is
O(N ).
Φ = DF R ∈ RM×N :
Dx is equal to randomly choose M indices from N
entries in x. F x = DCT (x), where DCT can be speeded
by FFT. Rx is equal to randomly permute the indices
of vector x. The memory cost is bounded by O(N ) in
order to store the sequence of random permutation and
the computation cost is O(N log N ).
Ψ ∈ RN ×N :
If Ψ is a deterministic matrix, it is not necessarily to
be stored. Thus, the memory cost of Ψx is O(N ).
The computational cost depends on whether Ψx can be
speeded up. In the worst case, it costs O(N 2 ).

Other operations such as rjT rj (Line 16) only involve multiplications between vectors. Thus, both memory cost and

7

computation cost are bounded by O(N ).
From the above analysis, one can see that the computation
cost of Algorithm 1 is mainly bounded by Ψ. We discuss some
applications below, where the computation cost involving Ψ
is low. For compressive sensing of images, wavelet transform
is often chosen as Ψ such that Ψx costs O(N ). For magnetic
resonance imaging (MRI), partial Fourier transform is selected
as the component F of the sensing matrix expressed as DF R.
In this case, Ψ is I in order to satisfy MIP or RIP, and costs
nothing. Spectrum sensing is another application, where Ψ is
a discrete Fourier transform matrix done with O(N log N ).
Now, the total cost of Algorithm 1 is discussed. Both the
computation and memory costs of Lines 1-10 except Line 6
(CG method) will be O(N log N ) and O(N ), respectively. As
for the memory cost of CG method, it needs O(N ). Therefore,
the total memory cost of Algorithm 1 is bounded by O(N ).
In addition, the computation cost is related to two factors,
i.e., the number of iterations to converge in CG and Ψ. They
are further discussed as follows.
Theorem 4. (Theorem 2.2.3 in [35]) Let H be symmetric and
positive-definite. Assume that there are exactly k < N distinct
eigenvalues in H. Then, CG terminates in at most k iterations.
In our case, H = WSTi AT AWSi is positive-semidefinite
instead of positive-definite. Thus, we derive the following
theorem.
Theorem 5. Given b = WSTi AT y and H = WSTi AT AWSi ,
solving Eq. (11) requires the number of iterations at most K
in CG, where K is the sparsity of an original signal.
Proof. Without loss of generality, let support set Si =
{1, 2, ..., i}. We start from another optimization problem:
s̄i = argmin kATSi (y − ASi ŝ)k2 .
ŝ

(12)

Following the same skill in Theorem 2, s̄i is a unique and
optimal solution to both Eq. (9) and Eq. (12). Thus, s̄i is also
the solution of Eq. (11) for the first i entries. Furthermore,
ATSi ASi is non-singular such that ATSi ASi is a positive-definite
matrix and has at most i distinct eigenvalues. From Theorem
4, solving Eq. (12) requires at most i iterations. Then, we
have b = WSTi#AT y = [y T ASi 0]T and H = WSTi AT AWSi =
"
ATSi ASi 0
in Eq. (11). When we only take the first
0
0
i entries of b and the left-top i × i submatrix of H into
consideration, it is equivalent to solving Eq. (12). This fact can
be checked trivially by comparing each step of CG for both
optimization problems in Eq. (11) and Eq. (12). Thus, the first
i entries of ŝ in Eq. (11) is updated in the same manner with
that of s̄ in Eq. (12). The remaining N − i entries of ŝ are
unrelated to convergence because the N − i entries of H ŝ are
zero. In sum, the required number of iterations to converge
in Eq. (11) is identical to that in Eq. (12). Hence, solving
Eq. (11) also requires at most i iterations. Since i ≤ K, the
number of iterations is at most K. We complete the proof.
Moreover, for each iteration in CG, the operation, Hdj ,
dominates the whole computation cost. It is obvious that if Ψ
can be executed with O(N log N ) or even lower computation

complexity, Hdj costs O(N log N ) since H involves Ψ,
which spends O(N log N ). Note that the total computation
complexity of CG-OMP will be O(N K 2 log N ), where K 2
comes from the outer loop in OMP, which needs O(K), and
solving Eq. (11) that requires the number of iterations at most
K in CG, as proved in Theorem 5. On the other hand, if, under
the worst case, Ψ costs O(N 2 ) operations, the computation
complexity of CG-OMP will be O(N 2 K 2 ). For applications
that accept finite-precision accuracy instead of exact precision,
CG [33] requires fewer steps (≤ K) to achieve approximation.
Under the circumstance, the complexity of CG-OMP nearly
approximates O(N K log N ) and O(N 2 K) operations for Ψ
with complexity O(N log N ) and O(N 2 ), respectively.
Consequently, a reformulation for solving a least square
problem is proposed based on CG such that the new matching
pursuit methodology can deal wtih large-scale signals quickly.
It should be noted that, in the future, CG may be substituted
with other first order methods that outperform CG. The
proposed idea can also be readily applied to other greedy
algorithms to enhance their performance.
E. Strategies for Reducing the Cost of Ψ
We further consider that if Ψ is a learned dictionary, it
will become a bottleneck for operator-based algorithms since
it requires O(N 2 ) for storage. To overcome this difficulty, Ψ
should √
be learned
in a tensor structure. Let x = vec(X), where
√
X ∈ R N × N and vec(·) is a vectorization operator. That is,
a two-dimensional vector is reshaped to a one-dimensional
vector. Then, we can learn a 2D dictionary such that x =
Ψs = vec(Ψ1 SΨT2 ) with s = vec(S) and Ψ = Ψ2 ⊗ Ψ1 (⊗ is
a Kronecker product). Under the circumstance, all operations
in CG involving ΦΨs can be replaced by Φvec(Ψ1 SΨT2 ).
Moreover, both Ψ1 and Ψ2 only require O(N ) in terms of
memory cost. In the literature, the existing algorithms for 2D
separable dictionary learning include [19][36][37].
IV. E XPERIMENTAL R ESULTS
In this section, we conduct comparisons among O-OMP,
M-OMP, and CG-OMP in terms of the memory cost and
computation cost. The code of OMP was downloaded from
SparseLab (http://sparselab.stanford.edu). We have also applied the proposed fast and cost-effective least square solver
to SP and OMPR in order to verify if our idea can speed the
family of matching pursuit algorithms. For SP and OMPR, we
implemented the corresponding original matrix-based versions
(M-OMPR and M-SP), original operator-based versions (OOMPR and O-SP), and proposed CG-based (CG-OMPR and
CG-SP). It should be noted that although both SP and OMPR
work well for increasingly adding a index to the support set
like OMP, they are not accelerated by Cholesky factorization.
A. Simulation Setting
The simulations were conducted in an Matlab R2012b
environment with an Intel CPU Q6600 and 4 GB RAM under
Microsoft Win7 (64 bits).
The model for the measurement vector in CS is y = As+η,
where η is an addictive Gaussian noise with standard deviation

8

ση . The input signal s was produced via a Gaussian models
as:

2
s ∼ pN 0, σon
,
(13)

which was also adopted in [38]. In Eq. (13), p is the probability
of the activity of a signal and controls the number of nonzero entries of x. Sparsity K is defined to be K = pN .
σon is standard deviation for input signal. Φ is designed from
SRM and Ψ is chosen to be a discrete cosine transform. In
the following experiments, M = N4 , K = M
4 , p = 0.0625,
σon = 1, and ση = 0.01.
B. Memory Cost Comparison
Fig. 2 shows the comparison in terms of memory cost vs.
signal length N . Since the matrix-based algorithms (M-OMP,
M-OMPR, and M-SP) run out of memory, their results are not
shown in Fig. 2 (note that the result regarding matrix-based
OMP can be found in Fig. 1). First, we can observe from
Fig. 2 that O-OMP, O-OMPR, and O-SP still require about
O(N 2 ) and fail to work when N > 218 . Second, in contrast
with O-OMP, O-OMPR, and O-SP although CG-OMP, CGOMPR, and CG-SP need more memory costs than the ideal
cost, which is O(N ), their slopes are nearly identical, which
seems to imply that the CG-based versions incur larger Big-O
constants. Thus, the proposed idea of fast and cost-effective
least square solver is readily incorporated with the existing
greedy algorithms to improve their capability of handling
large-scale signals.
4

Memory Cost (MegaBytes)

10

3

10

Figs. 3(a), (b), and (c) show the computation cost vs.
signal length for OMP, OMPR, and SP, respectively, under
the condition that the precision of CG was set to be exact.
In other words, we fix the same reconstruction quality for all
comparisons in Fig. 3 and discuss the computation costs for
these different versions of algorithms. It should be noted that
some curves are cut because of running out of memory.
From Fig. 3, it is observed that the operator-based strategy
(denoted with dash curves) or our proposed CG-based method
(denoted with solid curves) can effectively reduce the order of
computation cost in comparison with the matrix-based strategy
(denoted with solid-star curves). More specifically, in Fig.
3(a), it is noted that M-OMP is only fast than CG-OMP with
N ≤ 214 due to smaller Big-O constant. However, CG-OMP
outperforms M-OMP in the end since the order of computation
complexity of CG-OMP is lower than that of M-OMP. In
particular, such improvements are significant for large-scale
signals (with large N ). Moreover, in Figs. 3(b) and (c), OOMPR and O-SP have the same orders with M-OMPR and
M-SP because no Cholesky factorization is used.
On the other hand, when finite precision is considered, we
consider two cases of setting ξ = 10−5 and ξ = 10−10 to
verify that finite precision is adequate under the condition of
noisy interferences. Taking exact precision as the baseline, the
difference of SNR values for reconstruction between settings
for exact precision and ξ = 10−5 is about ±0.1 dB. Similarly,
the difference between exact precision and ξ = 10−10 is
about ±0.01 dB. Tables II, III, and IV further illustrate the
comparisons of computation costs under different precisions.
The precision setting ξ = 10−10 results in about four times
faster than the exact precision but only sacrifices ±0.1 dB
for the reconstruct quality, which is acceptable for many
applications. In fact, the performance occasionally is better
because exact precision may lead to over-fitting.

2

10

V. C ONCLUSIONS

Ideal Cost
CG−OMP
CG−OMPR
CG−SP
O−OMP
O−OMPR
O−SP

1

10

0

10

−1

10

16

17

18

19

20

21

22

23

24

log2 N

Fig. 2. Memory cost vs. Signal length, where M =
(p = 0.0625).

N
4

and P =

M
4

The bottleneck of greedy algorithms in the context of
compressive sensing is to solve the least square problem, in
particular, when facing large-scale data. In this paper, we
address this challenging issue and propose a fast but costeffective least square solver. Our solution has been theoretically proved and can be readily incorporated with the
existing greedy algorithms to improve their performance by
significantl reducing computation complexit and memory cost.
Case studies on combining our method and OMP, SP, and
OMPR have been conducted and shown promising results.

C. Computation Cost Comparison
Before illustrating the computation cost comparison, we
first discuss the convergence condition of CG in Algorithm
1 as follows: 1) For exact precision, krj k2 = 0 is set as the
stopping criterion. As described in Theorem 5, it costs at most
K iterations to converge. 2) With inexact or finite precision,
we set krj k2 ≤ ξ with ξ > 0. Though the precision is finite,
the signals, in practice, are interfered with noises and solutions
with finite precision are adequate. Under the circumstance, the
number of iterations required to converge will be significantly
decreased.

VI. ACKNOWLEDGMENT
This work was supported by Ministry of Science and
Technology, Taiwan, ROC, under grants MOST 104-2221-E001-019-MY3 and NSC 104-2221-E-001-030-MY3.
R EFERENCES
[1] D. L. Donoho, “Compressed sensing,” IEEE Transactions on Information Theory, vol. 52, no. 4, pp. 1289–1306, 2006.
[2] R. Baraniuk, “Compressive sensing,” IEEE Signal Processing Magazine,
vol. 24, no. 4, pp. 118–121, 2007.

9

TABLE II
T HE COMPUTATIONAL COST ( IN SECONDS ) OF CG-OMP UNDER DIFFERENT PRECISIONS .

o

N
Exact precision
krj k2 ≤ 10−10
krj k2 ≤ 10−5

211
1.2
0.6
0.4

212
4.7
2.3
1.8

213
16.2
8.7
7.0

214
59.8
35.5
26.9

215
241.0
131.3
103.6

216
980.3
564.7
432.1

217
3901.7
2143.9
1662.8

218
15301.4
8341.3
6479.2

TABLE III
T HE COMPUTATIONAL COST ( IN SECONDS ) OF CG-OMPR UNDER DIFFERENT PRECISIONS .

o

N
Exact precision
krj k2 ≤ 10−10
krj k2 ≤ 10−5

211
0.2
0.1
0.06

212
0.3
0.2
0.09

213
0.6
0.3
0.2

214
1.0
0.5
0.3

215
1.9
1.0
0.6

216
4.1
2.2
1.3

217
10.5
4.8
2.9

218
23.6
10.1
6.0

219
51.8
23.6
13.1

220
111.3
49.8
27.7

221
232.1
113.9
60.3

222
495.3
249.3
131.3

223
1109.1
532.3
293.4

224
2392.3
1130.4
643.6

223
1581.3
784.5
513.2

224
3195.4
1706.3
1096.8

TABLE IV
T HE COMPUTATIONAL COST ( IN SECONDS ) OF CG-SP UNDER DIFFERENT PRECISIONS .

o

N
Exact precision
krj k2 ≤ 10−10
krj k2 ≤ 10−5

211
0.3
0.2
0.09

212
0.5
0.3
0.2

213
0.9
0.5
0.4

214
1.5
0.8
0.6

215
3.0
1.6
0.9

216
6.2
3.6
2.1

[3] E. J. Candes and M. B. Wakin, “An introduction to compressive
sampling,” IEEE Signal Processing Magazine, vol. 25, no. 2, pp. 21–30,
2008.
[4] S. S. Chen, D. L. Donoho, and M. A. Saunders, “Atomic decomposition
by basis pursuit,” SIAM Journal of Scientific Computing, vol. 20, no. 1,
pp. 33–61, 1998.
[5] D. L. Donoho and M. Elad, “Maximal sparsity representation via
minimization,” Proceedings of the National Academy of Sciences, vol.
100, no. 5, pp. 2197–2202, 2003.
[6] J. A. Tropp and A. C. Gilbert, “Signal recovery from random measurements via orthogonal matching pursuit,” IEEE Transactions on
Information Theory, vol. 53, no. 12, pp. 4655–4666, 2007.
[7] D. Needella and J. A. Troppb, “Cosamp: Iterative signal recovery from
incomplete and inaccurate samplesstar,” Applied and Computational
Harmonic Analysis, vol. 26, no. 3, pp. 301–321, 2009.
[8] W. Dai, “Subspace pursuit for compressive sensing signal reconstruction,” IEEE Transactions on Information Theory, vol. 55, no. 5, pp.
2230–2249, 2009.
[9] C.F. Caiafa and A. Cichocki, “Computing sparse representations of
multidimensional signals using kronecker bases,” Neural Comput, vol.
25, no. 1, pp. 186 – 220, 2013.
[10] C.F. Caiafa and A. Cichocki, “Stable, robust and super fast reconstruction of tensors using multi-way projections,” IEEE Transactions on
Signal Processing, vol. 63, no. 3, pp. 780 – 793, 2015.
[11] L. Gan, “Block compressed sensing of natural images,” in Conf. on
Digital Signal Processing, 2007.
[12] S. Mun and J. E. Fowler, “Block compressed sensing of images using
directional transforms,” Proc. IEEE Int. Conf. Image Processing, pp.
3021–3024, 2009.
[13] A. Milzarek and M. Ulbrich, “A semismooth newton method with
multidimensional filter globalization for l1-optimization,” SIAM Journal
on Optimization, vol. 50, pp. 298–333, 2014.
[14] M. A. T. Figueiredo, R. D. Nowak, and S. J. Wright, “Gradient
projection for sparse reconstruction: Application to compressed sensing
and other inverse problems,” IEEE Journal of Selected Topics in Signal
Processing, vol. 1, no. 4, pp. 586–597, 2007.
[15] S. J. Wright, R. D. Nowak, and M. A. T. Figueiredo, “Sparse reconstruction by separable approximation,” Proceedings of IEEE International
Conference on Acoustics, Speech and Signal Processing, pp. 3373–3376,
2008.
[16] A. Beck and M. Teboulle, “A fast iterative shrinkage-thresholding
algorithm for linear inverse problems,” SIAM Journal on Imaging
Sciences, vol. 2, pp. 183–202, 2009.
[17] Z. Wen, W. Yin, D. Goldfarb, and Y. Zhang, “A fast algorithm for
sparse reconstruction based on shrinkage, subspace optimization, and
continuation,” SIAM Journal on Scientific Computing, vol. 32, pp. 1832–
1857, 2010.
[18] S. Mun and J. E. Fowler, “Motion-compensated compressed-sensing reconstruction for dynamic mri,” Proc. IEEE Int. Conf. Image Processing,
pp. 1006–1010, 2013.

217
16.6
7.8
5.2

218
37.9
17.3
10.0

219
79.4
36.1
22.5

220
173.2
78.3
49.1

221
359.2
177.6
109.2

222
740.6
382.9
241.6

[19] Y. Shi, X. Sun, J. Wang, and B. Yin, “Two dimensional synthesis sparse
model,” in IEEE ICME, 2013.
[20] Y. Rivenson and A. Stern, “Compressed imaging with a separable
sensing operator,” IEEE Signal Processing Letters, vol. 16, no. 6, pp.
449–452, 2009.
[21] Y. Rivenson and A. Ster, “Practical compressive sensing of large
images,” in 16th Int’l Conf. on Digital Signal Processing, 2009.
[22] E. H. Candes, J. Romberg, and T. Taio, “Robust uncertainty principles:
Exact signal reconstruction from highly incomplete frequency information,” IEEE Transactions on Information Theory, vol. 52, no. 2, pp.
489–507, 2006.
[23] T. T. Do, G. Lu, N. H. Nguyen, and T. D. Tran, “Fast and efficient compressive sensing using structurally random matrices,” IEEE Transations
on Signal Processing, vol. 60, pp. 139–154, 2012.
[24] N. Sidiropoulos and A. Kyrillidis, “Multi-way compressed sensing for
sparse low-rank tensors,” IEEE Signal Processing Letters, vol. 19, no.
11, pp. 757–760, 2012.
[25] Q. Li, D. Schonfeld, and S. Friedland, “Generalized tensor compressive
sensing,” IEEE International Conference on Multimedia and Expo, pp.
1–6, 2013.
[26] V. Cevher, S. Becker, and M. Schmidt, “Convex optimization for
big data: Scalable, randomized, and parallel algorithms for big data
analytics,” IEEE Signal Processing Magazine, vol. 31, no. 5, pp. 32–43,
2014.
[27] S. Becker, J. Bobin, and E. Candes, “Nesta: A fast and accurate firstorder method for sparse recovery,” SIAM Journal on Imaging Sciences,
vol. 4, pp. 1–39, 2011.
[28] Z. Wen, W. Yin, H. Zhang, and D. Goldfarb, “On the convergence of
an active-set method for ‘1 minimization,” Optimization Methods and
Software, vol. 27, pp. 1127–1146, 2012.
[29] M. F. Duarte, M. A. Davenport, D. Takhar, J. N. Laska, T. Sun, K. F.
Kelly, and R. G. Baraniuk, “Single-pixel imaging via compressive
sampling,” IEEE Signal Processing Magazine, vol. 25, no. 2, pp. 83–91,
2008.
[30] T. Blumensath and M. E. Davies, “Gradient pursuits,” IEEE Transactions on Signal Processing, vol. 56, no. 6, pp. 2370–2382, 2008.
[31] R. Chartrand and Yin Wotao, “Iteratively reweighted algorithms for
compressive sensing,” in Proceedings of IEEE International Conference
on Acoustics, Speech and Signal Processing, 2008, pp. 3869–3872.
[32] A. Tewari P. Jain and I. S. Dhillon, “Orthogonal matching pursuit with
replacement,” in NIPS, 2011.
[33] A. van der Sluis and H. A. van der Vorst, “The rate of convergence
of conjugate gradients,” Numerische Mathematik, vol. 48, pp. 543–560,
1986.
[34] E. F. Kaasschieter, “Preconditioned conjugate gradients for solving
singular system,” Journal of Computational and Applied Mathematics,
vol. 24, pp. 265–275, 1988.
[35] C. T. Kelley, “Iterative methods for linear and nonlinear equations,” in
Frontiers in Applied Mathematics, 1995.

10

5

10

4

Computational Cost (Sec)

10

3

10

2

10

1

10

0

10

M−OMP
CG−OMP
O−OMP

−1

10

−2

10

11

12

13

14

15

16

17

18

log N
2

(a) OMP
4

Computational Cost (Sec)

10

3

10

2

10

1

10

M−OMPR
CG−OMPR
O−OMPR

0

10

−1

10

10

12

14

16

18

20

22

24

log N
2

(b) OMPR
4

Computational Cost (Sec)

10

3

10

2

10

1

10

M−SP
CG−SP
O−SP

0

10

−1

10

10

12

14

16

18

20

22

24

log N
2

(c) SP
Fig. 3. Computation cost vs. signal length for OMP (a), OMPR (b), and SP
, and K = M
(p = 0.0625). The precision
(c), respectively, under M = N
4
4
of CG was set to be exact.

[36] S. Hawe, M. Seibert, and M. Kleinsteuber, “Separable dictionary
learning,” IEEE CVPR, pp. 438–445, 2013.
[37] S.-H. Hsieh, C.-S. Lu, and S.-C. Pei, “2d sparse dictionary learning
via tensor decomposition,” IEEE Global Conference on Signal and
Information Processing, pp. 492–496, 2014.
[38] H. Mohimani, M. Babaie-Zadeh, and C. Jutten, “A fast approach for
overcomplete sparse decomposition based on smoothed l0 norm,” IEEE
Transactions on Signal Processing, vol. 57, pp. 289–301, 2009.

