# Decentralized Expectation Propagation for Semi-Blind Channel Estimation in Cell-Free Networks

Zilu Zhao, Dirk Slock

► To cite this version:

Zilu Zhao, Dirk Slock. Decentralized Expectation Propagation for Semi-Blind Channel Estimation in Cell-Free Networks. ISIT 2024, IEEE International Symposium on Information Theory, IEEE, Jul 2024, Athens, Greece. pp.2152-2156, 10.1109/ISIT57864.2024.10619311. hal-05301307

HAL Id: hal-05301307

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

Submitted on 7 Oct 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>
# Decentralized Expectation Propagation for Semi-Blind Channel Estimation in Cell-Free Networks

Zilu Zhao, Dirk Slock

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

**Abstract**—This paper serves as a correction to the conference version [1]. In this work, we explore uplink communication in cell-free (CF) massive multiple-input multiple-output (MaMIMO) systems, employing semi-blind transmission structures to mitigate pilot contamination. We propose a simplified, decentralized method based on Expectation Propagation (EP) for semi-blind channel estimation. By utilizing orthogonal pilots, we preprocess the received signals to establish a simplified equivalent factorization scheme for the transmission process. Moreover, this study integrates Central Limit Theory (CLT) with EP, eliminating the need to introduce new auxiliary variables in the factorization scheme. We also refine the algorithm by assessing the variable scales involved. Finally, a decentralized approach is proposed to significantly reduce the computational demands on the Central Processing Unit (CPU).

## I. INTRODUCTION

One of the unique features of Cell-Free (CF) Massive MIMO (MaMIMO) networks is user terminals (UTs) in a given area are served by all the access points (APs) in the same area. This leads to the problem of pilot contamination. Addressing this issue, Semi-Blind approaches [2] have been explored to mitigate the effects of pilot contamination. In [2], the authors consider a deterministic approach treating the channel coefficients as deterministic unknown parameters and delve into the analysis of the Cramer-Rao bound and identifiability. In the context of Bayesian inference, the Semi-Blind approach is modeled as a bilinear inference problem, where the APs must jointly estimate both the channel coefficients and the user signals. Message-passing algorithms play a critical role here. (Loopy) Belief Propagation (BP) [3] is a widely used method in Bayesian inference. It is an iterative method exploiting the structure for a given factorization scheme (e.g., joint/posterior distribution) to lower the computational loads. One of the most potent message-passing algorithms is Expectation Propagation (EP) [4]. Besides exploiting the factorization scheme, EP also projects the beliefs (marginal posterior) to a family of simple distributions. We can consider BP a special case of EP, in which we assume the projection destination set to be the set of all distributions.

A centralized iterative algorithm has been explored based on variable level EP (VL-EP) [5] in which the authors assume Gaussian input and combine VL-EP with Expectation Maximization (EM). A more recent development is the distributed method proposed in [6]. This approach distributes the computational load at the Central Processing Unit (CPU) by enabling each Access Point (AP) to carry out part of the computation. In [7], the author introduces Decentralized Generalized Approximate Message Passing (D-GAMP). This method is a hybrid of Consensus Propagation [8] and Approximate Message Passing (AMP), effectively eliminating the need for CPU.

### A. Main Contributions

This paper presents a simplified, decentralized, EP-based method designed to address the Semi-Blind estimation problem in communication systems. By utilizing orthogonal pilots, we are able to decouple the channels for different users into mutually exclusive groups, which reduces computational complexity. To further decrease computational demands, we integrate Expectation Propagation (EP) with Central Limit Theory (CLT), treating the interference as noise. Drawing inspiration from [9], we introduce further simplifications through scale analysis. Additionally, to lessen the load on the central processing unit (CPU), we explore a decentralized scheme.

## II. SYSTEM MODEL

We consider a semi-blind signal model containing L APs. At the l-th AP,

$$[\mathbf{Y}_{p,l} \quad \mathbf{Y}_l] = \mathbf{H}_l [\mathbf{X}_p \quad \mathbf{X}] + [\mathbf{V}_{p,l} \quad \mathbf{V}_l]. \quad (1)$$

The received signals are composed of pilot part $\mathbf{Y}_{p,l} \in \mathbb{C}^{N \times P}$ and data part $\mathbf{Y}_l \in \mathbb{C}^{N \times T}$. The channels between different users are considered independent Gaussian i.e. $\text{vec}(\mathbf{H}_l) \sim \mathcal{CN}(0, \Xi_l)$ where $\Xi_l \in \mathbb{C}^{N \times K \times N K}$ is a block diagonal matrix of K blocks $\Xi_{h_{lk}} \in \mathbb{C}^{N \times N}$. The transmitted symbols can be decomposed as pilot symbols $\mathbf{X}_p \in \mathcal{S}^{K \times P}$ and data symbols $\mathbf{X} \in \mathcal{S}^{K \times T}$, where $\mathcal{S}$ is the constellation set. We assume that the elements $x_{kt}$ in $\mathbf{X}$ follow the categorical distribution $p(x_{kt})$. The signal power is denoted as $\sigma_x^2$. The noise is considered as i.i.d. Gaussian distribution, and thus, $\text{vec}([\mathbf{V}_{p,l} \quad \mathbf{V}_l]) \sim \mathcal{CN}(0, \sigma_v^2 \mathbf{I})$.

### A. Orthogonal Pilot sequences

If orthogonal pilot sequences are used, we can first preprocess the pilot observation by right multiplying it with $\mathbf{x}_{p,g}^*$ which is the conjugated g-th pilot sequence. This results in an equivalent observation $\mathbf{y}_{p,lg}$

$$\mathbf{y}_{p,lg} = \mathbf{Y}_{p,l} \mathbf{x}_{p,g}^* = \sum_{k \in G_g} P \sigma_x^2 \mathbf{h}_{lk} + \mathbf{v}_{p,lg} \quad (2)$$

where $\mathbf{v}_{p,lg} = \mathbf{V}_{p,l} \mathbf{x}_{p,g}^* \sim \mathcal{N}(\mathbf{v}_{p,lg} | 0, P \sigma_x^2 \sigma_v^2 \mathbf{I})$, $G_g$ denote the set of users using the g-th pilot sequence. We observe that every $\mathbf{h}_{lk}$ occurs only in one group $G_g$, and the cross-correlation $\mathbb{E}[\mathbf{v}_{p,lg} \mathbf{v}_{p,lg'}^H]$ is an all-zero matrix for all $g \neq g'$. Therefore, the observations $\mathbf{y}_{p,lg}$ and $\mathbf{y}_{p,lg'}$ are independent. With orthogonal pilots, the factorization scheme is derived as

$$\begin{aligned} & p(\{\mathbf{y}_{p,lg}\}, \{\mathbf{Y}_l\}, \{\mathbf{H}_l\}, \mathbf{X}, \{\mathbf{V}_l\}) \\ &= \prod_{k,t} p(x_{kt}) \prod_{l} \prod_{t_1=1}^{T} p(\mathbf{y}_{lt_1} | \mathbf{H}_l, \mathbf{x}_{:t_1}) \prod_{g} p(\mathbf{y}_{p,lg}, \mathbf{H}_{lg}) \end{aligned} \quad (3)$$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Partial factor graph</figcaption></figure>

where $H_{lg}$ is a matrix collecting all $k \in G_g$, $h_{lk}$ as its column vectors, and $x_{:t_1}$ denotes the $t_1$-th column of $X$. We will base our EP (BP) algorithm based on this factorization scheme.

### III. EXPECTATION PROPAGATION OVERVIEW

EP approximates the factors in a factorization scheme to simpler ones [10]. With a given factorization, the update algorithm in EP can be interpreted as message passing of two types of messages, i.e., the message $\mu_{\Psi;\theta_i}(\theta_i)$ from factor node $\Psi$ to variable node $\theta_i$ and the message $\mu_{\theta_i;\Psi}(\theta_i)$ from variable $\theta_i$ to factor $\Psi$: [11]

$$ \mu_{\theta_i;\Psi}(\theta_i) \propto \prod_{\Phi \neq \Psi} \mu_{\Phi;\theta_i}(\theta_i); \mu_{\Psi;\theta_i}(\theta_i) \propto \frac{\mathrm{proj}(b_{\Psi}(\theta_i))}{\mu_{\theta_i;\Psi}(\theta_i)}, \quad (4) $$

where $b_{\Psi}(\theta_i)$ is the belief of $\theta_i$ at node $\Psi$:

$$ b_{\Psi}(\theta_i) \propto \mu_{\theta_i;\Psi}(\theta_i) \int \Psi(\theta) \prod_{j \neq i} \mu_{\theta_j;\Psi}(\theta_j) d\theta_{\bar{i}}. \quad (5) $$

The notation $\theta_{\bar{i}}$ denotes all elements in $\theta$ except the i-th one. The operation $\mathrm{proj}(p)$ projects a given distribution $p$ into a target family $Q$ [11], i.e.,

$$ \mathrm{proj}(p) = \arg \min_{q \in Q} \mathrm{KLD}(p \| q), \quad (6) $$

where $\mathrm{KLD}(p \| q) = \int p(\theta) \ln \frac{p(\theta)}{q(\theta)} d\theta$ is the Kullback-Leibler divergence.

BP, on the other hand, can be considered a special form of EP. The only difference between BP and EP is that there is no projection step in BP. We assume all the messages in this paper are normalized to 1.

#### A. Expectation Propagation on Semi-Blind structure

For simplicity, we denote the factors in the factorization scheme (3) as

$$ \Psi_{1,kt} = p(x_{kt}); \Psi_{2,lt} = p(y_{lt}|H_l, x_{:t}); \Psi_{3,lg} = p(y_p,lg, H_{lg}). $$

The factor graph for (3) is illustrated in Fig. 1.

### IV. MESSAGE PASSING DERIVATIONS

This paper uses EP to estimate channel coefficients $h_{lk}$ and BP to estimate the data symbols $x_{kt}$. Furthermore, we specify the projection family of EP in this paper to be Gaussian distributions with diagonal covariance matrices. Now, we will examine each factor and derive its outbound message.

The message from $\Psi_{1,kt}$ to $x_{kt}$ can be computed directly since no projection is needed, i.e., $\mu_{\Psi_{1,kt};x_{kt}}(x_{kt}) = p(x_{kt})$.

#### A. Message from $\Psi_{2,lt}$ to $x_{kt}$

Following (4)-(5), the extrinsic at node $\Psi_{2,lt}$ is updated by

$$ \begin{aligned} \mu_{x_{kt};\Psi_{2,lt}}(x_{kt}) &\propto p(x_{kt}) \prod_{l' \neq l} \mu_{\Psi_{2,l't}';x_{kt}}(x_{kt}) \\ \mu_{h_{lk};\Psi_{2,lt}}(h_{lk}) &\propto \mu_{\Psi_{3,lg};h_{lk}}(h_{lk}) \prod_{t' \neq t} \mu_{\Psi_{2,lt}';h_{lk}}(h_{lk}), \end{aligned} \quad (7) $$

where the extrinsic of $h_{lk}$ can be computed as a Gaussian $\mu_{h_{lk};\Psi_{2,lt}}(h_{lk}) = CN(h_{lk}|m_{h_{lk};\Psi_{2,lt}}, C_{h_{lk};\Psi_{2,lt}})$ with

$$ \begin{aligned} C_{h_{lk};\Psi_{2,lt}} &= \left( C_{\Psi_{3,lg};h_{lk}}^{-1} + \sum_{t' \neq t} C_{\Psi_{2,lt}';h_{lk}}^{-1} \right)^{-1} \\ m_{h_{lk};\Psi_{2,lt}} &= C_{h_{lk};\Psi_{2,lt}} \left( C_{\Psi_{3,lg};h_{lk}}^{-1} m_{\Psi_{3,lg};h_{lk}} + \sum_{t' \neq t} C_{\Psi_{2,lt}';h_{lk}}^{-1} m_{\Psi_{2,lt}';h_{lk}} \right) \end{aligned} $$

According to the EP rule, the message from $\Psi_{2,lt}$ to $x_{kt}$ is

$$ \mu_{\Psi_{2,lt};x_{kt}}(x_{kt}) \propto \frac{\mathrm{proj}[b_{\Psi_{2,lt};x_{kt}}(x_{kt})]}{\mu_{x_{kt};\Psi_{2,lt}}(x_{kt})}, \quad (8) $$

where the belief (approximated posterior) at $\Psi_{2,lt}$ is defined as $b_{\Psi_{2,lt};x_{kt}}(x_{kt})$ with

$$ \begin{aligned} b_{\Psi_{2,lt};x_{kt}}(x_{kt}) &\propto \mu_{x_{kt};\Psi_{2,lt}}(x_{kt}) \sum_{\mathbf{x}_{\bar{k}t}} \int p(y_{lt}|x_{kt}h_{lk} + \sum_{i \neq k} x_{it}h_{li}) \\ &\qquad \cdot \mu_{h_{lk};\Psi_{2,lt}}(h_{lk}) \prod_{i \neq k} \mu_{h_{li};\Psi_{2,lt}}(h_{li}) \mu_{x_{it};\Psi_{2,lt}}(x_{it}) dH_{l}. \end{aligned} \quad (9) $$

We use the notation $\mathbf{x}_{\bar{k}t}$ to denote all the elements in $\mathbf{x}_{:t}$ except the k-th element. The integral (and summation) in (9) can be considered as a marginalization operation. Furthermore, we can view the extrinsic messages as hypothetical priors. Due to CLT, we approximate $\sum_{i \neq k} x_{it}h_{li}$ to a Gaussian where $x_{it} \sim \mu_{x_{it};\Psi_{2,lt}}(x_{it}), h_{li} \sim \mu_{h_{li};\Psi_{2,lt}}(h_{li})$ [11]. Therefore, (9) becomes

$$ b_{\Psi_{2,lt};x_{kt}}(x_{kt}) \propto \mu_{x_{kt};\Psi_{2,lt}}(x_{kt}) \quad (10) $$

$$ \iint p(y_{lt}|x_{kt}h_{lk} + z_{lkt}) \mu_{z_{lkt}}(z_{lkt}) dz_{lkt} \cdot \mu_{h_{lk};\Psi_{2,lt}}(h_{lk}) d\mathbf{h}_{lk}, $$

where $\mu_{z_{lkt}}(z_{lkt}) = CN(z_{lkt}|m_{z_{lkt}}, C_{z_{lkt}})$ with

$$ \mathbf{m}_{z_{lkt}} = \sum_{i \neq k} m_{x_{it};\Psi_{2,lt}} \mathbf{m}_{h_{li};\Psi_{2,lt}} \quad (11) $$

$$ \mathbf{C}_{z_{lkt}} = \sum_{i \neq k} r_{x_{it};\Psi_{2,lt}} \mathbf{C}_{h_{li};\Psi_{2,lt}} + \tau_{x_{it};\Psi_{2,lt}} \mathbf{m}_{h_{it};\Psi_{2,lt}} \mathbf{m}_{h_{it};\Psi_{2,lt}}^{H} $$

where $m_{x_{it};\Psi_{2,lt}}, \tau_{x_{it};\Psi_{2,lt}}$ and $r_{x_{it};\Psi_{2,lt}}$ are the mean, variance and second-order moment of the normalized message $\mu_{x_{it};\Psi_{2,lt}}$. By applying the Gaussian reproduction lemma [11] and the fact that Gaussian distribution integrates to one, the belief (10) becomes

$$ \begin{aligned} b_{\Psi_{2,lt};x_{kt}}(x_{kt}) &\propto CN(0|y_{lt}-\mathbf{m}_{z_{lkt}}-x_{kt}\mathbf{m}_{h_{lk}};\Psi_{2,lt}, \\ &\qquad \mathbf{C}_v+\mathbf{C}_{z_{lkt}}+|x_{kt}|^2\mathbf{C}_{h_{lk}};\Psi_{2,lt}) \cdot \mu_{x_{kt};\Psi_{2,lt}}(x_{kt}). \end{aligned} \quad (12) $$

Therefore, by BP rules, the outbound message is

$$ \begin{aligned} \mu_{\Psi_{2,lt};x_{kt}}(x_{kt}) &\propto CN(0|y_{lt}-\mathbf{m}_{z_{lkt}}-x_{kt}\mathbf{m}_{h_{lk}};\Psi_{2,lt}, \\ &\qquad \mathbf{C}_v+\mathbf{C}_{z_{lkt}}+|x_{kt}|^2\mathbf{C}_{h_{lk}};\Psi_{2,lt}) \end{aligned} \quad (13) $$
### B. Message from $\Psi_{2,lt}$ to $\mathbf{h}_{lk}$

Based on EP rules (4), the message to $\mathbf{h}_{lk}$ is

$$ \mu_{\Psi_{2,lt};\mathbf{h}_{lk}}(\mathbf{h}_{lk}) \propto \frac{\mathrm{proj}[b_{\Psi_{2,lt};\mathbf{h}_{lk}}(\mathbf{h}_{lk})]}{\mu_{\mathbf{h}_{lk};\Psi_{2,lt}}(\mathbf{h}_{lk})}, \quad (14) $$

where the belief is defined as

$$ \begin{aligned} b_{\Psi_{2,lt};\mathbf{h}_{lk}}(\mathbf{h}_{lk}) \propto & \sum_{\mathbf{x};t} \int p(y_{lt} | \sum_{i} x_{it} \mathbf{h}_{li}) \\ & \cdot \prod_{i} \mu_{\mathbf{h}_{li};\Psi_{2,lt}}(\mathbf{h}_{li}) \mu_{x_{it};\Psi_{2,lt}}(x_{it}) d\mathbf{h}_{l\bar{k}} \end{aligned} \quad (15) $$

We use $\mathbf{h}_{l\bar{k}}$ to denote all the column vectors in $\mathcal{H}_l$ except the $k$-th column. By using the same approach from (9) to (12), and separating the terms that contains only $x_{kt}$ [11] [6], the belief (15) becomes

$$ \begin{aligned} & b_{\Psi_{2,lt};\mathbf{h}_{lk}}(\mathbf{h}_{lk}) \\ & = \mathbb{E}_{b_{\Psi_{2,lt};x_{kt}}} \{ \mathcal{CN}[\mathbf{h}_{lk}|\mathbf{m}_{\hat{\mathbf{h}}_{lk}|x_{kt}}(x_{kt}), \mathbf{C}_{\hat{\mathbf{h}}_{lk}|x_{kt}}(x_{kt})] \} \end{aligned} \quad (16) $$

where $\mathbf{m}_{\hat{\mathbf{h}}_{lk}|x_{kt}}(\cdot)$ and $\mathbf{C}_{\hat{\mathbf{h}}_{lk}|x_{kt}}(\cdot)$ are defined as

$$ \begin{aligned} \mathbf{C}_{\hat{\mathbf{h}}_{lk}|x_{kt}}(x) &= [|x|^2(\mathbf{C}_v + \mathbf{C}_{z_{lk}})^{-1} + \mathbf{C}_{\mathbf{h}_{lk};\Psi_{2,lt}}^{-1}]^{-1} \\ \mathbf{m}_{\hat{\mathbf{h}}_{lk}|x_{kt}}(x) &= \mathbf{C}_{\hat{\mathbf{h}}_{lk}|x_{kt}}(x) \left[ \mathbf{C}_{\mathbf{h}_{lk};\Psi_{2,lt}}^{-1} \mathbf{m}_{\mathbf{h}_{lk};\Psi_{2,lt}} \right. \\ & \qquad \left. + |x|^2(\mathbf{C}_v + \mathbf{C}_{z_{lk}})^{-1} \frac{\mathbf{y}_{lt} - \mathbf{m}_{z_{lk}}}{x} \right], \end{aligned} \quad (17) $$

where $\mathbf{C}_v = \sigma_v^2 \mathbf{I}$. The mean $\mathbf{m}_{\hat{\mathbf{h}}_{lk}^2}$ and covariance $\mathbf{C}_{\hat{\mathbf{h}}_{lk}^2}$ of the belief distribution (16) are

$$ \begin{aligned} \mathbf{m}_{\hat{\mathbf{h}}_{lk}^2} &= \mathbb{E}_{b_{\Psi_{2,lt};x_{kt}}} [\mathbf{m}_{\hat{\mathbf{h}}_{lk}|x_{kt}}(x_{kt})] \\ \mathbf{C}_{\hat{\mathbf{h}}_{lk}^2}' &= \mathbb{E}_{b_{\Psi_{2,lt};x_{kt}}} [\mathbf{C}_{\hat{\mathbf{h}}_{lk}|x_{kt}}(x_{kt}) \\ & \quad + \mathbf{m}_{\hat{\mathbf{h}}_{lk}|x_{kt}}(x_{kt}) \mathbf{m}_{\hat{\mathbf{h}}_{lk}|x_{kt}}(x_{kt})^H] - \mathbf{m}_{\hat{\mathbf{h}}_{lk}^2} \mathbf{m}_{\hat{\mathbf{h}}_{lk}^2}^H. \end{aligned} \quad (18) $$

We project the belief at $\Psi_{2,lt}$ to a Gaussian with diagonal covariance matrix $\mathrm{proj}[b_{\Psi_{2,lt};\mathbf{h}_{lk}}(\mathbf{h}_{lk})] = \mathcal{CN}(\mathbf{h}_{lk}|\mathbf{m}_{\hat{\mathbf{h}}_{lk}^2}, \mathbf{C}_{\hat{\mathbf{h}}_{lk}^2})$, where $\mathbf{C}_{\hat{\mathbf{h}}_{lk}^2}$ is a diagonal matrix with the same diagonal elements as $\mathbf{C}_{\hat{\mathbf{h}}_{lk}^2}$. Finally, the message from $\Psi_{2,lt}$ to $\mathbf{h}_{lk}$ is

$$ \begin{aligned} \mu_{\Psi_{2,lt};\mathbf{h}_{lk}}(\mathbf{h}_{lk}) &= \mathcal{CN}(\mathbf{h}_{lk}|\mathbf{m}_{\Psi_{2,lt};\mathbf{h}_{lk}}, \mathbf{C}_{\Psi_{2,lt};\mathbf{h}_{lk}}) \\ & \propto \frac{\mathcal{CN}(\mathbf{h}_{lk}|\mathbf{m}_{\hat{\mathbf{h}}_{lk}^2}, \mathbf{C}_{\hat{\mathbf{h}}_{lk}^2})}{\mathcal{CN}(\mathbf{h}_{lk}|\mathbf{m}_{\mathbf{h}_{lk};\Psi_{2,lt}}, \mathbf{C}_{\mathbf{h}_{lk};\Psi_{2,lt}})}. \end{aligned} \quad (19) $$

### C. Message form $\Psi_{3,lg}$ to $\mathbf{h}_{lk}$

We assume $k \in G_g$. The extrinsic at $\Psi_{3,lg}$ is updated by

$$ \mu_{\mathbf{h}_{lk};\Psi_{3,lg}}(\mathbf{h}_{lk}) \propto \prod_{t} \mu_{\Psi_{2,lt};\mathbf{h}_{lk}}(\mathbf{h}_{lk}). \quad (20) $$

We denote this extrinsic message as $\mu_{\mathbf{h}_{lk};\Psi_{3,lg}}(\mathbf{h}_{lk}) = \mathcal{CN}(\mathbf{h}_{lk}|\mathbf{m}_{\mathbf{h}_{lk};\Psi_{3,lg}}, \mathbf{C}_{\mathbf{h}_{lk};\Psi_{3,lg}})$ with

$$ \begin{aligned} \mathbf{C}_{\mathbf{h}_{lk};\Psi_{3,lg}} &= \left( \sum_t \mathbf{C}_{\Psi_{2,lt};\mathbf{h}_{lk}}^{-1} \right)^{-1} \\ \mathbf{m}_{\mathbf{h}_{lk};\Psi_{3,lg}} &= \mathbf{C}_{\mathbf{h}_{lk};\Psi_{3,lg}} \left( \sum_t \mathbf{C}_{\Psi_{2,lt};\mathbf{h}_{lk}}^{-1} \mathbf{m}_{\Psi_{2,lt};\mathbf{h}_{lk}} \right). \end{aligned} $$

The belief of $\mathbf{h}_{lg}$ at the $\Psi_{3,lg}$ is

$$ b_{\Psi_{3,lg}}(\mathbf{h}_{lg}) \propto p(y_{p,lg}, \mathbf{h}_{lg}) \prod_{k \in G_g} p(\mathbf{h}_{lk}) \mu_{\mathbf{h}_{lk};\Psi_{3,lg}}(\mathbf{h}_{lk}). \quad (21) $$

All the factors appearing in (21) are Gaussian pdfs with diagonal covariance matrices. Therefore, the projection of $b_{\Psi_{3,lg}}(\mathbf{h}_{lk})$ results to itself. For simplicity, we define a hypothetical prior $q_{\mathbf{h}_{lk}|\mathbf{Y}_d}$ for $\mathbf{h}_{lk}$ in (21) as

$$ q_{\mathbf{h}_{lk}|\mathbf{Y}_d}(\mathbf{h}_{lk}) = \mathcal{N}(\mathbf{h}_{lk}|\mathbf{m}_{\mathbf{h}_{lk}|\mathbf{Y}_d}, \mathbf{C}_{\mathbf{h}_{lk}|\mathbf{Y}_d}) \\ \propto p(\mathbf{h}_{lk}) \mu_{\mathbf{h}_{lk};\Psi_{3,lg}}(\mathbf{h}_{lk}), \quad (22) $$

where

$$ \begin{aligned} \mathbf{C}_{\mathbf{h}_{lk}|\mathbf{Y}_d} &= (\Xi_{\mathbf{h}_{lk}}^{-1} + \mathbf{C}_{\mathbf{h}_{lk};\Psi_{3,lg}}^{-1})^{-1} \\ \mathbf{m}_{\mathbf{h}_{lk}|\mathbf{Y}_d} &= \mathbf{C}_{\mathbf{h}_{lk}|\mathbf{Y}_d} \mathbf{C}_{\mathbf{h}_{lk};\Psi_{3,lg}}^{-1} \mathbf{m}_{\mathbf{h}_{lk};\Psi_{3,lg}} \end{aligned} \quad (23) $$

The message from factor node $\Psi_{3,lg}$ to $\mathbf{h}_{lk}$ can be derived as

$$ \begin{aligned} \mu_{\Psi_{3,lg};\mathbf{h}_{lk}}(\mathbf{h}_{lk}) &\propto \frac{\int b_{\Psi_{3,lg}}(\mathbf{h}_{lg}) d\mathbf{h}_{l\bar{k}}}{\mu_{\mathbf{h}_{lk};\Psi_{3,lg}}(\mathbf{h}_{lk})} \\ &\propto p(\mathbf{h}_{lk}) \frac{\int p(\mathbf{y}_{p,lg}, \mathbf{h}_{lg}) \prod_{k \in G_g} q_{\mathbf{h}_{lk}|\mathbf{Y}_d}(\mathbf{h}_{lk}) d\mathbf{h}_{l\bar{k}}}{q_{\mathbf{h}_{lk}|\mathbf{Y}_d}(\mathbf{h}_{lk})} \end{aligned} \quad (24) $$

The fraction operation in the second line of (24) can be interpreted as component-wise conditionally-unbiased LMMSE estimation [12]. Therefore, the message from $\Psi_{3,lg}$ to $\mathbf{h}_{lk}$ is

$$ \mu_{\Psi_{3,lg};\mathbf{h}_{lk}}(\mathbf{h}_{lk}) = \mathcal{CN}(\mathbf{h}_{lk}|\mathbf{m}_{\Psi_{3,lg};\mathbf{h}_{lk}}, \mathbf{C}_{\Psi_{3,lg};\mathbf{h}_{lk}}), \quad (25) $$

where

$$ \begin{aligned} \mathbf{C}_{\Psi_{3,lg};\mathbf{h}_{lk}} &= \left[ \Xi_{\mathbf{h}_{lk}}^{-1} + \left( \frac{\sigma_v^2}{\sigma_x^2 P} \mathbf{I} + \sum_{k' \in G_g / \{k\}} \mathbf{C}_{\mathbf{h}_{lk}'|\mathbf{Y}_d} \right)^{-1} \right]^{-1} \\ \mathbf{m}_{\Psi_{3,lg};\mathbf{h}_{lk}} &= \Xi_{\mathbf{h}_{lk}} \left( \frac{\sigma_v^2}{\sigma_x^2 P} \mathbf{I} + \sum_{k' \in G_g / \{k\}} \mathbf{C}_{\mathbf{h}_{lk}'|\mathbf{Y}_d} + \Xi_{\mathbf{h}_{lk}} \right)^{-1} \\ & \quad \cdot \left( \frac{1}{\sigma_x^2 P} \mathbf{y}_{p,lg} - \sum_{k' \in G_g / \{k\}} \mathbf{m}_{\mathbf{h}_{lk}'|\mathbf{Y}_d} \right). \end{aligned} \quad (26) $$

## V. ASYMPTOTIC BEHAVIORS IN LARGE SYSTEMS

For scalable systems, $L \to \infty$ while $K = c_1 L$, $P = c_2 L$, $T = c_3 L$, where $c_1, c_2, c_3$ are some positive constants, we assume the channel coefficients $\forall l, n, k, \mathbb{E}[|\mathbf{h}_{lkn}|^2] = O(\frac{1}{L})$ data constellation symbols $\forall s \in S, s = O(1), 1/s = O(1)$. Furthermore, we assume the noise power scales as $\sigma_v^2 = O(1)$, $\sigma_v^{-2} = O(1)$. For simplicity, we define big-O-notations with matrix parameters to represent the element-wise asymptotic behavior, i.e., for matrices $\mathbf{A}, \mathbf{B}$ of the same size, we have $\mathbf{A} = O(\mathbf{B}) \Leftrightarrow \forall i, j, [\mathbf{A}]_{ij} = O([\mathbf{B}]_{ij})$.

**Assumption 1.** We assume that $\mathbb{E}(\mathbf{m}_{\mathbf{h}_{lk};\Psi_{2,lt}}) = 0$, $\mathbb{E}[m_{x_{it};\Psi_{2,lt}}\mathbf{m}_{\mathbf{h}_{lk};\Psi_{2,lt}}] = 0$, and $\forall i \neq j, \mathbb{E}([\mathbf{m}_{\mathbf{h}_{lk};\Psi_{2,lt}}\mathbf{m}_{\hat{\mathbf{h}}_{lk};\Psi_{2,lt}}^H]_{ij}) = 0$.

**Property 1.** For invertible matrices $\mathbf{A}, \mathbf{B}$, we have $(\mathbf{A}^{-1} + \mathbf{B}^{-1})^{-1} = \mathbf{A}(\mathbf{A} + \mathbf{B})^{-1}\mathbf{B} = \mathbf{B}(\mathbf{A} + \mathbf{B})^{-1}\mathbf{A}$.

**Lemma 1.** With proper initialization and Assumption 1, in each iteration, the updates satisfy $\mathbf{m}_{\mathbf{z}_{lkt}} = O(\mathbf{1})$, $\mathbf{C}_{\mathbf{z}_{lkt}} = O(\mathbf{I}) + O(\frac{1}{\sqrt{L}})$, $\mathbf{m}_{\hat{\mathbf{h}}_{lk}|x_{kt}}(x) = O(\frac{1}{\sqrt{L}})$, $\mathbf{C}_{\hat{\mathbf{h}}_{lk}|x_{kt}}(x) = O(\frac{1}{L}) + O(\frac{1}{L^2})$, $\mathbf{m}_{\Psi_{2,lt};\mathbf{h}_{lk}} = O(\frac{1}{\sqrt{L}})$, $\mathbf{C}_{\Psi_{2,lt};\mathbf{h}_{lk}} = O(\mathbf{I})$, $\mathbf{m}_{\Psi_{3,lg};\mathbf{h}_{lk}} = O(\frac{1}{\sqrt{L}})$, $\mathbf{C}_{\Psi_{3,lg};\mathbf{h}_{lk}} = O(\frac{1}{L})$, $\mathbf{m}_{\mathbf{h}_{lk};\Psi_{2,lt}} = O(\frac{1}{\sqrt{L}})$, $\mathbf{C}_{\mathbf{h}_{lk};\Psi_{2,lt}} = O(\frac{1}{L})$. Furthermore, $\mathbf{C}_{\Psi_{2,lt};\mathbf{h}_{lk}} = O(\mathbf{I})$.
*Proof.* We prove this lemma by mathematical induction. Due to a proper initialization, we can assume the messages $\mathbf{m}_{\hat{h}_{lk}}; \Psi_{2,lt}$, $\mathbf{C}_{\hat{h}_{lk}}; \Psi_{2,lt}$, $\mathbf{m}_{\Psi_{2,lt}}; \hat{h}_{lk}$, $\mathbf{C}_{\Psi_{2,lt}}; \hat{h}_{lk}$, $\mathbf{m}_{\Psi_{3,lg}}; \hat{h}_{lk}$, $\mathbf{C}_{\Psi_{3,lg}}; \hat{h}_{lk}$ are initialized with the above-mentioned scales. Then, we assume the lemma holds for the previous iterations and investigate the updates in the next iteration.

We first look at the update of $\mathbf{m}_{\mathbf{z}_{lkt}}$, $\mathbf{C}_{\mathbf{z}_{lkt}}$, which are updated according to (11). Similar to [13], we assume that $\forall i, m_{x_{it}}; \Psi_{2,lt}$ are weakly independent of $\mathbf{m}_{h_{li}}; \Psi_{2,lt}$. Since the elements in the constellation set scale with $O(1)$, we know $m_{x_{it}}; \Psi_{2,lt} = O(1)$. According to induction assumptions, the extrinsic mean $\mathbf{m}_{\hat{h}_{li}}; \Psi_{2,lt} = O(\frac{1}{\sqrt{L}})$. Due to Assumption 1, we use the results from [13, Lemma 1] to obtain $\mathbf{m}_{\mathbf{z}_{lkt}} = \sum_{i \neq k} m_{x_{it}}; \Psi_{2,lt} \mathbf{m}_{\hat{h}_{li}}; \Psi_{2,lt} = O(1)$ and the covariance matrix $\mathbf{C}_{\mathbf{z}_{lkt}} = O(\mathbf{I}) + O(\frac{1}{\sqrt{L}})$. For simplicity, we denote the diagonal terms of $\mathbf{C}_{\mathbf{z}_{lkt}}$ in (11) as $\mathbf{D}_{\mathbf{z}_{lkt}} = \sum_{i \neq k} r_{x_{it}}; \Psi_{2,lt} \mathbf{C}_{\hat{h}_{li}}; \Psi_{2,lt}$, and denote $\mathbf{B}_{\mathbf{z}_{lkt}} = \sum_{i \neq k} \tau_{x_{it}}; \Psi_{2,lt} \mathbf{m}_{\hat{h}_{li}}; \Psi_{2,lt} \mathbf{m}_{\hat{h}_{it}}; \Psi_{2,lt} = O(\mathbf{I}) + O(\frac{1}{\sqrt{L}})$. Thus, with these notations, $\mathbf{C}_{\mathbf{z}_{lkt}} = \mathbf{D}_{\mathbf{z}_{lkt}} + \mathbf{B}_{\mathbf{z}_{lkt}}$.

Now we investigate the update of $\mathbf{m}_{\hat{h}_{lk}}|_{x_{kt}}(x)$, $\mathbf{C}_{\hat{h}_{lk}}|_{x_{kt}}(x)$ in (17). By matrix inversion lemma,

$$ \mathbf{C}_{\hat{h}_{lk}|x_{kt}}(x) = \mathbf{C}_{\hat{h}_{lk}}; \Psi_{2,lt} - \mathbf{C}_{\hat{h}_{lk}}; \Psi_{2,lt} \mathbf{Q}_{lkt}^{\frac{H}{2}} \mathbf{W}_{lkt} (\mathbf{\Lambda}_{lkt} + \mathbf{I})^{-1} \mathbf{W}_{lkt}^{H} \mathbf{Q}_{lkt}^{-\frac{1}{2}} \mathbf{C}_{\hat{h}_{lk}}; \Psi_{2,lt}, \quad (27) $$

where we define the positive semi-definite diagonal matrix $\mathbf{Q}_{lkt} = \mathbf{C}_{\hat{h}_{lk}}; \Psi_{2,lt} + \frac{1}{|x|^2} (\mathbf{C}_v + \mathbf{D}_{\mathbf{z}_{lkt}}) = O(\mathbf{I})$ and $\mathbf{Q}_{lkt}^{\frac{1}{2}} \mathbf{Q}_{lkt}^{\frac{H}{2}} = \mathbf{Q}_{lkt}$. By eigendecomposition, we define $\frac{1}{|x|^2} \mathbf{Q}_{lkt}^{-\frac{H}{2}} \mathbf{B}_{\mathbf{z}_{lkt}} \mathbf{Q}_{lkt}^{\frac{H}{2}} = \mathbf{W}_{lkt} \mathbf{\Lambda}_{lkt} \mathbf{W}_{lkt}^{H}$, $\mathbf{I} = \mathbf{W}_{lkt} \mathbf{W}_{lkt}^{H}$. Therefore, $\mathbf{C}_{\hat{h}_{lk}|x_{kt}}(x) = O(\frac{1}{L}) + O(\frac{1}{L^2})$. We find the update of $\mathbf{m}_{\hat{h}_{lk}}|_{x_{kt}}$ in (17) is dominated by the first term $\mathbf{C}_{\hat{h}_{lk}|x_{kt}}(x) \mathbf{C}_{\hat{h}_{lk}}; \Psi_{2,lt} \mathbf{m}_{\hat{h}_{lk}}; \Psi_{2,lt}$. By neglecting higher order infinitesimal terms, we have

$$ \begin{aligned} \mathbf{m}_{\hat{h}_{lk}|x_{kt}} &\approx \mathbf{C}_{\hat{h}_{lk}|x_{kt}}(x) \mathbf{C}_{\hat{h}_{lk}}^{-1} \mathbf{m}_{\hat{h}_{lk}}; \Psi_{2,lt} = [\mathbf{I} + \mathbf{C}_{\hat{h}_{lk}}^{\frac{1}{2}}; \Psi_{2,lt} \cdot (\mathbf{C}_{\hat{h}_{lk}}^{\frac{H}{2}}; \Psi_{2,lt} \mathbf{C}_{\hat{h}_{lk}}^{-1} \mathbf{C}_{\hat{h}_{lk}}^{\frac{1}{2}} + \mathbf{I})^{-1} \mathbf{C}_{\hat{h}_{lk}}^{-\frac{H}{2}} \mathbf{C}_{\hat{h}_{lk}}^{-1} \mathbf{C}_{\hat{h}_{lk}}; \Psi_{2,lt} \\ &\approx \mathbf{m}_{\hat{h}_{lk}}; \Psi_{2,lt} = O(\frac{1}{\sqrt{L}}). \end{aligned} $$

To study the messages $\mathbf{m}_{\Psi_{2,lt}; \hat{h}_{lk}}$, $\mathbf{C}_{\Psi_{2,lt}; \hat{h}_{lk}}$, we first investigate the approximated (projected) belief of $\hat{h}_{lk}$ at $\Psi_{2,lt}$. From the previous discussion, $\mathbf{C}_{\hat{h}_{lk}}' \approx \mathbf{E}_{b_{\Psi_{2,lt};x_{kt}}[\mathbf{C}_{\hat{h}_{lk}|x_{kt}}(x_{kt})]}$, and thus,

$$ \mathbf{C}_{\hat{h}_{lk}^2} \approx \mathbf{C}_{\hat{h}_{lk}}; \Psi_{2,lt} - \mathbf{C}_{\hat{h}_{lk}}; \Psi_{2,lt} \cdot \mathbf{F} \cdot \mathbf{C}_{\hat{h}_{lk}}; \Psi_{2,lt}, \quad (28) $$

where $\mathbf{F} = \sum_{x \in S} b_{\Psi_{2,lt};x_{kt}}(x) [\frac{1}{|x|^2} \mathbf{B}_{\mathbf{z}_{lkt}} + \mathbf{Q}_{lkt}]^{-1}$. Thanks to the projection in EP, we are only interested in the diagonal elements of $\mathbf{C}_{\hat{h}_{lk}^2}'$. With the approximation $\mathbf{m}_{\hat{h}_{lk}|x_{kt}} \approx \mathbf{m}_{\hat{h}_{lk}}; \Psi_{2,lt}$, the n-th diagonal term reads

$$ \begin{aligned} [\mathbf{C}_{\hat{h}_{lk}^2}']_{nn} &= [\mathbf{C}_{\hat{h}_{lk}^2}]_{nn} \approx \tau_{h_{lnk}}; \Psi_{2,lt} - \tau_{h_{lnk}}^2; \Psi_{2,lt} [\mathbf{F}]_{nn} \\ &= [([\mathbf{F}]_{nn}^{-1} - \tau_{h_{lnk}}; \Psi_{2,lt})^{-1} + \tau_{h_{lnk}}^{-1}; \Psi_{2,lt}]^{-1} \end{aligned} \quad (29) $$

From the same analysis in (27), we know $\mathbf{Q}_{lkt}$ is asymptotic upper and lower bounded. Since $\mathbf{B}_{\mathbf{z}_{lkt}^1}$ is positive semi-definite, we have $[\mathbf{F}]_{nn} = O(1)$ and $[\mathbf{F}]_{nn}^1 = O(1)$. Substitute (29) into (19), and we obtain $[\mathbf{C}_{\Psi_{2,lt}; \hat{h}_{lk}}]_{nn} \approx [\mathbf{F}]_{nn}^{-1} - \tau_{h_{lnk}}; \Psi_{2,lt}$. Thus, $\mathbf{C}_{\Psi_{2,lt}; \hat{h}_{lk}} = O(\mathbf{I})$ and $\mathbf{C}_{\Psi_{2,lt}; \hat{h}_{lk}}^{-1} = O(\mathbf{I})$. Since $\mathbf{m}_{\hat{h}_{lk}|x_{kt}} \approx \mathbf{m}_{\hat{h}_{lk}}; \Psi_{2,lt}$, it is straightforward to see $\mathbf{m}_{\Psi_{2,lt}; \hat{h}_{lk}} = O(\frac{1}{\sqrt{L}})$.

The message covariance matrices $\mathbf{C}_{\Psi_{3,lg}; \hat{h}_{lk}}$ and $\mathbf{C}_{\hat{h}_{lk}; \Psi_{2,lt}}$ are both diagonal matrices. Due to the scalable-system assumption $T = O(L)$ and the elements in $\mathbf{C}_{\Psi_{2,lt}; \hat{h}_{lk}}$ being asymptotically upper and lower bounded, one can show $\mathbf{C}_{\Psi_{3,lg}; \hat{h}_{lk}} = O(\frac{1}{L})$ is upper bounded, and $\mathbf{m}_{\Psi_{3,lg}; \hat{h}_{lk}} = O(\frac{1}{\sqrt{L}})$ according to (23)-(26). We can then show $\mathbf{C}_{\hat{h}_{lk}; \Psi_{2,lt}} = O(\frac{1}{L})$ and $\mathbf{m}_{\hat{h}_{lk}; \Psi_{2,lt}} = O(\frac{1}{\sqrt{L}})$ according to (7). □

## VI. SIMPLICATION OF THE MESSAGES

We define beliefs at the variable nodes as

$$ \begin{aligned} b_{x_{kt}}(x_{kt}) &\propto p(x_{kt}) \prod_{l} \mu_{\Psi_{2,lt};x_{kt}}(x_{kt}) \\ b_{\hat{h}_{lk}}(\hat{h}_{lk}) &\propto \mu_{\Psi_{3,lg};\hat{h}_{lk}}(\hat{h}_{lk}) \prod_{t} \mu_{\Psi_{2,lt};\hat{h}_{lk}}(\hat{h}_{lk}). \end{aligned} \quad (30) $$

Compared to the extrinsic messages in (7), the beliefs in (30) only differ by one factor. Since Loopy BP is used for estimating $x_{kt}$, it has been shown in [9] that we can assume $b_{x_{kt}}(x_{kt}) \approx \mu_{x_{kt}}; \Psi_{2,lt}(x_{kt})$.

This work estimates the channel coefficients $\hat{h}_{lk}$ using EP. Therefore, a separate analysis from [9] is needed. We investigate the mean and covariance matrix difference between $b_{\hat{h}_{lk}}(\hat{h}_{lk})$ and $\mu_{\hat{h}_{lk}; \Psi_{2,lt}}$ based on the Lemma 2.

Substitute (7) into (30), and we obtain the belief at $\hat{h}_{lk}$ as

$$ b_{\hat{h}_{lk}}(\hat{h}_{lk}) \propto \mu_{\hat{h}_{lk}}; \Psi_{2,lt}(\hat{h}_{lk}) \mu_{\Psi_{2,lt}; \hat{h}_{lk}}(\hat{h}_{lk}). \quad (31) $$

Denote $\mathbf{m}_{\hat{h}_{lk}}$ and $\mathbf{C}_{\hat{h}_{lk}}$ as the mean and covariance matrix of $b_{\hat{h}_{lk}}(\hat{h}_{lk})$. From Lemma 2, we have

$$ \begin{aligned} \mathbf{C}_{\hat{h}_{lk}; \Psi_{2,lt}} - \mathbf{C}_{\hat{h}_{lk}} &= \mathbf{C}_{\hat{h}_{lk}; \Psi_{2,lt}}^2 (\mathbf{C}_{\hat{h}_{lk}; \Psi_{2,lt}} + \mathbf{C}_{\Psi_{2,lt}; \hat{h}_{lk}})^{-1} \\ \mathbf{m}_{\hat{h}_{lk}; \Psi_{2,lt}} - \mathbf{m}_{\hat{h}_{lk}} &= \mathbf{C}_{\hat{h}_{lk}; \Psi_{2,lt}} (\mathbf{C}_{\hat{h}_{lk}; \Psi_{2,lt}} + \mathbf{C}_{\Psi_{2,lt}; \hat{h}_{lk}})^{-1} \\ & \quad \cdot (\mathbf{m}_{\hat{h}_{lk}; \Psi_{2,lt}} - \mathbf{m}_{\Psi_{2,lt}; \hat{h}_{lk}}). \end{aligned} \quad (32) $$

It has been shown in the proof of Lemma 2 that the difference $\mathbf{m}_{\hat{h}_{lk}}; \Psi_{2,lt} - \mathbf{m}_{\Psi_{2,lt}; \hat{h}_{lk}}$ is a higher order infinitesimal relative to $\mathbf{m}_{\hat{h}_{lk}}; \Psi_{2,lt}$. Thus, the quotients $(\mathbf{C}_{\hat{h}_{lk}; \Psi_{2,lt}} - \mathbf{C}_{\hat{h}_{lk}}) \mathbf{C}_{\hat{h}_{lk}; \Psi_{2,lt}}$ and $(\mathbf{m}_{\hat{h}_{lk}; \Psi_{2,lt}} - \mathbf{m}_{\hat{h}_{lk}})/\mathbf{m}_{\hat{h}_{lk}; \Psi_{2,lt}}$ tend to zero as the system grows larger. Therefore, the difference in (32) are higher order infinitesimals relative to $\mathbf{C}_{\hat{h}_{lk}; \Psi_{2,lt}}$ and $\mathbf{m}_{\hat{h}_{lk}; \Psi_{2,lt}}$, respectively. Therefore, we have $b_{\hat{h}_{lk}}(\hat{h}_{lk}) \approx \mu_{\hat{h}_{lk}; \Psi_{2,lt}}(\hat{h}_{lk})$. Based on the above discussion, we propose to replace the extrinsic (7) at $\Psi_{2,lt}$ by (33) and (34) to reduce complexity further,

$$ \mu'_{\hat{h}_{lk}; \Psi_{2,lt}}(\hat{h}_{lk}) = b_{\hat{h}_{lk}}(\hat{h}_{lk}); \quad (33) $$

$$ \mu'_{x_{kt}; \Psi_{2,lt}}(x_{kt}) = b_{x_{kt}}(x_{kt}). \quad (34) $$

## VII. DECENTRALIZED METHOD

To obtain the belief of $x_{kt}$, we need to combine the message from all the AP. We consider the case where all the $L$ AP are connected, and the AP network has a tree structure. A decentralized message-passing method can be used based on the consensus propagation framework [8]. Define the normalized message from AP $l$ to AP $l'$:

$$ \nu_{l \to l'}(x_{kt}) \propto \mu_{\Psi_{2,lt};x_{kt}}(x_{kt}) \prod_{\overline{l' \in N(l)}/\{l'\}} \nu_{\overline{l} \to l}(x_{kt}), $$

where $N(l)$ denotes the set of connected neighbors of AP $l$. At convergence, the belief in (30) can be obtained by any AP $l$ as
**Algorithm 1 One Iteration of Decentralized EP**

Require: $\Xi_{h_{lk}}$, $y_{p,lg}$, $y_{lt}$, $p(x_{kt})$, $\sigma_x^2$, $\sigma_v^2$, $G_g$

1: Initialize $\mu_{\Psi_{3,lg};h_{lk}}$, $\mu_{\Psi_{2,lt};x_{kt}}$, $\mu_{\Psi_{2,lt};h_{lk}}$, $\nu_{l \to l'}(x_{kt})$

2: At all the APs, $\forall k, t$, update $b_{x_{kt}}$ according to (35)

3: **for** l=1:L **do**

4: &nbsp;&nbsp;&nbsp;&nbsp; $\forall k$, update $\mu_{h_{lk};\Psi_{3,lg}}$ based on (20)

5: &nbsp;&nbsp;&nbsp;&nbsp; $\forall k, t$, update $\mu'_{h_{lk};\Psi_{2,lt}}$ based on (30) and (34)

6: &nbsp;&nbsp;&nbsp;&nbsp; $\forall k, t$, update $\mu'_{x_{kt};\Psi_{2,lt}}$ based on (33)

7: &nbsp;&nbsp;&nbsp;&nbsp; $\forall k$, update $\mu_{\Psi_{3,lg};h_{lk}}$ based on (25)-(26)

8: &nbsp;&nbsp;&nbsp;&nbsp; $\forall k, t$, update $\mu_{\Psi_{2,lt};x_{kt}}$ based on (11)-(13)

9: &nbsp;&nbsp;&nbsp;&nbsp; $\forall k, t$, update $\mu_{\Psi_{2,lt};h_{lk}}$ based on (17)-(19)

10: &nbsp;&nbsp;&nbsp;&nbsp; $\forall l' \in N(l)$, $k, t$, update $\nu_{l \to l'}(x_{kt})$ based on (36)

11: **end for**

$$b'_{x_{kt}}(x_{kt}) \propto p(x_{kt})\mu_{\Psi_{2,lt};x_{kt}}(x_{kt}) \prod_{l' \in N(l)} \nu_{l' \to l}(x_{kt}). \quad (35)$$

Therefore, for a decentralized algorithm, we can replace the update of belief $b_{x_{kt}}$ in (30) by $b'_{x_{kt}}$ in (35). After updating the message $\mu_{\Psi_{2,lt};x_{kt}}$, we update the shared message by

$$\nu^{new}_{l \to l'}(x_{kt}) \propto \mu^{new}_{\Psi_{2,lt};x_{kt}}(x_{kt}) \prod_{\bar{v} \in N(l)/\{l'\}} \nu^{old}_{\bar{v} \to l}(x_{kt}), \quad (36)$$

where we use *new* and *old* to distinguish the message of different iterations. One possible ordering method is suggested in Algorithm 1.

## VIII. SIMULATION RESULTS

Our study simulates an environment within a 400 × 400 square meter area, equipped with 16 APs and 8 User Terminals (UTs). Each AP features $N = 2$ antennas and is positioned at coordinates $(\frac{400}{3}i, \frac{400}{3}j)$, $i, j \in \{0, 1, 2, 3\}$. The UTs are uniformly distributed throughout the area. We denote the distance between each UT $k$ and AP $l$ as $d_{lk}$. Channel covariances for each user $k$ at AP $l$ are modeled using $N \times N$ diagonal matrices, represented as $\sigma_{h_{lk}}^2 \mathbf{I}$, where $10 \log_{10}(\sigma_{h_{lk}}^2) = -30 - 36.7 \log_{10}(d_{lk})$.

All the neighboring APs within $\frac{400}{3}$ meters are connected and can exchange information of the estimated data symbols. Furthermore, as illustrated in Algorithm 1, a synchronized message-exchanging scheme is used.

The length of the orthogonal pilot sequences is set to $P = 6$ to introduce pilot contamination.

We employ a 4QAM constellation of length $T = 10$ for signal transmission and assume a noise power of $-96$ dBm. The signal-to-noise ratio (SNR) is adjusted by varying the transmitted power. We base our results on 100 different realizations, which are illustrated in Figure 2. The normalized mean squared error (NMSE) of the channel estimates is defined as $\text{NMSE} = \frac{\text{tr}[(\hat{H}-H)^2]}{\text{tr}[\hat{H}^2]}$, where $\hat{H}$ are synthesized from the mean of $b_{h_{lk}}(h_{lk})$ defined in (30) and the operation $| \cdot |^2$ is defined as $|\hat{H}|^2 = \hat{H}^H \hat{H}$.

In the VL-EP scenario, we generate data symbols drawn from i.i.d. Gaussian distribution and apply the VL-EP algorithm [5] for channel estimation. In the Genie-Aided scenario, we implement the proposed algorithm as if the data symbols are known. In the MMSE Genie-Aided scenario, all the APs

<figure><img src="image_3.png" /><figcaption>Fig. 2. NMSE vs SNR</figcaption></figure>

estimate the channel coefficients using the MMSE estimator with known channel coefficients.

## IX. CONCLUSIONS

This paper introduces a simplified, decentralized EP-based algorithm for bilinear joint estimation. To simplify the factorization scheme, we leverage orthogonal pilots and the CLT. Through asymptotic analysis, we further refine the message update scheme within the algorithm. Although originally developed for an acyclic network of APs, our simulation results confirm the algorithm's effectiveness even when the APs are interconnected in a cyclic network.

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

## REFERENCES

[1] Z. Zhao and D. Slock, "Decentralized Expectation Propagation for Semi-Blind Channel Estimation in Cell-Free Networks," in *IEEE International Symposium on Information Theory*, IEEE, Ed., 2024.

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

[3] M. J. Wainwright, M. I. Jordan et al., "Graphical Models, Exponential Families, and Variational Inference," *Foundations and Trends in Machine Learning*, 2008.

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

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

[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] K. Takeuchi, "Decentralized Generalized Approximate Message-Passing for Tree-Structured Networks," *arXiv preprint arXiv:2311.01032*, 2023.

[8] C. C. Moallemi and B. Van Roy, "Consensus Propagation," *IEEE Transactions on Information Theory*, 2006.

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

[10] K.-H. Ngo, M. Guillaud, A. Decurninge, S. Yang, and P. Schniter, "Multi-User Detection Based on Expectation Propagation for the Non-Coherent SIMO Multiple Access Channel," *IEEE Transactions on Wireless Communications*, 2020.

[11] Q. Zou, H. Zhang, C.-K. Wen, S. Jin, and R. Yu, "Concise Derivation for Generalized Approximate Message Passing Using Expectation Propagation," *IEEE Signal Processing Letters*, 2018.

[12] Z. Zhao, F. Xiao, and D. Slock, "Approximate message passing for not so large niid generalized linear models," in *International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)*, IEEE, Ed., 2023.

[13] P. Schniter, "A Simple Derivation of AMP and its State Evolution via First-Order Cancellation," *IEEE Transactions on Signal Processing*, 2020.