On Convergence of Approximate Message Passing

I Introduction

Approximate message passing is an algorithm derived from belief propagation that has been recently used with success in a number of sparse estimation problems, see e.g. . Highly non-trivial theoretical results were obtained on the performances of this algorithm . Based on these developments and the promising nature of their results we can anticipate that AMP based algorithms will become the state-of-the-art algorithms for many problems of practical interest.

Just as with any iterative algorithm the main question about AMP, besides its performance, is its convergence. This question is largely open except for the case of compressed sensing, i.e. estimation of a sparse ${\bf x}$ from noisy linear projections

with matrices $F$ having iid entries of zero mean, and ${\xi}$ a white Gaussian noise of variance $\Delta$. This last case has been treated in the rigorous proofs in the very large signal size limit of . However, for many other sparse estimation problems, or for slightly more general matrices $F$, the basic version of AMP fails to converge (and worst, can diverge violently). Attempts to fix these convergence issues were so far limited to rather basic and empirical strategies such as damping the iterations in various ways, or transforming the matrix by subtracting its mean. Such strategies are rarely discussed in the literature and often appear only in the associated implementations available online. Moreover, they are far from ensuring the convergence in all cases and some of these strategies (e.g. the mean removal) are not usable in more challenging signal processing settings where approximate message passing can be applied (e.g. the dictionary learning problem ). The main motivation of this work is to understand the origin of some of these convergence problems.

The structurally simplest case where AMP fails to converge appears to be when the measurement matrix $F$ has iid entries of non-zero mean. This problem was noticed by several authors, e.g. , and fixed in the implementations by removing the mean of the matrix. Indeed, the average of element of the measurement vector ${\bf y}$ reads

We denote $\overline{F}_{i}=\sum_{\mu}F_{\mu i}/M$ the average value of $F$ for column $i$. One can then work with the modified system $y_{\mu}-\overline{y}=\sum_{i}(F_{\mu i}-\overline{F}_{i})x_{i}$ where the mean of the new sensing matrix $F_{\mu i}-\overline{F}_{i}$ is zero. A similar (but different) trick is used in the implementation of . This ”remove mean” strategy is, however, not fully satisfactory because it is not understood why it is needed in the first place, nor under what conditions it restores the convergence. Moreover in some more general settings it is not applicable at all.

The goal of this paper is to analyze the origin of the non-convergence for non-zero mean matrices and discuss general strategies to prevent it. Such an understanding is a step towards the design of robustly convergent and hence more efficient AMP-based algorithms. We will hence consider matrices with entries generated as follows

For $\gamma=0$ this is the case that has been considered in the literature. To be specific and simple we will consider that the signal ${\bf x}$ was generated to have $\rho N$ non-zero entries that are iid normally distributed with zero mean and unit variance

We will consider the Bayesian version of the AMP algorithm that uses this prior information about the signal. A first observation is that AMP does not depend on $\gamma$ in an explicit way: this can be checked explicitly by repeating the detailed derivations of AMP present in the literature for $\gamma>0$ (follow e.g. the derivation in ).

On the other hand the asymptotic analysis of the performance of the algorithm — the state evolution — depends on $\gamma$ explicitly and hence we have to rederive it. The analysis of the state evolution for $\gamma>0$ will lead to an understanding of the origin of the convergence problems.

II The AMP algorithm

We consider the AMP algorithm in the form that was derived in . The main steps are a) going from belief propagation (BP) to a relaxed BP (r-BP) where only the two first moments of all messages are kept and b) using $N$ sites marginals instead of $N\times M$ messages and adding the compensating Onsager terms . Finally, AMP reads:

where $f_{k}(\Sigma^{2},R)$, here and in what follows, are the $k$-th connected cumulants w.r.t. the probability measure

where $Z(\Sigma^{2},R)$ is the normalization constant.

The variables $a_{i}$ and $v_{i}$ are the AMP estimators for the mean and variance of the component $i$ of the signal. The quality of the reconstruction can be evaluated by computing the mean squared error (MSE)

When $\gamma\!=\!0$, the performance of the AMP algorithm was analyzed rigorously in the limit of large system size via the state evolution $(E^{t+1},V^{t+1})=G(E^{t},V^{t})$, where $G$ is a function specified in . An important property of the Bayes optimal inference (i.e. when the signal was indeed generated from the assumed prior distribution) is that the two paramaters are equal in the large size limit, $E^{t}\!=\!V^{t}$, and the state evolution hence reduces to an iterative equation of a single real number, which is amenable to rigorous analysis . In statistical physics $E^{t}\!=\!V^{t}$ is called the Nishimori condition and is discussed in the context of compressed sensing in detail in . In general, when $\gamma\!=\!0$ we observed by analyzing the state evolution equations that even when at initial times $E^{t=0}\!\neq\!V^{t=0}$ the equality $E^{t}\!=\!V^{t}$ is restored after a sufficient number of iterations.

III State evolution with non-zero mean matrices

The state evolution of the AMP algorithm can be derived for measurement matrices with non-zero mean $\gamma>0$. Here we follow closely the derivation and notation from for zero mean matrices. Among the different variables, the statistical distribution of $R_{i}$ plays a crucial role in the determination of the state evolution it can be written as

where $s_{i}$ is the original signal component and

is a Gaussian random variable, and $a_{j\rightarrow\mu}^{t}$ is an auxiliary variable related closely to $a_{i}^{t}$ that appears in the derivation of the AMP algorithm. Assumptions used to derive AMP can be used to compute the mean and variance of $r_{i}^{t}$ over realizations of the problem. In the leading order we get

where we have defined a new order parameter

The parameter $D^{t}$ is not needed for zero mean matrices $\gamma=0$. For $\gamma>0$, however, the state evolution is written in terms of three parameters $E^{t}$, $V^{t}$ and $D^{t}$. The remaining steps in the derivation are basically identical to those for zero mean matrices and following we obtain

where ${\cal D}z$ is a Gaussian measure and

When the mean of the measurement matrix is zero, $\gamma=0$, these equations clearly reduce to those derived in .

Also for $\gamma>0$ we can identify the Nishimori condition, which reads $E^{t}=V^{t}$ (for the same reasons as for the previous case) and $D^{t}=0$ (since under Bayes optimal inference the mean of the estimator must be equal to the true mean of the signal). It is a question of simple algebraic verification to see that starting with $E^{t}=V^{t}$ and $D^{t}=0$ eqs. (19-21) lead to $E^{t+1}=V^{t+1}$ and $D^{t+1}=0$. Hence if we restrict ourselves to the space on which the Nishimori conditions hold (called the Nishimori line) there is no difference between the $\gamma=0$ and $\gamma>0$ case.

IV Instability of the Nishimori line

In this Section we analyze the dynamical stability of the Nishimori line (NL) under iterations of eqs. (19-21). We consider the space $(K,D)$ orthogonal to the NL, where $K=V-E$. We know that in this space $(K^{*}=0,D^{*}=0)$ is a fixed point. We can generically write

To analyze the stability we linearize around the fixed point considering the perturbations $\delta K^{t}=K^{t}-K^{*}$ and $\delta D^{t}=D^{t}-D^{*}$. The linearized formula reads

It follows from a straightforward algebraic verification that both the off-diagonal terms (the cross derivatives) are zero for the distribution $P(x)$ from eq. (4). The matrix ${\cal M}$ (25) is hence diagonal. For a more generic prior distribution the situation is slightly more involved, but qualitatively analogous to the one of (4). The diagonal terms read

where, as before, the functions $f_{k}(\Sigma^{2},R)$ are the $k$-th connected cumulants with respect to the measure ${\cal Q}(\Sigma^{2},R)$ (11), and where we denoted

The term $\partial_{K}f_{K}(V^{t})$ is independent of $\gamma$ and its module is always smaller than one. Hence the Nishimori line is stable in the direction $K=V-E$.

On the other hand the term $\lambda_{D}=\partial_{D}f_{D}(V^{t})$ has a non-trivial behavior that we illustrate in Fig. 1 for $\rho=0.1$, $\alpha=0.3$, $\Delta=10^{-10}$ and, respectively, $\gamma=1.9$, $\gamma=2.5$, $\gamma=2.9$ and $\gamma=3.6$. In the figure we identify three different regimes:

For $|\gamma|<\gamma_{c}^{(1)}$ the eigenvalue $\lambda_{D}$ is always less than $1$ in modulus.

For $\gamma_{c}^{(1)}<|\gamma|<\gamma_{c}^{(2)}$ the eigenvalue becomes greater than $1$ in modulus in a certain portion of the Nishimori line. In this region the evolution tends to make $|D|$ larger, while at the same time $V$ and $E$ decrease.

For $|\gamma|>\gamma_{c}^{(2)}$ the eigenvalue $\lambda_{D}$ is larger than $1$ in modulus in the whole range down to the fixed point. This tells us that any fluctuation of $D$ will be progressively enhanced.

We further realize that the expression used to calculate $\lambda_{D}$ depends on the value $V$ only through the variable $A$ (28) and not in an explicit way on the parameters $\alpha$ and $\Delta$. This means that the threshold value $\gamma_{c}^{(1)}$ is from its definition independent of $\alpha$ and $\Delta$. The threshold value $\gamma_{c}^{(2)}$ is also independent of $\alpha$ for $\Delta=0$ and only weakly dependent on both $\alpha$ and $\Delta$ for small values of $\Delta$. In Fig. 2 we hence plot the two threshold values for $\Delta=0$ (in which case they are both independent of the undersampling $\alpha$) as a function of the sparsity $\rho$.

V Comparing state evolution to AMP

We now discuss how does the instability of the Nishimori line translate into the behavior of the state evolution (SE) initialized usually as $E^{t=0}=V^{t=0}=\rho$ (corresponding to $a^{t=0}_{i}=0$ and $v^{t=0}_{i}=\rho$) and $D^{t=0}=0$. The SE was derived to correspond to the behavior of the AMP algorithm for sufficiently large system sizes $N$. We observe that

For $|\gamma|<\gamma_{c}^{(1)}$ the SE converges to the fixed point with monotonically decreasing $E=V$. There are really infinitesimal fluctuations in $D$ that are due to numerical precision but they are harmless.

For $\gamma_{c}^{(1)}<|\gamma|<\gamma_{c}^{(2)}$ the SE converges to the fixed point with monotonically decreasing $E=V$. In the region of $V$ in which $|\lambda_{D}|$ is larger than one, we observe that the numerical fluctuations of $D$ are slightly increased (especially if we are close to $\gamma_{c}^{(2)}$), without changing qualitatively the behavior of $V$, and when $|\lambda_{D}|$ becomes again smaller than 1 the fluctuations are reabsorbed.

For $|\gamma|>\gamma_{c}^{(2)}$ the fluctuations of $D$ are increased along the whole line $E=V$. At some point these fluctuations reach so large values that the difference $K=E-V$ grows and we observe a divergence of both $E$ and $V$.

Therefore, while with infinite numerical precision the SE should stay on the Nishimori line and converge whatever the value of $\gamma$ is, from the practical point of view the fluctuations due to numerical precision are sufficient to cause divergence in the third regime. Of course in the AMP algorithm the typical fluctuations are of order $1/\sqrt{N}$ hence relatively large and that is the reason why for $|\gamma|>\gamma_{c}^{(2)}$ AMP never converges. In fact these finite size fluctuations are so strong that even in the second regime $\gamma_{c}^{(1)}<|\gamma|<\gamma_{c}^{(2)}$ AMP might have problems. Therefore we observe a smooth transition in the success rate $R=$($\#$successes/$\#$failures) between $\gamma_{c}^{(1)}$ and $\gamma_{c}^{(2)}$ for finite $N$. When $N$ is increased this smooth transition becomes sharper. In the inset of Fig. 2 we show the success rate of the AMP algorithm averaged over $1000$ random instances of the measurement matrix for $N=1000,4000$ and over $500$ instances for $N=16000$. We see that, even if asymptotically, the reference value for the success/failure transition would be $\gamma_{c}^{(2)}$, for all practical system sizes, the right threshold to look at is rather $\gamma_{c}^{(1)}$.

VI Reducing the instability

There are at least two strategies that appear in the implementations of the AMP algorithm that improve its convergence. Let us discuss them now in the context of the above analysis.

A popular and generic strategy to improve convergence of iterative algorithms is “damping”, i.e. in every new iteration we update the variables only partially. Such damping (with different schemes) appears in basically every available implementation of AMP. In the view of the preceding analysis a dynamical instability is mitigated by such damping and the eigenvalue $|\lambda_{D}|$ becomes effectively smaller. Indeed AMP with damping converges well even for matrices with means slightly larger than those corresponding to $\gamma_{c}^{(2)}$ in Fig. 2.

Expectation maximization learning

In this paper so far we assumed the prior knowledge of the probability distribution of signal elements as well as of the measurement noise $\Delta$ and the sparsity $\rho$. A classical strategy of expectation maximization was suggested, tested and implemented in in order to learn these parameters when they are not known apriori. A careful investigation of the AMP algorithm with EM learning leads to a conclusion that with the learning the AMP has better convergence properties than without.

This can come as a surprise at first, but in the view of our above investigation it can now be easily explained. The EM update in a sense imposes (in an iterative way) the Nishimori condition, see the derivation of EM in , hence it should be expected that it also stabilizes the Nishimori line and consequently improves the convergence of AMP.

VII The sequential redemption

AMP being so sensitive to the mean of the matrix elements is surprising because the standard BP, when applied to discrete random problems, does not experience such problems. In this last section we argue that the convergence problems in the case of CS with non-zero mean measurement matrices are actually specific to the “parallel updates” (involving only matrix multiplications) performed naturally in the AMP algorithm that we presented in Sec. II. Let us recall the so-called relaxed-BP (r-BP) algorithm (for present notations see ) where messages are sent on the factor graph:

We intentionally wrote this algorithm without the time indices, because the update can be performed in two ways. First in the parallel one 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. For r-BP, this leads to the same computational complexity, however, it is important to realize that AMP is actually written assuming the r-BP with the parallel update. In Fig. 3 we compare the behavior of parallel and random sequential r-BP: as we see, the sequential update does not seem to be affected by the non-zero mean.

This observation of the parallel update being more problematic than the sequential one is actually not surprising a posteriori. In fact, such a lack of convergence is known to occur in parallel iterations in many problems due to instabilities just like the one we have studied here (see for instance the “modularity” instability in the hard-core model and coloring problems on random graphs). Using instead, when possible, a sequential r-BP update is therefore an interesting alternative. Nevertheless, it is not a universal solution since it by no means guarantees convergence for all matrices. Also, the disadvantage of the sequential r-BP update is that it looses the nice property of only involving matrix multiplication, a crucial property for scalability for operators, such as the fast Fourier transform, for which there exist efficient multiplication methods.

VIII Conclusions

We have analyzed the convergence problems of AMP in the specific case of compressed sensing with measurement matrices having iid entries of non-zero mean. Despite the fact that the AMP iterations are not modified w.r.t. the case of zero mean, the state evolution does contain an additional order parameter. The main result of the paper, contained in Sec. IV, is that the presence of this third parameter causes instabilities of the so-called Nishimori line and, therefore in the algorithm itself, if the mean of the matrix elements exceeds some critical value. In the last section we show that the convergence issue for matrices of non-zero mean are strongly mitigated when random sequential update is used in the message passing instead of the parallel one that is standard to AMP.

This analysis represents a step towards understanding of the nature of convergence issues in message passing algorithms that are ubiquitous in problems ranging from physics to information theory. More complete understanding of these issues is needed before message passing algorithms can become part of standard toolbox to solve a wide range of problems of practical interest.

Acknowledgment

This work has been supported by the ERC under the European Union’s 7th Framework Programme Grant Agreement 307087-SPARCS, and by the project TASC of the Labex PALM.

References