Sparse Estimation with the Swept Approximated Message-Passing Algorithm

I Introduction

Belief Propagation (BP) is a powerful iterative message passing algorithm for graphical models . However, it presents two main drawbacks when applied to highly connected continuous variable problems: first, the need to work with continuous probability distributions; and second, the necessity to iterate over one such probability distribution for each pair of variables.

The first problem can be addressed by projecting the distributions onto a finite number of moments and the second by utilizing the Thouless-Andreson-Palmer (TAP) approach where only single variable marginals are required. Approximate message passing (AMP), first introduced in , is one relaxation of BP that utilizes both of the aforementioned approximations in order to solve sparse estimation problems. In AMP’s more general setting, as is considered in Generalized AMP (GAMP) , the goal of the algorithm is the reconstruction of an $N$-dimensional sparse vector $\bf x$ given the knowledge of an $M$-dimensional vector $\bf y$ obtained via a possibly non-linear and/or probabilistic output function $h(z)$ performed on a set of linear projections. Specifically,

For example, if $h(z)=z+\xi$ where $\xi$ is a zero-mean i.i.d. Gaussian random variable, then $h(z)$ represents an additive white Gaussian noise (AWGN) channel. With this output function, in the setting $M\ll N$, (1) is simply the application of Compressed Sensing (CS) under noise. AMP is currently acknowledged as one of the foremost algorithms for such problems in terms of both its computational efficiency and in the number of measurements required for exact reconstruction of x. In fact, with properly chosen measurement matrices , one can achieve information-theoretically optimal reconstruction performance for CS, a hitherto unachievable bound with standard convex optimization approaches.

Just as with any iterative algorithm, the convergence properties of AMP are of chief analytical concern. Many rigorous results have been obtained on the performance of AMP in the case of i.i.d. and block i.i.d. matrices . Unfortunately, while AMP performs well for zero-mean i.i.d. projections, performance tends to drastically decline if one moves away from these simple scenarios. In fact, even for i.i.d. matrices with a small positive mean, the algorithm may violently diverge, leading to poor reconstruction results . This instability to slight variations from these strict assumptions on the projections is a serious problem for many practical applications of AMP.

The main theoretical reason for these convergence issues has been identified in . Namely, AMP’s use of a parallel update, instead of a sequential one, on the BP variables at each iteration. Three strategies have been proposed in recent literature to avoid this problem. First, one can highly damp the AMP iterations, as in . However, this often requires a damping factor so large that the cost, in terms of the number of iterations until convergence, is prohibitive. Additionally, it is not entirely clear how to determine an optimal damping factor to ensure convergence in general. Second, one can modify the problem a posteriori in order to come back to a more favorable situation. For instance, one might remove the mean of the matrix and of the measurements , or one might modify the algorithm according to the theoretical spectrum of the operator $\Phi$ , if it is known. This knowledge about the operator may be prohibitive and could therefore present a strong limitation in practice. Third, one might take one step backward in approximation from AMP to a BP-style iteration . This amounts to a huge cost in terms of both memory and computational efficiency as there are $O(N^{2})$ variables to update with BP as opposed to the the $O(N)$ utilized in AMP.

In this contribution, we solve these problems by deriving a modified and efficient AMP algorithm with greatly improved convergence properties while preserving the $O(N)$ iteration and memory cost of AMP. We accomplish this by a careful analysis of the relaxation leading from BP to AMP where we preserve the sequential, or swept, variable update pattern of BP in our AMP approach. This leads to a slightly modified set of update rules for the AMP and GAMP algorithms without affecting the fixed point in any way. The resulting algorithm, which we denote as Swept AMP (SwAMP), possesses impressive empirical convergence properties. The derivation of SwAMP is explained next in Sec. II. We then report, in Sec. III, numerical results for basic and 1-bit CS, as well as for group testing. In all of these cases, huge improvements over the state-of-the-art can be obtained while remaining robust to projections with troublesome properties.

II From Belief-Propagation to SwAMP for Signal Recovery

where we write $\propto$ as we neglect the normalization constant. The likelihood $P(\operatorname{\mathbf{y}}|\Phi,\operatorname{\mathbf{x}})$ is determined according to the constraints one wishes to enforce, which we consider to be of form $\operatorname{\mathbf{y}}=h(\Phi\operatorname{\mathbf{x}})$, with $h$ being, in general, any stochastic function. Here, we consider $h$ to be an AWGN channel One can generalize $h$ to be a more complicated output function. This generalization constitutes the change of AMP to GAMP . For example, we examine the case of 1-bit CS in Sec. III-C where $h$ is a non-linear sign function.,

where $\Delta$ is the variance of the AWGN and $\Phi_{\mu}$ is the $\mu^{th}$ row-vector of $\Phi$. Hence,

The prior $P_{0}(\operatorname{\mathbf{x}})$ is determined from the information we have on the structure of $\operatorname{\mathbf{x}}$. For CS, we are concerned with the recovery of sparse signals, i.e. ones with few non-zero values. Unstructured sparse signals can be modeled well by an i.i.d. Bernoulli sparse prior,

where $\psi(x_{i})$ can be any distribution, e.g. $\psi(x_{i})=\mathcal{N}(x_{i};\bar{x},\sigma^{2})$, and the degree of sparsity is controlled by the value $\rho\in$. Notice that, in this usual setting, both distributions are factorized, that is, the likelihood is in $M$ terms relative to the constraint over each $y_{\mu}$, and the prior is in $N$ terms relative to what is expected of each $x_{i}$. Factorized distributions such as these are well represented by graphical models , specifically, bipartite graphs in which the $M+N$ factors are represented by one type of node and the $N$ variables $x_{i}$ by another. Once the posterior distribution is written down, the estimate $\hat{\operatorname{\mathbf{x}}}$ may be assigned in different ways, according to what loss function one wishes to minimize. In this work, we are chiefly concerned with the minimum mean-squared error (MMSE) estimate, which can be shown to be the average of $x_{i}$ with respect to the posterior $P(\operatorname{\mathbf{x}}|\Phi,\operatorname{\mathbf{y}})$; if one were able to compute the posterior’s marginals, the MMSE estimate would read

The strategy employed by AMP is to infer the marginals of the posterior by using a relaxed version of the BP algorithm , and thus to arrive at the MMSE estimate of the unknown signal $\operatorname{\mathbf{x}}$.

II-B Relaxed Belief-Propagation

BP implements a message-passing scheme between nodes in a graphical model, ultimately allowing one to compute approximations of the posterior marginals. Messages $m_{i\to\mu}$ are sent from the variables nodes to the factor nodes and subsequent messages $m_{\mu\to i}$ are sent from factor nodes back to variable nodes that corresponds to algorithm’s current “beliefs” about the probabilistic distribution of the variables $x_{i}$. Since these distributions are continuous, the first relaxation step is to move to a projected version of these distributions, as described in . Here, we shall follow the notation of reference and use the following parametrization:

This leads (see ) to the following closed recursion sometimes called relaxed BP (r-BP):

where the functions $f$ are defined by the following prior-dependent integrals

After convergence, the single point marginals are given by

We intentionally write r-BP without specifying time indices since the updates can be performed in one of two ways. The first approach is to update in parallel, where all variables are updated at time $t$ given the state at time $t-1$. The second is the random sequential update where one picks a single index $i$ and updates all messages corresponding to it. A time-step is completed once all indices have been visited and updated once. As shown in , the sequential, or swept, iteration is much more stable for r-BP. We now turn our attention to AMP and to our proposed modification.

II-C Swept Approximate Message Passing

In the message-passing described in the previous section, $2(M\times N)$ messages are sent, one between each variable component and each measurement at each iteration. This creates a very large computational and memory burden for applications with large $N$, $M$. It is possible to rewrite the BP equations in terms of only $N+M$ messages by making the assumption that $\Phi$ is dense and that its elements are of magnitude $O(1/\sqrt{N})$. In statistical physics, this assumption leads to the TAP equations used in the study of spin glasses. For graphical models, such strategies have been discussed in . The use of TAP with r-BP provides the standard AMP iteration, as we now show. First we define

Next we expand around the marginals and disregard any $O(1)$ terms (see for details) to find:

Now let us investigate the expansion of the factor $\omega_{\mu}$ as we include the time, or iteration, indices $t$. First one has

which allows us to close the equations on the set of $a,v,R,\Sigma,V$ and $\omega$. Iterating all relations in parallel (i.e. updating all $R,\Sigma$’s, then $a,v$’s and then the $\omega,V$’s) provides the AMP iteration.

The implementation of the sequential update is not a straightforward task as many otherwise intuitive attempts lead to non-convergent algorithms. The key observation in the derivation of SwAMP is that (18) mixes different time indices: while the “$a$” and “$V$” are the “new ones”, the expression in the fraction is the “old” one, i.e. the one before the last iteration. The implication of this is that while $\sum_{i}\Phi_{\mu i}a_{i}$ and $V_{\mu}$ should be recalculated as the updates sweep over $i$ at a single time-step, the term $(y_{\mu}-\omega_{\mu})/(\Delta_{\mu}+V_{\mu})$ (which we denote as $g_{\mu}$ later on) should not. A corresponding bookkeeping then leads to the SwAMP algorithm for the evolution of $\omega_{\mu}$, $\Sigma^{2}_{i}$, $V_{\mu}$ and $R_{i}$ described in Alg. 1. At this point, the difference between AMP and SwAMP appears minimal, but, as we shall see, the differences in convergence properties turn out to be spectacular.

Finally, we note that this procedure can also be generalized, a la GAMP, for output channels other than the AWGN. The required change is minimal : one should replace the term $(y_{\mu}-\omega_{\mu})/(\Delta_{\mu}+V_{\mu})$ in the $R_{i}$ and $\omega_{\mu}$ updates with $g_{\text{out}}\left(\omega_{\mu},V_{\mu}\right)$, a generic function which depends on the channel. Specifically, $g_{\text{out}}(\omega,V)=\int dz\,P(y|z)\,e^{-\frac{(z-\omega)^{2}}{2V}}\left(\frac{z-\omega}{V}\right)$. Additionally, the $\frac{1}{\Delta_{\mu}+V_{\mu}}$ term in the $\Sigma_{i}^{2}$ update should be replaced by $-\frac{\partial g_{\text{out}}}{\partial\omega}$. Notice that all AWGN specific terms are recovered for $P(y|z)\propto e^{-\frac{(y-z)^{2}}{2\Delta}}$.

III Numerical Results

As discussed earlier, using projections of non-zero mean to sample $\operatorname{\mathbf{x}}$ is one of the simplest cases for which AMP can fail to converge. However, by using the proposed SwAMP approach, accurate estimates of $\operatorname{\mathbf{x}}$ can be obtained even when the mean of the projections is non-negligible. While it may be possible to use mean subtraction, our proposed approach does not require such preprocessing. Additionally, as we will show later, not all problems are amenable to such mean subtraction. To evaluate the effectiveness of SwAMP as compared to the standard parallel-update AMP iteration, we draw i.i.d. projections according to

We also considered an even more troublesome case for projections, namely, a set of projections which are strongly correlated. For these tests, we draw

III-B Group Testing

Group testing, also known as pooling in molecular biology, is an approach to designing experiments so as to reduce the number of tests required to identify rare events or faulty items. In the most naive approach to this problem, the number of tests is equal to the number of items, as each item is tested individually. However, since only a small fraction of the items may be faulty, the number of tests can be significantly reduced via pooling, i.e. testing many items simultaneously and allowing items to be included within multiple different tests. The nature of this linear combination of tests allows for a CS-type approach to faulty item detection, but with a few important caveats. First, the operator is extremely sparse since the number of pools, and the number of items in them, may be limited due to physical testing constraints. Second, the elements of this operator are commonly $0/1$. Group testing is therefore a very challenging application for AMP since the properties of the group testing operator do not match AMP’s assumptions.

In one recent work , the authors use both BP and AMP for group testing and found that while basic AMP would not converge, very good results—optimal ones, in fact—could be obtained by using a BP approach. This came at a large computational cost, however. Here, we have repeated the experiment of using the SwAMP approach instead of AMP and BP. In fact, for SwAMP, a sparse operator is a very advantageous situation in terms of computational efficiency. Since the projector is extremely sparse by construction, we may explicitly ignore operations involving null elements, thus considerably improving the algorithm’s speed, as seen in Fig. 2(b). Here, we also see that SwAMP’s computational complexity is on the order of $O(N^{2})$, as is AMP’s. Group testing experiments are shown in Fig. 2(a) where we use random $0/1$ projections, under the constraint that each projection should sum to $7$, to sample sparse $0/1$ signals with $K\ll N$ non-zero elements, where $N$ is the signal dimensionality. While AMP diverges when attempting to recover these signals, SwAMP converges to the correct solution in few iterations. Additionally, SwAMP very closely matches the BP transition, thus providing recovery performance better than convex optimization, just as BP does, but with much less computational complexity.

III-C 1-bit Compressed Sensing

One of the confounding factors regarding the practical implementation of CS in hardware devices is the treatment of measurement quantization. The original CS analysis provides recovery bounds based upon the assumption of real-valued measurements. However, in practice, hardware devices cannot capture such values with infinite precision, and so some kind of quantization on the measurements must be implemented. Specifically, if $Q(\cdot)$ is a uniform scalar quantizer, then $\operatorname{\mathbf{y}}=Q(\Phi\operatorname{\mathbf{x}},B),$ where $B$ is the number of bits used to represent the measurement. If signal recoverability is significantly impacted by small $B$, then the dimensionality reduction provided by CS may be lost by the requirement for many bits to encode each measurement.

Thankfully, recent works have shown CS recovery to be robust to quantization and the non-linear error it introduces. In fact, CS has been shown to be robust even in the extreme case $B=1$ known as 1-bit CS. In this case, the quantized measurements are given by

The non-linearity and severity of 1-bit CS requires special treatment from the CS recovery procedure. In , a renormalized fixed-point continuation (RFPC) algorithm was proposed. Later, analyzed the sensitivity of 1-bit CS to sign flips and proposed a noise-robust recovery algorithm, binary iterative hard thresholding (BIHT).

Both methods and show the effectiveness of algorithms grounded in statistical mechanics for quantized CS reconstruction. However, both assume an amenable set of projections. Even projections possessing small mean can cause large degradations in performance. While mean removal is occasionally effective in the usual CS setting, it cannot be used for 1-bit CS due to the nature of the sign operation in (21). An algorithm that can handle troublesome projectors can therefore be of great use. In Sec. II-C, we show how the SwAMP can be modified to the general-channel setting, as was done in GAMP. This generalization allows for 1-bit CS recovery with SwAMP under much more relaxed requirements for $\Phi$.

IV Conclusion

Exact analysis of the asymptotic state evolution of SwAMP, as well as a thorough analytical proof of its convergence, remains a challenging open problem for future work.

V Acknowledgments

This work has been supported in part by the ERC under the European Union’s 7th Framework Programme Grant Agreement 307087-SPARCS, by the Grant DySpaN of “Triangle de la Physique,” and by FAPESP under grant 13/01213-8.

References