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

Zilu Zhao, Dirk Slock

► To cite this version:

Zilu Zhao, Dirk Slock. Hierarchical Expectation Propagation for Semi-Blind Channel Estimation in Cell-Free Networks. ICASSP 2025, IEEE International Conference on Acoustics, Speech and Signal Processing, IEEE, Apr 2025, Hyderabad, India. 10.1109/ICASSP49660.2025.10889567. hal-05076627

HAL Id: hal-05076627

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

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
# Hierarchical 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**—In this work, we study uplink communication in cell-free (CF) massive multiple-input multiple-output (MaMIMO) systems, a promising architecture for next-generation networks. To address the challenge of pilot contamination, we employ semi-blind transmission structures that enable joint channel and data symbol estimation. However, Bayesian estimation in such semi-blind frameworks leads to intractable bilinear problems. To tackle this, we propose a simplified, distributed method based on Expectation Propagation (EP) for efficient semi-blind channel estimation. Notably, we identify that if the data constellation set can be decomposed into multiple sub-constellation sets with identical amplitudes, this structure can be leveraged to significantly reduce computational complexity. This approach is particularly advantageous for managing large constellation sizes, ensuring scalability and efficiency in practical systems. Additionally, approximations based on the Central Limit Theorem are incorporated to further simplify computations.

## I. INTRODUCTION

Unlike conventional mobile networks, Cell-Free (CF) networks no longer have the cellular topology. As a result, all the user terminals (UTs) will be served by all the access points (APs) in an area. A huge gain throughput gain can be obtained by transforming into the new CF network topology [1]. One of the challenges introduced by CF network is pilot contamination, which occurs when the number of UTs exceeds the length of pilot sequence. As a result, multiple UTs may use the same pilot sequence.

Semi-Blind approaches [2] have been explored to mitigate the effects of pilot contamination. In Semi-Blind approaches, the channel and data are jointly estimated with the presence of pilot contamination. The authors have also shown a performance gain by transforming from deterministic estimation to Bayesian estimation [2], [3]. However, Bayesian inference is intractable due to the high dimensional integrals.

A variety of Message Passing based algorithms has been proposed to simplify the high dimensional Bayesian inference problem into local low dimensional inference problems [4]–[8].

### A. Continuous Input

Variable Level Expectation Propagation (VL-EP) [9] is a variant of Expectation Propagation (EP) that assumes Gaussian input from UTs and combines EP with Expectation Maximization (EM) for joint bilinear estimation. In VL-EP, posterior distributions (referred to as beliefs in EP terminology) are used as extrinsic messages to reduce complexity.

In [10], the authors investigated the differences between VL-EP and EP algorithms and proposed a block coordinate descent method based on the Majorization-Minimization (MM) algorithm. However, in more realistic scenarios, UT inputs typically follow a discrete distribution, such as Quadrature Amplitude Modulation (QAM). Unfortunately, the aforementioned methods do not perform well for discrete input cases.

### B. Discrete Input

In [11], the authors proposed an EP-based distributed algorithm for bilinear estimation with discrete UT inputs. However, the complexity of this approach scales linearly with the size of the constellation set. In [12], the authors exploited the relationship between Bethe Free Energy (BFE) [4] and message-passing algorithms. By incorporating mean-field constraints into the BFE, the resulting low-complexity algorithm, proposed in [12], can be interpreted as a combination of Variational Bayes (VB) and EP.

### C. Main Contributions

In this paper, we propose an algorithm to mitigate the high complexity associated with large constellation sets. As observed in [13] and [11], differences in the amplitudes of constellation symbols further contribute to this complexity. Additionally, in decentralized or distributed algorithms, access points (APs) need to share information about data symbols through backhaul connections. This information is represented as a list of probabilities corresponding to each element in the constellation set, and its size grows linearly with the size of the constellation set. While the approach in [12] ensures that backhaul information exchanged between APs is independent of the constellation size, it still requires traversing the entire constellation set to determine the transmitted symbol.

To address these challenges, we exploit the hierarchical structure of the constellation set. Specifically, we assume that the original constellation set can be decomposed into multiple sub-constellation sets, each with identical amplitudes. A commonly used example is $4^M$-QAM, which can be decomposed into $M$ 4-QAM sub-constellation sets. By leveraging this decomposition, backhaul messages become Gaussian, requiring only the first and second-order moments to be shared between APs. Furthermore, this approach reduces the complexity of symbol determination from $4^M$ to $4M$ possible values per data symbol.

## II. SYSTEM MODEL

We consider a semi-blind signal model containing $L$ APs and $K$ UTs. At the $l$-th AP, it is represented by

$$[\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}^{M \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}^{NK \times NK}$ is a block diagonal matrix of $K$ blocks $\Xi_{hl_k} \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. 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} \ \mathbf{V}_l]) \sim \mathcal{CN}(0, \mathbf{C}_v)$ with $\mathbf{C}_v = \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 $\tilde{\mathbf{x}}_{p,g}^*$ which is the conjugated g-th pilot sequence. This results in an equivalent observation $\tilde{\mathbf{y}}_{p,lg}$

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

where $\tilde{\mathbf{v}}_{p,lg} = \frac{1}{P\sigma_x^2} \mathbf{V}_{p,l} \tilde{\mathbf{x}}_{p,g}^* \sim \mathcal{N}(\mathbf{v}_{p,lg}|\mathbf{0}, \mathbf{C}_{\tilde{\mathbf{v}}})$, $\mathbf{C}_{\tilde{\mathbf{v}}} = \frac{\sigma_v^2 \mathbf{I}}{P\sigma_x^2}$, and $G_g$ denotes 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}[\tilde{\mathbf{v}}_{p,lg} \tilde{\mathbf{v}}_{p,lg'}^H]$ is an all-zero matrix for all $g \neq g'$. Therefore, the observations $\tilde{\mathbf{y}}_{p,lg}$ and $\tilde{\mathbf{y}}_{p,lg'}$ are independent.

### B. Hierarchical Data Symbol

In this paper, we assume that the data symbols can be decomposed into

$$ x_{kt} = \sum_{m=1}^{M} x_{m,kt}, \quad (3) $$

where $|x_{m,kt}|^2 = \sigma_{xm}^2$. A common example of this constellation scheme is $4^M$-QAM. If we assume the transmission power of $4^M$-QAM to be $\sigma_x^2$, then we have

$$ \sigma_{xm}^2 = \frac{3 \sigma_x^2 \frac{1}{4}^{m-1}}{4 \left(1 - \frac{1}{4}^M\right)}. \quad (4) $$

We define an auxiliary variable

$$ \mathbf{z}_{m,lt} = \sum_k \mathbf{h}_{lk} x_{m,kt}. \quad (5) $$

The joint pdf can be factored as

$$ \begin{aligned} & p(\{\mathbf{y}_{p,lg}\}, \{\mathbf{Y}_l\}, \{\mathbf{z}_{m,lt}\}, \{\mathbf{H}_l\}, \mathbf{X}, \{\mathbf{V}_l\}) \\ & = \prod_l \prod_t f_{y_{lt}}(\mathbf{z}_{\{m\},lt}) \prod_m \prod_l \prod_t \delta_{m,lt}(\mathbf{z}_{m,lt}, \mathbf{h}_{l\{k\}}, \mathbf{x}_{m,\{k\}t}) \\ & \quad \cdot \prod_l \prod_g f_{H_{lG_g}}(\mathbf{H}_{lG_g}) \prod_m \prod_k \prod_t f_{x_{m,kt}}(x_{m,kt}) \end{aligned} \quad (6) $$

where

$$ \begin{aligned} f_{y_{lt}}(\mathbf{z}_{\{m\},lt}) &= p\left(\mathbf{y}_{lt} \middle| \sum_m \mathbf{z}_{m,lt}\right); \\ \delta_{m,lt}(\mathbf{z}_{m,lt}, \mathbf{h}_{l\{k\}}, \mathbf{x}_{m,\{k\}t}) &= p\left(\mathbf{z}_{m,lt} \middle| \sum_k \mathbf{h}_{lk} x_{m,kt}\right); \quad (7) \\ f_{H_{lG_g}}(\mathbf{H}_{lG_g}) &= p(\mathbf{y}_{p,lg}, \mathbf{H}_{lG_g}); \\ f_{x_{m,kt}}(x_{m,kt}) &= p(x_{m,kt}) \end{aligned} $$

where the curly bracket in the subscript $\theta_{\{n\}}$ denotes all $\theta_1, \dots, \theta_N$, and $\mathbf{H}_{lG_g}$ denotes a matrix collecting all $k \in G_g$, $\mathbf{h}_{lk}$ as its column vectors. 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 [14]. 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$: [15]

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

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 (9) $$

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

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

where $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 following (8). The only difference between BP and EP lies in the projection step: in BP, the argument is projected onto itself, i.e., in BP, we have $b_{\Psi} \in Q$. The projection step can be viewed as an approximation. To maintain high precision while ensuring low complexity, we use BP for estimating data symbol $x_{m,kt}$ and EP for estimating the Channel $\mathbf{h}_{lk}$. We exploit the Central Limit Theorem (CLT) when estimating $\mathbf{z}_{m,lt}$, resulting in BP being equivalent to EP for this estimation.

### A. Notations

Due to the relationship described in (8), the variable-to-factor messages are fully determined by the factor-to-variable messages. Therefore, in the following discussion, we will focus on deriving the factor-to-variable messages. In the EP procedure, we denote the mean and covariance of the variable-to-factor message $\mu_{\Psi;\theta_i}$ as $\mathbf{m}_{\Psi;\theta_i}$ and $\mathbf{C}_{\Psi;\theta_i}$, and those of the factor-to-variable message $\mu_{\theta_i;\Psi}$ as $\mathbf{m}_{\theta_i;\Psi}$ and $\mathbf{C}_{\theta_i;\Psi}$.

## IV. MESSAGE PASSING DERIVATIONS

We first look at the factor node $\delta_{m,lt}$. The (partially) marginalized belief at $\delta_{m,lt}$ is

$$ \begin{aligned} & b_{\delta_{m,lt}}(\mathbf{z}_{m,lt}, \mathbf{h}_{lk}, x_{m,kt}) \\ &= \int \delta(\mathbf{z}_{m,lt} - \sum_k \mathbf{h}_{lk} x_{m,kt}) \mu_{\mathbf{z}_{m,lt};\delta_{m,lt}}(\mathbf{z}_{m,lt}) \\ & \quad \cdot \prod_k \mu_{\mathbf{h}_{lk};\delta_{m,lt}}(\mathbf{h}_{lk}) \prod_k \mu_{x_{m,kt};\delta_{m,lt}}(x_{m,kt}) d\mathbf{h}_{l\bar{k}} dx_{m,l\bar{k}} \\ &= \mu_{\mathbf{z}_{m,lt};\delta_{m,lt}}(\mathbf{z}_{m,lt}) \mu_{\mathbf{h}_{lk};\delta_{m,lt}}(\mathbf{h}_{lk}) \mu_{x_{m,kt};\delta_{m,lt}}(x_{m,kt}) \\ & \quad \cdot \int d\mathbf{w}_{m,lkt} \delta(\mathbf{z}_{m,lt} - \mathbf{w}_{m,lkt} - \mathbf{h}_{lk} x_{m,kt}) \\ & \quad \cdot \int d\mathbf{h}_{l\bar{k}} dx_{m,l\bar{k}} \delta(\mathbf{w}_{m,lkt} - \sum_{k' \neq k} \mathbf{h}_{lk'} x_{m,k't}) \\ & \quad \cdot \prod_{k' \neq k} \mu_{\mathbf{h}_{lk'};\delta_{m,lt}}(\mathbf{h}_{lk'}) \prod_{k' \neq k} \mu_{x_{m,k't};\delta_{m,lt}}(x_{m,k't}), \end{aligned} \quad (11) $$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Factor Graph with Neighborhoods of both Factor (above) and Variable (below) Nodes. For simplicity, $\forall$ symbol is used to indicate that there are multiple nodes.</figcaption></figure>

where we define $w_{m,lkt} = \sum_{k' \neq k} h_{lk'} x_{m,k't}$ as the interference in (11). If we look at the second integral after the second equal sign in (11), this interference $w_{m,lkt}$ can be approximated as Gaussian by Central Limit Theory, with mean and covariance matrix

$$ 
\begin{aligned}
\mathbf{m}_{w_{m,lkt}} &= \sum_{k' \neq k} \mathbf{m}_{h_{lk'}} \delta_{m,lt} m_{x_{m,k'}t} \delta_{m,lt} \\
\mathbf{C}_{w_{m,lkt}} &= \sum_{k' \neq k} \mathbf{C}_{h_{lk'}} \delta_{m,lt} \sigma_{xm}^2 \\
&\quad + \mathbf{m}_{h_{lk'}} \delta_{m,lt} \mathbf{m}_{h_{lk'}}^H \delta_{m,lt} \tau_{x_{m,k'}t} \delta_{m,lt}.
\end{aligned} 
\quad (12) $$

With the above CLT approximation, the term $z_{m,lt} - w_{m,lkt}$ in (11) can be verified to be a Gaussian with mean and covariance matrix:

$$ 
\begin{aligned}
\mathbf{m}_{u_{m,lkt}} &= \mathbf{m}_{z_{m,lt}} \delta_{m,lt} - \mathbf{m}_{w_{m,lkt}} \\
\mathbf{C}_{u_{m,lkt}} &= \mathbf{C}_{z_{m,lt}} \delta_{m,lt} + \mathbf{C}_{w_{m,lkt}}
\end{aligned} 
\quad (13) $$

### A. Message from $\delta_{m,lt}$ to $x_{m,kt}$

From the derivation in [13] and the above CLT approximations, the message $\mu_{\delta_{m,lt};x_{m,kt}}$ according to (8) is Gaussian with mean and variance:

$$ 
\begin{aligned}
\tau_{\delta_{m,lt};x_{m,kt}} &= \left[ \mathbf{m}_{h_{lk}}^H \delta_{m,lt} \left( \mathbf{C}_{u_{m,lkt}} + \sigma_{xm}^2 \mathbf{C}_{h_{lk}} \delta_{m,lt} \right)^{-1} \right. \\
&\qquad \left. \cdot \mathbf{m}_{h_{lk}} \delta_{m,lt} \right]^{-1} \\
\mu_{\delta_{m,lt};x_{m,kt}} &= \tau_{\delta_{m,lt};x_{m,kt}} \mathbf{m}_{h_{lk}}^H \delta_{m,lt} \\
&\qquad \cdot \left( \mathbf{C}_{u_{m,lkt}} + \sigma_{xm}^2 \mathbf{C}_{h_{lk}} \delta_{m,lt} \right)^{-1} \mathbf{m}_{u_{m,lkt}}.
\end{aligned} 
\quad (14) $$

The belief of $x_{m,kt}$ at $\delta_{m,lt}$ factor node is:

$$ b_{\delta_{m,lt};x_{m,kt}}(x_{m,kt}) = \mu_{x_{m,kt};\delta_{m,lt}}(x_{m,kt})\mu_{\delta_{m,lt};x_{m,kt}}(x_{m,kt}) $$

Since BP is used for $x_{m,kt}$, we have:

$$ b_{x_{m,kt}}(x_{m,kt}) = b_{\delta_{m,lt};x_{m,kt}}(x_{m,kt}) \quad (15) $$

In most cases, $4^M$-QAM is used. To obtain the belief, we first look at the message $\delta_{m,lt} \rightarrow x_{m,kt}$:

$$ 
\begin{aligned}
\tau_{x_{m,kt}|\mathbf{y}} &= \left( \sum_l \tau_{\delta_{m,lt};x_{m,kt}}^{-1} \right)^{-1} \\
m_{x_{m,kt}|\mathbf{y}} &= \tau_{x_{m,kt}|\mathbf{y}} \left( \sum_l \tau_{\delta_{m,lt};x_{m,kt}}^{-1} m_{\delta_{m,lt};x_{m,kt}} \right).
\end{aligned} 
\quad (16) $$

Therefore,

$$ b_{x_{m,kt}}(x_{m,kt}) = \frac{p(x_{m,kt}) \mathcal{CN}(x_{m,kt} | m_{x_{m,kt}|\mathbf{y}}, \tau_{x_{m,kt}|\mathbf{y}})}{Z_{x_{m,kt}}}, \quad (17) $$

where:

$$ Z_{x_{m,kt}} = \sum_{x_{m,kt} \in S_m} p(x_{m,kt}) \mathcal{CN}(x_{m,kt} | m_{x_{m,kt}|\mathbf{y}}, \tau_{x_{m,kt}|\mathbf{y}}). \quad (18) $$

In the following, we denote the mean and variance of $b_{x_{m,kt}}$ as $m_{b_{x_{m,kt}}}$ and $\tau_{b_{x_{m,kt}}}$.

### B. Message from $\delta_{m,lt}$ to $h_{lk}$

With the CLT approximation, the joint belief pdf (11) can be written as:

$$ 
\begin{aligned}
b_{\delta_{m,lt}}(h_{lk}, x_{m,kt}) &= \mathcal{CN}(h_{lk} | \mathbf{m}_{h_{lk}} | x_{m,kt}, \mathbf{C}_{h_{lk}} | \sigma_{xm}^2) b_{\delta_{m,lt};x_{m,kt}}(x_{m,kt}),
\end{aligned} 
\quad (19) $$

where:

$$ 
\begin{aligned}
\mathbf{C}_{h_{lk}} | \sigma_{xm}^2 &= (\sigma_{xm}^2 \mathbf{C}_{u_{m,lkt}}^{-1} + \mathbf{C}_{h_{lk}}^{-1} \delta_{m,lt})^{-1} \\
\mathbf{m}_{h_{lk}} | x_{m,kt} &= \mathbf{C}_{h_{lk}} | \sigma_{xm}^2 \left( x_{m,kt}^* \mathbf{C}_{u_{m,lkt}}^{-1} \mathbf{m}_{u_{m,lkt}} + \mathbf{C}_{h_{lk}}^{-1} \delta_{m,lt} \mathbf{m}_{h_{lk}} \delta_{m,lt} \right)
\end{aligned} 
\quad (20) $$

The mean and covariance matrix of belief $b_{\delta_{m,lt};h_{lk}}$ for $h_{lk}$ at $\delta_{m,lt}$ can be obtained as [13]:

$$ 
\begin{aligned}
\mathbf{m}_{b_{\delta_{m,lt};h_{lk}}} &= \mathbf{C}_{h_{lk}} | \sigma_{xm}^2 \left( m_{b_{\delta_{m,lt};h_{lk}}}^* \mathbf{C}_{u_{m,lkt}}^{-1} \mathbf{m}_{u_{m,lkt}} + \mathbf{C}_{h_{lk}}^{-1} \delta_{m,lt} \mathbf{m}_{h_{lk}} \delta_{m,lt} \right) \\
\mathbf{C}_{b_{\delta_{m,lt};h_{lk}}}^{-1} &= \mathbf{C}_{h_{lk}}^{-1} | \sigma_{xm}^2 - \tau_{b_{x_{m,kt}}} \mathbf{C}_{u_{m,lkt}}^{-1} \mathbf{m}_{u_{m,lkt}} \cdot \left( 1 + \tau_{b_{x_{m,kt}}} \mathbf{m}_{u_{m,lkt}}^H \mathbf{C}_{u_{m,lkt}}^{-1} \mathbf{C}_{h_{lk}} | \sigma_{xm}^2 \mathbf{C}_{u_{m,lkt}}^{-1} \mathbf{m}_{u_{m,lkt}} \right)^{-1} \\
&\quad \cdot \mathbf{m}_{u_{m,lkt}}^H \mathbf{C}_{u_{m,lkt}}^{-1} \mathbf{m}_{u_{m,lkt}}.
\end{aligned} 
\quad (21) $$

The feedback message can be obtained by $\frac{\text{proj}[b_{\delta_{m,lt};h_{lk}}]}{\mu_{h_{lk}} \delta_{m,lt}}$:

$$ 
\begin{aligned}
\mathbf{C}_{\delta_{m,lt};h_{lk}} &= \frac{\mathbf{C}_{u_{m,lkt}}}{\sigma_{xm}^2} - \tau_{b_{x_{m,kt}}} \mathbf{m}_{u_{m,lkt}} \cdot \left[ \tau_{b_{x_{m,kt}}} \gamma_{aux_{m,lkt}} - \sigma_{xm}^4 \right]^{-1} \mathbf{m}_{u_{m,lkt}}^H \\
\mathbf{m}_{\delta_{m,lt};h_{lk}} &= \mathbf{C}_{\delta_{m,lt};h_{lk}} \mathbf{C}_{b_{\delta_{m,lt};h_{lk}}}^{-1} \mathbf{m}_{b_{\delta_{m,lt};h_{lk}}} - \mathbf{C}_{h_{lk}}^{-1} \delta_{m,lt} \mathbf{m}_{h_{lk}} \delta_{m,lt},
\end{aligned} 
\quad (22) $$

where:

$$ \gamma_{aux_{m,lkt}} = \mathbf{m}_{u_{m,lkt}}^H \left( \frac{\mathbf{C}_{u_{m,lkt}}}{\sigma_{xm}^2} + \mathbf{C}_{h_{lk}} \delta_{m,lt} \right)^{-1} \mathbf{m}_{u_{m,lkt}} \quad (23) $$
### C. Message from $\delta_{m,lt}$ to $\mathbf{z}_{m,lt}$

Similar to the computation of $w_{m,lkt}$ in (12), we approximate the marginalization over $\forall k$, $\mathbf{h}_{lk}$ and $x_{m,lkt}$ by CLT. Therefore, from (11), we have

$$ 
\begin{aligned}
b_{\delta_{m,lt}}(\mathbf{z}_{m,lt}) &= \mu_{\mathbf{z}_{m,lt};\delta_{m,lt}}(\mathbf{z}_{m,lt}) \\
&\quad \cdot \int \delta(\mathbf{z}_{m,lt} - \sum_k \mathbf{h}_{lk}x_{m,kt}) \prod_k \mu_{\mathbf{h}_{lk};\delta_{m,lt}}(\mathbf{h}_{lk}) \\
&\quad \cdot \prod_k \mu_{x_{m,kt};\delta_{m,lt}}(x_{m,kt})d\mathbf{h}_{lk}dx_{m,kt}
\end{aligned} 
\quad (24) $$

According to CLT, the integral in (24) can be approximated as Gaussian. Following (8), we can verify that the integral part is the feedback message, and thus,

$$ 
\begin{aligned}
\mathbf{m}_{\delta_{m,lt};\mathbf{z}_{m,lt}} &= \frac{1}{K-1} \sum_k \mathbf{m}_{\mathbf{w}_{m,l,k,t}} \\
\mathbf{C}_{\delta_{m,lt};\mathbf{z}_{m,lt}} &= \frac{1}{K-1} \sum_k \mathbf{C}_{\mathbf{w}_{m,l,k,t}}.
\end{aligned} 
\quad (25) $$

### D. Message from $f_{H_{lG_g}}$ to $\mathbf{h}_{lk}$

Since $f_{H_{lG_g}}$ is already Gaussian, the belief at $f_{H_{lG_g}}$ must be Gaussian,

$$ b_{f_{H_{lG_g}}} = f_{H_{lG_g}}(H_{lG_g}) \prod_{k \in G_g} \mu_{\mathbf{h}_{lk};f_{H_{lG_g}}}(\mathbf{h}_{lk}). \quad (26) $$

Following [13], we get

$$ 
\begin{aligned}
& \mu_{f_{H_{lG_g}};\mathbf{h}_{lk}}(\mathbf{h}_{lk}) \\
&= \mathcal{CN} \left( \mathbf{h}_{lk} | \tilde{\mathbf{y}}_{p,lg} - \sum_{k' \in G_g / \{k\}} \mathbf{m}_{\mathbf{h}_{lk'}|y}, \mathbf{C}_{\tilde{\mathbf{v}}} + \sum_{k' \in G_g / \{k\}} \mathbf{C}_{\mathbf{h}_{lk'}|y} \right) \\
&\quad \cdot \mathcal{CN}(\mathbf{h}_{lk} | \mathbf{0}, \Xi_{\mathbf{h}_{lk}}),
\end{aligned} 
\quad (27) $$

where

$$ 
\begin{aligned}
\mathbf{C}_{\mathbf{h}_{lk}|y} &= \left( \mathbf{C}_{\mathbf{h}_{lk};f_{H_{lG_g}}}^{-1} + \Xi_{\mathbf{h}_{lk}}^{-1} \right)^{-1} \\
\mathbf{m}_{\mathbf{h}_{lk}|y} &= \mathbf{C}_{\mathbf{h}_{lk}|y} \mathbf{C}_{\mathbf{h}_{lk};f_{H_{lG_g}}}^{-1} \mathbf{m}_{\mathbf{h}_{lk};f_{H_{lG_g}}}.
\end{aligned} 
\quad (28) $$

### E. Message from $f_{y_{lt}}$ to $\mathbf{z}_{m,lt}$

Since belief at $f_{y_{lt}}$ is

$$ b_{f_{y_{lt}}} = p(y_{lt} | \sum_m \mathbf{z}_{m,lt}) \prod_m \mu_{\mathbf{z}_{m,lt};f_{y_{lt}}}(\mathbf{z}_{m,lt}), \quad (29) $$

which is already a Gaussian, we can immediately derive the feedback message by Gaussian Reproduction Lemma [15]:

$$ 
\begin{aligned}
\mathbf{C}_{f_{y_{lt}};\mathbf{z}_{m,lt}} &= \mathbf{C}_{\mathbf{v}} + \sum_{m' \neq m} \mathbf{C}_{\mathbf{z}_{m,lt};f_{y_{lt}}} \\
\mathbf{m}_{f_{y_{lt}};\mathbf{z}_{m,lt}} &= \mathbf{y}_{lt} - \sum_{m' \neq m} \mathbf{m}_{\mathbf{z}_{m,lt};f_{y_{lt}}}
\end{aligned} 
\quad (30) $$

## V. SIMULATION RESULTS

We simulate an environment within a 400 × 400 square meter area, equipped with 16 APs and 8 User Terminals. 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

### Algorithm 1 One Iteration of Hierachical EP

Require: $\Xi_{\mathbf{h}_{lk}}, \tilde{\mathbf{y}}_{p,lg}, \mathbf{y}_{lt}, \sigma_x^2, \sigma_v^2, G_g, S$

1. Initialize $\mu_{\delta_{m,lt};\mathbf{h}_{lk}}, \mu_{\delta_{m,lt};\mathbf{z}_{m,lt}}, \mu_{\delta_{m,lt};x_{m,kt}}$

2. $\mathbf{C}_{f_{H_{lG_g}};\mathbf{h}_{lk}}$ and $\mathbf{m}_{f_{H_{lG_g}};\mathbf{h}_{lk}}$ via (28) → (27)

3. $\mathbf{C}_{f_{y_{lt}};\mathbf{z}_{m,lt}}$ and $\mathbf{m}_{f_{y_{lt}};\mathbf{z}_{m,lt}}$ via (30)

4. $\mathbf{C}_{\mathbf{u}_{m,lkt}}$ and $\mathbf{m}_{\mathbf{u}_{m,lkt}}$ via (13)

5. $\tau_{\delta_{m,lt};x_{m,kt}}$ and $\mu_{\delta_{m,lt};x_{m,kt}}$ via (14)

6. $b_{x_{m,kt}}(x_{m,kt})$ via (16) → (18) → (17)

7. $\mathbf{C}_{\mathbf{h}_{lk}|\sigma_x^2}$ via (20)

8. $\mathbf{C}_{b_{\delta_{m,lt};\mathbf{h}_{lk}}^{-1}}$ and $\mathbf{m}_{b_{\delta_{m,lt};\mathbf{h}_{lk}}}$ via (21)

9. $\gamma_{aux_{m,lkt}}$ via (23)

10. $\mathbf{C}_{\delta_{m,lt};\mathbf{h}_{lk}}^{-1} \mathbf{m}_{\delta_{m,lt};\mathbf{h}_{lk}}$ and $\mathbf{P}_{\delta_{m,lt};\mathbf{h}_{lk}}$ via (22)

11. $\mathbf{C}_{\delta_{m,lt};\mathbf{z}_{m,lt}}$ and $\mathbf{m}_{\delta_{m,lt};\mathbf{z}_{m,lt}}$ via (25)

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

matrices, represented as $\sigma_{\mathbf{h}_{lk}}^2 \mathbf{I}$, where $10 \log_{10}(\sigma_{\mathbf{h}_{lk}}^2) = -30 - 36.7 \log_{10}(d_{lk})$.

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

We employ $4^M$ QAM constellation of length $T = 16$ for signal transmission where $M \in \{1, 2\}$. Following [1], we employed power control to ensure the total received power at the APs is the same. We base our results on 50 different realizations, which are illustrated in Figure 2. 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 estimate the channel coefficients using the MMSE estimator with known channel coefficients which is the theoretical lower bound for all channel estimation methods.

## VI. CONCLUSIONS

This paper introduces a simplified, distributed EP-BP-based algorithm for bilinear joint estimation. By leveraging the structure of the data symbols, we reduce the complexity associated with the constellation size. The proposed algorithm is compatible with $4^M$-QAM data symbols. At each AP, the computational complexity is $O(MKT \cdot N^3)$, while the CPU handles the computation of data symbol beliefs with a complexity of $O(MLKT)$.

This approach can be extended to other message-passing algorithms to further reduce complexity. For example, the hybrid VB-EP algorithm in [12] treats the entire data sequence as an atomic variable. In contrast, our method considers each individual data symbol as an atomic variable.

**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] H. Q. Ngo, A. Ashikhmin, H. Yang, E. G. Larsson, and T. L. Marzetta, "Cell-free massive mimo: Uniformly great service for everyone," in *2015 IEEE 16th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)*, 2015, pp. 201–205.

[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] Z. Zhao and D. Slock, "Semi-blind sparse channel learning in cell-free massive mimo - a crb analysis," in *ICC, IEEE International Conference on Communications*, 2023.

[4] ——, "Bethe free energy and extrinsics in approximate message passing," *ASILOMAR, Asilomar Conference on Signals, Systems, and Computers*, 2023.

[5] 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 *ICASSP, IEEE International Conference on Acoustics, Speech and Signal Processing*, 2024.

[6] ——, "Approximate message passing for not so large niid generalized linear models," in *SPAWC, IEEE International Workshop on Signal Processing Advances in Wireless Communications*, 2023.

[7] Z. Zhao and D. Slock, "Variance predictions in vamp/uamp with right rotationally invariant measurement matrices for niid generalized linear models," in *EUSIPCO 2023, 31st European Signal Processing Conference*, 2023.

[8] ——, "Improved variance predictions in approximate message passing," in *MLSP 2023, IEEE 33rd International Workshop on Machine Learning for Signal Processing*, 2023.

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

[10] Z. Zhao and D. Slock, "Bilinear hybrid expectation maximization and expectation propagation for semi-blind channel estimation," in *EUSIPCO, European Signal Processing Conference*, 2024.

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

[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, Asilomar Conference on Signals, Systems, and Computers*, 2024.

[13] ——, "Decentralized expectation propagation for semi-blind channel estimation in cell-free networks," 2024. [Online]. Available: https://arxiv.org/abs/2410.01303

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

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