A Tutorial on Sparse Signal Acquisition and Recovery with Graphical Models
Volkan Cevher, Piotr Indyk, Lawrence Carin, Richard G. Baraniuk
I. I NTRODUCTION Many applications in digital signal processing, machine learning, and communications feature a linear regression problem in which unknown data points, hidden variables or codewords are projected into a lower dimensional space via

y = x + n.

In the signal processing context, we refer to x RN as the signal, y RM as measurements with
M < N , RM N as the measurement matrix, and n RM as the noise. The measurement matrix is a matrix with random entries in data streaming, an overcomplete dictionary of features in sparse
Bayesian learning, or a code matrix in communications [13]. Extracting x from y in (1) is ill-posed in general since M < N and the measurement matrix hence has a nontrivial null space; given any vector v in this null space, x + v denes a solution that produces the same observations y. Additional information is therefore necessary to distinguish the true x among the innitely many possible solutions [1, 2, 4, 5]. It is now well-known that sparse representations can provide crucial prior information in this dimensionality reduction; we therefore also refer to the problem of determining x in this particular setting as the sparse signal recovery. A signal x has a sparse representation x = in a basis RN N when K
N coefcients of can exactly represent
or well-approximate the signal x. Inspired by communications, coding and information theory problems, we often refer to the application of on x as encoding of x; the sparse signal recovery problem is then concerned with the decoding of x from y in the presence of noise. In the sequel, we assume the canonical sparsity basis, = I without loss of generality. The sparse signal recovery problem has been the subject of extensive research over the last few decades in several different research communities, including applied mathematics, statistics, and theoretical computer science [13, 6]. The goal of this research has been to obtain higher compression rates; stable recovery schemes; low encoding, update and decoding times; analytical recovery bounds; and resilience to noise. The momentum behind this research is well-justied: underdetermined linear regression problems in tandem with sparse representations underlie the paradigm of signal compression
and denoising in signal processing, the tractability and generalization of learning algorithms in machine learning, the stable embedding and decoding properties of codes in information theory, the effectiveness of data streaming algorithms in theoretical computer science, and neuronal information processing and interactions in computational neuroscience. An application du jour of the sparse signal recovery problem is compressive sensing (CS), which integrates the sparse representations with two other key aspects of the linear dimensionality reduction: information preserving projections and tractable recovery algorithms [1, 2, 46]. In CS, sparse signals are represented by a union of the
, K -dimensional subspaces, denoted as x K. We call the
set of indices corresponding to the nonzero entries the support of x. While the matrix is rank decient, it can be shown to preserve the information in sparse signals if it satises the so-called restricted isometry property (RIP). Intriguingly, a large class of random matrices have the RIP with high probability. Todays state-of-the-art CS systems can robustly and provably recover K -sparse signals from just M = O(K log(N/K)) noisy measurements using sparsity-seeking, polynomial-time optimization solvers or greedy algorithms. When x is compressible, in that it can be closely approximated as K sparse, then from the measurements y , CS can recover a close approximation to x. In this manner we can achieve sub-Nyquist signal acquisition, which requires uniform sampling rates at least two times faster than the signals Fourier bandwidth to preserve information. While such measurement rates based on sparsity are impressive and have the potential to impact a broad set of streaming, coding, and learning problems, sparsity is merely a rst-order description of signal structure; in many applications we have considerably more a priori information that previous approaches to CS fail to exploit. In particular, modern signal, image, and video coders directly exploit the fact that even compressible signal coefcients often sport a strong additional structure in the support of the signicant coefcients. For instance, the image compression standard JPEG2000 does not only use the fact that most of the wavelet coefcients of a natural image are small. Rather, it also exploits the fact that the values and locations of the large coefcients have a particular structure that is characteristic of natural images. Coding this structure using an appropriate model enables JPEG2000 and other similar algorithms to compress images close to the maximum amount possible, and signicantly better than a naive coder that just assigns bits to each large coefcient independently [1]. By exploiting a priori information on coefcient structure in addition to signal sparsity, we can make CS better, stronger, and faster. The particular approach we will focus in this tutorial is based on graphical models (GM) [3, 710]. As we will discover, GMs are not only useful for representing the prior information on x, but also lay the foundations for new kinds of measurement systems. GMs enable 2

Inference algorithms on graphical models typically exploit the factorization properties (2) and (3) of probability distributions. By manipulating the intermediate factors, it is often possible to compute the likelihood of a particular z value in an efcient, distributed manner. This strategy is exploited by the sum-product algorithm (also known as belief propagation) and max-product (also known as min-sum) for tree-structured DAGs with rigorous estimation guarantees when the random variables live in a discrete probability space [7]. These algorithms iteratively pass statistical information, denoted as messages, among neighboring vertices and converge in nite number of steps. Such guarantees, unfortunately, do not extend to arbitrary DAGs; inference in DAGs is typically NP-hard [11]. The sum-product and max-product algorithms are routinely applied to graphical models with cycles, leading to loopy belief propagation methods, even if their convergence and correctness guarantees for DAGs no longer hold in general [3, 11, 13]. Although there are certain local optimality guarantees associated with the xed points of loopy belief propagation, there are a number of natural inference problems arising in various applications in which loopy belief propagation either fails to converge, or provides poor results. Loopy belief propagation is therefore an approximate inference method. Surprisingly, stateof-the-art algorithms for decoding certain kinds of error-correcting codes are equivalent to loopy belief propagation. We revisit this topic in Section IV to discover a crucial role sparse representations play in providing theoretical guarantees for such algorithms. Approximate inference is also necessary in cases where the random variables are drawn from a continuous space, since corresponding marginal integrals needed for implementing Bayes rule cannot be analytically performed. Monte Carlo Markov chain sampling methods, such as importance sampling, Metropolis-Hastings and Gibbs sampling, provide computational means of doing approximate inference. The idea is to represent the pdf by a discrete set of samples, which is carefully weighted by some evidence likelihood so that the inference is consistent, with guarantees typically improving with larger sample sizes. See [3] for further background. Yet another approximate inference method is based on the calculus of variations, also known as variational Bayes (VB) [3, 8]. The quintessential VB example is the mean-eld approximation, which exploits the law of large numbers to approximate large sums of random variables by their means. In particular, the mean-eld approximation essentially decouples all the vertices in the graph, and then introduces a parameter, called a variational parameter, for each vertex. It then iteratively updates these variational parameters so as to minimize the cross-entropy between the approximate and true probability distributions. Updating the variational parameters then facilitates inference in an efcient and principled manner, often also providing bounds on the marginal likelihood to be calculated. Another special case of 5

the VB approach is the well-known expectation-maximization (EM) algorithm for MAP and maximum likelihood estimation, which has been extremely successful in a number of applications [3, 8, 14]. B. Compressive sensing (CS) 1) Sparse signal representations: No linear nonadaptive dimensionality reducing can preserve all of the information in all signals. Hence, researchers restrict the application of to not arbitrary signals
x RN but rather to some subset of RN , and in particular the set of sparse and compressible signals. A
sparsity/compressibility prior is then exploited as a tie-breaker to distinguish the true signal among the innitely many possible solutions of (1). While sparse signals live in K , the coefcients of a compressible signal x, when sorted in order of decreasing magnitude, decay according to the following power law:
xI(i) R i1/r , i = 1,. , N,
where I indexes the sorted coefcients, and R > 0 and r > 0 are constants. In the CS setting, we call a signal compressible when r 1. Thanks to the power-law decay of their coefcients, compressible signals are well-approximated by K -sparse signals in an appropriate norm. For instance, for all r < 1 and 0 <

N p i=1 |xi |

1, x xK
holds independent of N for any K (r/ ) 1r , where x

is the

p -norm,

xK = arg min x

is the best K -sparse approximation

of x (p 1), and x

is a pseudo-norm that counts the number of nonzeros of x [15].
2) Information preserving projections: The sparsity or compressibility of x is not sufcient alone for distinguishing x among all possible solutions to (1). The projection matrix must also work in tandem with the signal priors so that recovery algorithms can correctly identify the true signal. The restricted isometry property (RIP) assumption on achieves this by requiring to approximately preserve the distances between all signal pairs in the sparse signal set [5]. More formally, an M N matrix has the K -restricted isometry property (K -RIP) with constant

< 1 if, for all x K ,

x 2. 2
An alternative property is the restricted strong convexity (RSC) assumption, which is motivated by convex optimization arguments [16]. In general, the RSC assumption has an explicit dependence on the recovery algorithms objective function. For instance, if the recovery algorithms objective is to minimize the measurement error (e.g., y x 2 ), RSC requires t x 2
to be strictly positive for all x K. In
different contexts, other conditions on are also employed with varying levels of restrictions, such as null space property, spark, unique representation property, etc [6]. 6
While checking whether a measurement matrix satises the K -RIP, RSC, etc. has running time exponential in K and N , random matrices with iid subgaussian entries work with high probability provided M = O(K log(N/K)). Random matrices also have a so-called universality property in that, for any choice of orthonormal basis matrix , also has the K -RIP with high probability. This is useful when the signal is sparse in some basis other than the canonical basis. 3) Tractable recovery algorithms: To recover the signal x from y in (1), we exploit our a priori knowledge of its sparsity or compressibility. For example, to recover strictly sparse signals when there is no measurement noise, we can seek the sparsest x that agrees with the measurements y. While this optimization can recover a K -sparse signal from just M = 2K compressive measurements, it is not only a combinatorial, NP-hard problem, but also is not stable in the presence of noise [1]. Tractable recovery algorithms rely on conditions such as the RIP and RSC for stability and correctness and therefore require at least M = O(K log(N/K)) measurements. They can be grouped into two main camps: convex optimization and greedy approximation. The rst camp relies on a convex relaxation of seeking sparse solutions:

x = arg min x x

minimization as

s.t. y = x.

This optimization problem is known as basis pursuit and corresponds to a linear program that can be solved in polynomial time [1, 4, 5]. Adaptations to deal with additive noise in (1) have also been proposed; examples include basis pursuit with denoising (BPDN), and the least absolute shrinkage and selection operator (LASSO). The second camp nds the sparsest x agreeing with the measurements y through an iterative, greedy search over the coefcients. Example algorithms include iterative hard thresholding (IHT), compressive sampling matching pursuit (CoSaMP), and Subspace Pursuit (SP) [1]. Currently, convex optimization obtains the best recovery performance in theory, while its greedy counterparts, such as IHT, CoSaMP, and SP, offer desirable computational trade-offs, e.g., O N log2 N vs. O M 2 N 1.5 (interior point methods) [1]. Interestingly, algorithms in both camps have similar theoretical recovery guarantees from M = O(K log(N/K)) when the measurement matrix satises K -RIP with an algorithm dependent 2K (e.g., 2K < for basis pursuit [5]):

C1 K 1/2 x xK

+ C2 n 2 ,

which we refer to as an

guarantee, where the subscripts are matched to the norms adjacent to the is known as the irrecoverable energy.
inequality. In (7), x is the algorithm output, C1 and C2 are algorithm-dependent constants, and xK is the signal-dependent best K -sparse approximation. The term xxK 7
When the signals are compressible, we can never recover them fully under dimensionality reduction; we can only approximately recover them and only when r 1 in (4). III. G RAPHICAL MODELS FOR STRUCTURED SPARSITY While CS has the potential to revolutionize data acquisition, encoding, and processing in a number of applications, much work remains before it is ready for real-world deployment. In particular, for CS to truly live up its name, it is crucial that the theory and practice leverage concepts from state-of-the-art transform compression algorithms. In virtually all such algorithms, the key ingredient is a signal model for not just the coefcient sparsity but also the coefcient structure. As we saw in Section II, graphical models provide natural framework for capturing such dependencies among sparse signal coefcients. Hence, we review below some of the graphical models that are relevant for structured sparsity and the algorithms that result for sparse signal recovery. A. Sparsity priors Our rst departure from the standard, deterministic CS approach is to adopt a more general, Bayesian viewpoint as several others have pursued [1, 17, 18]. A probabilistic sparsity prior enforces sparsity or compressibility in each coefcient by heavily weighting zero or small-magnitude values in the probability distribution, while allowing a large range of values with low probabilities to obtain a small number of large coefcients. While this is expected in principle, we also require a generative aspect for such priors in order to build up structured sparsity models; that is, their statistical realizations result in sparse or compressible signals. This, in turn, is crucial for statistical consistency in high-dimensional scalings of the sparse signal recovery problem. Two-state mixture models are the canonical iid models for generating strictly sparse signals. Among mixture models, the spike-and-slab prior stipulates a density with a spike at zero surrounded symmetrically by a uniform distribution with specied boundaries. The mixture probability of the spike at zero approximately determines the percentage of zero coefcients of the signal, and the uniform distribution establishes the distribution of the signal coefcients on the signal support. To model compressible signals, we can use compressible priors, whose statistical realizations exhibit the power-law decay in (4) [15]. For algorithmic recovery guarantees for CS, we must have r 1 for such priors (c.f. (4)). For instance, Table III-A demonstrates that N -sample iid realizations of generalized Pareto, Students t, Frechet, and log-logistics distributions (each parameterized by a shape parameter

q > 0 and a scale parameter > 0) are compressible with parameter r = q. There also exist non-iid
compressible priors; multivariate Lomax distribution provides an elementary example whose pdf is given 8
TABLE I E XAMPLE DISTRIBUTIONS AND THE COMPRESSIBILITY PARAMETERS OF THEIR IID REALIZATIONS
Distribution Generalized Pareto Students t Fr chet e Log-Logistic Laplacian

pdf 1+

|x| (q+1) (q+1)/2 (x/)q

R N 1/q

2((q+1)/2) q(q/2) 1/q

r q N 1/q q q q log N

((q+1)/2) 2(q/2)

(q/) (x/)

1/q 1/q
(q/)(x/)q1 [1+(x/)q ]|x|/ e 2

log N

by MLD(x; q, ) 1 +

1 |xi |

[15]. The compressibility parameter MLD is simply r = 1,
irrespective of its shape parameter. To illustrate how we can exploit probabilistic priors in sparse signal recovery, we focus on two compressible priors on x: MLD and Students t. While MLD is relatively unknown, the Students t distribution has enjoyed tremendous attention in many fundamental problems, such as statistical modeling of natural images, relevance vector machines (RVM) and automatic relevance determination, and dictionary learning. We consider the case when there is no noise; the observations are then given by y = x with
M = O(K log(N/K)), which has innitely many solutions for x, as discussed in Section II. We know
that such projections preserve the information of the largest K -coefcients of signal realizations. We will assume that the parameters of MLD and the Students t priors are matched (i.e., the shape parameter of Students t is 1) so that the statistical realizations of both priors have the same decay prole. 1) MLD and basis pursuit: We rst exploit the MLD likelihood function for sparse signal recovery. For instance, when we ask for the solution that maximizes the MLD likelihood given y for i = , it is easy to see that we obtain the basis pursuit algorithm formulation in (6). In contrast, the conventional probabilistic motivation for the basis pursuit algorithm assumes that x is iid Laplacian, e.g., f (x) exp ( x 1 /) for some R+. We then face the optimization problem (6), when we seek the most likely signal from the Laplacian prior given y. Unfortunately, iid Laplacian realizations cannot be approximated as sparse since their compressibility parameter is r = log N , which is greater than 1 in general (c.f. Table III-A); it can be shown that x xK 1 x 1 requires K (1 )N with high probability [15]. Hence, while the
minimization correctly recovers the sparse vectors from M = O(K log(N/K)), it seems
that it does not correspond to an iid Laplacian prior on the sparse signal coefcients in a consistent Bayesian framework and may correspond to the MLD prior as this example demonstrates. 2) Students t and iterative reweighted least squares: We now exploit the Students t likelihood function 9

for sparse signal recovery, similar to MLD likelihood function above:
x = max f (x) = min x x log 1 + 2 xi2 , s.t. y = x.
Unfortunately, (8) is a non-convex problem. However, we can circumvent the non-convexity in (8) using a simple variational Bayes idea where we iteratively obtain a tractable upperbound on the log-term in (8) using the following inequality: u, v (0, ), log u log v + u/v 1. After some straightforward calculus, we obtain the iterative algorithm below, indexed by k , where x{k} is the k -th iteration estimate (x{0} = 0):

x{k} = min x

2 wi,{k} xi2 , s.t. y = x ; where wi,{k} = 2 + xi,{k} i 1
The decoding scheme in (9) is well-known as the iterative reweighted least squares (IRLS), where each iteration has an analytical solution [19]. Both priors in the above algorithms are trying to approximate their signal realizations, which have the same decay prole, and yet result in radically different algorithms. Both schemes have provable recovery guarantees; however, they differ in terms of their computational costs: O(M 2 N ) (IRLS) vs. O(M 2 N 1.5 ) (BP). B. Structured sparsity via GMs Compressible priors provide a natural launching point for incorporating further structure among the sparse signal coefcients. Here, by structure, we specically mean the probabilistic dependence and independence relationships among sparse signal coefcients as summarized by directed and undirected graphs. Such relationships are abundant in many diverse problems, such as speech recognition, computer vision, decoding of low-density parity-check codes, modeling of gene regulatory networks, gene nding and diagnosis of diseases. Graphical models also provide powerful tools for automatically learning these interactions from data. While learning GMs is an active research area, it is beyond the scope of this paper. 1) Tree graphs: Wavelet transforms sparsify piecewise smooth phenomena, including many natural and manmade signals and images. However, the signicant discrete wavelet transform (DWT) coefcients do not occur in arbitrary positions for such signals. Instead they exhibit a characteristic signal-dependent structure, with the large and small wavelet coefcients clustering along the branches of the wavelet tree. We highlight the persistency of the parent child relationship for several large DWT coefcients on the wavelet image in Figure 2. Hidden Markov trees (HMT) models and Gaussian scale mixtures (GSMs) on steerable pyramids succinctly and accurately captures this statistical structure [12, 20, 21]. The HMT models the probability 10

Original

Compressed

DWT-compressible image

DWT coefcients

Fig. 2.

10% of DWT coefcents

Shufed DWT coeff.

Most of the discrete wavelet transform (DWT) coefcients of natural images are small (blue in the bottom
row corresponds to zero amplitude). However, the signicant coefcients appear clustered on the wavelet trees [12, 20], as illustrated by the solid red lines. The dependence structure of the coefcients is crucial; by randomly shufing the DWT coefcients of the original image, we obtain a reconstruction that does not even remotely resemble a natural image, let alone the original.
density function of each wavelet coefcient as a mixture density with a hidden binary state that determines whether the coefcient is large or small. Wavelet coefcient persistence along the tree branches is captured by a tree-based Markov model that correlates the states of parent and children coefcients. GSMs exploit the Students t prior or its variants. Both models have been successfully applied to improve the performance of denoising, classication, and segmentation algorithms for wavelet sparse signals. For an explanation of the EM algorithm applied to GSMs, see the original EM paper [14]. Another Bayesian approach leverages the clustering of the DWT coefcients using the spike-and-slab prior along with VB inference. The main idea is to construct a tree-structured, conjugate-exponential family sparse model for wavelet coefcients where the mixture probability of the spike at an individual sparse coefcient on the wavelet tree is controlled by the size of its parent on the tree. Then, the VB updates can be efciently calculated with convergence guarantees. For instance, [17] demonstrates that the VB approach simultaneously improves both the recovery speed and performance over the state-of-the-art greedy and convex optimization based CS recovery approaches, discussed in Section II. Moreover, as 11

target

sparse recovery

LAMP recovery

Fig. 3.

LaMP Iter. #1

LaMP Iter. #2

[x, s] = arg max x ,s ij si sj +

i si + log(p(xi |si ))

1 y x 2 2
For the compressibility of the structured sparse signal x as signal dimensions vary, the pdfs p(xi |si ) should be chosen to result in a compressible prior for the signal coefcients, such as Students t, with proper scale parameters to differentiate the small and large coefcients. Although the objective function is non-convex, a local optimum can be efciently obtained via variational methods and alternating minimization techniques over the support and the signal coefcients. In the context of sparse signal recovery, recent work [18] exploits the Ising model and demonstrates that it naturally motivates a greedy search algorithm for background subtraction with clustered sparsity, dubbed Lattice matching pursuit (LAMP). Enabled by the Ising model, the LAMP algorithm (i) converges signicantly faster due to the reduction in the search space and (ii) recovers sparse signals from a much smaller number of measurements without sacricing stability. Similar optimization formulations can also be seen in applications involving face recognition in the presence of occlusions, where the occlusions are not only sparse but also contiguous. To solve (11), the LAMP algorithm relies on an approach inspired by the RIP assumption on the measurement matrices. The key observation for this iterative recovery scheme is that when the sampling matrix has K -RIP with constant
1 in (5), then the vector b = T y can serve as a rough
approximation of the original signal x. In particular, the largest K entries of b point toward the largest
K entries of the K -sparse signal x. Then, given the signal coefcient estimates, LAMP uses graph
cuts to obtain the MAP estimates of the latent support variables s. By clever book-keeping of the data 13

Fig. 4.

Incomplete

Lin. interp. (27.6dB)

IBP (35.2dB)

Learned dictionary

Natural images exhibit signicant self similarities that can be leveraged using the Indian buffet processes (IBP) in sparse
signal recovery even if a large number of pixels are missing. The IBP mechanism automatically infers a dictionary, in which the image patches are sparse, and the composition of these patches on the image to signicantly improve the recovery performance over linear interpolation.

(a) A bipartite graph. Fig. 5.

(b) Expansion.

An example of a bipartite graphs G = (A, B, E), over the left set A of size N and the right set B of size M. A
graph is an expander if any small subset S of A has many neighbors (S) in B.
Hence, in this section, we focus on another use of graphical models, to design high-performance measurement matrices that are sparse. Specically, we consider matrices whose entries are mostly equal to zero, while the non-zero entries are equal to 1.1 Each such matrix can be interpreted as bi-partite graph
G = (A, B, E) between a set A = {1. N } of N nodes, and a set B = {1. M } of M nodes: the
edge i j is present in the graph if and only if j,i = 1. Note that the nodes in A correspond to the coordinates of the signal x, while the nodes in B correspond to coordinates in the measurement vector
y. See Figure 5 for an illustration.
Sparse measurement matrices have several desirable features: (i) matrix-vector products during the encoding of the vector x into x can be performed very efciently, in time proportional to the number of non-zeros in ; (ii) the measurement vector x can be quickly updated if one or a few of its coordinates are modied crucial for processing massive data streams; (iii) the recovery process is quite efcient as well, since it relies on the matrix-vector product as a subroutine. Moreover, the graphical interpretation of such matrices enabled designing recovery algorithms using the belief propagation approach, leading to highly accurate recovery methods. Because of these reasons, various forms of sparse recovery using sparse matrices has been recently a subject of extensive research; see the recent survey [2] and the references therein. To complement the above, we also point out some disadvantages of sparse matrices. One of them is that they are directly applicable only to the case where the signal x is approximately sparse in the
Some of the constructions (e.g., [24]) allow the non-zero values to be equal to 1 as well.
canonical basis, i.e., x = (see the introduction for the notation). If the signal (as it is often the case) is sparse only after applying a linear transformation , then the actual measurement matrix is equal to

s = K. Let x be a K -sparse vector , and let S be a set of indices of non-zero entries in x. Then we
observe that for most nodes i S , only a small fraction of their neighbors (i) collide with any other non-zero entry. Such entries xi can be thus recovered correctly by the median-based procedure outlined earlier. To recover all entries, we can use an iterative renement approach, similar to the one from the previous section, and inspired by the bit-ipping algorithm for decoding low-density parity check codes. Consider rst the case where x is exactly K -sparse. The algorithm [30] starts by setting an initial approximation
x to 0 and iteratively renes x in order to achieve x = y. In each step, it tries to reduce x y 0 ,
by nding a pair (i, g) such that incrementing xi by g reduces the
difference. It is then shown that
if the graph is an (O(K), (1 )D)-expander for a sufciently small value of , then x converges to x in O(K) iterations. The running time of the algorithm is dominated by the preprocessing step, which takes time O(N D), after which each iteration can be performed in O(log N ) time or, in some cases, even faster. Since the graph G is an expander, it follows that the number of measurements is
M = O(KD) = O(K log(N/K)).
When x is not exactly K -sparse, then the
norm is no longer a suitable measure of progress, and
we need to use a more rened approach. To this end, we rst observe that the graph-theoretic notion of expansion has a natural and useful geometric interpretation. Specically, consider the following variant of the K -RIP in (5): we say that an M N matrix has the K -restricted isometry property in the norm (K -RIP1) with constant , if for all K -sparse vectors x, we have

x 1 (1 ) x

It has been shown in [31] that if is a matrix underlying a (K, D(1 /2))-expander, then /D satises K -RIP1 with constant. As a result, both
minimization and variants of the iterative algorithms
described in Section II can be used for sparse matrices. Unlike the methods using the standard RIP, the algorithms produce x which satises a somewhat weaker guarantee of the form

C x xK

Perhaps surprisingly, the simplest of the algorithms, called Sequential Sparse Matching Pursuit [32], is similar to the aforementioned algorithm of [30]. There are two key differences though. Firstly, each iteration reduces not the

error x y 0 , but the

error x y 1. This is because by the RIP1

property, the error x x

is small when x y
is small. Second, in order to be able to apply the
RIP1 property, we need to ensure that the vector x continues to be O(K)-sparse after we perform the updates. Since this property might cease to be true after some number of steps, we need to periodically re-sparsify the vector x by setting to zero all but the K largest (in absolute value) entries of x. This re-sparsication step is a standard tool in the design of iterative recovery algorithms for general vectors
See the survey [2] for a more detailed description of the algorithms for sparse matrices. V. C ONCLUSIONS A great deal of theoretic and algorithmic research has revolved around sparsity view of signals over the last decade to characterize new, sub-Nyquist sampling limits as well as tractable algorithms for signal recovery from dimensionality reduced measurements. Despite the promising advances made, real life applications require more realistic signal models that can capture underlying, application dependent order of sparse coefcients, better sampling matrices with information preserving properties that can be implemented in practical systems, and ever faster algorithms with provable recovery guarantees for real-time operation. On this front, we have seen that graphical models (GM) are emerging to effectively address the core of many of these desiderata. GMs provide a broad scaffold for automatically encoding the probabilistic 19
dependencies of sparse coefcients for sparse signal recovery. By exploiting the GM structure of signals beyond simple sparsity, we can radically reduce the number of measurements, increase noise robustness, and decrease recovery artifacts in signal acquisition. GMs are instrumental in constructing measurement matrices based on expander graphs. These matrices not only stably embed sparse signals into lower dimensions but also lead to faster recovery algorithms with rigorous guarantees. Moreover, the GMbased inference tools, such as variational methods, can estimate a posterior distribution for the sparse signal coefcients, providing condence bounds that are critical in many applications. To date, the sparse signal acquisition and recovery problemssurprisinglyhave been studied largely in isolation. Real progress in efcient signal recovery, processing and analysis requires that we unify probabilistic, structured sparsity models with sparse measurement matrices to simultaneously reduce sampling requirements and the computational complexity of recovery without compromising the recovery guarantees. This will in turn entail investigation of streaming algorithms, coding theory, and learning theory with a common, connecting element, which we expect to be graphical models. R EFERENCES

[1] R. G. Baraniuk, V. Cevher, and M. Wakin, Low-dimensional models for dimensionality reduction and signal recovery: A geometric perspective, Proceedings of the IEEE, 2010. [2] A. Gilbert and P. Indyk, Sparse recovery using sparse matrices, Proceedings of IEEE, 2010. [3] C. M. Bishop, Pattern recognition and machine learning. Springer, 2006. [4] D. L. Donoho, Compressed sensing, IEEE Trans. Info. Theory, vol. 52, pp. 12891306, Sept. 2006. [5] E. J. Cand` s, Compressive sampling, in Proc. International Congress of Mathematicians, vol. 3, (Madrid, Spain), e pp. 14331452, 2006. [6] A. Cohen, W. Dahmen, and R. DeVore, Compressed sensing and best k-term approximation, American Mathematical Society, vol. 22, no. 1, pp. 211231, 2009. [7] S. L. Lauritzen, Graphical Models. Oxford University Press, 1996. [8] M. I. Jordan, Z. Ghahramani, T. S. Jaakkola, and L. K. Saul, An introduction to variational methods for graphical models, Machine Learning, vol. 37, no. 2, pp. 183233, 1999. [9] A. Willsky, Multiresolution Markov models for signal and image processing, Proceedings of the IEEE, vol. 90, no. 8, pp. 13961458, 2002. [10] J. Pearl, Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference. Morgan Kaufmann Publishers, 1988. [11] D. Koller and N. Friedman, Probabilistic Graphical Models: Principles and Techniques. The MIT Press, 2009. [12] M. S. Crouse, R. D. Nowak, and R. G. Baraniuk, Wavelet-based statistical signal processing using Hidden Markov Models, IEEE Trans. Signal Processing, vol. 46, pp. 886902, Apr. 1998. [13] J. S. Yedidia, W. T. Freeman, and Y. Weiss, Generalized belief propagation, Advances in Neural Information Processing Systems, vol. 13, pp. 689695, 2001. [14] A. P. Dempster, N. M. Laird, and D. B. Rubin, Maximum likelihood from incomplete data via the EM algorithm, Journal of the Royal Statistical Society. Series B (Methodological), vol. 39, no. 1, pp. 138, 1977.
[15] V. Cevher, Learning with compressible priors, in NIPS, (Vancouver, B.C., Canada), 712 December 2008. [16] S. Negahban and M. J. Wainwright, Estimation of (near) low-rank matrices with noise and high-dimensional scaling, Arxiv preprint arXiv:0912.5100, 2009. [17] L. He and L. Carin, Exploiting structure in wavelet-based Bayesian compressive sensing, 2008. Preprint. Available at http://people.ee.duke.edu/ lcarin/Papers.html. [18] V. Cevher, M. F. Duarte, C. Hegde, and R. G. Baraniuk, Sparse signal recovery using Markov random elds, in Neural Information Processing Systems (NIPS), (Vancouver, B.C., Canada), 811 December 2008. [19] I. Daubechies, R. DeVore, M. Fornasier, and C. S. Gunturk, Iteratively Reweighted Least Squares Minimization for Sparse Recovery, Communications on Pure and Applied Mathematics, vol. 63, pp. 00010038, 2010. [20] J. K. Romberg, H. Choi, and R. G. Baraniuk, Bayesian tree-structured image modeling using wavelet-domain Hidden Markov Models, IEEE Trans. Image Processing, vol. 10, pp. 10561068, July 2001. [21] M. J. Wainwright and E. P. Simoncelli, Scale mixtures of Gaussians and the statistics of natural images, in Neural Information Processing Systems (NIPS) (S. A. Solla, T. K. Leen, and K.-R. M ller, eds.), vol. 12, (Cambridge, MA), u pp. 855861, MIT Press, Dec. 2000. [22] M. Zhou, H. Chen, J. Paisley, L. Ren, G. Sapiro, and L. Carin, Non-parametric bayesian dictionary learning for sparse image representations, in Neural Information Processing Systems (NIPS), 2009. [23] E. B. Sudderth and M. I. Jordan, Shared segmentation of natural scenes using dependent Pitman-Yor processes, in Advances in Neural Information Processing Systems, vol. 21, 2009. [24] M. Charikar, K. Chen, and M. Farach-Colton, Finding frequent items in data streams, Proceedings of International Colloquium on Automata, Languages and Programming (ICALP), 2002. [25] G. Cormode and S. Muthukrishnan, Improved data stream summaries: The count-min sketch and its applications, Proceedings of Latin American Theoretical Informatics Symposium (LATIN), 2004. [26] D. Baron, S. Sarvotham, and R. G. Baraniuk, Bayesian compressive sensing via belief propagation, to appear in IEEE Transactions on Signal Processing, 2010. [27] Y. Lu, A. Montanari, B. Prabhakar, S. Dharmapurikar, and A. Kabbani, Counter braids: a novel counter architecture for per-ow measurement, in SIGMETRICS 08: Proceedings of the 2008 ACM SIGMETRICS international conference on Measurement and modeling of computer systems, (New York, NY, USA), pp. 121132, ACM, 2008. [28] D. Donoho, A. Maleki, and A. Montanari, Message passing algorithms for compressed sensing, Proceedings of the National Academy of Sciences, 2009. [29] C. Estan and G. Varghese, New directions in trafc measurement and accounting: Focusing on the elephants, ignoring the mice, ACM Transactions on Computer Systems, 2003. [30] S. Jafarpour, W. Xu, B. Hassibi, and A. R. Calderbank, Efcient and robust compressed sensing using high-quality expander graphs, IEEE Transactions on Information Theory, vol. 33(9), 2009. [31] R. Berinde, A. Gilbert, P. Indyk, H. Karloff, and M. Strauss, Combining geometry and combinatorics: A unied approach to sparse signal recovery, in Proc. Allerton Conf. Communication, Control, and Computing, 2008. [32] R. Berinde and P. Indyk, Sequential sparse matching pursuit, Proceedings of Allerton Conference on Communication, Control, and Computing, 2009.