Statistical physics-based reconstruction in compressed sensing

Compressed sensing is triggering a major evolution in signal acquisition. It consists in sampling a sparse signal at low rate and later using computational power for its exact reconstruction, so that only the necessary information is measured. Currently used reconstruction techniques are, however, limited to acquisition rates larger than the true density of the signal. We design a new procedure which is able to reconstruct exactly the signal with a number of measurements that approaches the theoretical limit in the limit of large systems. It is based on the joint use of three essential ingredients: a probabilistic approach to signal reconstruction, a message-passing algorithm adapted from belief propagation, and a careful design of the measurement matrix inspired from the theory of crystal nucleation. The performance of this new algorithm is analyzed by statistical physics methods. The obtained improvement is confirmed by numerical studies of several cases.

The ability to recover high dimensional signals using only a limited number of measurements is crucial in many fields, ranging from image processing to astronomy or systems biology.Example of direct applications include speeding up magnetic resonance imaging without the loss of resolution, sensing and compressing data simultaneously [1] and the single-pixel camera [2].Compressed sensing is designed to directly acquire only the necessary information about the signal.This is possible when the signal is sparse in some known basis.In a second step one uses computational power to reconstruct the signal exactly [1,3,4].Currently, the best known generic method for exact reconstruction is based on converting the reconstruction problem into a convex optimization one, which can be solved efficiently using linear programming techniques [3,4].The 1 reconstruction is able to reconstruct accurately, provided the system size is large and the ratio of the number of measurements M to the number of non-zeros K exceeds a specic limit which can be proven by careful analysis [4,5].However, the limiting ratio is signicantly larger than 1.In this paper we improve on the performance of 1 minimization, and in the best possible way: we introduce a new procedure that is able to reach the optimal limit M/K → 1.This procedure, which we call seeded Belief Propagation (s-BP) is based on a new, carefully designed, measurement matrix.It is very powerful, as illustrated in Fig. 1.Its performance will be studied here with a joint use of numerical and analytic studies using methods from statistical physics [6].

Reconstruction in Compressed Sensing
The mathematical problem posed in compressed-sensing reconstruction is easily stated.Given an unknown signal which is a N -dimensional vector s, we make M measurements, where each measurement amounts to a projection of s on some known vector.The measurements are grouped into a M -component vector y, which is obtained from s by a linear transformation y = Fs.Depending on the application, this linear transformation can be for instance associated with measurements of Fourier modes or wavelet coefficients.The observer knows the M × N matrix F and the M measurements y, with M < N .His aim is to reconstruct s.This is impossible in general, but compressed sensing deals with the case where the signal s is sparse, in the sense that only K < N of its components are non-zero.We shall study the case where the non-zero components are real numbers and the measurements are linearly independent.In this case, exact signal reconstruction is possible in principle whenever M ≥ K + 1, using an exhaustive enumeration method which tries to solve y = Fx for all N K possible choices of locations of non-zero components of x: only one such choice gives a consistent linear system, which can then be inverted.However, one is typically interested in large instances where N 1, with M = αN and K = ρ 0 N .The enumeration method solves the compressed sensing problem in the regime where measurement rates are at least as large as the signal density, α ≥ ρ 0 , but in a time which grows exponentially with N , making it totally impractical.Therefore α = ρ 0 is the fundamental reconstruction limit for perfect reconstruction in the noiseless case, when the non-zero components of the signal are real numbers drawn from a continuous distribution.A general and detailed discussion of information-theoretically optimal reconstruction The first row is obtained with the 1-reconstruction algorithm [3,4] for a measurement matrix with iid elements of zero mean and variance 1/N .The second row is obtained with belief propagation, using exactly the same measurements as used in the 1 reconstruction.The third row is the result of the seeded belief propagation, introduced in this work, which uses a measurement matrix based on a chain of coupled blocks, see Fig. 4. The running time of all the three algorithms is comparable (asymptotically they are all quadratic in the size of the signal).In the second part of the figure we took as the sparse signal the relevant coefficients after two-step Haar transform of the picture of Lena.Again for this signal, with density ρ0 = 0.24, the s-BP procedure reconstructs exactly down to very low measurement rates (details are in Appendix G, the data is available online [7]).
In order to design practical 'low-complexity' reconstruction algorithms, it has been proposed [3,4] to find a vector x which has the smallest 1 norm, N i=1 |x i |, within the subspace of vectors which satisfy the constraints y = Fx, using efficient linear programming techniques.In order to measure the performance of this strategy one can focus on the measurement matrix F generated randomly, with independent Gaussian-distributed matrix elements of mean zero and variance 1/N , and the signal vector s having density 0 < ρ 0 < 1.The analytic study of the 1 reconstruction in the thermodynamic limit N → ∞ can be done using either geometric or probabilistic methods [5,11], or with the replica method [12][13][14].It shows the existence of a sharp phase transition at a value α 1 (ρ 0 ).When α is larger than this value, the 1 reconstruction gives the exact result x = s with probability going to one in the large N limit, when α < α 1 (ρ 0 ) the probability that it gives the exact result goes to zero.As shown in Fig. 2, α 1 (ρ 0 ) > ρ 0 and therefore the 1 reconstruction is suboptimal: it requires more measurements than would be absolutely necessary, in the sense that, if one were willing to do brute-force combinatorial optimization, no more than ρ 0 N measurements are necessary.
We introduce a new measurement and reconstruction approach, s-BP, that allows to reconstruct the signal by a practical method, which needs only ≈ ρ 0 N measurements.We shall now discuss its three ingredients: 1) a probabilistic approach to signal reconstruction, 2) a message-passing algorithm adapted from belief propagation [15], which is a procedure known to be efficient in various hard computational problems [16,17], and 3) an innovative design of the measurement matrix inspired from the theory of crystal nucleation in statistical physics and from recent developments in coding theory [18][19][20][21].Some previous works on compressed sensing have used these ingredients separately.In particular, adaptations of belief propagation have been developed for the compressed sensing reconstruction, both in the context of 1 reconstruction [11,22,23], and in a probabilistic approach [24].The idea of seeding matrices in compressed sensing was introduced in [25].It is however only the combined use of these three ingredients that allows us to reach the α = ρ 0 limit.αL1 αEM-BP α=ρ FIG. 2: Phase diagrams for compressed sensing reconstruction for two different signal distributions.On the left-hand side the ρ0N non-zero components of the signal are independent Gaussian random variables with zero mean and unit variance.On the right-hand side they are independent ±1 variables.The measurement rate is α = M/N .On both sides we show, from top to bottom: (a) The phase transition α 1 for 1 reconstruction [5,11,12] (which does not depend on the signal distribution).(b) The phase transition αEM−BP for EM-BP reconstruction based for both sides on the probabilistic model with Gaussian φ.(c) The data points which are numerical reconstruction thresholds obtained with the s-BP procedure with L = 20.The point gives the value of α where exact reconstruction was obtained in 50% of the tested samples, the top of the error bar corresponds to a success rate of 90%, the bottom of the bar to a success of 10%.The shrinking of the error bar with increasing N gives numerical support to the existence of the phase transition that we have studied analytically.These empirical reconstruction thresholds of s-BP are quite close to the α = ρ0 optimal line, and get closer to it when increasing N .The parameters used in these numerical experiments are detailed in Appendix E. (d) The line α = ρ0 that is the theoretical reconstruction limit for signals with continuous φ0.An alternative presentation of the same data using the convention of Donoho and Tanner [5] is shown in Appendix F.

A probabilistic approach
For the purpose of our analysis, we consider the case where the signal s has independent identically distributed (iid) components: , with 0 < ρ 0 < 1.In the large-N limit the number of non-zero components is ρ 0 N .Our approach handles general distributions φ 0 (s i ).
Instead of using a minimization procedure, we shall adopt a probabilistic approach.We introduce a probability measure P (x) over vectors x ∈ R N which is the restriction of the Gauss-Bernoulli measure P (x) = N i=1 [(1−ρ)δ(x i )+ ρφ(x i )] to the subspace |y − Fx| = 0 [26].In this paper, we use a distribution φ(x) which is a Gaussian with mean x and variance σ 2 , but other choices for φ(x) are possible.It is crucial to note that we do not require a priori knowledge of the statistical properties of the signal: we use a value of ρ not necessarily equal to ρ 0 , and the φ that we use is not necessarily equal to φ 0 .The important point is to use ρ < 1 (which reflects the fact that one searches a sparse signal).
Assuming that F is a random matrix, either where all the elements are drawn as independent Gaussian random variables with zero mean and the same variance, or of the carefully-designed type of 'seeding matrices' described below, we demonstrate in Appendix A that, for any ρ 0 -dense original signal s, and any α > ρ 0 the probability P (s) of the original signal goes to one when N → ∞.This result holds independently of the distribution φ 0 of the original signal, which does not need to be known.In practice, we see that s also dominates the measure when N is not very large.In principle, sampling configurations x proportionally to the restricted Gauss-Bernoulli measure P (x) thus gives asymptotically the exact reconstruction in the whole region α > ρ 0 .This idea stands at the roots of our approach, and is at the origin of the connection with statistical physics (where one samples with the Boltzmann measure).FIG.3: When sampling from the probability P (x), in the limit of large N , the probability that the reconstructed signal x is at a squared distance D = i (xi − si) 2 /N from the original signal s is written as ce N Φ(D) where c is a constant and Φ(D) is the free entropy.Left: Φ(D) for a Gauss-Bernoulli signal with ρ0 = ρ = 0.4 in the case of an unstructured measurement matrix F with independent random Gaussian-distributed elements.The evolution of the EM-BP algorithm is basically a steepest ascent in Φ(D) starting from large values of D. It goes to the correct maximum at D = 0 for large value of α but is blocked in the local maximum that appears for α < αEM−BP(ρ0 = 0.4) ≈ 0.59.For α < ρ0, the maximum is not at D = 0 and exact inference is impossible.The seeding matrix F, leading to the s-BP algorithm, succeeds in eliminating this local maximum.Right: Convergence time of the EM-BP and s-BP algorithms obtained through the replica analysis for ρ0 = 0.4.The EM-BP convergence time diverges as α → αEM−BP with the standard L = 1 matrices.The s-BP strategy allows to go beyond the threshold: using α1 = 0.7 and increasing the structure of the seeding matrix (here L = 2, 5, 10, 20), we approach the limit α = ρ0 (details of the parameters are given in Appendix E).

Sampling with Expectation Maximization Belief Propagation
The exact sampling from a distribution such as P (x), eq.(A2), is known to be computationally intractable [27].However, an efficient approximate sampling can be performed using a message-passing procedure that we now describe [24,[28][29][30].We start from the general belief-propagation formalism [17,31,32]: for each measurement µ = 1, . . ., M and each signal component i = 1, . . ., N , one introduces a 'message' m i→µ (x i ) which is the probability of x i in a modified measure where measurement µ has been erased.In the present case, the canonical belief propagation equations relating these messages can be simplified [11,[22][23][24]33] into a closed form that uses only the expectation a (t) i→µ and the variance v ).An important ingredient that we add to this approach is the learning of the parameters in P (x): the density ρ, and the mean x and variance σ 2 of the Gaussian distribution φ(x).These are three parameters to be learned using update equations based on the gradient of the so-called Bethe free entropy, in a way analogous to the expectation maximization [34][35][36].This leads to the Expectation Maximization Belief Propagation (EM-BP) algorithm that we will use in the following for reconstruction in compressed sensing.It consists in iterating the messages and the three parameters, starting from random messages a (0) i→µ and v (0) i→µ , until a fixed point is obtained.Perfect reconstruction is found when the messages converge to the fixed point a i→µ = s i and v i→µ = 0.
Like the 1 reconstruction, the EM-BP reconstruction also has a phase transition.Perfect reconstruction is achieved with probability one in the large-N limit if and only if α > α EM−BP .Using the asymptotic replica analysis, as explained below, we have computed the line α EM−BP (ρ 0 ) when the elements of the M × N measurement matrix F are independent Gaussian random variables with zero mean and variance 1/N and the signal components are iid.The location of this transition line does depend on the signal distribution, see Fig. 2, contrary to the location of the 1 phase transition.
Notice that our analysis is fully based on the case when the probabilistic model has a Gaussian φ.Not surprisingly, EM-BP performs better when φ = φ 0 , see the left-hand side of Fig. 2, where EM-BP provides a sizable improvement over 1 .In contrast, the right-hand side of Fig. 2 shows an adversary case when we use a Gaussian φ to reconstruct a binary φ 0 , in this case there is nearly no improvement over 1 reconstruction.

Designing seeding matrices
In order for the EM-BP message-passing algorithm to be able to reconstruct the signal down to the theoretically optimal number of measurements α = ρ 0 , one needs to use a special family of measurement matrices F that we call 'seeding matrices'.If one uses an unstructured F, for instance a matrix with independent Gaussian-distributed random elements, EM-BP samples correctly at large α, but at small enough α a metastable state appears in the measure P (x), and the EM-BP algorithm is trapped in this state, and is therefore unable to find the original signal (see Fig. 3), just as a supercooled liquid gets trapped in a glassy state instead of crystallizing.It is well known in crystallization theory that the crucial step is to nucleate a large enough seed of crystal.This is the purpose of the following design of F.
We divide the N variables into L groups of N/L variables, and the M measurements into L groups.The number of measurements in the p-th group is We then choose the matrix elements F µi independently, in such a way that, if i belongs to group p and µ to group q then F µi is a random number chosen from the normal distribution with mean zero and variance J q,p /N (see Fig. 4).The matrix J q,p is a L × L coupling matrix (and the standard compressed sensing matrices are obtained using L = 1 and α 1 = α).Using these new matrices, one can shift the BP phase transition very close to the theoretical limit.In order to get an efficient reconstruction with message passing, one should use a large enough α 1 .With a good choice of the coupling matrix J p,q , the reconstruction first takes place in the first block, and propagates as a wave in the following blocks p = 2, 3, . . ., even if their measurement rate α p is small.In practice, we use The whole reconstruction process is then analogous to crystal nucleation, where a crystal is growing from its seed (see Fig. 5).Similar ideas have been used recently in the design of sparse coding matrices for error-correcting codes [18][19][20][21].

Analysis of the performance of the seeded Belief Propagation procedure
The s-BP procedure is based on the joint use of seeding measurement matrices and of the EM-BP message-passing reconstruction.We have studied it with two methods: direct numerical simulations and analysis of the performance in the large N limit.The analytical result was obtained by a combination of the replica method and of the cavity method (also known as 'density evolution' or 'state evolution').The replica method is a standard method in statistical physics [6], which has been applied successfully to several problems of information theory [17,29,37,38] including compressed sensing [12][13][14].It can be used to compute the free entropy function Φ associated with the probability P (x) (see Appendix D), and the cavity method shows that the dynamics of the message-passing algorithm is a gradient dynamics leading to a maximum of this free-entropy.
When applied to the usual case of the full F matrix with independent Gaussian-distributed elements (case L = 1), the replica computation shows that the free-entropy Φ(D) for configurations constrained to be at a mean-squared distance D has a global maximum at D = 0 when α > ρ 0 , which confirms that the Gauss-Bernoulli probabilistic reconstruction is in principle able to reach the optimal compression limit α = ρ 0 .However, for α EM−BP > α > ρ 0 , where α EM−BP is a threshold that depends on the signal and on the distribution P (x), a secondary local maximum of Φ(D) appears at D > 0 (see Fig. 3).In this case the EM-BP algorithm converges instead to this secondary maximum and does not reach exact reconstruction.The threshold α EM−BP is obtained analytically as the smallest value of α such that Φ(D) is decreasing (Fig. 2).This theoretical study has been confirmed by numerical measurements of the number of iterations needed for EM-BP to reach its fixed point (within a given accuracy).This convergence time of BP to the exact reconstruction of the signal diverges when α → α EM−BP (see Fig. 3).For α < α EM−BP the EM-BP algorithm converges to a fixed point with strictly positive mean-squared error (MSE).This 'dynamical' transition is FIG. 4: Construction of the measurement matrix F for seeded compressed sensing.The elements of the signal vector are split into L (here L = 8) equally-sized blocks, the number of measurements in each block is Mp = αpN/L (here α1 = 1, αp = 0.5 for p = 2, . . ., 8).The matrix elements Fµi are chosen as random Gaussian variables with variance Jq,p/N if variable i is in the block p and measurement µ in the block q.In the s-BP algorithm we use Jp,q = 0, except for Jp,p = 1, Jp,p−1 = J1, and Jp−1,p = J2.Good performance is typically obtained with relatively large J1 and small J2.
similar to the one found in the cooling of liquids which go into a super-cooled glassy state instead of crystallizing, and appears in the decoding of error correcting codes [16,17] as well.
We have applied the same technique to the case of seeding-measurement matrices (L > 1).The cavity method allows to analytically locate the dynamical phase transition of s-BP.In the limit of large N , the MSE E p and the variance messages V p in each block p = 1, . . .L, the density ρ, the mean x, and the variance σ 2 of P (x) evolve according to a dynamical system which can be computed exactly (see Appendix E), and one can see numerically if this dynamical system converges to the fixed point corresponding to exact reconstruction (E p = 0 for all p).This study can be used to optimize the design of the seeding matrix F by choosing α 1 , L and J p,q in such a way that the convergence to exact reconstruction is as fast as possible.In Fig. 3 we show the convergence time of s-BP predicted by the replica theory for different sets of parameters.For optimized values of the parameters, in the limit of a large number of blocks L, and large system sizes N/L, s-BP is capable of exact reconstruction close to the smallest possible number of measurements, α → ρ 0 .In practice, finite size effects slightly degrade this asymptotic threshold saturation, but the s-BP algorithm nevertheless reconstructs signals at rates close to the optimal one regardless of the signal distribution, as illustrated in Fig. 2.
We illustrate the evolution of the s-BP algorithm in Fig. 5.The non-zero signal elements are Gaussian with zero mean, unit variance, density ρ 0 = 0.4 and measurement rate α = 0.5, which is deep in the glassy region where all other known algorithms fail.The s-BP algorithm first nucleates the native state in the first block and then propagates it through the system.We have also tested the s-BP algorithm on real images where the non-zero components of the signal are far from Gaussian, and the results are nevertheless very good, as shown in Fig. 1.This shows that the quality of the result is not due to a good guess of P (x).It is also important to mention that the gain in performance in using seeding-measurement matrices is really specific to the probabilistic approach: we have computed the phase diagram of 1 minimization with these matrices and found that, in general, the performance is slightly degraded with respect to the one of the full measurement matrices, in the large N limit.This demonstrates that it is the combination of the probabilistic approach, the message-passing reconstruction with parameter learning, and the seeding design of the measurement matrix that is able to reach the best possible performance.

Perspectives
The seeded compressed sensing approach introduced here is versatile enough to allow for various extensions.One aspect worth mentioning is the possibility to write the EM-BP equations in terms of N messages instead of the M × N parameters as described in Appendix C.This is basically the step that goes from rBP [24] to AMP [11] algorithm.FIG.5: Evolution of the mean-squared error at different times as a function of the block index for the s-BP algorithm.Exact reconstruction first appears in the left block whose rate α1 > αEM−BP allows for seeded nucleation.It then propagates gradually block by block driven by the free entropy difference between the metastable and equilibrium state.After little more than 300 iterations the whole signal of density ρ0 = 0.4 and size N = 50000 is exactly reconstructed well inside the zone forbidden for BP (see Fig. 2).Here we used L = 20, J1 = 20, J2 = 0.2, α1 = 1.0, α = 0.5.
It could be particularly useful when the measurement matrix has some special structure, so that the measurements y can be obtained in many fewer than M × N operations (typically in N log N operations).We have also checked that the approach is robust to the introduction of a small amount of noise in the measurements (see Appendix H).Finally, let us mention that, in the case where a priori information on the signal is available, it can be incorporated in this approach through a better choice of φ, and considerably improve the performance of the algorithm.For signal with density ρ 0 , the worst case, that we addressed here, is when the non-zero components of the signal are drawn from a continuous distribution.Better performance can be obtained with our method if these non-zero components come from a discrete distribution and one uses this distribution in the choice of φ.Another interesting direction in which our formalism can be extended naturally is the use of non-linear measurements and different type of noises.Altogether, this approach turns out to be very efficient both for random and structured data, as illustrated in Fig. 1, and offers an interesting perspective for fpractical compressed sensing applications.Data and code are available online at [7].
Appendix A: Proof of the optimality of the probabilistic approach Here we give the main lines of the proof that our probabilistic approach is asymptotically optimal.We consider the case where the signal s has iid components with 0 < ρ 0 < 1.And we study the probability distribution with a Gaussian φ(x) of mean zero and unit variance.We stress here that we consider general φ 0 , i.e. φ 0 is not necessarily equal to φ(x) (and ρ 0 is not necessarily equal to ρ).The measurement matrix F is composed of iid elements F µi such that if µ belongs to block q and i belongs to block p then F µi is a random number generated from the Gaussian distribution with zero mean and variance J q,p /N.The function δ (x) is a centered Gaussian distribution with variance 2 .We show that, with probability going to one in the large N limit (at fixed α = M/N ), the measure P (obtained with a generic seeding matrix F as described in the main text) is dominated by the signal if α > ρ 0 , α > ρ 0 (as long as φ 0 (0) is finite).
We introduce the constrained partition function: and the corresponding 'free entropy density': The notations E F and E s denote respectively the expectation value with respect to F and to s.I denotes an indicator function, equal to one if its argument is true, and equal to zero otherwise.The proof of optimality is obtained by showing that, under the conditions above, lim →0 Y(D, )/[(α − ρ 0 ) log(1/ )] is finite if D = 0 (statement 1), and it vanishes if D > 0 (statement 2).This proves that the measure P is dominated by D = 0, i.e. by the neighborhood of the signal x i = s i .The standard 'self-averageness' property, which states that the distribution (with respect to the choice of F and s) of log Y (D, )/N concentrates around Y(D, ) when N → ∞, completes the proof.We give here the main lines of the first two steps of the proof.
We first sketch the proof of statement 2. The fact that lim →0 Y(D, )/[(α − ρ 0 ) log(1/ )] = 0 when D > 0 can be derived by a first moment bound: where Y ann (D, ) is the 'annealed partition function' defined as In order to evaluate Y ann (D, ) one can first compute the annealed partition function in which the distances between x and the signal are fixed in each block.More precisely, we define By noticing that the M random variables a µ = i F µi (x i − s i ) are independent Gaussian random variables one obtains: where The behaviour of ψ(r) is easily obtained by standard saddle point methods.In particular, when r → 0, one has ψ(r) 1 2 ρ 0 log r.Using (A6), we obtain, in the small limit: where the maximum over r 1 , . . ., r L is to be taken under the constraint r 1 + • • • + r L > LD.Taking the limit of → 0 with a finite D, at least one of the distance r p must remain finite.It is then easy to show that lim where α is the fraction of measurements in blocks 2 to L. As α > ρ 0 , this is less singular than log(1/ )(α − ρ 0 ), which proves statement 2.
On the contrary, when D = 0, we obtain from the same analysis lim This annealed estimate actually gives the correct scaling at small , as can be shown by the following lower bound.When D = 0, we define V 0 as the subset of indices i where s i = 0, |V 0 | = N (1 − ρ 0 ), and V 1 as the subset of indices i where s i = 0, |V 1 | = N ρ 0 .We obtain a lower bound on Y (0, ) by substituting P (x) by the factors (1 − ρ)δ(x i ) when i ∈ V 0 and ρφ(x i ) when i ∈ V 1 .This gives: where M ij = αN µ=1 F µi F µj .The matrix M , of size ρ 0 N × ρ 0 N , is a Wishart-like random matrix.For α > ρ 0 , generically, its eigenvalues are strictly positive, as we show below.Using this property, one can show that, if i∈V1 φ(s i ) > 0, the integral over the variables u i in (A11) is strictly positive in the limit → 0. The divergence of Y(0, ) in the limit → 0 is due to the explicit term exp[N (ρ 0 − α) log ] in Eq. (A11).The fact that all eigenvalues of M are strictly positive is well known in the case of L = 1 where the spectrum has been obtained by Marcenko and Pastur.In general, the fact that all the eigenvalues of M are strictly positive is equivalent to saying that all the lines of the αN × ρ 0 N matrix F (which is the restriction of the measurement matrix to columns with non-zero signal components) are linearly independent.In the case of seeding matrices with general L, this statement is basically obvious by construction of the matrices, in the regime where in each block q, α q > ρ 0 and J qq > 0. A more formal proof can be obtained as follows.We consider the Gaussian integral This quantity is finite if and only if v is smaller than the smallest eigenvalue, λ min , of M .We now compute the annealed average Z ann (v) = E F Z(v).If Z ann (v) is finite, then the probability that λ min ≤ v goes to zero in the large N limit.Using methods similar to the one above, one can show that The saddle point equations have a solution at v = 0 (and by continuity also at v > 0 small enough), when 1 L L q=1 αq ρ0 = α ρ0 > 1 (it can be found for instance by iteration).Therefore, 2L n log E F Z(v) is finite for some v > 0 small enough, and therefore λ min > 0.

Appendix B: Derivation of Expectation maximization Belief Propagation
In this and the next sections we present the message-passing algorithm that we used for reconstruction in compressed sensing.In this section we derive its message-passing form, where O(N M ) messages are being sent between each signal component i and each measurement µ.This algorithm was used in [24], where it was called the relaxed belief propagation, as an approximate algorithm for the case of a sparse measurement matrix F. In the case that we use here of a measurement matrix which is not sparse (a finite fraction of the elements of F is non-zero, and all the non-zero elements scale as 1/ √ N ), the algorithm is asymptotically exact.We show here for completeness how to derive it.In the next section we then derive asymptotically equivalent equations that depend only on O(N ) messages.In statistical physics terms, this corresponds to the TAP equations [28] with the Onsager reaction term, that are asymptotically equivalent to the BP on fully connected models.In the context of compressed sensing this form of equations has been used previously [11] and it is called approximate message passing (AMP).In cases when the matrix F can be computed recursively (e.g. via fast Fourier transform), the running time of the AMP-type message passing is O(N log N ) (compared to the O(N M ) for the non-AMP form).Apart for this speed-up, both classes of message passing give the same performance.
We derive here the message-passing algorithm in the case where measurements have additive Gaussian noise, the noiseless case limit is easily obtained in the end.The posterior probability of x after the measurement of y is given by where Z is a normalization constant (the partition function) and ∆ µ is the variance of the noise in measurement µ.
The noiseless case is recovered in the limit ∆ µ → 0. The optimal estimate, that minimizes the MSE with respect to the original signal s, is obtained from averages of x i with respect to the probability measure P (x).Exact computation of these averages would require exponential time, belief propagation provides a standard approximation.The canonical BP equations for probability measure P (x) read where Z µ→i and Z i→µ are normalization factors ensuring that dx i m µ→i (x i ) = dx i m i→µ (x i ) = 1.These are integral equations for probability distributions that are still practically intractable in this form.We can, however, take advantage of the fact that after proper rescaling the linear system y = Fx is such a way that elements of y and x are of O(1), the matrix F µi has random elements with variance of O(1/N ).Using the Hubbard-Stratonovich transformation for ω = ( j =i F µj x j ) we can simplify eq.(B2) as Now we expand the last exponential around zero because the term F µj is small in N , we keep all terms that are of O(1/N ).Introducing means and variances as new messages we obtain Performing the Gaussian integral over λ we obtain where we introduced and the normalization Zµ→i contains all the x i -independent factors.The noiseless case corresponds to ∆ µ = 0.
To close the equations on messages a i→µ and v i→µ we notice that Messages a i→µ and v i→µ are respectively the mean and variance of the probability distribution m i→µ (x i ).For general φ(x i ) the mean and variance (B6-B7) will be computed using numerical integration over x i .Eqs. (B6-B7) together with (B10-B11) and (B12) then lead to closed iterative message-passing equations.In all the specific examples shown here and in the main part of the paper we used a Gaussian φ(x i ) with mean x and variance σ 2 .We define two functions Then the closed form of the BP update is where the a i and v i are the mean and variance of the marginal probabilities of variable x i .As we discussed in the main text the parameters ρ, x and σ are usually not known in advance.However, their values can be learned within the probabilistic approach.A standard way to do so is called expectation maximization [34].One realizes that the partition function is proportional to the probability of the true parameters ρ 0 , s, σ 0 , given the measurement y.Hence to compute the most probable values of parameters one searches for the maximum of this partition function.Within the BP approach the logarithm of the partition function is the Bethe free entropy expressed as [17] F where The stationarity conditions of Bethe free entropy (B18) with respect to ρ leads to where U i = γ A γ→i , and V i = γ B γ→i .Stationarity with respect to x and σ gives In statistical physics conditions (B22) are known under the name Nishimori conditions [35,37].In the expectation maximization eqs.(B22-B24) they are used iteratively for the update of the current guess of parameters.A reasonable initial guess is ρ init.= α.The value of ρ 0 s can also be obtained with a special line of measurement consisting of a unit vector, hence we assume that given estimate of ρ the x = ρ 0 s/ρ.In the case where the matrix F is random with Gaussian elements of zero mean and variance 1/N , we can also use for learning the variance: Appendix C: AMP-form of the message passing In the large N limit, the messages a i→µ and v i→µ are nearly independent of µ, but one must be careful to keep the correcting Onsager reaction terms.Let us define Then we have We now compute ω µ To express ω µ = i F µi a i→µ , we see that the first correction term has a contribution in F 3 µi , and can be safely neglected.On the contrary, the second term has a contribution in F 2 µi which one should keep.Therefore The computation of γ µ is similar, it gives: For a known form of matrix F these equations can be slightly simplified further by using the assumptions of the BP approach about independence of F µi and BP messages.This plus a law of large number implies that for matrix F with Gaussian entries of zero mean and unit variance one can effectively 'replace' every F 2 µi by 1/N in eqs.(C3,C4) and (C6,C7).This leads, for homogeneous or bloc matrices, to even simpler equations and a slightly faster algorithm.
Eqs. (C3,C4) and (C6,C7), with or without the later simplification, give a system of closed equations.They are a special form (P (x) and hence functions f a , f c are different in our case) of the approximate message passing of [11].
The final reconstruction algorithm for general measurement matrix and with learning of the P (x) parameters can hence be summarized in a schematic way: EM-BP(y µ , F µi , criterium, t max ) 1 Initialize randomly messages U i from interval [0, 1] for every component; 2 Initialize randomly messages V i from interval [−1, 1] for every component; 3 Initialize messages ω µ ← y µ ; 4 Initialize randomly messages γ µ from interval ]0, 1] for every measurement; 5 Initialize the parameters ρ ← α, x ← 0, σ 2 ← 1. 6 conv ← criterium + 1; t ← 0; 7 while conv > criterium and t < t max : 8 do t ← t + 1; 9 for each component i: Update U i according to eq. (C3).12 Update V i according to eq. (C4).13 for each measurement µ: 14 do Update ω µ according to eq. (C6).15 Update γ µ according to eq. (C7).16 Update ρ according to eq. (B22).17 Update x and σ 2 according to eq. (B24).18 for each component i: Note that in practice we use 'damping' (at each update, the new message is obtained as u times the old value plus 1 − u times the newly computed value, with a damping 0 < u < 1, typically u = 0.5) for both the update of messages and learning of parameters, empirically this speeds up the convergence.Note also that the algorithm is relatively robust with respect to the initialization of messages.The reported initialization was used to obtain the results in Fig. 1 of the main text.However, other initializations are possible.Note also that for specific classes of signals or measurement matrices the initial conditions may be adjusted to take into account the magnitude of the values F µi and y µ .
For a general matrix F one iteration takes O(N M ) steps, we observed the number of iterations needed for convergence to be basically independent of N , however, the constant depends on the parameters and the signal, see Fig. 3 in the main paper.For matrices that can be computed recursively (i.e.without storing all their N M elements) a speed-up is possible, as the message-passing loop takes only O(M + N ) steps.that the density evolution equations are the same for the message-passing and for the AMP equations.It turns out that the above equations close in terms of two parameters, the mean-squared error BP .From eqs. (D4-D7) easily gets a closed mapping In the main text we defined the function Φ(D) which is the free entropy restricted to configurations x for which D = N i=1 (x i − s i ) 2 /N is fixed.This is evaluated as the saddle point over Q, q, Q, q, m of the function Φ(Q, q, (Q − D + ρ 0 s 2 )/2, Q, q, m).This function is plotted in Fig. 3(a) of the main text.
In presence of Expectation Maximization learning of the parameters, the density evolution for the conditions (B22) and (B24) are The density evolution equations now provide a mapping EM−BP , ρ (t) , x (t) , σ (t) (D12) obtained by complementing the previous equations on E EM−BP with the update equations (D9,D10,D11).The next section gives explicitly the full set of equations in the case of seeding matrices, the ones for the full matrices are obtained by taking L = 1.These are the equations that we study to describe analytically the evolution of EM-BP algorithm and obtain the phase diagram for the reconstruction (see Fig. 2 in the main text).
Appendix E: Replica analysis and density evolution for seeding-measurement matrices Many choices of J 1 and J 2 actually work very well, and good performance for seeding-measurement matrices can be easily obtained.In fact, the form of the matrix that we have used is by no means the only one that can produce the seeding mechanism, and we expect that better choices, in terms of convergence time, finite-size effects and sensibility to noise, could be unveiled in the near future.
With the matrix presented in this work, and in order to obtain the best performance (in terms of phase transition limit and of speed of convergence) one needs to optimize the value of J 1 and J 2 depending on the type of signal.Fortunately, this can be analysed with the replica method.The analytic study in the case of seeding measurement matrices is in fact done using the same techniques as for the full matrix.The order parameters are now the MSE E p = q p − 2m p + ρ 0 s 2 and varianceV p = Q p − q p in each block p ∈ {1, . . ., L}.Consequently, we obtain the final dynamical system of 2L + 3 order parameters describing the density evolution of the s-BP algorithm.The order parameters at iteration t + 1 of the message-passing algorithm are given by: where: The functions f a (X, Y ), f c (X, Y ) were defined in (B13-B14), and the function g is defined as This is the dynamical system that we use in the paper in the noiseless case (∆ = 0) in order to optimize the values of α 1 , J 1 and J 2 .We can estimate the convergence time of the algorithm as the number of iterations needed in order to reach the successful fixed point (where all E p and V p vanish within some given accuracy).Figure 6 shows the convergence time of the algorithm as a function of J 1 and J 2 for Gauss-Bernoulli signals.
The numerical iteration of this dynamical system is fast.It allows to obtain the theoretical performance that can be achieved in an infinite-N system.We have used it in particular to estimate the values of L, α 1 , J 1 , J 2 that have good performance.For Gauss-Bernoulli signals, using optimal choices of J 1 , J 2 , we have found that perfect reconstruction can be obtained down to the theoretical limit α = ρ 0 by taking L → ∞ (with correction that scale as 1/L).Recent rigorous work by Donoho, Javanmard and Montanari [39] extends our work and proves our claim that s-BP can reach the optimal threshold asymptotically.
Practical numerical implementation of s-BP matches this theoretical performance only when the size of every block is large enough (few hundreds of variables).In practice, for finite size of the signal, if we want to keep the block-size reasonable we are hence limited to values of L of several dozens.Hence in practice we do not quite saturate the threshold α = ρ 0 , but exact reconstruction is possible very close to it, as illustrated in Fig. 2 in the main text, where the values that we used for the coupling parameters are listed in Table I.
In Fig. 3, we also presented the result of the s-BP reconstruction of the Gaussian signal of density ρ 0 = 0.4 for different values of L. We have observed empirically that the result is rather robust to choices of J 1 , J 2 and α 1 .FIG. 7: (color online): Same data as in Fig. 2 in the main paper, but using the convention of [5].The phase diagrams is now plot as a function of the under-sampling ratio ρDT = K/M = ρ/α and of the over-sampling ratio δDT = M/N = α.

Appendix G: Details on the phantom and Lena examples
In this section, we give a detailed description of the way we have produced the two examples of reconstruction in Fig. 1.It is important to stress that this figure is intended to be an illustration of the s-BP reconstruction algorithm.As such, we have used elementary protocols to produce true K-sparse signals and have not tried to optimize the sparsity nor to use the best possible compression algorithm; instead, we have limited ourselves to the simplest Haar wavelet transform, to make the exact reconstruction and the comparison between the different approaches more transparent.The Shepp-Logan example is a 128 2 picture that has been generated using the Matlab implementation.The Lena picture is a 128 2 crop of the 512 2 gray version of the standard test image.In the first case, we have worked with the sparse one-step Haar transform of the picture, while in the second one, we have worked with a modified picture where we have kept the 24 percent of largest (in absolute value) coefficients of the two-step Haar transform, while putting all others to zero.The datasets of the two images are available online [7].Compressed sensing here is done as follows: The original image is a vector o of N = L 2 pixels.The unknown vector x = Wo are the projections of the original image on a basis of one-or two-steps Haar wavelets.It is sparse by construction.We generate a matrix F as described above, and construct G = FW.The measurements are obtained by y = Go, and the linear system for which one does reconstruction is y = Fx.Once x has been found, the original image is obtained from o = W −1 x.We used EM-BP and s-BP with a Gauss-Bernoulli P (x).
On the algorithmic side, the EM-BP and 1 experiments were run with the same full gaussian random measurement matrices.The minimization of the 1 norm was done using the 1 -magic tool for Matlab, which is a free implementation that can be download at http://users.ece.gatech.edu/˜justin/l1magic/,using the lowest possible tolerance such that the algorithm outputs a solution.The coupling parameters of s-BP are given in Table II.A small damping was used in each iterations, we mixed the old and new messages, keeping 20% of the messages from the previous iteration.Moreover, since for both Lena and the phantom, the components of the signal are correlated, we have permuted randomly the columns of the sensing matrix F .This allows to avoid the (dangerous) situation where a given block contains only zero signal components.
Note that the number of iterations needed by the s-BP procedure to find the solution is moderate.The s-BP algorithm coded in Matlab finds the image in few seconds.For instance, the Lena picture at α = 0.5 requires about 500 iterations and about 30 seconds on a standard laptop.Even for the most difficult case of Lena at α = 0.3, we need around 2000 iterations, which took only about 2 minutes.In fact, s-BP (coded in Matlab) is much faster than the Matlab implementation of 1 -magic on the same machine.We report the MSE for all these reconstruction protocols in Table III.expectation maximization approach, which reads: . (H2) The dynamical system describing the density evolution has been written in Sect.E. It just needs to be complemented by the iteration on ∆ according to Eq. (H2).We have used it to study the evolution of MSE in the case where both P (x) and signal distribution are Gauss-Bernoulli with density ρ, zero mean and unit variance for the full measurement matrix and for the seeding one.Figure 8 shows the MSE as a function of α for a given ∆, using our algorithm, and the

1 FIG. 1 :
FIG. 1: Two illustrative examples of compressed sensing in image processing.Top: the original image, the Shepp-Logan phantom, of size N = 128 2 , is transformed via one step of Haar wavelets into a signal of density ρ0 ≈ 0.15.With compressed sensing one is thus in principle able to reconstruct exactly the image with M ≥ ρ0N measurements, but practical reconstruction algorithms generally need M larger than ρ0N .The five columns show the reconstructed figure obtained from M = αN measurements, with decreasing acquisition rate α.The first row is obtained with the 1-reconstruction algorithm[3,4] for a measurement matrix with iid elements of zero mean and variance 1/N .The second row is obtained with belief propagation, using exactly the same measurements as used in the 1 reconstruction.The third row is the result of the seeded belief propagation, introduced in this work, which uses a measurement matrix based on a chain of coupled blocks, see Fig.4.The running time of all the three algorithms is comparable (asymptotically they are all quadratic in the size of the signal).In the second part of the figure we took as the sparse signal the relevant coefficients after two-step Haar transform of the picture of Lena.Again for this signal, with density ρ0 = 0.24, the s-BP procedure reconstructs exactly down to very low measurement rates (details are in Appendix G, the data is available online[7]).