A robust parallel algorithm for combinatorial compressed sensing

TL;DR

This paper introduces the Robust-$\ell_0$ decoding algorithm, enhancing parallel combinatorial compressed sensing by effectively recovering sparse signals from noisy measurements, outperforming existing methods under Gaussian noise conditions.

Contribution

The paper develops a robust decoding algorithm that extends previous work by incorporating noise resilience through posterior distribution approximation.

Findings

Robust-$\ell_0$ outperforms existing algorithms in noisy environments.

The method effectively recovers sparse signals with Gaussian noise.

Analytic expressions for posteriors are derived under certain assumptions.

Abstract

In previous work two of the authors have shown that a vector $x \in \mathbb{R}^n$ with at most $k < n$ nonzeros can be recovered from an expander sketch $Ax$ in $\mathcal{O}(\mathrm{nnz}(A)\log k)$ operations via the Parallel-$\ell_0$ decoding algorithm, where $\mathrm{nnz}(A)$ denotes the number of nonzero entries in $A \in \mathbb{R}^{m \times n}$. In this paper we present the Robust-$\ell_0$ decoding algorithm, which robustifies Parallel-$\ell_0$ when the sketch $Ax$ is corrupted by additive noise. This robustness is achieved by approximating the asymptotic posterior distribution of values in the sketch given its corrupted measurements. We provide analytic expressions that approximate these posteriors under the assumptions that the nonzero entries in the signal and the noise are drawn from continuous distributions. Numerical experiments presented show that Robust-$\ell_0$ is superior to existing greedy and combinatorial compressed sensing algorithms in the presence of small to moderate signal-to-noise ratios in the setting of Gaussian signals and Gaussian additive noise.

Tables

Table 1. TABLE I: Extensions of scores used in Expander ℓ 0 subscript ℓ 0 \ell_{0} -decoding to identify candidate updates to the sparse signal x ^ ^ 𝑥 \hat{x} .

Table 2. TABLE II: Variants of Robust- ℓ 0 subscript ℓ 0 \ell_{0}

Algorithm	adaptive_k	quantised
Robust-${\ell }_0$	No	No
Robust-${\ell }_0$-adaptive	Yes	No
Robust-${\ell }_0$-quantised	No	Yes
Robust-${\ell }_0$-adaptive-quantised	Yes	Yes

Equations

$$\underset{x\in {\mathbb{R}}^n}{\min }{\Vert x\Vert }_0\text{ s.t. }Ax=y.$$

$$\Big{|}\Big{\{}i\in[m]:\sum_{j\in S}\operatorname{\mathds{1}}_{|A_{i,j}|>0}\Big{\}}\Big{|}>(1-\varepsilon)d|S|$$

$$\sum _{j\in S_1}{x}_j\ne \sum _{j\in S_2}{x}_j\forall S_1,S_2\subset \mathrm{supp}(x)\text{ s.t. }{S}_1\ne S_2.$$

$$n_e$$

$$n_z$$

$$p_z\left(\omega \right):=\mathbb{P}\left(y_i=0\mid {\hat{y}}_i=\omega \right).$$

$$p_e\left(\omega \right):=\mathbb{P}\left(y_{i_1}=y_{i_2}\mid {\hat{y}}_{i_1}-{\hat{y}}_{i_2}=\omega \right)$$

$$p_z\left(\omega \right)$$

$$p_e\left(\omega \right)$$

$${\breve{p}}_e(\omega )=\frac{p_e(\omega )}{{\max }_s{p}_e(s)},{\breve{p}}_z(\omega )=\frac{p_z(\omega )}{{\max }_s{p}_z(s)}.$$

$$n_e$$

$$n_z$$

$$\left|\{i\in \mathcal{N}(j):y_i=\omega \}\right|\ge (1-2\epsilon )d\forall (j,\omega )\in T$$

$${\mu }_0(E)$$

$${\mu }_1(E)$$

$${\mu }_{q+1}(E)$$

$${\mu }^{-}(E)=\mu (-E),\forall E\in \mathcal{B}(\mathbb{R})$$

$$\bar{\mu }(E)=\frac{\mu (E)+\mu (-E)}{2},\forall E\in \mathcal{B}(\mathbb{R})$$

$$\tilde{\mu }(E)=(\mu \ast {\mu }^{-})(E),\forall E\in \mathcal{B}(\mathbb{R}).$$

$$s_n:=\sum _{j=1}^n{b}_j{x}_j$$

$$s\sim \exp (-p)\sum _{q\ge 0}\frac{p^q}{q!}{\mu }_q.$$

$${\psi }_{s_n}(t)$$

$$=\prod _{j=1}^n\mathbb{E}\left[\exp (it{b}_j{x}_j)\right]$$

$$={\left(\left(1-\frac{p_n}{n}\right)+\left(\frac{p_n}{n}\right){\psi }_x(t)\right)}^n$$

$$={\left(1+\frac{p_n({\psi }_x(t)-1)}{n}\right)}^n.$$

$$\underset{n\to \infty }{\lim }{\psi }_{s_n}(t)=\exp (-p)\exp (p{\psi }_x(t)).$$

$${\left[{\psi }_x(t)\right]}^q={\psi }_{w_q}(t).$$

$$z\sim \exp (-p)\sum _{q\ge 0}\frac{p^q}{q!}{\mu }_q.$$

$${\psi }_z(t)$$

$$=\exp (-p)\sum _{q\ge 0}\frac{p^q}{q!}{\psi }_{w_q}(t)$$

$$=\exp (-p)\sum _{q\ge 0}\frac{p^q}{q!}{\left({\psi }_x(t)\right)}^q$$

$$=\exp (-p)\sum _{q\ge 0}\frac{{\left(p{\psi }_x(t)\right)}^q}{q!}$$

$$=\exp (-p)\exp (p{\psi }_x(t)).$$

$${\hat{y}}_i\stackrel{(d)}{\to }{\hat{y}}_i^{\ast }\text{mbox}and{\hat{y}}_{i_1}-{\hat{y}}_{i_2}\stackrel{(d)}{\to }{\hat{g}}^{\ast }$$

$${\hat{y}}_i^{\ast }$$

$${\hat{g}}^{\ast }$$

$$y_i=\sum _{j=1}^n{A}_{i,j}{x}_j.$$

Full text

I Introduction

Compressed sensing is a well studied method by which a sparse or compressible vector can be acquired by a number of measurements proportional to the number of its dominant entries [2], [3]. To fix notation, let $\chi_{k}^{n}$ be the set of vectors in $\mathbb{R}^{n}$ that have at most $k$ nonzero entries, let $x\in\chi_{k}^{n}$ and let $A\in\mathbb{R}^{m\times n}$ be a matrix with $m<n$. We will refer to $A$ as the measurement matrix, $x$ as the signal and $y=Ax$ as the measurements. The goal of compressed sensing is to recover the sparsest, most parsimonious, $x\in\mathbb{R}^{n}$ from the measurements $y$ and the matrix $A$. Letting $\|\cdot\|_{0}$ denote the number of non-zeros in $x$, the problem of finding $x$ can be written as

[TABLE]

Many algorithms have been developed to solve this problem or equivalent formulations and there are good theoretical results on when and how fast recovery of a signal is possible given certain types of measurement matrix $A$ and signal $x$. These algorithms can be broadly categorized into convex optimization based algorithms like those implemented in [4, 5, 6, 7] and greedy algorithms [8, 9, 10, 11, 12, 13, 14], and were designed and analysed for the setting of dense sensing matrices; e.g. independent (sub-)Gaussian entries or randomly subsampled Fourier matrices. For a more detailed introduction to compressed sensing see [15].

Here we extend an algorithm proposed in [1] which can be used to recover exactly a sparse signal from its expander sketch (see Section II for details). Specifically, [1] proposed Parallel-$\ell_{0}$ (Algorithm 2), for noiseless combinatorial compressed sensing which is guaranteed to converge in $\mathcal{O}(\operatorname{nnz}(A)\log(k))$ where the sensing matrix $A$ is an expander matrix (Definition I.1) and the signal $x\in\chi_{k}^{n}$ is dissociated in the sense of Definition I.2 or the signal is drawn independently of $A$. For alternative combinatorial compressed sensing algorithms see, for example, [16, 17, 18, 19, 20, 21]. We borrow notation from combinatorics and use the shorthands $[n]:=\{1,\dots,n\}$, $[n]^{(k)}:=\{S\subset[n]:|S|=k\}$ where $|S|$ denotes the cardinality of the set $S$, and $[n]^{(\leq k)}:=\cup_{\ell\leq k}[n]^{(\ell)}$ for $n,k\in\mathbb{N}$ and $k<n$.

Definition I.1** (Expander matrices [1]).**

The matrix $A\in\{0,1\}^{m\times n}$ is a $(k,\varepsilon,d)$-expander matrix if $\sum_{i=1}^{m}\operatorname{\mathds{1}}_{|A_{i,j}|>0}=d$ for all $j\in[n]$ and

[TABLE]

for all $S\in[n]^{(\leq k)}$. We denote by $\mathbb{E}_{k,\varepsilon,d}^{m\times n}$ the set of $(k,\varepsilon,d)$-expander matrices of dimension $m\times n$.

Definition I.2** (Dissociated signals [1]).**

A signal $x\in\mathbb{R}^{n}$ is said to be dissociated if

[TABLE]

An example of (almost surely) dissociated signals are those drawn from a continuous distribution. It is shown in [1] that if $y=Ax$ is an expander sketch and $x\in\chi_{k}^{n}$ is dissociated, then there exists a subset $T\subset[n]$ such that, for each $j\in T$, $|\{i\in[m]:y_{i}=x_{j}\}|$ is bounded below by a positive constant depending on $d$ and $\varepsilon$. This guarantees that if $|\{i\in\operatorname{supp}(a_{j}):y_{i}=y_{\ell}\}|>d/2$ then $x_{j}=y_{\ell}$. Parallel-$\ell_{0}$ (Algorithm 2) implements this observation by letting $\hat{x}=0$ and estimating the decrease in $\|y\|_{0}$ when performing the update $\hat{x}_{j}\leftarrow\hat{x}_{j}+y_{\ell}$. We denote for $j\in[n]$ its neighbours by $\mathcal{N}(j)\vcentcolon=\{i\in[m]:|A_{ij}|>0\}$; to estimate the decrease in $\|y\|_{0}$, Parallel-$\ell_{0}$ computes

[TABLE]

We extend their approach to the additive noise signal model of $\hat{y}=y+\eta$ with $y=Ax$ and $\eta\in\mathbb{R}^{m}$ by replacing (1)-(2) with scores estimating the distribution of $n_{e}$ and $n_{z}$, e.g. (8)-(9). That is, we follow a Bayesian approach to the computation of these scores and estimate:

The probability of $y_{i}=0$ given that we observe $\hat{y}_{i}$.

[TABLE] 2. 2.

The probability of $y_{i_{1}}=y_{i_{2}}$ given that we observe $\hat{y}_{i_{1}}-\hat{y}_{i_{2}}$.

[TABLE]

Among our contributions are series approximations for (3)-(4) when the signals and measurements are generated according to the generating model given in Definition I.3 and illustrated in Figure 1. In what follows, we let $\mathbb{D}(\mathbb{R})$ be the set of distributions supported on $\mathbb{R}$. If $\mu\in\mathbb{D}(\mathbb{R})$, we write $z\sim\mu$ to denote that $z$ was drawn according to the distribution $\mu$. We also use the notation $v_{i}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\mu$ to denote that each $v_{i}$ is drawn independently at random from $\mu$. Finally, we use $U(S)$ to denote the uniform distribution over a set $S$.

Definition I.3** (Generating model $\operatorname{GM}(n,m,k,d,\mu,\nu)$).**

Let $n,m,k,d\in\mathbb{N}$ be such that $k<m<n$ and $d\ll m$. Let $\mu,\nu\in\mathbb{D}(\mathbb{R})$. Then, the problem $(A,\hat{y})$ is drawn from the model $\operatorname{GM}(n,m,k,d,\mu,\nu)$ if $A\in\{0,1\}^{m\times n}$ and $\hat{y}\in\mathbb{R}^{m}$ are such that

each column of $A$ has a support drawn uniformly at random from $[m]^{(d)}$; 2. 2.

$\operatorname{supp}(x)$ is drawn uniformly at random from $[n]^{(k)}$; 3. 3.

$x_{j}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\mu$ for each $j\in\operatorname{supp}(x)$; 4. 4.

$\eta_{i}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\nu$ for each $i\in[m]$; 5. 5.

$\eta_{i}$ is independent of $x_{j}$ for all $i\in[m],j\in\operatorname{supp}(x)$; 6. 6.

$y=Ax$ and $\hat{y}=y+\eta$.

We write $(A,\hat{y})\sim\operatorname{GM}(n,m,k,d,\mu,\nu)$ to denote problem instances drawn from this signal model.

Remark**.**

It is important to note that a matrix $A\in\mathbb{R}^{m\times n}$ generated under the model presented in Definition I.3 and Figure 1, is a $(k,\varepsilon,d)$-expander matrix with high probability, see [22], [15, Theorem 13.6].

Moreover, the generating model in Definition I.3 also allows us to define robust estimates for (3)-(4) for general noise and signal distributions and to any degree of accuracy under the assumption that these probability measures are available. From there we can define noisy analogues to the values $n_{e}$ and $n_{z}$ in (1)-(2) used in Parallel-$\ell_{0}$ but which are robust to additive noise. The contributions of this work are two-fold: (i) to present principled ways to compute (3)-(4) in the case where the nonzeros in $\eta$ and $x$ are drawn from continuous probability distributions; (ii) to provide a variation of Parallel-$\ell_{0}$ that is robust to noise. While other similar generating models can be considered using the techniques presented here, for ease of exposition and clarity, we restrict our description to this model.

Theorem I.4** (Probabilities for general signal and noise distributions).**

Let $\delta,\rho\in(0,1)$. For each $n>1$, let $m=\delta n$, $k=\rho m$ and $d<m$. If $\mu,\nu\in\mathbb{D}(\mathbb{R})$ and $(A,\hat{y})\sim\operatorname{GM}(n,m,k,d,\mu,\nu)$. Then as $n\rightarrow\infty$,

[TABLE]

Where $\mu_{q},\bar{\mu}_{q},\tilde{\nu}$ are probability measures constructed as in Definition III.1.

Equations (5) and (6) allows us to quantify the uncertainty associated with computing the score for Parallel-$\ell_{0}$ under the presence of additive noise. Note that equations (5) and (6) can be easily adapted to alternative generative models, such as where the expected density of nonzeros per row varies, but for expository clarity we restrict our discussion to this somewhat generic model. It will be discussed in Section III-C that equations (5) and (6) should not be used directly, but instead be scaled by considering normalised functions $\breve{p}_{e}$ and $\breve{p}_{z}$ defined as,

[TABLE]

In the most general case $n_{e}$ and $n_{z}$ can be written as the sum of individual scores $q_{e}$ and $q_{z}$ as follows

[TABLE]

A confidence threshold $t>0$ needs to be given for some variants of our algorithm, so to simplify the exposition we include this parameter in the scores $q_{e}(\cdot\mid t)$ and $q_{z}(\cdot\mid t)$ regardless of whether it is used or not. A summary of the score functions are given in Table I.

I-A Outline of the manuscript

The structure of this paper is as follows: Section II comprises a review of combinatorial compressed sensing and the Parallel-$\ell_{0}$ decoding algorithm extended here. Robust-$\ell_{0}$ decoding and the associated scores (8)-(9) are presented in Section III. In Section IV we present numerical experiments which demonstrate Robust-$\ell_{0}$ to perform superior to a number of leading greedy and combinatorial compressed sensing algorithms.

II Background: Combinatorial Compressed Sensing and $\ell_{0}$-decoding

As mentioned in the previous section, the branch of combinatorial compressed sensing measures $x\in\chi_{k}^{n}$ with the adjacency matrix of an expander graph. These matrices are of very low complexity in terms of generation and storage, and also promise faster encoding and decoding than their dense counterparts, see Theorem II.2. In this section we review the basic elements of expander graphs and combinatorial compressed sensing.

There have been various algorithms proposed to reconstruct a sparse vector $x$ from measurements $y=Ax$ when $A$ is an expander matrix, see [17], [18], [20] and [19]; the work presented here starts with the Parallel-$\ell_{0}$ algorithm and to improve this algorithm by making it robust to noisy measurements, results in Robust-$\ell_{0}$ (Algorithm 1). The key observation for the Parallel-$\ell_{0}$ algorithm is given by the following Lemma.

Lemma II.1**.**

Let $y=Ax$, $x$ dissociated, $A\in\mathbb{E}_{k,\varepsilon,d}^{m\times n}$ with $\varepsilon<\frac{1}{4}$. Then there exists a nonempty set $T\subset[n]\times\mathbb{R}$ such that

[TABLE]

and for every tuple in $T$ that satisfies this property, we have $w=x_{j}$.

This means that at each iteration, if the residual $r$ is non-zero, i.e. if we have not yet found the correct $x$, then there is a set of entries in $x$ that we can change so that we reduce the number of non-zeros in $r$ by at least $|T|(1-2\varepsilon)d$.

Theorem II.2** (Convergence of Algorithm 2 [1]).**

Let $A\in\mathbb{E}_{k,\varepsilon,d}^{m\times n}$ and let $\varepsilon\leq\frac{1}{4}$, and $x\in\chi_{k}^{n}$ be dissociated. Then, Parallel-$\ell_{0}$ with $\alpha=(1-2\varepsilon)d$ can recover $x$ from $y=Ax$ in $\mathcal{O}(\log k)$ iterations of complexity $\mathcal{O}(dn)$.

To put this result into context and show its applicability, we recall the remark after Definition I.3, stating that random matrices as considered in this work are indeed expander matrices with high probability.

We furthermore emphasize the fact that the algorithm is designed in a way that allows for massively parallel implementations.

III Main contributions: $\ell_{0}$-decoding for Noisy Measurements

We now consider the case where the measurements $y$ are subject to additive noise, i.e. instead of $y=Ax$, we measure $\hat{y}=y+\eta$ where $\eta$ is a realization of a random variable with $\eta_{i}\sim\nu$.

Parallel-$\ell_{0}$ is not able to cope with additive noise, as it needs to make decisions whether a value in $\hat{y}$ is zero and whether two values in $\hat{y}$ are equal to each other. While for very small noise levels we could consider two values as equal if they are within a certain number of standard deviations, for larger noise levels the decision becomes more challenging. A discussed in Section I we need to know $p_{z}(\hat{y}_{i})$ and $p_{e}(\hat{y}_{i_{1}}-\hat{y}_{i_{2}})$ which correspond, respectively, to the probability of $y_{i}=0$ given that we observe $\hat{y}_{i}$ and the probability of $y_{i_{1}}=y_{i_{2}}$ given that we observe $\hat{y}_{i_{1}}-\hat{y}_{i_{2}}$. The functions $p_{z}$ and $p_{e}$ depend on the parameters fed into the generative model in Definition I.3 and in particular on the distribution of $y$ and $\hat{y}$. Hence, since $y$ is a vector of sparse inner products we need to understand the limiting behaviour of sparse sums.

We include Definition III.1 in order to remind the reader of some actions on measures used in this manuscript as so as to be relatively self contained.

Definition III.1** (Measures [23]).**

Let $\mathcal{B}(\mathbb{R})$ denote the Borel $\sigma$-algebra over $\mathbb{R}$. If $E\in\mathcal{B}(\mathbb{R})$ let $-E:=\{-x:x\in E\}$. Let $\mu\in\mathbb{D}(\mathbb{R})$, we define the following measures.

The $q$-convolution,

[TABLE] 2. 2.

The negative measure,

[TABLE] 3. 3.

The symmetrized measure,

[TABLE] 4. 4.

The measure associated with the difference of two random variables,

[TABLE]

Lemma III.2** (Limiting distribution for sparse sums of random variables).**

Let $p\in(0,1)$, let $\mu\in\mathbb{D}(\mathbb{R})$ and let $\mu_{q}\in\mathbb{D}(\mathbb{R})$ be its the $q$-fold convolution. For each $n\geq 1$, let

[TABLE]

be such that,

$x_{j}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\mu$ for each $j\in[n]$, 2. 2.

$b_{j}\stackrel{{\scriptstyle\text{i.i.d.}}}{{\sim}}\mathrm{Ber}(\frac{p_{n}}{n})$ for each $j\in[n]$ with $p_{n}\to p\in\mathbb{R}$ as $n\to\infty$.

Then, as $n\to\infty$ it holds that $s_{n}\overset{(d)}{\to}s$ where

[TABLE]

Proof.

Let $\psi_{s_{n}}(t)$ be the characteristic function of $s_{n}$. Let $x\sim\mu$, then

[TABLE]

Taking the limit $n\to\infty$ we see that

[TABLE]

Letting $w_{q}=\sum_{j=1}^{q}x_{j}$ it holds by the independence of $\{x_{1},\dots,x_{q}\}$ that $w_{q}\sim\mu_{q}$ and

[TABLE]

Now, consider a random variable $z$ distributed according to

[TABLE]

The characteristic function of $z$ is given by

[TABLE]

Therefore (15) equals (12). By Lévy’s continuity Theroem pointwise convergence of the characteristic functions implies weak convergence of the random variables (cf. [23, Theorem 15.23]) and hence the statement follows. ∎

Theorem III.3** (Distribution of $\hat{y}_{i}$ and $\hat{y}_{i_{1}}-\hat{y}_{i_{2}}$).**

Fix $\delta,\rho\in(0,1)$ and for $n\in\mathbb{N}$ let $m=\delta n$, $k=\rho m$ and $d\ll m$. Furthermore, let $\mu$ and $\nu$ be measures and assume that $\hat{y}=y+\eta=Ax+\eta$ is drawn from the model $\mathrm{GM}(n,m,k,d,\mu,\nu)$. Then as $n\to\infty$

[TABLE]

where

[TABLE]

Proof.

To show (16), let $i\in[m]$ and

[TABLE]

By our assumptions on $A$ and $x$,

[TABLE]

where for two events $E_{1}$ and $E_{2}$ we let $E_{1}\wedge E_{2}$ be the conjunction of the events. Note that if $A_{i,j}x_{j}\neq 0$, then $j\in\operatorname{supp}(x)$ so $A_{i,j}x_{j}=x_{j}\sim\mu$. Hence, letting $b_{j}\sim\mathrm{Ber}\left(\frac{d\rho}{n}\right)$,

[TABLE]

Invoking Lemma III.2 with $p_{n}=p=d\rho$ we obtain that $y_{i}\rightarrow y_{i}^{*}$ as $n\to\infty$ where

[TABLE]

By the independence of $y_{i}^{*}$ and $\eta_{i}$, since $\hat{y}_{i}=y_{i}+\eta_{i}$, the distribution of $\hat{y}_{i}^{*}$ is given by

[TABLE]

Equation (16) follows from (18) and the distributivity of the convolution operator.

To show (17), let $i_{1},i_{2}\in[m]$ be such that $i_{1}\neq i_{2}$ and let

[TABLE]

Similarly to the previous case, we compute

[TABLE]

Note that if $(A_{i_{1},j}-A_{i_{2},j})x_{j}\neq 0$, then $j\in\operatorname{supp}(x)$ and $(A_{i_{1},j}-A_{i_{2},j})$ is either $+1$ with probability $\frac{1}{2}$ or $-1$ with probability $\frac{1}{2}$. Then, Hence,

[TABLE]

then,

[TABLE]

Letting $b^{\prime}_{j}\sim\mathrm{Ber}\left(\frac{2d\rho}{n}(1-o(1))\right)$,

[TABLE]

Again, invoking Lemma III.2 with $p_{n}={2d\rho}(1-o(1))$ we obtain that $p=2d\rho$ and also that as $n\rightarrow\infty$, $y_{i_{1}}-y_{i_{2}}\to g^{*}$ with

[TABLE]

Therefore, Given that $\eta_{i_{1}}-\eta_{i_{2}}\sim\nu\ast\nu^{-}$ and that

[TABLE]

we convolve (19) with $\nu\ast\nu^{-}$ to recover (17). ∎

We are now ready to prove Theorem I.4,

Proof. Theorem I.4.

Using Bayes rule we write

[TABLE]

From (16) we can deduce that

[TABLE]

and using equation (16) from Theorem III.3 we obtain that as $n\to\infty$,

[TABLE]

Coupling (20), (22) and (21) yields (5).

Again, by Bayes rule,

[TABLE]

Noting that $\hat{y}_{i_{1}}-\hat{y}_{i_{2}}=\left(y_{i_{1}}-y_{i_{2}}\right)+\left(\eta_{i_{1}}-\eta_{i_{2}}\right)$,

[TABLE]

By (16) from Theorem III.3 we obtain that as $n\to\infty$,

[TABLE]

Coupling (23), (24), and (25) yields (6).

∎

III-A Explicit formulas for the centred Gaussian case

We further elucidate (5)-(6) from Theorem I.4 in the case when $\mu$ and $\nu$ are Gaussian with mean zero, which are the distributions considered in Section IV. To this end, let $\mu=\mathcal{N}(0,\sigma_{s}^{2})$ and $\nu=\mathcal{N}(0,\sigma_{n}^{2})$. We observe that in this case $\bar{\mu}=\mu$, $\mu_{q}=\bar{\mu}_{q}=\mathcal{N}(0,q\sigma_{s}^{2})$ and $\tilde{\nu}=\mathcal{N}(0,2\sigma_{n}^{2})$. Hence, denoting by $\varphi(\cdot\mid\sigma^{2})$ the probability density function of a centred Gaussian random variable with variance $\sigma^{2}$, we obtain

[TABLE]

III-A1 Estimating the tails in the Gaussian case

We approximate the infinite sum in in the denominators of (26) and (27) by doing an approximation to the tail of this summation.

Lemma III.4** (Sums of centred density functions).**

Let $\mu_{i}$ be the density function of a random varible with mean zero and variance $\sigma_{i}^{2}$ and let $\alpha_{i}>0$ be such that $\sum_{i}\alpha_{i}=1$. Then, $\mu=\sum_{i}\alpha_{i}\mu_{i}$ is the density function of a random variable with mean zero and variance $\sum_{i}\alpha_{i}\sigma_{i}^{2}$.

Proof.

Let $x\sim\mu$ and $x_{i}\sim\mu_{i}$ be such that $\mathbb{E}\left[{x_{i}}\right]=0$ and $\mathbb{V}ar\left[{x_{i}}\right]=\sigma_{i}^{2}$.

[TABLE]

∎

To simplify notation let $\sigma_{q}^{2}=q\sigma_{s}^{2}+\sigma_{n}$ for $q\in\mathbb{N}\cup\{0\}$. Consider the series in the denominator of (22)

[TABLE]

Let,

[TABLE]

and write

[TABLE]

Note that $S_{b}/R_{z}$ satisfies the conditions of Lemma III.4 so it corresponds to the density function of a centred random variable with variance

[TABLE]

Therefore,

[TABLE]

A similar argument shows that

[TABLE]

Where $\sigma_{q,e}^{2}=q\sigma_{s}^{2}+\sigma_{n}^{2}$, and

[TABLE]

III-B Comparison with empirical probabilities

We test the approximations given in (28) and (29) by randomly generating $\hat{y}$ and $y$ according the generating models $\operatorname{GM}(n,\delta n,\rho\delta n,7,\mathcal{N}(0,1),\mathcal{N}(0,\sigma^{2}))$, for $\delta=0.3$, $\rho\in\{0.1,0.3\}$ and $\sigma\in\{10^{-3},10^{-2}\}$. The results can be seen in Figure 2.

Overall the analytical expressions fit the empirical probabilities very well, indicating that the approximations we made in the calculations above are justified. However, we observe that as $\rho$ and especially $\sigma$ increase, both functions drop significantly. This means that for large values of these parameters, the noise eventually dominates and it is difficult to decided whether values are zero or equal.

III-C Scaled probabilities $\breve{p}_{z}$ and $\breve{p}_{e}$

We mentioned in Section I that our algorithms do not implement the functions $p_{z}$ and $p_{e}$ exactly, but a scaled version of these. As observed in Figure 2, the value of $\max_{s}p_{e}(s)$ and $\max_{s}p_{z}(s)$ varies significantly as $\sigma$ and $\rho$ change. Algorithm 1 evaluates whether a given score is large or not by implementing a sweeping parameter $t$ that is set to one at the beginning of the algorithm and decreased by a constant $c$ after every iteration. In order to use a fixed initial $t$ we consider the scaled probabilities $\breve{p}_{e}$ and $\breve{p}_{z}$ in $\eqref{eq:breve}$; otherwise the initial value of $t$ would depend on $\sigma$ and $\rho$.

III-D Adaptive $k$

Algorithm 1 can optionally account for the sparsity of the current estimate via the flag adaptive_k. In the noiseless model of $r=A\hat{x}$, an update of the form $\hat{x}\leftarrow\hat{x}+\omega e_{j}$ with Parallel-$\ell_{0}$ guarantees that $j\in\operatorname{supp}(x)$ so that at the next iteration the problem with residual $r-a_{j}\omega$ and $(k-1)$-sparse signal is considered. The adaptive_k flag updates the sparsity prior in $\hat{x}$ after every update in hope of having more reliable estimates of $p_{e}$ and $p_{z}$. We will see in the numerical experiments that under the data generating model that we tested this strategy does not bring substantial benefits. We don’t rule out the posssibility that there are other signal and noise distributions for which this flag becomes especially useful, but we leave that as future work.

IV Numerical experiments

In this section we present numerical experiments which validate the efficacy of Robust-$\ell_{0}$ decoding. In particular, we contrast Robust-$\ell_{0}$ with other state-of-the-art greedy algorithms for compressed sensing in terms of their ability to recover the measured signal for varying problem sizes $(k,m,n)$ as well as their computational complexity. To facilitate reproducibility we begin by describing the stopping conditions and measures used to denote successful recovery in the presence of noise in Section IV-A, along with how the parameter $c$ is varied in Section IV-B. We then present the main numerical results in Section IV-C where the algorithms phase transitions and runtime are presented, along with Sections IV-D and IV-E which show further details on Robust-$\ell_{0}$ decoding’s performance as a function of noise variance and for extreme subsampling respectively.

IV-A Stopping conditions

We are interested in the signal model $y=Ax+\eta$ where $\eta\sim\mathcal{N}(0,\sigma^{2}I_{m\times m})$. If $\hat{x}$ is an approximation to $x$ the residual is $r=y-A\hat{x}$. Note that if $\hat{x}=x$, then

[TABLE]

and we should not seek reductions in the residual below $\|\eta\|_{1}$ which would result in fitting to the additive noise. We further account for the variance of $\|\eta\|_{1}$ and denote the algorithm to have successfully recovered $x$ if $\hat{x}$ satisfies

[TABLE]

for some $C_{1}\geq 0$. We should be aware that the right hand side of (30) might be greater than 1 for some choices of $k$, $m$ and $\sigma$. When this happens, the stopping condition (30) becomes invalid since we expect $\hat{x}$ to have captured a proportion of the $\ell_{1}$-energy of $\|x\|_{1}$. Hence, if the right hand side of (30) is greater than $\frac{1}{10}$ we clip the upper bound at this value and use the stopping condition

[TABLE]

For the numerical experiments conducted in this section we consider nonzero entries in $x$ drawn as $x_{i}\sim\mathcal{N}(0,1)$, and noise $\eta_{i}\sim\mathcal{N}(0,\sigma^{2})$ for which $\mathbb{E}\left[\|x\|_{1}\right]=k\sqrt{\frac{2}{\pi}}$ and $\mathbb{E}\left[\|\eta\|_{1}\right]=m\sigma\sqrt{\frac{2}{\pi}}$ and moreover $\mathrm{Var}\left[\|\eta\|_{1}\right]=m\sigma^{2}\left(1-\frac{2}{\pi}\right)$, see e.g. [24].

IV-B Selection of parameter $c$ in Algorithm 1

The sweeping parameter $t$ in Algorithm 1 is initialised at 1 and updated by decreasing it by a constant value $c$. We observed in our experiments that the quality of the phase transitions are sensitive on the parameter $c$ especially for low $\delta$ and $\rho$. We do not provide a way to fine-tune $c$, but we run our phase transitions with the following choices:

If the algorithm is quantised,

[TABLE] 2. 2.

If the algorithm is continuous,

[TABLE]

The values in (32) and (33) were chosen heuristically for $\nu$ and $\mu$ Gaussian.

IV-C Phase transitions and runtime

We benchmark the variants of Robust-$\ell_{0}$ against other greedy algorithms via their phase-transitions and runtime. The user can supply two binary flags, adaptive_k and quantised which yield four different variants of Robust-$\ell_{0}$. We assigned a unique label to each of these variants as described in Table II.

The phase transition of a compressed-sensing algorithm [25] is the largest value of $k/m$ for which the algorithm is typically able recovery all $k$ sparse vectors with sparsity less than $k$ for a fixed $m/n$. Hence, for a fixed value of $\delta=m/n$ the phase transition of an algorithm is the largest value $\rho^{*}(\delta)$ for which the algorithm converges for all $\rho(\delta)<\rho^{*}(\delta)$. The value $\rho^{*}(m/n)$ often converges to a fixed value as $n\rightarrow\infty$, so phase transitions often partition the $\delta\times\rho$ space into two regions: One in which the algorithm converges with high probability and another in which the algorithm doesn’t converge with high probability. We benchmark Robust-$\ell_{0}$ against the algorithms presented in [18, 19, 26]. Specifically, our tests include the following algorithms,

{Robust-$\ell_{0}$, Robust-$\ell_{0}$-adaptive, Robust-$\ell_{0}$-quantised, Robust-$\ell_{0}$-adaptive-quantised, SSMP, SMP, CGIHT}.

In the deterministic case Parallel-$\ell_{0}$ was compared against a range of combinatorial compressed sensing algorithms in [1]; out of those we have selected SSMP and SMP as these perform best and are similar in nature to Robust-$\ell_{0}$. We also compare with CGIHT, as this algorithm was shown to be the fastest among the greedy algorithms compared in [13, 26]. Figures 3a, 3e and 3i show the phase transition curves for these algorithms with $\sigma=10^{-3},\;10^{-2},\;10^{-1}$ respectively. The curves were computed by setting $n=2^{18}$, $d=7$, and using the stopping condition $\|r\|_{1}\leq\mathbb{E}\left[\|\eta\|_{1}\right]=m\sigma\sqrt{\frac{2}{\pi}}$ and a success condition (30) with $C_{1}=1$. The testing is done at $m=\delta_{p}n$ for

[TABLE]

For each $\delta_{p}$, we set $\rho=0.01$ and generate 10 synthetic problems to apply the algorithms to, with a problem generated as in $\operatorname{GM}$ with the given parameters and $\mu$ and $\nu$ being normal Gaussian. If at least one such problem was recovered successfully, the sparsity ratio $\rho$ is increased by $0.01$ and the experiment is repeated. Following the testing framework in [27], the recovery data is fitted using a logistic function and finally the 50% recovery transition function is computed and presented in the phase transition plots contained herein. Figures 3b, 3f and 3j show a selection map for these algorithms. Namely, these plots indicate which algorithm requires the least computational time111All the numerical results presented here were performed using a Linux machine with Intel Xeon E5-2643 CPUs 3.30 GHz, NVIDIA Tesla K10 GPUs and executed from Matlab R2016b. The code was added to the GAGA library available at http://www.gaga4cs.org/, and described in [28], so as to facilitate large scale benchmarking against the other algorithms presented here which are also contained in the aforementioned library. at each point in the $\delta\times\rho$ space where the algorithm converges. Finally, Figures 3c, 3g and 3k show the total time for convergence in milliseconds for the fastest algorithm at each point in the $\delta\times\rho$ space. The ratio of the time for the second fastest algorithm over the time for the fastest algorithm is given in Figures 3d, 3h, and 3l.

We can see from Figure 3 that CGIHT [26] dominates the upper region of the phase transition space, while the Robust-$\ell_{0}$ algorithms only converge for $\rho\lessapprox 0.3$ which is consistent with the observed phase transitions for Parallel-$\ell_{0}$ [1]. In terms of speed, Robust-$\ell_{0}$ seems to be the most competitive for $\sigma\in\{10^{-3},10^{-2}\}$ and $\rho\lessapprox 0.2$. However, for large noise, $\sigma=10^{-1}$, CGIHT becomes the fastest algorithm of all. We remark that while CGIHT performs very well in our numerical tests, the current theory developed for it does not hold in the setting considered here, as it requires zero-mean columns in $A$.

Figure 4 shows the widths for the Robust-$\ell_{0}$ algorithms. The widths measure how sharp the phase transition of an algorithm is; namely, how thin the boundary between the region of recovery with high-probability and the region of recovery where combinatorial search is needed. It has been shown that the widths of a compressed sensing algorithm tend to zero as $n\rightarrow\infty$ when decoding with linear programming [29], and we usually expect the same behaviour for other algorithms [26]. Figure 4 show that the widths for the Robust-$\ell_{0}$ algorithms indeed decrease with $n$ and with $\delta$. The observed smoothness of the phase transition widths signal also suggest that the stopping conditions of the algorithm are consistent for the problem under consideration.

IV-D Dependence on noise variance, $\sigma$

We now investigate the extent to which the phase transitions of the algorithm decrease as we increase the noise level $\sigma$. To do this, we consider $\delta\in\{0.01,0.02,0.1,0.2\}$ and for each value of $\delta$ we compute we define the grid

[TABLE]

Then at each value of $\sigma$ we let $\rho=0.01$ and draw ten problem instances from $\operatorname{GM}$ with signal distribution $\mathcal{N}(0,1)$ and noise distribution $\mathcal{N}(0,\sigma^{2})$. If at least one of the problems was recovered successfully, then we set $\rho\leftarrow\rho+0.01$ and repeat the experiment. We do this procedure for each of the algorithms considered and record the largest $\rho$ having at least $50\%$ success rate. The results are shown in Figure 5. In order to show where the clipping in (31) becomes active, the figures also show a TOL-curve which partitions the space into the region where the right hand side of (31) equals $\frac{1}{10}$ (bottom region) and the region where it equals the right hand side of (30) (top region). We can appreciate from Figure 5 that for small $\delta$ both Robust-$\ell_{0}$ and SSMP have the best recovery capabilities, with Robust-$\ell_{0}$ being preferable for noise levels $\sigma\lessapprox 10^{-2}$ and SSMP being better suited for larger noise levels. For larger $\delta$, CGIHT is preferable except for very low noise levels.

IV-E Phase transitions for extreme subsampling, $\delta\ll 1$

The numerical experiments of Parallel-$\ell_{0}$ in [1] showed flat phase transitions; that is, it was observed that $\rho^{*}(\delta)$ remained approximately $0.3$ as $\delta\rightarrow 0$ provided $n$ was sufficiently large. While Robust-$\ell_{0}$ does not exhibit precisely the same behaviour in the presence of noise, we do observe that $\rho^{*}(\delta)$ remains nontrivial even for $\delta$ as small as $10^{-3}$, again provided $n$ is sufficiently large. We provide numerical evidence for this in Figure 6. For each $(\delta,\sigma)\in\{0.001,0.01\}\times\{0.001,0.01\}$ we let $\rho=0.01$ and solve ten problems drawn from $\operatorname{GM}$ with nonzero distribution $\mathcal{N}(0,1)$ and noise distribution $\mathcal{N}(0,\sigma^{2})$. If at least one problem instance converges, we average the run-time of the problems that converged and repeat the process with $\rho\leftarrow\rho+0.01$. We plot the timing at each $\rho$ for parameter values for which at which at least $50\%$ of the problems were successfully recovered under the criteria (31).

It can be seen in Figure 6a-6d that the phase transition either remains nearly unchanged or increases as $n$ increases from $2^{22}$ to $2^{24}$. In particular, Figure 6a shows that for $\sigma=0.001$ and $\delta=0.001$, the phase transitions for all variants of the Robust-$\ell_{0}$ algorithms in fact increase from $\rho\approx 0.08$ to $\rho\approx 0.1$ as $n$ increases from $2^{22}$ to $2^{24}$. Additionally, contrasting Figures 6a and 6b or 6c and 6d shows that a ten-fold increase in $\sigma$ has the expected effect of reducing the phase transition and increasing the computational time. Figures 6a and 6b show the phase transition for Robust-$\ell_{0}$ and Robust-$\ell_{0}$-adaptive remain significant even for $\delta$ as small as $10^{-3}$. Figures 6c and 6d show results for the same set of experiments, but for $\delta=0.01$ which corresponds to a ten-fold increase in $\delta$ over the value used in Figures 6a and 6b. For $\sigma=0.001$ the phase transition of Robust-$\ell_{0}$ reaches $\rho\approx 0.17$, while for $\sigma=0.01$ the phase transitions drop to $\rho\approx 0.12$ and there is an increase in the computational time.

V Conclusions

We have shown that the decoding framework presented in [1] can be extended to the case where the measurements are corrupted by additive noise. This framework is extended by deriving the posterior distribution of an entry in the residual being zero or being equal to another residual entry given the corrupted measurements. This Bayesian approach to decoding was implemented in Robust-$\ell_{0}$ and its four variants. We show that the resulting algorithms inherits some desirable properties from Parallel-$\ell_{0}$ like high phase transitions for low $\delta$ and large $n$ and low-latency. However, these qualities are weakened by the corruption of the measurements. Our numerical experiments show that Robust-$\ell_{0}$ should be considered in cases of moderate noise and $\rho\lessapprox 0.3$.

Acknowledgments

The authors would like to thank William Carson for many helpful discussions on imaging applications of this work. This work was supported by The Alan Turing Institute under the EPSRC grant EP/N510129/1.