Cell-Free Massive MIMO Detection: A Distributed Expectation Propagation Approach

As illustrated in Fig. 1, we consider a cell-free massive MIMO network with $L$ distributed APs, each equipped with $N$ antennas to serve single-antenna users. All APs are connected to a CPU that has abundant computing resources. Denote ${\mathbf{h}}_{kl}\in{\mathbb{C}}^{N\times 1}\sim\mathcal{N}_{{\mathbb{C}}}(\mathbf{0},{\mathbf{R}}_{kl})$ as the channel between the $k$-th user and the $l$-th AP, where ${\mathbf{R}}_{kl}\in{\mathbb{C}}^{N\times N}$ is the spatial correlation matrix, and $\beta_{k,l}=\mathrm{tr}({\mathbf{R}}_{kl})/N$ as the large-scale fading coefficient involving the geometric path loss and shadowing. In the uplink data transmission phase, we consider $\mathcal{M}_{k}\subset\{1,\ldots,L\}$ as the subset of APs that serve the $k$-th user and define the DCC matrices ${\mathbf{D}}_{kl}\in{\mathbb{C}}^{N\times N}$ based on $\mathcal{M}_{k}$ as

Specifically, ${\mathbf{D}}_{kl}$ is an identity matrix if the $l$-th AP serves the $k$-th user, and is a zero matrix, otherwise. Furthermore, we define $\mathcal{D}_{l}$ as the set of user indices that are served by the $l$-th AP

where the cardinality of $\mathcal{D}_{l}$ is denoted as $|\mathcal{D}_{l}|$. When each AP serves all users, we have $|\mathcal{D}_{l}|=K$ accordingly. Assuming that perfect CSI is available at the local APs, the received signal ${\mathbf{y}}_{l}\in{\mathbb{C}}^{N\times 1}$ at the $l$-th AP is then given by

where $x_{k}\in{\mathbb{C}}$ is the transmitted symbol drawn from an $M$-QAM constellation, and $p_{k}>0$ is the transmit power at the $k$-th user. The additive noise at the $l$-th AP is denoted as ${\mathbf{n}}_{l}\in{\mathbb{C}}^{N\times 1}\sim\mathcal{N}_{{\mathbb{C}}}(\mathbf{0}_{N},\sigma^{2}{\mathbf{I}}_{N})$. Let ${\mathbf{x}}=[x_{1},\ldots,x_{K}]^{T}$ denote the transmitted vector from all users and ${\mathbf{H}}_{l}=[{\mathbf{h}}_{1l},\ldots,{\mathbf{h}}_{Kl}]\in{\mathbb{C}}^{N\times K}$ is the channel from all users to AP $l$. If the $k$-th user is not associated with the $l$-th AP, the channel vector is ${\mathbf{h}}_{kl}=\mathbf{0}$. Furthermore, we denote ${\mathbf{H}}=[{\mathbf{H}}_{1}^{T},\ldots,{\mathbf{H}}_{L}^{T}]^{T}\in{\mathbb{C}}^{LN\times K}$ as the channel matrix between all users and APs. The uplink detection problem for cell-free massive MIMO is to detect the transmitted data ${\mathbf{x}}$ based on the received signals ${\mathbf{y}}_{l}\,(l=1,2,\ldots,L)$, channel matrix ${\mathbf{H}}$, and noise power $\sigma^{2}$.

For the cell-free massive MIMO detection, there are four commonly-adopted linear receivers [27] with different levels of cooperation among APs. This is elaborated as follows.

• •

  Fully centralized receiver: The pilot and data signals received at all APs are sent to the CPU for channel estimation and data detection. The CPU may perform the MMSE or MRC detection.

• •

  Partially distributed receiver with large-scale fading decoding: First, each AP estimates the channels and uses the linear MMSE detector to detect the received signals. Then, the detected signals are collected at the CPU for joint detection for all users by utilizing the large-scale fading decoding (LSFD) method. Compared to the fully centralized receiver, only the channel statistics are utilized at the CPU but the pilot signals are not required to be sent to the CPU.

• •

  Partially distributed receiver with average decoding: It is a special case of the partially distributed receiver with LSFD. The CPU performs joint detection for all users by simply taking the average of the local estimates. Thus, no channel statistics are required to be transmitted to the CPU via the fronthaul.

• •

  Fully distributed receiver: It is a fully distributed approach in which the data detection is performed at the APs based on the local channel estimates. No information is required to be transferred to the CPU. Each AP may perform MMSE detection.

All of the above-mentioned four receivers achieve performance that is far from the optimal one due to the linear processing. In contrast, non-linear receivers have shown great advantages in terms of BER in cellular massive MIMO systems but at the cost of higher computational complexity [31]. Thanks to the relatively high computing abilities at the CPU and the large number of APs in cell-free massive MIMO systems, we can offload some of the computational-intensive operations to the CPU and distribute the partial computation tasks to APs. Following this idea, we will propose a distributed non-linear detector for cell-free massive MIMO systems.

We first apply the Bayesian inference to recover the signals ${\mathbf{x}}$ from the received signal ${\mathbf{y}}$ in the data detection stage with the linear model ${\mathbf{y}}={\mathbf{H}}{\mathbf{x}}+{\mathbf{n}}$, where ${\mathbf{y}}=[{\mathbf{y}}_{1}^{T},\ldots,{\mathbf{y}}_{L}^{T}]^{T}\in{\mathbb{C}}^{LN\times 1}$, ${\mathbf{H}}=[{\mathbf{H}}_{1}^{T},\ldots,{\mathbf{H}}_{L}^{T}]^{T}\in{\mathbb{C}}^{LN\times K}$ and ${\mathbf{x}}=[x_{1},\ldots,x_{K}]^{T}{\mathbb{C}}^{K\times 1}$. Based on Bayes’ theorem, the posterior probability is given by

where $\mathtt{P}({\mathbf{y}}|{\mathbf{x}},{\mathbf{H}})$ is the likelihood function with the known channel matrix and $\mathtt{P}({\mathbf{x}})$ is the prior distribution of ${\mathbf{x}}$. Given the posterior probability $\mathtt{P}({\mathbf{x}}|{\mathbf{y}},{\mathbf{H}})$, the Bayesian MMSE estimate is obtained by

However, the Bayesian MMSE estimator is not tractable because the marginal posterior probability in (5) involves a high-dimensional integral, which motivates us to develop an advanced method to approximate (4) effectively. Furthermore, different from conventional cellular MIMO systems[31], the posterior probability in (4) has to be rewritten in a distributed way because of the inherent distributed characteristic of cell-free massive MIMO systems, given by

As a result, we can develop an efficient distributed MIMO detector to approximate the posterior probability (6) by resorting to message passing on factor graph and take the characteristic of cell-free massive MIMO into consideration. The factor graph and proposed detector will be elaborated in next sections.

Factor graph is a powerful tool to develop efficient algorithms for solving inference problems by performing different message-passing rules. We first introduce several concepts for the factor graph [37], as illustrated in Fig. 2, where the hollow circles represent the variable nodes and the solid squares represent the factor nodes. For a factor graph with factor nodes $\{f_{0},f_{1},\ldots,f_{L}\}$, all connected to a set of variable nodes $\{{\mathbf{x}}_{i}\}$, the messages are computed and updated iteratively by performing the following rules.

1) Approximate beliefs: The approximate belief $b_{\mathrm{app}}({\mathbf{x}})$ on variable node ${\mathbf{x}}$ is $\mathcal{N}_{\mathbb{C}}({\mathbf{x}},{\mathbf{x}}_{\mathrm{A}}^{\mathrm{ext}},v_{\mathrm{A}}^{\mathrm{ext}})$, where ${\mathbf{x}}_{\mathrm{A}}^{\mathrm{ext}}=\mathtt{E}[{\mathbf{x}}|b_{\mathrm{sp}}]$ and $v_{\mathrm{A}}^{\mathrm{ext}}=\mathrm{mean}(\mathrm{diag}(\mathrm{Cov}[{\mathbf{x}}|b_{\mathrm{sp}}]))$ are the mean and average variance of the corresponding beliefs $b_{\mathrm{sp}}({\mathbf{x}})\propto\prod_{l}\mu_{f_{l}\rightarrow{\mathbf{x}}}({\mathbf{x}})$.

2) Variable-to-factor messages: The message from a variable node ${\mathbf{x}}$ to a connected factor node $f_{l}$ is given by

which is the product of the message from other factors.

3) Factor-to-variable messages: The message from a factor node $f_{l}$ to a connected variable node ${\mathbf{x}}_{i}$ is $\mu_{f_{l}\rightarrow{\mathbf{x}}_{i}}({\mathbf{x}})\propto\int f({\mathbf{x}}_{i},\{{\mathbf{x}}_{j}\}_{j\neq i})\prod_{j\neq i}\mu_{{\mathbf{x}}_{j}\rightarrow f_{l}}{\mathrm{d}}{\mathbf{x}}_{j}$.

If we treat the $L$ APs and the CPU as the factor nodes, we can construct the factor graph for cell-free massive MIMO detection. By further applying the above message-passing rules to the constructed factor graph, we develop the distributed EP detector for cell-free massive MIMO systems. The detailed derivation is given in Appendix A-A.

The distributed EP-based detector is illustrated in Algorithm 1. Compared with the centralized EP detector, it deploys partial calculation modules at the AP based on the local information and then sends the estimated posterior mean and variance to the CPU for combining. In particular, the input of the algorithm is the received signal ${\mathbf{y}}_{l}$, channel matrix ${\mathbf{H}}_{l}$, and noise level $\sigma^{2}$, while the output is the recovered signal ${\mathbf{x}}_{\mathrm{B}}^{\mathrm{post},T}$ in the $T$-th iteration. The initial parameters are $\boldsymbol{\gamma}_{l}^{(0)}=\mathbf{0}$, and $\lambda_{l}^{(0)}=\frac{1}{E_{x}}$, where

Note that $\boldsymbol{\gamma}_{l}^{(0)}$ and $\lambda_{l}^{(0)}$ are the initialized extrinsic information, and $E_{x}$ is the power of the transmitted symbol $x_{k}$. The block diagram of the proposed distributed EP detector is illustrated in Fig. 3, which is composed of two modules, A and B. Each module uses the turbo principle in iterative decoding, where each module passes the extrinsic messages to the other module and this process repeats until convergence. It is derived from the factor graph with the messages updated and passed between different pairs of nodes that are assumed to follow Gaussian distributions. As the Gaussian distribution can be fully characterized by its mean and variance, only the mean and variance need to be calculated and passed between different modules. To better understand the distributed EP detection algorithm, we elaborate the details of each equation in Algorithm 1. Specifically, module A is the linear MMSE (LMMSE) estimator performed at the APs according to the following linear model

In the $t$-iteration, the explicit expressions for the posterior covariance matrix $\boldsymbol{\Sigma}_{l}^{t}$ and mean vector $\boldsymbol{\mu}_{l}^{t}$ are given by (14) and (15), respectively. Note that each AP only uses the local channel ${\mathbf{H}}_{l}$ to detect the transmitted signal ${\mathbf{x}}$. For the ease of notation, we omit the iteration index $t$ for all estimates for the mean and variance in Algorithm 1. Then, the variance estimate $v_{\mathrm{A},l}^{\mathrm{post}}=\mathrm{tr}(\boldsymbol{\Sigma}_{l}^{t})/|\mathcal{D}_{l}|$ and mean estimate ${\mathbf{x}}_{\mathrm{A},l}^{\mathrm{post}}=\boldsymbol{\mu}_{l}^{t}$ are transferred to the CPU to compute the extrinsic information $v_{\mathrm{A},l}^{\mathrm{ext}}$ (16) and ${\mathbf{x}}_{\mathrm{A},l}^{\mathrm{ext}}$ (17), respectively. The extrinsic information ${\mathbf{x}}_{\mathrm{A},l}^{\mathrm{\mathrm{ext}}}$ can be regarded as the AWGN observation given by

where ${\mathbf{n}}^{\mathrm{ext}}_{\mathrm{A},l}\sim\mathcal{N}_{{\mathbb{C}}}(0,v_{\mathrm{A},l}^{\mathrm{ext}}{\mathbf{I}})$ [33]. As a result, the linear model in (3) is decoupled into $K$ parallel and independent AWGN channels with equivalent noise $v_{\mathrm{A},l}^{\mathrm{ext}}$ for each AP. Subsequently, the CPU collects all extrinsic means $\{{\mathbf{x}}_{\mathrm{A},l}^{\mathrm{ext}}\}_{l=1}^{L}$ and variances $\{v_{\mathrm{A},l}^{\mathrm{ext}}\}_{l=1}^{L}$ and performs MRC¹¹1In our paper, we focus on linear combining of the local estimate $\mathbf{x}_{\mathrm{A},l}^{\mathrm{ext}}$ and $v_{\mathrm{A},l}^{\mathrm{ext}}$ from APs as $\mathbf{x}_{\mathrm{A}}^{\mathrm{ext}}=\sum_{l=1}^{L}a_{l}\mathbf{x}_{\mathrm{A},l}^{\mathrm{ext}},\quad l=1,\ldots,L$. The optimal combining is MRC, which can maximize the post-combination SNR of the final AWGN observation $\mathrm{x}_{\mathrm{A}}^{\mathrm{ext}}$ with the linear combination of $\mathbf{x}_{\mathrm{A},l}^{\mathrm{ext}}$ where $\mathbf{x}_{\mathrm{A}}^{\mathrm{ext}}=\mathbf{x}+\mathbf{n}^{\mathrm{eq}}$ and $\mathbf{n}_{\mathrm{A}}^{\mathrm{ext}}\sim\mathcal{N}_{\mathbb{C}}(0,v_{\mathrm{A}}^{\mathrm{ext}}\mathbf{I})$. The final expressions are shown in (18) and (19) and the detailed derivation is given in Appendix A-B.. The MRC expressions (18) and (19) are obtained by maximizing the post-combination signal-to-noise ratio (SNR) of the AWGN observation ${\mathbf{x}}_{\mathrm{A}}^{\mathrm{ext}}$ at the CPU, given by

where ${\mathbf{n}}^{\mathrm{ext}}_{\mathrm{A}}\sim\mathcal{N}_{{\mathbb{C}}}(0,v_{\mathrm{A}}^{\mathrm{ext}}{\mathbf{I}})$. The CPU uses the posterior mean estimator to detect the signal ${\mathbf{x}}$ from the equivalent AWGN model (11). Then, the posterior mean and variance are computed by the posterior MMSE estimator for the equivalent AWGN model in (20) and (21). As the transmitted symbol is assumed to be drawn from the $M$-QAM set $\mathcal{S}=\{s_{1},s_{2},\ldots,s_{M}\}$, the corresponding expressions for each element in (20) and (21) are given by

where $[x_{\mathrm{B}}^{\mathrm{post}}]_{k}$, $[v_{\mathrm{B}}^{\mathrm{post}}]_{k}$, and $[x_{\mathrm{A}}^{\mathrm{ext}}]_{k}$ are the $k$-th element in ${\mathbf{x}}_{\mathrm{B}}^{\mathrm{post}}$, ${\mathbf{v}}_{\mathrm{B}}^{\mathrm{post}}$, and ${\mathbf{x}}_{\mathrm{A}}^{\mathrm{ext}}$, respectively. The posterior mean ${\mathbf{x}}_{\mathrm{B}}^{\mathrm{post}}$ and variance ${\mathbf{v}}_{\mathrm{B}}^{\mathrm{post}}$ are then utilized to compute the extrinsic information $\lambda_{l}^{(t)}$ and $\boldsymbol{\gamma}_{l}^{(t)}$ for each AP in (22) and (23), where the function $\mathrm{mean}(\cdot)$ is used to compute the mean of a vector. Finally, the extrinsic information $\lambda_{l}^{(t)}$ and $\boldsymbol{\gamma}_{l}^{(t)}$ are transferred to each AP in the next iteration. The whole procedure is executed iteratively until it is terminated by a stopping criterion or a maximum number of iterations.

In the following, we provide the complexity analysis and fronthaul overhead for different detectors in Tables I and II, respectively. For the proposed distributed EP detector, the computational complexity at each AP is dominated by the LMMSE estimator for estimating the signal ${\mathbf{x}}_{\mathrm{A},l}^{\mathrm{post}}$, which is $\mathcal{O}(N|\mathcal{D}_{l}|^{2})$ because of the matrix inversion required in (14), while the computational complexity at the CPU is $O(L|\mathcal{D}_{l}|)$ in each iteration. Furthermore, if the number of antennas $N$ is less than $|\mathcal{D}_{l}|$, we can use the matrix inversion lemma to carry out the matrix inversion in (14) as follows

where ${\mathbf{D}}=\lambda_{l}^{(t-1)}{\mathbf{I}}$ and the computational complexity is reduced to $\mathcal{O}(|\mathcal{D}_{l}|N^{2})$. Therefore, the overall computational complexity at each AP is $\mathcal{O}(T|\mathcal{D}_{l}|N^{2})$ while the overall computational complexity at the CPU is $\mathcal{O}(TL|\mathcal{D}_{l}|)$ for $T$ iterations. As observed in Table I, the distritbuted EP detector mainly increases the computational complexity at the CPU when compared with other distributed linear receivers.

We further compare the number of complex scalars that need to be transmitted from the APs to the CPU via the fronthauls of different detectors in Table II. We assume that $\tau_{c}$ and $\tau_{p}$ are the coherence time and pilot length, respectively. The results for the 4 baseline detectors are from [27]. With the proposed detector, $\sum\limits_{l=1}^{L}(\tau_{c}-\tau_{p})2T(|\mathcal{D}_{l}|+1)$ scalars need to be passed from the APs to the CPU and no statistical parameters are required to be passed, where $T$ denotes the total number of iterations. The fronthaul overhead of the proposed distributed EP detector is similar to those of the centralized MMSE and EP detectors. The detailed comparison shall be determined by the value of system parameters.

ICD has been explored for MIMO [38, 39], OFDM [40], and massive MIMO [41] systems with a low-resolution analog-to-digital converter (ADC) [42]. Recently, it has been investigated in cell-free massive MIMO to improve the BER performance and reduce the pilot overhead. Different from [43] that solved the complex biconvex optimization problem with high complexity, we develop a low-complexity turbo-like ICD architecture. The iterative procedure is summarized in Algorithm 2.

As illustrated in Fig. 3, the proposed ICD scheme employs a similar idea as iterative decoding. At each AP, the channel estimator and Module A in the proposed distributed EP detector are performed. They exchange extrinsic information iteratively. In the ICD processing stage, the channel estimation is first performed based on the received pilot signal and transmitted pilot. Then, the data detector performs signal detection by taking both the channel estimation error and channel statistics into consideration. It will then feed back the detected data and detection error to the channel estimator. Finally, data-aided channel estimation is employed with the help of the detected data. The whole ICD procedure is executed iteratively until it is terminated by a stopping criterion or a maximum number of ICD iterations.

Similar to Algorithm 1, the input of the ICD scheme is the pilot signal matrix, $\mathbf{X}^{{\mathrm{p}}}$, the received signal matrix corresponding to the pilot matrix, ${\mathbf{Y}}^{{\mathrm{p}}}$, and that corresponding to the data matrix, ${\mathbf{Y}}^{{\mathrm{d}}}$, in each time slot for each AP. In the $r$-th ICD iteration²²2 Here we use $r$ to denote the number of data feedback from the data detector to the channel estimator. For example, $r=1$ means no data feedback while $r=4$ means 3 data feedback should be performed., $\hat{{\mathbf{H}}}^{(r)}$ is the estimated channel matrix, $\hat{{\mathbf{X}}}^{({\mathrm{d}},r)}$ is the estimated data matrix, and $\hat{{\mathbf{V}}}_{\mathrm{est},l}$ and $\hat{{\mathbf{V}}}_{\mathrm{det},l}$ are used to compute the covariance matrix for the equivalent noise in the signal detector and channel estimator, respectively. The final output of the signal detector is the detected data matrix $\hat{{\mathbf{X}}}^{(\mathrm{d},R)}$, where $R$ is the total number of ICD iterations. Compared with the conventional receiver where the channel estimator and signal detector are designed separately, the ICD scheme can improve the BER performance by considering the characteristics of the channel estimation error and channel statistics. After performing signal detection, the detected payload data will be utilized for channel estimation. Next, we will elaborate the whole procedure in detail.

We consider adopting the classical LMMSE estimator in the channel estimation stage at each AP. To facilitate the representation of the channel estimation problem, we consider the matrix vectorization to (9) and rewrite it as

where ${\mathbf{A}}^{{\mathrm{p}}}=({\mathbf{X}}^{{\mathrm{p}}})^{T}\otimes{\mathbf{I}}_{N}\in\mathbb{C}^{\tau_{p}N\times KN}$, ${\mathbf{y}}^{{\mathrm{p}}}_{l}=\mathrm{vec}({\mathbf{Y}}^{{\mathrm{p}}}_{l})\in\mathbb{C}^{\tau_{p}N\times 1}$, $\bar{{\mathbf{h}}}_{l}=\mathrm{vec}({\mathbf{H}}_{l})\in\mathbb{C}^{NK\times 1}$, and ${\mathbf{n}}^{{\mathrm{p}}}_{l}=\mathrm{vec}(\mathbf{N}^{{\mathrm{p}}}_{l})\in\mathbb{C}^{\tau_{p}N\times 1}$. We denote $\otimes$ as the matrix Kronecker product and $\mathrm{vec}(\cdot)$ as the vectorization operation. In the pilot-only based channel estimation stage, the LMMSE estimate of $\bar{{\mathbf{h}}}_{l}$ is given by

where ${\mathbf{R}}_{\bar{{\mathbf{h}}}_{l}\bar{{\mathbf{h}}}_{l}}$ is the channel covariance matrix and depends on the spatial correlation matrix ${\mathbf{R}}_{kl}$. Based on the property of the LMMSE estimator, $\hat{\bar{{\mathbf{h}}}}_{l}^{{\mathrm{p}}}$ is a Gaussian random vector, and the channel estimation error vector $\Delta{{\mathbf{h}}^{{\mathrm{p}}}_{l}}=\hat{\bar{{\mathbf{h}}}}^{{\mathrm{p}}}_{l}-{\mathbf{h}}_{l}$ is also a Gaussian random vector with zero-mean and covariance matrix ${\mathbf{R}}_{\Delta{{\mathbf{h}}_{l}^{\mathrm{p}}}}$ given by

After performing channel estimation, the estimated channel is used to perform data detection at each AP. In the data transmission stage, the received data signal vector ${\mathbf{y}}^{{\mathrm{d}}}[n]$ corresponding to the $n$-th data in each coherence time can be expressed by ${\mathbf{y}}^{{\mathrm{d}}}[n]={\mathbf{H}}{\mathbf{x}}^{{\mathrm{d}}}[n]+{\mathbf{n}}^{{\mathrm{d}}}[n]$, where ${\mathbf{n}}^{{\mathrm{d}}}[n]\sim\mathcal{N}_{\mathbb{C}}(0,\sigma^{2}{\mathbf{I}}_{N})$ is the AWGN vector³³3As each AP has similar signal models, we remove the index $l$ for the AP and use ${\mathbf{y}}^{{\mathrm{d}}}[n]$, ${\mathbf{H}}$, ${\mathbf{x}}^{{\mathrm{d}}}[n]$, and ${\mathbf{n}}^{{\mathrm{d}}}[n]$ to denote ${\mathbf{y}}_{l}^{{\mathrm{d}}}[n]$, ${\mathbf{H}}_{l}$, ${\mathbf{x}}_{l}^{{\mathrm{d}}}[n]$, and ${\mathbf{n}}_{l}^{{\mathrm{d}}}[n]$, respectively. To simplify the notation, we have the same operation for the following symbols.. We can rewrite it in a matrix from as ${\mathbf{Y}}^{{\mathrm{d}}}={\mathbf{H}}{\mathbf{X}}^{{\mathrm{d}}}+{\mathbf{N}}^{{\mathrm{d}}}$, where ${\mathbf{N}}^{{\mathrm{d}}}=\left[{\mathbf{n}}^{{\mathrm{d}}}[1],\ldots,{\mathbf{n}}^{{\mathrm{d}}}[\tau_{d}]\right]\in{\mathbb{C}}^{N\times\tau_{d}}$ is the AWGN matrix in the data transmission stage. Denote the estimated channel $\hat{{\mathbf{H}}}\in\mathbb{C}^{N\times K}$ as $\hat{{\mathbf{H}}}={\mathbf{H}}+\Delta{\mathbf{H}}$, where $\Delta{\mathbf{H}}$ is the channel estimation error. The signal detection problem can be formulated as

where $\hat{{\mathbf{n}}}^{{\mathrm{d}}}[n]={\mathbf{n}}^{{\mathrm{d}}}[n]-\Delta{\mathbf{H}}{\mathbf{x}}^{{\mathrm{d}}}[n]$ is the equivalent noise in the signal detector and the statistical characterization of $\hat{{\mathbf{n}}}^{{\mathrm{d}}}[n]$ is given in Appendix B-A.

To further improve the BER performance, a data-aided channel estimation approach is adopted in the channel estimation stage. First, conventional pilot-only based channel estimation is performed and the transmitted symbols are detected. Then, the detected symbols are fed back to the channel estimator as additional pilot symbols to refine the channel estimation. In the channel training stage, pilot matrix ${\mathbf{X}}^{{\mathrm{p}}}$ is transmitted. The received signal matrix corresponding to the pilot matrix ${\mathbf{Y}}^{{\mathrm{p}}}\in\mathbb{C}^{N\times\tau_{p}}$ is expressed as ${\mathbf{Y}}^{{\mathrm{p}}}={\mathbf{H}}{\mathbf{X}}^{{\mathrm{p}}}+{\mathbf{N}}^{{\mathrm{p}}}$, where ${\mathbf{N}}^{{\mathrm{p}}}=\left[{\mathbf{n}}^{{\mathrm{p}}}[1],\ldots,\mathbf{n}^{{\mathrm{p}}}[\tau_{p}]\right]\in\mathbb{C}^{N\times\tau_{p}}$ is the AWGN matrix and each column $\mathbf{n}^{{\mathrm{p}}}[n]\sim\mathcal{N}_{\mathbb{C}}(0,\sigma^{2}\mathbf{I}_{N})$ for $n=1,\ldots,\tau_{p}$. The estimated data matrix $\hat{\mathbf{X}}^{{\mathrm{d}}}$ can be expressed as $\hat{{\mathbf{X}}}^{{\mathrm{d}}}={\mathbf{X}}^{{\mathrm{d}}}+{\mathbf{E}}^{{\mathrm{d}}}$, where $\mathbf{E}^{{\mathrm{d}}}$ is the signal detection error matrix. In the data-aided channel estimation stage, the estimated $\hat{{\mathbf{X}}}^{{\mathrm{d}}}$ are fed back to the channel estimator as additional pilots. Then, the received signal matrix ${\mathbf{Y}}^{{\mathrm{d}}}$ corresponding to $\hat{{\mathbf{X}}}^{{\mathrm{d}}}$ can be expressed as

where $\hat{{\mathbf{N}}}^{{\mathrm{d}}}={\mathbf{N}}^{{\mathrm{d}}}-{\mathbf{H}}{\mathbf{E}}^{{\mathrm{d}}}$ is the equivalent noise for the data-aided pilot $\hat{{\mathbf{X}}}^{{\mathrm{d}}}$. The statistical information of the $n$-th column of $\hat{{\mathbf{N}}}^{{\mathrm{d}}}$ is $\hat{\mathbf{n}}^{{\mathrm{d}}}[n]\sim\mathcal{N}_{{\mathbb{C}}}(\mathbf{0},\hat{\mathbf{V}}^{{\mathrm{d}}}[n])$ for $n=1,\ldots,\tau_{d}$, where $\hat{\mathbf{V}}^{{\mathrm{d}}}[n]$ will be utilized for data-aided channel estimation. We denote ${\mathbf{Y}}=[{\mathbf{Y}}^{{\mathrm{p}}},{\mathbf{Y}}^{{\mathrm{d}}}]$ as the received signal matrix corresponding to the overall transmitted signal. By concatenating the received pilot and data signals, we have

where ${\mathbf{X}}=[{\mathbf{X}}^{{\mathrm{p}}},\hat{{\mathbf{X}}}^{{\mathrm{d}}}]$ and ${\mathbf{N}}=[{\mathbf{N}}^{{\mathrm{p}}},\hat{{\mathbf{N}}}^{{\mathrm{d}}}]$ can be interpreted as the equivalent pilot signal and noise in the data-aided channel estimation stage, respectively. The statistical characteristic for the equivalent noise and corresponding LMMSE estimator are derived in Appendix B-B.

We first provide a state evolution analysis framework to predict the asymptotic performance of the distributed EP detector in the large-system limit⁴⁴4The large-system limit means that the cell-free massive MIMO system with $L,K\rightarrow\infty$..

The function $\mathrm{MSE}(\cdot)$ is given by

where the expectation is with respect to $x$. $\tau_{k,l}$ are the eigenvalues of ${\mathbf{H}}^{H}_{l}{\mathbf{H}}_{l}$.

The state equations illustrated in Proposition 1 can be proved rigorously in a similar way as [44]. To better understand the state evolution analysis, we give an intuitive derivation here. In Algorithm 1, we have explicit expressions for computing the mean and variance for each module. By substituting (14) into (16), we have the following expression

The asymptotic expressions for $v_{\mathrm{A}}^{\mathrm{ext},t}$ and $\lambda_{l}^{(t)}$ can then be derived from (18) and (22). Finally, the asymptotic MSE can be interpreted as the MSE of the decoupled scalar AWGN channels (11) and is related to the distribution of the transmitted signal ${\mathbf{x}}$. Next, we will give a specific example for Proposition 1.

Example 1: If the channel matrix ${\mathbf{H}}_{l}\,(l=1,2,\ldots,L)$ is an i.i.d. Gaussian distributed matrix, then $v_{\mathrm{A},l}^{\mathrm{ext}}$ will converge to the following asymptotic expression by adopting the result of the $\mathcal{R}$-transform⁵⁵5The concept of $\mathcal{R}$-transform can be found in [45], which enables the characterization of the asymptotic spectrum of a sum of suitable matrices from their individual asymptotic spectra. of the average empirical eigenvalue distribution given by,

Furthermore, if the data symbol is drawn from a quadrature phase-shift keying (QPSK) constellation, the $\mathrm{MSE}$ is expressed by

Furthermore, the BER w.r.t. ${\mathbf{x}}$ can be evaluated through the equivalent AWGN channel (11) with an equivalent $\mathrm{SNR}=1/v_{\mathrm{A},l}^{\mathrm{ext},T}$ and is given by

where $Q(x)=\int_{x}^{\infty}\mathrm{D}z$ is the $Q$-function. In fact, the MSE and BER are determined given the knowledge of the AWGN channel (11) with $\mathrm{SNR}=1/v_{\mathrm{A},l}^{\mathrm{ext},T}$, which is known as the decoupling principle. Thus, if the data symbol is drawn from other $M$-QAM constellations, the corresponding BER can be easily obtained using a closed-form BER expression [46].

The fixed points of the distributed EP detection algorithm are the stationary points of a relaxed version of the Kullback-Leibler (KL) divergence to minimize the posterior probability (4), which is given by

The KL divergence between two distributions is defined as,

Generally, it is difficult to approximate the complex posterior distribution (4) by a tractable approach. The moment matching method is utilized to exploit the relaxations of the KL divergence minimization problem [47]. The process for solving the KL divergence minimization problem is equivalent to optimizing the Bethe free energy subject to the moment-matching constraint. That is,

where $J\left(b,b_{1},\ldots,b_{L},q\right)=\mathrm{KL}\left(b\|e^{\log f({\mathbf{x}})}\right)+\sum\limits_{l=1}^{L}\mathrm{KL}\left(b_{l}\|e^{\log f_{L}({\mathbf{x}})}\right)+H(q)$ and $H(\cdot)$ is the differential entropy. In addition, $b({\mathbf{x}})$, $\{b_{c}({\mathbf{x}})\}_{l=1}^{L}$, and $q({\mathbf{x}})$ are given by

respectively, where $Z_{0}\left({\mathbf{x}}_{\mathrm{A}}^{\mathrm{ext}}\right)$, $\{Z_{c}\left({\mathbf{x}}_{\mathrm{B},l}^{\mathrm{ext}}\right)\}_{c=1}^{C}$, and $Z_{q}\left(\hat{{\mathbf{x}}}\right)$ are the normalization factors corresponding to their density functions. Based on the above-mentioned optimization process, we have

for any fixed points $v$ and $\hat{{\mathbf{x}}}$ in Algorithm 1. This fixed point characterization is helpful to perform the convergence analysis and understand the algorithm from the optimization perspective.

In this subsection, we consider the BER performance of the proposed detector in the conventional cell-free massive MIMO where each AP serves all users. We assume that perfect CSI can be obtained at each AP. The simulation parameters are set as $N=8$, $K=32$, and $L=8$. Fig. 4 compares the achievable BER of the proposed distributed EP, centralized EP detectors [31], and other baselines [27] using QPSK and $16$-QAM modulation symbols. The results are obtained by Monte Carlo simulations with $10,000$ independent channel realizations. We denote “deEP” as the distributed EP detector proposed in Section III. It can be observed that the proposed distributed EP detector outperforms the distributed MMSE detector with only one EP iteration. This is because deEP incorporates the prior information of the symbol ${\mathbf{x}}$ into the posterior mean estimator for the equivalent AWGN channel (10). Furthermore, it outperforms the centralized MMSE detector with $T=5$ iterations for different modulation symbols. Compared with the centralized EP detectors [31], the distributed EP detector achieves comparable BER performance.

We next verify the accuracy of the analytical framework in Fig. 5 with different system settings and modulation schemes. As shown in the figure, the BERs of the proposed detector match well with the derived analytical results, which demonstrates the accuracy of the analytical framework developed in Section V-A. Although the analytical framework is derived in the large-system limit, it still can predict the BER performance in small-sized MIMO systems illustrated in Fig. 5(b). Therefore, instead of performing time-consuming Monte Carlo simulations to obtain the corresponding performance metrics, we can predict the theoretical behavior by the state evolution equations. Furthermore, the analytical framework can be further utilized to optimize the system design.