# Reconciling AMP Algorithms derived from Belief Propagation or the Large System Limit Bethe Free Energy

Zilu Zhao, Fangqing Xiao, Christo Thomas, Dirk Slock

► To cite this version:

Zilu Zhao, Fangqing Xiao, Christo Thomas, Dirk Slock. Reconciling AMP Algorithms derived from Belief Propagation or the Large System Limit Bethe Free Energy. ICASSP 2025, IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), IEEE, Apr 2025, Hyderabad, India. pp.1-5, 10.1109/ICASSP49660.2025.10889454. hal-04991590

HAL Id: hal-04991590

https://hal.science/hal-04991590v1

Submitted on 14 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Reconciling AMP Algorithms derived from Belief Propagation or the Large System Limit Bethe Free Energy

Zilu Zhao<fnref n="1" />, Fangqing Xiao<fnref n="1" />, Christo K. Thomas<fnref n="2" /> and Dirk Slock<fnref n="1" />

**Abstract**—When derived from the Bethe Free Energy (BFE) of the Generalized Linear Model (GLM), Approximate Message Passing (AMP) algorithms combine two asymptotic Large System Limit (LSL) simplifications which are asymptotic Gaussianity of extrinsics and large random matrix theory based asymptotic variance computations. In the provably convergent AMBGAMP algorithm, a LSL version of the BFE is derived. In Expectation Propagation (EP) style minimization, the LSL BFE cost function is augmented with Lagrangian terms for mean and variance consistency constraints, augmented with a quadratic version of the mean constraints as in the Method of Multipliers (MM). The mean Lagrange multipliers then get updated ADMM-style (Alternating Direction of MM). In this approach, the weights of the MM terms need to be carefully chosen, which is not part of the MM philosophy, and the Lagrange multipliers have no particular meaning. On the other hand, AMP can be derived by directly introducing LSL simplifications in the Belief Propagation (BP) algorithm that minimizes the original GLM BFE. This allows to relate extrinsic messages to posterior pdfs by first-order Taylor series expansion based perturbations. We also apply LSL approximations to the variances of the various Gaussians involved, which in fact leads to a rederivation of a fundamental LSL theorem describing the deterministic limit of posterior variances. We show that this LSL version of BP leads to BFE modifications that correspond to the augmented Lagrangian of the LSL BFE, explaining its weights and Lagrange Multipliers. These insights should facilitate the extension of AMP to more complex settings such as bilinear models.

## I. INTRODUCTION

Sparse signal recovery is a fundamental problem in signal processing with a wide range of applications. Many of these problems can be framed as the task of estimating a latent vector $x$ based on a correlated observation vector $y$ [1]. In the Bayesian framework, the complexity of Canonical Methods such as MMSE and MAP experiences exponential growth as the dimension of the problem grows.

By exploiting the structure of the models, graphical model based methods prove to be effective. Belief Propagation (BP) transforms the global inference problem into a local inference problem as outlined by [2]. Expectation Propagation (EP) was introduced in [3] and has been shown to share a similar updating scheme as BP, but for computational efficiency, the messages in BP are projected into a suitable member of the family of exponential distributions [3].

In both [1] and [4], the authors unify EP and BP within the framework of minimizing variational free energy. They demonstrate the close relationship between the fixed points of various message-passing algorithms and the stationary points of Bethe Free Energy (BFE).

EP can serve as an inference method in the generalized linear model (GLM). However, the computational cost corresponds to propagating $2MN$ messages as in Fig. 1 when the data matrix $A$ is of size $M \times N$. Generalized Approximate Mes-

sage Passing (GAMP) [5] builds upon EP, but through the application of large system approximations (LSA), it effectively reduces the number of messages to $M+N$ extrinsics and (marginal) posteriors, providing a more computationally efficient approach.

In [6], the authors investigated the fixed points of the Generalized AMP (GAMP) algorithm for GLMs. They discovered that GAMP shares the same fixed point as the stationary points of the Large System Limit Bethe Free Energy (LSL BFE). In [7] we then proposed AMBGAMP which is guaranteed to converge. Building upon the works of [1], [8], [7], [9], [10], and [11], we present the contributions described in the abstract.

## II. BETHE FREE ENERGY OF THE GENERALIZED LINEAR MODEL

### A. Bethe Free Energy (BFE)

Consider a pdf factorization

$$p(x, y) \propto \prod_{\alpha} f_{x_{\alpha}}(x_{\alpha}), \quad (1)$$

where $x_\alpha$ is a subvector of $x$. In case of a tree-structured factor graph, an alternative equivalent form is [2]

$$p(x|y) = \frac{\prod_{\alpha} p(x_{\alpha})}{\prod_{i} p(x_i)^{M_i-1}}, \quad (2)$$

where $M_i$ is the number of subvectors $x_\alpha$ that contain $x_i$. In (2), the $p(x_\alpha)$ and $p(x_i)$ are the exact factor (subvector) resp. variable marginals.

The concept of variational free energy suggests that to infer the marginals from a tree structured $p(x, y)$ given in (1), we can use as trial distribution

$$q_x(x) = \frac{\prod_{\alpha} q_{x_\alpha}(x_\alpha)}{\prod_{i} q_{x_i}(x_i)^{M_i-1}}. \quad (3)$$

The true marginals can be obtained by [1]

$$ \begin{gathered} \min_{q_{x\alpha}(x_\alpha), q_{x_i}(x_i)} F = D[q(x) || \prod_\alpha f_{x_\alpha}(x_\alpha)]; \\ \text{s.t. } \forall \alpha, \forall i \in \mathcal{I}_\alpha, q_{x_i}(x_i) = \int q_{x_\alpha}(x_\alpha) dx_{\bar{x}_i}, \end{gathered} \quad (4) $$

where we define the shorthand notation (for arbitrary nonnegative functions $q, p$) $D(q||p) = \int q(x) \ln \frac{q(x)}{p(x)} dx$ (which is the Kullback-Leibler Divergence (KLD) in case of normalized $q, p$) and $\bar{x}_i$ denotes all $x$ except $x_i$. The free energy can be expanded as

$$F = \sum_{\alpha} D[q_{x\alpha}(x_{\alpha}) || f_{x\alpha}(x_{\alpha})] + \sum_{i} (M_i - 1) H[q_{x_i}(x_i)], \quad (5)$$

<footnote n="1">¹Communication Systems Department, EURECOM, France</footnote>

<footnote n="2">²NEWS Lab, Virginia Tech, VA, USA</footnote>

<footnote>{Zilu.Zhao, Fangqing.Xiao, Dirk.Slock}@eurecom.fr, ChristoKT@vt.edu</footnote>
where $H(.)$ denotes entropy in nats. Note that this representation only holds for a tree structured distribution. For general graphs that contain loops, (2) no longer holds. Thus, in cases with loops, (5) is only an approximation of the variational free energy. The expression (5) is instead called Bethe free energy.

## B. BFE of the GLM for BP

We consider a GLM with

$$p(\mathbf{x}) = \prod_{i=1}^{N} p(x_i), \quad \mathbf{z} = \mathbf{A}\mathbf{x}, \quad p(\mathbf{y}|\mathbf{z}) = \prod_{j=1}^{M} p(y_j|z_j), \quad (6)$$

where the ratio $N/M$ is a constant for large system considerations. We interpret the linear mixing as a conditional probability
$$p(\mathbf{z}|\mathbf{x}) = \delta(\mathbf{z} - \mathbf{A}\mathbf{x}). \quad (7)$$

From this general linear model, a joint (loopy) factorization scheme comes up naturally:

$$p(\mathbf{x}, \mathbf{z}|\mathbf{y}) \propto p(\mathbf{x}, \mathbf{y}, \mathbf{z}) = p(\mathbf{y}|\mathbf{z})\delta(\mathbf{z} - \mathbf{A}\mathbf{x})p(\mathbf{x}). \quad (8)$$

According to the definition of BFE (5), the associated BFE based on the joint factorization scheme (8) is calculated [1] as

$$F = D[q_x(\mathbf{x})\|p(\mathbf{x})] + D[q_z(\mathbf{z})\|p(\mathbf{y}|\mathbf{z})] + \sum_{i} H[q_{x_i}(x_i)] \\ + D[b_{x,z}(\mathbf{x}, \mathbf{z})\|\delta(\mathbf{z} - \mathbf{A}\mathbf{x})] + \sum_{j} H[q_{z_j}(z_j)], \quad (9)$$

where $q_x$, $q_z$, $b_{x,z}$, $q_{x_i}$ and $q_{z_j}$ are only approximate posteriors because of the loops in the factor graph. Since we need to minimize the BFE given by (9), the distribution function $b_{x,z}(\mathbf{x}, \mathbf{z})$ must be of the form

$$b_{x,z}(\mathbf{x}, \mathbf{z}) = b_x(\mathbf{x})\delta(\mathbf{z} - \mathbf{A}\mathbf{x}), \quad (10)$$

to avoid an infinite value of the KLD, leading to $D[b_{x,z}(\mathbf{x}, \mathbf{z})\|\delta(\mathbf{z} - \mathbf{A}\mathbf{x})] = -H[b_x]$. For BP, the BFE (9) needs to be minimized w.r.t. marginal consistency constraints $q_x(x_i) = b_x(x_i) = q_{x_i}(x_i)$, $q_z(z_j) = q_{z_j}(z_j)$. Given the independent priors for $\mathbf{x}$, $\mathbf{z}$, minimization of the BFE leads to $q_x(\mathbf{x}) = \prod_i q_{x_i}(x_i)$, $q_z(\mathbf{z}) = \prod_j q_{z_j}(z_j)$. Furthermore, the maximization of $H[b_x]$ under marginal constraints leads to $b_x(\mathbf{x}) = \prod_i b_{x_i}(x_i)$. Together with the marginal constraints, this leads to the cancellation of the entropy terms in $\mathbf{x}$ in the BFE, which becomes

$$\sum_{i} D[q_{x_i}(x_i)\|p(x_i)] + \sum_{j} D[q_{z_j}(z_j)\|p(y_j|z_j)] + \sum_{j} H[q_{z_j}(z_j)] \quad (11)$$

which needs to be minimized under the constraint $\mathbf{z} = \mathbf{A}\mathbf{x}$.

## III. GAMP FROM LSL BELIEF PROPAGATION

In reGVAMP [12], [10], extrinsics in the GLM are built from the equivalent Gaussian linear model, which introduces equivalent Gaussian priors from Gaussian posterior approximations and Gaussian extrinsics.

GAMP exploits LSL simplifications of reGVAMP for a random $\mathbf{A}$ with i.i.d. signs which leads to

(i) Gaussianity of extrinsics (also in reGVAMP), and

(ii) independence of marginals (extra w.r.t. reGVAMP).

(ii) leads to the large system simplifications of the variances, avoiding covariance matrix inverses. But also posterior and extrinsic estimates $\hat{\mathbf{x}}$, $\hat{\mathbf{z}}$ and $\mathbf{r}$, $\mathbf{p}$ that are constructed by combining decoupled pieces of information. These estimates are non-linear MMSE and CWCU MMSE estimates in

<figure><img src="image_2.png" /><figcaption>Fig. 1. Factor Graph for the GLM used by GAMP.</figcaption></figure>

general. Extrinsics are not obtained as linear perturbations of corresponding MMSE estimates because those are not necessarily close to each other. Rather the interplay between $\mathbf{x}$ and $\mathbf{z}$ is exploited with perturbations due to the small effect of a single term in $\mathbf{A}$ in the LSL. In both reGVAMP and GAMP, we have:

Gaussian extrinsics: $e_x(\mathbf{x}) = \mathcal{N}(\mathbf{x}; \mathbf{r}, \tau_r)$, $e_z(\mathbf{z}) = \mathcal{N}(\mathbf{z}; \mathbf{p}, \tau_p)$ and

Posterior marginals proportional to: $q_x(\mathbf{x}) \sim p_x(\mathbf{x}) e_x(\mathbf{x})$, $q_z(\mathbf{z}) \sim p_{y|\mathbf{z}}(\mathbf{y}|\mathbf{z}) e_z(\mathbf{z})$ with Gaussian approximations $\mathcal{N}(\mathbf{x}; \hat{\mathbf{x}}, \tau_{\mathbf{x}})$, $\mathcal{N}(\mathbf{z}; \hat{\mathbf{z}}, \tau_{\mathbf{z}})$ resp. (where $\mathbf{r}, \tau_{\mathbf{r}}, \mathbf{p}, \tau_{\mathbf{p}}$ etc. are vectors, e.g. $\mathcal{N}(\mathbf{x}; \mathbf{r}, \tau_{\mathbf{r}})$ is short for $\mathcal{N}(\mathbf{x}; \mathbf{r}, \mathbf{D}_{\tau_{\mathbf{r}}})$, and $\mathbf{D}_{\tau_{\mathbf{r}}} = \text{diag}(\tau_r)$). In [11], it is shown that LSL simplifications of BP lead to the GAMP algorithm (see e.g. [7]) which computes the indicated marginal posteriors that minimize the cost function [8] (to be minimized w.r.t. $q_x$, $q_z$ and $\mathbf{u}$ later)

$$D(q_x \| p_x e_x / Z_x) + D(q_z \| p_z e_z / Z_z). \quad (12)$$

We get per component

$$\min_{q_{x_i}} D(q_{x_i} \| g_{x_i} / Z_{x_i}) \Rightarrow q_{x_i} = g_{x_i} / Z_{x_i}, \quad Z_{x_i} = \int g_{x_i}(x_i) dx_i, \\ - \ln g_{x_i}(x_i) = f_{x_i}(x_i) + \frac{1}{2\tau_{r_i}} [(x_k - r_i)^2 - r_i^2]. \quad (13)$$

The partition function $Z_{x_i}$ acts as cumulant generating function:

$$\tau_{r_i} \frac{\partial \ln Z_{x_i}}{\partial r_i} = \mathbb{E}(x_i|q_{x_i}) = \mathbb{E}(x_i|r_i, \tau_{r_i}) = \hat{x}_i \quad (14)$$

$$\tau_{r_i}^2 \frac{\partial^2 \ln Z_{x_i}}{\partial r_i^2} = \text{var}(x_i|r_i, \tau_{r_i}) = \tau_{x_i}. \quad (14)$$

We also get per component

$$\min_{q_{z_k}} D(q_{z_k} \| g_{z_k} / Z_{z_k}) \Rightarrow q_{z_k} = g_{z_k} / Z_{z_k} \\ Z_{z_k} = \int g_{z_k}(z_k) dz_k, \quad -\ln g_{z_k}(z_k) = \\ f_{z_k}(z_k) + \frac{1}{2\tau_{p_k}^2} [(z_k - p_k)^2 - p_k^2]. \quad (15)$$

The partition function $Z_{z_k}$ acts again as cumulant generating function:

$$-\frac{\partial \ln Z_{z_k}}{\partial s_k} = \mathbb{E}(z_k|q_{z_k}) = \mathbb{E}(z_k|p_k, \tau_{p_k}) = \hat{z}_k \quad (16)$$

$$\frac{\partial^2 \ln Z_{z_k}}{\partial s_k^2} = \text{var}(z_k|p_k, \tau_{p_k}) = \tau_{z_k}. \quad (16)$$

The LSL BP derivation also leads to the following identities

$$Z_z(p, y, \tau_p) = \int p_{y|z}(y|z) e^{-\frac{1}{2\tau_p}(z-p)^2} dz \\ \frac{\partial \ln Z_z}{\partial p} = \frac{Z_z'}{Z_z} = s = \frac{\hat{z}_p - p}{\tau_p}, \quad \hat{z}_p = \frac{1}{Z_z} \int z p_{y|z}(y|z) e^{-\frac{1}{2\tau_p}(z-p)^2} dz \\ \frac{\partial^2 \ln Z_z}{\partial p^2} = -\tau_s = \frac{Z_z''}{Z_z} - (\frac{Z_z'}{Z_z})^2 = -(1 - \tau_z / \tau_p) / \tau_p$$
and updates of the following quantities

$$
\begin{align}
\mathbf{p} &= \mathbf{A} \hat{\mathbf{x}} - \tau_k^p \cdot \mathbf{s}, \quad \tau_p - \mathbf{S} \tau_x \nonumber \\
\mathbf{r} &= \hat{\mathbf{x}} + \tau_r \cdot \mathbf{A}^T \mathbf{s}, \quad \tau_r = 1 \cdot (\mathbf{S}^T \tau_s) \tag{17}
\end{align}
$$

where $\mathbf{S} = \mathbf{A} \cdot \mathbf{A}$ and we use the notations: $\|\mathbf{u}\|_{\tau}^{2} = \sum_{i} u_{i}^{2} / \tau_{i}$, element-wise multiplication as in $\mathbf{s} \cdot \mathbf{\tau}$ and element-wise division as in $\mathbf{1} \cdot / \mathbf{\tau}$, and $\mathbf{1}$ is a vector of ones.

IV. LSL BFE AND EP

After the LSL simplifications [11], the BFE from (11) with
marginal pdf consistency constraints can be seen to become
equivalent to the following LSL-BFE [7], [9] :

$$
\begin{align}
\min_{q_x, q_z, \tau_p, \mathbf{u}} & D[q_x \| p_x] + D[q_z \| p_y | \mathbf{z}] + \frac{1}{2} \sum_k \left[ \frac{\operatorname{var}(z_k | q_z)}{\tau_{p_k}} + \ln(\tau_{p_k}) \right] \nonumber \\
\text{s.t.} & E[\mathbf{z} | q_z] = \mathbf{A} \mathbf{u} \tag{18}
\end{align}
$$

$$
\begin{align*}
E[x|q_x] &= \mathbf{u} \\
\tau_p &= \mathbf{S} \operatorname{var}(x|q_x).
\end{align*}
$$

We will exploit some useful relations

$$
\begin{align}
\forall \tau, \mathbf{c}^T \operatorname{var}(\mathbf{x}|q_x) &= \int \|\mathbf{x} - \mathbf{u}\|_{1/\tau}^2 q_x(\mathbf{x}) \, d\mathbf{x} \tag{19} \\
\sum_k \frac{\operatorname{var}(z_k|q_z)}{\tau_{p_k}} &= \int \|\mathbf{z} - \mathbf{A}\mathbf{u}\|_{\tau_p}^2 q_z(\mathbf{z}) \, d\mathbf{z}. \nonumber
\end{align}
$$

The Lagrangian of (18) becomes

$$
\begin{align}
L = D[q_x \| p_x] &+ D[q_z \| p_y|\mathbf{z}] + \frac{1}{2} \sum_k \left[ \frac{\operatorname{var}(z_k|q_z)}{\tau_{p_k}} + \ln(\tau_{p_k}) \right] \nonumber \\
&+ \lambda_{\mu_z}^T \left( \mathbf{A}\mathbf{u} - \int \mathbf{z} q_z(\mathbf{z}) d\mathbf{z} \right) + \lambda_{\mu_x}^T \left( \mathbf{u} - \int \mathbf{x} q_x(\mathbf{x}) d\mathbf{x} \right) \nonumber \\
&- \frac{1}{2} \lambda_{\tau}^T \left( \tau_p - \mathbf{S} \operatorname{var}(\mathbf{x}|q_x) \right) \tag{20}
\end{align}
$$

The derivatives w.r.t. $q_x$, $q_z$, $\tau_p$, $\mathbf{u}$ become

$$
\begin{align}
\frac{\partial L}{\partial q_x} &= \ln(q_x) - \ln(p_x) - \lambda_{\mu_x}^T \mathbf{x} + \frac{1}{2} \|\mathbf{x} - \mathbf{u}\|_{1/\cdot}^2 \cdot (\mathbf{S}^T \lambda_\tau) \nonumber \\
\frac{\partial L}{\partial q_z} &= \ln(q_z) - \ln(p_y|\mathbf{z}) - \lambda_{\mu_z}^T \mathbf{z} + \frac{1}{2} \|\mathbf{z} - \mathbf{A}\mathbf{u}\|_{\tau_p}^2 \nonumber \\
\frac{\partial L}{\partial \tau_{p_k}} &\propto -\frac{\operatorname{var}(z_k|q_z)}{\tau_{p_k}^2} + \frac{1}{\tau_{p_k}} - \lambda_{\tau_k} \nonumber \\
\frac{\partial L}{\partial \mathbf{u}} &= -\mathbf{A}^T \mathbf{D}_{\tau_p}^{-1} (\hat{\mathbf{z}} - \mathbf{A}\mathbf{u}) + \mathbf{A}^T \lambda_{\mu_z} + \lambda_{\mu_x} \nonumber \\
&\quad - \mathbf{D}_{\mathbf{S}^T \lambda_\tau} (\hat{\mathbf{x}} - \mathbf{u}), \tag{21}
\end{align}
$$

where $\hat{\mathbf{z}} = \mathbb{E}(\mathbf{z}|q_z)$ and $\hat{\mathbf{x}} = \mathbb{E}(\mathbf{x}|q_x)$. Zeroing derivatives:

$$
q_x(\mathbf{x}) \propto p_x(\mathbf{x}) e^{-\frac{1}{2} \|\mathbf{x} - \mathbf{u}\|_{1/\cdot}^2} \cdot (\mathbf{S}^T \lambda_\tau) e^{\lambda_{\mu_x}^T \mathbf{x}} \quad (22)
$$

$$
q_z(\mathbf{z}) \propto p_{y|z}(\mathbf{y}|\mathbf{z}) e^{-\frac{1}{2} \|\mathbf{z} - \mathbf{A}\mathbf{u}\|_{\tau_p}^2} e^{\lambda_{\mu_z}^T \mathbf{z}} \quad (23)
$$

$$
\lambda_{\tau_k} = \frac{1}{\tau_{p_k}} - \frac{\tau_{z_k}}{\tau_{p_k}^2} \quad (24)
$$

$$
\begin{align}
& \left[ \mathbf{A}^T \mathbf{D}_{\tau_p}^{-1} \mathbf{A} + \mathbf{D}_{\mathbf{S}^T \lambda_\tau} \right] \mathbf{u} = \mathbf{A}^T \mathbf{D}_{\tau_p}^{-1} \hat{\mathbf{z}} \tag{25} \\
& \quad + \mathbf{D}_{\mathbf{S}^T \lambda_\tau} \hat{\mathbf{x}} - \mathbf{A}^T \lambda_{\mu_z} - \lambda_{\mu_x} \nonumber
\end{align}
$$

where $\tau_{z_k} = \mathbb{E}[(z_k - \hat{z}_k)^2 | q_z]$. By satisfying the two mean
constraints in (18), the equation (25) becomes

$$
\mathbf{A}^T \lambda_{\mu_z} = -\lambda_{\mu_x} \quad (26)
$$

A solution can be obtained by solving the system of 7 equa-
tions containing (22), (23), (24), (26) and the three constraint
equations in (18).

V. ITERATIVE SOLUTION LEADING TO GAMP

We ignore pdf normalization for simplicity. Furthermore, we
use red symbols to indicate parameters to be updated.

A. Update of $\boldsymbol{\lambda}_{\mu_z}$

Consider (23) and the two mean constraints in (18)

$$
\mathbb{E} \left[ \mathbf{z} | p_{y|z}(\mathbf{y}|\mathbf{z}) e^{-\frac{1}{2} \|\mathbf{z}-\hat{\mathbf{z}}\|_{\tau_p}^2} e^{\lambda_{\mu_z}^{(\text{new})T} \mathbf{z}} \right] \quad (27)
$$

$$
= \mathbb{E} \left[ \mathbf{z} | p_{y|z}(\mathbf{y}|\mathbf{z}) e^{-\frac{1}{2} \|\mathbf{z}-\mathbf{A}\hat{\mathbf{x}}\|_{\tau_p}^2} e^{\lambda_{\mu_z}^T \mathbf{z}} \right] = \hat{\mathbf{z}} \quad (28)
$$

We first use (28) to obtain $\hat{\mathbf{z}}$. Then we use this newly obtained $\hat{\mathbf{z}}$ to update $\lambda_{\mu_z}^{(\text{new})}$, since we need to keep the exponential factor identical in order not to change the mean, i.e.,

$$
e^{-\frac{1}{2} \|\mathbf{z} - \hat{\mathbf{z}}\|_{\tau_p}^2} e^{\lambda_{\mu_z}^{(\text{new})T} \mathbf{z}} = e^{-\frac{1}{2} \|\mathbf{z} - \mathbf{A}\hat{\mathbf{x}}\|_{\tau_p}^2} e^{\lambda_{\mu_z}^T \mathbf{z}}. \quad (29)
$$

If we want to bridge the GAMP from [7] and BFE, we can
denote

$$
p = \mathbf{A}\hat{\mathbf{x}} + \mathbf{D}(\tau_p)\lambda_{\mu_z}. \quad (30)
$$

With definition (30), the expression (28) can be written as

$$
\hat{\mathbf{z}} = \mathbb{E} \left[ \mathbf{z} | p_{y|z}(\mathbf{y}|\mathbf{z}) e^{-\frac{1}{2} \|\mathbf{z} - \mathbf{p}\|_{\tau_p}^2} \right]. \quad (31)
$$

Therefore, the updating for $\lambda_{\mu_z}^{(\text{new})}$ according to (29) becomes

$$
\lambda_{\mu_z}^{(\text{new})} = \mathbf{D}(\tau_p^{-1})(\mathbf{p} - \hat{\mathbf{z}}). \quad (32)
$$

It is now clear that we can relate to the LSL BP GAMP above
(or [7]) if we define

$$
\mathbf{s} = -\lambda_{\mu_z} \cdot \quad (33)
$$

For further use, we also state the computation of $\tau_{\hat{\mathbf{z}}}$ explicitly:

$$
\tau_{\hat{\mathbf{z}}} = \mathbb{E} \left[ (\mathbf{z} - \hat{\mathbf{z}})^2 \middle| p_{y|z}(\mathbf{y}|\mathbf{z}) e^{-\frac{1}{2} \|\mathbf{z} - \mathbf{p}\|_{\tau_p}^2} \right] \quad (34)
$$

$$
= \mathbb{E} \left[ (\mathbf{z} - \hat{\mathbf{z}})^2 \middle| p_{y|z}(\mathbf{y}|\mathbf{z}) e^{-\frac{1}{2} \|\mathbf{z} - \hat{\mathbf{z}}\|_{\tau_p}^2} e^{\lambda_{\mu_z}^{(\text{new})T} \mathbf{z}} \right] \quad (35)
$$

where $\mathbf{z}^{.2}$ denotes element-wise square of vector $\mathbf{z}$. (34) and (35) result in the same solution.

B. Update of $\boldsymbol{\lambda}_{\mu_x}$

According to (26), we can update $\lambda_{\mu_x}$ by

$$
\lambda_{\mu_x} = -\mathbf{A}^T \lambda_{\mu_x}. \quad (36)
$$

To show the relation between this paper and [7], we define

$$
\begin{align}
\tau_r &= 1 \cdot (\mathbf{S}^T \lambda_\tau) \nonumber \\
\mathbf{r} &= \hat{\mathbf{x}}^{\text{old}} + \mathbf{D}_{\tau_r} \mathbf{A}^T \mathbf{s} \tag{37}
\end{align}
$$

Then the updated posterior mean and variance becomes

$$
\hat{\mathbf{x}} = \mathbb{E} \left[ \mathbf{x} | p_{\mathbf{x}}(\mathbf{x}) e^{-\frac{1}{2} \|\mathbf{x} - \mathbf{r}\|_{\tau_r}^2} \right] \quad (38)
$$

$$
\tau_{\hat{\mathbf{x}}} = \mathbb{E} \left[ \mathbf{x} - \hat{\mathbf{x}} \right]^2 |p_{\mathbf{x}}(\mathbf{x}) e^{-\frac{1}{2} \|\mathbf{x} - \mathbf{r}\|_{\tau_r}^2} ,
$$

where we also used the mean constraint for $\mathbf{x}$ in (18).
### C. The update of $\lambda_{\tau}$ and $\tau_p$

The updates of these two variables are quite straightforward. They are already explicitly given by (24) and the variance constraint in (18). To show the relation with GAMP in [7] explicitly, we can define $\tau_s$ and then get

$$ \tau_s = \lambda_{\tau} \text{ from which } \tau_p = S\tau_x, \tau_{s_k} = \frac{1}{\tau_{p_k}} - \frac{\tau_{z_k}}{\tau_{p_k}^2}. \quad (39) $$

## VI. ITERATIVE SOLUTION LEADING TO AMBGAMP

GAMP does not use the extra variable $u$ in (18) (hence uses $u = \hat{x}$) and as result is an algorithm that does not correspond to alternating optimization of a BFE, with the resulting convergence issues. For AMBGAMP, we keep variable $u$, and use (22)-(26) along with the three constraints in (18) as a system of 8 equations to be solved.

### A. Update of $\lambda_{\mu_z}$

Consider (23) and the mean constraint in (18)

$$ \mathbb{E}[\mathbf{z}|p_{y|z}(\mathbf{y}|\mathbf{z})e^{-\frac{1}{2}\|\mathbf{z}-\hat{\mathbf{z}}\|^2_{\tau_p}}e^{\lambda_{\mu_z}^{(\text{new})T}\mathbf{z}}] = \hat{\mathbf{z}} \quad (40) $$

$$ = \mathbb{E}[\mathbf{z}|p_{y|z}(\mathbf{y}|\mathbf{z})e^{-\frac{1}{2}\|\mathbf{z}-\mathbf{A}\mathbf{u}\|_{\tau_p}^2}e^{\lambda_{\mu_z}^T\mathbf{z}}\mathbf{z}] = \hat{\mathbf{z}} \quad (41) $$

To make the connection with AMBGAMP in [7], we define

$$ p = \mathbf{A}u + \tau_p \lambda_{\mu_z} \quad (42) $$

Similar to the previous section, we have the update

$$ \lambda_{\mu_z}^{(\text{new})} = (p - \hat{z}) / \tau_p. \quad (43) $$

Substitute (42) into (43), and we have

$$ \lambda_{\mu_z}^{(\text{new})} = \lambda_{\mu_z} + (\mathbf{A}u - \hat{z}) / \tau_p. \quad (44) $$

If we define $s = -\lambda_{\mu_z}$, (45)

it then follows

$$ s^{(\text{new})} = s + (\hat{z} - \mathbf{A}u) / \tau_p. \quad (46) $$

For the convenience of the further discussion, we write the update for the posterior mean and variance of $\mathbf{z}$

$$ \hat{\mathbf{z}} = \mathbb{E}[\mathbf{z}|p_{y|z}(\mathbf{y}|\mathbf{z})e^{-\frac{1}{2}\|\mathbf{z}-\mathbf{p}\|^2_{\tau_p}}] \quad (47) $$

$$ \tau_{\hat{\mathbf{z}}} = \mathbb{E}\left[(\mathbf{z} - \hat{\mathbf{z}})^2 | p_{y|z}(\mathbf{y}|\mathbf{z})e^{-\frac{1}{2}\|\mathbf{z}-\mathbf{p}\|^2_{\tau_p}}\right]. \quad (48) $$

### B. Update of $\lambda_{\mu_{\infty}}$

We can use (26) and (45) to obtain the update

$$ \lambda_{\mu_{\infty}} = -\mathbf{A}^T \lambda_{\mu_z} = \mathbf{A}^T s \quad (49) $$

If we define (and note that $\lambda_{\tau} = \tau_s$)
$$ \tau_r = 1/(\mathbf{S}^T \lambda_{\tau}), \quad \mathbf{r} = \mathbf{u} + \tau_r \cdot \lambda_{\mu_{\infty}} = \mathbf{u} + \tau_r \cdot (\mathbf{A}^T \mathbf{s}) \lambda_{\mu_{\infty}}, \quad (50) $$

we have the explicit update for $\hat{x}$ and $\tau_{\hat{x}}$:

$$ \hat{x} = \mathbb{E}[\mathbf{x}|p_{\mathbf{x}}(\mathbf{x})e^{-\frac{1}{2}\|\mathbf{x}-\mathbf{r}\|^2_{\tau_r}}] \quad (51) $$

$$ \tau_{\hat{x}} = \mathbb{E}\left[(\mathbf{x} - \hat{x})^2 | p_{\mathbf{x}}(\mathbf{x})e^{-\frac{1}{2}\|\mathbf{x}-\mathbf{r}\|^2_{\tau_r}}\right]. \quad (51) $$

### C. Update of $u$

By combining (25) and (26), we get the solution

$$ u = [\mathbf{A}^T \mathbf{D}_{\tau_p}^{-1} \mathbf{A} + \mathbf{D}_{\tau_r}^{-1}]^{-1} (\mathbf{A}^T \mathbf{D}_{\tau_p}^{-1} \hat{\mathbf{z}} + \mathbf{D}_{\tau_r}^{-1} \hat{\mathbf{x}}). \quad (52) $$

AMBGAMP actually updates $u$ by applying SGD with step-size by linesearch to the quadratic cost function that (52) is the solution of. The update of $\lambda_{\tau}$ and $\tau_p$ are identical to the updates in GAMP in (39).

## VII. CONCLUDING REMARKS

In this paper, we have shown that it is possible to derive the convergent AMBGAMP algorithm by analyzing the KKT conditions for optimizing the LSL BFE. And this while avoiding the quadratic augmentation terms of the Method of Moments, which require a very particular choice in their weights, and circumventing the ADMM-style update of a Lagrange multiplier. This is thanks to the introduction of the auxiliary variable $u$ in the mean consistency constraints. This $u$ is optimized to minimize the BFE equivalent in (12) and can be interpreted to be a MMSE estimate of an equivalent underlying Gaussian linear model. On the other hand, another solution to the LSL BFE, which eliminates $u$ via $u = \hat{x}$, leads to GAMP and corresponds to the original LSL BP based derivation, optimizing BFE with LSL approximations. Hence we have reconciled these seemingly different approaches. The variance predictions in (AMB)GAMP are based on a sign i.i.d. model for $\mathbf{A}$, which leads to decorrelation and Gaussianity after multiplication of a vector with $\mathbf{A}$ or $\mathbf{A}^T$, similar to spreading and despreading in CDMA. Another somewhat popular model for $\mathbf{A}$ is the Right Rotationally Invariant class, in which (only) the right singular vectors of $\mathbf{A}$ are modeled as random, and in particular as Haar distributed. This is the motivation for Vector AMP (VAMP) [13]. To keep complexity low however, VAMP has to restrict diagonal covariances to multiples of identity, which e.g. is not useful for Sparse Bayesian Learning [14]. GAMP-style low complexity algorithms can be derived also, but they require some correction terms in the variance predictions, stemming from the Haar distribution [15], [16].

## VIII. ACKNOWLEDGEMENTS

EURECOM's research is partially supported by its industrial members: ORANGE, BMW, SAP, iABG, Norton LifeLock, by the Franco-German projects CellFree6G and 5G-OPERA, by the EU H2030 project CONVERGE, and by a Huawei France funded Chair towards Future Wireless Networks.

## REFERENCES

[1] D. Zhang, X. Song, W. Wang, G. Fettweis, and X. Gao, "Unifying Message Passing Algorithms Under the Framework of Constrained Bethe Free Energy Minimization," *IEEE Transactions on Wireless Communications*, vol. 20, no. 7, 2021.

[2] M. J. Wainwright, M. I. Jordan et al., "Graphical Models, Exponential Families, and Variational Inference," *Foundations and Trends® in Machine Learning*, vol. 1, no. 1–2, 2008.

[3] T. Minka et al., "Divergence Measures and Message Passing," Citeseer, Tech. Rep., 2005.

[4] T. Heskes, M. Opper, W. Wiegerinck, O. Winther, and O. Zoeter, "Approximate Inference Techniques with Expectation Constraints," *Journal of Statistical Mechanics: Theory and Experiment*, vol. 2005, no. 11, 2005.
[5] Q. Zou and H. Yang, "A Concise Tutorial on Approximate Message Passing," arXiv preprint arXiv:2201.07487, 2022.

[6] S. Rangan, P. Schniter, E. Riegler, A. K. Fletcher, and V. Cevher, "Fixed Points of Generalized Approximate Message Passing with Arbitrary Matrices," *IEEE Transactions on Information Theory*, vol. 62, no. 12, 2016.

[7] C. Kurisummoottil Thomas, Z. Zhao, and D. Slock, "Towards Convergent Approximate Message Passing by Alternating Constrained Minimization of Bethe Free Energy," in *IEEE Information Theory Workshop (ITW)*, Saint Malo, France, 2023.

[8] S. Rangan, A. Fletcher, P. Schniter, and U. Kamilov, "Inference for Generalized Linear Models via Alternating Directions and Bethe Free Energy Minimization," *IEEE Trans. Info. Theory*, Jan. 2017.

[9] Z. Zhao and D. Slock, "Bethe Free Energy and Extrinsics in Approximate Message Passing," in *IEEE Asilomar Conf. Signals, Systems and Computers*, 2023.

[10] Z. Zhao, F. Xiao, and D. Slock, "Vector approximate message passing for not so large N.I.I.D. generalized I/O linear models," in *IEEE Int'l Conf. Acoustics, Speech and Sig. Proc. (ICASSP)*, Seoul, 2024.

[11] ——, "Extrinsics and Linearized Component-Wise Conditionally Unbiased MMSE Estimation as in GAMP," in *IEEE Asilomar Conf. Signals, Systems and Computers*, 2024.

[12] ——, "Approximate Message Passing for Not So Large niid Generalized Linear Models," in Proc. Int'l Workshop Signal Processing Advances in Wireless Comm's (SPAWC), Sept. 2023.

[13] S. Rangan, P. Schniter, and A. K. Fletcher, "Vector Approximate Message Passing," *IEEE Trans. On Info. Theo.*, Oct. 2019.

[14] C. K. Thomas and D. Slock, "SAVE - Space alternating variational estimation for sparse Bayesian learning," in *IEEE Data Science Workshop*, June 2018.

[15] Z. Zhao and D. Slock, "Variance Predictions in VAMP/UAMP with Right Rotationally Invariant Measurement Matrices for niid Generalized Linear Models," in *European Sig. Proc. Conf. (EUSIPCO)*, Helsinki, Finland, 2023.

[16] ——, "Improved Variance Predictions in Approximate Message Passing," in *IEEE Int'l Workshop Machine Learning and Sig. Proc. (MLSP)*, Rome, Italy, 2023.