Massive multiple-input multiple-output (mMIMO) is a key technology for improving spectral efficiency and interference management in fifth-generation (5G) and beyond wireless networks. By exploiting spatial diversity, mMIMO supports multiple users within the same frequency-time block [1]- [5]. However, its performance may degrade under certain conditions: i.) in highly correlated channels, where rank deficiency S. Nassirpour, T. V. Nguyen, T. Ratnarajah, and D. H. N. Nguyen are with the Department of Electrical and Computer Engineering, San Diego State University, San Diego, CA 92182, USA. Emails: snassirpour@sdsu.edu, tnguyen58@sdsu.edu, tratnarajah@sdsu.edu, and duy.nguyen@sdsu.edu.
H. Q. Ngo is with the School of Electronics, Electrical Engineering and Computer Science, Queen's University, Belfast, UK. Email: hien.ngo@qub.ac.uk.
L. N. Tran is with the School of Electrical and Electronic Engineering, University College Dublin, Dublin, Ireland. Email: nam.tran@ucd.ie. limits coverage; ii.) with widely distributed users, where large user-base station (BS) distances cause significant path loss; and iii.) at cell edges, where users experience strong inter-cell interference from neighboring mMIMO BSs.
To tackle the above challenges, cell-free mMIMO (CF-mMIMO) has attracted growing interest from both academia and industry [6]- [11]. Unlike traditional mMIMO, which uses a large number of collocated antenna elements at a single BS, CF-mMIMO employs many distributed access points (APs). Each AP may be equipped with a single or multiple antenna elements, and they are positioned far apart from one another. This distributed architecture reduces channel correlation and lowers average path loss between users and APs compared to the collocated antenna setup in traditional mMIMO.
Moreover, CF-mMIMO adopts a user-centric design, unlike the cell-centric approach of traditional mMIMO where a single BS serves all users in a cell. In CF-mMIMO, distributed APs coordinate with a central processing unit (CPU) via fronthaul links, enhancing connectivity and overall system performance.
Initial CF-mMIMO studies assumed fully centralized processing at the CPU [9], [10], showing notable gains in median and 95%-likely spectral efficiency over traditional small-cell mMIMO, where each BS antenna serves its own users. However, this approach demands high fronthaul bandwidth, as each AP must employ a high-resolution quantizer (e.g., 10+ bits) to generate nearly continuous signals for transmission to the CPU. This becomes increasingly unsustainable as the number of users in 5G and beyond networks continues to grow. To mitigate this, the authors of [11] proposed local processing at APs, enabling partial to fully decentralized architectures. Their approach used a linear filter (i.e., the linear minimum mean squared error (LMMSE) filter) to enhance spectral efficiency under limited fronthaul bandwidth.
While linear filters offer low computational complexity, their performance declines when the number of AP antennas is small relative to the number of users or when channels are highly correlated. To address these limitations, nonlinear methods such as approximate message passing (AMP), generalized AMP (GAMP), expectation propagation (EP), and variational Bayesian (VB) inference have been proposed [12]- [16]. AMP is effective for data detection in mMIMO under independent and identically distributed (i.i.d.) Rayleigh fading by decoupling the system into parallel additive white Gaussian noise (AWGN) channels [12], [13]. However, it may diverge under ill-conditioned or non-zero-mean channels (e.g., correlated or Rician fading) [14]. Moreover, GAMP extends AMP to support non-Gaussian priors and quantized measurements [15]; however, it suffers instability with noni.i.d. measurement matrices or small system sizes. On the other hand, EP minimizes Kullback-Leibler (KL) divergence via message passing and moment matching [16], but it lacks convergence guarantees and may diverge in high-dimensional settings.
The performance of CF-mMIMO networks strongly depends on accurate channel estimation (CE), motivating several studies on CE techniques [17], [18]. For instance, [17] proposed a subspace-based semi-blind CE scheme to mitigate pilot contamination, which arises when the number of users exceeds the available pilot duration. To further address this issue, [19] developed a graph coloring-based pilot assignment method that models user interference through AP selection, improving pilot reuse and reducing contamination. Beyond traditional CE, recent works have explored joint channel estimation and data detection (JED), where data detected during the data transmission phase is used to refine CE performance. In [20], a generalized superimposed pilot scheme was introduced for JED, distributing data symbols across coherence time slots via linear precoders. Then, [21] proposed an iterative JED algorithm to solve a relaxed maximum a posteriori (MAP) problem using the forward-backward splitting technique.
Motivated by the above, in this paper, we focus on the JED problem in the uplink of CF-mMIMO networks and propose an approach based on VB inference. The VB method offers a robust framework for approximating posterior distributions by optimizing a simpler, tractable distribution to closely match the intractable true posterior. Originally developed for machine learning, VB techniques have recently gained attraction in wireless communications [14], [22]- [24]. Unlike other nonlinear methods like AMP, GAMP, and EP, VB is well-suited for ill-conditioned channels (e.g., correlated channels) and has solid convergence properties [25], [26].
It is worth noting that recent studies in [27]- [29] have explored joint user activity detection and channel estimation in CF-mMIMO networks, formulating the problem as a compressed sensing (CS) task by leveraging user sparsity. Specifically, [27] employed a single measurement vector (SMV)based minimum mean squared error (MMSE) approach, while [28] used the GAMP algorithm to tackle the CS problem. In [29], the authors extended the problem to include data detection and proposed a grant-free scheme based on bilinear inference, using the bilinear Gaussian belief propagation (Bi-GaBP) algorithm. Their method achieved efficient joint estimation without data spreading and mitigated pilot contamination through a low-coherence pilot design.
To address the limited bandwidth of fronthaul links, in addition to the previously discussed approach of local processing at the APs, two alternative methods can be considered. The first, called the quantization-and-estimation (Q-E) scenario, uses low-bit quantizers at each AP to send quantized received signals to the CPU. The second, referred to as the estimationand-quantization (E-Q) scenario, performs local CE at each AP during the pilot transmission phase and forwards quantized versions of the local CE and received signals during the data transmission phase to the CPU. In both cases, final JED is performed at the CPU. Unlike the perfect fronthaul link (PFL) scenario with unquantized received signals, the CPU in Q-E and E-Q relies on quantized inputs, introducing challenges due to the nonlinearity of the quantization function and the resulting quantization noise. In [30], the authors studied spectral and energy efficiency under Q-E and E-Q, using uniform quantization modeled by the Max algorithm and Bussgang decomposition.
Next, the work in [31] focused on joint user activity detection and channel estimation in a CF-mMIMO network under Q-E with mixed-quantization APs, where some APs receive unquantized signals while others use low-bit quantizers. This problem was modeled as a CS task, and the authors proposed an approach based on multiple measurement vectors and GAMP to tackle it. Then, the study in [32] addressed the JED problem in the CF-mMIMO network under Q-E. The authors proposed a two-stage solution: first, a de-quantization step based on scalar Gaussian approximation and Bussgang decomposition to estimate the statistics of the unquantized signals; second, they applied a BiGaBP algorithm for JED. Nevertheless, the Bussgang decomposition discussed in [31], [32], which linearizes the relationship between quantized and unquantized signals, has two key limitations: it requires access to signal statistics before and after quantization, and it is applicable only when the input to the quantizer follows a Gaussian distribution. These constraints limit its applicability and leave room for improvements beyond Bussgang-based methods [33]. Building on this idea, [34] studied the JED problem in mMIMO systems with quantized signals. The authors used a variant of belief propagation (BP) within a GAMP-based framework to approximate the marginal distributions of data and channel components, and derived analytical expressions for quantized observations using the Gaussian cumulative distribution function (CDF). However, this method may diverge under correlated channels, a known limitation of AMP-type algorithms. To address this, [35] proposed a VB-based JED framework for mMIMO systems and used the expectation maximization (EM) technique to estimate the precision terms. The VB framework guarantees convergence to at least a local optimum [25], [26].
Motivated by the results in [35], in this work, we propose a VB-based approach to model the nonlinear relationship between unquantized and quantized signals in both Q-E and E-Q scenarios. Furthermore, we adopt Gamma priors for the precision parameters, generalizing the VB framework and showing that the VB-EM approach is a special case. We evaluate the performance of our proposed VB-based method in terms of symbol error rate (SER), normalized mean squared error (NMSE) for channel estimation, computational complexity, and fronthaul signaling overhead.
Contribution: Our main contributions are as follows:
• We consider uplink communications in a CF-mMIMO network and propose a method based on VB inference to address the JED problem in both Q-E and E-Q scenarios. • Unlike previous works in [31], [32] that utilized Bussgang decomposition to linearize the relationship between quantized and unquantized signals, we leverage the VB framework to approximate the posterior distribution of the unquantized signal given its quantized counterpart. • We model uncertainty in the JED process via residual inter-user interference, assumed to be Gaussian, capturing both noise and detection/estimation errors. Its precision is also treated as a random variable and estimated within the VB framework. • We assess the performance of our proposed VB-based approach under the PFL, Q-E, and E-Q scenarios in terms of SER, channel NMSE, computational complexity, and fronthaul signaling overhead. We compare its performance with LMMSE(PFL), GAMP(PFL), GAMP(Q-E) [34], VB-EM(PFL), and VB-EM(Q-E) [35], where VB-EM is built upon the VB framework and estimates the precision of the residual inter-user interference using EM. Numerical results show that our VB-based methods consistently outperform LMMSE(PFL) by capturing nonlinear relationships, and also surpass the GAMP-and VB-EM-based methods, thanks to a unified, fully VB-driven design that ensures convergence even under correlated channels. In contrast, GAMP may diverge, and EM-based precision estimation within VB-EM can degrade performance. Finally, VB(Q-E) slightly outperforms VB(E-Q) due to local CE errors in the latter, aligning with the findings in [30].
The rest of this paper is structured as follows. Section II introduces the system model. Section III details the proposed VB framework for the JED problem. Section IV presents the simulation results, and Section V provides the conclusion.
Notation: In this paper, scalars are represented by italic letters, vectors by bold lowercase letters, and matrices by bold uppercase letters. The notation CN (m, C) denotes a complex Gaussian random vector with mean m and covariance C. Γ(a, b) implies a Gamma distribution parametrized by a and b. The space of x × y complex-valued matrices is denoted as C x×y . We use diag{v} to represent a diagonal matrix formed from the vector v. Further, we indicate the determinant, transpose, conjugate transpose, Euclidean norm, and Frobenius norm of matrix A by |A|, A ⊤ , A H , ∥A∥, and ∥A∥ F , respectively. The symbols ∼ and ∝ represent "distributed according to" and "proportional to," respectively. The element in the i th row and j th column of matrix X is denoted as [X] ij . We use [X] :m to indicate the m th column of X. Additional notations include Tr(•) for the trace function, exp{•} for the exponential function, and sign(•) for the signum function. The identity matrix of size M is represented by I M . The complex conjugate of x is indicated as x * , and the real and imaginary parts of x are denoted by ℜ{x} and ℑ{x}, respectively. The probability density function (PDF) of a length-K complex-valued random vector x ∼ CN (m, C) is given by: CN (x; m, C)
The expected value and variance of x under the distribution p(x) are denoted by E p(x) [x] and Var p(x) [x], respectively. We use ⟨x⟩, ⟨|x| 2 ⟩, and τ x to represent the mean, second moment, and variance of x under a variational distribution q(x). We denote the indicator function by 1(•), which takes the value of one if the given condition is true and zero otherwise.
In this study, we consider an uplink scenario in a CF-mMIMO network, depicted in Fig. 1, where a CPU utilizes L distributed APs to concurrently serve K users. Each AP is equipped with M antennas, while each user has a single antenna. Our focus is on addressing a JED problem within this network. To this end, we develop a method leveraging VB inference to approximate the true posterior distributions. This section begins by describing the channel model, provides a brief introduction to VB inference, and then presents the problem formulation.
We use h i,ℓ = [h [1] i,ℓ , h [2] i,ℓ , . . . , h
to denote the communication channel between the i th user and the ℓ th AP, where h
represents the channel between the m th antenna element of the ℓ th AP and the i th user. In this paper, we consider 1 ≤ K ≤ M L and assume communication over a coherence interval during which the channels remain constant.
To characterize h i,ℓ , we assume that h i,ℓ ∼ CN (0, Σ i,ℓ ), where Σ i,ℓ is the channel covariance matrix, given by:
where β i,ℓ denotes the large-scale fading between the i th user and the ℓ th AP, and Σℓ is the spatial correlation between the antenna elements at the ℓ th AP. We assume that the channel covariance matrix Σ i,ℓ is known, which is a common assumption in the literature [32], [35]. We consider that the channels are independent across users, which leads to E[h i,ℓ h H j,ℓ ] = 0 for 1 ≤ ℓ ≤ L, 1 ≤ i, j ≤ K, and i ̸ = j.
Received signal: In the considered CF-mMIMO network, the received signal at duration t at the m th antenna element of the ℓ th AP is denoted by r m,ℓ,t and is modeled as:
where x i,t represents the transmitted signal from the i th user at duration t, while n m,ℓ,t ∼ CN (0, N 0 ) is the i.i.d. AWGN at the m th antenna element of the ℓ th AP at duration t. We then organize the received signals, channels, transmitted signals, and noise elements into vectors to formulate a linear model for the received signals at the ℓ th AP at duration t as below:
where
In this work, we examine an uplink CF-mMIMO network with the goal of jointly estimating
To do this, we propose an approach relying on VB inference. Prior to exploring the methodology in-depth, we provide a concise introduction to the core principles of VB.
Variational Bayesian (VB) inference is a powerful statistical approach rooted in machine learning, designed to approximate the posterior distribution of latent variables efficiently. Consider the set of observed variables r and the set of V latent variables x. To detect x, it is necessary to evaluate the posterior distribution p(x|r), which is often computationally intractable. To address this, VB approximates p(x|r) using a distribution q(x) parameterized by variational variables, chosen from a predefined family Q. The goal is to ensure that q(x) is as close as possible to p(x|r). This is achieved by formulating an optimization problem that minimizes the KL divergence from q(x) to p(x|r):
where q ⋆ (x) represents the optimal variational approximation, and the KL divergence is given by:
The KL divergence reaches its minimum when q(x) matches the true posterior p(x|r). However, since deriving the exact posterior is usually infeasible, a practical approach involves restricting q(x) to a simplified family of distributions. A common choice in the literature [25] is the mean-field variational family, where the latent variables are assumed to be independent, leading to q(x) = V i=1 q i (x i ). Within this framework, the optimal solution for q i (x i ) is given by [31]
where ⟨•⟩ denotes the expectation with respect to all latent variables except x i , using the current variational densities q -i (x -i ) = V j=1,j̸ =i q j (x j ). To solve the optimization problem in (4), the Coordinate Ascent Variational Inference (CAVI) algorithm is commonly employed. This iterative method sequentially updates each q ⋆ i (x i ), as shown in (6), while keeping the others fixed, thereby ensuring a monotonic improvement in the objective function in (4). The CAVI algorithm is guaranteed to converge to a local optimum [25], [26].
Taking into account the received signals at all APs results in the following system model at duration t:
where
. In this study, we partition each coherence interval, consisting of T symbols, into two phases: a pilot transmission phase of length T p symbols and a data transmission phase of length T d symbols, satisfying T p + T d = T . Here, we use R p and R d to denote the received signal matrices during the pilot and data transmission phases, respectively, which are given by:
where
In this paper, we use ζ p and ζ d to denote the residual interuser interference during the pilot and data transmission phases, respectively. These terms capture uncertainty in the JED process, which includes noise, CE error, and data detection error, and are defined as follows:
We represent the estimated values of H and X d by Ĥ = H + e h and Xd = X d + e x , respectively, where e h is the CE error, and e x denotes the data detection error. By substituting these estimated values into (10) and (11), we obtain:
To effectively capture the uncertainty, we model ζ p and ζ d as i.i.d. zero-mean Gaussian random variables. We then define γ p and γ d = [γ d,Tp+1 , γ d,Tp+2 , . . . , γ d,T ] ⊤ to represent the precision of the combined effect of noise and CE error during the pilot transmission phase, and the precision of the combined effect of noise, CE error, and data detection error during the data transmission phase, respectively. Since the underlying errors in (12) and ( 13) are assumed to be unknown Gaussian variables, we treat the precision parameters γ p and γ d as unknown i.i.d. Gamma random variables and estimate them within the VB framework.
Notably, in the ideal scenario where both CE and data detection are error-free, the residual inter-user interference simplifies to the noise term alone. As a result, γ p and γ d represent the noise precision.
Based on (6), in order to apply the VB method, it is necessary to compute the joint distribution,
, which can be derived as follows:
In this study, our objective is to compute the Bayesian optimal estimates for both the data symbol X d and the channel matrix H. Achieving this requires the posterior distribution p(X d , H, γ p , γ d |R p , R d , X p , Σ), which is often difficult to determine. Hence, we use the VB framework to approximate this posterior distribution, with detailed explanations provided in the subsequent section.
IN CF-MMIMO In this section, we propose a VB-based method to address the JED problem in the uplink of a CF-mMIMO network by approximating the intractable posterior distribution p(X d , H, γ p , γ d |R p , R d , X p , Σ). To this end, we adopt a mean-field variational distribution q(X d , H, γ p , γ d ) within the VB framework, as described in the following.
Based on (6), computing the optimal variational densities in (15) requires the joint distribution p(R p , R d , X d , H, γ p , γ d ; X p , Σ), given by:
To use (16) within the VB framework, we need to specify the prior distributions for p(x i,t ) and p(γ d,t ) for T p + 1 ≤ t ≤ T , p(γ p ), and p(h i,ℓ ; Σ i,ℓ ). We assume the prior distribution of x i,t as p(x i,t ) = a∈S p a δ(x i,t -a), where p a represents the probability of the constellation point a ∈ S, and S denotes the signal constellation. Moreover, we assume
JED processing is performed at the CPU and can operate on either unquantized or quantized signals forwarded by the APs. While forwarding unquantized signals requires high fronthaul bandwidth, transmitting quantized signals from the APs to the CPU reduces the bandwidth demand. To analyze these situations, we consider the following three scenarios.
In the PFL scenario, we assume that all APs forward their unquantized received signals to the CPU. We then use the CAVI algorithm to find the local optimal solutions: q ⋆ (h i,ℓ ), q ⋆ (x i,t ), q ⋆ (γ d,t ), and q ⋆ (γ p ). In the following subsections, we describe the update procedure for each of these variables.
1) Updating h i,ℓ : By evaluating the expectation of ( 16) over all latent variables except h i,ℓ , we derive the variational distribution q(h i,ℓ ) as shown in (17) (at the top of the next page). Given that h i,ℓ is assumed to follow a Gaussian prior distribution, it is reasonable to guess that q(h i,ℓ ) will also be Gaussian. Our derivations in (17) confirm this point. Specifically, q(h i,ℓ ) is Gaussian with the following covariance matrix and mean:
We use the following lemma to compute the variational posterior mean of x i,t , γ p , and γ d . Lemma 1. [35] Let y ∈ C m×1 , A ∈ C m×n , and x ∈ C n×1 be three independent random matrices (vectors) with respect to a variational density q y,A,x (y, A, x) = q y (y)q A (A)q x (x). Suppose A is column-wise independent and ⟨a i ⟩ and Σ ai are the variational mean and covariance matrix of the i th column of A. Let ⟨y⟩ and Σ y be the variational mean and covariance matrix of y and ⟨x⟩ and Σ x be the variational mean and covariance matrix of x. Consider F is an arbitrary Hermitian matrix. Here, (y -Ax) H F(y-Ax) with respect to q y,A,x (y, A, x), is given by:
where B = diag Tr{FΣ a1 }, . . . , Tr{FΣ an } .
Proof: The proof of this lemma is provided in [35].
2) Updating x i,t : We apply this update exclusively during the data transmission phase. We compute the expectation of ( 16) over all latent variables, excluding x i,t , to derive the variational distribution q i (x i,t ), as shown in (21) on the following page, where
and
serves as a linear approximation for x i,t . To evaluate ∥h i ∥ 2 in (23), we utilize Lemma 1, giving:
where Σi = [ Σi,1 , Σi,2 , . . . , Σi,L ]. Since the prior distribution p(x i,t ) is discrete, the variational distribution q i (x i,t ) is also discrete. Hence, normalization is required such that
As a result, the expected value and variance of x i,t under the variational distribution are given by:
3) Updating γ p : Here, we determine the variational distribution q(γ p ) by taking the expectation of ( 16) with respect to all latent variables except γ p , leading to the following expression:
Based on (28), q(γ p ) is Gamma distribution with mean
where
4) Updating γ d,t : Similar to the VB process for γ p , we have:
which results in q(γ d,t ) being Gamma distribution with the following mean:
where
If we set a p = b p = a d,t = b d,t = 0 in ( 29) and ( 31), the expectations ⟨γ p ⟩ and ⟨γ d,t ⟩ become equivalent to the estimated values of the precision of the residual inter-user interference described in [35] using the EM technique. Hence, the VB-EM method can be seen as a special case of our proposed VB method, which employs Gamma priors for the residual terms.
It is essential to note that, in our VB framework, γ p is assumed to be a scalar due to the single-shot nature of the CE process over the entire pilot block, whereas γ d is modeled as a vector to better capture the dynamic and iterative nature of the JED process during the data transmission phase.
As discussed earlier, the PFL scenario requires APs to forward unquantized signals to the CPU, which typically demands high-resolution quantizers to approximate continuous signals [32]. However, this results in a high fronthaul bandwidth requirement, making the approach impractical in real-world
systems. To address this challenge, we introduce the Q-E and E-Q scenarios, explained in the following subsections.
In this case, each AP quantizes its received signal to b bits and then forwards it to the CPU. Then, the CPU employs the quantized received signals to perform JED. Here, y ℓ,t is the quantized received signal vector at the ℓ th AP at duration t, which is given by: (32) where Q b (•) represents a quantization operator with b bits, utilizing a uniform scalar quantization approach. It operates based on a set of 2 b -1 thresholds, denoted as {d 1 , d 2 , . . . , d 2 b -1 }. For convenience, we define the thresholds to satisfy
The step size for quantization, denoted as ∆, is used to express the quantization thresholds as d i = (-2 b-1 + i)∆ for i = 1, 2, . . . , 2 b -1. This results in a quantized output q as follows:
Further, we use q low = d i-1 and q up = d i to represent the lower and upper bounds of the quantization bin that contains q. Then, in order to apply the VB framework, we need to formulate the joint probability as (34), shown at the top of the next page, where y t is the stacked vector of y ℓ,t from all APs at duration t,
Here, our goal is to use the VB framework to get the variational distribution
We then iteratively apply the CAVI algorithm to obtain the local optimal solutions q ⋆ (r t ), q ⋆ (h i,ℓ ), q ⋆ (x i,t ), q ⋆ (γ p ), and q ⋆ (γ d,t ). The derivations for these variational distributions, except for q ⋆ (r t ), follow directly from the procedures established for the PFL scenario in Section III-A. Therefore, we focus on detailing the steps for obtaining q ⋆ (r t ) given the quantized received signal y t in the following part.
1) Updating r t : This update contains two parts: Pilot Transmission Phase: By computing the expectation of ( 34) for all latent variables except for r t when 1 ≤ t ≤ T p , we express the variational distribution q(r t ) as:
∝ 1(r t ∈ [y low t , y up t ]) × exp -⟨γ p ⟩∥r t -⟨H⟩x t ∥ 2 . Notice the variational distribution q(r t ) is inherently separable as M m=1 L ℓ=1 q(r m,ℓ,t ) and the variational distribution q(r m,ℓ,t ) is the complex Gaussian distribution obtained from bounding r m,ℓ,t ∼ CN (⟨[H ℓ ] :m ⟩x t , γ -1 p ) to the interval (y low m,ℓ , y up m,ℓ ). Data Transmission Phase: Similar to the pilot transmission phase, we have: q(r t ) ∝ exp {⟨ln p(y t |r t ) + ln p(r t |x t , H, γ d,t )⟩} (37)
∝ 1(r t ∈ [y low t , y up t ]) × exp -⟨γ d,t ⟩∥r t -⟨H⟩⟨x t ⟩∥ 2 , when T p + 1 ≤ t ≤ T . Here, q(r m,ℓ,t ) is also a complex Gaussian distribution obtained from bounding r m,ℓ,t ∼ CN (⟨[H ℓ ] :m ⟩⟨x t ⟩, γ -1 d,t ) to the interval (y low m,ℓ,t , y up m,ℓ,t ). Later, in the Appendix, we will use U and V functions to represent the variational mean ⟨r m,ℓ,t ⟩ and variance τ r m,ℓ,t , respectively. Algorithm 1 presents a pseudocode for performing the JED process using our proposed VB method in the Q-E scenario, where I tr denotes the number of iterations in the CAVI algorithm. We will discuss the fronthaul signaling overhead required in the Q-E scenario in detail later in Section IV-D.
It is important to mention that the PFL scenario follows the steps outlined in Algorithm 1, except for steps 7 to 9, 12, and 13, which are skipped.
In the two previous scenarios, the CPU is responsible for the entire process. However, it could be possible to leverage APs to partially process the received signals, then forward the processed signals to the CPU, which would complete the remaining processing. To achieve this, we introduce the E-Q scenario, which includes two parts: i.) AP pre-processing in the pilot transmission phase, and ii.) CPU processing during the data transmission phase. In particular, the ℓ th AP performs
Algorithm 1: VB(Q-E) Method Find ⟨r m,ℓ,t ⟩, ∀m∀ℓ, using ( 36) and (55). 9 Get τ r m,ℓ,t , ∀m∀ℓ, using ( 36) and (56). Attain ⟨r m,ℓ,t ⟩, ∀m∀ℓ, using (37) and (55).
13
Obtain τ r m,ℓ,t , ∀m∀ℓ, using (37) and (56).
14 Apply ( 19) and ( 18) to get h i,ℓ and Σi,ℓ , ∀i∀ℓ, respectively.
15 for t = Tp + 1, . . . , T do 16 for i = 1, 2, . . . , K do 17 Get qi(xi,t) as in (21). 18 Calculate ⟨xi,t⟩ and τ x i,t based on (26) and (27). 19 Set r m,ℓ,t = ⟨r m,ℓ,t ⟩, h i,ℓ = ⟨h i,ℓ ⟩, Σ i,ℓ = Σi,ℓ . 20 Set xi,t = ⟨xi,t⟩, γp = ⟨γp⟩, γ d,t = ⟨γ d,t ⟩. 21 until convergence; 22 Compute xi,t = arg max a∈S qi(a), ∀i and t ∈ [Tp + 1, T ]. 23 Calculate Ĥ = H.
local processing relying only on its received signal, r ℓ,t , 1 ≤ t ≤ T p , to find h
as a local estimate of h i,ℓ . Then, it computes h q i,ℓ = Q(h loc i,ℓ ) as a quantized version of h loc i,ℓ and forwards it to the CPU. Next, in the data transmission phase, the ℓ th AP sends h q i,ℓ and y ℓ,t to the CPU and lets the CPU follow the VB framework to do JED based on H q and Y d , where H q is the stacked matrix containing all h q i,ℓ , ∀i, ∀ℓ. In the following, we explain the details of the VB approach for AP pre-processing and CPU processing tasks.
AP pre-processing: In this part, the ℓ th AP performs CE using its unquantized received signal r ℓ,t . To do so, it follows the CAVI algorithm to update unknown variables h loc i,ℓ and γ loc p,ℓ , where γ loc p,ℓ is the local estimate of γ p at the ℓ th AP. 1) Updating h loc i,ℓ : By computing the expectation of ( 16) for all latent variables except for h i,ℓ , and by taking the point into account that only r ℓ,t is available the ℓ th AP, we express the variational distribution q(h loc i,ℓ ) as (38) (at the top of the next page), which is Gaussian with the subsequent covariance matrix and mean:
2) Updating γ loc p,ℓ : We take the expectation of ( 16) with respect to all latent variables except for γ p , to derive the variational distribution q(γ loc p,ℓ ) as follows:
which leads to
where
. CPU processing: In this part, the CPU utilizes the information sent by all APs (i.e., H q and Y d ) and applies the VB method to obtain the variational distribution q
Next, we leverage the CAVI algorithm to find the local optimal solutions for the variational distributions q(r t ), q(h i,ℓ ), and q(x i,t ), and q(γ d,t ). To achieve this, we need to compute the joint probability, which is given by (44) on the next page, where H loc is the stacked matrix containing all h loc i,ℓ , ∀i∀ℓ. In this part, the updating procedures are analogous to what we presented in Section III-B, except for h i,ℓ , by ignoring
the information corresponding to the pilot transmission phase. Thus, we only provide the details about finding q(h i,ℓ ).
3) Updating h i,ℓ : We compute the expectation of ( 44) for all latent variables except for h i,ℓ . We express the variational distribution q(h i,ℓ ) as (45) (on the next page), where hi,ℓ and Σi,ℓ are given by:
hi,ℓ = Σi,ℓ
It is essential to note that since h loc i,ℓ is estimated from local observations at the ℓ th AP, it typically differs from the true channel h i,ℓ . We model this mismatch as h loc i,ℓ = h i,ℓ + e h,i,ℓ , where e h,i,ℓ is the corresponding CE error, and each element of e h,i,ℓ is assumed to be i.i.d. zero-mean complex Gaussian with known variance N e i,ℓ . Under this model, (45) describes the distribution of a quantized version of a noisy Gaussian variable. Moreover, the variational distribution q(h i,ℓ ) is inherently separable as M m=1 K i=1 L ℓ=1 q(h m,i,ℓ ). Therefore, the corresponding mean and variance of ℜ{h m,i,ℓ }, denoted by ⟨ℜ{h m,i,ℓ }⟩ and [ Σ [1] i,ℓ ] mm , respectively, can be computed using the approach proposed in [34] as (48) and (49), shown on the next page, where u(x) and U (x) represent the PDF and CDF of the random variable x, respectively.
Similar to (48), ( 49), (50), and (51), by replacing the real parts with imaginary parts, we can derive ⟨ℑ{h m,i,ℓ }⟩ and [ Σ [2] i,ℓ ] mm as the corresponding mean and variance of ℑ{h m,i,ℓ }, respectively. Finally,
Algorithm 2 provides a pseudocode for implementing the VB(E-Q) method to carry out the JED process. The AP preprocessing is detailed in steps 3 through 10, while the CPU processing is performed in steps 11 through 28. In the next section, we will provide a detailed discussion on the fronthaul signaling overhead requirements for the E-Q scenario.
To assess the performance of our proposed VB methods, in this section, we focus on the uplink scenario in a CF-mMIMO network. We conduct a comparison between our proposed VB methods (i.e., VB(PFL), VB(Q-E), and VB(E-Q)) and both linear and nonlinear JED methods in terms of SER, channel NMSE, computational complexity, and signaling overhead over the fronthaul links.
We evaluate the performance of our VB methods using quadrature phase shift keying (QPSK) modulation. Unless stated otherwise, we set L = 8, M = 4, K = 16, I tr = 50, T p = 32, T d = 128, and p a = 1/|S|, where |S| denotes the cardinality of S. Moreover, we normalize the covariance matrix Σ i,ℓ , ∀i, ∀ℓ, so that all diagonal elements are equal to 1. This ensures that E[∥h i,ℓ ∥ 2 ] = M . We then determine the noise variance N 0 based on the signal-to-noise ratio (SNR), which is expressed as:
In this paper, to model the spatial correlation, we adopt the exponential correlation model [36], applied independently to each column of H ℓ . The corresponding covariance matrix Σ i,ℓ is defined as:
, ∀i∀ℓ, prior distribution of pa, ∀a ∈ S, ap, bp, and
12 Initialize τ r m,ℓ,t = 0, ∀m∀ℓ∀t. Attain ⟨r m,ℓ,t ⟩, ∀m∀ℓ, using (37) and (55). 17 Obtain τ r m,ℓ,t , ∀m∀ℓ, using (37) and (56).
18
Apply (48) to get h i,ℓ , ∀i∀ℓ. Get qi(xi,t) as in (21).
Set r m,ℓ,t = ⟨r m,ℓ,t ⟩, h i,ℓ = ⟨h i,ℓ ⟩, Σ i,ℓ = Σi,ℓ . by averaging over 100 trials.
In this part, we evaluate the SER performance of our proposed VB methods across various cases, including different SNR levels, number of users, number of APs, pilot transmission length, and data transmission length.
Fig. 2 shows the SER as a function of SNR for VB(PFL), LMMSE(PFL), and VB-DD(PFL), as well as VB(Q-E) and VB(E-Q) with i.i.d. channels, evaluated with 1-bit, 2-bit, and 3-bit quantizer configurations, as the SNR ranges from 0 to 20 dB. In our VB-based methods, we set a p = b p = a d,t = b d,t = 0. The LMMSE(PFL) approach performs JED by first estimating the channels during the pilot phase using known pilot signals and then utilizing the estimated channels to detect data during the data transmission phase. Furthermore, the VB-DD(PFL) benchmark is a data detection scheme that has access to perfect channel knowledge and PFL, and applies VB inference solely for data detection.
It is evident that the performance of all methods improves with higher SNRs, as the signal becomes stronger in comparison to the noise. Among these, VB-DD(PFL) demonstrates the best performance due to its nonlinear Bayesian model and access to perfect channel knowledge and unquantized received signals from all APs. VB(PFL) ranks second, as it uses unquantized signals but is affected by CE errors inherent to the JED process. The VB(Q-E, 3bits) and VB(E-Q, 3bits) methods take rank third due to their use of the quantized information. Finally, the LMMSE(PFL) method ranks fourth. Despite the CPU having access to received signals from all APs, this method relies on a linear model, resulting in worse performance compared to the VB-based methods with unquantized or 3-bit quantized signals. Notably, VB(PFL), VB(Q-E, 3bits), and 3bits) demonstrate gains of about 4 dB, 2 dB, and 2 dB, respectively, at an SER of 10 -3 , compared to the LMMSE(PFL) method.
Additionally, Fig. 2 illustrates that as the resolution of the quantized signals decreases, the performance of both the VB(Q-E) and VB(E-Q) methods deteriorates. For lowresolution cases, the performance curves for these methods saturate at high SNRs. This occurs because, at high SNRs, quantization noise becomes dominant over background noise, rendering further increases in SNR ineffective for performance improvement. Furthermore, as we can see, the performance of VB(E-Q, 2bits) and VB(E-Q, 1bit) is slightly lower than that of their counterparts based on VB(Q-E), consistent with the results reported in [30]. This is due to CE errors during the estimation phase, where each AP estimates the channels using only its locally received signals.
Fig. 3 presents the SER performance of the LMMSE(PFL), VB(PFL), VB(Q-E, 3bits), and VB(E-Q, 3bits) methods for T p ∈ [10, 64] with i.i.d channels when a p = b p = a d,t = b d,t = 0. It is worth noting that VB methods do not strictly require orthogonal pilot signals (i.e., T p ≥ K). The results show that the performance of all methods improves as T p increases. Initially, the SER drops rapidly, but the rate of improvement slows down as T p becomes larger. This happens because, in the JED process during the data transmission phase, the SER performance is influenced by the prior channel knowledge obtained during the pilot phase. When T p is small, the initial CE is poor due to limited pilot observations, resulting in higher SER. As T p increases, the prior CE improves, leading to better detection performance. However, once T p becomes sufficiently large (i.e., T p ≥ 48), the performance gain becomes marginal since the JED process iteratively refines the CE during the data transmission phase. Thus, further increasing T p yields limited additional benefits.
Moreover, VB(Q-E, 3bits) and VB(E-Q, 3bits) reach the saturation point earlier. This is because, after achieving a reasonable CE, their performance becomes limited by quantization noise rather than estimation accuracy. Further, as anticipated, the VB-based methods outperform the LMMSE(PFL) method, benefiting from the nonlinear VB framework.
Next, in Fig. 4, we evaluate the SER performance of LMMSE(PFL), VB(PFL), VB(Q-E, 3bits), and VB(E-Q, 3bits) under i.i.d channels with respect to T d , where T d varies from 16 to 192 and a p = b p = a d,t = b d,t = 0. This figure shows that VB(PFL) outperforms the others due to access to unquantized signals from all APs and the nonlinear VB framework. In contrast, the performance of LMMSE(PFL) remains constant because it only performs CE during the pilot transmission phase and does not involve a recursive optimization process. Consequently, any detection errors that occur during short data transmission periods persist even in scenarios with long T d values.
Then, Fig. 5 shows the SER as a function of the number of users for LMMSE(PFL) and three variations of our VB- based methods with i.i.d channels, when K ∈ [12,32] and a p = b p = a d,t = b d,t = 0. As observed, the performance of all methods degrades as K increases. This can be explained to the fact that a network with K 1 users can be approximated as one with K 2 users, where K 2 < K 1 , by deactivating K 1 -K 2 users. However, when these deactivated users become active, they introduce additional interference to the JED process. The trend in Fig. 5 aligns with the previous figures, with VB(PFL) providing the best performance and the VB-based methods consistently outperforming the LMMSE(PFL) method.
We also provide an SER analysis of the LMMSE(PFL), VB(PFL), VB(Q-E, 3bits), and VB(E-Q, 3bits) methods with respect to different numbers of APs (i.e., L ∈ [4, 10]) in Fig. 6 under i.i.d channels with a p = b p = a d,t = b d,t = 0. This figure shows that the SER performance of all methods improves as L increases, due to the availability of a larger number of received signals. Furthermore, our proposed VB-based methods achieve a lower SER compared to LMMSE(PFL), owing to their use of the nonlinear VB framework.
Given that our proposed VB-based methods are nonlinear, it is important to evaluate their performance against other nonlinear benchmarks. To this end, we consider the GAMP(PFL) and GAMP(Q-E) methods introduced in [34], along with the VB-EM(PFL) and VB-EM(Q-E) techniques from [35]. Note that while these benchmarks were originally developed for mMIMO systems, they can be extended to CF-mMIMO sys-tems. Moreover, since this work is the first to introduce an E-Q approach for the JED problem and no direct counterpart exists in the literature, we limit our comparison to the VB(PFL) and VB(Q-E) methods against the aforementioned nonlinear baselines. Fig. 7 presents the SER performance under correlated channels with α = 0.35 + j0. 35, where the bisection method is used to determine the optimal values of a p , b p , a d,t , and b d,t . Other settings are configured analogously to Fig. 2. Here, VB(PFL) outperforms the other PFL-based methods, while VB(Q-E, 3bits) achieves superior performance than other Q-E-based approaches. This improvement is attributed to the fact that GAMP(PFL) and GAMP(Q-E) are not robust under correlated channels and may diverge. Furthermore, while VB-EM(PFL) and VB-EM(Q-E) in [35] use the EM algorithm to estimate precision, our approach generalizes this by modeling the precision parameters using a Gamma distribution, of which the EM-based solution is a special case.
In this part, we focus on the channel NMSE metric to assess the performance of our VB-based methods. We calculate the channel NMSE using the following formula in dB scale.
In Fig. 8, we evaluate the channel NMSE performance of our VB-based methods compared to LMMSE(PFL) and VB-CE(PFL), where the simulation settings are the same as in Fig. 2. The VB-CE(PFL) method is a CE approach that uses VB inference while treating the entire transmission block as pilot data. Our VB-based methods outperform LMMSE(PFL) for the reasons discussed in Fig. 2. Further, at low SNRs, a noticeable performance gap exists between our methods and VB-CE(PFL) due to the limited informativeness of data symbols in noisy conditions. However, this gap narrows as the SNR increases, since the data becomes more informative and contributes more effectively to the CE process.
Next, in Fig. 9 and Fig. 10, we study the channel NMSE performance across different values of T p and T d , respectively, where the simulation settings are the same as in Fig. 3 and Fig. 4. In these figures, we compare LMMSE(PFL) with VB(PFL), VB(Q-E, 3bits), and VB(E-Q, 3bits). The results are similar to those in Fig. 3 and Fig. 4, respectively, where VB(PFL) outperforms the others. This observation shows that the CE and data detection procedures interactively affect each other, and improving one leads to improvements in the other.
In this subsection, we focus on the comparative analysis of computational complexity between LMMSE(PFL), VB(PFL), GAMP(PFL), VB-EM(PFL), VB(Q-E), GAMP(Q-E), VB-EM(Q-E), and VB(E-Q). In particular, the complexity of LMMSE(PFL) is O(M LK 2 + |S|K) [14]. Moreover, the complexity of VB(PFL) is expressed as O(I tr M 3 L 3 K + (T d + 1)M LK + |S|K ), due to the matrix inversion in (18) and matrix multiplications in (29) and (31). The complexity of GAMP(PFL) and VB-EM(PFL) matches that of VB(PFL) [34], [35]. Then, the complexity of VB(Q-E) is O(I tr M 3 L 3 K + (T d + 1)M LK + |S|K ). Furthermore, the complexity of GAMP(Q-E) and VB-EM(Q-E) is equal to the complexity of VB(Q-E). Finally, the complexity of VB(E-Q) is given by O(I tr 2M 3 LK +(T d +2)M K +|S|K ), due to the matrix inversions in (39) and (46), and matrix multiplications in (42) and within the CPU process. As a result, LMMSE(PFL) has the lowest complexity, followed by VB(E-Q), which has the second-lowest complexity, while the other methods based on PFL and Q-E exhibit the highest complexity.
In this part, we analyze the fronthaul signaling overhead requirements of the methods based on PFL, Q-E, and E-Q, as discussed in the previous subsection. In real-world scenarios, the unquantized signals at the APs must first be quantized (with high resolution) before being forwarded to the CPU. For this purpose, we assume a 10-bit resolution to represent the unquantized signals. In this case, the PFL-based methods require 10M L(T p + T d ) bits for fronthaul signaling overhead. In contrast, the VB(Q-E, 3bits), GAMP(Q-E, 3bits), and VB-EM(Q-E, 3bits) methods use 3-bit quantizers at the APs, where each AP quantizes its received signals and forwards them to the CPU, requiring 3M L(T p + T d ) bits. On the other hand, VB(E-Q, 3bits) involves local CE at each AP, followed by quantization and forwarding of the quantized channels to the CPU, requiring 3M L(K + T d ) bits. This is lower than the methods based on (Q-E, 3bits) as K ≤ T p .
Based on the previous two subsections, we observe that, from a computational complexity perspective, among the VBbased methods, VB(PFL) and VB(Q-E) have the highest complexity. On the other hand, in terms of fronthaul signaling overhead, VB(PFL) requires the highest fronthaul load, and VB(Q-E, 3bits) imposes a higher load than VB(E-Q, 3bits). Therefore, there is a trade-off among the different VB-based methods. Specifically, if fronthaul bandwidth is not a concern, VB(PFL) is the optimal choice. If computational complexity is not a constraint but fronthaul bandwidth is limited, the VB(Q-E, 3bits) method is more favorable. Finally, if both computational complexity and fronthaul bandwidth are limited, the VB(E-Q, 3bits) method is the most suitable option.
In this work, we studied the JED problem in CF-mMIMO networks. We proposed a VB-based approach to address the limitations of fronthaul bandwidth by using low-bit quantizers at the APs. Specifically, we considered two scenarios: Q-E and E-Q. In the Q-E scenario, each AP sends a quantized version of its received signals to the CPU, while in the E-Q scenario, each AP initially performs local CE during the pilot transmission phase and then shares a quantized version of the estimated channels, along with a quantized version of the received signals during the data transmission phase, with the CPU. We evaluated our VB-based approach under PFL, Q-E, and E-Q scenarios and compared its performance with LMMSE(PFL), GAMP(PFL), GAMP(Q-E), VB-EM(PFL), and VB-EM(Q-E). Notably, VB(Q-E) and VB(E-Q) outperform LMMSE(PFL) due to the nonlinear VB framework. In addition, our VBbased methods outperform the nonlinear benchmarks thanks to their efficient convergence to a local optimum and the superior performance of the pure VB framework compared to approaches that rely on the EM technique. Future research can proceed in several promising directions. One direction is to explore the impact of near-field communication effects on the JED process in CF-mMIMO systems. Moreover, in any Bayesian-based statistical CE approach, including VB methods, knowledge of the channel covariance matrix is essential for effective estimation. This information can, in principle, be