Performance Bounds for Expander-based Compressed Sensing in the presence of Poisson Noise

TL;DR

This paper develops performance bounds for expander graph-based compressed sensing under Poisson noise, demonstrating their advantages over traditional dense sensing matrices through theoretical analysis and experiments.

Contribution

It introduces a novel expander graph-based sensing paradigm and a MAP recovery algorithm tailored for Poisson noise, with proven error bounds and improved performance.

Findings

Expander graphs outperform dense matrices in Poisson noise scenarios.

The proposed MAP algorithm effectively recovers sparse signals.

Experimental results confirm theoretical bounds and advantages.

Abstract

This paper provides performance bounds for compressed sensing in the presence of Poisson noise using expander graphs. The Poisson noise model is appropriate for a variety of applications, including low-light imaging and digital streaming, where the signal-independent and/or bounded noise models used in the compressed sensing literature are no longer applicable. In this paper, we develop a novel sensing paradigm based on expander graphs and propose a MAP algorithm for recovering sparse or compressible signals from Poisson observations. The geometry of the expander graphs and the positivity of the corresponding sensing matrices play a crucial role in establishing the bounds on the signal reconstruction error of the proposed algorithm. The geometry of the expander graphs makes them provably superior to random dense sensing matrices, such as Gaussian or partial Fourier ensembles, for the Poisson noise model. We support our results with experimental demonstrations.

Figures

Equations

$$(1-2ϵ)d\Vert x{\Vert }_1\le \Vert Ax{\Vert }_1\le d\Vert x{\Vert }_1$$

$$\Vert u{\Vert }_1\ge \Vert v{\Vert }_1-\Delta$$

$$\frac{(1-2ϵ)}{(1-6ϵ)}\left(2\Vert u_{\overline{S}}{\Vert }_1+\Delta \right)+\frac{2}{d(1-6ϵ)}\Vert Au-Av{\Vert }_1.$$

$$\Vert Au-Av{\Vert }_1+2dϵ\Vert y{\Vert }_1\ge (1-2ϵ)d\Vert y_S{\Vert }_1.$$

$$\Vert u{\Vert }_1$$

$$\ge \Vert u{\Vert }_1-2\Vert u_{\overline{S}}{\Vert }_1$$

$$+\Vert (u-v){\Vert }_1-2\Vert {(u-v)}_S{\Vert }_1-\Delta$$

$$\ge \Vert u{\Vert }_1-2\Vert u_{\overline{S}}{\Vert }_1-\Delta +\Vert u-v{\Vert }_1$$

$$-\frac{2\Vert Au-Av{\Vert }_1+4dϵ\Vert u-v{\Vert }_1}{(1-2ϵ)d}.$$

$$(\bm{I}_{\Lambda})_{i}=\begin{cases}1&\mbox{ if }i\in{\Lambda}\\ $0$&\mbox{otherwise}\end{cases}$$

$$\mathbf{y}\sim \text{mbox}Poisson\left(\mathbf{\Phi }{\mathbf{\alpha }}^{\ast }\right).$$

$$\mathrm{Pr}\left[Y_j|(\mathbf{\Phi }{\mathbf{\alpha }}^{\ast }{)}_j\right]=\{\begin{array}{ll}\frac{(\mathbf{\Phi }{\mathbf{\alpha }}^{\ast }{)}_j^{Y_j}}{Y_j!}{e}^{-(\mathbf{\Phi }{\mathbf{\alpha }}^{\ast }{)}_j}& \text{mbox}if(\mathbf{\Phi }{\mathbf{\alpha }}^{\ast }{)}_j\ne 0\\ \mathbf{\delta }(Y_j)& \text{mbox}else\end{array}$$

$$\mathbf{\delta }(Y_j)=\{\begin{array}{ll}1& \text{mbox}if{Y}_j=0\\ 0& \text{mbox}else\end{array}$$

$$\underset{(\mathbf{\Phi }{\mathbf{\alpha }}^{\ast }{)}_j\to 0}{\lim }\frac{(\mathbf{\Phi }{\mathbf{\alpha }}^{\ast }{)}_j^{Y_j}}{Y_j!}{e}^{-(\mathbf{\Phi }{\mathbf{\alpha }}^{\ast }{)}_j}=\delta (Y_j).$$

$$\Theta =\left\{f_1,\cdots ,f_{|\Theta |}\right\}$$

$$\forall f_i\in \Theta :\Vert f_i{\Vert }_1=1,f_i⪰0.$$

$$\Gamma =\left\{x_i=f_i+\lambda {\mathbf{I}}_{\Lambda }\right\}.$$

$$\sum _{\mathbf{x}\in \Gamma }{e}^{-\text{pen}(\mathbf{x})}\le 1.$$

$$\log \mathcal{L}(\mathbf{y}|\mathbf{x})$$

$$\hat{\mathbf{x}}\doteq \underset{\mathbf{x}\in \Gamma }{\mathrm{argmin}}\left\{\sum _{j=1}^m(\mathbf{\Phi }\mathbf{x}{)}_j-y_j\log {\left(\mathbf{\Phi }\mathbf{x}\right)}_j+2\text{pen}(\mathbf{x})\right\}$$

$$\forall \mathbf{x}\in \Gamma :\frac{m\lambda }{d}\le \Vert \mathbf{\Phi }\mathbf{x}{\Vert }_1\le 1+\frac{m\lambda }{d}.$$

$$\Vert \mathbf{\Phi }({\mathbf{\alpha }}^{\ast }-\hat{\mathbf{x}}){\Vert }_1^2$$

$$\Vert {\mathbf{\beta }}^{\ast }-\hat{\mathbf{\beta }}{\Vert }_1^2=\sum _{i=1}^m{\left|({\mathbf{\beta }}^{\ast }{)}_i-(\hat{\mathbf{\beta }}{)}_i\right|}^2$$

$$\le \sum _{i,j=1}^m{\left|({\mathbf{\beta }}^{\ast }{)}_i^{1/2}-(\hat{\mathbf{\beta }}{)}_i^{1/2}\right|}^2.{\left|({\mathbf{\beta }}^{\ast }{)}_j^{1/2}+(\hat{\mathbf{\beta }}{)}_j^{1/2}\right|}^2$$

$$\le 2\sum _{i=1}^m{\left|({\mathbf{\beta }}^{\ast }{)}_i^{1/2}-(\hat{\mathbf{\beta }}{)}_i^{1/2}\right|}^2.\sum _{j=1}^m\left|({\mathbf{\beta }}^{\ast }{)}_j+(\hat{\mathbf{\beta }}{)}_j\right|$$

$$\le 2\left(2+\frac{m\lambda }{d}\right)\sum _{i=1}^m{\left|({\mathbf{\beta }}^{\ast }{)}_i^{1/2}-(\hat{\mathbf{\beta }}{)}_i^{1/2}\right|}^2.$$

$$2\log \frac{1}{\left(\int \sqrt{p(\mathbf{Y}|\mathbf{g})p(\mathbf{Y}|\mathbf{h})}d\nu (\mathbf{y})\right)}=\sum _{j=1}^m\left((\mathbf{g}{)}_j^{1/2}-(\mathbf{h}{)}_j^{1/2}\right).$$

$${\mathbb{E}}_{\mathbf{Y}|\mathbf{\Phi }{\mathbf{\alpha }}^{\ast }}\left[\sum _{i=1}^m{\left|(\mathbf{\Phi }{\mathbf{\alpha }}^{\ast }{)}_i^{1/2}-(\mathbf{\Phi }\hat{\mathbf{x}}{)}_i^{1/2}\right|}^2\right]$$

$$\text{mbox}KL\left(p(\mathbf{Y}|\mathbf{\Phi }{\mathbf{\alpha }}^{\ast })\parallel p(\mathbf{Y}|\mathbf{\Phi }\mathbf{x})\right)\le \frac{d\Vert {\mathbf{\alpha }}^{\ast }-\mathbf{x}{\Vert }_1^2}{\lambda }.$$

$$\text{mbox}KL\left(p(\mathbf{Y}|\mathbf{\Phi }{\mathbf{\alpha }}^{\ast })\parallel p(\mathbf{Y}|\mathbf{\Phi }\mathbf{x})\right)$$

$$\mathbb{E}\left[\Vert \mathbf{\Phi }({\mathbf{\alpha }}^{\ast }-\hat{\mathbf{x}}){\Vert }_1\right]\le \sqrt{6}\underset{\tilde{\mathbf{x}}\in \Gamma }{\min }\sqrt{\frac{d}{\lambda }}\Vert {\mathbf{\alpha }}^{\ast }-\tilde{\mathbf{x}}{\Vert }_1+\sqrt{2\text{pen}(\tilde{\mathbf{x}})}$$

$$\mathbb{E}\left[\Vert \mathbf{\Phi }({\mathbf{\alpha }}^{\ast }-\hat{\mathbf{x}}){\Vert }_1^2\right]$$

$$\mathbb{E}\left[\Vert \mathbf{\Phi }({\mathbf{\alpha }}^{\ast }-\hat{\mathbf{x}}){\Vert }_1^2\right]\le 6\underset{\tilde{\mathbf{x}}\in \Gamma }{\min }\left(\frac{d}{\lambda }\Vert {\mathbf{\alpha }}^{\ast }-\tilde{\mathbf{x}}{\Vert }_1^2+2\text{pen}(\tilde{\mathbf{x}})\right).$$

$$\mathbb{E}\left[\Vert \mathbf{\Phi }({\mathbf{\alpha }}^{\ast }-\hat{\mathbf{x}}){\Vert }_1\right]\le \sqrt{6}\underset{\tilde{\mathbf{x}}\in \Gamma }{\min }\left(\sqrt{\frac{d}{\lambda }}\Vert {\mathbf{\alpha }}^{\ast }-\tilde{\mathbf{x}}{\Vert }_1+\sqrt{2\text{pen}(\tilde{\mathbf{x}})}\right).$$

Full text

I Introduction

The goal of compressive sampling or compressed sensing (CS) [1, 2] is to replace conventional sampling by a more efficient data acquisition framework, requiring fewer measurements whenever the measurement or compression is costly. This paradigm is particularly enticing in the context of photon-limited applications (such as low-light imaging) and digital fountain codes, since photo-multiplier tubes used in photon-limited imaging are large and expensive, and the number of packets transmitted via a digital fountain code is directly tied to coding efficiency. In these and other settings, however, we cannot directly apply standard methods and analysis from the CS literature, since these are based on assumptions of bounded, sparse, or Gaussian noise. Therefore, very little is known about the validity or applicability of compressive sampling to photon-limited imaging systems and streaming data communication.

The Poisson model is often used to model images acquired by photon-counting devices, particularly when the number of photons is small and a Gaussian approximation is inaccurate [3]. Another application is data streaming, in which streams of data are transmitted through a channel with Poisson statistics.

The Poisson model, commonly used to describe photon-limited measurements and discrete-time memoryless Poisson communication channels, pose significant theoretical and practical challenges in the context of CS. One of the key challenges is the fact that the measurement error variance scales with the true intensity of each measurement, so that we cannot assume uniform noise variance across the collection of measurements. The approach considered in this paper hinges, like most CS methods, on reconstructing a signal from compressive measurements by optimizing a sparsity-regularized data-fitting expression. In contrast to many CS approaches, however, we measure the fit of an estimate to the data using the Poisson log likelihood instead of a squared error term.

In previous work [4, 5], we showed that a Poisson noise model combined with conventional dense CS sensing matrices (properly scaled) yielded performance bounds which were somewhat sobering relative to bounds typically found in the literature. In particular, we found that if the number of photons (or packets) available to sense were held constant, and if the number of measurements, $m$, was above some critical threshold, then larger $m$ in general led to larger bounds on the error between the true and the estimated signals. This can intuitively be understood as resulting from the low signal-to-noise ratio of each of the $m$ measurements, which decays with $m$ when the number of photons (packets) is held constant.

This paper demonstrates that the bounds developed in previous work can be improved by considering alternatives to dense sensing matrices formed by making iid draws from a given probability distribution. In particular, we show that sensing matrices given by scaled adjacency matrices of expander graphs have important theoretical characteristics (especially an $\ell_{1}$ version of the restricted isometry property) which are ideally suited to controlling the performance of Poisson CS.

Expander graphs have been recently proposed as an alternative to dense random matrices within the compressed sensing framework, leading to computationally efficient recovery algorithms [6, 7, 8]. The approach described in this paper consists of the following key elements:

• •

expander sensing matrices and the RIP-1 associated with them;

• •

a reconstruction objective function which explicitly incorporates the Poisson likelihood;

• •

a collection of candidate estimators; and

• •

a penalty function defined over the collection of candidates which satisfies the Kraft inequality and which can be used to promote sparse reconstructions.

II Compressed Sensing using Expander Graphs

We start by defining an expander graph.

Definition II.1** (Expander Graph)**

A $(k,\epsilon)$-expander graph is a bipartite graph $V=(A,B),|A|=$ $n,~{}|B|=$ $m$, where $A$ is the set of variable nodes and $B$ is the set of parity (or check) nodes, which is unbalanced, i.e $m=o(n)$, and is left regular with left degree $d$, such that for any $S\subset A$ with $|S|\leq k$ the set of neighbors ${\cal N}(S)$ of $S$ has size $|{\cal N}(S)|>(1-\epsilon)d|S|$ .

Expander graphs have been recently proposed as a means of constructing efficient compressed sensing algorithms [6, 7, 8]. Figure 1 illustrates such a graph. The following proposition, proved using probabilistic methods, states that expander graphs are optimal in terms of the number of measurements required for compressive sampling:

Proposition II.1.1

For any $1\leq k\leq\frac{n}{2}$ and any positive $\epsilon$, there exists a $(k,\epsilon)$-expander graph with left degree $d=O\left(\frac{\log(\frac{n}{k})}{\epsilon}\right),$ and right set size $m=O\left(\frac{k\log(\frac{n}{k})}{\epsilon^{2}}\right).$

One reason why expander graphs are good sensing candidates is that the adjacency matrix of any expander graph almost preserves the $\ell_{1}$ norm of any sparse vector (RIP-1). Berinde et al have shown that the RIP-1 property can be derived from the expansion property [7]. In Section IV we exhibit the role this property plays in the performance of the maximum a posteriori (MAP) estimation algorithm for recovering sparse vectors in the presence of the Poisson noise.

Proposition II.1.2** (RIP-1 property of the expander graphs)**

Let $A$ be the $m\times n$ adjacency matrix of a $(k,\epsilon)$ expander graph $G$. Then for any $k$-sparse vector $x\in\mathbb{R}^{n}$ we have:

[TABLE]

The following theorem is a direct consequence of the RIP-1 property. This theorem states that, for any almost $k$-sparse vector111By “almost sparsity” we mean that the vector has at most $k$ significant entries. $u$, if there exists a vector $v$ whose $\ell_{1}$ norm is close to that of $u$, and if $v$ approximates $u$ in the measurement domain, then $v$ properly approximates $u$. In Section IV we show that the proposed MAP decoding algorithm outputs a vector satisfying the two conditions above, and hence approximately recovers the desired signal.

Theorem II.2

Let $A$ be the adjacency matrix of a $(k,\epsilon)$-expander and $u,v$ be two vectors in $\mathbb{R}^{n}$, such that

[TABLE]

for some positive $\Delta$. Let $S$ be the set of $k$ largest (in magnitude) coefficients of $u$, and $\overline{S}$ be the set of remaining coefficients. Then $\|u-v\|_{1}$ is upper-bounded by

[TABLE]

Proof:

Let $y=u-v$, and $\langle S_{1},\cdots,S_{t}\rangle$ be a decreasing partitioning of $\overline{S}$ (with respect to coefficient magnitudes), such that all sets but (possibly) $S_{t}$ have size $k$. Note that $S_{0}=S$. Let $\tilde{A}$ be a submatrix of $A$ containing rows from ${\cal N}(S)$. Then, following the argument of Berinde et al. [7], we have the following inequality:

[TABLE]

Now, using the triangle inequality and Eq. (2), we obtain

[TABLE]

Rearranging the inequality completes the proof. ∎

Finally, note that, since the graph is regular, there exists a minimal set ${\Lambda}$ of variable (left) nodes with size at most $m$, such that its neighborhood covers all of the check nodes, i.e ${\cal N}({\Lambda})={B}$. Let $\bm{I}_{\Lambda}$ be an index vector such that

[TABLE]

where $(\bm{I}_{\Lambda})_{i}$ denotes the $i$th entry of $\bm{I}_{\Lambda}$. Then $A\bm{I}_{\Lambda}\succeq I_{m\times 1}.$ The role of ${\Lambda}$ is to guarantee that recovery candidates are non-zero vectors in the measurement domain. This is crucial in compressed sensing with Poisson noise, and we will explain this issue in detail in the next sections.

III Compressed Sensing in the presence of Poisson Noise

Recall that a signal is defined to be “almost $k$-sparse” if it has at most $k$ significant entries, while the remaining entries have near-zero values. Let $\bm{\alpha}^{*}_{k}$ be the best $k$-term approximation of $\bm{\alpha}^{*}$, and $\bm{\Phi}_{m\times n}$ be the sensing matrix. Let $\mathbb{Z}_{+}=\{0,1,2,\cdots\}$. We assume that each entry of the measured vector $\bm{y}\in\mathbb{Z}_{+}^{m}$ is sensed independently according to a Poisson model:

[TABLE]

That is, for each index $j$ in $\{1,\cdots,m\}$, the random variable $Y_{j}$ is sampled from a Poisson distribution with mean $(\bm{\Phi}\bm{\alpha}^{*})_{j}$:

[TABLE]

where

[TABLE]

Note that

[TABLE]

We use MAP (maximum *a posteriori *probability) decoding for recovering a good estimate for $\bm{\alpha}^{*}$, given measurements $\bm{y}$ in the presence of the Poisson noise. Let

[TABLE]

be a set of candidate estimates for $\bm{\alpha}^{*}$ such that

[TABLE]

We would like to find the best possible *a posteriori *estimate, given the observation vector $\bm{y}$. Moreover, to maintain consistency between the maximum likelihood and the MAP decoding, we impose the requirement that no candidate MAP estimator can have a zero coordinate if the corresponding measurement is non-zero. To guarantee this, let $\lambda\ll\frac{1}{k\log n}$ be a small parameter. We define

[TABLE]

Then since $\lambda$ is strictly positive, we will have $\Phi x\succ 0$ for any estimate $x$ in $\Gamma$. This allows us to run the MAP decoding over the set $\Gamma$ and output the (one-to-one) corresponding estimate from $\Theta$. We show this precisely in the next section. This relaxation allows the MAP decoding to work properly and guarantees recovering an estimate from $\Theta$ with expected $\ell_{1}$ error close to the error of the best estimate in $\Theta$.

Let $\text{pen}(x)$ be a nonnegative penalty function based on our prior knowledge about the estimates in $\Gamma$ (or equivalently let $\widehat{\text{pen}}(\theta)$ be a penalty function over $\Theta$). The only constraint that we impose on the penalty function is the Kraft inequality

[TABLE]

For instance, we can impose less penalty on sparser signals or construct a penalty based on any other prior knowledge about the underlying signal. The log-likelihood of the measurement, according to Eq. (3), is

[TABLE]

We will show that the maximum *a posteriori *estimate

[TABLE]

has error close to the error of the best estimate in $\Gamma$. The decoding in (6) is a MAP algorithm over the set of estimates $\Gamma$, where the likelihood is computed according to the Poisson model (3) and the penalty function corresponds to a negative log prior on the candidate estimators in $\Gamma$.

IV Performance of MAP Recovery on Almost Sparse Signals

Let $A$ be the $m\times n$ adjacency matrix of a $(2k,1/16)$-expander with left degree $d$. Also let $\bm{\Phi}=\frac{A}{d}$ be the sensing matrix. From definition of $\bm{I}_{\Lambda}$ and $\Gamma$, and since the adjacency matrix of any graph only consists of zeros and ones, for any estimate $\bm{x}\in\Gamma$ we have $\bm{\Phi}\bm{x}\succeq\frac{\lambda}{d}$. Moreover, from the RIP-1 property of the expander graphs stated in Lemma (II.1.2) we know that for any signal $\bm{x}$, $\|\bm{\Phi}\bm{x}\|_{1}\leq\|\bm{x}\|_{1}$, and $(1-2\epsilon)\|\bm{x}\|_{1}\leq\|\bm{\Phi}\bm{x}\|_{1}$ for any $k$-sparse signal $\bm{x}$. Hence by definition of $\Gamma$

[TABLE]

Lemma IV.1

Let $\bm{\Phi}$ be the normalized expander sensing matrix, $\bm{\alpha}^{*}$ be the original $k$-sparse signal and $\hat{\bm{x}}$ be the minimizer of the Equation (6). Then

[TABLE]

Proof:

Let ${\bm{\beta^{*}}}=\bm{\Phi}\bm{\alpha}^{*}$ and ${\bm{\hat{\beta}}}=\bm{\Phi}\hat{\bm{x}}$. Then

[TABLE]

The first and the second inequalities are by Cauchy–Schwarz, while the third inequality is a consequence of the RIP-1 property of the expander graphs (Lemma II.1.2) and Eq. (7). ∎

Lemma IV.2

Given two Poisson parameter vectors $\bm{g},\bm{h}\in\mathbb{R}^{m}_{+}$, the following equality holds:

[TABLE]

Proof:

The proof follows from expanding the term $\int{\sqrt{p(\bm{Y}|\bm{g})p(\bm{Y}|\bm{h})}d\nu(\bm{y})}$, and is provided in [4]. ∎

Lemma IV.3

Let $\bm{\Phi}$ be the expander sensing matrix, $\bm{\alpha}^{*}$ be the original almost $k$-sparse signal, and $\hat{\bm{x}}$ be a minimizer in Eq. (6). Finally let $\bm{y}$ be the compressive measurements of $\bm{\alpha}^{*}$ in Poisson model. Then

[TABLE]

Proof:

The proof exploits techniques from Li and Baron [9], and Kolaczyk and Nowak [10]. ∎

Now we show that in Poisson setting for all estimates $\bm{x}$ in $\Gamma$, the relative entropy term $\mbox{KL}\left(p(\bm{Y}|\bm{\Phi}\bm{\alpha}^{*})\parallel p(\bm{Y}|\bm{\Phi}\bm{x})\right)$ is upper bounded by the squared $\ell_{1}$ norm of $\bm{\alpha}^{*}-\bm{x}$:

Lemma IV.4

For any estimate $\bm{x}\in\Gamma$ the following inequality holds:

[TABLE]

Proof:

[TABLE]

The first inequality is $\log t\leq t-1$, and the second inequality is by the RIP-1 property of $\bm{\Phi}$ (Eq. (1)) and definition of $\Gamma$ (Eq. (4)). ∎

Lemma IV.5

Let $\bm{\Phi}$ be the expander sensing matrix, $\bm{\alpha}^{*}$ be the original almost $k$-sparse signal, and $\hat{\bm{x}}$ be a minimizer in Eq. (6). Then

[TABLE]

Proof:

Lemmas IV.1, IV.3, and IV.4 together imply

[TABLE]

Since $\lambda\ll\frac{1}{k\log(n/k)}$, and $m=O\left(k\log(n/k)\right)$, the ratio $\frac{m\lambda}{d}$ is much less than $1$. So $2+\frac{m\lambda}{d}\ll 3$, and

[TABLE]

Now since the function $f(\bm{x})=\|A\bm{x}+b\|_{1}^{2}$ is convex and the square root function is strictly increasing, by applying Jensen’s inequality we get

[TABLE]

∎

Theorem IV.6

Let $\bm{\Phi}$ be the expander sensing matrix, $\lambda\ll\frac{1}{k\log(n/k)}$ be a small positive value, $\bm{\alpha}^{*}$ be the original almost $k$-sparse signal compressively sampled in the presence of Poisson noise, $\hat{\bm{x}}$ be a minimizer in Eq. (6), and $\hat{\bm{f}}$ be the corresponding estimate in $\Theta$, i.e $\hat{\bm{x}}=\hat{\bm{f}}+\lambda\bm{I}_{\Lambda}$. Then

[TABLE]

Proof:

In Lemma IV.5, we have bounded $\|\bm{\Phi}(\bm{\alpha}^{*}-\hat{\bm{x}})\|_{1}$. Now we can use Theorem II.2 to bound $\|\bm{\alpha}^{*}-\hat{\bm{x}}\|_{1}$. We have used a $(2k,1/16)$-expander. Also since $\|\bm{\alpha}^{*}\|_{1}=1$, and any $x$ in $\Gamma$ has the form $\theta+\lambda\bm{I}_{\Lambda}$ where $\|\theta\|_{1}=1$, and $\|\bm{I}_{\Lambda}\|_{1}=m$, and since $\hat{\bm{x}}\in\Gamma$, we get $\|\hat{\bm{x}}\|_{1}\leq\|\hat{\bm{f}}\|_{1}+\lambda\|\bm{I}_{\Lambda}\|_{1}$ and hence

[TABLE]

As a result, by Theorem II.2 and Lemma IV.5 we get

[TABLE]

Consequently, we have derived a bound on how much $\hat{\bm{x}}$ differs from $\bm{\alpha}^{*}$. Since any $\bm{x}$ in $\Gamma$ has the form $\bm{f}+\lambda\bm{I}_{\Lambda}$ for some estimate $\bm{f}$ in $\Theta$, using the triangle inequality we get

[TABLE]

and so

[TABLE]

∎

By substituting the values $m=O\left(k\log(n/k)\right)$, and $d=O\left(\log(n/k)\right)$, and choosing

[TABLE]

we can guarantee that $\mathbb{E}\left[\|\bm{\alpha}^{*}-\hat{\bm{f}}\|_{1}\right]$ is of order

[TABLE]

Remark IV.7

It has been shown by Willett et.al. [4, 5] that, using random dense matrices, the MAP reconstruction algorithm can reconstruct a signal $\bm{\alpha}^{*}$ satisfying $\|\bm{\alpha}^{*}\|_{1}=1$ with the expected error of

[TABLE]

Hence, for random dense matrices there is an $O\left(m^{-1}\right)$ min-max approximation error. This error cannot be made arbitrarily small by increasing the number of measurements as the first term in (11) also depends on $m$. However, as stated earlier, the bounds of [4, 5] are not restricted to signals that are sparse in the canonical basis.

V Experimental Results

To validate our results via simulation, we generated random sparse signals, simulated Poisson observations of the signal multiplied by the proposed expander graph sensing matrix, and reconstructed the signal using the proposed objective function in (6).

Each signal was a length $n=100,000$ signal with $k$ non-zero elements, where $k$ ranged from $1$ to $4,000$. Each of the non-zero elements was assigned intensity $I$, where $I$ was $10$, $100$, $1,000$, or $10,000$. The locations of the non-zero elements were selected uniformly at random for each trial. The sensing matrix was a scaled adjacency matrix of an expander graph, as described earlier, with $d=16$ and the number of rows $m=40,000$.

Reconstruction was performed using a method described in [11] for reconstruction of sparse signals from indirect Poisson measurements, precisely the situation encountered here. The penalty function used in this implementation is proportional to $\|\bm{x}\|_{1}$; constructing a penalty function of this form which satisfies the Kraft inequality is a subject of ongoing work. (The authors would like to thank Mr. Zachary Harmany for his assistance with the implementation of this algorithm.) After each trial, the normalized $\ell_{1}$ error was computed as $\|\bm{\alpha}^{*}-\hat{\bm{x}}\|_{1}/\|\bm{\alpha}^{*}\|_{1}$, and the errors were averaged over $50$ trials. The results of this experiment are presented in Figure 2.

VI Conclusions

In this paper we investigated the advantages of expander-based sensing over dense random sensing in the presence of Poisson noise. Even though Poisson model is essential in some applications, dealing with this noise model is challenging as the noise is not bounded, or even as concentrated as Gaussian noise, and is signal-dependent.

Here we proposed using normalized adjacency matrices of expander graphs as an alternative construction of sensing matrices, and we showed that the binary nature and the RIP-1 property of these matrices yield provable consistency for a MAP reconstruction algorithm.

Acknowledgements

The authors would like to thank Zachary Harmany for his assistance with the implementation of the reconstruction algorithm, and Piotr Indyk for his insightful comments on the performance of the expander graphs.