Compressed sensing reconstruction using Expectation Propagation

TL;DR

This paper introduces a Bayesian inference approach using Expectation Propagation for compressed sensing, demonstrating superior accuracy over existing methods, especially with correlated measurement matrices.

Contribution

The paper presents a novel application of Expectation Propagation to compressed sensing within a Bayesian framework, outperforming existing algorithms in accuracy.

Findings

EP-based method outperforms state-of-the-art algorithms

Method is more accurate with correlated measurement matrices

Bayesian approach effectively enforces sparsity constraints

Abstract

Many interesting problems in fields ranging from telecommunications to computational biology can be formalized in terms of large underdetermined systems of linear equations with additional constraints or regularizers. One of the most studied ones, the Compressed Sensing problem (CS), consists in finding the solution with the smallest number of non-zero components of a given system of linear equations $\boldsymbol y = \mathbf{F} \boldsymbol{w}$ for known measurement vector $\boldsymbol{y}$ and sensing matrix $\mathbf{F}$. Here, we will address the compressed sensing problem within a Bayesian inference framework where the sparsity constraint is remapped into a singular prior distribution (called Spike-and-Slab or Bernoulli-Gauss). Solution to the problem is attempted through the computation of marginal distributions via Expectation Propagation (EP), an iterative computational scheme originally developed in Statistical Physics. We will show that this strategy is comparatively more accurate than the alternatives in solving instances of CS generated from statistically correlated measurement matrices. For computational strategies based on the Bayesian framework such as variants of Belief Propagation, this is to be expected, as they implicitly rely on the hypothesis of statistical independence among the entries of the sensing matrix. Perhaps surprisingly, the method outperforms uniformly also all the other state-of-the-art methods in our tests.

Equations

$$\mathbf{y}=\mathbf{F}\mathbf{w},$$

$$E(\mathbf{w})=\frac{1}{2}{\Vert \mathbf{y}-\mathbf{F}\mathbf{w}\Vert }_2:=\frac{1}{2}{\left(\mathbf{y}-\mathbf{F}\mathbf{w}\right)}^T\cdot \left(\mathbf{y}-\mathbf{F}\mathbf{w}\right),$$

$$P(\mathbf{y}|\mathbf{F},\mathbf{w})$$

$$P\left(\mathbf{w}|\mathbf{F},\mathbf{y}\right)$$

$$\psi (w_i)=(1-\rho )\delta (w_i)+\frac{\rho }{\sqrt{2\pi \lambda }}{e}^{-\frac{w_i^2}{2\lambda }},$$

$$Q(\mathbf{w}|\mathbf{F},\mathbf{y})$$

$$Z_Q={\left(2\pi \right)}^{\frac{N}{2}}{\left(\det \mathbf{\Sigma }\right)}^{\frac{1}{2}},$$

$${\mathbf{\Sigma }}^{-1}:=\beta {\mathbf{F}}^T\mathbf{F}+\mathbf{D},\bar{\mathbf{w}}:=\mathbf{\Sigma }(\beta {\mathbf{F}}^T\mathbf{y}+\mathbf{D}\mathbf{a}),$$

$$Q^{(n)}(\mathbf{w}|\mathbf{F},\mathbf{y})$$

$${\mathbf{\Sigma }}_{(n)}^{-1}=\beta {\mathbf{F}}^T\mathbf{F}+{\mathbf{D}}_{(n)},{\bar{\mathbf{w}}}_{(n)}={\mathbf{\Sigma }}_{(n)}\left(\beta {\mathbf{F}}^T\mathbf{y}+{\mathbf{D}}_{(n)}\mathbf{a}\right),$$

$${⟨w_n⟩}_{Q^{\left(n\right)}}$$

$$⟨w_n{⟩}_Q$$

$$⟨w_n^2{⟩}_Q$$

$${\left({\bar{w}}_{(n)}\right)}_n$$

$$d_n$$

$$a_n$$

$${ϵ}_t=\underset{n}{\max }\left|{⟨w_n⟩}_{Q_t^{\left(n\right)}}-{⟨w_n⟩}_{Q_{t-1}^{\left(n\right)}}\right|+\left|{⟨w_n^2⟩}_{Q_t^{\left(n\right)}}-{⟨w_n^2⟩}_{Q_{t-1}^{\left(n\right)}}\right|,$$

$$r=\frac{{\sum }_{k=1}^N(w_k-w_{sm})({\hat{w}}_k-{\hat{w}}_{sm})}{\sqrt{{\sum }_{k=1}^N(w_k-w_{sm}{)}^2}\sqrt{{\sum }_{k=1}^N({\hat{w}}_k-{\hat{w}}_{sm}{)}^2}},$$

$$MSE=\frac{1}{N}\sum _{k=1}^N(w_k-{\hat{w}}_k{)}^2.$$

$${\alpha }_1({\rho }_0)=2(1-{\rho }_0)H({\hat{\chi }}^{-1/2})+{\rho }_0.$$

$$\hat{\chi }={\alpha }^{-1}\left[2(1-{\rho }_0)\left((\hat{\chi }+1)H({\hat{\chi }}^{-1/2})-{\hat{\chi }}^{1/2}\frac{e^{-1/(2\hat{\chi })}}{\sqrt{2\pi }}\right)+{\rho }_0(\hat{\chi }+1)\right],$$

$$MSE=\rho MS{E}_1+(1-\rho )MS{E}_2,$$

$$\mathbf{F}=({\mathbf{f}}_1,...,{\mathbf{f}}_M{)}^T,$$

$${\mathbf{f}}_i\sim \mathcal{N}(0,\mathbf{S}),i=1,...,M.$$

$$\mathbf{S}={\mathbf{Y}}^T\mathbf{Y}+\mathbf{\Delta },$$

$$\mathbf{F}^{\prime}=\left[\begin{array}[]{ccc}\mathbf{I}&|&\mathbf{G}\end{array}\right]=\left[\begin{array}[]{cccccccc}1&&&&G_{1,1}&\dots&G_{1,N-M}\\ &1&&&G_{2,1}&\dots&G_{2,N-M}\\ &&\dots&&\dots&\dots&\dots\\ &&&1&G_{M,1}&\dots&G_{M,N-M}\end{array}\right]$$

$$\mathbf{w}=:({\mathbf{w}}^{(d)},{\mathbf{w}}^{(i)})=w_1^{(d)},\dots ,w_M^{(d)},w_1^{(i)},\dots ,w_{N-M}^{(i)}$$

$${\mathbf{w}}^{(d)}+\mathbf{G}{\mathbf{w}}^{(i)}={\mathbf{y}}^{\prime }$$

$${\stackrel{ˉ}{}\mathbf{w}}^{(d)}=-\mathbf{G}{\bar{\mathbf{w}}}^{(i)}+{\mathbf{y}}^{\prime }$$

$${\mathbf{\Sigma }}_{{\mathbf{w}}^{(i)}}={\left({\mathbf{D}}_{{\mathbf{w}}^{(i)}}+{\mathbf{G}}^T{\mathbf{D}}_{{\mathbf{w}}^{(d)}}\mathbf{G}\right)}^{-1},{\mathbf{\Sigma }}_{{\mathbf{w}}^{(d)}}=\mathbf{G}{\mathbf{\Sigma }}_{{\mathbf{w}}^{(i)}}{\mathbf{G}}^T$$

$$\psi (w_n)=(1-\rho )\delta (w_n)+\frac{\rho }{\sqrt{2\pi \lambda }}{e}^{-\frac{w_n^2}{2\lambda }},$$

$$Q^{(n)}(w_n)=\frac{1}{Z_{Q^{(n)}}}{\tilde{Q}}_n(w_n){\psi }_n(w_n)$$

Full text

1 Introduction

The problem of Compressed Sensing (CS) [1, 2, 3] has led to significant developments in the field of sparse approximation and representation [4], together with the more traditional framework of optimal signal processing [4, 5].

In general, CS deals with the reconstruction of a sparse $N-$dimensional signal from $M$, often noisy, measurements. In this context, a sparse signal is characterized by a large (possibly the largest) number of components equal to zero. In practical cases, the regime of interest is given by the $M\ll N$, $K\ll N$ limit of the problem, where $K$ is the number of non-zero components. A brute-force exhaustive search of the $K-$dimensional minimal support of the $N-$dimensional signal would lead to an exponential explosion in the search space as there are $N\choose K$ possible base supports to be explored. In the noiseless limit of CS, even if the $K$-dimensional minimal base were given, at least $M=K$ measurements are necessary to uniquely identify the solution.

CS can be easily formulated: let $\bm{w}\in{\mathbb{R}}^{N}$ be compressed into a vector $\bm{y}\in{\mathbb{R}}^{M}$ (from now on we will assume $M<N$) through a linear transformation:

[TABLE]

where $\mathbf{F}\in{\mathbb{R}}^{M\times N}$ is a linear operator of maximal rank often referred to as the measurement or sensing matrix. The problem is to determine the vector $\bm{w}$ from the knowledge of the measurement matrix $\mathbf{F}$ and of the compressed vector $\bm{y}$.

The relevance of CS in statistical physics stems from the link with the statistical mechanics of disordered systems, in analogy with many other combinatorial optimization problems [6, 7, 8]. In particular CS is an optimization problem with quenched disorder (the measurement matrix) and, as such, it is amenable to analytic treatment using replica theory [9] as developed in the field of spin-glasses [10].

From equation (1), it is clear that, as long as $M<N$, there are infinitely many solutions to the system of equations. However, if one imposes the supplementary condition on the sparsity of the signal $\bm{w}$, solving the problem may still lead to a unique retrieved signal. Thus, a problem in CS is to determine the set of conditions under which it is possible to find a solution of equation (1) that is as sparse as possible. A practical way to enforce sparsity is through the minimization of the $L_{p}$-norm of the $\bm{w}$ vector in the space of solutions which is defined in equation (1). A particularly successful line of research has been pursued through the minimization of the $L_{1}$ norm [1, 11, 3, 12, 2, 13, 9], for which the following prediction was obtained in the case of independent and identically distributed (i.i.d.) random Gaussian measurement matrices: in the large $N$ limit, there is a non-trivial region $\rho=K/N$, $\alpha=M/N$ where an exact reconstruction of the original signal is indeed possible. The parameters $\rho$ and $\alpha$ that characterize the signal are called the signal density and the measurement rate, respectively.

However, in many applications, the sensing matrix may not be random [14] or the i.i.d. assumption might not hold and, in these cases, several algorithms could fail to decode the signal. Examples of deterministic matrices include chirp sensing matrices [15], which have been applied to image reconstruction [16], and second-order Reed-Muller matrices [17], whereas examples of correlated random matrices include random partial Fourier matrices [18, 19], which are encountered in MRI [20] as well as in other applications [21, 22], and partial random circulant and Toeplitz matrices [23], which arise in the presence of convolutions.

Our approach here, is to treat CS as a Bayesian inference problem, i.e. we focus our attention on determining the probability of observing, a posteriori, the signal $\bm{w}$ when we have observed a set of measurements $\bm{y}$ and given the sensing matrix $\mathbf{F}$. While equation (1) is easily cast in a Gaussian-like likelihood function, the prior over the variables $\bm{w}$, enforcing the sparsity constraint, plays a key role on the goodness of the solution and on the complexity of the problem. A prior that corresponds to the minimization of a $L_{p}$ norm, with $p={1,2}$, allows an easy marginalization of the corresponding posterior probability with the drawback of softening the sparsity requirement. Similarly to [9, 24, 25], we consider the so-called spike-and-slab prior [26] which is understood as the minimization of the $L_{0}$ norm of the vector $\bm{w}$. The posterior distribution results to be intractable in practice and thus some approximations need to be sought.

Almost twenty years ago Expectation Propagation (EP), an iterative scheme to approximate intractable distributions, has been introduced first in the field of statistical physics [27, 28] and shortly after in the field of theoretical computer science [29]. Recently, EP inspired inference strategies – similar to the one we present here – have been proposed to solve other underdetermined linear constraint problems such as the problem of sampling solutions from the reconstruction of large scale metabolic networks [30] and of tomographic images [31]. Here we propose an efficient and accurate EP-based reconstruction strategy for CS which, moreover, does not require i.i.d. measurement matrices.

Other attempts to solve sparse linear models using the EP approximation can be found in [32] – where the authors use a similar EP implementation to the one we adopt here and make use of a Laplace prior – and in [33] – where, in addition to the spike-and-slab, a Bernoulli prior on the components of the signal is introduced and the EP update scheme involves only three approximating factors of the original posterior distribution.

All these approaches have the same computational complexity, which is dominated by a single matrix inversion per EP cycle. In our implementation, we integrate the inference of the signal with a variational learning of the density parameter of the spike-and-slab prior. Besides, we show in A that the original EP scheme is equivalent to an alternative formulation of the EP update equations that takes into account the linear constraints (1) exactly. This implementation has the advantage of reducing the size of the matrix to be inverted, thus decreasing the computing time. Throughout the work presented in this paper, we have numerically checked that both procedures can be used interchangeably in the case of noiseless measurements, as they lead to the same results.

In conclusion, we believe that our EP algorithm shows some original features that have not been proposed by other methods, namely: (i) the quality of EP reconstruction on correlated measurement matrices is the same as in the case of uncorrelated measurements whereas all other methods we tested fail; (ii) to reconstruct signals in the noiseless case, we introduced an EP formulation performing analytically the zero temperature limit; (iii) we learn the density parameter of the spike and slab prior from the data by minimizing iteratively the EP free energy.

The outline is as follows: we introduce the EP approach to CS and set the basic notations in section 2, we show a thorough comparison with other state of the art algorithms in section 3 and, finally, we present our conclusions in section 4.

2 Expectation Propagation

In this section we introduce the so-called finite temperature formulation of EP, where we allow the measurement vector $\bm{y}$ to be noisy. We consider the undersampling regime $M<N$, although there is no technical limitation in considering the $M\geq N$ case. The linear constraints in equation (1) can be alternatively mapped into a minimization problem of the following quadratic form:

[TABLE]

From a Bayesian perspective, the probability of observing the vector ${\bm{y}}$ given the matrix $\mathbf{F}$ and the vector ${\bm{w}}$ is:

[TABLE]

where we introduced a fictitious inverse temperature $\beta$ that we can take as large as we wish in order to enforce the linear bounds expressed by equation (1). Alternatively we can interpret equation (3) as the probability of observing an additive noise vector $\bm{y}-\mathbf{F}\bm{w}$ whose elements are distributed according to a Gaussian density of zero mean and variance $\beta^{-1}$. The posterior distribution for the vector ${\bm{w}}$ reads:

[TABLE]

where in the last step we restricted the structure of the prior $P({\bm{w}}):=\prod_{i=1}^{N}\psi_{i}(w_{i})$ to a factorized form, although more general structures can be considered, e.g. as in [34]. In this work, we have considered the so-called spike-and-slab prior [26]:

[TABLE]

where $\delta$ denotes the Dirac delta function, in order to model any prior knowledge about the sparsity of the signal.

We seek a solution vector $\hat{\bm{w}}$ whose components are the first moments of the marginal densities of equation (4). Contrarily to the maximum a posteriori estimate, it can be proven that this strategy minimizes the mean squared error between the true and the recovered signal. Unfortunately, due to the non-convex nature of the spike-and-slab prior, there exists no technique able to perform the marginalization of the posterior probability in equation (4). In the following we introduce the EP approximation scheme which relies on an adaptive Gaussian approximation of the marginal probabilities of interest.

2.1 The approximate posterior distribution

EP [29] is an efficient approximation to compute posterior probabilities. EP was first introduced as an improved mean-field method [27, 28] and further developed in the framework of Bayesian inference problems in the seminal work of Minka [29]. The approximated distribution consists in substituting the typically analytically intractable $\psi_{i}$ priors with univariate Gaussian distributions $\phi_{i}(w_{i})={\cal N}(w_{i};a_{i},d_{i})$ of mean $a_{i}$ and variance $d_{i}.$ The approximated posterior thus reads:

[TABLE]

where:

[TABLE]

[TABLE]

and $\bm{D}$ is a diagonal matrix having diagonal elements $d_{1}^{-1},\dots,d_{N}^{-1}$. We now need a way to fix the parameters ${\bm{a},\bm{d}}$ of the prior. To do so, we focus on the $n^{th}$ variable $w_{n}$ (with $1\leq n\leq N$), and in particular on its approximated prior $\phi_{n}$. We can define the tilted distribution $Q^{(n)}$ as:

[TABLE]

where:

[TABLE]

and, in analogy with equation (8), $\bm{D}_{(n)}$ is a diagonal matrix of elements $d_{m}^{-1}$ for all diagonal elements $m\neq n$ and zero for $m=n$. The *tilted *distribution differs from the approximated posterior $Q$ as it contains the original prior $\psi_{n}$ instead of the $n$-th approximated prior $\phi_{n}$. Intuitively, we expect that the ${n}$-th tilted distribution provides a better approximation of the expectation values related to the $n$-th variable than the multivariate Gaussian approximation. From a computational point of view, the presence of a single intractable prior in the tilted distribution does not prejudice the efficiency of the algorithm.

A natural way of fixing the optimal $\bm{a},\bm{d}$ parameters consists in requiring the approximated $Q$ distribution to be as similar as possible to the tilted distribution $Q^{(n)}.$ To do so, we minimize the Kullback-Leibler divergence $D_{KL}\left(Q^{(n)}||Q\right)$. Perhaps not surprisingly, this procedure is equivalent to equating the first two moments of the two distributions:

[TABLE]

where $\langle\cdot\rangle_{Q}$ and $\langle\cdot\rangle_{Q^{(n)}}$ denote averages w.r.t. $Q(\bm{w})$ and $Q^{(n)}(\bm{w})$, respectively.

Notice that the computation of the moments of the tilted distribution on the left-hand side of equation (11) depends on the functional form of the prior considered. We refer to B for the expression of the moments of the tilted distributions used in the case of a spike-and-slab prior.

Thanks to the multivariate Gaussian form of the approximated distribution, it is a simple exercise to compute the moments of $Q$:

[TABLE]

From equations (8) and (10), it is clear that the two matrices $\Sigma^{-1}$, and $\Sigma_{(n)}^{-1}$ differ only in a diagonal term. We can thus exploit a low-rank update property to relate the two inverses. It turns out that the tilted parameters are related to the approximated ones:

[TABLE]

Upon imposing the moment matching condition (11), we eventually get an explicit equation for the prior parameters $a_{n},d_{n}$:

[TABLE]

The ${\bm{a}},{\bm{d}}$ parameters are sequentially updated until a fixed point is eventually reached: numerically, we need to set a threshold below which the algorithm stops. To this purpose, for each iteration $t$ (i.e. for each update of the ${\bm{a}},{\bm{d}}$ vectors), we can define an error $\epsilon_{t}$ as:

[TABLE]

where $Q_{t}^{\left(n\right)}$ is the tilted distribution with parameters computed at iteration $t$. In practice, the algorithm stops when $\epsilon_{t}<10^{-6}$. At convergence, the tilted distributions provide an approximation to the marginal densities of the posterior in equation (4) and their first moments $\left\langle w_{n}\right\rangle_{Q^{\left(n\right)}}$ provide the estimate $\hat{\bm{w}}$ of the unknown vector $\bm{w}$.

For the sake of convenience, the EP procedure with low rank update that we have just presented is summarized in Algorithm 1.

3 Experimental results with synthetic data

In this section, we present our empirical results obtained by means of numerical simulations. We will first consider compressed sensing with i.i.d. random sensing matrices and then we will discuss some results related to the case of correlated random matrices.

3.1 Uncorrelated measurements

We consider $M\times N$ measurement matrices having i.i.d. entries sampled from a standard normal distribution $\mathcal{N}(0,1)$. The signal vector has $K=\rho N$ nonzero components, which are also sampled from a Gaussian distribution with zero mean and unit variance. The measurements are assumed to be noiseless. Note that, in general, the parameter $\rho$ in equation (5) is unknown and needs to be estimated. We show in C how one can infer $\rho$ within the framework provided by EP.

We have run EP throughout the $\rho$-$\alpha$ plane. The parameters used in our EP simulations are $\lambda=1$ and, in the finite temperature formulation, $\beta=10^{9}$.

In order to measure the quality of the reconstruction, we consider the sample Pearson correlation coefficient $r$ of the true vector and of the reconstructed vector:

[TABLE]

where $w_{sm}$ and $\hat{w}_{sm}$ are the sample means of the signal and of the inferred vector, respectively. We also consider the within-sample mean squared error as a measure of the reconstruction error:

[TABLE]

We plot the sample correlation coefficient in figure 1 for a single simulation at each given $\rho$ and $\alpha$ and progressively larger values of $N$, namely $N=200$ (figure 1(a)), $N=400$ (figure 1(b)), $N=800$ (figure 1(c)) and $N=1600$ (figure 1(d)). The $L_{0}$ line represents the theoretical limit $M=K$, under which a perfect reconstruction is impossible, whereas the $L_{1}$ line was obtained in [9, 35] using the replica method in the limit $N\rightarrow\infty$ and $M\rightarrow\infty$, with $\alpha$ finite. The white region is the one in which EP does not converge.

The plots suggest that there exists a phase transition line $\{(\rho,\alpha^{EP}(\rho)),0\leq\rho\leq 1\}$, which is located under the $L_{1}$ line. In order to obtain the coordinates of the points along the transition line, one can proceed numerically by using a bisection-like algorithm. After discretizing the interval $0\leq\rho\leq 1$, one can select two starting values of $\alpha$ for each discretized value $\rho_{0}$. For instance, a reasonable choice, which is the one we adopt, is taking $\alpha_{0}$ on the $L_{0}$ minimization line, namely $\alpha_{0}(\rho_{0})=\rho_{0}$, and $\alpha_{1}(\rho_{0})$ on the $L_{1}$ minimization line, where $\alpha_{1}(\rho_{0})$ is expressed as [9]:

[TABLE]

In the last equation, $\hat{\chi}$ is given by the solution of [35]:

[TABLE]

where $H(x)=\int_{x}^{+\infty}\exp(-\frac{t^{2}}{2})/(2\pi)dt$. Then, one performs the EP inference for configurations corresponding to those points and to the point $(\rho^{*},\alpha^{*})$, where $\rho^{*}=\rho_{0}$ and $\alpha^{*}=(\alpha_{0}+\alpha_{1})/2$ and computes their mean squared error. If the difference $|MSE(\alpha_{1})-MSE(\alpha^{*})|$ is negligible (i. e. smaller than a certain threshold $\delta$), then we set $\alpha_{1}=\alpha^{*}$. Otherwise, we set $\alpha_{0}=\alpha^{*}$. We recompute the middle point and repeat the procedure until we reach the desired accuracy $\Delta\alpha_{min}$ on the points located at the boundary between the successful and unsuccessful reconstruction regions. We summarize the procedure in Algorithm 2. The resulting transition line is shown in figure 2.

We define probability of convergence as the empirical frequency that for a random instance of the signal $\bm{w}$ and of the measurement matrix $\mathbf{F}$, the algorithmic error $\epsilon_{t}$ becomes arbitrarily small after some iteration, as the maximum number of iterations $t$ is increased. In practice, we set a threshold for the error $\epsilon_{t}$ equal to $10^{-6}$ and we estimated the probability of convergence as the fraction of times the algorithm fulfilled the convergence criterion. Empirically, it turns out that the probability of convergence of EP, increases with the number of variables $N$ (not shown). Moreover, the fluctuations of the Pearson correlation coefficient $r$ and of the MSE beyond the transition line decrease as the number of variables $N$ becomes larger, whereas their average values do not seem to depend on the size of the system. We show this in the case of the MSE in figure 3, for $N=400$ and $N=1600$.

Finally, we note that the mean squared error can be expressed as follows:

[TABLE]

where $MSE_{1}$ is the mean squared error between the vector of the $K$ non zero components of the signal and the corresponding vector extracted from the inferred signal and $MSE_{2}$ is the mean squared error between the original and inferred vectors having the remaining $N-K$ elements as components. The latter corresponds to the squared norm of the last $N-K$ components of the inferred vector (divided by $N-K$). Beyond the CS threshold of EP, the dominant contribution to the reconstruction error comes from the estimate of the $K$ non-zero components of the reconstructed vector, implying that, overall, EP is still quite accurate in discriminating the zero entries of the signal. This is shown in figures 4 and 4, in which, respectively, $\rho MSE_{1}$ and $(1-\rho)MSE_{2}$ are compared to the total mean squared error.

3.2 Correlated measurement matrices

We consider the case of correlated measurement matrices $\mathbf{F}$:

[TABLE]

whose rows $\bm{f}_{i}\in\mathbb{R}^{N}$ are correlated but linearly independent samples drawn from a multivariate Gaussian distribution:

[TABLE]

The covariance matrix $\mathbf{S}$ is designed according to the following functional form:

[TABLE]

where $\mathbf{Y}$ is a $k\times N$ matrix with random i.i.d. Gaussian $\mathcal{N}(0,1)$ entries and having a controllable rank $k$ and $\bm{\Delta}$ is a diagonal $N\times N$ matrix with positive Gaussian i.i.d. eigenvalues. Notice that the product $\mathbf{Y}^{T}\mathbf{Y}$ is symmetric and positive semi-definite by construction. Adding the matrix $\bm{\Delta}$ ensures that $\mathbf{S}$ has maximum rank.

We first study the retrieval performance of EP and of Expectation Maximization Belief Propagation (EMBP), a similar message passing reconstruction algorithm [24, 25], implemented in MATLAB and available at http://aspics.krzakala.org.

The BP approximation lies on the independence of the entries of the sensing matrix, a condition that is not generally fulfilled in the matrices considered in this section. In particular, for small values of $k$ the covariances and the variances of $\mathbf{S}$ are of the same order of magnitude. However, as $k$ increases, these variances become dominant with respect to the off-diagonal entries of $\mathbf{S}$ and the associated multivariate Gaussian measure becomes more and more similar to the product of independent univariate Gaussian distributions. In figures 5 and 5, we compare the MSE associated with EMBP and with EP when using correlated matrices and for $\rho=0.3$, $\rho=0.5$ respectively. The signal density is learned in both cases, using Expectation Maximization in the case of EMBP and using the free energy-based variational method described in C in the case of EP. Only converged trials have been taken into account. While BP fails to correctly reconstruct the signal and to infer the signal density in the presence of the correlations we introduced in the measurements, the performance of the EP reconstruction is unaffected. As we expected, EMBP performances improve as $k$ increases. However, we note that at low enough values of $k$, such as those we have considered, EMBP never achieves zero MSE, not even when $\alpha=1$ (that is when we have as many equations as variables).

The fraction of converged trials is generally far lower in the case of EMBP than in the case of EP and decreases as $N$ is increased. For example, in the case of $k=1$, for $N=50$, it is of the order of one in a thousand of simulated trials, two orders of magnitude lower than the fraction of converged EP simulations at the same values of $N$ and $k$ (not shown).

We also plot the MSE resulting from the EP reconstruction for i.i.d. measurement matrices and correlated measurement matrices that are constructed using equations (20), (21) and (22) for various values of $k$. By considering $k=50$, $k=20$, $k=10$, $k=5$ and $k=1$, we obtain that the associated mean squared errors do not exhibit significant discrepancies, as shown in figures 5 and 5. This confirms that the EP inference of the signal is not altered with respect to the case of i.i.d. sensing matrices and retains its correct-incorrect reconstruction threshold.

Finally, we tested several algorithms for sparse reconstruction on linear systems with the same type of correlated measurement matrices considered so far. More precisely, these algorithms are Basis Pursuit [36], Orthogonal Matching Pursuit (OMP) [37], Regularized Orthogonal Matching Pursuit [38], Compressive Sampling Matching Pursuit (CoSaMP) [39], Subspace Pursuit [40], Smoothed L0 (SL0) [41], Approximate Message Passing (AMP) [13] and, again, Expectation Maximization Belief Propagation [24]. These algorithms are implemented in the C++ library KL1p [42]. This specific implementation makes use of the linear algebra library Armadillo [43, 44].

In order to compare the performance of these algorithms, we generated $N_{t}=100$ random gaussian i.i.d. signals of length $N=100$ and as many random correlated sensing matrices, with $k=5$. For any given pair of signal $\bm{w}$ and measurement matrix $\mathbf{F}$, we attempted to recover the original signal by means of EP and of the algorithms included in KL1p. The results are presented in figure 6. As we can see in figure 6 and as further highlighted in the semi-logarithmic plot in figure 6, EP is the only algorithm exhibiting an incorrect-correct reconstruction phase transition, whereas all the other methods that we considered fail to retrieve the signal regardless of the value of $\alpha$. In terms of running time, EP appears to be comparable to most of the other reconstruction techniques, as shown in figure 6.

4 Conclusions

We have proposed an EP-based scheme for efficient CS reconstruction whose computational complexity is dominated by a matrix inversion per iteration which requires $O(N^{3})$ operations

By analyzing the reconstruction achieved by EP in the case of undersampled linear systems with random i.i.d. measurement matrices, we showed that EP exhibits a phase transition, in analogy to other message passing inspired algorithms. We numerically computed the threshold and the related phase transition diagram and found that signal reconstruction is possible below the $L_{1}$ minimization line (see figures 1, and 2).

Finally, we investigated the case of correlated measurement matrices and found that the EP threshold persists, implying that EP is still capable of accurately retrieving the signal beyond a critical $\alpha$ and that, contrary to the case of other reconstruction algorithms, it is robust against the presence of statistical structure in the measurements.

5 Acknowledgements

We acknowledge the kind hospitality of the Physics Department of the University of Havana where part of this work has been completed. This work has received support by INFERNET – Marie Skłodowska-Curie grant agreement number 734439 – and by the SmartData@PoliTO interdepartmental center of Politecnico di Torino.

Appendix A Expectation Propagation in the zero temperature limit

In section 2, we implemented the linear constraints in equation (1) as the $\beta\rightarrow\infty$ limit of the multivariate Gaussian measure. We are going to show here how one can compute analytically this limit. Let us rewrite the previous formalism in a slightly different matrix notation. We will assume that the rank of the measurement matrix is maximum and equal to $M$, with $M<N$ (if this is not the case, we can easily remove the linearly dependent rows from the measurement matrix $\mathbf{F}$). Through Gaussian elimination, we can transform the matrix $\mathbf{F}$ to a row echelon form:

[TABLE]

The structure of the linear constraint induced by the row echelon representation suggests to split the $\bm{w}$ variable into two sets of variables: the first $M$ variables (dependent) and a second set of $N-M$ variables (independent). To do so, we define:

[TABLE]

where ${\bm{w}}^{(d)}\in\mathbb{R}^{M}$, and ${\bm{w}}^{(i)}\in\mathbb{R}^{N-M}$. The linear constraint in equation (1) now reads:

[TABLE]

where ${\bm{y}}^{\prime}\in\mathbb{R}^{M}$ is the transformed measurement vector. Assuming, as in the previous subsection, a Gaussian prior for the $\bm{w}^{(i)}$ variables, there follows a Gaussian statistics on the ${\bm{w}}^{(d)}$ variables with consistent moments:

[TABLE]

where the $\bm{D}_{{\bm{w}}^{(d)}},\bm{D}_{{\bm{w}}^{(i)}}$ are the diagonal matrices whose entries are the inverses of the variances of the Gaussian priors associated with the dependent and independent variables, respectively. We note that from the previous relations the moments of the dependent ${\bm{w}}^{(d)}$ variables can be directly related to those of the ${\bm{w}}^{(i)}$ variables, which allows us to compute everything in terms of the inverse of a smaller matrix of size $\left(N-M\right)\times(N-M)$ compared to the finite temperature case.

At this point, the parameters $\bm{a}$ and $\bm{d}$ can be updated by moment matching as described in section 2.

Appendix B Moments of the spike-and-slab prior

We use the spike-and-slab prior, defined as:

[TABLE]

where $\rho=K/N$ is the density of the signal $\bm{w}$. For spike-and-slab priors, the $n$-th marginal of the tilted distribution (9) is given by:

[TABLE]

in which:

[TABLE]

where we have denoted $\left(\bar{w}_{(n)}\right)_{n}$ and $\left(\Sigma_{(n)}\right)_{n,n}$ by $\bar{w}^{\prime}_{n}$ and $\Sigma^{\prime}_{n}$, respectively, in order to simplify the notation.

The partition function in (24) reads:

[TABLE]

and the first and second moment of $w_{n}$ with respect to $Q^{(n)}$ are given by:

[TABLE]

and by:

[TABLE]

respectively.

Appendix C Learning of the density parameter of the prior

C.1 The EP Free Energy function

Let us consider equation (9) and rewrite it through equation (6) as:

[TABLE]

We define:

[TABLE]

from which the so-called EP free energy follows:

[TABLE]

The converged means $\bm{a}$ and variances $\bm{d}$ of EP, which fulfill the moment matching conditions (11) for all $n=1,\dots,N$, are fixed points of the EP free energy, where the latter are obtained from:

[TABLE]

for $n=1,\dots,N$. In order to show this, we shall prove that the moment matching conditions imply that (32) and (33) are satisfied. We have for $\langle w_{n}\rangle_{Q^{(l)}}$ and $\langle w_{n}^{2}\rangle_{Q^{(l)}}$:

[TABLE]

In the last equality, for $\alpha=1$ ($\alpha=2$), the integral being averaged w.r.t. ${Q^{(l)}(w_{l})}$ is the first (second) moment of $w_{n}$, conditioned on $w_{l}$ and computed w.r.t. $Q$. These moments depend on $w_{l}$ through the mean (squared mean) of $Q(w_{n}|w_{l})$. As such mean depends linearly on $w_{l}$, $\langle w_{n}\rangle_{Q(w_{n}|w_{l})}$ and $\langle w^{2}_{n}\rangle_{Q(w_{n}|w_{l})}$ depend on $w_{l}$ linearly and quadratically, respectively. Therefore, for $\alpha=1,2$, by the moment matching conditions, we have that:

[TABLE]

implying that $\langle w_{n}^{\alpha}\rangle_{Q^{(l)}}=\langle w_{n}^{\alpha}\rangle_{Q}$ and thus that the conditions (32) and (33) are identically fulfilled.

C.2 Learning of the density

We are interested in finding the maximum likelihood value of the density parameter $\rho$ which appears in the prior factors $\psi_{n}(w_{n})$. The likelihood of the parameters of the prior is given by:

[TABLE]

and maximizing this likelihood corresponds to minimizing the associated free energy $F(\rho,\lambda)=-\log Z(\rho,\lambda)$.

At the fixed point of EP, the free energy is approximated by $F_{EP}$ and the parameters can be learned by gradient descent. In particular, we have for the signal density $\rho$:

[TABLE]

where $t$ denotes the current iteration and $\eta$ is a learning rate. In the simulations of this paper, we have taken $\eta=5\times 10^{-4}$.

The parameters of the prior enter in the EP free energy through the contributions associated with each of the tilted distributions. Such contributions read:

[TABLE]

for

[TABLE]

Therefore, we have for $\frac{\partial F_{EP}}{\partial\rho}$:

[TABLE]

where:

[TABLE]

and:

[TABLE]

yielding:

[TABLE]

By taking the derivative w.r.t. $\rho$ one more time, we see that $F_{EP}$ is a strictly convex function of $\rho$ for $\lambda>0$:

[TABLE]

which guarantees that the sought value of $\rho$ is unique at fixed $\bar{w}_{n}$ and $\Sigma_{n,n}$.