Covariance-Based Cooperative Activity Detection for Massive Grant-Free Random Access

For the problem in model (1), the unknown device state vectors for different APs are different. Moreover, there are some common characteristics among the neighboring APs. To enhance the detection performance, we first associate a local estimator with each AP. Incorporating the estimates of the neighboring APs, i.e., sparsity-promoting and the similarity-promoting terms [5, 11], we can modify the local estimator to associate a regularized local cost function with each AP.

Firstly, we design the local estimator. It is well known that the covariance-based massive activity detection is equivalent to recovering the device state vector $\bm{\gamma}_{b}$ from the noisy measures $\mathbf{Y}_{b}$ with the knowledge of the pre-defined pilot sequence matrix $\mathbf{S}$. In general, the estimation of the device state vector $\bm{\gamma}_{b}$ can be formulated as a maximum likelihood estimation problem [6]. In particular, for a given $\mathbf{Y}_{b}$, each column of $\mathbf{Y}_{b}$, denoted as $\mathbf{y}_{bm}$, $1\leq m\leq M$, can be termed as an independent sample having the following multivariate complex Gaussian distribution:

where the covariance matrix is calculated by $\mathbb{E}[\mathbf{y}_{bm}\mathbf{y}_{bm}^{H}]$ and $\mathbf{I}$ denotes the identity matrix. For convenience, we define $\bm{\Sigma}_{b}=\mathbf{S}\bm{\Gamma}_{b}\mathbf{S}^{H}+\sigma^{2}\mathbf{I}$. Then, the likelihood of $\mathbf{Y}_{b}$ given $\bm{\gamma}_{b}$ can be represented as

where $\det(\cdot)$ and $\text{tr}(\cdot)$ are operators that return the determinant and the trace of a matrix, respectively. By exploiting the Gaussianity, we can obtain the Maximum Likelihood (ML) estimator of $\bm{\gamma}_{b}$ at the $b$th AP as follows:

where $\hat{\bm{\Sigma}}_{b\mathbf{y}}=\frac{1}{M}\mathbf{Y}_{b}\mathbf{Y}_{b}^{H}$ denotes the sample covariance matrix of the received signal of the $b$th AP averaged over different antennas. Based on (4), the maximum likelihood estimation problem can be formulated as $\arg\min_{\bm{\gamma}_{b}\in\mathbb{R}_{+}}f(\bm{\gamma}_{b})$.

Secondly, since the activity detection is a typical sparse signal processing problem, we propose a sparsity-promoting term to facilitate cooperative detection. The specific sparsity pattern can be simultaneously observed at different APs, namely the indices of nonzero entries of $\bm{\gamma}_{b}$ are consistent for $b=1,2,\cdots,B$. Because each AP only communicates with its neighbor APs, it cannot obtain the global information about the sparsity pattern. Moreover, it is quite challenging to split this global quantity into several local quantities consisting of components only from the neighboring nodes. In this case, for the $b$th AP, we define a local parameter matrix consisting of the parameter vectors of all its neighbors, which can be directly obtained as follows:

where $l_{i}\in\mathcal{N}_{b}^{-}$ is the index set of neighbors of the $b$th AP except itself, $\mathcal{N}_{b}$ denotes the index set of the neighbors of the $b$ AP including itself, and $|\mathcal{N}_{b}^{-}|_{c}$ and $|\mathcal{N}_{b}|_{c}$ denote the cardinality of the set $\mathcal{N}_{b}^{-}$ and $\mathcal{N}_{b}$, respectively. Consequently, we aim to impose sparsity constraints on the row vectors of matrix $\mathbf{R}_{b}$ to exploit the joint sparsity. To this end, this paper designs a novel sparsity-promoting term, which is given by

where $\theta>0$ is the penalty parameter, $\mathbf{R}_{b}(n,:)$ is the $n$th row of $\mathbf{R}_{b}$, and $\left\|\cdot\right\|_{2}$ denotes the $l_{2}$ norm of a matrix. Herein, $g(\bm{\gamma}_{b})$ is the logarithmic smooth function which can promote row sparsity [5], where the nonzero rows are penalized by minimizing $g(\bm{\gamma}_{b})$. In this way, a common sparsity profile across the columns of the local parameter matrix $\mathbf{R}_{b}$ is promoted. Although the sparsity-promoting term is imposed on the local parameter matrix $\mathbf{R}_{b}$, the cooperative nature promotes a common sparsity profile across all columns of the global device state vectors $\{\bm{\gamma}_{b}\}_{b=1}^{B}$.

Thirdly, we design a similarity-promoting term to improve the detection performance. The supports of the global device state vector $\{\bm{\gamma}_{b}\}_{b=1}^{B}$ for all APs are the same, but the amplitudes of the nonzero entries at the APs are different from each other due to the effects of different path loss. In particular, the device state vectors of neighboring APs have a large number of similar entries and only a relatively small number of distinct entries. Motivated by these observations, we design a similarity-promoting function as follows

where $c_{lb}$ are linear weights satisfying the conditions: $\sum\limits_{l\in{\mathcal{N}_{b}}}{c_{lb}}=1,~{}~{}{c_{lb}}=0~{}~{}\forall l\notin{\mathcal{N}_{b}}$. $\Psi_{l}(\bm{\gamma}_{b}-\bm{\gamma}_{l})$ is a convex penalty function, minimized at $\Psi_{l}(\mathbf{0})$, which encourages similarity between $\bm{\gamma}_{b}$ and $\bm{\gamma}_{l}$. Note that the log-likelihood $f(\bm{\gamma}_{b})$ depends on the empirical covariance $\hat{\bm{\Sigma}}_{b\mathbf{y}}$. In high-dimensional settings, where the length of pilot sequences $L$ is larger than the number of AP antennas $M$, $\hat{\bm{\Sigma}}_{b\mathbf{y}}$ will be relatively different from the covariance matrix $\bm{\Sigma}_{b}$. By enforcing structural similarity, each $\bm{\Sigma}_{b}$ can exploit from the fact that neighboring AP estimates should be similar to each other.

The specific expressions in the penalty function $\Psi_{l}(\cdot)$ form can be set different. In general, it dependents on the assumptions imposed on the problem, one may choose the most appropriate penalty for the data at hand. For example, $l_{1}$-norm penalty $\Psi_{l}(\mathbf{x})=\sum_{n=1}^{N}|x_{n}|$, where $x_{n}$ and $|\cdot|$ denote the $n$th element of the vector $\mathbf{x}$ and the absolute value, respectively. This penalty function encourages the changes of limited number of values between neighbor APs, while the rest of the structure remains the same. In other words, it borrows information aggressively across neighbors, encouraging not only similar structure but also similar values. As a result, this penalty is suitable for massive access where only a small fraction of potential devices to change their states at a time slot and is adopted in this paper.

After defining the similarity-promoting term and sparsity-promoting term, accumulating them into (4) leads to the following novel regularized local cost function at the $b$th AP:

where $\beta>0$ and $\tau>0$ are the penalty parameters used to enforce sparsity and similarity, respectively. In the following, we design a massive activity detection algorithm to minimize the local cost function at each AP.

Note that the first term of (8), i.e., $f(\bm{\gamma}_{b})$ is differentiable and geodesically convex [12]. However, as stated in the above subsection, $\Psi_{l}(\bm{\gamma}_{b}-\bm{\gamma}_{l})$ is discontinuous, i.e., the third term of the local cost function could be a sum of non-smooth functions, and the second term is also potentially non-differentiable. In addition, the unknown variables $\bm{\gamma}_{b}$ for neighboring APs are coupled with each other. These obstacles make the problem intractable to solve and existing algorithms are not applicable to such a problem. In the following, we design a decentralized approximate separating strategy for minimizing the cost function in (8) based on the forward-backward splitting strategy [13], which can handle the non-smooth problem and is especially amenable to solve the high-dimensional activity detection problem due to its fast convergence rate and its conceptual and mathematical simplicity.

Before proceeding, we recall the forward-backward splitting approach for minimizing (8), which is given by the iteration

where $\bigtriangledown f(\cdot)$ is the gradient of function $f$, $\eta_{b}$ is the step size for the $b$th AP, and $\bm{\gamma}_{b}^{t}$ denotes the value of $\bm{\gamma}_{b}$ in the $t$th iteration. The gradient descent step is the forward step and the proximal step is the backward step [13]. Note that the proximal operator of a function $h$ is a mapping function given by: $\text{prox}_{\eta h}(\mathbf{y})=\arg\min_{\mathbf{u}}h(\mathbf{u})+\frac{1}{2\eta}\left\|\mathbf{u}-\mathbf{y}\right\|_{2}^{2}$ with variables $\mathbf{y}$ and $\mathbf{u}$, and a step-size $\eta>0$ [14].

Unfortunately, it is prohibitively challenging to directly evaluate the proximal operators with respect to similarity-promoting function $\Psi(\bm{\gamma}_{b})$ and the sum of $\beta g(\bm{\gamma}_{b})+\tau\Psi(\bm{\gamma}_{b})$. Moreover, the calculation of $\Psi(\bm{\gamma}_{b})$ over all the number of neighborhood, $|\mathcal{N}_{b}|_{c}$, in each iteration is expensive. Motivated by Douglas Rachford splitting in [13], where two of the proximal operators can be updated alternately, this paper aims to handle the proximal operator of function $\Psi(\bm{\gamma}_{b})$ and $g(\bm{\gamma}_{b})$ separately. Specifically, we first calculate an estimator $\mathbf{x}_{b}^{t}$ of the subgradient $\partial\Psi(\bm{\gamma}_{b}^{t})$ and then incorporate the gradient descent step into the proximal step with respect to sparsity-promoting term $g(\cdot)$ for the $b$th AP, which is given by

where $\mathbf{z}_{b}^{t}$ is a intermediate variable. Then, according to the update rule of Douglas Rachford splitting, we incorporate $\mathbf{z}_{b}^{t}$ into proximal operator with respect to similarity-promoting function $\Psi(\bm{\gamma}_{b})$:

The intermediate variable $\mathbf{z}_{b}^{t}$ and device state vector $\bm{\gamma}_{b}^{t}$ iterate alternately and their values approach to each other. When converging to optimality, their values are identical. Afterwards, in order to overcome the difficulty in processing the non-smooth finite sum term and to reduce the computational overhead, this paper proposes a splitting strategy that uses the proximal operator of a single function $\Psi_{l}$ in each iteration to approximate the proximal operator of the average of $|\mathcal{N}_{b}|_{c}$ non-smooth functions $\Psi_{l}$. In mathematical terms, we first choose $l$ randomly from the set $\mathcal{N}_{b}$ with probabilities $\{p_{1},p_{2},\cdots,p_{\left|\mathcal{N}_{b}\right|_{c}}\}$. Then, utilizing the proximal operator, a specific step is introduced as follows

where $\mathbf{x}_{b}^{l,t}$ is the estimator of subgradient $\partial\Psi_{l}(\bm{\gamma}_{b}^{t+1})$ for the randomly selected $l$th neighbor of the $b$th AP in the $t$th iteration. Let $c_{lb}^{t}$ denotes the combiner at $t$th iteration. In the sequel, $\eta_{b}^{l}$ can be set to $\eta_{b}^{l}=\frac{c_{lb}^{t}\eta_{b}}{p_{l}}$, which is a stochastic approximation of $\eta_{b}$ controlled by the combiner and the probability of being selected. In this way, we are able to treat the difficult term in (8) with non-smooth finite sum term for any size of cardinality $|\mathcal{N}_{b}|_{c}$.

Since $\mathbf{z}_{b}^{t}$ and $\bm{\gamma}_{b}^{t}$ converge to the same value, (12) is an accurate approximation of (9) if $\mathbf{x}_{b}^{l,t}=\partial\Psi_{l}(\bm{\gamma}_{b}^{t+1})$ and $\mathbf{x}_{b}^{t}=\partial\Psi(\bm{\gamma}_{b}^{t+1})$ hold. Thus, we must ensure that $\partial\Psi_{l}(\bm{\gamma}_{b}^{t+1})$ is close to $\mathbf{x}_{b}^{l,t}$ to obtain an accurate estimator. According to the definition of proximal operator, equation (12) satisfies

Hence, we can arrive at the following subgradient estimator $\mathbf{x}_{b}^{l,t+1}$ such that (12) holds:

where the right hand side (RHS) of (14) is obtained by replacing the proximal step in (13) by $\bm{\gamma}_{b}^{t+1}$ and further reorganizing the left hand side (LHS) of formula (13). Consequently, the subgradient estimator $\mathbf{x}_{b}^{t}$ in (10) can be updated as

which exploits the fact that $\mathbf{x}_{b}^{t}=\sum_{l=1}^{|\mathcal{N}_{b}|_{c}}c_{lb}^{t}\mathbf{x}_{b}^{l,t}$ and as stated in (12), only a single selected $\mathbf{x}_{b}^{l,t}$ is updated in each iteration.

Since the proximal operator needs to be calculated at each iteration in (10) and (12), it is important to derive closed form expressions for evaluating them exactly. According to the well-known Sherman-Morrison rank-1 update identity [15], we obtain

with $\bm{\Sigma}_{bn}=\bm{\Sigma}_{b}-{\gamma}_{bn}\mathbf{s}_{n}\mathbf{s}_{n}^{H}$, where ${\gamma}_{bn}$ is the $n$th element of $\bm{\gamma}_{b}$. Applying the well-known determinant identity yields

Then, substituting (III-C) and (17) into (4) and taking the derivative of $f(\bm{\gamma}_{b})$ with respect to ${\gamma}_{bn}$, we have

Correspondingly, the gradient $\bigtriangledown f(\bm{\gamma}_{b}^{t})$ can be derived by computing the following derivative: $\bigtriangledown f(\bm{\gamma}_{b}^{t})=\text{col}\{\bigtriangledown f(\gamma_{b1}^{t}),\cdots,\bigtriangledown f(\gamma_{bN}^{t})\}$, where $\text{col}(\cdot)$ denotes a column vector. By substituting it into (10) and calculating the closed form expression of the proximal operator of $g(\bm{\gamma}_{b})$, we can obtain the following intermediate recursion

with $\bm{\varsigma_{b}}^{t}=\bm{\gamma}_{b}^{t}-\eta_{b}\bigtriangledown f(\bm{\gamma}_{b}^{t})-\tau\eta_{b}\mathbf{x}_{b}^{t}$, where $\varsigma_{bn}^{t}$ is the $n$th element of $\bm{\varsigma}_{b}^{t}$.

Now, we turn to derive the recursion of $\bm{\gamma}_{b}^{t}$ in (12). Since $\Psi_{l}(\cdot)$ in (12) is fully separable, its proximal operator can be evaluated component-wise:

where the minimization operator is to preserve the positivity of $\gamma_{bn}$. Here, $z_{bn}^{t}$ and $x_{bn}^{l,t}$ are the $n$th entry of vectors $\mathbf{z}_{b}^{t}$ and $\mathbf{x}_{b}^{l,t}$, respectively. For (15) and (23), the estimation performance depends, to a great extent, on the cooperation strategy specified by the combiner $c_{lb}^{t}$. This paper adopts the following adaptive combiner

where $\rho$ is a large constant set beforehand. Note that the term $\|\bm{\gamma}_{b}^{t-1}-\bm{\gamma}_{l}^{t-1}\|_{2}$ in (27) accounts for the distance between the local estimates of the $b$th AP and its $l$th neighbor. The combiner $c_{lb}^{t}$ is inversely proportional to such a distance. When the distance defined above between two APs is large, the $b$th AP tends to decrease the combination weight, or even discard the information from this neighbor. Conversely, the $b$th AP will increase the combination weight when the distance of estimation between two APs is small. For clarity, the pseudo-code of the CMD algorithm is summarized in Algorithm 1.

Once an estimate $\{\bm{\gamma}_{b}\}_{b=1}^{B}$ is obtained, we employ the element-wise thresholding at each AP to determine $\chi_{n}$ from $\gamma_{bn}$ which is the $n$-th entry of $\bm{\gamma}_{b}$. In specific, $\chi_{n}=1$ if $\gamma_{b^{0}n}>\imath\sigma^{2}$ for a pre-specified threshold $\imath>0$, and $\chi_{n}=0$ otherwise [16]. Herein, $b^{0}$ is the AP closest to device $n$. Note that there exists a particular fixed point $\bm{\gamma}_{b}^{*}$ of the problem $\arg\min_{\bm{\gamma}_{b}\in\mathbb{R}_{+}}F(\bm{\gamma}_{b})$ for (10) and (12). Based on the information across the 6G wireless network, it can be shown that the iterates $\bm{\gamma}_{b}^{t}$ converge to this particular fixed point under certain step-size conditions.

In what follows, the computational complexity and communication cost of the proposed CMD algorithm is analyzed. In each iteration, for an arbitrary AP, the computational complexity mainly arises from the matrix multiplication, and the overall computational complexity of the CMD algorithm is $\mathcal{O}(L^{2}N)$. Although the computational complexity of sample covariance $\hat{\bm{\Sigma}}_{b\mathbf{y}}$ is $\mathcal{O}(L^{2}M)$, it only needs to be calculated once at each time slot before the iteration. For the communication cost, in each iteration, each AP needs to transmit $N$-dimensional intermediate $\bm{\gamma}_{b}^{t}$ to its neighboring APs. Thus, for all APs, the CMD algorithm needs to exchange $N\sum_{b=1}^{B}\left|\mathcal{N}_{b}^{-}\right|_{c}$ parameters. Since the APs exchange intermediate variables instead of the received signal matrix $\mathbf{Y}_{b}$, the communication cost does not grow as the number $M$ of AP antennas increases.

Remark 1: It is interesting to emphasize that the computational complexity and the communication cost at each iteration of the CMD algorithm do not grow as the number of each AP antennas, $M$, increases. Note that the CMD algorithm can also be adopted for data detection in unsourced random access [17], where an arbitrary AP only needs to know which messages are sent without identifying which message belongs to which device. Suppose that each active device has $q$ bits to send, then the total number of potential devices $N$ in the device activity detection problem is replaced by $2^{J}$, where the small size $J$ is obtained by divided $q$-bit message into $\mathcal{Z}$ blocks. The details please refer to our journal paper. Thus, the computational cost of CMD algorithm reduces to $2^{J}\sum_{b=1}^{B}\left|\mathcal{N}_{b}^{-}\right|_{c}$. This implies that the communication cost of the proposed CMD algorithm does not grow by increasing the total number of potential devices, which is an appealing feature for reality IoT networks.