On the Convergence of Approximate Message Passing with Arbitrary Matrices

I Introduction

for separable $F(\mathbf{z})=\sum_{i}F_{i}(z_{i})$ and $G(\mathbf{x})=\sum_{j}G_{j}(x_{j})$. Such problems arise in a range of applications including statistical regression, inverse problems, and compressed sensing.

Most current numerical methods for solving the constrained optimization problem (2) attempt to exploit the separable structure of the objective function (2) using approaches like iterative shrinkage and thresholding (ISTA) , the alternating direction method of multipliers (ADMM) , or primal-dual approaches .

In recent years, however, there has also been considerable interest in approximate message passing (AMP) methods that apply Gaussian and quadratic approximations to loopy belief propagation (BP) in graphical models . AMP applied to max-sum loopy BP produces a sequence of estimates that approximate $\widehat{\mathbf{x}}_{\text{\sf MAP}}$, while AMP applied to sum-product loopy BP produces a sequence of estimates that approximate $\widehat{\mathbf{x}}_{\text{\sf MMSE}}$. For zero-mean i.i.d. sub-Gaussian $\mathbf{A}$ in the large-system limit (i.e., $m,n\rightarrow\infty$ with fixed $m/n$), AMP methods are characterized by a state evolution whose fixed points, when unique, coincide with $\widehat{\mathbf{x}}_{\text{\sf MAP}}$ or $\widehat{\mathbf{x}}_{\text{\sf MMSE}}$ . In addition, for large but finite-sized i.i.d. Gaussian matrices, recent work shows that AMP is close to Bayes-optimal.

Unfortunately, a rigorous characterization of AMP for generic $\mathbf{A}$ remains lacking. The recent papers studied the fixed-points of the generalized AMP (GAMP) algorithm from for generic $\mathbf{A}$. In , it was established that the fixed points of max-sum GAMP coincide with the critical points of the optimization objective in (2). Similarly, established that the fixed points of sum-product GAMP are critical points of a large-system version of the Bethe free energy from . However, the papers did not discuss the convergence of the algorithm to those fixed points. Indeed, similar to other loopy BP algorithms, GAMP may diverge, as demonstrated for mildly ill-conditioned $\mathbf{A}$ in . Likewise, showed that AMP can diverge with non-zero-mean i.i.d. Gaussian $\mathbf{A}$ and the divergence can, in fact, be predicted via a state-evolution analysis.

For general loopy BP, a variety of methods have been proposed to improve convergence, including coordinate descent, tree re-weighting, and double loop methods . In this paper, we propose and analyze a “damped” modification of GAMP that is similar to the technique used in Gaussian belief propagation —a closely related algorithm. We also point out connections between damped GAMP and the primal-dual hybrid-gradient (PDHG) algorithm popular in convex optimization. This connection enhances the interpretability of AMP methods, especially for those who are less familiar with belief propagation.

Our first main result establishes a necessary and sufficient condition on the global convergence of damped GAMP for arbitrary $\mathbf{A}$ in the special case of Gaussian $P(x_{j})$ and $P(y_{i}|z_{i})$ (i.e., quadratic $F$ and $G$) and fixed scalar stepsizes. This condition (see Theorem 2 below) shows that, with sufficient damping, the Gaussian GAMP algorithm can be guaranteed to converge. However, the amount of damping grows with the peak-to-average ratio of the squared singular values of $\mathbf{A}$. This result explains why Gaussian GAMP converges (with high probability) for large i.i.d. Gaussian $\mathbf{A}$, but it also explains why it needs to be damped significantly for non-zero-mean, low-rank, or otherwise ill-conditioned $\mathbf{A}$.

Our second result establishes the local convergence of GAMP for strictly convex $F$ and $G$ and arbitrary, but fixed, vector-valued stepsizes. This sufficient condition is similar to the Gaussian case, but involves a certain row-column normalized version of $\mathbf{A}$. (See Theorem 3 below.)

Finally, we present numerical experiments that verify the tightness of the sufficient conditions from Theorems 2 and 3.

II Damped GAMP

The GAMP algorithm was introduced in and rigorously analyzed in . The procedure (see Algorithm 1) produces a sequence of estimates $\widehat{\mathbf{x}}^{t},t=1,2,\dots$, that, in max-sum mode, approximate $\widehat{\mathbf{x}}_{\text{\sf MAP}}$ and, in sum-product mode, approximate $\widehat{\mathbf{x}}_{\text{\sf MMSE}}$. The two modes differ only in the definition of the scalar estimation functions ${g}_{s}$ and $g_{x}$ used in lines 8, 9, 12, and 13 of Algorithm 1:

using \bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{r}=[\tau_{r_{1}},\dots,\tau_{r_{n}}]^{\text{\sf T}}, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{p}=[\nu_{p_{1}},\dots,\nu_{p_{m}}]^{\text{\sf T}}, and

Note (3) implements scalar MAP denoising under prior $P(x_{j})\!\propto\!\exp(-G(x_{j}))$ and variance-$\tau_{r_{j}}$ Gaussian noise.

and so (6) is the scalar MMSE denoiser under $P(x_{j})\!\propto\!\exp(-G(x_{j}))$ and variance-$\tau_{r_{j}}$ Gaussian noise.

Note that, in Algorithm 1 and the sequel, $\mathbf{a}.\mathbf{b}$ and $\mathbf{a}./\mathbf{b}$ denote component-wise multiplication and division, respectively, between vectors $\mathbf{a}$ and $\mathbf{b}$.

Algorithm 1 reveals the computational efficiency of GAMP: the vector-valued MAP and MMSE estimation problems are reduced to a sequence of scalar estimation problems in Gaussian noise. Specifically, each iteration involves multiplications by $\mathbf{S}$, $\mathbf{S}^{\text{\sf T}}$, $\mathbf{A}$ and $\mathbf{A}^{\text{\sf H}}$ along with simple scalar estimations on the components $x_{j}$ and $z_{i}$; there are no vector-valued estimations or matrix inverses.

We note that Algorithm 1 writes GAMP in a “symmetrized” form, where the steps in lines 6-9 mirror those in lines 10-13. This differs from the way that GAMP is presented in most other publications, such as , which is obtained by replacing the variables $\mathbf{s}$, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{p}, and $\mathbf{p}$ in Algorithm 1 by $-\mathbf{s}$, \mathbf{1}./\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{p}, and \mathbf{p}.\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{p}, respectively. Note that, thoughout this paper, we use $\tau$ for variance quantities and $\nu$ for precision (i.e., inverse variance) quantities.

II-B Damped GAMP

Algorithm 1 includes a small but important modification to the original GAMP from : lines 9 and 13 perform damping using constants $\theta_{s},\theta_{x}\in(0,1]$ that slow the updates of $\mathbf{s}^{t},\mathbf{x}^{t}$ when $\theta_{s},\theta_{x}<1$, respectively. The original GAMP implicitly uses $\theta_{s}=1=\theta_{x}$. In the sequel, we establish—analytically—that damping facilitates the convergence of GAMP for general $\mathbf{A}$, a fact that has been empirically observed in past works (e.g., ).

II-C GAMP with Scalar Stepsizes

The computational complexity of Algorithm 1 is dominated by the matrix-vector multiplications involving $\mathbf{A}$, $\mathbf{A}^{\text{\sf H}}$, $\mathbf{S},$ and $\mathbf{S}^{\text{\sf T}}$. In , a scalar-stepsize simplification of GAMP was proposed to avoid the multiplications by $\mathbf{S}$ and $\mathbf{S}^{\text{\sf T}}$, roughly halving the per-iteration complexity. The meaning of “stepsize” will become clear in the sequel. Algorithm 2 shows the scalar-stepsize version of Algorithm 1.

For use in the sequel, we now show that scalar-stepsize GAMP is equivalent to vector-stepsize GAMP under a different choice of $\mathbf{S}$. While Algorithm 1 uses $\mathbf{S}=\mathbf{A}.\mathbf{A}$, Algorithm 2 effectively uses

i.e., a constant matrix having the same average value as $\mathbf{A}.\mathbf{A}$. Thus, the two algorithms coincide when $|A_{ij}|$ is invariant to $i$ and $j$. To see the equivalence, we first note that, under $\mathbf{S}$ from (7), line 6 in Algorithm 1 would produce a version of 1/{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{p}^{t} containing identical elements $1/\nu_{p}^{t}$, where

for \tau_{x}^{t}=(1/n)\mathbf{1}^{\text{\sf T}}\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x}. Similarly, line 10 would produce a vector 1/\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{r}^{t} with identical elements $1/\tau_{r}^{t}$, where

for \nu_{s}^{t}=(1/m)\mathbf{1}^{\text{\sf T}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{s}^{t}. Furthermore, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{p}^{t}=\nu_{p}^{t}\mathbf{1} and line 8 imply that \nu_{s}^{t}=(\nu_{p}^{t}/m)\mathbf{1}^{\text{\sf T}}{g}_{s}({\mathbf{p}}^{t},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{p}^{t}), while \bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{r}^{t}=\tau_{r}^{t}\mathbf{1} and line 12 imply that \tau_{x}^{t\!+\!1}=(\tau_{r}^{t}/n)\mathbf{1}^{\text{\sf T}}g_{x}(\mathbf{r}^{t},\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{r}^{t}). Applying these modifications to Algorithm 1, we arrive at Algorithm 2.

II-D Relation to Primal-Dual Hybrid Gradient Algorithms

An important case of (2) is when $F$ and $G$ are closed proper convex functionals and the solution $\widehat{\mathbf{x}}_{\text{\sf MAP}}$ exists. Recently, there has been great interest in solving this problem from the primal-dual perspective , which can be described as follows. Consider $F^{*}$, the convex conjugate of $F$, as given by the Legendre-Fenchel transform

For closed proper convex $F$, we have $F^{**}=F$, and so

which gives the equivalent saddle-point formulation of (2),

The so-called primal-dual hybrid-gradient (PDHG) algorithm recently studied in is defined by the iteration

where $\theta\in$ is a relaxation parameter. Line (11) can be recognized as proximal gradient ascent in the dual variable ${\mathbf{s}}$ using stepsize $\nu_{p}$, while line (12) is proximal gradient descent in the primal variable $\mathbf{x}$ using stepsize $\tau_{r}$.

PDHG can be related to damped scalar-stepsize GAMP as follows. Since $F$ is proper, closed, and convex, we can apply the Moreau identity

to (4), after which the assumed separability of $F$ implies that

Thus, under $\theta_{s}=1$, scalar GAMP’s update of ${\mathbf{s}}$ (in line 8 of Algorithm 2) matches PDHG’s in (11). Similarly, noting the connection between (3) and (12), it follows that, under $\theta_{x}=1$, scalar GAMP’s update of $\mathbf{x}$ (in line 12 of Algorithm 2)) matches the PDHG update (13) under $\theta=0$.

In summary, PDHG under $\theta=0$ (the Arrow-Hurwicz case) would be equivalent to non-damped scalar GAMP if the stepsizes $\nu_{p}^{t}$ and $\tau_{r}^{t}$ were fixed over the iterations. GAMP, however, adapts these stepsizes. In fact, under the existence of the second derivative $f^{\prime\prime}$, it can be shown that

implying that, for smooth $F$ and $G$, GAMP updates $\tau_{x}^{t}$ according to the average local curvature of $G$ at the point $\mathbf{x}=\operatorname{prox}_{\tau_{r}^{t}G}(\mathbf{r}^{t})$ and updates $\nu_{s}^{t}$ according to the average local curvature of $F^{*}$ at the point ${\mathbf{s}}=\operatorname{prox}_{\nu_{p}^{t}F^{*}}({\mathbf{p}}^{t})$. A different form of PDHG stepsize adaptation has been recently considered in , one that is not curvature based.

Meanwhile, PDHG under $\theta\neq 0$ is similar to fixed-stepsize damped scalar GAMP with $\theta_{s}=1$ and $\theta_{x}=1+\theta$, although not the same. Note that PDHG uses the damped version of $\mathbf{x}$ only in the dual update (11) whereas GAMP uses the damped version of $\mathbf{x}$ in both primal and dual updates. Also, PDHG relaxes only the primal variable $\mathbf{x}$, whereas damped GAMP relaxes (or damps) both primal and dual variables.

III Damped Gaussian GAMP

Although Algorithms 1 and 2 apply to generic distributions $P(x_{j})$ and $P(y_{i}|z_{i})$, we find it useful to at first consider the simple case of Gaussian distributions, and in particular

where $\tau_{0_{j}}$ are variances and $\nu_{w_{i}}$ are precisions (i.e., inverse variances). In this case, the scalar estimation functions used in max-sum mode are identical to those in sum-product mode, and are linear :

Henceforth, we use “Gaussian GAMP” (GGAMP) when referring to GAMP under the estimation functions (17).

III-B Convergence of GGAMP Stepsizes

We first establish the convergence of the GGAMP stepsizes in the case of an arbitrary matrix $\mathbf{A}$. For the vector-stepsize case in Algorithm 1, lines 8 and 12 become

and, combining these with lines 6 and 10, we get

which are invariant to $\theta_{s},\theta_{x},{\mathbf{s}}^{t}$, and $\mathbf{x}^{t}$. The scalar-stepsize case in Algorithm 2 is similar, and in either case, the following theorem shows that the GGAMP stepsizes always converge.

Consider Algorithms 1 or 2) with Gaussian estimation functions (17) defined for any vectors {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{w} and \bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{0}>\bm{0}. Then, as $t\rightarrow\infty$, the stepsizes {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{p}^{t},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{s}^{t},\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{r}^{t},\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x}^{t} (or their scalar versions) converge to unique fixed points that are invariant to $\theta_{s}$ and $\theta_{x}$.

IV Scalar-Stepsize GGAMP Convergence

An important special case that we now consider is scalar-stepsize GGAMP from Algorithm 2 under identical variances, i.e.,

for some $\nu_{w}$ and $\tau_{0}>0$. In this case, lines 7 and 11 give

and, combining these with lines 5 and 9, we get

IV-B Convergence

We now investigate the convergence of the primal and dual variables $\mathbf{x}^{t}$ and ${\mathbf{s}}^{t}$ for scalar GGAMP. Since, for this algorithm, the previous section established that, as $t\rightarrow\infty$, the stepsizes $\nu_{p}^{t}$ and $\tau_{r}^{t}$ converge independently of $\theta_{s},\theta_{x},{\mathbf{s}}^{t}$, and $\overline{\mathbf{x}}^{t}$, we henceforth consider GGAMP with fixed stepsizes $\nu_{p}^{t}=\nu_{p}$ and $\tau_{r}^{t}=\tau_{r}$, where $\nu_{p}$ and $\tau_{r}$ are the fixed points of (22) for Algorithm 2. (A generalization to arbitrary fixed stepsizes will be given in Section V.)

Under Gaussian priors (i.e., (17)) with identical variances (20), scalar-stepsize GAMP from Algorithm 2 converges for any $\nu_{w}$ and $\tau_{0}>0$ when

Conversely, it diverges for large enough $\tau_{0}\nu_{w}$ when

Theorem 2 provides a simple necessary and sufficient condition on the convergence of scalar GGAMP. To better interpret this condition, recall that $\|\mathbf{A}\|_{2}^{2}$ is the maximum squared singular value of $\mathbf{A}$ and that $\|\mathbf{A}\|_{F}^{2}$ is the sum of the squared singular values of $\mathbf{A}$ (i.e., $\|\mathbf{A}\|^{2}_{F}=\sum_{i=1}^{\min\{m,n\}}\sigma^{2}_{i}(\mathbf{A})$). Thus

is the peak-to-average ratio of the squared singular values of $\mathbf{A}$. Convergence condition (24) can then be rewritten as

meaning that, for GGAMP convergence, it is necessary and sufficient to choose $\kappa_{\max}(\theta_{s},\theta_{x})$ above the peak-to-average ratio of the squared singular values.

When there is no damping (i.e., $\theta_{s}=1=\theta_{x}$), the definitions in (23) and (27) can be combined to yield

More generally, for $\theta_{s},\theta_{x}\in(0,1]$, it can be shown that

so that the necessary and sufficient GGAMP convergence condition (27) can be rewritten as

which implies that, by choosing sufficiently small damping constants $\theta_{s}$ and $\theta_{x}$, scalar-stepsize GGAMP can always be made to converge.

Condition (30) also helps to understand the effect of $\kappa(\mathbf{A})$ on the GGAMP convergence rate. For example, if we equate $\theta_{s}=\theta_{x}=\theta$ for simplicity, then (30) implies that

Thus, if GGAMP converges at rate $\theta$, then after $\theta$ is adjusted to ensure convergence, GGAMP will converge at a rate below $\sqrt{C/\kappa(\mathbf{A})}$. So larger peak-to-average ratios $\kappa(\mathbf{A})$ will result in slower convergence.

IV-C Examples of Matrices

To illustrate how the level of damping is affected by the nature of the matrix $\mathbf{A}$, we consider several examples.

with equality when $m=n$, and where the approximation becomes exact in the large-system limit. Because this Marcenko-Pastur bound coincides with the $\theta_{s}=1=\theta_{x}$ case (28) of the convergence condition (27), our analysis implies that, for large i.i.d. matrices, scalar stepsize GGAMP will converge without damping, thereby confirming the state evolution analysis. Note that we require that the asymptotic value of $m/n\neq 1$ so that the inequality in (32) is strict; when $m=n$, (32) becomes an equality and we obtain a condition $\Gamma(\theta_{s},\theta_{x})=\|\mathbf{A}\|_{2}^{2}/\|\mathbf{A}\|_{F}^{2}$ right on the boundary between convergence and divergence, where Theorem 2 does not make any statements.

Subsampled unitary matrices

Suppose that $\mathbf{A}$ is constructed by removing either columns or rows, but not both, from a unitary matrix. Then, $\kappa(\mathbf{A})=1$, so, from (29), $\kappa(\mathbf{A})<\kappa_{\max}(\theta_{s},\theta_{x})$ for any $\theta_{s},\theta_{x}\in(0,1]$. Hence, scalar GGAMP will converge with or without damping.

Linear filtering

where $H(e^{j\omega})$ is the DTFT of $\mathbf{h}$. Equation (33) implies that more damping is needed as the filter becomes more narrowband. For example, if $H(e^{j\omega})$ has a normalized bandwidth of $B\in(0,1]$, then $\kappa(\mathbf{A})\approx 1/B$ and, relative to an allpass filter, GGAMP will need to slow by a factor of $O(\sqrt{B})$.

Low-rank matrices

which, from (31), implies the need to choose a damping constant $\theta<\sqrt{Cr/\min\{m,n\}}$, slowing the algorithm by a factor of $\sqrt{\min\{m,n\}/r}$ relative to a full-rank matrix. Hence, more damping is needed as the relative rank decreases.

Walk-summable matrices

Closely related to Gaussian GAMP is Gaussian belief propagation , which performs a similar iterative algorithm to minimize a general quadratic function of the form $f(\mathbf{x})=\mathbf{x}^{\text{\sf H}}\mathbf{J}\mathbf{x}+\mbox{Real}\{\mathbf{c}^{\text{\sf H}}\mathbf{x}\}$ for some positive definite matrix $\mathbf{J}$. Sufficient conditions for the convergence of Gaussian belief propagation were first shown in , but those conditions are difficult to verify. In a now classic result, showed that Gaussian belief propagation will converge when

where $|\mathbf{I}-\mathbf{J}|$ is the component-wise magnitude. The condition (34) is called walk summability, with the constraints $J_{ii}=1$ being for normalization.

A quadratic function $f$ is said to be convex decomposable if it can be written in the form $f(\mathbf{x})=\sum_{i}f_{i}(x_{i})+\sum_{i,j}f_{ij}(x_{i},x_{j})$ where $\{f_{i}\}$ are strictly convex quadratic functions and $\{f_{ij}\}$ are convex quadratic functions. Moallemi and Van Roy showed that if a quadratic objective function is convex decomposable then min-sum message passing converges to the global minimum. In , it was shown that a function is convex decomposable if and only if it is walk-summable (i.e., the two properties are equivalent).

To compare walk summability with GGAMP, first observe that, in the identical-variance case (20), GGAMP performs the same quadratic minimization with a particular $\mathbf{c}$ and with

Now, consider the high-SNR case, where $\tau_{0}=1$ and $\nu_{w}^{-1}\approx 0$, so that $\mathbf{J}\approx\mathbf{A}^{\text{\sf H}}\mathbf{A}$. Then the walk-summability condition (34) reduces to

where the normalizations $J_{ii}=1$ imply that the columns of $\mathbf{A}$ have unit norm, i.e., that $\|\mathbf{A}\|^{2}_{F}=n$. Note that, if (35) is satisfied, then

Applying these results to the $\kappa(\mathbf{A})$ definition (26), we find

where the latter inequality follows from inspection of (28). We conclude that, in the high-SNR regime, walk summability is sufficient for GGAMP to converge with or without damping.

V Local Stability for Strictly Convex Functions

We next consider the convergence with a more general class of scalar estimation functions ${g}_{s}$ and $g_{x}$: those that are twice continuously differentiable with first derivatives bounded as

for all ${\mathbf{p}}$, $\mathbf{r}$, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{p} and \bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{r}. This condition arises in the important case of minimizing strictly convex functions. Specifically, if GAMP is used in max-sum mode so that the scalar estimation functions are given by (3) and (4) with strictly convex, twice differentiable functions $G_{i}$ and $F_{j}$, then (3), (4), and (16) show that the conditions in (37) will be satisfied.

Let $\mathbf{x}^{t+1}=\bm{f}_{t}(\mathbf{x}^{t})$ for $t=0,1,2,\cdots$ be a dynamical system with a fixed point $\mathbf{x}^{*}$ (i.e., $\bm{f}_{t}(\mathbf{x}^{*})=\mathbf{x}^{*}~{}\forall t$). We say that the system is locally stable at $\mathbf{x}^{*}$ if $\exists\delta>0$ such that, if $\|\mathbf{x}^{0}-\mathbf{x}^{*}\|<\delta$, then $\lim_{t\rightarrow\infty}\mathbf{x}^{t}=\mathbf{x}^{*}$.

Outside of the Gaussian scenario, we have not yet established conditions on the global convergence of GAMP for general scalar estimation functions.Interestingly, it was shown by Moallemi and Van Roy that, for a certain class of convex optimization problems characterized by “scaled diagonal dominance”, max-sum BP converges. As future work, it would be interesting to study whether max-sum GAMP also converges for this class of problems. Instead, we now establish conditions on local stability, as defined in . To simplify the analysis, we will assume that the GAMP algorithm uses arbitrary but fixed stepsize vectors {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{p} and \bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{r}.

Under these assumptions, consider any fixed point $({\mathbf{p}},\mathbf{r})$ of the GAMP method, and define the matrices

evaluated at that fixed point. Note that, under assumption (37), the components of $\mathbf{q}_{s}$ and $\mathbf{q}_{x}$ lie in $(0,1)$. Define the matrix

Then (38)-(39), together with lines 8 and 10 of Algorithm 1, imply

Hence, the column norms of $\widetilde{\mathbf{A}}$ in (39) are less than one. Similar arguments can be use to establish that, for any $i$,

so that $\widetilde{\mathbf{A}}$ also has row norms less than one. We will thus call $\widetilde{\mathbf{A}}$ the row-column normalized matrix.

Consider any fixed point $({\mathbf{s}},\mathbf{x})$ of GAMP Algorithm 1 or Algorithm 2 with fixed vector or scalar stepsizes {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{p} and \bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{r}, respectively, and scalar estimation functions ${g}_{s}$ and $g_{x}$ satisfying the above conditions. Then, the fixed point is locally stable if

for $\widetilde{\mathbf{A}}$ defined in (39). For the Gaussian GAMP algorithm, the same condition implies the algorithm is globally stable.

To relate this condition to Theorem 2, consider the case when {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{s} and \bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x} are fixed points of (19) with $\mathbf{S}=\mathbf{A}.\mathbf{A}$, i.e., the component-wise magnitude square of $\mathbf{A}$. From (41) and (42), we have that

Thus, the peak-to-average ratio of $\widetilde{\mathbf{A}}$ as defined in (26) is bounded below as

Hence, a sufficient condition to satisfy (43) is given by

In comparison, (27) and (29) show that a Gaussian GAMP with scalar step sizes converges is $\kappa(\mathbf{A})<C/(\theta_{s}\theta_{x})$. We conclude that the sufficient condition for the vector-stepsize GAMP algorithm to converge is similar to the scalar-stepsize GAMP algorithm, but where the peak-to-average ratio is measured on a certain normalized matrix.

VI Numerical Results

In this section, we present some numerical simulations to verify Theorems 2 and 3. This section is divided into two parts: the first part is on the global convergence of damped GGAMP (Theorem 2) and the second part is on the local stability of damped GAMP (Theorem 3).

In this experiment, the elements of $\mathbf{x}$ were drawn i.i.d. ${\mathcal{N}}(0,1)$ and the measurements were generated using the AWGN model as discussed above. For each choice of damping factor $\theta_{s}=\theta_{x}$, scalar stepsize GGAMP was run from the fixed initialization $\{\mathbf{x}^{0}\!=\!\bm{0}$, ${\mathbf{s}}^{-1}\!=\!\bm{0}$, $\tau_{x}\!=\!1\}$ and the MSE after $5000$ iterations was recorded. This experiment was then repeated for $100$ realizations of $\{\mathbf{A},\mathbf{x},\mathbf{w}\}$. The damping factors $\theta_{s}=\theta_{x}$ were varied from $0.7$ to $1$ in steps of $0.005$. To test the validity of Theorem 2, we present the results in term of the “excess MSE,” defined as the ratio of the MSE achieved by GAMP to the MMSE, which was computed in closed form. To enhance the readability of the plots, the excess MSE was clipped at $100$ dB.

Figures 1 and 2 show the excess MSE versus $\kappa_{\max}(\theta_{s},\theta_{x})$, which—according to Theorem 2—is the maximum allowed value of $\kappa(\mathbf{A})$ under which GGAMP will converge with damping factors $(\theta_{s},\theta_{x})$, as defined in (27). In both figures, the dimensions of $\mathbf{A}$ were $200\times 100$, and the excess MSE from each realization is plotted as a dot. The figures show that the excess MSE was zero dB whenever $\kappa_{\max}(\theta_{s},\theta_{x})>\kappa(\mathbf{A})$, and conversely the excess MSE was greater than zero dB whenever $\kappa_{\max}(\theta_{s},\theta_{x})<\kappa(\mathbf{A})$, which verifies the claim of Theorem 2.

VI-B Local Convergence of GAMP

To test the local stability of damped GAMP, we used the following procedure. For each realization of $\{\mathbf{A},\mathbf{x},\mathbf{y}\}$, the parameters \{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{p},\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{r},\theta_{s},\theta_{x}\} were chosen and vector-stepsize GAMP was run from the initialization $\{\mathbf{x}^{0}=\bm{0}$, ${\mathbf{s}}^{-1}=\bm{0}$, \bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x}=\bm{1}\} with the stepsizes fixed at the chosen \{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{p},\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{r}\}. The values of \{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{p},\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{r},\theta_{s},\theta_{x}\} were chosen so that GAMP converged to some fixed point $\{\mathbf{x},{\mathbf{s}},{\mathbf{p}},\mathbf{r}\}$; more details are provided below. Next, GAMP was initialized near to the fixed point and tested for local convergence (under the same fixed stepsizes \{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{p},\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{r}\}.) In particular, it was initialized at $\{\mathbf{x}^{0}=\mathbf{x}+\mathbf{x}_{\varepsilon}$, ${\mathbf{s}}^{-1}={\mathbf{s}}\}$, where the elements of $\mathbf{x}_{\varepsilon}$ were drawn i.i.d. ${\mathcal{N}}(0,1)$, with $\mathbf{x}_{\varepsilon}$ subsequently normalized such that the initial MSE was $15$ dB above the MSE at the fixed point. This test was repeated $20$ times for each fixed point. If $\theta_{s}\theta_{x}\|\widetilde{\mathbf{A}}\|_{2}^{2}<1$ then, according to Theorem 3, GAMP should converge to the fixed point. Each dot in Figures 3-6 represents the excess MSE, now defined as the ratio of the maximum MSE among all local runs of GAMP to the MSE at the fixed point. The above procedure was repeated for a range of $\theta_{s}=\theta_{x}$ and many realizations of $\{\mathbf{A},\mathbf{x},\mathbf{y}\}$, as detailed below. As before, the excess MSE values were clipped at $100$ dB before plotting.

Figures 3 and 4 show the excess MSE versus $\theta_{s}\theta_{x}\|\widetilde{\mathbf{A}}\|_{2}^{2}$ for Bernoulli-Gaussian $\mathbf{x}$ with sparsity rate 0.1 and AWGN measurements. Figure 3 investigates the case where $\kappa(\mathbf{A})=4$ and Figure 4 investigates the case where $\kappa(\mathbf{A})=10$. For each plot, the dimensions of $\mathbf{A}$ were $200\times 100$, the stepsizes were $\nu_{p_{i}}=(\sum_{j=1}^{n}A_{ij}^{2})^{-1}~{}\forall i$, the damping factors $\theta_{s}=\theta_{x}$ were varied from $0.45$ to $0.95$ in steps of $0.05$, and $50$ realizations of $\{\mathbf{A},\mathbf{x},\mathbf{y}\}$ were tested. Also, $\tau_{r_{j}}=(\sum_{i=1}^{m}A_{ij}^{2})^{-1}$ in Figure 3 and $\tau_{r_{j}}=(10\sum_{i=1}^{m}A_{ij}^{2})^{-1}$ in Figure 4, for all $j$. This particular choice of \bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{r} was used to ensure that the fixed-stepsized GAMP converged to a fixed point for the chosen range of $\theta_{s},\theta_{x}$.

Figure 5 and 6 show the excess MSE versus $\theta_{s}\theta_{x}\|\widetilde{\mathbf{A}}\|_{2}^{2}$ for Bernoulli-Gaussian $\mathbf{x}$ with sparsity rate 0.1 and binary measurements. Figure 5 investigates the case where $\kappa(\mathbf{A})=4$ and Figure 6 investigates the case where $\kappa(\mathbf{A})=10$. For each plot, the dimensions of $\mathbf{A}$ were $400\times 100$, the stepsizes were $\nu_{p_{i}}=10~{}\forall i$ and $\tau_{r_{j}}=1~{}\forall j$, the damping factors $\theta_{s}=\theta_{x}$ were varied from $0.45$ to $0.95$ in steps of $0.05$, and $50$ realizations of $\{\mathbf{A},\mathbf{x},\mathbf{y}\}$ were tested.

Figures 3-6 show an excess MSE of $\approx 0$ dB whenever $\theta_{s}\theta_{x}\|\widetilde{\mathbf{A}}\|_{2}^{2}<1$, hence verifying Theorem 3.

Conclusions

A key outstanding issue for the adoption of AMP-related methods is their convergence for generic finite-dimensional linear transforms. Similar to other loopy BP-based methods, standard forms of AMP may diverge. In this paper, we presented a damped version of the generalized AMP algorithm that, when used with fixed stepsizes, can guarantee global convergence for Gaussian distributions and local convergence for the minimization of strictly convex functions (i.e., strictly concave log-priors). The required amount of damping is related to the peak-to-average ratio of the squared singular values of the transform matrix. However, much remains unanswered: Most importantly, we have yet to derive a condition for global convergence even in the case of strictly convex functions. Secondly, our analysis assumes the use of fixed stepsizes. Third, short of computing the peak-to-average singular-value ratio, we proposed no method to compute the damping constants. Hence, an adaptive method may be useful in practice. One such method, , has been proposed, but it comes without convergence guarantees. Thus, future work might aim to analyze the convergence of such methods. Also, a more recent algorithm, Vector AMP (VAMP) , has improved convergence on larger classes of random matrices. Another line of future work could seek conditions for convergence of VAMP on deterministic matrices.

Appendix A Proof of Theorem 1

The variance updates of both Algorithms 1 and 2 are both of the form (19) with different choices of $\mathbf{S}$. So, the theorem will be proven by showing that the updates (19) converge for any non-negative matrix $\mathbf{S}\geq 0$. To this end, we use the results in . Specifically, for any {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{w} and \bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{0}>0, define the functions

so that the updates (19) can be written as

It is easy to check that, for any $\mathbf{S}\geq 0$,

\Phi_{s}(\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x})>0,

\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x}\geq\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x}^{\prime}\Rightarrow\Phi_{s}(\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x})\leq\Phi_{s}(\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x}^{\prime}), and

For all $\alpha>1$, \Phi_{s}(\alpha\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x})>(1/\alpha)\Phi_{s}(\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x}).

with the analogous properties being satisfied by \Phi_{x}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{s}). Now let $\Phi:=\Phi_{x}\circ\Phi_{s}$ be the composition of the two functions so that \bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}^{t\!+\!1}_{x}=\Phi(\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}^{t}_{x}). Then, $\Phi$ satisfies the three properties:

\Phi(\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x})>0,

\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x}\geq\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x}^{\prime}\Rightarrow\Phi(\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x})\geq\Phi(\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x}^{\prime}), and

For all $\alpha>1$, \Phi(\alpha\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x})<\alpha\Phi(\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x}).

Also, for any {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{s}\geq 0, we have \Phi_{x}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{s})\leq\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{0} and therefore, \Phi(\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x})\leq\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{0} for all \bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x}\geq 0. Hence, taking any \bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x}\geq\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{0}, we obtain:

Using Theorem 2 in , it can be shown that the updates \bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}^{t\!+\!1}_{x}=\Phi(\bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x}^{t}) converge to a unique fixed point. A similar argument shows that {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{s}^{t} also converges to a unique fixed point.

Appendix B Linear System Stability Condition

The proofs of both Theorems 2 and 3 are based on analyzing the GAMP algorithm via an equivalent linear system and then applying results from linear stability theory. For both results we will show that the condition of the theorem is equivalent to an eigenvalue test on a certain matrix.

First consider the Gaussian GAMP algorithm with fixed vector stepsizes. With fixed stepsizes and Gaussian estimation functions (17), Algorithm 1 reduces to a linear system:

Note that the components of $\mathbf{q}_{s}$ and $\mathbf{q}_{x}$ are in $(0,1)$. We can write the system (45) in matrix form as

for an appropriate matrix $\mathbf{G}$ and vector $\mathbf{b}$. The matrix $\mathbf{G}$ is given by

Note that both $\mathbf{D}_{x}$ and $\mathbf{D}_{s}$ are diagonal matrices with entries in the interval $(0,1)$.

Now, consider the case of the more general scalar estimation functions satisfying (37) and other assumptions in Section V. Due to the differentiability assumptions, to prove the local stability, we only have to look at the linearization of the system around the fixed points . With fixed stepsizes, the linearization of the updates in Algorithm 1 around any fixed point is given by

where the matrices $\mathbf{Q}_{s}$ and $\mathbf{Q}_{x}$ in (46) are replaced by the derivatives (38). This linear system is also of the form (47) with the same matrix (48). Also, under the assumptions of the theorem, $\mathbf{q}_{s}$ and $\mathbf{q}_{x}$ are vectors with components in $(0,1)$.

Hence, we conclude that to prove the global stability of Gaussian GAMP, or the local stability of GAMP under the assumptions of Theorem 3, it suffices to show that the linear system (47) with a matrix $\mathbf{G}$ of the form (48) is stable. The matrices $\mathbf{D}_{s}$ and $\mathbf{D}_{x}$ are given in (49) where $\mathbf{Q}_{s}$ and $\mathbf{Q}_{x}$ are diagonal matrices with elements in $(0,1)$.

To evaluate this condition, first recall that the linear system (47) is stable when the eigenvalues of $\mathbf{G}$ are in the unit circle. However, if we define

the eigenvalues of $\mathbf{G}$ are identical to those of $\mathbf{H}$ given by

Expanding the matrix product in (52), we get

For stability, we need to show that for any $|\lambda|\geq 1$, $\mathbf{H}_{\lambda}$ is invertible. We simplify this condition as follows: Consider any $\lambda$ with $|\lambda|\geq 1$. Now, $\mathbf{D}_{s}$ in (49b) is a diagonal matrix with entries in $[0,1)$. Hence $\lambda\mathbf{I}-\mathbf{D}_{s}$ is invertible since $|\lambda|\geq 1$. Therefore, taking a Schur complement, we see that $\mathbf{H}_{\lambda}$ is invertible if and only if the matrix

is invertible. We can summarize the result as follows.

Consider the GAMP Algorithm 1 for any scalar estimation functions satisfying the conditions in Section V including (37). The GAMP algorithm is locally stable around a fixed point if and only if $\mathbf{J}_{\lambda}$ is invertible for all $|\lambda|\geq 1$, where

and $\mathbf{F}$ is given in (53). In the special case of Gaussian estimation functions (17), the above condition implies the GAMP Algorithm 1, will be globally stable.

A similar calculation can be performed for the GAMP algorithm with scalar stepsizes. In this case, the vector stepsizes such as \bm{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\tau}}_{x} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bm{\nu}}_{s} are replaced with the scalar quantities $\tau_{x}$ and $\nu_{s}$. For the case of Gaussian estimation functions (17) and identical variances (20) we obtain the following:

Consider the GAMP Algorithm 2 with scalar stepsizes, Gaussian scalar estimation functions (17) and identical variances (20). Then, the algorithm is globally stable if and only if $\mathbf{J}_{\lambda}$ is invertible for all $|\lambda|\geq 1$, where

Appendix C Proof of Theorem 2

Our first step in the proof is to simplify the condition in Lemma 2.

Consider the GAMP algorithm with scalar stepsizes, Algorithm 2, with the Gaussian scalar estimation functions (17) and fixed stepsizes. Then the system is stable if and only if

From Lemma 2, we know that the system is stable if and only if $\mathbf{J}_{\lambda}$ in (57) is invertible for all $|\lambda|\geq 1$. To evaluate this condition, suppose that $\mathbf{J}_{\lambda}$ is not invertible for some $|\lambda|\geq 1$. Then, there exists an $\mathbf{v}\neq 0$ such that $\mathbf{J}_{\lambda}\mathbf{v}=0$, which implies that

Using the expression for $\mathbf{F}$ in (58), this is equivalent to

Thus, $\mathbf{v}$ is an eigenvector of $\mathbf{A}^{\text{\sf H}}\mathbf{A}$. But, $\sigma^{2}$ is an eigenvalue of $\mathbf{A}^{\text{\sf H}}\mathbf{A}$ if and only if $\sigma$ is a singular value of $\mathbf{A}$. Hence, we conclude that $\mathbf{J}_{\lambda}$ is non-invertible if and only if there exists a singular value $\sigma$ of $\mathbf{A}$ such that

Equivalently, we have shown that the system is stable if and only if the the second-order polynomial

has stable roots for all singular values of $\mathbf{A}$, $\sigma$. Now recall that $d_{s}$ and $d_{x}\in(0,1)$. By the Jury stability condition, the $p(\lambda)$ has stable roots if and only $p(1)>0$ and $p(-1)>0$. Now, the first condition is always satisfied since

So, the polynomial is stable if and only if

For this to be true for all singular values of $\mathbf{A}$, we need

Thus, the system is stable if and only if (60) is satisfied with

So, we simply need to prove that (62) matches the definition in (61). To this end, first note that

where (a) follows from the definition $q_{x}=\tau_{x}/\tau_{r}$ in (49a) and (b) follows from the fixed-point equation (22b). Similarly using (49b) and (22a), we obtain that

Substituting (63) and (64) into (62), we obtain (61) and the lemma is proven. $\Box$

where $\gamma$ is defined in (61) and the minimization is over $\nu_{w}$ with the other parameters, $\|\mathbf{A}\|^{2}_{F}$, $\tau_{0}$, $m$ and $n$, being fixed. It follows that if

then the system is stable for all $\nu_{w}$. Conversely, if

then there exists at least one $\nu_{w}$ such that the system is unstable. So, the theorem will be proven if we can show that $\Gamma$ defined in (65) matches the expression in (23).

To calculate the minima in (65), it is useful to write a scaled version of the updates. Let

Then, the fixed points of (22) are given by

Moreover, the minimization in (65) is equivalent to

since minimizing over $\nu_{w}$ is equivalent to minimizing over $u$ in the scaled system. To evaluate the minima (69), we first prove the following.

That is, the minima is achieved as $u\rightarrow 0$.

Substituting (67) into (68) and applying (71), we obtain

Now let $s^{\prime}$, $x^{\prime}$ and $A^{\prime}(s,u)$ denote the derivatives with respect to $u$. From (67) we have

Therefore, $(sx)^{2}>\beta$ and hence, from (75), $s^{\prime}>0$. It follows that

since both $2-\theta_{x}-\theta_{s}\geq 0$ and $\theta_{x}\theta_{x}>0$. Hence, from (72), we have

and it follows that the $\gamma$ is minimized by taking $u$ as small as possible. Therefore,

We conclude by evaluating the limit in (70). The following lemma shows that value of the minimization agrees with (23), and hence completes the proof of the theorem.

For any damping constants $\theta_{s}$, $\theta_{x}$, the limit in (70) is given by (23).

First consider the case when $\beta\geq 1$ (i.e. $m\geq n$). In this case, as $u\rightarrow 0$ the solutions to the fixed points (67) will satisfy $s\rightarrow 0$ and $x\rightarrow\infty$. Hence, the limit of $A(s,u)$ in (73) is

where (a) used (72); (b) used (73) and (c) used the fact that $\beta=m/n$. This proves the $m\geq n$ case of (23).

For the case when $\beta<1$ (i.e. $m<n$) and $u=0$, the solutions to fixed point in (67) are

Substituting $s=1-\beta$ and $u=0$ into (72),

where again we have used the fact that $\beta=m/n$. Therefore,

and this proves the $m<n$ case of (23). $\Box$

Appendix D Proof of Theorem 3

Then $0\not\in\mbox{conv}(P)$, the convex hull of $P$.

Write $\lambda$ in polar coordinates, $\lambda=re^{i\theta}$. We first consider the case where $\theta\in(0,\pi)$. Under this assumption, we claim for all $z\in P$,

Since $P$ is compact, this would imply that (77) holds for all $z\in\mbox{conv}(P)$. In particular, $0\not\in\mbox{conv}(P)$. So, we need to show that (77) holds for all $z\in P$.

Now, since $\theta\in(0,\pi)$, $\sin\theta>0$. Also, since $|\lambda|\geq 1$, $r\geq 1$. Therefore, $r^{2}>d_{s}d_{x}$ since $d_{s},d_{x}<1$. Hence, (79) shows that (77) holds for all $z\in P$.

Similarly, for the case when $\theta\in(-\pi,0)$, (79) shows that

for all $z\in P$. The same argument then shows that $0\not\in\mbox{conv}(P)$.

It remains to consider the cases when $\theta=0$ or $\theta=\pi$. For $\theta=0$, $\lambda=r$ and any $z\in P$ is of the form,

where (a) follows from the fact that $r>d_{s}$ and (b) follows from the fact that $r>d_{x}$. So, for all $z\in P$, $z$ is real and positive. Hence, $0\not\in\mbox{conv}(P)$. Similarly, when $\theta=\pi$, $\lambda=-r$ and

where (a) follows since $d_{x}>0$ and (b) follows since $r\geq 1$ and $\sigma^{2}<1$. Therefore, for all $z\in P$, $z$ is real and negative. Hence, $0\not\in\mbox{conv}(P)$. We have thus shown that $0\not\in\mbox{conv}(P)$ for all values of $\theta$. $\Box$

We can now prove the main result. Suppose that (43) is satisfied. By the definition of $\mathbf{F}$ in (53) and $\widetilde{\mathbf{A}}$ in (39), we have that

Suppose that $\mathbf{J}_{\lambda}$ in (56) is not invertible for some $\lambda$ with $|\lambda|\geq 1$. Then, there exists an $\mathbf{x}$ with $\|\mathbf{x}\|^{2}=1$ such that $\mathbf{x}^{\text{\sf H}}\mathbf{J}_{\lambda}\mathbf{x}=0$. Therefore, if we define $\mathbf{y}=\mathbf{F}\mathbf{x}$, the definition of $\mathbf{J}_{\lambda}$ in (56) shows that

Since $\mathbf{D}_{x}$ and $\mathbf{D}_{s}$ are diagonal, we have

Since $\|\mathbf{x}\|^{2}=1$, we have $\sum_{j}|x_{j}|^{2}=1$. Also, since $\|\mathbf{F}\|_{2}^{2}=\sigma^{2}_{\max}(\mathbf{F})<1$,

for some $\sigma^{2}<1$. Therefore, (82) shows that

Now, from (38) and the contractivity assumption (37), the elements of the diagonal matrices $\mathbf{Q}_{x}$ and $\mathbf{Q}_{s}$ must be in the interval $(0,1)$. Hence, from (49), the elements $d_{x_{j}}$ and $d_{s_{j}}\in(0,1)$. Therefore, $d_{x},d_{s,\max}$ in (84) are in $(0,1)$. From Lemma 6, $0\not\in\mbox{conv}(P_{\lambda})$ which is a contradiction of (83). Hence, the assumption that $\mathbf{J}_{\lambda}$ is not invertible must be false, and the theorem is proven.

References