# Non-Circular Combined Expectation Propagation and Variational Bayes for Semi-Blind Channel Estimation in Cell-Free Systems

Zilu Zhao, Dirk Slock

► To cite this version:

Zilu Zhao, Dirk Slock. Non-Circular Combined Expectation Propagation and Variational Bayes for Semi-Blind Channel Estimation in Cell-Free Systems. ICNC 2025, International Conference on Computing, Networking and Communications, IEEE, Feb 2025, Honolulu, United States. pp.775-779, 10.1109/ICNC64010.2025.10993718 . hal-05076605

HAL Id: hal-05076605

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

Submitted on 21 May 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
# Non-Circular Combined Expectation Propagation and Variational Bayes for Semi-Blind Channel Estimation in Cell-Free Systems

Zilu Zhao, Dirk Slock

Communication Systems Department, EURECOM, France
zilu.zhao@eurecom.fr, dirk.slock@eurecom.fr

**Abstract**—In this work, we investigate uplink communication in Semi-Blind Cell-Free (CF) Massive Multiple-Input Multiple-Output (MaMIMO) systems. A key challenge in CF MaMIMO systems is pilot contamination, which arises when multiple user terminals (UTs) share the same pilot sequence due to an imbalance between the number of UTs and the length of the pilot sequence. Semi-blind approaches have been proposed to address this issue by enabling access points (APs) to jointly estimate both the channel and user data. This joint estimation, however, leads to a bilinear problem. To develop a tractable algorithm that exploits the finite alphabet of the user data, we analyze the constrained Bethe Free Energy (BFE) of the bilinear system and propose a message-passing algorithm based on minimizing the constrained BFE. We decompose complex quantities into in-phase and quadrature components, removing circular Gaussian constraints. Additionally, we find that replacing normal covariance constraints with average covariance constraints significantly reduces computational complexity. Simulation results indicate that this simplification incurs negligible performance loss.

## I. INTRODUCTION

In Cell-Free (CF) Massive Multiple-Input Multiple-Output (MaMIMO) systems, user terminals (UTs) are simultaneously served by all access points (APs) in a given region. A significant challenge in CF MaMIMO systems is pilot contamination, which occurs when the number of users exceeds the length of the pilot sequences. Consequently, APs cannot estimate the channel solely based on pilot sequences. To address this issue, semi-blind channel estimation is employed [1]. In Semi-Blind settings, APs jointly estimate the channel and user data based on received signals and limited pilot sequences.

### A. Prior Work

Bayesian estimation in semi-blind structures holds significant potential [1], but it also presents challenges due to high-dimensional and intractable integrals. Message-passing algorithms, particularly Expectation Propagation (EP) [2] and Belief Propagation (BP) [3], are widely used in Bayesian estimation. Both EP and BP assume a factored joint probability density function and simplify high-complexity global inference problems into manageable local inference tasks. EP further reduces complexity by approximating the factors of the joint pdf with simpler forms, such as Gaussian distributions.

Variable-Level EP (VL-EP) was introduced for Gaussian input data by combining Expectation-Maximization (EM) with EP [4]. To improve its convergence properties, hybrid EM-EP and loop-free EM-EP algorithms were proposed in [5]. However, these approaches are not designed to handle user symbols from finite alphabets.

1) **Expectation Propagation for Gaussian Mixture Models:** The bilinear combination of a Gaussian distribution (e.g., channel distribution) and a discrete distribution (e.g., input data distribution) leads to a Gaussian Mixture Model. To address the limitations of VL-EP, a distributed bilinear-EP algorithm was proposed in [6], which adopts a brute-force approach to inference over finite alphabets, avoiding high-dimensional computations by considering only one data symbol at a time. Inspired by [6] and [7], the authors of [8] proposed a simplified decentralized bilinear-EP algorithm.

2) **Bethe Free Energy:** The Bethe Free Energy (BFE) is another powerful tool for Bayesian inference. It represents the variational energy of a factored joint pdf under a specific trial distribution, whose form is determined by the factorization scheme of the joint pdf. It has been demonstrated [9] that various message-passing algorithms, such as EP and BP, can be derived by optimizing BFE under different constraints on the trial distribution.

Hybrid Vector Message Passing (HVMP) [10] was proposed based on BFE optimization, introducing a mean-field constraint for the bilinear factor. However, this method does not account for finite alphabets and entirely neglects the correlation between the channel and data.

### B. Main Contributions

We propose a low-complexity algorithm for semi-blind channel and data estimation, leveraging a framework based on BFE constrained optimization. To effectively handle posterior interference, we introduce an auxiliary variable, which enables a more tractable optimization process. Unlike prior works such as [6] and [8], our method treats the entire data sequence of a single user as a single atomic variable, significantly simplifying the estimation procedure. To address non-analytical integrals that arise during posterior estimation, we incorporate mean-field assumptions into the belief factors by including delta functions. This approximation streamlines the derivations and reduces computational complexity. Furthermore, we derive the algorithm by separating complex quantities into their in-phase and quadrature components.

Lastly, we introduce a covariance averaging operation to further minimize complexity while maintaining robust performance. This combination of methods ensures a computationally efficient and scalable approach to semi-blind estimation.

## II. SYSTEM MODEL

We examine the uplink cell-free semi-blind network containing K single-antenna user terminals (UTs) and L access points (APs). Each AP is equipped with M antennas. The received signals of the l-th AP is

$$[\hat{\mathbf{Y}}_{p,l} \quad \hat{\mathbf{Y}}_l] = \hat{\mathbf{H}}_l [\hat{\mathbf{X}}_p^\mathrm{T} \quad \hat{\mathbf{X}}^\mathrm{T}] + [\hat{\mathbf{V}}_{p,l} \quad \hat{\mathbf{V}}_l] \in \mathbb{C}^{M \times (P+T)}, \quad (1)$$

where $\hat{\mathbf{H}}_l \in \mathbb{C}^{M \times K}$ models the channel matrix, $\hat{\mathbf{X}}_p \in \mathbb{C}^{P \times K}$ models the pilot sequences of all the users, $\hat{\mathbf{X}} \in \mathbb{C}^{T \times K}$ models the transmitted data sequences.

The complex model can be transformed into in-phase and quadrature representation

$$\begin{bmatrix} \Re[\hat{\mathbf{Y}}_{p,l}] & \Re[\hat{\mathbf{Y}}_l] \\ \Im[\hat{\mathbf{Y}}_{p,l}] & \Im[\hat{\mathbf{Y}}_l] \end{bmatrix} = \begin{bmatrix} \Re[\hat{\mathbf{H}}_l] & -\Im[\hat{\mathbf{H}}_l] \\ \Im[\hat{\mathbf{H}}_l] & \Re[\hat{\mathbf{H}}_l] \end{bmatrix} \begin{bmatrix} \Re[\hat{\mathbf{X}}_p^\mathrm{T}] & \Re[\hat{\mathbf{X}}^\mathrm{T}] \\ \Im[\hat{\mathbf{X}}_p^\mathrm{T}] & \Im[\hat{\mathbf{X}}^\mathrm{T}] \end{bmatrix} + \mathbf{V}, \quad (2)$$

$$\text{where } \mathbf{V} = \begin{bmatrix} \Re[\hat{\mathbf{V}}_{p,l}] & \Re[\hat{\mathbf{V}}_l] \\ \Im[\hat{\mathbf{V}}_{p,l}] & \Im[\hat{\mathbf{V}}_l] \end{bmatrix}.$$
For simplicity, we define the following notations

$$
\begin{align*}
\mathbf{Y}_{p,l} &= \begin{bmatrix} \Re[\hat{\mathbf{Y}}_{p,l}] \\ \Im[\hat{\mathbf{Y}}_{p,l}] \end{bmatrix}; \quad \mathbf{Y}_l = \begin{bmatrix} \Re[\hat{\mathbf{Y}}_l] \\ \Im[\hat{\mathbf{Y}}_l] \end{bmatrix}; \quad \mathbf{H}_l = \begin{bmatrix} \Re[\hat{\mathbf{H}}_l] \\ \Im[\hat{\mathbf{H}}_l] \end{bmatrix}; \\
\mathbf{x}_{k,\Re} &= \Re[\hat{\mathbf{x}}_k]; \quad \mathbf{x}_{k,\Im} = \Im[\hat{\mathbf{x}}_k]; \quad \mathbf{X}_k^\mathrm{T} = \begin{bmatrix} \mathbf{x}_{k,\Re}^\mathrm{T} \\ \mathbf{x}_{k,\Im}^\mathrm{T} \end{bmatrix}; \quad \mathbf{x}_k = \mathrm{vec}[\mathbf{X}_k^\mathrm{T}]
\end{align*}
$$

We assume that $\mathbf{H}_l$ is column-wise independent, with the distri-
bution of its $k$-th column given by $\mathbf{h}_{lk} \sim \mathcal{N}(\mathbf{h}_{lk}|\mathbf{0}, \Xi \mathbf{h}_{lk})$. Let
$\mathbf{S} = \{\mathbf{s}_1, \dots, \mathbf{s}_{|\mathbf{S}|}\}$ denote the set of symbol constellations. Each data
symbol is assumed to be independent of others, with a distribution
$p(\mathbf{x}_{kt})$, where $\mathbf{x}_{kt}$ represents the $t$-th column of $\mathbf{X}_k$.

The power of each pilot sequence (complex) is assumed to be $P\sigma_x^2$.
The noise is modeled as additive white Gaussian noise (AWGN),
with each entry in $\mathbf{V}$ having power $\sigma_v^2/2$. For simplicity, we define
$\mathbf{C}_v = \frac{\sigma_v^2}{2}\mathbf{I}$.

A. Orthogonal Pilots

When orthogonal pilots are used, we correlate the received pilot
signals $Y_{p,l}$ with the $g$-th pilot sequence $\tilde{\mathbf{x}}_{p,g}$ (not to confuse with
the pilot sequence of the $g$-th user) to obtain the correlated version
of the received pilot signals $\tilde{\mathbf{y}}_{p,lg}$:

$$
\tilde{\mathbf{y}}_{p,lg} = \hat{\mathbf{Y}}_{p,l} \tilde{\mathbf{x}}_{p,g}^* = P \sigma_x^2 \hat{\mathbf{H}}_{lG_g} \mathbf{1}_{|G_g|} + \tilde{\mathbf{v}}_{p,lg}, \quad (3)
$$

where we use $G_g$ to denote the UTs groups using the g-th pilot
sequence. The columns of $\hat{\mathbf{H}}_{lG_g}$ are composed of the complex
channel coefficients corresponding to the users using the g-th pilot,
i.e., $\hat{\mathbf{h}}_{lk}$ is a column of $\hat{\mathbf{H}}_{lG_g}$ if $\hat{\mathbf{x}}_{p,k} = \hat{\mathbf{x}}_{p,g}$. We denote $\tilde{\mathbf{v}}_{p,lg} =$
$\hat{\mathbf{V}}_{p,l} \tilde{\mathbf{x}}_{p,g}^*$ which is the transformed noise following a distribution
$\mathcal{CN}(\mathbf{0}, \sigma_x^2 \sigma_v^2 P \mathbf{I}_M)$. The complex expression is transformed into in-
phase and quadrature expression. Denote

$$
\tilde{\mathbf{y}}_{p,lg} = \begin{bmatrix} \Re[\tilde{\mathbf{y}}_{p,lg}] \\ \Im[\tilde{\mathbf{y}}_{p,lg}] \end{bmatrix}; \quad \tilde{\mathbf{v}}_g = \begin{bmatrix} \Re[\tilde{\mathbf{v}}_g] \\ \Im[\tilde{\mathbf{v}}_g] \end{bmatrix} \quad (4)
$$

Thus, we have

$$
\tilde{\mathbf{y}}_{p,lg} = P \sigma_x^2 \sum_{k \in G_g} \mathbf{h}_{lk} + \tilde{\mathbf{v}}_g, \quad (5)
$$

where $\tilde{\mathbf{v}}_g \sim \mathcal{N}(\mathbf{v}_g | \mathbf{0}_{2M}, \frac{\sigma_v^2 \sigma_x^2 P}{2} \mathbf{I}_{2M})$. For simplicity, denote $\mathbf{C}_{v,p} = \frac{\sigma_v^2 \sigma_x^2 P}{2} \mathbf{I}_{2M}$ and define $\mathbf{H}_{lG_g}$ to be a matrix whose columns are $\mathbf{h}_{lk}$ such that $k \in G_g$ and its vectorization $\mathbf{h}_{lG_g} = \mathrm{vec}(\mathbf{H}_{lG_g})$.

B. Factored Joint Distribution

Introduce an auxiliary variable

$$
\mathbf{Z}_{lk} = \tilde{\mathbf{H}}_{lk} \mathbf{X}_{k}^{\mathrm{T}} = \mathbf{h}_{lk} \mathbf{x}_{\Re,k}^{\mathrm{T}} + \mathbf{S} \mathbf{h}_{lk} \mathbf{x}_{\Im,k}^{\mathrm{T}} \quad (6)
$$

where

$$
\mathbf{S} = \begin{bmatrix} \mathbf{I}_M & -\mathbf{I}_M \end{bmatrix}; \quad \tilde{\mathbf{H}}_{lk} = \begin{bmatrix} \mathbf{h}_{lk} & \mathbf{S} \mathbf{h}_{lk} \end{bmatrix}. \quad (7)
$$

The vectorization $\mathbf{Z}_{lk}$ can be represented as $\mathbf{z}_{lk} = \mathrm{vec}[\mathbf{Z}_{lk}]$.
Therefore, the likelihood of $\mathbf{Z}_{lk}$ is captured by Dirac function
$p(\mathbf{Z}_{lk}|\mathbf{h}_{lk}, \mathbf{x}_k) = \delta(\mathbf{Z}_{lk} - \mathbf{h}_{lk}\mathbf{x}_{\Re,k}^\mathrm{T} - \mathbf{S}\mathbf{h}_{lk}\mathbf{x}_{\Im,k}^\mathrm{T})$. The joint proba-
bility density function (PDF) can be derived as

$$
\begin{align}
& p(\hat{\mathbf{Y}}_{p,\{l\}}, \hat{\mathbf{Y}}_{\{l\}}, \hat{\mathbf{Z}}_{\{l\}\{k\}}, \hat{\mathbf{H}}_{\{l\}}, \hat{\mathbf{X}}) \nonumber \\
& = \prod_{l} p(\mathbf{Y}_l | \mathbf{Z}_{\{l\}\{k\}}) \prod_{l} \prod_{k} p(\mathbf{Z}_{lk} | \mathbf{h}_{lk}, \mathbf{x}_k) \tag{8} \\
& \prod_{l} \prod_{g} p(\tilde{\mathbf{y}}_{p,lg}, \mathbf{H}_{lG_g}) \prod_{k} p(\mathbf{x}_k). \nonumber
\end{align}
$$

For simplicity, we define

$$
\begin{align}
f_{z_l}(\mathbf{z}_{l\{k\}}) &\propto p(\mathbf{Y}_l|\mathbf{Z}_{l\{k\}}); \quad f_{h_{lG_g}}(\mathbf{h}_{lG_g}) \propto p(\tilde{\mathbf{y}}_{p,lg}, \mathbf{H}_{lG_g}) \tag{9} \\
f_{x_k}(\mathbf{x}_k) &= p(\mathbf{x}_k); \quad f_{\delta_{lk}}(\mathbf{z}_{lk}, \mathbf{h}_{lk}, \mathbf{x}_k) \propto p(\mathbf{Z}_{lk}|\mathbf{h}_{lk}, \mathbf{x}_k). \nonumber
\end{align}
$$

The factorization given by (8) admits a factor graph [3]. We denote
$\mathbb{F} = \{f_{z_l}, f_{h_{lG_g}}, f_{x_k}, \delta_{lk}\}$ as the set of all factor nodes and $\mathbb{V} = \{\mathbf{z}_{lk}, \mathbf{h}_{lk}, \mathbf{x}_k\}$ as the set of all variable nodes.

III. BETHE FREE ENERGY OPTIMIZATION FRAMEWORK

Bethe free energy is the approximated variational free energy between
the true probability (8) and a constrained Bethe approximation trial
function. For a given factored pdf p, its trial pdf b is obtained by:

$$
p(\boldsymbol{\theta}) \propto \prod_{\alpha} f_{\alpha}(\boldsymbol{\theta}_{\alpha}) \Rightarrow b(\boldsymbol{\theta}) = \frac{\prod_{\alpha} b_{f_{\alpha}}(\boldsymbol{\theta}_{\alpha})}{\prod_{i} b_{\theta_{i}}(\boldsymbol{\theta}_{i})^{|N_{i}|-1}}, \quad (10)
$$

s.t.

$$
\forall \alpha, \theta_i \in \theta_{\alpha}, \quad \int b_{f_{\alpha}}(\theta_{\alpha}) d\theta_{\tilde{i}} = b_{\theta_i}(\theta_i) \quad (11)
$$

where $|N_i|$ denotes the number of factors $f_\alpha$ that contain $\theta_i$ and $\theta_{\tilde{i}}$
denotes all the variables except $\theta_i$.

With (10)-(11), the BFE can be obtained as

$$
\text{BFE} = D[b(\boldsymbol{\theta})] \prod_{\alpha} f_{\alpha}(\boldsymbol{\theta}_{\alpha}) = \sum_{\alpha} D(b_{f_{\alpha}} \| f_{\alpha}) + \sum_{i} (|N_i| - 1) H(b_{\theta_i}), \quad (12)
$$

where we define $D(b||q) = \int b(\theta) \ln \frac{b(\theta)}{q(\theta)} d\theta$, and $H(\cdot)$ as entropy.
It is worth noticing that (12) only holds if the factorization (10) is
loop-free and strict constraints (11) are applied. Otherwise, (12) is
only an approximation.

A. Bethe Approximation with Constraints

Following [9], the BFE of (8) is:

$$
\begin{align}
\text{BFE} = & \sum_l D[b_{f_{z_l}}(\mathbf{z}_{l\{k\}}) \| f_{z_l}(\mathbf{z}_{l\{k\}})] \nonumber \\
& + \sum_{l,g} D[b_{f_{h_{lG_g}}}(\mathbf{h}_{lG_g}) \| f_{h_{lG_g}}(\mathbf{h}_{lG_g})] + \sum_k D[b_{f_{x_k}}(\mathbf{x}_k) \| f_{x_k}(\mathbf{x}_k)] \nonumber \\
& + \sum_{l,k} D[b_{\delta_{lk}}(\mathbf{z}_{lk}, \mathbf{h}_{lk}, \mathbf{x}_k) \| f_{\delta_{lk}}(\mathbf{z}_{lk}, \mathbf{h}_{lk}, \mathbf{x}_k)] + \sum_{l,k} H[b_{\mathbf{z}_{lk}}(\mathbf{z}_{lk})] \nonumber \\
& + \sum_{l,k} H[b_{\mathbf{h}_{lk}}(\mathbf{h}_{lk})] + \sum_k L \cdot H[b_{\mathbf{x}_k}(\mathbf{x}_k)]. \tag{13}
\end{align}
$$

where all the factor-level beliefs $b_{f_{z_l}}, b_{\delta_{z_{lk}}}, b_{f_{h_{lG_g}}}, b_{f_{x_k}}$, and
variable-level beliefs $b_{h_{lk}}, b_{z_{lk}}, b_{x_k}$ are proper distributions normal-
ized to one. Furthermore, to make all these factors consistent, the
variable-level beliefs must be the marginal distribution of the factor-
level beliefs. For all $l \in [1, L], k \in [1, K]$, the constraints for the $\mathbf{x}_k$
are

$$
\int b_{\delta_{lk}}(\mathbf{z}_{lk}, \mathbf{h}_{lk}, \mathbf{x}_k) d\mathbf{z}_{lk} d\mathbf{h}_{lk} = b_{\mathbf{x}_k}(\mathbf{x}_k) \quad (14)
$$

$$
b_{f_{x_k}}(\mathbf{x}_k) = b_{x_k}(\mathbf{x}_k). \quad (15)
$$

However, satisfying the strict constraints of $\mathbf{h}_{lk}$ and $\mathbf{z}_{lk}$ will lead
to an intractable problem. Therefore, we relax the strict constraints
to first and second-order moment constraints (specifically, mean and
covariance constraints). W.l.o.g., we denote those sufficient statistics
as $\phi_{\mathbf{h}_{lk}}(\mathbf{h}_{lk})$, $\phi_{\mathbf{z}_{lk}}(\mathbf{z}_{lk})$

$$
\mathbb{E}_{b_{f_{z_l}}} [\phi_{z_{lk}}(\mathbf{z}_{lk})] = \mathbb{E}_{b_{z_{lk}}} [\phi_{z_{lk}}(\mathbf{z}_{lk})] \quad (16)
$$

$$
\mathbb{E}_{\delta_{lk}}[\phi_{z_{lk}}(\mathbf{z}_{lk})] = \mathbb{E}_{b_{z_{lk}}}[\phi_{z_{lk}}(\mathbf{z}_{lk})] \quad (17)
$$

$$
\mathbb{E}_{b_{f_{h_{lG_g}}}}[\phi_{h_{lk}}(\mathbf{h}_{lk})] = \mathbb{E}_{b_{h_{lk}}}[\phi_{h_{lk}}(\mathbf{h}_{lk})] \quad (18)
$$

$$
\mathbb{E}_{b_{\delta_{lk}}}[\phi_{h_{lk}}(\mathbf{h}_{lk})] = \mathbb{E}_{b_{h_{lk}}}[\phi_{h_{lk}}(\mathbf{h}_{lk})] \quad (19)
$$
Moreover, to make the further derivation tractable with finite input **X**, we only consider the average covariance constraints of elements within every size-2M block $\forall t \in [1, T]$, $[\mathbf{z}_{lk}]_{2M(t-1)+1:2Mt}$.

### B. Bethe Free Energy Optimization

The optimization criteria can be concluded by

$$ \begin{gathered} \min_{b} \text{BFE} \\ \text{s.t.} \quad (14) \sim (19). \end{gathered} \quad (20) $$

We observe the term $D[b_{\delta_{lk}}(\mathbf{z}_{lk}, \mathbf{h}_{lk}, \mathbf{x}_k)\|\delta(\mathbf{Z}_{lk} - \mathbf{h}_{lk}\mathbf{x}_k^T)]$ in (13). Since we need to minimize the BFE, the posterior factor $b_{\delta_{lk}}$ must contain the factor $\delta(\mathbf{Z}_{lk} - \mathbf{h}_{lk}\mathbf{x}_k^T)$ to avoid infinity BFE value. In order to have an analytical algorithm, we use the following mean-field approximation for the joint belief $b_{\delta_{lk}}$:

$$ b_{\delta_{lk}}(\mathbf{z}_{lk}, \mathbf{h}_{lk}, \mathbf{x}_k) = b_{\delta_{h,lk}}(\mathbf{h}_{lk}) b_{\delta_{\omega,lk}}(\mathbf{x}_k) p(\mathbf{Z}_{lk}|\mathbf{h}_{lk}, \mathbf{x}_k), \quad (21) $$

where the belief $b_{\delta_{h,lk}}$ and $b_{\delta_{\omega,lk}}$ are beliefs normalized to one. By using Lagrangian method, we can obtain the following message-passing style system of equations along with (21):

$$ b_{f\mathbf{z}_l}(\mathbf{z}_l\{k\}) = p(\mathbf{Y}_l|\mathbf{z}_l\{k\}) \prod_k \mu_{\mathbf{z}_{lk}; f\mathbf{z}_l}(\mathbf{z}_{lk}) \quad (22) $$

$$ b_{f\mathbf{h}_{lG_g}}(\mathbf{h}_{lG_g}) = p(\tilde{\mathbf{y}}_{p,lg}, \mathbf{h}_{lG_g}) \prod_{k \in G_g} \mu_{\mathbf{h}_{lk}; f\mathbf{h}_{lG_g}}(\mathbf{h}_{lk}) \quad (23) $$

$$ b_{f\mathbf{x}_k}(\mathbf{x}_k) = p(\mathbf{x}_k)\mu_{\mathbf{x}_k; f\mathbf{x}_k}(\mathbf{x}_k) \quad (24) $$

$$ b_{\delta_{h,lk}}(\mathbf{h}_{lk}) = \mu_{\mathbf{h}_{lk}; \delta_{lk}}(\mathbf{h}_{lk}) e^{\int b_{\delta_{\omega,lk}}(\mathbf{x}_k) \ln \mu_{\mathbf{z}_{lk}; \delta_{lk}}(\mathrm{vec}(\mathbf{h}_{lk}\mathbf{x}_k^T)) d\mathbf{x}_k} \quad (25) $$

$$ b_{\delta_{\omega,lk}}(\mathbf{x}_k) = \mu_{\mathbf{x}_k; \delta_{lk}}(\mathbf{x}_k) e^{\int b_{\delta_{h,lk}}(\mathbf{h}_{lk}) \ln \mu_{\mathbf{z}_{lk}; \delta_{lk}}(\mathrm{vec}(\mathbf{h}_{lk}\mathbf{x}_k^T)) d\mathbf{h}_{lk}} \quad (26) $$

$$ b_{\mathbf{z}_{lk}}(\mathbf{z}_{lk}) = \mu_{\mathbf{z}_{lk}; f\mathbf{z}_l}(\mathbf{z}_{lk}) \mu_{\mathbf{z}_{lk}; \delta_{lk}}(\mathbf{z}_{lk}) \quad (27) $$

$$ b_{\mathbf{h}_{lk}}(\mathbf{h}_{lk}) = \mu_{\mathbf{h}_{lk}; f\mathbf{h}_{lG_g}}(\mathbf{h}_{lk}) \mu_{\mathbf{h}_{lk}; \delta_{lk}}(\mathbf{h}_{lk}) \quad (28) $$

$$ b_{\mathbf{x}_k}(\mathbf{x}_k) = [\mu_{\mathbf{x}_k; f\mathbf{x}_k}(\mathbf{x}_k) \prod_l \mu_{\mathbf{x}_k; \delta_{lk}}(\mathbf{x}_k)]^{1/L}, \quad (29) $$

The equations (21)~(26) describes the factor level beliefs while (27)~(29) are variable level beliefs. For all $f \in \mathbb{F}$, $\theta \in \mathbb{V}$, we interpret $\mu_{\theta;f}$ as the variable to factor message. Furthermore, we can define the factor to variable messages such that the following relation holds [11]

$$ \forall f \in N(\theta), \mu_{\theta;f}(\theta) = \prod_{f' \in N(\theta)/\{f\}} \mu_{f';\theta}(\theta), \quad (30) $$

where $N(\theta)$ denotes the neighborhood around $\theta$ in the corresponding factor graph. Thus, (29) can be rewritten into the message passing form

$$ b_{\mathbf{x}_k}(\mathbf{x}_k) = \mu_{f\mathbf{x}_k}; \mathbf{x}_k(\mathbf{x}_k) \prod_l \mu_{\delta_{lk}}; \mathbf{x}_k(\mathbf{x}_k) \quad (31) $$

Since the sufficient statistics we consider here are first and second-order moments, the messages $\mu_{f\mathbf{h}_{lG_g}}; \mathbf{h}_{lk}$, $\mu_{\delta_{lk}; \mathbf{h}_{lk}}$, $\mu_{f\mathbf{z}_l}; \mathbf{z}_{lk}$ and $\mu_{\delta_{lk}; \mathbf{z}_{lk}}$ are all (unnormalized) Gaussian distributions.

## IV. ALGORITHM DERIVATION

### A. Message from $\delta_{lk}$ to $\mathbf{x}_k$

By satisfying the constraint (14) between (26) and (31), we obtain the iterative updating expression for the message from $\delta_{lk}$ to $\mathbf{x}_k$:

$$ \mu_{\delta_{lk}; \mathbf{x}_k}(\mathbf{x}_k) = e^{\int b_{\delta_{h,lk}}(\mathbf{h}_{lk}) \ln \mu_{\mathbf{z}_{lk}; \delta_{lk}}[\mathbf{A}(\mathbf{h}_{lk})\mathbf{x}_k] d\mathbf{h}} \quad (32) $$

where $\mu_{\mathbf{z}_{lk}; \delta_{lk}}(\mathbf{z}_{lk}) = \mu_{f\mathbf{y}_l}; \mathbf{z}_{lk}(\mathbf{z}_{lk})$, and we define the linear transformation matrix

$$ \mathbf{A}(\mathbf{h}_{lk}) = \mathbf{I}_T \otimes [\mathbf{h}_{lk} \quad \mathbf{S}\mathbf{h}_{lk}]. \quad (33) $$

For simplicity, denote

$$ \tilde{\mathbf{H}}_{lk} = [\mathbf{h}_{lk} \quad \mathbf{S}\mathbf{h}_{lk}]. \quad (34) $$

The message $\mu_{\mathbf{z}_{lk}; \delta_{lk}}$ can be explicitly described as

$$ \mu_{\mathbf{z}_{lk}; \delta_{lk}}(\mathbf{z}_{lk}) = \mathcal{N}(\mathbf{z}_{lk} | \mathbf{m}_{\mathbf{z}_{lk}; \delta_{lk}}, \mathbf{C}_{\mathbf{z}_{lk}; \delta_{lk}}), \quad (35) $$

where $\mathbf{C}_{\mathbf{z}_{lk}; \delta_{lk}}$ is a $2TM \times 2TM$ block diagonal matrix

$$ \mathbf{C}_{\mathbf{z}_{lk}; \delta_{lk}} = \text{blkdiag}(\mathbf{C}_{\mathbf{z}_{lk1}; \delta_{lk1}}, \dots, \mathbf{C}_{\mathbf{z}_{lkT}; \delta_{lkT}}, \dots, \mathbf{C}_{\mathbf{z}_{lkT}; \delta_{lkT}}), \quad (36) $$

with blocks of size $2M \times 2M$ and the operation $\text{blkdiag}(\cdot)$ forms a block diagonal matrix the same way as defined in MatLab. The expression in (32) expands as a Gaussian distribution of $\mathbf{x}_k$. We denote the first and second-order moments of the belief of $\mathbf{h}_{lk}$ at factor node $\delta_{lk}$ as

$$ \mathbb{E}_{b_{\delta_{lk}}}[\mathbf{h}_{lk}] = \mathbf{m}_{b_{\delta_{hlk}}}; \quad \mathbb{E}_{b_{\delta_{lk}}}[\mathbf{h}_{lk}\mathbf{h}_{lk}^{\top}] = \mathbf{R}_{b_{\delta_{hlk}}}, \quad (37) $$

The corresponding covariance matrix of $\mathbf{x}_k$ can be computed from (32) as

$$ \mathbf{C}_{\mathbf{z}_{lk}; \mathbf{x}_k} = \text{blkDiag}(\mathbf{C}_{\mathbf{z}_{lk1}; \mathbf{x}_{k1}}, \dots, \mathbf{C}_{\mathbf{z}_{lkT}; \mathbf{x}_{kT}}, \dots, \mathbf{C}_{\mathbf{z}_{lk1}; \mathbf{x}_{k1}}), \quad (38) $$

where

$$ \mathbf{C}_{\delta_{lk}; \mathbf{x}_{kt}} = \begin{bmatrix} \gamma_{\delta_{lk}; \mathbf{x}_{kt}, \Re} & \gamma_{\delta_{lk}; \mathbf{x}_{kt}, \Im} \\ \gamma_{\delta_{lk}; \mathbf{x}_{kt}, \Im} & \gamma_{\delta_{lk}; \mathbf{x}_{kt}, \Re} \end{bmatrix}^{-1}, \quad (39) $$

$$ \begin{aligned} \gamma_{\delta_{lk}; \mathbf{x}_{kt}, \Re} &= \text{tr}[\mathbf{C}_{\mathbf{z}_{lk}; \delta_{lk}}^{-1} \mathbf{R}_{b_{\delta_{hlk}}}] \\ \gamma_{\delta_{lk}; \mathbf{x}_{kt}, \Im} &= \text{tr}[\mathbf{C}_{\mathbf{z}_{lk}; \delta_{lk}}^{-1} \mathbf{S}_{R_{b_{\delta_{hlk}}}}] \end{aligned} \quad (40) $$

$$ \begin{aligned} \gamma_{\delta_{lk}; \mathbf{x}_{kt}, \Re} &= \text{tr}[\mathbf{C}_{\mathbf{z}_{lk}; \delta_{lk}}^{-1} \mathbf{R}_{b_{\delta_{hlk}}}] \\ \gamma_{\delta_{lk}; \mathbf{x}_{kt}, \Im} &= \text{tr}[\mathbf{C}_{\mathbf{z}_{lk}; \delta_{lk}}^{-1} \mathbf{S}_{R_{b_{\delta_{hlk}}}}^{\top}] \end{aligned} \quad (40) $$

Similarly, the mean of the feedback Gaussian message is obtained by the integral

$$ \mathbf{m}_{\delta_{lk}; \mathbf{x}_k} = [\mathbf{m}_{\delta_{lk}; \mathbf{x}_{k1}}^{\top}, \dots, \mathbf{m}_{\delta_{lk}; \mathbf{x}_{kT}}^{\top}]^{\top}, \quad (41) $$

where

$$ \mathbf{m}_{\delta_{lk}; \mathbf{x}_{kt}} = \mathbf{C}_{\delta_{lk}; \mathbf{x}_{kt}} \mathbf{M}_{b_{\delta_{\bar{H}_{lk}}}}^{\top} \mathbf{C}_{\mathbf{z}_{lk}; \delta_{lk}}^{-1} \mathbf{m}_{\mathbf{z}_{lk}; \delta_{lk}}, \quad (42) $$

$$ \text{with } \mathbf{M}_{b_{\delta_{\bar{H}_{lk}}}} = [\mathbf{m}_{b_{\delta_{hlk}}}, \mathbf{S}\mathbf{m}_{b_{\delta_{hlk}}}] $$

### B. Message from $\mathbf{x}_k$ to $\delta_{lk}$ and the Belief $b_{\mathbf{x}_k}$

Due to (14), (26) and (31), the belief of $\mathbf{x}_k$ at $\delta_{lk}$ is

$$ \begin{aligned} b_{\delta_{\omega,lk}}(\mathbf{x}_k) &= b_{\mathbf{x}_k}(\mathbf{x}_k) = \mu_{\mathbf{x}_k; \delta_{lk}}(\mathbf{x}_k) \mu_{\delta_{lk}; \mathbf{x}_k}(\mathbf{x}_k) \\ &= p(\mathbf{x}_k) \prod_l \mu_{\delta_{lk}; \mathbf{x}_k}(\mathbf{x}_k) \propto p(\mathbf{x}_k) \mathcal{N}(\mathbf{x}_k | \mathbf{m}_{\mathbf{y}|\mathbf{x}_k}, \mathbf{C}_{\mathbf{y}|\mathbf{x}_k}) \end{aligned} \quad (43) $$

where

$$ \mathbf{C}_{\mathbf{y}|\mathbf{x}_k} = \left( \sum_l \mathbf{C}_{\delta_{lk}; \mathbf{x}_k}^{-1} \right)^{-1} ; \quad \mathbf{m}_{\mathbf{y}|\mathbf{x}_k} = \mathbf{C}_{\mathbf{y}|\mathbf{x}_k} \left( \sum_l \mathbf{C}_{\delta_{lk}; \mathbf{x}_k}^{-1} \mathbf{m}_{\delta_{lk}; \mathbf{x}_k} \right) \quad (44) $$

Due to the block diagonal structure of the covariance matrix, we denote the t-th block as

$$ \mathbf{C}_{\mathbf{y}|\mathbf{x}_{kt}} = [\mathbf{C}_{\mathbf{y}|\mathbf{x}_{k}}]_{2(t-1)+1:2t, 2(t-1)+1:2t} \quad (45) $$

$$ \mathbf{m}_{\mathbf{y}|\mathbf{x}_{kt}} = [\mathbf{m}_{\mathbf{y}|\mathbf{x}_{k}}]_{2(t-1)+1:2t} $$
The symbol-wise posterior mean and covariance matrix for $x_{kt}$ according to (43) are

$$
\begin{align*}
\mathbf{m}_{b_{\mathbf{x},kt}} &= \frac{1}{Z} \sum_{i} \mathbf{c}_i p(\mathbf{x}_{kt} = \mathbf{s}_i) \mathcal{N}(\mathbf{s}_i | \mathbf{m}_{\mathbf{y}|\mathbf{x}_{kt}}, \mathbf{C}_{\mathbf{y}|\mathbf{x}_{kt}}) \\
\mathbf{C}_{b_{\mathbf{x},kt}} &= \frac{1}{Z} \sum_{i} \mathbf{c}_i \mathbf{s}_i^\top p(\mathbf{x}_{kt} = \mathbf{s}_i) \mathcal{N}(\mathbf{s}_i | \mathbf{m}_{\mathbf{y}|\mathbf{x}_{kt}}, \mathbf{C}_{\mathbf{y}|\mathbf{x}_{kt}}) \\
&\quad - \mathbf{m}_{b_{\mathbf{x},kt}} \mathbf{m}_{b_{\mathbf{x},kt}}^\top,
\end{align*}
$$

where $Z = \sum_i p(x_{kt} = s_i) \mathcal{N}(s_i | m_{y|x_{kt}}, C_y|x_{kt})$.

Combine the symbol-wise means and covariance matrices into
sequence-wise mean and covariance matrix

$$
\begin{align}
\mathbf{m}_{b_{\mathbf{x},k}} &= \left[ \mathbf{m}_{b_{\mathbf{x},k,1}}^{\mathrm{T}}, \dots, \mathbf{m}_{b_{\mathbf{x},k,T}}^{\mathrm{T}} \right]^{\mathrm{T}} \nonumber \\
\mathbf{C}_{b_{\mathbf{x},k}} &= \mathrm{blkdiag}(\mathbf{C}_{b_{\mathbf{x},k,1}}, \dots, \mathbf{C}_{b_{\mathbf{x},k,T}}) \tag{46}
\end{align}
$$

Thus, the posterior means and covariance matrices of the real and
imaginary parts become

$$
\begin{align}
\mathbf{m}_{b_{\mathbf{x},k,\Re}} &= \mathbf{F}_1 \mathbf{m}_{b_{\mathbf{x},k}}; \quad \mathbf{m}_{b_{\mathbf{x},k,\Im}} = \mathbf{F}_2 \mathbf{m}_{b_{\mathbf{x},k}} \nonumber \\
\mathbf{C}_{b_{\mathbf{x},k,\Re}} &= \mathbf{F}_1 \mathbf{C}_{b_{\mathbf{x},k}} \mathbf{F}_1^\mathrm{T}; \quad \mathbf{C}_{b_{\mathbf{x},k,\Im}} = \mathbf{F}_1 \mathbf{C}_{b_{\mathbf{x},k}} \mathbf{F}_2^\mathrm{T} \tag{47} \\
\mathbf{C}_{b_{\mathbf{x},k,\Im}} &= \mathbf{F}_2 \mathbf{C}_{b_{\mathbf{x},k}} \mathbf{F}_1^\mathrm{T}; \quad \mathbf{C}_{b_{\mathbf{x},k,\Im}} = \mathbf{F}_2 \mathbf{C}_{b_{\mathbf{x},k}} \mathbf{F}_2^\mathrm{T} \nonumber
\end{align}
$$

where $\mathbf{F}_1 = \mathbf{I}_T \otimes [1 \ 0]$, $\mathbf{F}_2 = \mathbf{I}_T \otimes [0 \ 1]$. It is easy to see that $\mathbf{C}_{b_{\mathbf{x},k,\Re}}$ is a $T \times T$ diagonal matrix where the $T$-th diagonal entry corresponds to the top-left elements of the $T$-th block matrix in the block diagonal matrix $\mathbf{C}_{b_{\mathbf{x},k}}$. Similarly, $\mathbf{C}_{b_{\mathbf{x},k,\Im}}$, $\mathbf{C}_{b_{\mathbf{x},k,\Re}}$ and $\mathbf{C}_{b_{\mathbf{x},k,\Im}}$ correspond to the top-right, bottom left and bottom right elements of the $T$-th block matrix in $\mathbf{C}_{b_{\mathbf{x},k}}$ respectively.

C. Message from $\delta_{lk}$ to $\mathbf{h}_{lk}$

We consider the feedback message

$$
\mu_{\delta_{lk}; \mathbf{h}_{lk}}(\mathbf{h}_{lk}) = e^{\int b_{\mathbf{x},k}(\mathbf{x}_k) \ln \mu_{\mathbf{z}_{lk}; \delta_{lk}}[\mathbf{B}(\mathbf{x}_k)\mathbf{h}_{lk}]d\mathbf{x}_k}, \quad (48)
$$

where $\mathbf{B}(\mathbf{x}_k) = \mathbf{x}_{\Re,k} \otimes \mathbf{I}_{2M} + \mathbf{x}_{\Im,k} \otimes \mathbf{S}$. The expression (48) result to a Gaussian distribution of $\mathbf{h}_{lk}$. Due to the expectation (integral) operation, we get

$$
\begin{align}
\mathbf{C}_{\delta_{lk}; \mathbf{h}_{lk}} &= \left( \sum_t r_{b_{x,kt},\Re} \mathbf{C}_{\mathbf{z}_{lkt};\delta_{lk}}^{-1} + r_{b_{x,kt},\Im} \mathbf{S}^\mathrm{T} \mathbf{C}_{\mathbf{z}_{lkt};\delta_{lk}}^{-1} \mathbf{S} \right) \tag{49} \\
&\quad \left( r_{b_{x,kt},\Re} \mathbf{C}_{\mathbf{z}_{lkt};\delta_{lk}}^{-1} \mathbf{S} + r_{b_{x,kt},\Im} \mathbf{S}^\mathrm{T} \mathbf{C}_{\mathbf{z}_{lkt};\delta_{lk}}^{-1} \right)^{-1}, \nonumber
\end{align}
$$

and the mean

$$
\begin{equation}
\begin{split}
\mathbf{m}_{\delta_{lk}; \mathbf{h}_{lk}} = \mathbf{C}_{\delta_{lk}; \mathbf{h}_{lk}} & \cdot \sum_t (m_{b_{x,kt},\Re} \mathbf{I}_{2M} + m_{b_{x,kt},\Im} \mathbf{S}^\mathrm{T}) \mathbf{C}_{\mathbf{z}_{lkt};\delta_{lk}}^{-1} \mathbf{m}_{\mathbf{z}_{lkt};\delta_{lk}}.
\end{split}
\tag{50}
\end{equation}
$$

D. Message from $\delta_{lk}$ to $\mathbf{z}_{lk}$

The vectorization of $\mathbf{Z}_{lk}$ can be expressed as

$$
\mathbf{z}_{lk} = \mathbf{x}_{\Re,k} \otimes \mathbf{h}_{lk} + \mathbf{x}_{\Im,k} \otimes (\mathbf{S}\mathbf{h}_{lk}). \quad (51)
$$

The posterior (belief) mean according to the relation (21) becomes

$$
\mathbf{m}_{b_{\delta, \mathbf{z}, \mathbf{l}, \mathbf{k}}} = \mathbf{m}_{b_{\mathbf{x}, k, \Re}} \otimes \mathbf{m}_{b_{\delta, \mathbf{h}, \mathbf{l}, \mathbf{k}}} + \mathbf{m}_{b_{\mathbf{x}, k, \Im}} \otimes (\mathbf{S} \mathbf{m}_{b_{\delta, \mathbf{h}, \mathbf{l}, \mathbf{k}}}) \quad (52)
$$

The correlation matrix of $\mathbf{z}_{lk}$ is:

$$
\begin{align}
\mathbf{R}_{b_{\delta, \mathbf{z}, \mathbf{l}, \mathbf{k}}} ={}& \mathbf{R}_{b_{\mathbf{x}, k, \Re}} \otimes \mathbf{R}_{b_{\delta, \mathbf{h}, \mathbf{l}, \mathbf{k}}} + \mathbf{R}_{b_{\mathbf{x}, k, \Re}} \otimes (\mathbf{R}_{b_{\delta, \mathbf{h}, \mathbf{l}, \mathbf{k}}} \mathbf{S}^{\mathrm{T}}) \tag{53} \\
& + \mathbf{R}_{b_{\mathbf{x}, k, \Im}} \otimes (\mathbf{S} \mathbf{R}_{b_{\delta, \mathbf{h}, \mathbf{l}, \mathbf{k}}}) + \mathbf{R}_{b_{\mathbf{x}, k, \Im}} \otimes (\mathbf{S} \mathbf{R}_{b_{\delta, \mathbf{h}, \mathbf{l}, \mathbf{k}}} \mathbf{S}^{\mathrm{T}}) \nonumber
\end{align}
$$

To make the algorithm trackable, we project $b_{\delta_{z,lk}}$ to Gaussian with
block diagonal covariance matrix (block size $2M$).

$$
\mathbf{C}_{b_{\delta, z, l_k}} = \mathrm{blkdiag}_{2M} [\mathbf{R}_{b_{\delta, z, l_k}} - \mathbf{m}_{b_{\delta, z, l_k}} \mathbf{m}_{b_{\delta, z, l_k}}^{\mathrm{T}}], \quad (54)
$$

where the operation `blkdiag2M`[·] : ℝ²M×2M → ℝ²M×2M keeps only the block matrices of size 2M along the diagonal while setting all the other elements to zero. Another treatment is to project to Gaussian with identical block diagonal matrices:

$$
\mathbf{C}_{b_{\delta, z, l_k}} = \mathbf{I}_T \otimes \frac{1}{T} \mathrm{blktrace}_{2M} [\mathbf{R}_{b_{\delta, z, l_k}} - \mathbf{m}_{b_{\delta, z, l_k}} \mathbf{m}_{b_{\delta, z, l_k}}^{\mathrm{T}}], \quad (55)
$$

where `blktrace2M`[·] : ℝ²M×2M → ℝ²M×2M sums over all the
block matrices of size 2M along the diagonal. Thus, the feedback
message becomes

$$
\mu_{\delta_{lk}; \mathbf{z}_{lk}}(\mathbf{z}_{lk}) = \frac{\mathcal{N}(\mathbf{z}_{lk} | \mathbf{m}_{b_{\delta_{z},lk}}, \mathbf{C}_{b_{\delta_{z},lk}})}{\mathcal{N}(\mathbf{z}_{lk} | \mathbf{m}_{\mathbf{z}_{lk};\delta_{lk}}, \mathbf{C}_{\mathbf{z}_{lk};\delta_{lk}})}. \quad (56)
$$

Note that the above division may lead to a message that cannot be
normalized to a distribution.

Define the eigenvalue matrix $\Lambda_{\delta_{lk};z_{lk}}$ and unitary eigenvector $\mathbf{U}_{\delta_{lk};z_{lk}}$ such that $\mathbf{C}_{\delta_{lk};z_{lk}} = \mathbf{U}_{\delta_{lk};z_{lk}} \Lambda_{\delta_{lk};z_{lk}} \mathbf{U}_{\delta_{lk};z_{lk}}^\mathrm{T}$. We propose the following correction: for all $\lambda \in \mathrm{diag}[\Lambda_{\delta_{lk};z_{lk}}]$, we clip $\lambda^{-1}$ to the range $[10^{-8}, 10^8]$. Since we are using the iterative algorithm to find the fixed point of the BFE, resetting the value will not change the final result.

E. Message from $f_{\mathbf{h}_{lG_g}}$ to $\mathbf{h}_{lk}$ and Belief at $f_{\mathbf{h}_{lG_g}}$

The posterior belief at $f_{\mathbf{h}_{lG_g}}$ is

$$
\begin{equation}
\begin{split}
b_{f_{\mathbf{h}_{lG_g}}}(\mathbf{h}_{lG_g}) ={}& \mathcal{N}(\tilde{\mathbf{y}}_{p,lg} | P\sigma_x^2 \sum_{k' \in G_g} \mathbf{h}_{lk'}, \mathbf{C}_{v,p}) \\
& \qquad \prod_{k' \in G_g} \mu_{\mathbf{h}_{lk'}}:f_{\mathbf{h}_{lG_g}}(\mathbf{h}_{lk'}) p(\mathbf{h}_{lk'}),
\end{split}
\tag{57}
\end{equation}
$$

where $\mu_{\mathbf{h}_{lk}}:f_{\mathbf{h}_{lG_g}}(\mathbf{h}_{lk}) = \mu_{\delta_{lk}}:\mathbf{h}_{lk}(\mathbf{h}_{lk})$. The feedback message is

$$
\begin{align}
\mu_{f_{\mathbf{h}_{lG_g}}; \mathbf{h}_{lk}}(\mathbf{h}_{lk}) &= \int b_{f_{\mathbf{h}_{lG_g}}} (\mathbf{h}_{lG_g}) d\mathbf{h}_{l\bar{k}} / \mu_{f_{\mathbf{h}_{lG_g}}; \mathbf{h}_{lk}}(\mathbf{h}_{lk}) \nonumber \\
&= \int \mathcal{N}(\mathbf{h}_{lk} | \frac{\tilde{\mathbf{y}}_{p,lg}}{P\sigma_x^2} - \sum_{k' \in G_g / \{k\}} \mathbf{h}_{lk'}, \frac{\mathbf{C}_{v,p}}{(P\sigma_x^2)^2}) \nonumber \\
&\quad \cdot \prod_{k' \in G_g / \{k\}} \mu_{\mathbf{h}_{lk'}}:f_{\mathbf{h}_{lG_g}}(\mathbf{h}_{lk'}) p(\mathbf{h}_{lk'}) d\mathbf{h}_{l\bar{k}} \tag{58}
\end{align}
\[
\text{Gaussian interference distribution}
\]
$$

We denote the interference plus noise mean and covariance as

$$
\begin{align}
\mathbf{C}_{\mathbf{h}_{\bar{k}}} &= \frac{\mathbf{C}_{v,p}}{(P\sigma_x^2)^2} + \sum_{k' \in G_g / \{k\}} \mathbf{C}_{\mathbf{h}_{lk'} | y} \tag{59} \\
\mathbf{m}_{\mathbf{h}_{\bar{k}}} &= \sum_{k' \in G_g / \{k\}} \mathbf{m}_{\mathbf{h}_{lk'} | y}, \nonumber
\end{align}
$$

where

$$
\begin{align}
\mathbf{C}_{\mathbf{h}_{lk}|y} &= (\Xi_{\mathbf{h}_{lk}}^{-1} + \mathbf{C}_{\mathbf{h}_{lk};f_{\mathbf{h}_{lG_g}}}^{-1})^{-1} \nonumber \\
\mathbf{m}_{\mathbf{h}_{lk}|y} &= \mathbf{C}_{\mathbf{h}_{lk}|y}(\mathbf{C}_{\mathbf{h}_{lk};f_{\mathbf{h}_{lG_g}}}^{-1} \mathbf{m}_{\mathbf{h}_{lk};f_{\mathbf{h}_{lG_g}}}). \tag{60}
\end{align}
$$

Due to the Gaussian reproduction lemma [ ], the feedback message
can be obtained as

$$
\begin{align}
\mathbf{C}_{f_{\mathbf{h}_{lG_g}}; \mathbf{h}_{lk}} &= (\mathbf{C}_{\mathbf{h}_{l\bar{k}}}^{-1} + \Xi_{\mathbf{h}_{lk}}^{-1})^{-1} \nonumber \\
\mathbf{m}_{f_{\mathbf{h}_{lG_g}}; \mathbf{h}_{lk}} &= \mathbf{C}_{f_{\mathbf{h}_{lG_g}}; \mathbf{h}_{lk}} \left[ \mathbf{C}_{\mathbf{h}_{l\bar{k}}}^{-1} \left( \frac{\tilde{\mathbf{y}}_{p,lg}}{P\sigma_x^2} - \mathbf{m}_{\mathbf{h}_{l\bar{k}}} \right) \right]. \tag{61}
\end{align}
$$
Algorithm 1 Proposed Method

Require: $\forall l, g, k, \tilde{y}_{p,tg}, y_t, p(x_k), p(h_{lk}), p(y_t|z_{l1}, \dots, z_{lK})$

1: Initialize: All the factor-to-variable messages are set to zero mean and unit covariance matrices.

2: $[\forall l, \text{At AP } l, \text{ execute the following loop}]$

3: repeat $[\forall l' \in N(l)k, g]$

4: $C_{f_{h_{lG_g}};h_{lk}}$ and $m_{f_{h_{lG_g}};h_{lk}}$ by (60) $\rightarrow$ (59) $\rightarrow$ (61)

5: $C_{f_{y_l};z_{lk}}$ and $m_{f_{y_l};z_{lk}}$ by (65)

6: $C_{b_{\delta_{h,lk}}}$ and $m_{b_{\delta_{h,lk}}}$ by (62)

7: $R_{b_{\delta_{h,lk}}} = C_{b_{\delta_{h,lk}}} + m_{b_{\delta_{h,lk}}} m_{b_{\delta_{h,lk}}}^T$

8: $C_{\delta_{l_k};x_k}$ and $m_{\delta_{l_k};x_k}$ by (40) $\rightarrow$ (39) $\rightarrow$ (42) $\rightarrow$ (41) $\rightarrow$ (38)

9: Compute $m_{b_{\Phi_k,k},\mathbb{R}}$, $m_{b_{\Phi_k,k},\mathbb{R}}$, $C_{b_{\Phi_k,k},\mathbb{R}}$, $C_{b_{\Phi_k,k},\mathbb{R}}$, $C_{b_{\Phi_k,k},\mathbb{R}}$, $C_{b_{\Phi_k,k},\mathbb{R}}$ via (45)-(47)

10: $C_{\delta_{l_k};h_{lk}}$ and $m_{\delta_{l_k};h_{lk}}$ by (49) and (50)

11: $C_{\delta_{l_k};z_{lk}}$ and $m_{\delta_{l_k};z_{lk}}$ via (52) $\rightarrow$ (53) $\rightarrow$ (54) or (55) $\rightarrow$ (56)

12: until Convergence

According to relation (58), the mean and covariance matrix of the belief are

$$ 
\begin{aligned}
C_{b_{\delta_{h,lk}}} &= \left( C_{f_{h_{lG_g}};h_{lk}}^{-1} + C_{\delta_{l_k};h_{lk}}^{-1} \right)^{-1} \\
m_{b_{\delta_{h,lk}}} &= C_{b_{\delta_{h,lk}}} \left( C_{f_{h_{lG_g}};h_{lk}}^{-1} m_{f_{h_{lG_g}};h_{lk}} + C_{\delta_{l_k};h_{lk}}^{-1} m_{\delta_{l_k};h_{lk}} \right)
\end{aligned}
\quad (62) $$

F. Message from $f_{y_l}$ to $z_{lk}$

The belief at $f_{y_l}$ is

$$ b_{f_{y_l}}(z_{l\{k\}}) = \mathcal{N}(y_l | \sum_k z_{lk}, C_v) \prod_k \mu_{z_{lk};f_{y_l}}(z_{lk}), \quad (63) $$

where $\mu_{z_{lk};f_{y_l}}(z_{lk}) = \mu_{\delta_{l_k};z_{lk}}(z_{lk})$. The feedback message is then obtained by

$$ \mu_{f_{y_l};z_{lk}}(z_{lk}) = \frac{\int b_{f_{y_l}}(z_{l\{k\}}) d\bar{z}_{l\bar{k}}}{\mu_{z_{lk};f_{y_l}}(z_{lk})} \quad (64) $$

which can be computed to be a Gaussian with mean and covariance matrix:

$$ m_{f_{y_l};z_{lk}} = y_l - \sum_{k' \neq k} m_{z_{lk};f_{y_l}}; \quad C_{f_{y_l};z_{lk}} = C_v + \sum_{k' \neq k} C_{z_{lk};f_{y_l}}. \quad (65) $$

For the algorithm with covariance averaging, the overall complexity at every AP is $O(M^3K+M^2KT)$. While the CPU has a complexity of $O(LKT+KT|S|)$. For algorithms using (54) instead of (55), the complexity at the APs will be $O(M^3KT)$.

V. SIMULATION RESULTS

In this section, we verify the algorithm using numerical simulations. We consider a $400m \times 400m$ area with $M = 16$ APs and $K = 8$ UTs. The APs are located at the coordinates $(\frac{400}{3}i, \frac{400}{3}j)$, where $i, j \in \{0, \dots, 3\}$. The UTs are uniformly randomly distributed over this area. The fading model we use is [6],

$$ \sigma_{l,k}^2[\text{dB}] = -30.5 - 36.7 \log_{10}(d_{lk}), \quad (66) $$

where $d_{lk}$ is the distance between AP $l$ and UT $k$. To induce pilot contamination, the pilot sequence length is set to $P = 4$. Furthermore, the pilots are randomly assigned to the users. We utilize 4-QAM and 16-QAM modulation schemes to generate $X$, with the data length set to $T = 16$. To ensure fairness in the simulations, power control is applied for each UT, maintaining equal total received power across

<figure><img src="image_2.png" /><figcaption>Fig. 1. NMSE vs SNR</figcaption></figure>

all users. Our proposed method is compared to an alternative BFE-based approach derived under circular complex Gaussian constraints. In Genie-Aided scenarios, we assume the user data to be known for performance benchmarking.

Compared to the complex version of EP [12] with a complexity of $O[(M^3 + |S|)KT]$ at each AP, the complexity of our proposed method with covariance averaging is $O(M^3K + M^2KT)$.

VI. CONCLUSIONS

We propose a BFE-constrained optimization algorithm that separates complex quantities into in-phase and quadrature components, enabling it to handle a wider range of messages for approximating the true posterior. Furthermore, the complexity of the algorithm is effectively managed by employing averaged covariance matrices.

**Acknowledgements** EURECOM's research is partially supported by its industrial members: ORANGE, BMW, SAP, iABG, Norton LifeLock, and by the Franco-German project CellFree6G.

REFERENCES

[1] R. Gholami, L. Cottatellucci, and D. Slock, “Tackling Pilot Contamination in Cell-Free Massive MIMO by Joint Channel Estimation and Linear Multi-User Detection,” in IEEE International Symposium on Information Theory (ISIT), 2021.

[2] T. P. Minka, “A Family of Algorithms for Approximate Bayesian Inference,” Ph.D. dissertation, Massachusetts Institute of Technology, 2001.

[3] 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, pp. 1–305, 2008.

[4] R. Gholami, L. Cottatellucci, and D. Slock, “Message Passing for a Bayesian Semi-Blind Approach to Cell-Free Massive MIMO,” in Asilomar Conference on Signals, Systems, and Computers. IEEE, 2021.

[5] Z. Zhao and D. Slock, “Bilinear Hybrid Expectation Maximization and Expectation Propagation for Semi-Blind Channel Estimation,” in European Signal Processing Conference (EUSIPCO), 2024.

[6] A. Karataev, C. Forsch, and L. Cottatellucci, “Bilinear Expectation Propagation for Distributed Semi-Blind Joint Channel Estimation and Data Detection in Cell-Free Massive MIMO,” IEEE Open Journal of Signal Processing, 2024.

[7] J. T. Parker, P. Schniter, and V. Cevher, “Bilinear Generalized Approximate Message Passing—Part I: Derivation,” IEEE Transactions on Signal Processing, 2014.

[8] Z. Zhao and D. Slock, “Decentralized Expectation Propagation for Semi-Blind Channel Estimation in Cell-Free Networks,” 2024. [Online]. Available: https://arxiv.org/abs/2410.01303

[9] 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, 2021.

[10] H. Jiang, X. Yuan, and Q. Guo, “Hybrid Vector Message Passing for Generalized Bilinear Factorization,” arXiv preprint arXiv:2401.03626, 2024.

[11] 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.

[12] Z. Zhao and D. Slock, “Decentralized Message-Passing for Semi-Blind Channel Estimation in Cell-Free Systems Based on Bethe free energy Optimization,” in Asilomar Conference on Signals, Systems, and Computers, 2024.