Distributed Compressive Sensing

Introduction

A core tenet of signal processing and information theory is that signals, images, and other data often contain some type of structure that enables intelligent representation and processing. The notion of structure has been characterized and exploited in a variety of ways for a variety of purposes. In this paper, we focus on exploiting signal correlations for the purpose of compression.

Current state-of-the-art compression algorithms employ a decorrelating transform such as an exact or approximate Karhunen-Loève transform (KLT) to compact a correlated signal’s energy into just a few essential coefficients . Such transform coders exploit the fact that many signals have a sparse representation in terms of some basis, meaning that a small number $K$ of adaptively chosen transform coefficients can be transmitted or stored rather than $N\gg K$ signal samples. For example, smooth signals are sparse in the Fourier basis, and piecewise smooth signals are sparse in a wavelet basis ; the commercial coding standards MP3 , JPEG , and JPEG2000 directly exploit this sparsity.

While the theory and practice of compression have been well developed for individual signals, distributed sensing applications involve multiple signals, for which there has been less progress. Such settings are motivated by the proliferation of complex, multi-signal acquisition architectures, such as acoustic and RF sensor arrays, as well as sensor networks. These architectures sometimes involve battery-powered devices, which restrict the communication energy, and high aggregate data rates, limiting bandwidth availability; both factors make the reduction of communication critical.

Fortunately, since the sensors presumably observe related phenomena, the ensemble of signals they acquire can be expected to possess some joint structure, or inter-signal correlation, in addition to the intra-signal correlation within each individual sensor’s measurements. In such settings, distributed source coding that exploits both intra- and inter-signal correlations might allow the network to save on the communication costs involved in exporting the ensemble of signals to the collection point . A number of distributed coding algorithms have been developed that involve collaboration amongst the sensors . Note, however, that any collaboration involves some amount of inter-sensor communication overhead.

In the Slepian-Wolf framework for lossless distributed coding , the availability of correlated side information at the decoder (collection point) enables each sensor node to communicate losslessly at its conditional entropy rate rather than at its individual entropy rate, as long as the sum rate exceeds the joint entropy rate. Slepian-Wolf coding has the distinct advantage that the sensors need not collaborate while encoding their measurements, which saves valuable communication overhead. Unfortunately, however, most existing coding algorithms exploit only inter-signal correlations and not intra-signal correlations. To date there has been only limited progress on distributed coding of so-called “sources with memory.” The direct implementation for sources with memory would require huge lookup tables . Furthermore, approaches combining pre- or post-processing of the data to remove intra-signal correlations combined with Slepian-Wolf coding for the inter-signal correlations appear to have limited applicability, because such processing would alter the data in a way that is unknown to other nodes. Finally, although recent papers provide compression of spatially correlated sources with memory, the solution is specific to lossless distributed compression and cannot be readily extended to lossy compression settings. We conclude that the design of constructive techniques for distributed coding of sources with both intra- and inter-signal correlation is a challenging problem with many potential applications.

2 Compressive sensing (CS)

A new framework for single-signal sensing and compression has developed under the rubric of compressive sensing (CS). CS builds on the work of Candès, Romberg, and Tao and Donoho , who showed that if a signal has a sparse representation in one basis then it can be recovered from a small number of projections onto a second basis that is incoherent with the first.Roughly speaking, incoherence means that no element of one basis has a sparse representation in terms of the other basis. This notion has a variety of formalizations in the CS literature . CS relies on tractable recovery procedures that can provide exact recovery of a signal of length $N$ and sparsity $K$, i.e., a signal that can be written as a sum of $K$ basis functions from some known basis, where $K$ can be orders of magnitude less than $N$.

The implications of CS are promising for many applications, especially for sensing signals that have a sparse representation in some basis. Instead of sampling a $K$-sparse signal $N$ times, only $M=O(K\log N)$ incoherent measurements suffice, where $K$ can be orders of magnitude less than $N$. Moreover, the $M$ measurements need not be manipulated in any way before being transmitted, except possibly for some quantization. Finally, independent and identically distributed (i.i.d.) Gaussian or Bernoulli/Rademacher (random $\pm 1$) vectors provide a useful universal basis that is incoherent with all others. Hence, when using a random basis, CS is universal in the sense that the sensor can apply the same measurement mechanism no matter what basis sparsifies the signal .

While powerful, the CS theory at present is designed mainly to exploit intra-signal structures at a single sensor. In a multi-sensor setting, one can naively obtain separate measurements from each signal and recover them separately. However, it is possible to obtain measurements that each depend on all signals in the ensemble by having sensors collaborate with each other in order to combine all of their measurements; we term this process a joint measurement setting. In fact, initial work in CS for multi-sensor settings used standard CS with joint measurement and recovery schemes that exploit inter-signal correlations . However, by recovering sequential time instances of the sensed data individually, these schemes ignore intra-signal correlations.

3 Distributed compressive sensing (DCS)

In this paper we introduce a new theory for distributed compressive sensing (DCS) to enable new distributed coding algorithms that exploit both intra- and inter-signal correlation structures. In a typical DCS scenario, a number of sensors measure signals that are each individually sparse in some basis and also correlated from sensor to sensor. Each sensor separately encodes its signal by projecting it onto another, incoherent basis (such as a random one) and then transmits just a few of the resulting coefficients to a single collection point. Unlike the joint measurement setting described in Section 1.2, DCS requires no collaboration between the sensors during signal acquisition. Nevertheless, we are able to exploit the inter-signal correlation by using all of the obtained measurements to recover all the signals simultaneously. Under the right conditions, a decoder at the collection point can recover each of the signals precisely.

The DCS theory rests on a concept that we term the joint sparsity — the sparsity of the entire signal ensemble. The joint sparsity is often smaller than the aggregate over individual signal sparsities. Therefore, DCS offers a reduction in the number of measurements, in a manner analogous to the rate reduction offered by the Slepian-Wolf framework . Unlike the single-signal definition of sparsity, however, there are numerous plausible ways in which joint sparsity could be defined. In this paper, we first provide a general framework for joint sparsity using algebraic formulations based on a graphical model. Using this framework, we derive bounds for the number of measurements necessary for recovery under a given signal ensemble model. Similar to Slepian-Wolf coding , the number of measurements required for each sensor must account for the minimal features unique to that sensor, while at the same time features that appear among multiple sensors must be amortized over the group. Our bounds are dependent on the dimensionality of the subspaces in which each group of signals reside; they afford a reduction in the number of measurements that we quantify through the notions of joint and conditional sparsity, which are conceptually related to joint and conditional entropies. The common thread is that dimensionality and entropy both quantify the volume that the measurement and coding rates must cover. Our results are also applicable to cases where the signal ensembles are measured jointly, as well as to the single-signal case.

While our general framework does not by design provide insights for computationally efficient recovery, we also provide interesting models for joint sparsity where our results carry through from the general framework to realistic settings with low-complexity algorithms. In the first model, each signal is itself sparse, and so we could use CS to separately encode and decode each signal. However, there also exists a framework wherein a joint sparsity model for the ensemble uses fewer total coefficients. In the second model, all signals share the locations of the nonzero coefficients. In the third model, no signal is itself sparse, yet there still exists a joint sparsity among the signals that allows recovery from significantly fewer than $N$ measurements per sensor. For each model we propose tractable algorithms for joint signal recovery, followed by theoretical and empirical characterizations of the number of measurements per sensor required for accurate recovery. We show that, under these models, joint signal recovery can recover signal ensembles from significantly fewer measurements than would be required to recover each signal individually. In fact, for two of our three models we obtain best-possible performance that could not be bettered by an oracle that knew the the indices of the nonzero entries of the signals.

4 Paper organization

Section 2 overviews the single-signal CS theories and provides a new result on CS recovery. While some readers may be familiar with this material, we include it to make the paper self-contained. Section 3 introduces our general framework for joint sparsity models and proposes three example models for joint sparsity. We provide our detailed analysis for the general framework in Section 4; we then address the three models in Section 5. We close the paper with a discussion and conclusions in Section 6. Several appendices contain the proofs.

Compressive Sensing Background

The standard procedure for compressing sparse and nearly-sparse signals, known as transform coding, is to (i) acquire the full $N$-sample signal $x$; (ii) compute the complete set of transform coefficients $\{\vartheta(n)\}$; (iii) locate the $K$ largest, significant coefficients and discard the (many) small coefficients; (iv) encode the values and locations of the largest coefficients. This procedure has three inherent inefficiencies: First, for a high-dimensional signal, we must start with a large number of samples $N$. Second, the encoder must compute all of the $N$ transform coefficients $\{\vartheta(n)\}$, even though it will discard all but $K$ of them. Third, the encoder must encode the locations of the large coefficients, which requires increasing the coding rate since the locations change with each signal.

We will focus our theoretical development on exactly $K$-sparse signals and defer discussion of the more general situation of compressible signals where the coefficients decay rapidly with a power law but not to zero. Section 6 contains additional discussion on real-world compressible signals, and presents simulation results.

2 Incoherent projections

These inefficiencies raise a simple question: For a given signal, is it possible to directly estimate the set of large $\vartheta(n)$’s that will not be discarded? While this seems improbable, Candès, Romberg, and Tao and Donoho have shown that a reduced set of projections can contain enough information to recover sparse signals. A framework to acquire sparse signals, often referred to as compressive sensing (CS) , has emerged that builds on this principle.

In CS, we do not measure or encode the $K$ significant $\vartheta(n)$ directly. Rather, we measure and encode $M<N$ projections $y(m)=\langle x,\phi_{m}^{T}\rangle$ of the signal onto a second set of functions $\{\phi_{m}\},m=1,2,\ldots,M$, where $\phi_{m}^{T}$ denotes the transpose of $\phi_{m}$ and $\langle\cdot,\cdot\rangle$ denotes the inner product. In matrix notation, we measure $y=\Phi x$, where $y$ is an $M\times 1$ column vector and the measurement matrix $\Phi$ is $M\times N$ with each row a measurement vector $\phi_{m}$. Since $M<N$, recovery of the signal $x$ from the measurements $y$ is ill-posed in general; however the additional assumption of signal sparsity makes recovery possible and practical.

The CS theory tells us that when certain conditions hold, namely that the basis $\{\psi_{n}\}$ cannot sparsely represent the vectors $\{\phi_{m}\}$ (a condition known as incoherence ) and the number of measurements $M$ is large enough (proportional to $K$), then it is indeed possible to recover the set of large $\{\vartheta(n)\}$ (and thus the signal $x$) from the set of measurements $\{y(m)\}$ . This incoherence property holds for many pairs of bases, including for example, delta spikes and the sine waves of a Fourier basis, or the Fourier basis and wavelets. Signals that are sparsely represented in frames or unions of bases can be recovered from incoherent measurements in the same fashion. Significantly, this incoherence also holds with high probability between any arbitrary fixed basis or frame and a randomly generated one. In the sequel, we will focus our analysis to such random measurement procedures.

with overwhelming probability. (Thanks to the incoherence between the two bases, if the original signal is sparse in the $\vartheta$ coefficients, then no other set of sparse signal coefficients $\vartheta^{\prime}$ can yield the same projections $y$.)

Let $\Phi$ be an $M\times N$ measurement matrix, where $M\leq K$. Then, aside from pathological cases (specified in the proof), no signal $x=\Psi\vartheta$ with $\|\vartheta\|_{0}=K$ can be uniquely recovered from the $M$-dimensional measurement vector $y=\Phi x$.

The second statement of the theorem differs from the first in the following respect: when $K<M<2K$, there will necessarily exist $K$-sparse signals $x$ that cannot be uniquely recovered from the $M$-dimensional measurement vector $y=\Phi x$. However, these signals form a set of measure zero within the set of all $K$-sparse signals and can safely be avoided with high probability if $\Phi$ is randomly generated independently of $x$.

Comparing the second and third statements of Theorem 1, we see that one measurement separates the achievable region, where perfect recovery is possible with probability one, from the converse region, where with overwhelming probability recovery is impossible. Moreover, Theorem 1 provides a strong converse measurement region in a manner analogous to the strong channel coding converse theorems of information theory .

This optimization problem, also known as Basis Pursuit, is significantly more approachable and can be solved with traditional linear programming techniques whose computational complexities are polynomial in $N$.

There is no free lunch, however; according to the theory, more than $K+1$ measurements are required in order to recover sparse signals via Basis Pursuit. Instead, one typically requires $M\geq cK$ measurements, where $c>1$ is an overmeasuring factor. As an example, we quote a result asymptotic in $N$. For simplicity, we assume that the sparsity scales linearly with $N$; that is, $K=SN$, where we call $S$ the sparsity rate.

Set $K=SN$ with $0<S\ll 1$. Then there exists an overmeasuring factor $c(S)=O(\log(1/S))$, $c(S)>1$, such that, for a $K$-sparse signal $x$ in basis $\Psi$, the following statements hold:

In an illuminating series of papers, Donoho and Tanner have characterized the overmeasuring factor $c(S)$ precisely. In our work, we have noticed that the overmeasuring factor is quite similar to $\log_{2}(1+S^{-1})$. We find this expression a useful rule of thumb to approximate the precise overmeasuring ratio. Additional overmeasuring is proven to provide robustness to measurement noise and quantization error .

Throughout this paper we use the abbreviated notation $c$ to describe the overmeasuring factor required in various settings even though $c(S)$ depends on the sparsity $K$ and signal length $N$.

5 Signal recovery via greedy pursuit

Iterative greedy algorithms have also been developed to recover the signal $x$ from the measurements $y$. The Orthogonal Matching Pursuit (OMP) algorithm, for example, iteratively selects the vectors from the matrix $\Phi\Psi$ that contain most of the energy of the measurement vector $y$. The selection at each iteration is made based on inner products between the columns of $\Phi\Psi$ and a residual; the residual reflects the component of $y$ that is orthogonal to the previously selected columns. The algorithm has been proven to successfully recover the acquired signal from incoherent measurements with high probability, at the expense of slightly more measurements, . Algorithms inspired by OMP, such as regularized orthogonal matching pursuit , CoSaMP , and Subspace Pursuit have been shown to attain similar guarantees to those of their optimization-based counterparts. In the following, we will exploit both Basis Pursuit and greedy algorithms for recovering jointly sparse signals from incoherent measurements.

6 Properties of random measurements

In addition to offering substantially reduced measurement rates, CS has many attractive and intriguing properties, particularly when we employ random projections at the sensors. Random measurements are universal in the sense that any sparse basis can be used, allowing the same encoding strategy to be applied in different sensing environments. Random measurements are also future-proof: if a better sparsity-inducing basis is found for the signals, then the same measurements can be used to recover a more accurate view of the environment. Random coding is also robust: the measurements coming from each sensor have equal priority, unlike Fourier or wavelet coefficients in current coders. Finally, random measurements allow a progressively better recovery of the data as more measurements are obtained; one or more measurements can also be lost without corrupting the entire recovery.

7 Related work

Several researchers have formulated joint measurement settings for CS in sensor networks that exploit inter-signal correlations . In their approaches, each sensor $n\in\{1,2,\ldots,N\}$ simultaneously records a single reading $x(n)$ of some spatial field (temperature at a certain time, for example).Note that in Section 2.7 only, $N$ refers to the number of sensors, since each sensor acquires a signal sample. Each of the sensors generates a pseudorandom sequence $r_{n}(m),m=1,2,\ldots,M$, and modulates the reading as $x(n)r_{n}(m)$. Each sensor $n$ then transmits its $M$ numbers in sequence to the collection point where the measurements are aggregated, obtaining $M$ measurements $y(m)=\sum_{n=1}^{N}x(n)r_{n}(m)$. Thus, defining $x=[x(1),x(2),\ldots,x(N)]^{T}$ and $\phi_{m}=[r_{1}(m),r_{2}(m),\ldots,r_{N}(m)]$, the collection point automatically receives the measurement vector $y=[y(1),y(2),\ldots,y(M)]^{T}=\Phi x$ after $M$ transmission steps. The samples $x(n)$ of the spatial field can then be recovered using CS provided that $x$ has a sparse representation in a known basis. These methods have a major limitation: since they operate at a single time instant, they exploit only inter-signal and not intra-signal correlations; that is, they essentially assume that the sensor field is i.i.d. from time instant to time instant. In contrast, we will develop signal models and algorithms that are agnostic to the spatial sampling structure and that exploit both inter- and intra-signal correlations.

Recent work has adapted DCS to the finite rate of innovation signal acquisition framework and to the continuous-time setting . Since the original submission of this paper, additional work has focused on the analysis and proposal of recovery algorithms for jointly sparse signals .

Joint Sparsity Signal Models

In this section, we generalize the notion of a signal being sparse in some basis to the notion of an ensemble of signals being jointly sparse.

We denote by $\Phi_{j}$ the measurement matrix for signal $j$; $\Phi_{j}$ is $M_{j}\times N$ and, in general, the entries of $\Phi_{j}$ are different for each $j$. Thus, $y_{j}=\Phi_{j}x_{j}$ consists of $M_{j}<N$ random measurements of $x_{j}$. We will emphasize random i.i.d. Gaussian matrices $\Phi_{j}$ in the following, but other schemes are possible, including random $\pm 1$ Bernoulli/Rademacher matrices, and so on.

with 0 denoting a matrix of appropriate size with all entries equal to 0. We then have $Y=\Phi X$. Equation (4) shows that separate measurement matrices have a characteristic block-diagonal structure when the entries of the sparse vector are grouped by signal.

Below we propose a general framework for joint sparsity models (JSMs) and three example JSMs that apply in different situations.

2 General framework for joint sparsity

We now propose a general framework to quantify the sparsity of an ensemble of correlated signals $x_{1},x_{2},\ldots,x_{J}$, which allows us to compare the complexities of different signal ensembles and to quantify their measurement requirements. The framework is based on a factored representation of the signal ensemble that decouples its location and value information.

The joint sparsity level $D$ of the signal ensemble $X$ is the number of columns of the smallest matrix $P\in\mathcal{P}_{F}(X)$.

In contrast to the single-signal case, there are several natural choices for what matrices $P$ should be members of a joint sparsity model $\mathcal{P}$. We restrict our attention in the sequel to what we call common/innovation component JSMs. In these models each signal $x_{j}$ is generated as a combination of two components: ($i$) a common component $z_{C}$, which is present in all signals, and ($ii$) an innovation component $z_{j}$, which is unique to each signal. These combine additively, giving

Note, however, that the individual components might be zero-valued in specific scenarios. We can express the component signals as

The span of this location matrix (i.e., the set of signal ensembles $X$ that it can generate) remains unchanged if we remove any one of the columns, i.e., if we drop any entry of the value vector $\Theta$. This provides us with a lower-dimensional representation $\Theta^{\prime}$ of the same signal ensemble $X$ under the JSM $\mathcal{P}$; the joint sparsity of $X$ is $D=2$.

3 Example joint sparsity models

Since different real-world scenarios lead to different forms of correlation within an ensemble of sparse signals, we consider several possible designs for a JSM $\mathcal{P}$. The distinctions among our three JSMs concern the differing sparsity assumptions regarding the common and innovation components.

In this model, we suppose that each signal contains a common component $z_{C}$ that is sparse plus an innovation component $z_{j}$ that is also sparse. Thus, this joint sparsity model (JSM-1) $\mathcal{P}$ is represented by the set of all matrices of the form (5) with $K_{C}$ and all $K_{j}$ smaller than $N$. Assuming that sparsity reduction is not possible, the joint sparsity $D=K_{C}+\sum_{j\in\Lambda}K_{j}$.

A practical situation well-modeled by this framework is a group of sensors measuring temperatures at a number of outdoor locations throughout the day. The temperature readings $x_{j}$ have both temporal (intra-signal) and spatial (inter-signal) correlations. Global factors, such as the sun and prevailing winds, could have an effect $z_{C}$ that is both common to all sensors and structured enough to permit sparse representation. More local factors, such as shade, water, or animals, could contribute localized innovations $z_{j}$ that are also structured (and hence sparse). A similar scenario could be imagined for a network of sensors recording light intensities, air pressure, or other phenomena. All of these scenarios correspond to measuring properties of physical processes that change smoothly in time and in space and thus are highly correlated .

3.2 JSM-2: Common sparse supports

The JSM-2 model is immediately applicable to acoustic and RF sensor arrays, where each sensor acquires a replica of the same Fourier-sparse signal but with phase shifts and attenuations caused by signal propagation. In this case, it is critical to recover each one of the sensed signals. Another useful application for this framework is MIMO communication .

Similar signal models have been considered in the area of simultaneous sparse approximation . In this setting, a collection of sparse signals share the same expansion vectors from a redundant dictionary. The sparse approximation can be recovered via greedy algorithms such as Simultaneous Orthogonal Matching Pursuit (SOMP) or MMV Order Recursive Matching Pursuit (M-ORMP) . We use the SOMP algorithm in our setting (Section 5.2) to recover from incoherent measurements an ensemble of signals sharing a common sparse structure.

3.3 JSM-3: Nonsparse common component + sparse innovations

A practical situation well-modeled by this framework is where several sources are recorded by different sensors together with a background signal that is not sparse in any basis. Consider, for example, a verification system in a component production plant, where cameras acquire snapshots of each component to check for manufacturing defects. While each image could be extremely complicated, and hence nonsparse, the ensemble of images will be highly correlated, since each camera is observing the same device with minor (sparse) variations.

JSM-3 can also be applied in non-distributed scenarios. For example, it motivates the compression of data such as video, where the innovations or differences between video frames may be sparse, even though a single frame may not be very sparse. In this case, JSM-3 suggests that we encode each video frame separately using CS and then decode all frames of the video sequence jointly. This has the advantage of moving the bulk of the computational complexity to the video decoder. The PRISM system proposes a similar scheme based on Wyner-Ziv distributed encoding .

There are many possible joint sparsity models beyond those introduced above, as well as beyond the common and innovation component signal model. Further work will yield new JSMs suitable for other application scenarios; an example application consists of multiple cameras taking digital photos of a common scene from various angles . Extensions are discussed in Section 6.

Theoretical Bounds on Measurement Rates

In this section, we seek conditions on $\mathcal{M}=(M_{1},M_{2},\ldots,M_{J})$, the tuple of number of measurements from each sensor, such that we can guarantee perfect recovery of $X$ given $Y$. To this end, we provide a graphical model for the general framework provided in Section 3.2. This graphical model is fundamental in the derivation of the number of measurements needed for each sensor, as well as in the formulation of a combinatorial recovery procedure. Thus, we generalize Theorem 1 to the distributed setting to obtain fundamental limits on the number of measurements that enable recovery of sparse signal ensembles.

Based on the models presented in Section 3, recovering $X$ requires determining a value vector $\Theta$ and location matrix $P$ such that $X=P\Theta$. Two challenges immediately present themselves. First, a given measurement depends only on some of the components of $\Theta$, and the measurement budget should be adjusted between the sensors according to the information that can be gathered on the components of $\Theta$. For example, if a component $\Theta(d)$ does not affect any signal coefficient $x_{j}(\cdot)$ in sensor $j$, then the corresponding measurements $y_{j}$ provide no information about $\Theta(d)$. Second, the decoder must identify a location matrix $P\in\mathcal{P}_{F}(X)$ from the set $\mathcal{P}$ and the measurements $Y$.

We introduce a graphical representation that captures the dependencies between the measurements in $Y$ and the value vector $\Theta$, represented by $\Phi$ and $P$. Consider a feasible decomposition of $X$ into a full-rank matrix $P\in\mathcal{P}_{F}(X)$ and the corresponding $\Theta$; the matrix $P$ defines the sparsities of the common and innovation components $K_{C}$ and $K_{j}$, $1\leq j\leq J$, as well as the joint sparsity $D=K_{C}+\sum_{j=1}^{J}K_{j}$. Define the following sets of vertices: ($i$) the set of value vertices $V_{V}$ has elements with indices $d\in\{1,\ldots,D\}$ representing the entries of the value vector $\Theta(d)$, and ($ii$) the set of measurement vertices $V_{M}$ has elements with indices $(j,m)$ representing the measurements $y_{j}(m)$, with $j\in\Lambda$ and $m\in\{1,\ldots,M_{j}\}$. The cardinalities for these sets are $|V_{V}|=D$ and $|V_{M}|=\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{8.30453pt}{0.43057pt}}{M}$, respectively.

We now introduce a bipartite graph $G=(V_{V},V_{M},E)$, that represents the relationships between the entries of the value vector and the measurements (see for details). The set of edges $E$ is defined as follows:

For every $d\in\{1,2,\dots,K_{C}\}\subseteq V_{V}$ and $j\in\Lambda$ such that column $d$ of $P_{C}$ does not also appear as a column of $P_{j}$, we have an edge connecting $d$ to each vertex $(j,m)\in V_{M}$ for $1\leq m\leq M_{j}$.

For every $d\in\{K_{C}+1,K_{C}+2,\dots,D\}\subseteq V_{V}$, we consider the sensor $j$ associated with column $d$ of $P$, and we have an edge connecting $d$ to each vertex $(j,m)\in V_{M}$ for $1\leq m\leq M_{j}$.

In words, we say that $y_{j}(m)$, the $m^{th}$ measurement of sensor $j$, measures $\Theta(d)$ if the vertex $d\in V_{V}$ is linked to the vertex $(j,m)\in V_{M}$ in the graph $G$. An example graph for a distributed sensing setting is shown in Figure 1.

2 Quantifying redundancies

In order to obtain sharp bounds on the number of measurements needed, our analysis of the measurement process must account for redundancies between the locations of the nonzero coefficients in the common and innovation components. To that end, we consider the overlaps between common and innovation components in each signal. When we have $z_{c}(n)\neq 0$ and $z_{j}(n)\neq 0$ for a certain signal $j$ and some index $1\leq n\leq N$, we cannot recover the values of both coefficients from the measurements of this signal alone; therefore, we will need to recover $z_{c}(n)$ using measurements of other signals that do not feature the same overlap. We thus quantify the size of the overlap for all subsets of signals $\Gamma\subset\Lambda$ under a feasible representation given by $P$ and $\Theta$, as described in Section 3.2.

The overlap size for the set of signals $\Gamma\subset\Lambda$, denoted $K_{C}(\Gamma,P)$, is the number of indices in which there is overlap between the common and the innovation component supports at all signals $j\notin\Gamma$:

We also define $K_{C}(\Lambda,P)=K_{C}(P)$ and $K_{C}(\emptyset,P)=0$.

For $\Gamma\subset\Lambda$, $K_{C}(\Gamma,P)$ provides a penalty term due to the need for recovery of common component coefficients that are overlapped by innovations in all other signals $j\notin\Gamma$. Intuitively, for each entry counted in $K_{C}(\Gamma,P)$, some sensor in $\Gamma$ must take one measurement to account for that entry of the common component — it is impossible to recover such entries from measurements made by sensors outside of $\Gamma$. When all signals $j\in\Lambda$ are considered, it is clear that all of the common component coefficients must be recovered from the obtained measurements.

3 Measurement bounds

Converse and achievable bounds for the number of measurements necessary for DCS recovery are given below. Our bounds consider each subset of sensors $\Gamma\subseteq\Lambda$, since the cost of sensing the common component can be amortized across sensors: it may be possible to reduce the rate at one sensor $j_{1}\in\Gamma$ (up to a point), as long as other sensors in $\Gamma$ offset the rate reduction. We quantify the reduction possible through the following definition.

The conditional sparsity of the set of signals $\Gamma$ is the number of entries of the vector $\Theta$ that must be recovered by measurements $y_{j}$, $j\in\Gamma$:

The joint sparsity gives the number of degrees of freedom for the signals in $\Lambda$, while the conditional sparsity gives the number of degrees of freedom for signals in $\Gamma$ when the signals in $\Lambda\setminus\Gamma$ are available as side information. Note also that Definition 1 for joint sparsity can be extended to a subset of signals $\Gamma$ by considering the number of entries of $\Theta$ that affect these signals:

The bipartite graph introduced in Section 4.1 is the cornerstone of Theorems 3, 4, and 5, which consider whether a perfect matching can be found in the graph; see the proofs in Appendices B, D, and E, respectively, for detail.

(Achievable, known $P$) Assume that a signal ensemble $X$ is obtained from a common/innovation component JSM $\mathcal{P}$. Let $\mathcal{M}=(M_{1},M_{2},\ldots,M_{J})$ be a measurement tuple, let $\{\Phi_{j}\}_{j\in\Lambda}$ be random matrices having $M_{j}$ rows of i.i.d. Gaussian entries for each $j\in\Lambda$, and write $Y=\Phi X$. Suppose there exists a full rank location matrix $P\in\mathcal{P}_{F}(X)$ such that

for all $\Gamma\subseteq\Lambda$. Then with probability one over $\{\Phi_{j}\}_{j\in\Gamma}$, there exists a unique solution $\widehat{\Theta}$ to the system of equations $Y=\Phi P\widehat{\Theta}$; hence, the signal ensemble $X$ can be uniquely recovered as $X=P\widehat{\Theta}$.

(Achievable, unknown $P$) Assume that a signal ensemble $X$ and measurement matrices $\{\Phi_{j}\}_{j\in\Lambda}$ follow the assumptions of Theorem 3. Suppose there exists a full rank location matrix $P^{*}\in\mathcal{P}_{F}(X)$ such that

for all $\Gamma\subseteq\Lambda$. Then $X$ can be uniquely recovered from $Y$ with probability one over $\{\Phi_{j}\}_{j\in\Gamma}$.

(Converse) Assume that a signal ensemble $X$ and measurement matrices $\{\Phi_{j}\}_{j\in\Lambda}$ follow the assumptions of Theorem 3. Suppose there exists a full rank location matrix $P\in\mathcal{P}_{F}(X)$ such that

for some $\Gamma\subseteq\Lambda$. Then there exists a solution $\widehat{\Theta}$ such that $Y=\Phi P\widehat{\Theta}$ but $\widehat{X}:=P\widehat{\Theta}\neq X$.

4 Discussion

The bounds in Theorems 3–5 are dependent on the dimensionality of the subspaces in which the signals reside. The number of noiseless measurements required for ensemble recovery is determined by the dimensionality $\dim(\mathcal{S})$ of the subspace $\mathcal{S}$ in the relevant signal model, because dimensionality and sparsity play a volumetric role akin to the entropy $H$ used to characterize rates in source coding. Whereas in source coding each bit resolves between two options, and $2^{NH}$ typical inputs are described using $NH$ bits , in CS we have $M=\dim(\mathcal{S})+O(1)$. Similar to Slepian-Wolf coding , the number of measurements required for each sensor must account for the minimal features unique to that sensor, while at the same time features that appear among multiple sensors must be amortized over the group.

Practical Recovery Algorithms and Experiments

Although we have provided a unifying theoretical treatment for the three JSM models, the nuances warrant further study. In particular, while Theorem 4 highlights the basic tradeoffs that must be made in partitioning the measurement budget among sensors, the result does not by design provide insight into tractable algorithms for signal recovery. We believe there is additional insight to be gained by considering each model in turn, and while the presentation may be less unified, we attribute this to the fundamental diversity of problems that can arise under the umbrella of jointly sparse signal representations. In this section, we focus on tractable recovery algorithms for each model and, when possible, analyze the corresponding measurement requirements.

We first characterize the sparse common signal and innovations model JSM-1 from Section 3.3.1. For simplicity, we limit our description to $J=2$ signals, but describe extensions to multiple signals as needed.

Assume the measurement matrices $\{\Phi_{j}\}_{j\in\Lambda}$ contain i.i.d. Gaussian entries. Then the signal ensemble $X$ can be recovered with probability one if the following conditions hold:

Our joint recovery scheme provides a significant savings in measurements, because the common component can be measured as part of any of the $J$ signals.

1.2 Stochastic signal model for JSM-1

To give ourselves a firm footing for analysis, in the remainder of Section 5.1 we use a stochastic process for JSM-1 signal generation. This framework provides an information theoretic setting where we can scale the size of the problem and investigate which measurement rates enable recovery. We generate the common and innovation components as follows. For $n\in\{1,\ldots,N\}$ the decision whether $z_{C}(n)$ and $z_{j}(n)$ is zero or not is an i.i.d. Bernoulli process, where the probability of a nonzero value is given by parameters denoted $S_{C}$ and $S_{j}$, respectively. The values of the nonzero coefficients are then generated from an i.i.d. Gaussian distribution. The outcome of this process is that $z_{C}$ and $z_{j}$ have sparsities $K_{C}\sim\mbox{Binomial}(N,S_{C})$ and $K_{j}\sim\mbox{Binomial}(N,S_{j})$. The parameters $S_{j}$ and $S_{C}$ are sparsity rates controlling the random generation of each signal. Our model resembles the Gaussian spike process , which is a limiting case of a Gaussian mixture model.

Likelihood of sparsity reduction and overlap: This stochastic model can yield signal ensembles for which the corresponding generating matrices $P$ allow for sparsity reduction; specifically, there might be overlap between the supports of the common component $z_{C}$ and all the innovation components $z_{j}$, $j\in\Lambda$. For $J=2$, the probability that a given index is present in all supports is $S_{R}:=S_{C}S_{1}S_{2}$. Therefore, the distribution of the cardinality of this overlap is $K_{R}\sim\mbox{Binomial}(N,S_{R})$. We must account for the reduction obtained from the removal of the corresponding number of columns from the location matrix $P$ when the total number of measurements $M_{1}+M_{2}$ is considered. In the same way we can show that the distributions for the number of indices in the overlaps required by Corollary 1 are $K_{C}(\{1\})\sim\mbox{Binomial}(N,S_{C,\{1\}})$ and $K_{C}(\{2\})\sim\mbox{Binomial}(N,S_{C,\{2\}})$, where $S_{C,\{1\}}:=S_{C}(1-S_{1})S_{2}~{}~{}~{}\mbox{and}~{}~{}~{}S_{C,\{2\}}:=S_{C}S_{1}(1-S_{2}).$

Measurement rate region: To characterize DCS recovery performance, we introduce a measurement rate region. We define the measurement rate $R_{j}$ in an asymptotic manner as

Thus, we also set $S_{X_{j}}=S_{C}+S_{j}-S_{C}S_{j}$, $j\in\{1,2\}$. For a measurement rate pair $(R_{1},R_{2})$ and sources $X_{1}$ and $X_{2}$, we evaluate whether we can recover the signals with vanishing probability of error as $N$ increases. In this case, we say that the measurement rate pair is achievable.

Anticipated converse: As shown in Corollary 1, for indices $n$ such that $x_{1}(n)$ and $x_{2}(n)$ differ and are nonzero, each sensor must take measurements to account for one of the two coefficients. In the case where $S_{1}=S_{2}=S_{I}$, the joint sparsity rate is $S_{C}+2S_{I}-S_{C}S_{I}^{2}$. We define the measurement function $c^{\prime}(S):=S\cdot c(S)$ based on Donoho and Tanner’s oversampling factor $c(S)$ (Theorem 2). It can be shown that the function $c^{\prime}(\cdot)$ is concave; in order to minimize the sum rate bound, we “explain” as many of the sparse coefficients in one of the signals and as few as possible in the other. From Corollary 1, we have $R_{1},R_{2}\geq S_{I}+S_{C}S_{I}-S_{C}S_{I}^{2}$. Consequently, one of the signals must “explain” this sparsity rate, whereas the other signal must explain the rest:

Unfortunately, the derivation of $c^{\prime}(S)$ relies on Gaussianity of the measurement matrix, whereas in our case $\Phi$ has a block matrix form. Therefore, the following conjecture remains to be proved rigorously.

Let $J=2$ and fix the sparsity rate of the common component $S_{C}$ and the innovation sparsity rates $S_{1}=S_{2}=S_{I}$. Then the following conditions on the measurement rates are necessary to enable recovery with probability one:

Let $J=2$, $N\to\infty$ and fix the sparsity rate of the common component $S_{C}$ and the innovation sparsity rates $S_{1}=S_{2}=S_{I}$. If the measurement rates satisfy the following conditions:

The achievable measurement rate region of Theorem 6 is loose with respect to the region of the anticipated converse Conjecture 1 (see Figure 2). We leave for future work the characterization of a tight measurement rate region for computationally tractable (polynomial time) recovery techniques.

1.6 Simulations for JSM-1

Recovering two signals with symmetric measurement rates: Our simulation setting is as follows. The signal components $z_{C}$, $z_{1}$, and $z_{2}$ are assumed (without loss of generality) to be sparse in $\Psi=I_{N}$ with sparsities $K_{C}$, $K_{1}$, and $K_{2}$, respectively. We assign random Gaussian values to the nonzero coefficients. We restrict our attention to the symmetric setting in which $K_{1}=K_{2}$ and $M_{1}=M_{2}$, and consider signals of length $N=50$ where $K_{C}+K_{1}+K_{2}=15$.

Recovering two signals with asymmetric measurement rates: In Figure 2, we compare separate CS recovery with the anticipated converse bound of Conjecture 1, the achievable bound of Theorem 6, and numerical results.

We use $J\in\{1,2,\ldots,10\}$ signals and choose the same sparsity rates $S_{C}=0.2$ and $S_{I}=0.05$ as the asymmetric rate simulations; here we use symmetric measurement rates and let $N=500$. The results of Figure 4 describe the smallest symmetric measurement rates for which we always succeeded recovering the signal over $100$ simulation runs. As $J$ increases, lower measurement rates can be used; the results compare favorably with the lower bound from Conjecture 1, which gives $R_{j}\approx 0.232$ as $J\to\infty$.

2 Recovery strategies for common sparse supports (JSM-2)

The algorithms we propose are inspired by conventional greedy pursuit algorithms for CS (such as OMP ). In the single-signal case, OMP iteratively constructs the sparse support set $\Omega$; decisions are based on inner products between the columns of $\Phi$ and a residual. In the multi-signal case, there are more clues available for determining the elements of $\Omega$.

When there are many correlated signals in the ensemble, a simple non-iterative greedy algorithm based on inner products will suffice to recover the signals jointly. For simplicity but without loss of generality, we assume that an equal number of measurements $M_{j}=M$ are taken of each signal. We write $\Phi_{j}$ in terms of its columns, with $\Phi_{j}=\left[\phi_{j,1},\phi_{j,2},\dots,\phi_{j,N}\right]$.

Get greedy: Given all of the measurements, compute the test statistics

and estimate the elements of the common coefficient support set by

When the sparse, nonzero coefficients are sufficiently generic (as defined below), we have the following surprising result, which is proved in Appendix H.

In words, with fewer than $K$ measurements per sensor, it is actually possible to recover the sparse support set $\Omega$ under the JSM-2 model. One can also show the somewhat stronger result that, as long as $\sum_{j}M_{j}\gg N$, TP recovers $\Omega$ with probability approaching one. We have omitted this additional result for brevity. Of course, this approach does not recover the $K$ coefficient values for each signal; at least $K$ measurements per sensor are required for this.

Assume that the nonzero coefficients in the $\theta_{j}$ are i.i.d. Gaussian random variables. Then the following statements hold:

Let the measurement matrices $\Phi_{j}$ contain i.i.d. Gaussian entries, with each matrix having an overmeasuring factor of $c=1$ (that is, $M_{j}=K$ for each measurement matrix $\Phi_{j}$). Then TP recovers all signals from the ensemble $\{x_{j}\}$ with probability approaching one as $J\rightarrow\infty$.

Let $\Phi_{j}$ be a measurement matrix with overmeasuring factor $c<1$ (that is, $M_{j}<K$), for some $j\in\Lambda$. Then with probability one, the signal $x_{j}$ cannot be uniquely recovered by any algorithm for any value of $J$.

The first statement is an immediate corollary of Theorem 7; the second statement follows because each equation $y_{j}=\Phi_{j}x_{j}$ would be underdetermined even if the nonzero indices were known. Thus, under the JSM-2 model, the TP algorithm asymptotically performs as well as an oracle decoder that has prior knowledge of the locations of the sparse coefficients. From an information theoretic perspective, Corollary 2 provides tight achievable and converse bounds for JSM-2 signals. We should note that the theorems in this section have a slightly different flavor than Theorem 4 and 5, which ensure recovery of any sparse signal ensemble, given a suitable set of measurement matrices. Theorem 7 and Corollary 2 above, in contrast, rely on a random signal model and do not guarantee simultaneous performance for all sparse signals under any particular measurement ensemble. Nonetheless, we feel this result is worth presenting to highlight the strong subspace concentration behavior that enables the correct identification of the common support.

In the technical reports , we derive an approximate formula for the probability of error in recovering the common support set $\Omega$ given $J$, $K$, $M$, and $N$. While theoretically interesting and potentially practically useful, these results require $J$ to be large. Our numerical experiments show that the number of measurements required for recovery using TP decreases quickly as $J$ increases. However, in the case of small $J$, TP performs poorly. Hence, we propose next an alternative recovery technique based on simultaneous greedy pursuit that performs well for small $J$.

2.2 Recovery via iterative greedy pursuit

In practice, the common sparse support among the $J$ signals enables a fast iterative algorithm to recover all of the signals jointly. Tropp and Gilbert have proposed one such algorithm, called Simultaneous Orthogonal Matching Pursuit (SOMP) , which can be readily applied in our DCS framework. SOMP is a variant of OMP that seeks to identify $\Omega$ one element at a time. A similar simultaneous sparse approximation algorithm has been proposed using convex optimization . We dub the DCS-tailored SOMP algorithm DCS-SOMP.

To adapt the original SOMP algorithm to our setting, we first extend it to cover a different measurement matrix $\Phi_{j}$ for each signal $x_{j}$. Then, in each DCS-SOMP iteration, we select the column index $n\in\{1,2,\dots,N\}$ that accounts for the greatest amount of residual energy across all signals. As in SOMP, we orthogonalize the remaining columns (in each measurement matrix) after each step; after convergence we obtain an expansion of the measurement vector $y_{j}$ on an orthogonalized subset of the columns of basis vectors. To obtain the expansion coefficients in the sparse basis, we then reverse the orthogonalization process using the QR matrix factorization. Finally, we again assume that $M_{j}=M$ measurements per signal are taken.

Select the dictionary vector that maximizes the value of the sum of the magnitudes of the projections of the residual, and add its index to the set of selected indices

Orthogonalize the selected basis vector against the orthogonalized set of previously selected dictionary vectors

Iterate: Update the estimate of the coefficients for the selected vector and residuals

De-orthogonalize: Consider the relationship between $\Gamma_{j}=[\gamma_{j,1},\gamma_{j,2},\ldots,\gamma_{j,M}]$ and the $\Phi_{j}$ given by the QR factorization $\Phi_{j,{\widehat{\Omega}}}=\Gamma_{j}R_{j}$, where $\Phi_{j,{\widehat{\Omega}}}=[\phi_{j,n_{1}},\phi_{j,n_{2}},\ldots,\phi_{j,n_{M}}]$ is the so-called mutilated basis.We define a mutilated basis $\Phi_{\Omega}$ as a subset of the basis vectors from $\Phi=[\phi_{1},\phi_{2},\ldots,\phi_{N}]$ corresponding to the indices given by the set $\Omega=\{n_{1},n_{2},\ldots,n_{M}\}$, that is, $\Phi_{\Omega}=[\phi_{n_{1}},\phi_{n_{2}},\ldots,\phi_{n_{M}}]$. This concept can be extended to vectors in the same manner. Since $y_{j}=\Gamma_{j}\beta_{j}=\Phi_{j,{\widehat{\Omega}}}x_{j,{\widehat{\Omega}}}=\Gamma_{j}R_{j}x_{j,{\widehat{\Omega}}}$, where $x_{j,{\widehat{\Omega}}}$ is the mutilated coefficient vector, we can compute the signal estimates $\{\widehat{x}_{j}\}$ as

where $\widehat{x}_{j,{\widehat{\Omega}}}$ is the mutilated version of the sparse coefficient vector $\widehat{x}_{j}$.

In practice, we obtain $\widehat{c}K$ measurements from each signal $x_{j}$ for some value of $\widehat{c}$. We then use DCS-SOMP to recover the $J$ signals jointly. We orthogonalize because as the number of iterations approaches $M$ the norms of the residues of an orthogonal pursuit decrease faster than for a non-orthogonal pursuit; indeed, due to Step 3 the algorithm can only run for up to $M$ iterations. The computational complexity of this algorithm is $O(JNM^{2})$, which matches that of separate recovery for each signal while reducing the required number of measurements.

Thanks to the common sparsity structure among the signals, we believe (but have not proved) that DCS-SOMP will succeed with $\widehat{c}<c(S)$. Empirically, we have observed that a small number of measurements proportional to $K$ suffices for a moderate number of sensors $J$. Based on our observations, described in Section 5.2.3, we conjecture that $K+1$ measurements per sensor suffice as $J\rightarrow\infty$. Thus, this efficient greedy algorithm enables an overmeasuring factor $\widehat{c}=(K+1)/K$ that approaches $1$ as $J$, $K$, and $N$ increase.

2.3 Simulations for JSM-2

We now present simulations comparing separate CS recovery versus joint DCS-SOMP recovery for a JSM-2 signal ensemble. Figure 5 plots the probability of perfect recovery corresponding to various numbers of measurements $M$ as the number of sensors varies from $J=1$ to $32$, over 1000 trials in each case. We fix the signal lengths at $N=50$ and the sparsity of each signal to $K=5$.

With DCS-SOMP, for perfect recovery of all signals the average number of measurements per signal decreases as a function of $J$. The trend suggests that for large $J$ close to $K$ measurements per signal should suffice. On the contrary, with separate CS recovery, for perfect recovery of all signals the number of measurements per sensor increases as a function of $J$. This occurs because each signal experiences an independent probability $p\leq 1$ of successful recovery; therefore the overall probability of complete success is $p^{J}$. Consequently, each sensor must compensate by making additional measurements. This phenomenon further motivates joint recovery under JSM-2.

Finally, we note that we can use algorithms other than DCS-SOMP to recover the signals under the JSM-2 model. Cotter et al. have proposed additional algorithms (such as M-FOCUSS) that iteratively eliminate basis vectors from the dictionary and converge to the set of sparse basis vectors over which the signals are supported. We hope to extend such algorithms to JSM-2 in future work.

3 Recovery strategies for nonsparse common component + sparse innovations (JSM-3)

The JSM-3 signal ensemble model from Section 3.3.3 provides a particularly compelling motivation for joint recovery. Under this model, no individual signal $x_{j}$ is sparse, and so recovery of each signal separately would require fully $N$ measurements per signal. As in the other JSMs, however, the commonality among the signals makes it possible to substantially reduce this number. Again, the potential for this savings is evidenced by specializing Theorem 4 to the context of JSM-3.

If $\Phi_{j}$ is a random Gaussian matrix for all $j\in\Lambda$, $\mathcal{P}$ is defined by JSM-3, and

then the signal ensemble $X$ can be uniquely recovered from $Y$ with probability one.

This suggests that the number of measurements of an individual signal can be substantially decreased, as long as the total number of measurements is sufficiently large to capture enough information about the nonsparse common component $z_{C}$. The term $K_{R}$ denotes the number of indices where the common and all innovation components overlap, and appears due to the sparsity reduction that can be performed at the common component before recovery. We also note that when the supports of the innovations are independent, as $J\to\infty$, it becomes increasingly unlikely that a given index will be included in all innovations, and thus the terms $K_{C}(\Gamma,P)$ and $K_{R}$ will go to zero. On the other hand, when the supports are completely matched (implying $K_{j}=K$, $j\in\Lambda$), we will have $K_{R}=K$, and after sparsity reduction has been addressed, $K_{C}(\Gamma,P)=0$ for all $\Gamma\subseteq\Lambda$.

Successful recovery of the signal ensemble $\{x_{j}\}$ requires recovery of both the nonsparse common component $z_{C}$ and the sparse innovations $\{z_{j}\}$. To help build intuition about how we might accomplish signal recovery using far fewer than $N$ measurements per sensor, consider the following thought experiment.

Estimate common component: Define the matrix $\widehat{\Phi}$ as the vertical concatenation of the regularized individual measurement matrices $\widehat{\Phi}_{j}=\frac{1}{M_{j}\sigma_{j}^{2}}\Phi_{j}$, that is, $\widehat{\Phi}=[\widehat{\Phi}_{1}^{T},\widehat{\Phi}_{2}^{T},\ldots,\widehat{\Phi}_{J}^{T}]^{T}$. Calculate the estimate of the common component as $\widehat{z}_{C}=\frac{1}{J}\widehat{\Phi}^{T}Y$.

Estimate measurements generated by innovations: Using the previous estimate, subtract the contribution of the common part from the measurements and generate estimates for the measurements caused by the innovations for each signal: $\widehat{y}_{j}=y_{j}-\Phi_{j}\widehat{z}_{C}$.

Obtain signal estimates: Estimate each signal as the sum of the common and innovations estimates; that is, $\widehat{x}_{j}=\widehat{z}_{C}+\widehat{z}_{j}$.

The following theorem, proved in Appendix I, shows that asymptotically, by using the TECC algorithm, each sensor needs to only measure at the rate dictated by the sparsity $K_{j}$.

Assume that the nonzero expansion coefficients of the sparse innovations $z_{j}$ are i.i.d. Gaussian random variables and that their locations are uniformly distributed on $\{1,2,\ldots,N\}$. Let the measurement matrices $\Phi_{j}$ contain i.i.d. $\mathcal{N}(0,\sigma_{j}^{2})$ entries with $M_{j}\geq K_{j}+1$. Then each signal $x_{j}$ can be recovered using the TECC algorithm with probability approaching one as $J\rightarrow\infty$.

For large $J$, the measurement rates permitted by Theorem 8 are the best possible for any recovery strategy for JSM-3 signals, even neglecting the presence of the nonsparse component. These rates meet the minimum bounds suggested by Corollary 3, although again Theorem 8 is of a slightly different flavor, as it does not provide a uniform guarantee for all sparse signal ensembles under any particular measurement matrix collection. The CS technique employed in Theorem 8 involves combinatorial searches that estimate the innovation components; we have provided the theorem simply as support for our intuitive development of the TECC algorithm. More efficient techniques could also be employed (including several proposed for CS in the presence of noise ). It is reasonable to expect similar behavior; as the error in estimating the common component diminishes, these techniques should perform similarly to their noiseless analogues.

3.2 Recovery via Alternating Common and Innovation Estimation (ACIE)

The preceding analysis demonstrates that the number of required measurements in JSM-3 can be substantially reduced through joint recovery. While Theorem 8 shows theoretical gains as $J\rightarrow\infty$, practical gains can also be realized with a moderate number of sensors. In particular, suppose in the TECC algorithm that the initial estimate $\widehat{z}_{C}$ is not accurate enough to enable correct identification of the sparse innovation supports $\{\Omega_{j}\}$. In such a case, it may still be possible for a rough approximation of the innovations $\{z_{j}\}$ to help refine the estimate $\widehat{z}_{C}$. This in turn could help to refine the estimates of the innovations.

where $\Phi_{j,{\widehat{\Omega}}_{j}}$ is the mutilated matrix $\Phi_{j}$ corresponding to the indices in ${\widehat{\Omega}}_{j}$, and the $M_{j}\times(M_{j}-|{\widehat{\Omega}}_{j}|)$ matrix $Q_{j}$ has orthonormal columns that span the orthogonal complement of $\Phi_{j,{\widehat{\Omega}}_{j}}$.

This construction allows us to remove the projection of the measurements into the aforementioned span to obtain measurements caused exclusively by vectors not in ${\widehat{\Omega}}_{j}$:

Thus, the modified measurements $\widetilde{Y}=\left[\widetilde{y}_{1}^{T}~{}\widetilde{y}_{2}^{T}~{}\ldots~{}\widetilde{y}_{J}^{T}\right]^{T}$ and modified measurement matrix $\widetilde{\Phi}=\left[\widetilde{\Phi}_{1}^{T}~{}\widetilde{\Phi}_{2}^{T}~{}\ldots~{}\widetilde{\Phi}_{J}^{T}\right]^{T}$ can be used to refine the estimate of the common component of the signal,

where $A^{\dagger}=(A^{T}A)^{-1}A^{T}$ denotes the pseudoinverse of matrix $A$.

When the innovation support estimate is correct (${\widehat{\Omega}}_{j}=\Omega_{j}$), the measurements $\widetilde{y_{j}}$ will describe only the common component $z_{C}$. If this is true for every signal $j$ and the number of remaining measurements $\sum_{j=1}^{J}(M_{j}-K_{j})\geq N$, then $z_{C}$ can be perfectly recovered via (18). However, it may be difficult to obtain correct estimates for all signal supports in the first iteration of the algorithm, and so we find it preferable to refine the estimate of the support by executing several iterations.

Estimate common component: Update estimate $\widetilde{z_{C}}$ according to (17)–(18).

Estimate innovation supports: For each sensor $j$, after subtracting the contribution $\widetilde{z_{C}}$ from the measurements, $\widehat{y_{j}}=y_{j}-\Phi_{j}\widetilde{z_{C}}$, estimate the support of each signal innovation ${\widehat{\Omega}}_{j}$.

Estimate innovation coefficients: For each $j$, estimate the coefficients for the indices in ${\widehat{\Omega}}_{j}$,

where $\widehat{z}_{j,{\widehat{\Omega}}_{j}}$ is a mutilated version of the innovation’s sparse coefficient vector estimate $\widehat{z}_{j}$.

Recover signals: Compute the estimate of each signal as $\widehat{x}_{j}=\widetilde{z_{C}}+\widehat{z}_{j}$.

Estimation of the supports in Step 3 can be accomplished using a variety of techniques. We propose to run a fixed number of iterations of OMP; if the supports of the innovations are known to match across signals — as in JSM-2 — then more powerful algorithms like SOMP can be used. The ACIE algorithm is similar in spirit to other iterative estimation algorithms, such as turbo decoding .

3.3 Simulations for JSM-3

We now present simulations of JSM-3 recovery for the following scenario. Consider $J$ signals of length $N=50$ containing a common white noise component $z_{C}(n)\sim\mathcal{N}(0,1)$ for $n\in\{1,\ldots,N\}$. Each innovations component $z_{j}$ has sparsity $K=5$ (once again in the time domain), resulting in $x_{j}=z_{C}+z_{j}$. The signals are generated according to the model used in Section 5.1.6.

We study two different cases. The first is an extension of JSM-1: we select the supports for the various innovations separately and then apply OMP to each signal in Step 3 of the ACIE algorithm in order to estimate its innovations component. The second case is an extension of JSM-2: we select one common support for all of the innovations across the signals and then apply the DCS-SOMP algorithm (Section 5.2.2) to estimate the innovations in Step 3. In both cases we use $L=10$ iterations of ACIE. We test the algorithms for different numbers of signals $J$ and calculate the probability of correct recovery as a function of the (same) number of measurements per signal $M$.

Figure 6(a) shows that, for sufficiently large $J$, we can recover all of the signals with significantly fewer than $N$ measurements per signal. As $J$ grows, it becomes more difficult to perfectly recover all $J$ signals. We believe this is inevitable, because even if $z_{C}$ were known without error, then perfect ensemble recovery would require the successful execution of $J$ independent runs of OMP. Second, for small $J$, the probability of success can decrease at high values of $M$. We believe this behavior is due to the fact that initial errors in estimating $z_{C}$ may tend to be somewhat sparse (since $\widehat{z}_{C}$ roughly becomes an average of the signals $\{x_{j}\}$), and these sparse errors can mislead the subsequent OMP processes. For more moderate $M$, it seems that the errors in estimating $z_{C}$ (though greater) tend to be less sparse. We expect that a more sophisticated algorithm could alleviate such a problem; the problem is also mitigated at higher $J$.

Figure 6(b) shows that when the sparse innovations share common supports we see an even greater savings. As a point of reference, a traditional approach to signal acquisition would require $1600$ total measurements to recover these $J=32$ nonsparse signals of length $N=50$. Our approach requires only approximately $10$ random measurements per sensor — a total of $320$ measurements — for high probability of recovery.

Discussion and Conclusions

In this paper we have extended the theory and practice of compressive sensing to multi-signal, distributed settings. The number of noiseless measurements required for ensemble recovery is determined by the dimensionality of the subspace in the relevant signal model, because dimensionality and sparsity play a volumetric role akin to the entropy used to characterize rates in source coding. Our three example joint sparsity models (JSMs) for signal ensembles with both intra- and inter-signal correlations capture the essence of real physical scenarios, illustrate the basic analysis and algorithmic techniques, and indicate the significant gains to be realized from joint recovery. In some sense, distributed compressive sensing (DCS) is a framework for distributed compression of sources with memory, which has remained a challenging problem for some time.

In addition to offering substantially reduced measurement rates, the DCS-based distributed source coding schemes we develop here share the properties of CS mentioned in Section 2. Two additional properties of DCS make it well-matched to distributed applications such as sensor networks and arrays . First, each sensor encodes its measurements separately, which eliminates the need for inter-sensor communication. Second, DCS distributes its computational complexity asymmetrically, placing most of it in the joint decoder, which will often have more computational resources than any individual sensor node. The encoders are very simple; they merely compute incoherent projections with their signals and make no decisions.

Second, (random) measurements are real numbers; quantization gradually degrades the recovery quality as the quantization becomes coarser . Moreover, in many practical situations some amount of measurement noise will corrupt the $\{x_{j}\}$, making them not exactly sparse in any basis. While characterizing these effects and the resulting rate-distortion consequences in the DCS setting are topics for future work, there has been work in the single-signal CS literature that we should be able to leverage, including variants of Basis Pursuit with Denoising , robust iterative recovery algorithms , CS noise sensitivity analysis , the Dantzig Selector , and one-bit CS .

Third, in some applications, the linear program associated with some DCS decoders (in JSM-1 and JSM-3) could prove too computationally intense. As we saw in JSM-2, efficient iterative and greedy algorithms could come to the rescue, but these need to be extended to the multi-signal case. Recent results on recovery from a union of subspaces give promise for efficient, model-based algorithms .

Finally, we focused our theory on models that assign common and innovation components to the signals in the ensemble. Other models tailored to specific applications can be posed; for example, in hyperspectral imaging applications, it is common to obtain strong correlations only across spectral slices within a certain neighborhood. It would be then appropriate to pose a common/innovation model with separate common components that are localized to a subset of the spectral slices obtained. Results similar to those obtained in Section 4 are simple to derive for models with full-rank location matrices.

Appendix A Proof of Theorem 1

Statement 2 is an application of the achievable bound of Theorem 4 to the case of $J=1$ signal. It remains then to prove Statements 1 and 3.

Statement 3 (Converse, $M\leq K$): If $M<K$, there is insufficient information in the vector $y$ to recover the $K$ nonzero coefficients of $\theta$; thus we assume $M=K$. In this case, there is a single explanation for the measurements only if there is a single set $\Omega$ of $K$ linearly independent columns and the nonzero indices of $\theta$ are the elements of $\Omega$. Aside from this pathological case, the rank of subsets $\Phi_{\widehat{\Omega}}$ will generally be less than $K$ — which would prevent robust recovery of signals supported on ${\widehat{\Omega}}$, or will be equal to $K$ — which would give ambiguous solutions among all such sets ${\widehat{\Omega}}$. $\Box$

Appendix B Proof of Theorem 3

has rank $D$, and thus $\Theta$ is the unique solution to the equation $Y=\Phi P\Theta=\Upsilon\Theta$.

We recall that, under our common/innovation model, $P$ has the form (5), where $P_{C}$ is an $N\times K_{C}$ submatrix of the $N\times N$ identity, and each $P_{j}$, $j\in\Lambda$, is an $N\times K_{j}$ submatrix of the $N\times N$ identity. To prove that $\Upsilon$ has rank $D$, we will require the following lemma, which we prove in Appendix C.

If (7) holds, then there exists a mapping $\mathcal{C}:\{1,2,\dots,K_{C}\}\rightarrow\Lambda$, assigning each element of the common component to one of the sensors, such that for each $\Gamma\subseteq\Lambda$,

Intuitively, the existence of such a mapping suggests that ($i$) each sensor has taken enough measurements to cover its own innovation (requiring $K_{j}$ measurements) and perhaps some of the common component, ($ii$) for any $\Gamma\subseteq\Lambda$, the sensors in $\Gamma$ have collectively taken enough extra measurements to cover the requisite $K_{C}(\Gamma,P)$ elements of the common component, and ($iii$) the extra measurements are taken at sensors where the common and innovation components do not overlap. Formally, we will use the existence of such a mapping to prove that $\Upsilon$ has rank $D$.

We proceed by noting that $\Upsilon$ has the form

where each $\Phi_{j}P_{C}$ (respectively, $\Phi_{j}P_{j}$) is an $M_{j}\times K_{C}$ (respectively, $M_{j}\times K_{j}$) submatrix of $\Phi_{j}$ obtained by selecting columns from $\Phi_{j}$ according to the nonzero entries of $P_{C}$ (respectively, $P_{j}$). In total, $\Upsilon$ has $D$ columns (19). To argue that $\Upsilon$ has rank $D$, we will consider a sequence of three matrices $\Upsilon_{0}$, $\Upsilon_{1}$, and $\Upsilon_{2}$ constructed from small modifications to $\Upsilon$.

We begin by letting $\Upsilon_{0}$ denote the “partially zeroed” matrix obtained from $\Upsilon$ using the following construction. We first let $\Upsilon_{0}=\Upsilon$ and then make the following adjustments:

For each $j$ such that $P_{j}$ has a column that matches column $k$ of $P_{C}$ (note that by Lemma 2 this cannot happen if $\mathcal{C}(k)=j$), let $k^{\prime}$ represent the column index of the full matrix $P$ where this column of $P_{j}$ occurs. Subtract column $k^{\prime}$ of $\Upsilon_{0}$ from column $k$ of $\Upsilon_{0}$. This forces to zero all entries of $\Upsilon_{0}$ formerly corresponding to column $k$ of the block $\Phi_{j}P_{C}$.

If $k<K_{C}$, add one to $k$ and go to step 2.

The matrix $\Upsilon_{0}$ is identical to $\Upsilon$ everywhere except on the first $K_{C}$ columns, where any portion of a column overlapping with a column of $\Phi_{j}P_{j}$ to its right has been set to zero. Thus, $\Upsilon_{0}$ satisfies the next two properties, which will be inherited by matrices $\Upsilon_{1}$ and $\Upsilon_{2}$ that we subsequently define:

Each entry of $\Upsilon_{0}$ is either zero or a Gaussian random variable.

All Gaussian random variables in $\Upsilon_{0}$ are i.i.d.

Finally, because $\Upsilon_{0}$ was constructed only by subtracting columns of $\Upsilon$ from one another,

We now let $\Upsilon_{1}$ be the matrix obtained from $\Upsilon_{0}$ using the following construction. For each $j\in\Lambda$, we select $K_{j}+\sum_{k=1}^{K_{C}}1_{\mathcal{C}(k)=j}$ arbitrary rows from the portion of $\Upsilon_{0}$ corresponding to sensor $j$. Using (19), the resulting matrix $\Upsilon_{1}$ has

rows. Also, because $\Upsilon_{1}$ was obtained by selecting a subset of rows from $\Upsilon_{0}$, it has $D$ columns and satisfies

We now let $\Upsilon_{2}$ be the $D\times D$ matrix obtained by permuting columns of $\Upsilon_{1}$ using the following construction:

Let $\Upsilon_{1}^{\prime}$ denote the columns of $\Upsilon_{1}$ corresponding to the entries of $\Phi_{j}P_{j}$ (the innovation components of sensor $j$), and concatenate $\Upsilon_{1}^{\prime}$ to $\Upsilon_{2}$, i.e., let $\Upsilon_{2}\leftarrow[\Upsilon_{2}~{}\Upsilon_{1}^{\prime}]$. There are $K_{j}$ such columns.

If $j<J$, let $j\leftarrow j+1$ and go to Step 2.

Because $\Upsilon_{1}$ and $\Upsilon_{2}$ share the same columns up to reordering, it follows that

Based on its dependency on $\Upsilon_{0}$, and following from Lemma 2, the square matrix $\Upsilon_{2}$ meets properties P1 and P2 defined above in addition to a third property:

All diagonal entries of $\Upsilon_{2}$ are Gaussian random variables.

Having identified these three properties satisfied by $\Upsilon_{2}$, we will prove by induction that, with probability one over $\Phi$, such a matrix has full rank.

Let $\Upsilon^{(d-1)}$ be a $(d-1)\times(d-1)$ matrix having full rank. Construct a $d\times d$ matrix $\Upsilon^{(d)}$ as follows:

Proof of Lemma 3: When $d=1$, $\Upsilon^{(d)}=[\omega]$, which has full rank if and only if $\omega\neq 0$, which occurs with probability one.

When $d>1$, using expansion by minors, the determinant of $\Upsilon^{(d)}$ satisfies

where $C=C(\Upsilon^{(d-1)},v_{1},v_{2})$ is independent of $\omega$. The matrix $\Upsilon^{(d)}$ has full rank if and only if $\det(\Upsilon^{(d)})\neq 0$, which is satisfied if and only if

By assumption, $\det(\Upsilon^{(d-1)})\neq 0$ and $\omega$ is a Gaussian random variable that is independent of $C$ and $\det(\Upsilon^{(d-1)})$. Thus, $\omega\neq\frac{-C}{\det(\Upsilon^{(d-1)})}$ with probability one. $\Box$

Appendix C Proof of Lemma 2

To prove this lemma, we apply tools from graph theory.

We seek a matching within the graph $G=(V_{V},V_{M},E)$ from Figure 1, i.e., a subgraph $(V_{V},V_{M},\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{5.61287pt}{0.43057pt}}{E})$ with $\makebox[0.0pt][l]{\hskip 2.05pt\rule[8.12498pt]{5.61287pt}{0.43057pt}}{E}\subseteq E$ that pairs each element of $V_{V}$ with a unique element of $V_{M}$. Such a matching will immediately give us the desired mapping $\mathcal{C}$ as follows: for each $k\in\{1,2,\dots,K_{C}\}\subseteq V_{V}$, we let $(j,m)\in V_{M}$ denote the single node matched to $k$ by an edge in $E$, and we set $\mathcal{C}(k)=j$.

To prove the existence of such a matching within the graph, we invoke a version of Hall’s marriage theorem for bipartite graphs . Hall’s theorem states that within a bipartite graph $(V_{1},V_{2},E)$, there exists a matching that assigns each element of $V_{1}$ to a unique element of $V_{2}$ if for any collection of elements $\Pi\subseteq V_{1}$, the set $E(\Pi)$ of neighbors of $\Pi$ in $V_{2}$ has cardinality $|E(\Pi)|\geq|\Pi|$.

In the context of our lemma, Hall’s condition requires that for any set of entries in the value vector, $\Pi\subseteq V_{V}$, the set $E(\Pi)$ of neighbors of $\Pi$ in $V_{M}$ has size $|E(\Pi)|\geq|\Pi|$. We will prove that if (7) is satisfied, then Hall’s condition is satisfied, and thus a matching must exist.

We would now like to show that $\sum_{j\in S_{\Pi}}K_{j}+K_{C}(S_{\Pi},P)\geq|\Pi|$, and thus if (7) is satisfied for all $\Gamma\subseteq\Lambda$, then (24) is satisfied in particular for $S_{\Pi}\subseteq\Lambda$.

In general, the set $\Pi$ may contain vertices for both common components and innovation components. We write $\Pi=\Pi_{I}\cup\Pi_{C}$ to denote the disjoint union of these two sets.

Thus, $\sum_{j\in S_{\Pi}}K_{j}+K_{C}(S_{\Pi},P)\geq|\Pi_{I}|+|\Pi_{C}|=|\Pi|$, and so (7) implies (24) for any $\Pi$, and so Hall’s condition is satisfied, and a matching exists. Because in such a matching a set of vertices in $V_{M}$ matches to a set in $V_{V}$ of lower or equal cardinality, we have in particular that (20) holds for each $\Gamma\subseteq\Lambda$. $\Box$

Appendix D Proof of Theorem 4

Given the measurements $Y$ and measurement matrix $\Phi$, we will show that it is possible to recover some $P\in\mathcal{P}_{F}(X)$ and a corresponding vector $\Theta$ such that $X=P\Theta$ using the following algorithm: