Deep-Unfolded Massive Grant-Free Transmission
in Cell-Free Wireless Communication Systems

This paper proposes a novel framework for joint active user detection, channel estimation, and data detection (JACD) in cell-free systems with grant-free access. We start by formulating JACD as an optimization problem that fully exploits the sparsity of the wireless channel and data matrices and then approximately solve it using forward-backward splitting (FBS) [31, 32]. To enable FBS, we relax the discrete constellation constraint to its convex hull and employ JACD methods with the incorporation of either a regularizer or a posterior mean estimator (PME), guiding symbols toward discrete constellation points, resulting in the box-constrained FBS and PME-based JACD algorithms, respectively. To reduce complexity, we replace the exact proximal operators with approximate shrinkage operations and approximate PME computations. To improve convergence and JACD performance, we include per-iteration step sizes and a momentum strategy. To avoid tedious manual parameter tuning, we employ deep unfolding (DU) to jointly tune all of the algorithm hyper-parameters using machine learning tools. To improve active user detection (AUD) performance, we include a novel soft-output AUD module that jointly considers the estimated data and channel matrix. Based on the aforementioned modifications, we have developed the deep unfolding versions of the box-constrained FBS and PME-based JACD algorithms, referred to as DU-ABC and DU-POEM. We use Monte–Carlo simulations to demonstrate the superiority of our framework compared to existing methods.

Matrices, column vectors, and sets are denoted by uppercase boldface letters, lowercase boldface letters, and uppercase calligraphic letters, respectively. $\mathbf{A}(m,n)$ stands for the element of matrix $\mathbf{A}$ at the $m$th row and $n$th column, and $\mathbf{a}(m)$ stands for the $m$th entry of the vector $\mathbf{a}$. The $M\times N$ all-ones matrix is $\mathbf{1}_{M\times N}$ and the $M\times M$ identity matrix is $\mathbf{I}_{M}$. The unit vector, which is zero except for the $m$th entry, is $\mathbf{e}_{m}$, and the zero vector is $\mathbf{0}$; the dimensions of these vectors will be clear from the context. We use hat symbols to refer to the estimated values of a variable, vector, or matrix. The superscripts $\left(\cdot\right)^{T}$ and $\left(\cdot\right)^{H}$ represent transpose and conjugate transpose, respectively, and the superscript $(k)$ denotes the $k$th iteration. The Frobenius norm is denoted by $\|\cdot\|_{F}$. In addition, $\operatorname{diag}\{x_{1},\ldots,x_{N}\}$ stands for a diagonal matrix with entries $x_{1},\ldots,x_{N}$ on the main diagonal; $\operatorname{det}(\mathbf{A})$ stands for the determinant of $\mathbf{A}$. The cardinality of a set $\mathcal{Q}$ is $|\mathcal{Q}|$. The operators $\odot$ and $\propto$ denote Hadamard product and proportional relationships. For $\mathbf{x}\in\mathbb{C}^{N}$, $\text{Re}\{\mathbf{x}\}\in\mathbb{R}^{N}$ and $\text{Im}\{\mathbf{x}\}\in\mathbb{R}^{N}$ represent its real and imaginary part, respectively. $\mathbb{P}\{\cdot\}$ and $\mathbb{E}\{\cdot\}$ denote probability and expectation, respectively; the indicator function $\mathbb{I}\left\{\cdot\right\}$ returns $1$ for valid conditions and $0$ otherwise. The multivariate complex Gaussian probability distribution with mean vector $\mathbf{m}$ and covariance matrix $\bm{\Sigma}$ evaluated at $\mathbf{x}\in\mathbb{C}^{M}$ is

and the symbol $\triangleq$ means “is defined as.” Besides, the shrinkage operation is defined as

and the element-wise clamp function is defined as

The rest of this paper is organized as follows. Section II presents the necessary prerequisites for the subsequent derivations. Section III introduces the system model and formulates the JACD optimization problem. Section IV introduces two distinct JACD algorithms: (i) box-constrained FBS and (ii) PME-based JACD. Section V deploys DU techniques and details the AUD module. Section VI investigates the efficacy of the proposed algorithms via Monte–Carlo simulations. Section VII concludes the paper. Proofs are relegated to the appendices.

FBS, also known as proximal gradient methods, is widely used for solving a wide variety of convex optimization problems [31, 32]. FBS splits the objective function into two components: a smooth function, denoted as $f(\mathbf{S})$, and another not necessarily smooth function, $g(\mathbf{S})$, and solves the following optimization problem:

FBS iteratively performs a gradient step in the smooth function $f(\mathbf{S})$ (denoted as the forward step) and a proximal operation to find a solution in the vicinity of the minimizer of the function $g(\mathbf{S})$ (denoted as the backward step). The forward step proceeds as

where the superscript $(k)$ denotes the $k$th iteration, $\tau^{(k)}$ represents the step size at iteration $k$, and $\nabla f\left(\mathbf{S}\right)$ is the gradient of $f(\mathbf{S})$. The backward step proceeds as

This process is iterated for $k=1,2,\ldots$ until a predefined convergence criterion is met or a maximum number of $K$ iterations has been reached.

Here, we introduce Definitions 1 and 2 to specify the probability distributions that the random vector $\mathbf{x}$ may follow.

Definition 1: For a random vector $\mathbf{x}$ taken from the discrete set $\mathcal{S}$ ($\mathbf{0}\in\mathcal{S}$), we call $\mathbf{x}$ follows the $\theta$-mixed discrete uniform distribution on $\mathcal{S}$, denoted as $\mathbf{x}\sim U_{\theta,\mathcal{S}}\{\mathbf{x}\}$, if $\mathbb{P}\{\mathbf{x}\}=U_{\theta,\mathcal{S}}\{\mathbf{x}\}=\frac{\theta}{|\mathcal{S}|-1}\mathbb{I}\{\mathbf{x}\neq\mathbf{0}\}+(1-\theta)\mathbb{I}\{\mathbf{x}=\mathbf{0}\},\forall\mathbf{x}\in\mathcal{S}$, where $\theta\in(0,1)$ is the probability of $\mathbf{x}$ being non-zero.

Definition 2: For a random vector $\mathbf{x}$ in a discrete set $\mathcal{S}$ ($\mathbf{0}\notin\mathcal{S}$), we call $\mathbf{x}$ follows a discrete uniform distribution on $\mathcal{S}$, denoted as $\mathbf{x}\sim U_{\mathcal{S}}\{\mathbf{x}\}$, if $\mathbb{P}\{\mathbf{x}\}=U_{\mathcal{S}}\{\mathbf{x}\}=\frac{1}{|\mathcal{S}|},\forall\mathbf{x}\in\mathcal{S}$.

We also define the PME of a vector $\mathbf{x}$ under the specific conditions as follows:

Definition 3: Given the observation vector $\hat{\mathbf{x}}=\mathbf{x}+\mathbf{e}$ with $\mathbf{x}\in\mathcal{S}$ following the prior distribution $\mathbb{P}\{\mathbf{x}\}$ and Gaussian estimation error $\mathbf{e}\sim\mathcal{CN}\left(\mathbf{0},N_{\text{e}}\mathbf{I}\right)$, we call

the PME of $\mathbf{x}$ under the prior $\mathbb{P}\{\mathbf{x}\}$.

Following our previous work in [1, 2], we model the input-output relation of this scenario as follows:

Here, the received signal matrix from all APs is denoted by $\mathbf{Y}\in\mathbb{C}^{MP\times R}$, and $\xi_{n}\in\{0,1\}$ is $n$th UE’s activity indicator with $\xi_{n}=1$ indicating the $n$th UE is active and $\xi_{n}=0$ otherwise. We assume that the activity between UEs is independent and all UEs have the same activity probability $\alpha$, i.e., $\mathbb{P}\{\xi_{n}=1\}=\alpha,\forall n$. The channel vector between the $n$th UE and all APs is $\mathbf{h}_{n}\triangleq[\mathbf{h}_{n,1}^{T},\ldots,\mathbf{h}_{n,P}^{T}]^{T}\in\mathbb{C}^{MP}$ with $\mathbf{h}_{n,p}\in\mathbb{C}^{M}$ being the channel vector between the $n$th UE and the $p$th AP. The $n$th UE’s signal vector $\bar{\mathbf{x}}_{n}\triangleq[\mathbf{x}_{\text{P},n}^{T},\bar{\mathbf{x}}_{\text{D},n}^{T}]^{T}\in\mathbb{C}^{R}$ consists of the pilot vector $\mathbf{x}_{\text{P},n}\in\mathbb{C}^{R_{\text{P}}}$ and the data vector $\bar{\mathbf{x}}_{\text{D},n}\in\mathcal{Q}^{R_{\text{D}}}$ with data entries independently and uniformly sampled from the constellation set $\mathcal{Q}$. Entries in the noise matrix $\mathbf{N}\in\mathbb{C}^{MP\times R}$ are assumed to be i.i.d. circularly-symmetric complex Gaussian variables with variance $N_{0}=1$.

For ease of notation, we rewrite the system model in (8) as

where $\mathbf{Y}=[\mathbf{Y}_{\text{P}},\mathbf{Y}_{\text{D}}]$ with $\mathbf{Y}_{\text{P}}\in\mathbb{C}^{MP\times R_{\text{P}}}$ and $\mathbf{Y}_{\text{D}}\in\mathbb{C}^{MP\times R_{\text{D}}}$ being the received pilot matrix and received data matrix, respectively. The channel matrix is given by $\mathbf{H}\triangleq\left[\xi_{1}\mathbf{h}_{1},\ldots,\xi_{N}\mathbf{h}_{N}\right]\in\mathbb{C}^{MP\times N}$. The signal matrix $\mathbf{X}=[\mathbf{X}_{\text{P}},\mathbf{X}_{\text{D}}]\in\mathbb{C}^{N\times R}$ contains the pilot matrix $\mathbf{X}_{\text{P}}\triangleq[\mathbf{x}_{\text{P},1},\ldots,\mathbf{x}_{\text{P},N}]^{T}\in\mathbb{C}^{N\times R_{\text{P}}}$ and the data matrix $\mathbf{X}_{\text{D}}\triangleq[\mathbf{x}_{\text{D},1},\ldots,\mathbf{x}_{\text{D},N}]^{T}\in\mathbb{C}^{N\times R_{\text{D}}}$, where $\mathbf{x}_{\text{D},n}\triangleq\xi_{n}\bar{\mathbf{x}}_{\text{D},n}\in\bar{\mathcal{Q}}$ with $\bar{\mathcal{Q}}\triangleq\{\mathcal{Q}^{R_{\text{D}}},\mathbf{0}\}$.

A typical signal matrix $\mathbf{X}$ is depicted in Fig. 2. Here, the data matrix $\mathbf{X}_{\text{D}}$ is a row-sparse matrix due to sporadic UE activity, and the pilot matrix $\mathbf{X}_{\text{P}}$ is characterized as a non-sparse matrix, retaining all known pilot information for subsequent optimization. A typical channel matrix $\mathbf{H}$ is depicted in Fig. 3, where sparsity is due to two reasons: (i) the UEs’ sporadic activity results in column sparsity and (ii) the inherent discrepancies in large-scale fading between UEs and different APs, caused by their varying distances, lead to the block sparsity within each non-zero column.

While row/column sparsity and block sparsity have been widely studied in compressed sensing, the intricate interplay of these sparsities in our setting remains largely unexplored. In our system model (9), the estimation of both $\mathbf{H}$ and $\mathbf{X}$ poses a significant challenge, as $\mathbf{H}$ exhibits both column and row sparsity, while $\mathbf{X}_{\text{D}}$ displays row sparsity. In the subsequent problem formulation and algorithm derivation, we will effectively leverage the sparsity of both the channel matrix $\mathbf{H}$ and the data matrix $\mathbf{X}_{\text{D}}$ to enhance JACD performance.

Using the system model (9), we formulate the maximum-a-posteriori JACD problem for massive grant-free transmission in cell-free wireless communication systems as [1, 2]

where the channel law $P\!\left(\mathbf{Y}|\mathbf{H},\mathbf{X}_{\text{D}}\right)$ is given by

Here, we employ the complex-valued block-Laplace model and Laplace model for block sparsity in $\mathbf{H}$ and column sparsity in $\mathbf{X}_{\text{D}}$, respectively, as follows [1, 2]

with $\mu_{h}$ and $\mu_{x}$ indicating the sparsity levels of $\mathbf{H}$ and $\mathbf{X}_{\text{D}}$, respectively.

By inserting (11), (12), and (13) into (10), we can rewrite the JACD problem as problem $\mathcal{P}_{1}$ with $N_{0}=1$, expressed in problem (14) above. This problem aims to estimate the channel matrix $\mathbf{H}$ and the data matrix $\mathbf{X}_{\text{D}}$ using the received signal matrix $\mathbf{Y}$ and the pilot matrix $\mathbf{X}_{\text{P}}$, in which UE activity is indicated by the column sparsity of $\mathbf{H}$ and row sparsity of $\mathbf{X}_{\text{D}}$.

The discrete set $\bar{\mathcal{Q}}$ renders $\mathcal{P}_{1}$ a discrete-valued optimization problem for which a naïve exhaustive search is infeasible. To circumvent this limitation, as in [1, 2], we relax the discrete set $\bar{\mathcal{Q}}$ to its convex hull $\mathcal{B}^{R_{\text{D}}}$, thereby transforming $\mathcal{P}_{1}$ into a continuous-valued optimization problem. The set $\mathcal{B}$ is given by¹¹1For quadrature phase shift keying (QPSK), the set $\mathcal{B}$ can be specified as $\mathcal{B}=\{x\in\mathbb{C}:-B\leq\text{Re}\{x\}\leq B,-B\leq\text{Im}\{x\}\leq B\}$.

To push the entries in the data matrix $\mathbf{X}_{\text{D}}$ towards points in discrete set $\bar{\mathcal{Q}}$, two distinct methodologies can be used: Method (i) introduces a regularizer (also called penalty term) into the objective function of $\mathcal{P}_{1}$ [1, 64] and method (ii) leverages the PME to denoise the estimated data matrix $\mathbf{X}_{\text{D}}$ [2]. As such, we can transform the original discrete-valued optimization problem $\mathcal{P}_{1}$ into problem $\mathcal{P}_{2}$ in (16) above. The penalty parameter $\lambda$ is $\lambda>0$ for method (i) and $\lambda=0$ for method (ii). For the regularizer $\mathcal{C}\left(\mathbf{X}_{\text{D}}\right)$, many alternatives are possible, such as $\mathcal{C}\left(\mathbf{X}_{\text{D}}\right)=-\|\mathbf{X}_{\text{D}}\odot\mathbf{X}_{\text{D}}^{*}-B^{2}\mathbf{1}_{N\times R_{\text{D}}}\|_{F}^{2}$ utilized in [1]. In the next section, we will develop the FBS algorithm based on method (i) and method (ii) respectively to efficiently compute approximate solutions to problem $\mathcal{P}_{2}$.

For this method, we utilize the FBS with the incorporation of the regularizer $\mathcal{C}(\mathbf{X}_{\text{D}})$ to approximately solve the non-convex problem $\mathcal{P}_{2}$. Using the definition $\mathbf{S}\triangleq[\mathbf{H}^{H},\mathbf{X}_{\text{D}}]^{H}\in\mathbb{C}^{(MP+R_{\text{D}})\times N}$ in [1, 2, 54], we can split the objective function of $\mathcal{P}_{2}$ into the two functions

where $\mathcal{X}\left(\mathbf{X}_{\text{D}}\right)=+\infty\cdot\mathbb{I}\left\{\mathbf{X}_{\text{D}}\notin\mathcal{B}^{N\times R_{\text{D}}}\right\}$ represents the constraint ${\mathbf{X}}_{\text{D}}\in\mathcal{B}^{N\times R_{\text{D}}}$ in problem $\mathcal{P}_{2}$. Following the FBS framework outlined in Section II-A, the corresponding forward and backward steps can be specified as follows.

An alternative technique that considers the discrete constellation constraint in $\mathbf{X}_{\text{D}}$ uses a PME, which denoises the data matrix $\mathbf{X}_{\text{D}}$; for this approach, we set $\lambda=0$. We now only discuss the differences to the box-constrained FBS algorithm from Section IV-A with $\lambda=0$, which are in the computation of the proximal operator for $\mathbf{X}_{\text{D}}$.

In our paper, the PME assumes that one observes a signal of interest through noisy Gaussian observations (see Definition 3 in Section II-B) and obtains the best estimates in terms of minimizing the mean square error by leveraging the known prior probability of the signal. As in [65, 66], instead of calculating a proximal operator for the data matrix $\mathbf{X}_{\text{D}}$, we simply denoise the output of the forward step using a carefully designed PME. To this end, we model the $n$th UE’s estimated data vector $\hat{\mathbf{x}}_{\text{D},n}^{(k)}$ in the $k$th iteration as

Here, $\mathbf{x}_{\text{D},n}$ is unknown and $\mathbf{e}^{(k)}\in\mathbb{R}^{R_{\text{D}}}$ is the estimation error at iteration $k$, which we assume to be complex Gaussian following the distribution $\mathbf{e}^{(k)}\sim\mathcal{CN}(\mathbf{0},N_{\text{e}}^{(k)}\mathbf{I}_{R_{\text{D}}})$. Here, the vector $\mathbf{x}_{\text{D},n}$ from (26) follows $\alpha$-mixed discrete uniform distribution on $\bar{\mathcal{Q}}$, i.e., $\mathbf{x}_{\text{D},n}\sim U_{\alpha,\bar{\mathcal{Q}}}\{\mathbf{x}_{\text{D},n}\}$ (see Definition 1 in Section II-B), where $\alpha$ is the UE activity probability. Accordingly, we can employ the PME of $\mathbf{x}_{\text{D},n}$ under the prior $U_{\alpha,\bar{\mathcal{Q}}}\{\mathbf{x}_{\text{D},n}\}$, i.e., $\text{PME}(\hat{\mathbf{x}}_{\text{D},n}^{(k)},\bar{\mathcal{Q}},U_{\alpha,\bar{\mathcal{Q}}}\{\mathbf{x}\},N_{\text{e}}^{(k)})$ (refer to Definition 3 in Section II-B) as the estimate of the data matrix $\mathbf{X}_{\text{D}}$ [2], expressed as in (27). Algorithm 2 outlines the pseudocode for the proposed PME-based JACD algorithm.

After estimating the channel and data matrices over $K$ iterations using box-constrained FBS or PME-based JACD, we proceed with AUD and DD for both algorithms. As in [1], we determine UE activity based on the channel energy. Specifically, if the channel energy of the $n$th UE surpasses a threshold $T_{\text{AUD}}$, then the UE is deemed active; otherwise, we consider it inactive. Mathematically, we describe AUD as follows:

The estimated UE activity indicators $\{\hat{\xi}_{n}\}_{\forall n}$ can now be used to update the estimated data matrix as

where $\tilde{\mathbf{X}}_{\text{D}}\in\mathcal{Q}^{N\times R_{\text{D}}}$ is obtained by mapping the entries in $\mathbf{X}_{\text{D}}^{(K+1)}$ to the nearest symbols in $\mathcal{Q}$ as follows:

In Section V-C, we will introduce a trainable soft-output AUD module that leverages information from both the estimated channel and data matrices to improve AUD performance.

We assess the computational complexity of our algorithms using the number of complex-valued multiplications, which are expressed in big-O notation. For each iteration, the complexity of the forward steps in these algorithms is $O(MPNR+MPNR_{\text{D}}+\mathbb{I}\{\lambda\neq 0\}NR_{\text{D}}\mathcal{C}_{1})$. Here, $\mathcal{C}_{1}$ varies depending on the chosen regularizer²²2For instance, $\mathcal{C}_{1}=1$ for the regularizer $\mathcal{C}\left(\mathbf{X}_{\text{D}}\right)=-\|\mathbf{X}_{\text{D}}\odot\mathbf{X}_{\text{D}}^{*}-B^{2}\mathbf{1}_{N\times R_{\text{D}}}\|_{F}^{2}$.. The per-iteration complexities of the backward steps in the box-constrained FBS and PME-based JACD algorithms are $O(MPN+NR_{\text{D}})$ and $O(MPN+NR_{\text{D}}|\mathcal{Q}|^{R_{\text{D}}})$, respectively. In addition, the computational complexity of the AUD and DD module is $O(MPN+NR_{\text{D}})$.

The summation terms in the PME expression (27) lead to exponential complexity with the data length $R_{\text{D}}$ serving as the exponent, which results in higher complexity for the PME-based JACD algorithm compared to the box-constrained FBS algorithm. In Section V, we will leverage DU to tune hyper-parameters in these algorithms, improving their effectiveness automatically. Besides that, we introduce approximations to replace the costly computations in (25) and (27) in the backward steps of these algorithms to further reduce complexity.

The forward step and the proximal operations on $\mathbf{H}$ in the backward step of the DU-POEM algorithm align with those of the DU-ABC algorithm with $\lambda^{(k)}=0,\forall k$, as detailed in (31)-(35). We now only focus on the proximal operation for $\mathbf{X}_{\text{D}}$ in DU-POEM, which is different from that of the DU-ABC.

As for the PME of $\mathbf{X}_{\text{D}}$ in the PME-based JACD algorithm, the summation of numerous exponential terms in (27) can result in high computational complexity and lead to numerical stability issues. To mitigate these issues, we show the following proposition that reveals the linear relationship between the PME of $\mathbf{x}_{\text{D},n}$ under two specific prior distributions: $U_{\alpha,\bar{\mathcal{Q}}}\{\mathbf{x}_{\text{D},n}\}$ and $U_{\mathcal{Q}^{R_{\text{D}}}}\{\mathbf{x}_{\text{D},n}\}$ (refer to Definition 2 in Section II-B).

Proposition 2: Given observation vector $\hat{\mathbf{x}}=\mathbf{x}+\mathbf{e}$ with Gaussian estimation error $\mathbf{e}\sim\mathcal{CN}\left(\mathbf{0},N_{\text{e}}\mathbf{I}\right)$, there is a linear relationship between $\text{PME}(\hat{\mathbf{x}},\bar{\mathcal{S}},U_{\alpha,\bar{\mathcal{S}}}\{\mathbf{x}\},N_{\text{e}})$ and $\text{PME}(\hat{\mathbf{x}},\mathcal{S},U_{\mathcal{S}}\{\mathbf{x}\},N_{\text{e}})$, i.e.,

where $\bar{\mathcal{S}}=\{\mathcal{S},\mathbf{0}\}$, and the coefficient $C_{\text{PME}}(\hat{\mathbf{x}},\mathcal{S},\alpha,N_{\text{e}})$ is defined as

The proof is given in Appendix B.

According to Proposition 2, we can reformulate the PME of $\mathbf{x}_{\text{D},n}$ under the prior $U_{\alpha,\bar{\mathcal{Q}}}\{\mathbf{x}\}$, $\text{PME}(\hat{\mathbf{x}}_{\text{D},n}^{(k)},\bar{\mathcal{Q}},U_{\alpha,\bar{\mathcal{Q}}}\{\mathbf{x}\},N_{\text{e}}^{(k)})$, in equation (27) as the product of a coefficient and $\text{PME}(\hat{\mathbf{x}}_{\text{D},n}^{(k)},\mathcal{Q}^{R_{\text{D}}},U_{\mathcal{Q}^{R_{\text{D}}}}\{\mathbf{x}\},N_{\text{e}}^{(k)})$. The main advantage of this reformulation is that we can further decouple each element in $\text{PME}(\hat{\mathbf{x}}_{\text{D},n}^{(k)},\mathcal{Q}^{R_{\text{D}}},U_{\mathcal{Q}^{R_{\text{D}}}}\{\mathbf{x}\},N_{\text{e}}^{(k)})$ and compute them independently, which is explained by Proposition 3. The proof is given in Appendix C.

Proposition 3: Given observation vector $\hat{\mathbf{x}}=\mathbf{x}+\mathbf{e}\in\mathbb{C}^{S}$ with $\mathbf{x}\sim U_{\mathcal{S}}\{\mathbf{x}\}$ and Gaussian estimation error $\mathbf{e}\sim\mathcal{CN}\left(\mathbf{0},N_{\text{e}}\mathbf{I}_{S}\right)$, we can express the $s$-th entry of $\text{PME}(\hat{\mathbf{x}},\mathcal{S},U_{\mathcal{S}}\{\mathbf{x}\},N_{\text{e}})$ as