S-AMP: Approximate Message Passing for General Matrix Ensembles

Burak Çakmak, Ole Winther, Bernard H. Fleury

I Introduction

Consider an $N\times K$ linear observation model described by

The scalar functions $\eta_{t}$ , $t\geq 0$ , in (2) are obtained by applying an additional optimization procedure based upon the so-called state evolution formula for the underlying measurement matrix ensemble . In (3), $\eta^{\prime}_{t}(x)=d\eta_{t}(x)/dx$ , $t\geq 0$ . Moreover for a vector ${\mathchoice{\mbox{\boldmath$ \displaystyle u $}}{\mbox{\boldmath$ \textstyle u $}}{\mbox{\boldmath$ \scriptstyle u $}}{\mbox{\boldmath$ \scriptscriptstyle u $}}}\triangleq(u_{1},\dots,u_{K})$ , $\left<{\mathchoice{\mbox{\boldmath$ \displaystyle u $}}{\mbox{\boldmath$ \textstyle u $}}{\mbox{\boldmath$ \scriptstyle u $}}{\mbox{\boldmath$ \scriptscriptstyle u $}}}\right>\triangleq\sum_{k=1}^{K}u_{k}/K$ and $\alpha\triangleq N/K$ . The vectors ${\mathchoice{\mbox{\boldmath$ \displaystyle\mu $}}{\mbox{\boldmath$ \textstyle\mu $}}{\mbox{\boldmath$ \scriptstyle\mu $}}{\mbox{\boldmath$ \scriptscriptstyle\mu $}}}^{t}$ and ${\mathchoice{\mbox{\boldmath$ \displaystyle z $}}{\mbox{\boldmath$ \textstyle z $}}{\mbox{\boldmath$ \scriptstyle z $}}{\mbox{\boldmath$ \scriptscriptstyle z $}}}^{t}$ are referred to as the current estimate of $\textstyle x$ and the corresponding residual, respectively. Finally $(\cdot)^{\dagger}$ denotes transposition.

AMP has two appealing properties. Firstly, when the entries of $\textstyle A$ are independent identically distributed (iid) Gaussian with zero mean and variance $1/N$ , AMP yields the minimum mean square error (MMSE) estimator in the large system limit . Secondly, AMP includes a so-called Onsager reaction term, i.e, $\alpha^{-1}\left<\eta_{t-1}^{\prime}(\cdot)\right>{\mathchoice{\mbox{\boldmath$ \displaystyle z $}}{\mbox{\boldmath$ \textstyle z $}}{\mbox{\boldmath$ \scriptstyle z $}}{\mbox{\boldmath$ \scriptscriptstyle z $}}}^{t-1}$ in (3), that corrects the naive mean field approximation. In statistical physics such a technique is known as the Thouless-Anderson-Palmer (TAP) correction .

The adaptive TAP (ADATAP) mean field theory was introduced in . In ADATAP the form of Onsager reaction term depends on the measurement matrix, see [4, Eq. (20) & (51)]. Indeed, a connection between ADATAP and AMP has been recently realized in . The connection is based on some approximations of the Gibbs free energy, which are derived using the replica method, see [5, Eq. (10) & (11)] and the references therein.

Inference techniques based on the free energy optimization have become popular in the literature of information theory and in machine learning , and references therein. The important results exploited in this contribution is that the fixed points of belief propagation (BP) and expectation propagation (EP) are the stationary points of the Bethe Free energy (BFE) under a set of marginalization consistency constraints and moment consistency constraints , respectively.

The conventional approximate message passing methods presented in the literature are based on a Gaussian approximation of loopy BP on a dense graph, . By contrast, the method presented in this paper is based on probabilistic inference on a tree graph. Specifically we consider an exact Gibbs free energy formulation (i.e. a BFE formulation on a tree probabilistic graph) under first and second-moment consistency constraints. Our analysis relies on the stationary point equations of the constrained Gibbs free energy. In particular we propose a novel algorithm whose fixed points are the stationary points of the constrained Gibbs free energy in the large system limit. This algorithm – we coin it S-AMP – executes the following iteration steps:

with $\rm S_{{\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}}$ denoting the S-transform of the asymptotic eigenvalue distribution (AED) of ${\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}^{\dagger}{\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}$ (see, e.g. ). Later in the paper we will show that the optimality of S-AMP follows by its design rather than based upon an optimization procedure as in .

To show that AMP is a special case of S-AMP, let the entries of $\textstyle A$ be iid with zero mean variance $1/N$ . Then, as $N,K\to\infty$ with the ratio $\alpha=N/K$ fixed, $\rm S_{{\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}}(\omega)=1/(1+\omega/\alpha)$ [13, Eq. (2.87)]. Inserting this expression in (6) we obtain the iteration steps (2)-(3) of AMP.

II Gibbs Free Energy with Moment Constraints

Consider the $N\times K$ linear observation model (1). For Bayesian inference, we assign a prior $p_{k}(x_{k})$ for all $k\in\mathcal{K}$ . Hence the joint posterior pdf can be written as

with $p({\mathchoice{\mbox{\boldmath$ \displaystyle y $}}{\mbox{\boldmath$ \textstyle y $}}{\mbox{\boldmath$ \scriptstyle y $}}{\mbox{\boldmath$ \scriptscriptstyle y $}}}|{\mathchoice{\mbox{\boldmath$ \displaystyle x $}}{\mbox{\boldmath$ \textstyle x $}}{\mbox{\boldmath$ \scriptstyle x $}}{\mbox{\boldmath$ \scriptscriptstyle x $}}})$ and $Z$ denoting the likelihood given by (1) and a normalization constant, respectively. The factor graph representation of (7) is a tree. Thus the BFE for (7) is equal to the Gibbs free energy [6, Theorem 3], which is given by

When we define a Lagrangian for (8) that accounts for the set of marginalization consistency constrains, then at its stationary point, the belief $b_{k}(x_{k})$ is equal to $p(x_{k}|{\mathchoice{\mbox{\boldmath$ \displaystyle y $}}{\mbox{\boldmath$ \textstyle y $}}{\mbox{\boldmath$ \scriptstyle y $}}{\mbox{\boldmath$ \scriptscriptstyle y $}}})$ for all $k\in\mathcal{K}$ . We consider the Gibbs energy formulation with a set of moment consistency constraints, instead of marginalization constraints. Specifically, following the arguments of we define the Lagrangian

The term $\mathcal{Z}$ accounts for the set of the normalization constraints for the beliefs:

We formulate the estimation procedure for $x_{k},k\in\mathcal{K}$ as

where $b_{k}^{\star}(x_{k})$ represents the belief of $x_{k}$ at a stationary point of (9).

In the sequel we derive the stationary points equations of the Lagrangian (9). For the sake of notational compactness we define

In (13) we have introduced the $K\times K$ diagonal matrix and the $K\times 1$ vector $\textstyle\bar{\gamma}$ whose entries are respectively $\bar{\Lambda}_{kk}=\bar{\lambda}_{k}$ and $\bar{\gamma}_{k}$ , $k\in\mathcal{K}$ .

The stationary points of the Lagrangian (9) are obtained to be of the form

Furthermore we define for any $k\in\mathcal{K}$

It is shown in [12, Eq. (31)-(35)] that $\eta(\kappa_{k};\lambda_{k})$ and $\eta^{\prime}(\kappa_{k};\lambda_{k})/\lambda_{k}$ give the mean and the variance of the belief (17), respectively. With these definitions, the identities resulting from the moment consistency constraints are given by

We now derive a simple expression for (11). By making use of the identities in (16) and (20), we write first

Furthermore by the definitions in (13) we have

Let us introduce the $K\times K$ diagonal matrix $\textstyle\Lambda$ and the $K\times 1$ vector $\textstyle\gamma$ whose entries are respectively $\Lambda_{kk}=\lambda_{k}$ and $\gamma_{k}$ , $k\in\mathcal{K}$ . Then, making use of the identity in (23) we can write

where $\text{diag}({\mathchoice{\mbox{\boldmath$ \displaystyle\Sigma $}}{\mbox{\boldmath$ \textstyle\Sigma $}}{\mbox{\boldmath$ \scriptstyle\Sigma $}}{\mbox{\boldmath$ \scriptscriptstyle\Sigma $}}})$ is the $K\times K$ diagonal matrix with $\text{diag}({\mathchoice{\mbox{\boldmath$ \displaystyle\Sigma $}}{\mbox{\boldmath$ \textstyle\Sigma $}}{\mbox{\boldmath$ \scriptstyle\Sigma $}}{\mbox{\boldmath$ \scriptscriptstyle\Sigma $}}})_{kk}=\Sigma_{kk}$ , $k\in\mathcal{K}$ . Then, by invoking the identities (20) and(22) we arrive at the sought explicit form for (11):

As a matter of fact equations (27)–(29) coincide with the fixed point equations of ADATAP that are obtained by applying the cavity approach-new in statistical physics, see [4, Eq. (20), (25) and (26)].

The step in (29) requires a matrix inversion, which is desirable to avoid in order to keep the complexity of fixed point algorithms devised from (27)–(29) low. In the authors circumvent this complexity problem by using the so-called self-averaging method [4, Section 3.1] in the large system limit. The following theorem restates a result presented in [4, Section 3.1] in terms of the function $\eta^{\prime}$ and the R-transform in free probability (see e.g. ).

[4, Section 3.1] Let ${\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}^{\dagger}{\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}$ have an AED as $N,K\to\infty$ with the ratio $\alpha=N/K$ fixed. Let $\left<\eta^{\prime}({\mathchoice{\mbox{\boldmath$ \displaystyle\kappa $}}{\mbox{\boldmath$ \textstyle\kappa $}}{\mbox{\boldmath$ \scriptstyle\kappa $}}{\mbox{\boldmath$ \scriptscriptstyle\kappa $}}};{\mathchoice{\mbox{\boldmath$ \displaystyle\Lambda $}}{\mbox{\boldmath$ \textstyle\Lambda $}}{\mbox{\boldmath$ \scriptstyle\Lambda $}}{\mbox{\boldmath$ \scriptscriptstyle\Lambda $}}})\right>\triangleq\frac{1}{K}\sum_{k\in\mathcal{K}}\eta^{\prime}(\kappa_{k};\lambda_{k})$ . Then, as $N,K\to\infty$ with the ratio $\alpha=N/K$ fixed, for all $k\in\mathcal{K}$ $\lambda_{k}$ converges almost surely to the macroscopic quantity $\lambda$ that is the solution ofBy abusing the notation we define $\eta^{\prime}({\mathchoice{\mbox{\boldmath$ \displaystyle\kappa $}}{\mbox{\boldmath$ \textstyle\kappa $}}{\mbox{\boldmath$ \scriptstyle\kappa $}}{\mbox{\boldmath$ \scriptscriptstyle\kappa $}}};\lambda)\triangleq\eta^{\prime}({\mathchoice{\mbox{\boldmath$ \displaystyle\kappa $}}{\mbox{\boldmath$ \textstyle\kappa $}}{\mbox{\boldmath$ \scriptstyle\kappa $}}{\mbox{\boldmath$ \scriptscriptstyle\kappa $}}};\lambda{\mathchoice{\mbox{\boldmath$ \displaystyle I $}}{\mbox{\boldmath$ \textstyle I $}}{\mbox{\boldmath$ \scriptstyle I $}}{\mbox{\boldmath$ \scriptscriptstyle I $}}})$ , with $\textstyle I$ denoting the identity matrix of appropriate dimension.

with $\rm R_{\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}$ denoting the R-transform of the AED of ${\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}^{\dagger}{\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}$ .

Making use of the relation between the R-transform and the S-transform [14, Table 6] in (30) we obtain the following corollary.

Let the random matrix $\textstyle A$ be defined as in Theorem 1. Then, we have

with $\rm S_{{\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}}$ denoting the S-transform of the AED of ${\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}^{\dagger}{\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}$ .

III Fixed Point Algorithms

In this section we use the stationary point equations obtained in the previous section to introduce three fixed point iterative algorithms. Firstly we will present the classical EP scheme for (1) and the ADATAP scheme . Secondly we derive the S-AMP algorithm mentioned in the introduction.

All three recovery schemes have the following basis update step in common, which results by time-indexing the first identity in (13):

Since only one element of ${\mathchoice{\mbox{\boldmath$ \displaystyle\bar{\Lambda} $}}{\mbox{\boldmath$ \textstyle\bar{\Lambda} $}}{\mbox{\boldmath$ \scriptstyle\bar{\Lambda} $}}{\mbox{\boldmath$ \scriptscriptstyle\bar{\Lambda} $}}}^{t}$ is updated in each iteration the matrix inversion lemma can be applied to reduce the complexity of this step to $O(K^{2})$ , e.g. see [9, Eq. (37)]. This makes (32) suitable for applications with moderately large dimensions of $\textstyle A$ .

In the following we present the compact form of the EP scheme for (1) (e.g. see ) and the ADATAP scheme . First we start with defining update steps common to both algorithms. They follow by merely time indexing (29) for $k\in\mathcal{K}$ :

EP updates $\mu_{k}^{t},k\in\mathcal{K}$ based on the second identity in (13), (21) and (28):

ADATAP updates $\mu^{t}_{k},k\in\mathcal{K}$ based on the stationary points identities in (27)–(28):

Depending on the system model, EP and ADATAP may exhibit a poor convergence behavior, and may even diverge. A procedure to improve the convergence behavior consists in introducing a damping factor, say $\epsilon$ , when updating e.g. $\mu^{t}_{k}$ in (36) and (38) as $(1-\epsilon)\mu_{k}^{t}+\epsilon\eta(\kappa_{k}^{t};\lambda_{k}^{t})\to\mu_{k}^{t}$ . However this approach leads to very slow convergence which might require thousands of iterations, e.g. see [5, Section V]. Regarding more advanced damping schemes we refer the reader to .

III-B S-AMP

In the sequel we derive a new fixed point algorithm from the stationary points identities (27)–(29). The algorithm yields S-AMP in the large system limit.

From this definition we “devise” the following identity:

Making use of (27), (28) (with definition (39)), and (40) we obtain the new fixed point algorithm

where $\lambda_{k}^{t}$ satisfies the system of equations

Like AMP, this scheme includes by design a natural damping factor $(1-\sigma_{w}^{2}\lambda_{k}^{t-1})$ for the contribution $z_{n,k}^{t-1}$ . Specifically in this scheme just like in AMP, we do not need a step-size parameter. However,at each iteration solving $\lambda_{k}^{t}$ from (43) is non-trivial in general. In this respect, the scheme in (33) can be considered as an approximation of (43).

By the design of $\lambda_{k}^{t}$ through (43), and Theorem 1, for all $k\in\mathcal{K}$ , $\lambda^{t}_{k}$ converges almost surely to a macroscopic quantity $\lambda^{t}$ as $N,K\to\infty$ with the ratio $\alpha$ fixed. Furthermore invoking Corollary 1 the quantity $\lambda^{t}$ is the solution of the identity

Consequently we obtain the iteration steps (4)-(6) of S-AMP in their scalar form:

We note that by the definition, ${\mathchoice{\mbox{\boldmath$ \displaystyle\kappa $}}{\mbox{\boldmath$ \textstyle\kappa $}}{\mbox{\boldmath$ \scriptstyle\kappa $}}{\mbox{\boldmath$ \scriptscriptstyle\kappa $}}}^{t}={\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}^{\dagger}{\mathchoice{\mbox{\boldmath$ \displaystyle z $}}{\mbox{\boldmath$ \textstyle z $}}{\mbox{\boldmath$ \scriptstyle z $}}{\mbox{\boldmath$ \scriptscriptstyle z $}}}^{t}+{\mathchoice{\mbox{\boldmath$ \displaystyle\mu $}}{\mbox{\boldmath$ \textstyle\mu $}}{\mbox{\boldmath$ \scriptstyle\mu $}}{\mbox{\boldmath$ \scriptscriptstyle\mu $}}}^{t}$ .

In , the function $\eta_{t}({\mathchoice{\mbox{\boldmath$ \displaystyle\kappa $}}{\mbox{\boldmath$ \textstyle\kappa $}}{\mbox{\boldmath$ \scriptstyle\kappa $}}{\mbox{\boldmath$ \scriptscriptstyle\kappa $}}}^{t})$ in AMP updates is referred as “an appropriate sequence of non-linear functions”. By contrast, by design of the iterative process of S-AMP, we have the definition of $\eta_{t}({\mathchoice{\mbox{\boldmath$ \displaystyle\kappa $}}{\mbox{\boldmath$ \textstyle\kappa $}}{\mbox{\boldmath$ \scriptstyle\kappa $}}{\mbox{\boldmath$ \scriptscriptstyle\kappa $}}}^{t})$ via the fixed point equation (44). Note that, $\lambda^{t}$ must be solved at each iteration $t$ from this equation. Depending on the prior pdf’s, obtaining closed form expression for $\lambda^{t}$ is often non-trivial. In fact this shows how S-AMP (or AMP in particular) can be a very advanced estimator as (44) directly relates the asymptotic stationary point identity in (31). In order to better comprehend this aspect, in the following we examine $\lambda^{t}$ for the linear estimation problem.

The optimality of AMP for the linear estimation problem with the zero mean iid matrix ensemble, was proven in [2, Section 2.1] via a minimization procedure upon the state evolution formula. We, by contrast, have the definition of S-AMP of which we can show the optimality for the general matrix ensembles.

Consider the linear observation model $\eqref{system}$ . Let the entries of $\textstyle x$ be iid Gaussian with zero mean and variance one, i.e. $p_{k}(x_{k})=N(x_{k}|0,1)$ , $k\in\mathcal{K}$ . Then the asymptotic MMSE of (1) reads In , the notation $\eta$ is used for (48). For convenience we adopt the notation $\tau$ .

with ${\rm P}_{{\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}}(x)$ denoting the AED of ${\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}^{\dagger}{\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}$ . Recall that the fixed points of S-AMP are the stationary points of the Gibbs free energy under the moment consistency constraints in the large system limit. Therefore, for the given a Gaussian prior, S-AMP must be a MMSE estimator in the large system limit. Namely the following relation must be satisfied:

We show next that actually for any $t\geq 0$ , $\left<\eta_{t}^{\prime}({\mathchoice{\mbox{\boldmath$ \displaystyle\kappa $}}{\mbox{\boldmath$ \textstyle\kappa $}}{\mbox{\boldmath$ \scriptstyle\kappa $}}{\mbox{\boldmath$ \scriptscriptstyle\kappa $}}}^{t})\right>/\lambda^{t}=\tau_{{\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}}(\sigma_{w}^{2})$ . First notice that with the choice of the prior we have $\left<\eta_{t}^{\prime}({\mathchoice{\mbox{\boldmath$ \displaystyle\kappa $}}{\mbox{\boldmath$ \textstyle\kappa $}}{\mbox{\boldmath$ \scriptstyle\kappa $}}{\mbox{\boldmath$ \scriptscriptstyle\kappa $}}}^{t})\right>=\lambda^{t}/(1+\lambda^{t})$ . From the definition in (44), we have

The S-transform can be formulated in terms of $\tau_{{\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}}(\sigma_{w}^{2})$ [13, Definition 2.15]. Using this formula we write

Thus $\lambda^{t}=1/\tau_{{\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}}(\sigma_{w}^{2})-1$ , which confirms the optimality of S-AMP for the linear estimation problem.

IV A Sub-Optimal Variant of S-AMP

In the previous subsection we derived the explicit expression for $\lambda^{t}$ when the prior pdf’s are Gaussian. However solving $\lambda^{t}$ from (44) for the other prior pdf’s is often non-trivial. A direct approach consists in including an inner loop to solve (44) iteratively at each iteration. That would, however, create an overhead that we would like to avoid. Instead, we propose in the following a sub-optimal scheme for $\lambda_{t}$ that does not require any inner loop. We approximate the optimal $\lambda^{t}$ defined by (44) with $\lambda_{s}^{t}$ that satisfies

Here we note that, the sub-optimal scheme coincides with the same fixed point equations of the optimal scheme.

When the entries of $\textstyle A$ are iid with zero mean and variance $1/N$ , the sub-optimal scheme coincides with the classical recursion of AMP used in the literature e.g. . In fact, from (52), it is easy to obtain the so-called state-evolution formula for the iid zero mean matrix ensemble. Finally, we note that it is possible to introduce more advanced recovery schemes. But this is out of the scope of this contribution.

In the sequel we assess the performance of the sub-optimal variant of S-AMP. Due to the space limitation we only consider the system model used in [5, Section 5] for Bayesian inference in compressed sensing. Accordingly, the prior pdf’s are Bernoulli-Gaussian: $p_{k}(x_{k})=(1-\rho)\delta(x_{k})+\rho N(x_{k}|0,1)$ , $k\in\mathcal{K}$ . We refer the reader to [12, Eq. (67) & (68)] for the closed form expressions of $\eta_{t}(\cdot)$ and $\eta_{t}^{\prime}(\cdot)$ in this case.

We consider the sub-optimal variant of S-AMP for two scenarios: i) the random row-orthogonal matrix ensemble, i.e. ${\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}=\alpha^{-\frac{1}{2}}{\mathchoice{\mbox{\boldmath$ \displaystyle P $}}{\mbox{\boldmath$ \textstyle P $}}{\mbox{\boldmath$ \scriptstyle P $}}{\mbox{\boldmath$ \scriptscriptstyle P $}}}_{\alpha}{\mathchoice{\mbox{\boldmath$ \displaystyle O $}}{\mbox{\boldmath$ \textstyle O $}}{\mbox{\boldmath$ \scriptstyle O $}}{\mbox{\boldmath$ \scriptscriptstyle O $}}},\alpha\leq 1$ , where ${\mathchoice{\mbox{\boldmath$ \displaystyle P $}}{\mbox{\boldmath$ \textstyle P $}}{\mbox{\boldmath$ \scriptstyle P $}}{\mbox{\boldmath$ \scriptscriptstyle P $}}}_{\alpha}$ is the $N\times K$ matrix with entries $[{\mathchoice{\mbox{\boldmath$ \displaystyle P $}}{\mbox{\boldmath$ \textstyle P $}}{\mbox{\boldmath$ \scriptstyle P $}}{\mbox{\boldmath$ \scriptscriptstyle P $}}}_{\alpha}]_{ij}=\delta_{ij},\forall ij$ , with $\delta_{ij}$ denoting the Kronecker delta, and $\textstyle O$ is the $K\times K$ Haar matrix ; ii) iid zero mean Gaussian matrix ensemble. Note that in the latter case, the sub-optimal variant coincides with the classical AMP recursion as in . In the former case, with a straightforward calculus in free probability we obtain that ${\rm S}_{{\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}}(\omega)=(1+\omega)/(1+\omega/\alpha)$ and

where $\chi^{t}\triangleq\left<\eta_{t-1}^{\prime}({\mathchoice{\mbox{\boldmath$ \displaystyle A $}}{\mbox{\boldmath$ \textstyle A $}}{\mbox{\boldmath$ \scriptstyle A $}}{\mbox{\boldmath$ \scriptscriptstyle A $}}}^{\dagger}{\mathchoice{\mbox{\boldmath$ \displaystyle z $}}{\mbox{\boldmath$ \textstyle z $}}{\mbox{\boldmath$ \scriptstyle z $}}{\mbox{\boldmath$ \scriptscriptstyle z $}}}^{t-1}+{\mathchoice{\mbox{\boldmath$ \displaystyle\mu $}}{\mbox{\boldmath$ \textstyle\mu $}}{\mbox{\boldmath$ \scriptstyle\mu $}}{\mbox{\boldmath$ \scriptscriptstyle\mu $}}}^{t-1})\right>/(\alpha\sigma_{w}^{2}\lambda_{s}^{t-1})$ .

In [5, Section 5], the authors report the estimated normalized mean square estimation error (nmsee) of the damped-ADATAP scheme for the setting (i). For each trial up to 3000 iterations are executed. In Figure 1 we report the nmsee of the suboptimal variant of S-AMP applied in the same context versus the number of iterations. Details are reported in the caption of Fig. 1. Note that no divergence behavior was observed in all performed trials. A comparison of the curves in Fig. 1 with the corresponding curves reported in [5, Fig. 1] show that both recovery schemes achieve the same performance, but with a significantly smaller number of iterations for the sub-optimal variant.

V Conclusion

We developed a novel low-complexity fixed-point algorithm for linear observation systems from the equations of the stationary points of the exact Gibbs free energy under first- and second-moment consistency constraints in the large system limit. The algorithm that we call S-AMP extends AMP for general matrix ensembles. Specifically, AMP is a special case of S-AMP when the measurement matrix has iid zero mean entries. The optimality of S-AMP follows by its design. Furthermore, we define a sub-optimal variant of S-AMP, which is easy to implement. This sub-optimal recovery scheme shows excellent performance for the row-orthogonal matrix ensemble in compressed sensing and it converges in around 40 iterations without showing any divergence behavior.