# Cooperative Activity Detection: Sourced and Unsourced Massive Random Access
Paradigms

Xiaodan Shao, Xiaoming Chen, Derrick Wing Kwan Ng, Caijun Zhong, and Zhaoyang
Zhang Part of this paper has been accepted for presentation at the IEEE Global
Communications Conference (GLOBECOM), Dec. 2020 [1].Xiaodan Shao
(shaoxiaodan@zju.edu.cn), Xiaoming Chen (chen_xiaoming@zju.edu.cn), Caijun
Zhong (caijunzhong@zju.edu.cn), and Zhaoyang Zhang<EMAIL_ADDRESS>are
with the College of Information Science and Electronic Engineering, Zhejiang
University, Hangzhou 310016, China. Derrick Wing Kwan Ng<EMAIL_ADDRESS>is with the School of Electrical Engineering and Telecommunications,
University of New South Wales, Sydney, NSW 2052, Australia.

###### Abstract

This paper investigates the issue of cooperative activity detection for grant-
free random access in the sixth-generation (6G) cell-free wireless networks
with sourced and unsourced paradigms. First, we propose a cooperative
framework for solving the problem of device activity detection in sourced
random access. In particular, multiple access points (APs) cooperatively
detect the device activity via exchanging low-dimensional intermediate
information with their neighbors. This is enabled by the proposed covariance-
based algorithm via exploiting both the sparsity-promoting and similarity-
promoting terms of the device state vectors among neighboring APs. A
decentralized approximate separating approach is introduced based on the
forward-backward splitting strategy for addressing the formulated problem.
Then, the proposed activity detection algorithm is adopted as a decoder of
cooperative unsourced random access, where the multiple APs cooperatively
detect the list of transmitted messages regardless of the identity of the
transmitting devices. Finally, we provide sufficient conditions on the step
sizes that ensure the convergence of the proposed algorithm in the sense of
Bregman divergence. Simulation results show that the proposed algorithm is
efficient for addressing both sourced and unsourced massive random access
problems, while requires a shorter signature sequence and accommodates a
significantly larger number of active devices with a reasonable antenna array
size, compared with the state-of-art algorithms.

###### Index Terms:

Cooperative activity detection, sourced random access, unsourced random
access, 6G cell-free wireless networks, covariance-based detection.

##  I Introduction

Massive machine-type communication (mMTC) or massive access, which is expected
to be a typical application scenario for 6G wireless networks, aims to meet
the demand of massive connectivity for the Internet-of-Things (IoT) [2]-[4].
Unfortunately, applying conventional grant-based random access schemes for
massive access leads to an exceedingly long access latency and a prohibitive
signaling overhead. As a remedy, grant-free random access schemes have been
proposed and considered as a promising technical candidate for realizing
massive access [5]-[7], where active devices transmit their data signals
without obtaining a grant from their home base stations (BSs) after sending
pre-assigned signature sequences.

In general, grant-free random access includes two different paradigms, namely
sourced and unsourced random access [8]-[21]. For sourced random access, the
BS is interested in both the messages and the identities of the devices that
generated them. In particular, sourced random access is mainly applied to
device-oriented applications. For example, in the application of status
monitoring, the BS is required to know which sensors send the message and
updates their status periodically. Hence, each device is preassigned to a
unique signature sequence and the BS identifies the device activity based on
the received signals by detecting which sequences are transmitted. Since the
active device detection for sourced random access is a typical sparse signal
processing problem, numerous compressed sensing (CS)-based approaches have
been employed to handle the detection problems. For instance, in [8]-[10],
approximate message propagation (AMP) algorithms were designed for joint
activity detection and channel estimation (JADCE) in different scenarios by
exploiting the statistics of wireless channels. Then, the authors in [11]
proposed an orthogonal AMP-multiple measurement vectors (OAMP-MMV) algorithm
that applies Hadamard pilot matrice and reduces the computational complexity
by using fast Fourier transform. In addition, the authors in [12] proposed a
low-complexity dimension reduction-based JADCE algorithm, which projects the
original device state matrix to a low-dimensional space by exploiting its
sparse and low-rank structure. Note that the aforementioned approaches in
[8]-[12] performed activity detection based on the instantaneous received
signals, which require exceedingly long signature sequences in the scenario of
massive access. To tackle this problem, the covariance-based algorithms were
proposed to improve the performance of device activity detection in [13] and
[15], where the detection problem was handled by a coordinate descent
algorithm. After that, a low computational complexity covariance-based
algorithm with constant-modulus pilots was proposed for JADCE, which shows low
computational complexity and appealing activity detection performance [16]. In
fact, the covariance-based algorithm can outperform the AMP algorithm with the
same length of signature. However, the superior system performance comes from
the expense of the use of a relatively large number of receive antennas
compared with what traditional CS-based algorithms needed. As a result, active
device detection in massive access has emerged as a challenging problem due to
a large number of devices and the limited radio resources in $6$G wireless
networks.

As for unsourced random access, the AP is interested in the transmitted
messages only and not the identity of the devices. In fact, unsourced random
access is motivated by practical IoT scenarios, where millions of low-cost
devices have their codebook hardwired at the moment of production [18]. Such
that all the devices share a common codebook in unsourced random access. Since
the same codebook is exploited, the AP can only decode the list of transmitted
messages irrespectively of the identity of the active devices. In this
context, the computational complexity can be significantly reduced compared
with the case that assigns a unique individual codebook to each device [5].
Therefore, unsourced random access is mainly applied to content-oriented
applications. For example, in quality inspections process of smart factories,
a fraction of devices send their current reliability information to the AP,
while the AP detects all received messages, computes a weighted average of
current informations, then generates a performance index. Recently, a low-
complexity unsourced random access algorithm based on the coupled CS framework
was proposed in [19]. Specifically, the transmission slot is partitioned into
subslots and each active device sends a codeword from a common codebook across
different subslots. To this connection, an inner encoder is needed which maps
each submessage into one column of a given coding matrix. Then, at the BS, the
inner decoder must identify which columns of the matrix have been transmitted.
On the other hand, an outer tree-based decoder is applied to stitch the
decoded sequences. However, this work assumed only a single-antenna receiver
equipped at the BS and the results are not applicable to the case of multiple
antennas. Afterward, the authors in [21] extended the model in [19] to a case
of large-scale receive antenna arrays with Rayleigh fading, where a maximum
likelihood (ML)-based activity detection scheme was adopted as the inner
decoder of a concatenated coding scheme in [18]. Specifically, the algorithm
in [21] avoided the use of a signature sequence longer than the number of
active devices. However, it usually requires a large number of antennas at the
BS in order to accumulate more active devices [20]. Moreover, for a short
subslot length, the number of active devices that can be accommodated is
limited by the detection capabilities of the inner decoder [21].

To overcome these challenges in sourced and unsourced random access,
cooperative activity detection with multiple APs can be applied to detect the
active device and the active columns of the coding matrix. Inspired by this,
the authors in [22] considered the active device detection in multi-cell
massive multiple-input-multiple-output (MIMO) systems and cooperative MIMO
systems using the AMP algorithm. For multi-cell massive MIMO, the authors
assumed that each AP operates independently to detect the active devices of
its own cell. In this case, inter-cell interference is treated as noise, which
is a severe limiting factor for achieving reliable activity detection. As for
cooperative MIMO, each AP performs active device detection locally and detects
the devices from neighboring cells as well. Then, the detection results in the
form of log-likelihood ratio are forwarded to a central unit where final
decisions on device activity are carried out. Then, the authors in [23]
proposed an active device identification method in a multi-cell network
adopting second-order Reed-Muller sequences, where each AP identifies all in-
cell devices and the transmissions from out-of-the-cell devices are regarded
as interference. Compared with [22] and [23], the proposed cell-free framework
in this paper can avoid the inter-cell interference due to the fact that the
concepts of cell and cell boundary do not exist in our system.

Motivated by these works, this paper designs a unified framework for both
sourced and unsourced activity detection in 6G cell-free wireless networks
[24]-[26], where multiple APs are deployed in a vast area to serve all devices
located in this area via a fronthaul network [27]. For instance, the authors
in [26] studied grant-free massive access in cell-free massive MIMO based IoT,
where multiple APs cooperate in the network to serve massive devices. Herein,
the connections between APs can be set differently according to the
communication radius. Unlike the conventional cooperative MIMO where the final
decisions on device activity are only carried out in a central unit, the
proposed framework does not need a centralized fusion center and the
cooperative activity detection among the APs only needs to exchange low-
dimensional intermediate information, i.e., device state vectors, thereby all
the APs in the system can obtain the detection results. Therefore, it is
reliable and robust to AP and/or fronthaul link failure. More importantly, to
reduce the amount of associated signaling overhead under this framework, this
paper designs a scalable computationally efficient algorithm to detect the
activity. In summary, our main contributions are listed as follows:

  1. 1.

We propose a unified cooperative activity detection framework for sourced and
unsourced random access based on the covariance of the received signals in 6G
cell-free wireless networks to support massive IoT with heterogeneous
application requirements.

  2. 2.

We develop a novel low-complexity cooperative activity detection (CAD)
algorithm that exploits the special characteristic of the device state vectors
of interest among neighboring APs, namely joint similarity and sparsity.

  3. 3.

We provide a theoretical analysis of the convergence property of the proposed
CAD algorithm and the result shows that it enjoys a convergence rate of
$\mathcal{O}(1/t)$ under certain general conditions with an iteration index
$t$. Moreover, extensive simulation results confirm the effectiveness of the
proposed CAD algorithm in both sourced and unsourced random access paradigms.

The rest of this paper is organized as follows. Section II gives a brief
introduction of 6G cell-free wireless networks in a sporadic device activity
pattern scenario. Then, Section III proposes the cooperative sourced random
access scheme. Section IV designs the cooperative unsourced random access
scheme. Next, Section V analyzes the performance of the proposed CAD
algorithm. Afterward, Section VI provides extensive simulation results to
illustrate the performance of the proposed algorithm. Finally, Section VII
concludes the paper.

_Notations_ : We use bold letters to denote matrices or vectors, non-bold
letters to denote scalars, $\mathbf{I}$ to denote the identity matrix,
$\mathbb{C}^{A\times B}$ to denote the space of complex matrices of size
$A\times B$, $|\cdot|$ to denote the absolute value of a complex number,
$(\cdot)^{H}$ and $(\cdot)^{T}$ to denote conjugate transpose and transpose
respectively, $\left\|\cdot\right\|_{F}$ to denote Frobenius norm of a matrix,
$\left\|\cdot\right\|_{2}$ to denote the $l_{2}$-norm of an input vector.
$\|\cdot\|_{0}$ denotes the $l_{0}$-norm defined as the number of nonzero
elements of an input vector. $\mathbf{A}(n,:)$ denotes the $n$th row of matrix
$\mathbf{A}$. $\|\mathbf{A}\|_{2,0}$ denotes the $l_{20}$-norm defined as the
number of nonzero elements of vector
$\left[\|\mathbf{A}(1,:)\|_{2},\cdots,\|\mathbf{A}(n,:)\|_{2}\right]$.
$\det(\cdot)$ and $\text{tr}(\cdot)$ are operators that return the determinant
and the trace of an input matrix, respectively.
$\left\langle\mathbf{a},\mathbf{b}\right\rangle$ denotes the inner product of
the vectors $\mathbf{a}$ and $\mathbf{b}$. $\text{col}(\cdot)$ denotes a
column vector. $|\cdot|_{c}$ denotes the cardinality of a set. $\odot$ denotes
element wise multiplication. $\bigtriangledown f(\cdot)$ denotes the gradient
of a function $f(\cdot)$. $\partial f(\cdot)$ denotes the subgradient of a
function $f(\cdot)$. $\mathcal{O}(\cdot)$ stands for the big-O notation.

##  II System Model

Consider a 6G cell-free wireless network comprising $B$ APs. The APs are
equipped with $M$ antennas each, serving $N$ uniformly distributed single-
antenna IoT devices in a vast area. Each AP is connected to several adjacent
APs via fronthaul links and can only communicate with its one-hop neighbors
for reducing the communication load, as shown in Fig. 1. In 6G wireless
networks, the density of IoT devices is usually huge, e.g., $10$ devices per
$m^{2}$. However, due to the bursty communication characteristic of IoT
applications, only a fraction of IoT devices are active at any given time
slot. We adopt $\mathcal{K}$ to denote the set of active devices with
cardinality $K=\left|\mathcal{K}\right|_{c}\ll N$ being the number of the
active devices, which is a random variable. For convenience, we define
$\chi_{n}$ as a binary activity indicator with ${\chi_{n}}=1$ if the $n$th
device is active, and ${\chi_{n}}=0$ otherwise. Moreover, we represent the
$M$-dimensional channel vector from the $n$th device to the $b$th AP as
$\sqrt{g_{b,n}}\mathbf{h}_{b,n}$, where $g_{b,n}$ is the large-scale fading
depending on the location of the $n$th device and
$\mathbf{h}_{b,n}\in\mathbb{C}^{M}$ is the corresponding small-scale fading
following independent and identically distributed (i.i.d.) complex Gaussian
distribution with zero mean and unit variance.

Figure 1: Illustration of a 6G cell-free wireless network with multiple APs.

In order to reduce the access latency and the system signaling overhead, a
grant-free random access protocol is advocated for 6G cell-free wireless
networks [5]. Specifically, at the beginning of each time slot, the active
devices transmit the corresponding signature sequences over uplink channels
simultaneously and then the APs perform the activity detection based on the
received signal cooperatively. Herein, it is assumed that all signature
sequences $\mathbf{s}_{n}\in\mathbb{C}^{L}$ are generated following i.i.d.
complex Gaussian distribution with zero mean and unit variance which are known
at the APs in advance. Thus, the received signal
$\mathbf{Y}_{b}\in\mathbb{C}^{L\times M}$ at the $b$th AP can be expressed as

$\displaystyle\mathbf{Y}_{b}=\sum_{n=1}^{N}\chi_{n}\mathbf{s}_{n}\sqrt{g_{b,n}}\mathbf{h}_{b,n}^{T}+\mathbf{W}_{b}=\mathbf{S}\bm{\Gamma}_{b}^{\frac{1}{2}}\mathbf{H}_{b}+\mathbf{W}_{b},$
(1)

where
$\mathbf{H}_{b}=[\mathbf{h}_{b,1},\cdots,\mathbf{h}_{b,N}]^{T}\in\mathbb{C}^{N\times
M}$ denotes the small-scale fading channel matrix,
$\mathbf{S}=[\mathbf{s}_{1},\cdots,\mathbf{s}_{N}]\in\mathbb{C}^{L\times N}$
denotes the horizontal stack of all signature sequences, and
$\mathbf{W}_{b}\in\mathbb{C}^{L\times M}$ is the additive white Gaussian noise
(AWGN) marix with i.i.d. entries $\thicksim\mathcal{CN}(0,\sigma^{2})$, where
$\sigma^{2}$ denotes the noise power at each antenna. Define
$\bm{\gamma}_{b}=[\gamma_{b,1},\cdots,\gamma_{b,N}]^{T}\in\mathbb{R}^{N}$ as
the diagonal entries of diagonal matrix $\bm{\Gamma}_{b}$, representing the
device state vector of the $b$th AP with $\gamma_{b,n}=\chi_{n}g_{b,n}$.

In this paper, we aim to provide two grant-free random access schemes
according to the characteristics and requirements of IoT applications. The
first kind of IoT applications have to know which device sends the messages.
In this case, we should detect the identity of active device, namely sourced
random access. The second kind of IoT applications only need to know what
messages are sent. In this case, we directly detect the messages, namely
unsourced random access. In the following, we will design sourced and
unsourced random access schemes via exchanging some intermediate variables
based on a limited cooperation among multiple APs.

##  III Cooperative Sourced Random Access

For sourced random access, each device is assigned to an unique sequence.
Based on the received sequences, the APs perform the cooperative device
activity detection. In this section, we first propose a cooperative detection
framework for 6G cell-free wireless networks with a massive number of IoT
devices. Then, we design a corresponding cooperative detection algorithm for
such a sourced random access scenario.

###  III-A Cooperative Massive Detection Framework

For the activity detection problem based on model in (1), the unknown device
state vectors for different APs are generally different. Moreover, there are
some common characteristics among the neighboring APs. To enhance the
detection performance, we first associate a local estimator for each AP due to
the fact that different AP estimates a different device state vector. Then, to
incorporate the estimates of neighboring APs, i.e., sparsity-promoting and the
similarity-promoting terms [12, 34], we can modify the local estimator to
associate a regularized local cost function with each AP.

####  III-A1 Covariance Based Local Estimator

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 sequence matrix $\mathbf{S}$. In general, the
estimation of the device state vector $\bm{\gamma}_{b}$ can be formulated as a
ML estimation problem [13]. 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:

$\mathbf{y}_{bm}\sim\mathcal{CN}(\mathbf{0},\mathbf{S}\bm{\Gamma}_{b}\mathbf{S}^{H}+\sigma^{2}\mathbf{I}),$
(2)

where covariance matrix is calculated by
$\mathbb{E}\left\\{\mathbf{y}_{bm}\mathbf{y}_{bm}^{H}\right\\}$. 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

$\displaystyle
P(\mathbf{Y}_{b}|\bm{\gamma}_{b})=\frac{1}{\det(\pi\bm{\Sigma}_{b})^{M}}\exp(-\text{tr}(\bm{\Sigma}_{b}^{-1}\mathbf{Y}_{b}\mathbf{Y}_{b}^{H})).$
(3)

By exploiting the Gaussianity, we can obtain the ML estimator of
$\bm{\gamma}_{b}$ at the $b$th AP as follows:

$\displaystyle
f(\bm{\gamma}_{b})=-P(\mathbf{Y}_{b}|\bm{\gamma}_{b})=\ln\text{det}(\bm{\Sigma}_{b})+\text{tr}(\bm{\Sigma}_{b}^{-1}\hat{\bm{\Sigma}}_{b\mathbf{y}}),$
(4)

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 ML estimation problem can
be formulated as
$\arg\min_{\bm{\gamma}_{b}\in\mathbb{R}_{+}}f(\bm{\gamma}_{b})$.

####  III-A2 Sparsity-Promoting Term Design

Secondly, since the focused device state vectors have inherently structured
sparsity, we propose a sparsity-promoting term to facilitate cooperative
detection among multiple APs. In [12], the rows of the interested matrix share
the common support, thus appropriate co-regularizer was added to the mean-
square error criterion to promote row sparsity. In this paper, 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 neighboring APs,
it cannot obtain the global information of the system-level 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 one-hop neighbors, which can be directly
obtained as follows:

$\displaystyle\mathbf{R}_{b}\\!=\\!\left[\bm{\gamma}_{l_{1}},\bm{\gamma}_{l_{2}},\bm{\gamma}_{l_{i}},\cdots,\bm{\gamma}_{l_{|\mathcal{N}_{b}^{-}|_{c}}},\bm{\gamma}_{b}\right]$
$\displaystyle\in\mathbb{C}^{N\times(|\mathcal{N}_{b}|_{c})},$
$\displaystyle\forall b\in\\{1,2,\cdots,B\\},$ (5)

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$th AP including itself. Consequently, we aim to promote row sparsity of
matrix $\mathbf{R}_{b}$ to exploit the joint sparsity. In general, the
$l_{20}$-norm which is a combination of $l_{0}$-norm and $l_{2}$-norm can
promote the row sparsity of an input matrix [37]. However, the $l_{20}$-norm
is non-convex and non-differentiable, and its minimization is non-
deterministic polynomial-time hardness. To this end, this paper adopts the
following sparsity-promoting term proposed in our work [12]

$g(\bm{\gamma}_{b})=\sum_{n=1}^{N}\underbrace{\left(\left\|\mathbf{R}_{b}(n,:)\right\|_{2}-\frac{1}{\theta}\ln(1+\theta\left\|\mathbf{R}_{b}(n,:)\right\|_{2})\right)}_{\Phi(\mathbf{R}_{b}(n,:))},$
(6)

Calculating the Maclaurin series of
$\ln(1+\theta\left\|\mathbf{R}_{b}(n,:)\right\|_{2})$, we note that when
$\|\mathbf{R}_{b}(n,:)\|_{2}\rightarrow 0$,
$\Phi(\mathbf{R}_{b}(n,:))\rightarrow\frac{\theta}{2}\|\mathbf{R}_{b}(n,:)\|_{2}^{2}$,
however, when $\|\mathbf{R}_{b}(n,:)\|_{2}\rightarrow\infty$,
$\Phi(\mathbf{R}_{b}(n,:))\rightarrow\|\mathbf{R}_{b}(n,:)\|_{2}$. Herein,
$\theta>0$ is the penalty parameter used to guarantee that the curve of
$\Phi(\mathbf{R}_{b}(n,:))$ is as close as possible to the curve of
$\|\mathbf{R}_{b}(n,:)\|_{2}$. In other words, the function
$\Phi(\mathbf{R}_{b}(n,:))$ is used to approximates the $l_{2}$-norm of
$\mathbf{R}_{b}(n,:)$, since $\Phi(\mathbf{R}_{b}(n,:))$ is differentiable
with respect to the row vector $\mathbf{R}_{b}(n,:)$ at $\mathbf{0}$.
Moreover, the $g(\bm{\gamma}_{b})$ is convex since it relaxes the non-convex
$l_{0}$-norm by the summation operator. 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}$. Note that depends on the assumptions
imposed on the problem, one may design other appropriate penalties for
sparsity-promoting term containing smooth or non-smooth functions.

####  III-A3 Similarity-Promoting Term Design

Thirdly, we design a similarity-promoting term to improve the detection
performance. In [28] and [34], to promote similarity among optimum parameter
vectors, a $l_{2}$-norm based co-regularizer has been introduced to the
diffusion least-mean-square algorithm. In the considered problem, the device
state vectors of interest among neighboring APs also exhibit similarity.
Specifically, the indices of nonzero entries of the global device state vector
$\\{\bm{\gamma}_{b}\\}_{b=1}^{B}$ for all APs should be the same, but the
amplitudes of the nonzero entries at the APs are different from each other due
to the effects of different path losses. In particular, since the neighboring
APs have short distances, the corresponding device state vectors consist of
some similar nonzero entries. Moreover, the ML estimator $f(\bm{\gamma}_{b})$
depends on the empirical covariance $\hat{\bm{\Sigma}}_{b\mathbf{y}}$. In
high-dimensional settings, where the length of 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, the update of each $\bm{\Sigma}_{b}$ can
exploit from the fact that the estimates of neighboring APs should be similar
to each other. Motivated by these observations, we design a similarity-
promoting function as follows

$\Psi(\bm{\gamma}_{b})=\sum_{l\in\mathcal{N}_{b}}c_{lb}\Psi_{l}(\bm{\gamma}_{b}-\bm{\gamma}_{l}),~{}\forall
b\in\\{1,2,\cdots,B\\},$ (7)

where $\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}$. In this paper, the specific
expression in the penalty function $\Psi_{l}(\bm{\gamma}_{b}-\bm{\gamma}_{l})$
form is set to $l_{1}$-norm penalty
$\Psi_{l}(\mathbf{x})=\sum_{n=1}^{N}|x_{n}|$, where $x_{n}$ denotes the $n$th
element of the vector $\mathbf{x}$. This penalty function can encourage a
large number of elements of device state vector to be identical across APs. In
other words, it borrows information aggressively across neighbors, encouraging
not only similar structure but also similar values. The adopted penalty is
suitable for massive access where only a small fraction of potential devices
are active at a time slot. Herein, $c_{lb}$ are linear weights satisfying the
conditions:

$\sum\limits_{l\in{\mathcal{N}_{b}}}{c_{lb}}=1,~{}\textrm{and}~{}~{}{c_{lb}}=0~{}~{}\forall
l\notin{\mathcal{N}_{b}}.$ (8)

After defining the similarity-promoting term and sparsity-promoting term, the
amalgamation of (4), (6), and (7) leads to the following novel regularized
local cost function to be adopted at the $b$th AP:

$\displaystyle F(\bm{\gamma}_{b})\\!=\\!f(\bm{\gamma}_{b})+\beta
g(\bm{\gamma}_{b})+\tau\Psi(\bm{\gamma}_{b}),~{}\forall
b\\!\in\\!\\{1,2,\cdots,B\\},$ (9)

where $\beta>0$ and $\tau>0$ are the penalty parameters used to enforce
sparsity and similarity, respectively [29]. By varying the values of these
penalty parameters, one can strike a balance between the ML estimator and
promoting terms. In (9), in addition to exploiting the local received signal
through certain covariance matrix, the structured sparsity and similarity
pattern observed in neighboring APs are leveraged to improve activity
detection performance. To ensure that each AP can solve the problem in an
autonomous and adaptive manner using only local interactions, in Section III.
B and Section III. C, we will design a cooperative activity detection (CAD)
algorithm to jointly minimize the costs involved (9) at each AP.

###  III-B A Decentralized Approximate Separating Strategy

Note that the first term of (9), i.e., $f(\bm{\gamma}_{b})$ is differentiable
and geodesically convex, which is a generalized form of classical convexity
and can guarantee that all local minima of geodesic convexity functions are
globally minimum [30]. 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, e.g., [31] and [34], are not applicable to such a problem. In this
paper, we aim to design a scheme which computes device state vectors of all
APs independently to facilitate decentralized implementation. In the
following, we design a decentralized approximate separating strategy for
minimizing the cost function in (9) based on the forward-backward splitting
strategy [32], 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 (9), which is given by the iteration

$\displaystyle\bm{\gamma}_{b}^{t+1}$ $\displaystyle=$
$\displaystyle\bm{\gamma}_{b}^{t}-\eta_{b}^{t}\bigtriangledown
f(\bm{\gamma}_{b}^{t})-\tau\eta_{b}^{t}\partial\Psi(\bm{\gamma}_{b}^{t+1})-\beta\eta_{b}^{t}\partial
g(\bm{\gamma}_{b}^{t+1})$ (10) $\displaystyle=$
$\displaystyle\overbrace{\text{prox}_{\eta_{b}^{t}(\tau\Psi+\beta
g)}(\underbrace{\bm{\gamma}_{b}^{t}-\eta_{b}^{t}\bigtriangledown
f(\bm{\gamma}_{b}^{t})}_{\textrm{Forward step}})}^{\textrm{Backward step}},$

where $\eta_{b}^{t}$ is the step size for the $b$th AP at iteration $t$, and
$\bm{\gamma}_{b}^{t}$ denotes the value of $\bm{\gamma}_{b}$ in the $t$th
iteration. Note that $\Psi(\bm{\gamma}_{b}^{t+1})$ and $\partial
g(\bm{\gamma}_{b}^{t+1})$ have the effect of evaluating the gradient at
$\bm{\gamma}^{t+1}$. This form is known to have better approximation
properties than conventional gradient descent algorithm which evaluates the
gradient at $\bm{\gamma}^{t}$ [32]. However, it cannot be rewritten as an
iteration that gives $\bm{\gamma}_{b}^{t+1}$ explicitly in terms of
$\bm{\gamma}_{b}^{t}$. For this reason, the proximal operator of a function
$h$ is introduced, which is 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},$
(11)

with variables $\mathbf{y}$ and $\mathbf{u}$, and a step-size $\eta>0$ [39].
As shown in the last equation of (10), the gradient descent step is the
forward step and the proximal step is the backward step. Herein, the proximal
step can be viewed as an implicit discretization of the non-differentiable
function, which can be applied to overcome the non-smoothness of
$\Psi(\bm{\gamma}_{b})$ and the proximal algorithm has a faster convergence
rate than the approach based on the subgradient strategy [32].

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 [32], 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

$\displaystyle\mathbf{z}_{b}^{t}$ $\displaystyle=$
$\displaystyle\bm{\gamma}_{b}^{t}-\eta_{b}^{t}\bigtriangledown
f(\bm{\gamma}_{b}^{t})-\tau\eta_{b}^{t}\mathbf{x}_{b}^{t}-\beta\eta_{b}^{t}\partial
g(\mathbf{z}_{b}^{t})$ (12) $\displaystyle=$
$\displaystyle\text{prox}_{\beta\eta_{b}^{t}g}(\bm{\gamma}_{b}^{t}-\eta_{b}^{t}\bigtriangledown
f(\bm{\gamma}_{b}^{t})-\tau\eta_{b}^{t}\mathbf{x}_{b}^{t}),$

where $\mathbf{z}_{b}^{t}$ is an 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})$:

$\displaystyle\bm{\gamma}_{b}^{t+1}=\text{prox}_{\tau\eta_{b}^{l,t}\Psi}(\mathbf{z}_{b}^{t}+\tau\eta_{b}^{t}\mathbf{x}_{b}^{t}).$
(13)

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. In order to
overcome the difficulty in processing the non-smooth finite sum term and to
reduce the computational overhead of (13), this paper proposes a decentralized
separating strategy that approximately separates the proximal operator of the
similarity-term of the $b$th AP, $\Psi(\bm{\gamma}_{b})$, into the proximal
operator of one of the neighbors of AP $b$,
$\Psi_{l}(\bm{\gamma}_{b}-\bm{\gamma}_{l})$, in each iteration. 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, we
incorporate $\mathbf{z}_{b}^{t}$ into proximal operator with respect to the
sub-function of similarity-promoting
$\Psi_{l}(\bm{\gamma}_{b}-\bm{\gamma}_{l})$:

$\displaystyle\bm{\gamma}_{b}^{t+1}$ $\displaystyle=$
$\displaystyle\mathbf{z}_{b}^{t}-\tau\eta_{b}^{l,t}\left[\partial\Psi_{l}(\bm{\gamma}_{b}^{t+1})-\mathbf{x}_{b}^{l,t}\right]$
(14) $\displaystyle=$
$\displaystyle\text{prox}_{\tau\eta_{b}^{l,t}\Psi_{l}}(\mathbf{z}_{b}^{t}+\tau\eta_{b}^{l,t}\mathbf{x}_{b}^{l,t}),$

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}$ denote the
combiner at $t$th iteration and set
$\eta_{b}^{l,t}=\frac{c_{lb}^{t}\eta_{b}^{t}}{p_{l}}$, which is a stochastic
approximation of $\eta_{b}^{t}$ controlled by the combiner and the probability
of being selected. In this manner, we are able to treat the difficult terms in
(9) 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, (14) is an accurate approximation of (10) 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 in (11), equation (14) satisfies

$\displaystyle\frac{(\mathbf{z}_{b}^{t}\\!+\\!\tau\eta_{b}^{l,t}\mathbf{x}_{b}^{l,t}\\!-\\!\text{prox}_{\tau\eta_{b}^{l,t}\Psi_{l}}(\mathbf{z}_{b}^{t}\\!+\\!\tau\eta_{b}^{l,t}\mathbf{x}_{b}^{l,t}))}{\tau\eta_{b}^{l,t}}\in\partial\Psi_{l}(\bm{\gamma}_{b}^{t+1}).$
(15)

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

$\mathbf{x}_{b}^{l,t+1}=\mathbf{x}_{b}^{l,t}+\frac{1}{\tau\eta_{b}^{l,t}}(\mathbf{z}_{b}^{t}-\bm{\gamma}_{b}^{t+1}),$
(16)

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

$\mathbf{x}_{b}^{t+1}=\mathbf{x}_{b}^{t}+c_{lb}^{t}(\mathbf{x}_{b}^{l,t+1}-\mathbf{x}_{b}^{l,t}),$
(17)

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 (14), only a single selected $\mathbf{x}_{b}^{l,t}$ is
updated in each iteration $t$.

###  III-C Derivation of Algorithm Recursions

Since the proximal operator needs to be calculated at each iteration in (12)
and (14), it is important to derive closed-form expressions for evaluating
$\bm{\gamma}_{b}^{t+1}$ and $\mathbf{z}_{b}^{t}$ exactly. We start by
calculating the gradient of local estimator $f(\bm{\gamma}_{b})$ in (4).
According to the well-known Sherman-Morrison rank-1 update identity [33], we
obtain

$\displaystyle\left(\bm{\Sigma}_{bn}+{\gamma}_{bn}\mathbf{s}_{n}\mathbf{s}_{n}^{H}\right)^{-1}=\bm{\Sigma}_{bn}^{-1}-\frac{{\gamma}_{bn}\bm{\Sigma}_{bn}^{-1}\mathbf{s}_{n}\mathbf{s}_{n}^{H}\bm{\Sigma}_{bn}^{-1}}{1+{\gamma}_{bn}\mathbf{s}_{n}^{H}\bm{\Sigma}_{bn}^{-1}\mathbf{s}_{n}},$
(18)

with

$\bm{\Sigma}_{bn}=\bm{\Sigma}_{b}-{\gamma}_{bn}\mathbf{s}_{n}\mathbf{s}_{n}^{H},$
(19)

where ${\gamma}_{bn}$ is the $n$th element of $\bm{\gamma}_{b}$. Applying the
well-known determinant identity111Determinant identity [35] states that if
$\mathbf{A}$ and $\mathbf{B}$ are matrices of sizes $m\times n$ and $n\times
m$, then
$\text{det}(\mathbf{I}+\mathbf{A}\mathbf{B})=\text{det}(\mathbf{I}+\mathbf{B}\mathbf{A})$.
to $\bm{\Sigma}_{bn}+{\gamma}_{bn}\mathbf{s}_{n}\mathbf{s}_{n}^{H}$ yields

$\displaystyle\text{det}(\bm{\Sigma}_{bn}+{\gamma}_{bn}\mathbf{s}_{n}\mathbf{s}_{n}^{H})=(1+{\gamma}_{bn}\mathbf{s}_{n}^{H}\bm{\Sigma}_{bn}^{-1}\mathbf{s}_{n})\text{det}(\bm{\Sigma}_{bn}).$
(20)

Then, substituting (18) and (20) into (4) and taking the derivative of
$f(\bm{\gamma}_{b})$ with respect to ${\gamma}_{bn}$ leads to

$\displaystyle\bigtriangledown
f(\gamma_{bn})=\frac{\mathbf{s}_{n}^{H}\bm{\Sigma}_{bn}^{-1}\mathbf{s}_{n}}{1+{\gamma}_{bn}\mathbf{s}_{n}^{H}\bm{\Sigma}_{bn}^{-1}\mathbf{s}_{n}}-\frac{\mathbf{s}_{n}^{H}\bm{\Sigma}_{bn}^{-1}\hat{\bm{\Sigma}}_{b\mathbf{y}}\bm{\Sigma}_{bn}^{-1}\mathbf{s}_{n}}{(1+{\gamma}_{bn}\mathbf{s}_{n}^{H}\bm{\Sigma}_{bn}^{-1}\mathbf{s}_{n})^{2}}.$
(21)

Correspondingly, the gradient $\bigtriangledown f(\bm{\gamma}_{b}^{t})$ can be
derived as $\bigtriangledown
f(\bm{\gamma}_{b}^{t})=\text{col}\\{\bigtriangledown
f(\gamma_{b1}^{t}),\cdots,\bigtriangledown f(\gamma_{bN}^{t})\\}$. Observe
that the full gradient of $f(\bm{\gamma}_{b})$ over $\bm{\gamma}_{b}$ is
calculated in each iteration, resulting in high computational complexity. To
perform as few gradient calculations as possible for striking a trade-off
between algorithm exploration and accuracy, this paper adopts a non-dense
gradient update rule where the full gradient is computed with a certain
probability $0<\bar{p}\leq 1$ resulting in a cheap per-iteration cost.
Specifically, we construct a function of the current update in the $t$
iteration for computing the gradient as follow

$\displaystyle\widetilde{\bigtriangledown
f}(\bm{\gamma}_{b}^{t})=\left\\{\begin{array}[]{l}\bigtriangledown
f(\bm{\gamma}_{b}^{t}),~{}\text{with probability}~{}\bar{p},\\\
\mathbf{0},~{}~{}~{}~{}\text{with probability}~{}1-\bar{p},\end{array}\right.$
(24)

Indeed, since each AP updates the gradient with probability $\bar{p}$ at each
iteration, all these updates add up to a dense one.

By substituting (24) into (12) and calculating the proximal step of
$g(\bm{\gamma}_{b})$, we can obtain the following intermediate recursion

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mathbf{z}_{b}^{t}=\bm{\varsigma}_{b}^{t}-\eta_{b}^{t}\beta\text{col}\left\\{\frac{\varsigma_{b1}^{t}}{\left\|\mathbf{R}_{b}^{t}(1,:)\right\|_{2}},\cdots,\frac{\varsigma_{bN}^{t}}{\left\|\mathbf{R}_{b}^{t}(N,:)\right\|_{2}}\right\\},$
(25)

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

Note that the step size $\eta_{b}^{t}$ in (25) is an essential hyper-
parameter, which should be designed for improving the detection performance.
In particular, a large step size leads to a faster convergence initially but
early saturation, while small step size slows the convergence, but postpone
the saturation. In fact, it is a difficult task in practice as one approximate
step size is determined by the Lipschitz constant which is hard to estimate.
The authors in [36] proposed the adaptive step size for stochastic gradient
based on the quasi-Newton property:

$\displaystyle\eta_{b}^{t}=\frac{\left\|\bm{\gamma}_{b}^{t}-\bm{\gamma}_{b}^{t-1}\right\|_{2}^{2}}{(\bm{\gamma}_{b}^{t}-\bm{\gamma}_{b}^{t-1})^{H}(\widetilde{\bigtriangledown
f}(\bm{\gamma}_{b}^{t})-\widetilde{\bigtriangledown
f}(\bm{\gamma}_{b}^{t-1}))},$ (26)

which is the local estimation of Lipschitz constant and approximates the
inverse of Hessian matrix of $f(\cdot)$ at $\bm{\gamma}_{b}^{t}$. However, the
denominator of (26) is potentially small due to the fact that
$\bm{\gamma}_{b}^{t}$ is highly sparse. Also, the gradient even remains
constant between two consecutive iteration indices due to the application of
gradient in a probabilistic manner in (24). In consequence,
$\tilde{\eta}_{b}^{t}$ may approach infinity. To avoid this case, we introduce
the term $\epsilon\left\|\bm{\gamma}_{b}^{t}-\bm{\gamma}_{b}^{t-1}\right\|$ to
the denominator of (26) to keep the its denominator positive and control the
lower bound of the denominator in each iteration. The resulting step size is
explicitly given by

$\displaystyle\eta_{b}^{t}=\frac{\left\|\bm{\gamma}_{b}^{t}-\bm{\gamma}_{b}^{t-1}\right\|_{2}^{2}}{\left|(\bm{\gamma}_{b}^{t}\\!-\\!\bm{\gamma}_{b}^{t-1})^{H}(\widetilde{\bigtriangledown
f}(\bm{\gamma}_{b}^{t})\\!-\\!\widetilde{\bigtriangledown
f}(\bm{\gamma}_{b}^{t-1}))\right|\\!+\\!\epsilon\left\|\bm{\gamma}_{b}^{t}\\!-\\!\bm{\gamma}_{b}^{t-1}\right\|_{2}^{2}},$

where $\epsilon\geq 0$ is the adjustment parameter. Noting that (III-C) can
adaptive tune the step size without imposing much extra computational
complexity burden.

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

$\displaystyle\gamma_{bn}^{t+1}=\left\\{\begin{array}[]{l}\min\left(\tau\eta_{b}^{l,t}\frac{z_{bn}^{t}+\tau\eta_{b}^{l,t}x_{bn}^{l,t}-{\gamma}_{ln}^{t}}{|z_{bn}^{t}\\!+\\!\tau\eta_{b}^{l,t}x_{bn}^{l,t}-{\gamma}_{ln}^{t}|},z_{bn}^{t}+\tau\eta_{b}^{l,t}x_{bn}^{l,t}\right)\\\
+z_{bn}^{t}+\tau\eta_{b}^{l,t}x_{bn}^{l,t},~{}\text{if}~{}~{}z_{bn}^{t}+\tau\eta_{b}^{l,t}x_{bn}^{l,t}\neq
0,\\\
0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\text{if}~{}~{}z_{bn}^{t}+\tau\eta_{b}^{l,t}x_{bn}^{l,t}=0,\end{array}\right.$
(31)

where the minimum selection 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. After
getting $\gamma_{bn}^{t+1}$, the covariance matrix $\bm{\Sigma}_{b}$ in (19)
can be updated according to the distance between $\gamma_{bn}^{t+1}$ and
$\gamma_{bn}^{t}$, which is given by

$\bm{\Sigma}_{b}^{t+1}=\bm{\Sigma}_{b}^{t}+(\gamma_{bn}^{t+1}-\gamma_{bn}^{t})\mathbf{s}_{n}\mathbf{s}_{n}^{H}.$
(32)

From (17) and (31), the estimation performance depends, to a great extent, on
the cooperation strategy specified by the combiner $c_{lb}^{t}$. Similar to
work [34], this paper adopts the following adaptive combiner

$\displaystyle
c_{lb}^{t}=\left\\{\begin{array}[]{l}\frac{2}{\left|\mathcal{N}_{b}^{-}\right|_{c}}\frac{1}{1+\exp(\rho\|\bm{\gamma}_{b}^{t-1}-\bm{\gamma}_{l}^{t-1}\|_{2})},~{}l\in\mathcal{N}_{b}^{-},\\\
1-\sum_{l\in\mathcal{N}_{b}^{-}}c_{lb}^{t},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}l=b,\\\
0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}l\notin\mathcal{N}_{b},\end{array}\right.$
(36)

where $\rho>0$ is a large constant set beforehand. Note that the term
$\|\bm{\gamma}_{b}^{t-1}-\bm{\gamma}_{l}^{t-1}\|_{2}$ in (36) 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 value of combiner, or even discard the information from this
neighbor. Conversely, the $b$th AP would increase the combiner when the
distance of estimation between two APs is small.

Once the estimate $\\{\hat{\bm{\gamma}}_{b}\\}_{b=1}^{B}$ of the device state
vectors is obtained, we employ the element-wise thresholding at each AP to
determine $\chi_{n}$ from $\hat{\gamma}_{bn}$, i.e., $\chi_{n}=1$, if
$\hat{\gamma}_{bn}>\imath_{b}$, and $\chi_{n}=0$ otherwise, where the
threshold $\imath_{b}=\omega_{b}\sigma^{2}$ for a tunable parameter
$\omega_{b}>0$. For clarity, the pseudo-code of the cooperative activity
detection (CAD) algorithm for sourced random access is summarized in Algorithm
1. Different from the cooperative MIMO [22], where each BS performs signal
processing locally, the vectors of log-likelihood ratio are then forwarded to
the central unit for decision making, the proposed CAD algorithm exchanges
device state vector iteratively among neighboring APs. In this context, the
access latency in cooperative MIMO is longer than that of the proposed CAD
algorithm. On the other hand, the proposed CAD algorithm is more reliable and
robust to AP and/or fronthaul link failure compared with cooperative MIMO.

Algorithm 1 CAD for Sourced Random Access

1: Input: $\\{\mathbf{Y}_{b}\\}_{b=1}^{B}$, $\mathbf{S}$,
$\\{\hat{\bm{\Sigma}}_{b\mathbf{y}}=\frac{1}{M}\mathbf{Y}_{b}\mathbf{Y}_{b}^{H}\\}_{b=1}^{B}$,
step size $\\{\eta_{b}^{0}\\}_{b=1}^{B}$, and total iterations $T$

2: Initialization: $\\{\bm{\gamma}_{b}^{0}=\mathbf{0}\\}_{b=1}^{B}$,
$\\{\bm{\Sigma}_{b}^{0}=\sigma^{2}\mathbf{I}\\}_{b=1}^{B}$,
$\\{\mathbf{x}_{b}^{l,0},l\in\mathcal{N}_{b}\\}_{b=1}^{B}$,
$\\{c_{lb}^{0}\\}_{b=1}^{B}$,
$\\{\mathbf{x}_{b}^{0}=\sum_{l\in\mathcal{N}_{b}}c_{lb}^{0}\mathbf{x}_{b}^{l,0}\\}_{b=1}^{B}$

3: for $t=1:T$ do

4: for each AP $b$:

5: Adaptation:

6: Compute $\widetilde{\bigtriangledown f}(\bm{\gamma}_{b}^{t})$ based on (24)

7: Compute $\mathbf{z}_{b}^{t}$ based on (25)

8: Compute adaptive step size $\eta_{b}^{t}$ based on (III-C)

9: Choose $l$ randomly from the set $\mathcal{N}_{b}$ with probabilities
$\\{p_{1},p_{2},\cdots,p_{\left|\mathcal{N}_{b}\right|_{c}}\\}$

10: Compute adaptive combiner $c_{lb}^{t}$ based on (36)

11: Compute $\eta_{b}^{l,t}=\frac{c_{lb}^{t}\eta_{b}^{t}}{p_{l}}$

12: Randomly select a permutation $t_{1},t_{2},\cdots,t_{N}$ of the coordinate
indices $\\{1,2,\cdots,N\\}$ of $\bm{\gamma}_{b}^{t}$

13: for $n=1:N$ do

14: Update $\gamma_{bt_{n}}^{t+1}$ based on (31)

15:
$\bm{\Sigma}_{b}^{t+1}=\bm{\Sigma}_{b}^{t}+(\gamma_{bt_{n}}^{t+1}-\gamma_{bt_{n}}^{t})\mathbf{s}_{t_{n}}\mathbf{s}_{t_{n}}^{H}$

16: end for

17: Compute $\mathbf{x}_{b}^{l,t+1}$ based on (16)

18: Compute $\mathbf{x}_{b}^{t+1}$ based on (17)

19: Communication:

20: Transmit $\bm{\gamma}_{b}^{t}$ to its one-hop neighbor AP

21: end for

22: Output: $\\{\bm{\gamma}_{b}^{t+1}\\}_{b=1}^{B}$

##  IV Cooperative Unsourced Random Access

In this section, we design an unsourced random access scheme for 6G cell-free
wireless networks. Indeed, most works on unsourced random access comprise an
inner code for recovering the submessage and an outer code for stitching the
submessages [19, 21]. In this context, the proposed CAD algorithm in Section
III can be employed as an inner decoder of inner code for unsourced random
access. Hence, we can provide a unified algorithm for both sourced and
unsourced random access paradigms.

For the proposed cooperative unsourced random access scheme, the message of
each device is split into several subblocks, which are independently
transmitted and cooperatively recovered at all APs, and then to be integrated
to reconstruct the original message. Herein, the transmission and the recovery
of the subblock messages are called cooperative inner code, while the split
and integration of the message is called outer code. In what follows, we
provide the details about the cooperative inner and outer code by exploiting
the characteristics of 6G cell-free wireless networks.

###  IV-A Cooperative Inner Code

It is assumed that each active device has a message of length $\bar{b}=uR$
bits to send with $u=\mathcal{Z}L$, where $u$ denotes the length of signal
samples, $\mathcal{Z}$ denotes the number of subblocks, $R$ denotes the sample
rate, and $L$ denotes the subslot length. The active devices synchronously
transmit their messages over the coherence time block to all the APs in 6G
cell free wireless network. Then, each AP $b$ produces a list
$\mathcal{L}_{b}$ of the transmitted messages $\\{m_{b,k}:k\in\mathcal{K}\\}$
of the active devices. For some integer $J>0$, the $\bar{b}$-bit message is
divided into $\mathcal{Z}$ subblocks of size
$q_{1},q_{2},\cdots,q_{\mathcal{Z}}$ satisfying the following conditions:
$\sum_{i}q_{i}=\bar{b}$, $q_{1}=J$, and $q_{i}<J$ for all
$i=2,\cdots,\mathcal{Z}$. Herein, all subblocks $i=2,\cdots,\mathcal{Z}$ are
augmented to size $J$ by appending the parity bits and the pseudo-random
parity-check equations are identical for all devices. Such a systematic linear
block code based on random parity checks is known as the tree code proposed in
[19].

Then, the inner code is used to transmit in sequence the $\mathcal{Z}$
subblocks under the subslot length $L$. Let $\mathbf{S}\in\mathbb{C}^{L\times
2^{J}}=[\mathbf{s}_{1},\mathbf{s}_{2},\cdots,\mathbf{s}_{2^{J}}]$ represent a
common codebook for all devices with columns being the normalized codewords
such that $\|\mathbf{s}_{e}\|_{2}^{2}=L,\forall e\in\\{1,2,\cdots,2^{J}\\}$.
For subblock $i$, let $e_{k}(i)$ denote the $J$-bit messages produced by the
$k$th active device characterized by integers in $\\{1,2,\cdots,2^{J}\\}$. The
active device $k$ simply sends the column $\mathbf{s}_{e_{k}(i)}$ of the
coding matrix $\mathbf{S}$ to all APs. As such, the received signal at the
$b$th AP during the $i$th subblock can be expressed as

$\displaystyle\mathbf{Y}_{b,i}$
$\displaystyle=\sum_{k\in\mathcal{K}}\sqrt{g_{b,k}}\mathbf{s}_{e_{k}(i)}\mathbf{h}_{b,k}^{T}+\mathbf{W}_{b,i}$
$\displaystyle=\mathbf{S}\bm{\Xi}_{b,i}\mathbf{G}_{b}^{1/2}\mathbf{H}_{b}+\mathbf{W}_{b,i},~{}\forall
i\in\\{1,2,\cdots,\mathcal{Z}\\},$ (37)

where $g_{b,k}$, $\mathbf{h}_{b,k}\in\mathbb{C}^{M}$, and
$\mathbf{H}_{b}\in\mathbb{C}^{N\times M}$ are the large-scale fading
component, the small-scale fading vector, and the corresponding small-scale
fading matrix (Gaussian i.i.d. entries $\thicksim\mathcal{CN}(0,1)$) as
described in Section III, $\mathbf{G}_{b}=\text{diag}(g_{b,1},\cdots,g_{b,N})$
denotes the diagonal large-scale fading matrix, $\mathbf{W}_{b,i}$ is the AWGN
matrix with i.i.d. entries $\thicksim\mathcal{CN}(0,\sigma^{2})$, and
$\bm{\Xi}_{b,i}\in\\{0,1\\}^{2^{J}\times N}$ is a binary selection matrix
where for the active device $k$ the corresponding column is all-zero but has a
single one in position $e_{k}(i)$, while for all inactive devices the
corresponding columns contain all zeros.

The elements of the $r$th row of the matrix
$\bm{\Xi}_{b,i}\mathbf{G}_{b}^{1/2}\mathbf{H}_{b}$ follow the distribution of
$\mathcal{CN}(0,\sum_{k\in\mathcal{K}}g_{b,k}\xi_{r,k}^{b,i})$, where
$\xi_{r,k}^{b,i}$ denotes the $(r,k)$-th element of $\bm{\Xi}_{b,i}$, which
equals to one if $r=e_{k}(i)$ and zero otherwise. Let
$\gamma_{bi}^{r}=\sum_{k\in\mathcal{K}}g_{b,k}\xi_{r,k}^{b,i}$ and
${\bm{\Gamma}}_{b,i}=\text{diag}(\gamma_{bi}^{1},\cdots,\gamma_{bi}^{2^{J}})\in\mathbb{C}^{2^{J}\times
2^{J}}$, then, the received signal can be transformed as

$\displaystyle\mathbf{Y}_{b,i}$ $\displaystyle=$
$\displaystyle\mathbf{S}\bm{\Gamma}_{b,i}^{\frac{1}{2}}\bar{\mathbf{H}}_{b}+\mathbf{W}_{b,i},i=1,2,\cdots,\mathcal{Z},$
(38)

where $\bar{\mathbf{H}}_{b}$ has the i.i.d. elements
$\thicksim\mathcal{CN}(0,1)$. Note that for a given subblock $i$, (38) has the
same form as (1). Moreover, it is seen that the dimensions of the variable
matrices in (38) are independent of the number of total devices $N$.

The decoding of the inner code is equivalent to recovering the non-zero
elements of the diagonal matrix $\bm{\Gamma}_{b,i}$ from the noisy
measurements $\mathbf{Y}_{b,i}$ for all APs and subblocks. Under this context,
given subblock $i$, the codeword state vectors
$\\{\bm{\gamma}_{bi}=[\gamma_{bi}^{1},\cdots,\gamma_{bi}^{2^{J}}]\\}_{b=1}^{B}$
has the same characteristics as the device state vectors
$\\{\bm{\gamma}_{b}\\}_{b=1}^{B}$. Specifically, they are jointly sparse and
the amplitudes of the nonzero entries at the same position of neighboring APs
are similar to each other. It is interesting to find that the decoding of
inner code is the same as the cooperative activity detection for sourced
random access. Hence, the proposed CAD algorithm in the last section can be
employed as the cooperative decoding algorithm of the inner code.

###  IV-B Outer Code in Cooperative Systems

Once the estimation of the codeword state vector,
$\hat{\bm{\gamma}}_{b}^{i}=[\hat{\gamma}_{b1}^{i},\hat{\gamma}_{b2}^{i},\cdots,\hat{\gamma}_{b2^{J}}^{i}]^{T}$,
is obtained, it is then passed to the outer decoder for further identify all
possible subblock messages by separating data and parity bits. In this paper,
we utilize the concatenated code as the outer code. In what following, we
provide a decoding approach for the outer code.

For the $b$th AP, the set of active columns of codebook $\mathbf{S}$ at
subblock $i$ is given by

$\displaystyle\mathcal{C}_{b,i}=\\{r\in\\{1,2,\cdots,2^{J}\\}:\hat{\gamma}_{br}^{i}\geq\imath_{b,i}\\},$
(39)

where the threshold is defined as $\imath_{b,i}=\nu_{b,i}\sigma^{2}$ with a
tuneable parameter $\nu_{b,i}>0$ for AP $b$ and subblock $i$. The set
$\mathcal{C}_{b,i}$ forms the list of submessages at the $b$th AP and subblock
$i$. For the $b$th AP, since the subblocks contain parity bits, not all
message sequences in
$\mathcal{C}_{b,1}\times\mathcal{C}_{b,2}\times\cdots\times\mathcal{C}_{b,\mathcal{Z}}$
are needed. Therefore, the outer decoder is used to identify all possible
message sequences. In specific, starting from $i=1$ and proceeding in order,
the integer indices $\mathcal{C}_{b,i}$ is converted back to their binary
representation, and then the sequences satisfying the parity are output as the
original messages. A description of the CAD for unsourced random access is
summarized in Algorithm 2.

Algorithm 2 CAD for Unsourced Random Access

1: Input: $\\{\mathbf{Y}_{b}\\}_{b=1}^{B}$ and
$\mathbf{S}\in\mathbb{C}^{L\times 2^{J}}$

2: Outer encoder: The $\mathcal{Z}$ subblocks of messages are encoded via a
tree code

3: Inner encoder: Transmit in sequence the $\mathcal{Z}$ subblocks to all APs
based on $\mathbf{S}$

4: for $i=1:\mathcal{Z}$ do

5: for each AP $b$:

6: Cooperative Inner decoder: Inner decoding with the Algorithm 1 to obtain
the estimation $\hat{\bm{\gamma}}_{b}^{i}$

7: end for

8: Outer decoder: Find the set $\\{\mathcal{C}_{b,i}\\}_{b=1}^{B}$ by hard
threshold and stitch together the sequence of submessages

9: Output: The list $\\{\mathcal{L}_{b}\\}_{b=1}^{B}$ of the transmitted
messages

##  V Performance Analysis

As mentioned earlier, the CAD algorithm is the key of both cooperative sourced
and unsourced random access. In this section, we provide the analysis of
computational complexity, communication cost, and the convergence of the
proposed CAD algorithm.

###  V-A Computational Complexity and Communication Cost

In what follows, the computational complexity and communication cost of the
proposed CAD algorithm is analyzed. In each iteration of Algorithm 1, for an
arbitrary AP, the computational complexity mainly arises from the matrix
multiplication, and the overall computational complexity of the CAD algorithm
is $\mathcal{O}(L^{2})$. 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 start of iteration.
For an arbitrary AP, we compare the proposed algorithm with two detection
algorithms from the perspective of the computational complexity, including the
AMP algorithm [22] and the OAMP-MMV algorithm [11]. It can be seen in Table I
that the computational complexity of the proposed CAD algorithm is superior to
other two algorithms in mMTC applications.

TABLE I: Computation Complexity of the Considered Schemes. | Schemes
---
| Computational
---
complexity
Proposed CAD | $\mathcal{O}(TL^{2}+L^{2}M)$
AMP | $\mathcal{O}(T(MNL+LM^{2}+NM^{3}))$
OAMP-MMV | $\mathcal{O}(T(MN\log_{2}N+NM^{3}))$

For the communication cost, in each iteration, each AP needs to transmit
$N$-dimensional intermediate $\bm{\gamma}_{b}^{t}$. Thus, for all APs, the CAD
algorithm needs to exchange
$N\sum_{b=1}^{B}\left|\mathcal{N}_{b}^{-}\right|_{c}$ number of parameters.
For the activity detection in cooperative unsourced random access, the role of
the total number of potential devices $N$ in the device activity detection
problem is replaced by the number of messages $2^{J}\ll N$ with the small size
$J$. Thus the communication cost of the CAD algorithm reduce from
$N\sum_{b=1}^{B}\left|\mathcal{N}_{b}^{-}\right|_{c}$ to
$2^{J}\sum_{b=1}^{B}\left|\mathcal{N}_{b}^{-}\right|_{c}$.

_Remark 1_ : It is interesting to emphasize that the computational complexity
and the communication cost at each iteration of the CAD algorithm do not grow
as the number of antennas at each AP, $M$, increases. Moreover, cooperative
unsourced random access implies that both the computational complexity and the
communication cost of the proposed CAD algorithm do not grow by increasing the
total number of potential devices. Instead, the overhead is determined by the
length of subblock $J$. In practical applications, the wireless fronthaul
links among APs usually have limited capacity in 6G cell-free wireless
network. Thus, $J$ can be flexibly designed to strike a balance between the
detection accuracy and required signaling overhead. In addition, a small
number of APs can be set for cooperation, as will be verified by simulations
in Section VI, where the activity detection error can approach zero with the
existence of limited fronthaul links.

###  V-B Convergence Analysis

In this subsection, we establish the convergence (to a stationary point of
$\bm{\gamma}_{b}^{t}$) of the proposed CAD algorithm. Before giving a more
detailed convergence result, let us first provide the following lemma which is
instrumental in analyzing the convergence of the proposed algorithm.

_Lemma 1_ : Let
$\bm{\zeta}_{b}^{t}=\bm{\gamma}_{b}^{t}-\bar{p}\eta_{b}^{t}\triangledown
f(\bm{\gamma}_{b}^{t})$ and
$\bm{\zeta}_{b}^{*}=\bm{\gamma}_{b}^{*}-\bar{p}\eta_{b}^{t}\bigtriangledown
f(\bm{\gamma}_{b}^{*})$, where $\bm{\gamma}_{b}^{*}$ is the optimal solution
of minimizing the problem in (9). When $f(\cdot)$ is $\mathcal{L}_{f}$-smooth
(i.e., $\|\bigtriangledown f(\bm{\gamma}_{b}^{t})-\bigtriangledown
f(\bm{\gamma}_{b}^{t-1}\|\leq\mathcal{L}_{f}\|\bm{\gamma}_{b}^{t}-\bm{\gamma}_{b}^{t-1}\|$)
defined in [37], we have

$\displaystyle\mathbb{E}\left\\{\left\|\bm{\zeta}_{b}^{t}-\bm{\zeta}_{b}^{*}\right\|_{2}^{2}\right\\}\leq\mathbb{E}\left\\{\left\|\bm{\gamma}_{b}^{t}-\bm{\gamma}_{b}^{*}\right\|_{2}^{2}\right\\}-\bar{p}\varpi\eta_{b}^{t}d_{f}(\bm{\gamma}_{b}^{t},\bm{\gamma}_{b}^{*}),$
(40)

for an arbitrary AP $b$, the smoothness constant $\mathcal{L}_{f}$, and some
constant $\varpi>0$, if the step sizes are chosen to satisfy

$\displaystyle\frac{1}{\mathcal{L}+\epsilon}\leq\eta_{b}^{t}<\min\left(\frac{2}{\mathcal{L}_{f}},\frac{1}{\epsilon}\right),b=1,2,\cdots,B,$
(41)

where $d_{f}(\bm{\gamma}_{b}^{t},\bm{\gamma}_{b}^{*})$ denotes the Bregman
divergence of a function $f(\cdot)$ and is given by
$d_{f}(\bm{\gamma}_{b}^{t},\bm{\gamma}_{b}^{*}):=f(\bm{\gamma}_{b}^{t})-f(\bm{\gamma}_{b}^{*})-\left\langle\triangledown
f(\bm{\gamma}_{b}^{t}),\bm{\gamma}_{b}^{t}-\bm{\gamma}_{b}^{*}\right\rangle$.

###### Proof:

Please refer to Appendix A. ∎

As stated, the bound (40) in Lemma 1, which depends on the Bregman divergence,
will play an important role in verifying the convergence of the proposed CAD
algorithm. We then turn to prove the existence and the optimality of the fixed
points of step 7 and step 13 in Algorithm 1.

_Lemma 2_ : For an arbitrary AP $b$, the solution $\bm{\gamma}_{b}^{*}$ of
minimizing the problem in (9) exists for recursions step 7 and step 13 in
Algorithm 1, i.e., it holds that

$\displaystyle\bm{\gamma}_{b}^{*}=\text{prox}_{\beta\eta_{b}^{t}g}(\bm{\gamma}_{b}^{*}-\eta_{b}^{t}\widetilde{\bigtriangledown
f}(\bm{\gamma}_{b}^{*})-\tau\eta_{b}^{t}\mathbf{x}_{b}^{*}),$ (42)
$\displaystyle\bm{\gamma}_{b}^{*}=\text{prox}_{\tau\eta_{b}^{l,t}\Psi_{l}}(\mathbf{z}_{b}^{*}+\tau\eta_{b}^{l,t}\mathbf{x}_{b}^{l,*}),$
(43)

for any $\eta_{b}^{t}$ and $\eta_{b}^{l,t}$, if

$\displaystyle\bigtriangledown
f(\bm{\gamma}_{b}^{*})+\tau\sum_{l\in\mathcal{N}_{b}}c_{lb}^{t}\mathbf{x}_{b}^{l,*}+\beta\bm{\varrho}_{b}^{*}=\mathbf{0},$
(44)

holds for the vectors
$\mathbf{x}_{b}^{*}\in\partial\Psi(\bm{\gamma}_{b}^{*})$,
$\mathbf{x}_{b}^{l,*}\in\partial\Psi_{l}(\bm{\gamma}_{b}^{*})$, and
$\bm{\varrho}_{b}^{*}\in\partial g(\bm{\gamma}_{b}^{*})$. Condition (44) means
that when the optimality is obtained, the subgradient of local cost function
in (9) reaches $0$.

###### Proof:

Please refer to Appendix B. ∎

Note that there exists a particular fixed point $\bm{\gamma}_{b}^{*}$ and we
will show that the iterations of $\bm{\gamma}_{b}^{t}$ converge to this
particular fixed point in the following. We collect the information from
across all APs into block vectors and matrices. In particular, for the $b$th
AP, we define the mean of the variables $\bm{\gamma}_{b}^{t}$ with respect to
time series $t$ as
$\bar{\bm{\gamma}}_{b}^{t}=\frac{1}{t}\sum_{j=0}^{t-1}\bm{\gamma}_{b}^{t}$.
The global stochastic quantities are defined as

$\displaystyle\bm{\gamma}^{*}=\text{col}\\{\bm{\gamma}_{1}^{*},\bm{\gamma}_{2}^{*},\cdots,\bm{\gamma}_{B}^{*}\\},$
$\displaystyle\bar{\bm{\gamma}}^{t}=\text{col}\\{\bar{\bm{\gamma}}_{1}^{t},\bar{\bm{\gamma}}_{2}^{t},\cdots,\bar{\bm{\gamma}}_{B}^{t}\\},$
$\displaystyle\mathbf{h}^{t}=\text{col}\\{d_{f}(\bar{\bm{\gamma}}_{1}^{t},\bm{\gamma}_{1}^{*}),d_{f}(\bar{\bm{\gamma}}_{2}^{t},\bm{\gamma}_{2}^{*}),\cdots,d_{f}(\bar{\bm{\gamma}}_{B}^{t},\bm{\gamma}_{B}^{*})\\}.~{}~{}~{}~{}~{}$
(45)

The next result provides a convergence rate of the proposed CAD algorithm.

_Theorem 1:_ If the step size satisfies
$\frac{1}{\mathcal{L}_{f}+\epsilon}\leq\eta_{b}^{t}<\min(\frac{2}{\mathcal{L}_{f}},\frac{1}{\epsilon})$,
$b=1,2,\cdots,B$, the inequality

$\displaystyle\mathbb{E}\left\\{\mathbf{h}^{t}\right\\}\leq\frac{1}{\bar{p}\varpi
t}\frac{1}{\bm{\eta}^{t}}\odot\bm{\mathcal{W}}^{0},$ (46)

holds where
$\frac{1}{\bm{\eta}^{t}}=\text{col}\\{\frac{1}{\eta_{1}^{t}},\frac{1}{\eta_{2}^{t}},\cdots,\frac{1}{\eta_{B}^{t}}\\}$,
the expectation is taken over random choice of neighbor AP $l$ of the $b$th
AP, and
$\bm{\mathcal{W}}^{0}=\text{col}\\{\mathcal{W}_{1}^{0},\mathcal{W}_{2}^{0},\cdots,\mathcal{W}_{N}^{0}\\}$
with

$\displaystyle\mathcal{W}_{b}^{0}=\|\bm{\gamma}_{b}^{0}-\bm{\gamma}_{b}^{*}\|_{2}^{2}+\tau^{2}\sum_{l\in\mathcal{N}_{b}}(\eta_{b}^{l,0})^{2}\|\mathbf{x}_{b}^{l,0}-\mathbf{x}_{b}^{l,*}\|_{2}^{2}.$
(47)

###### Proof:

Similar to the Lyapunov function adopted in [38] and [39], which can be used
to prove the stability of an equilibrium, we define the Lyapunov function for
each AP of the CAD algorithm at iteration $t$ as

$\displaystyle\mathcal{W}_{b}^{t+1}$ $\displaystyle=$
$\displaystyle\mathbb{E}\left\\{\|\bm{\gamma}_{b}^{t+1}-\bm{\gamma}_{b}^{*}\|_{2}^{2}\right\\}$
(48) $\displaystyle+$
$\displaystyle\sum_{l\in\mathcal{N}_{b}}(\tau\eta_{b}^{l,t+1})^{2}\mathbb{E}\left\\{\|\mathbf{x}_{b}^{l,t+1}-\mathbf{x}_{b}^{l,*}\|_{2}^{2}\right\\},$

which measures both distance between $\bm{\gamma}_{b}^{t}$ and
$\bm{\gamma}_{b}^{*}$ and the sum of the distance between
$\mathbf{x}_{b}^{l,t}$ and $\mathbf{x}_{b}^{l,*}$ for all the neighbors of the
$b$th AP.

For establishing the upper bound on the first term of $\mathcal{W}_{b}^{t}$,
we start by applying the non-expansiveness property of the proximal map [32],
which implies

$\displaystyle\\!\\!\\!\\!\\!\\!(1+\frac{1}{\tau\eta_{b}^{l,t}\mathcal{L}_{\Psi_{l}}})\|\mathbf{x}-\text{prox}_{\tau\eta_{b}^{l,t}\Psi_{l}}(\mathbf{x})-(\mathbf{y}-\text{prox}_{\tau\eta_{b}^{l,t}\Psi_{l}}(\mathbf{y}))\|_{2}^{2}$
$\displaystyle\\!\\!\\!\\!\\!\\!+\|\text{prox}_{\tau\eta_{b}^{l,t}\Psi_{l}}(\mathbf{x})-\text{prox}_{\tau\eta_{b}^{l,t}\Psi_{l}}(\mathbf{y}))\|_{2}^{2}\leq\|\mathbf{x}-\mathbf{y}\|_{2}^{2},$
(49)

when the function $\Psi_{l}$ is non-smooth with the smoothness constant
$\mathcal{L}_{\Psi_{l}}=+\infty$ for any $\mathbf{x}$ and $\mathbf{y}$ [40].
Considering step 13 and step 16 in Algorithm 1, we substitute
$\mathbf{x}=\mathbf{z}_{b}^{t}+\tau\eta_{b}^{l,t}\mathbf{x}_{b}^{l,t}$ and
$\mathbf{y}=\bm{\gamma}_{b}^{*}+\tau\eta_{b}^{l,t}\mathbf{x}_{b}^{l,*}$ into
the inequality in (V-B) with the optimality condition of Lemma 2:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\|\bm{\gamma}_{b}^{t+1}-\bm{\gamma}_{b}^{*}\|_{2}^{2}+(1+\frac{1}{\tau\eta_{b}^{l,t}\mathcal{L}_{\Psi_{l}}})$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\cdot\|\mathbf{z}_{b}^{t}+\tau\eta_{b}^{l,t}\mathbf{x}_{b}^{l,t}-\bm{\gamma}_{b}^{t+1}-(\bm{\gamma}_{b}^{*}+\tau\eta_{b}^{l,t}\mathbf{x}_{b}^{l,*}-\bm{\gamma}_{b}^{*})\|_{2}^{2}=$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\|\bm{\gamma}_{b}^{t+1}\\!-\\!\bm{\gamma}_{b}^{*}\|_{2}^{2}\\!+\\!(1\\!+\\!\frac{1}{\tau\eta_{b}^{l,t}\mathcal{L}_{\Psi_{l}}})(\tau\eta_{b}^{l,t})^{2}\|\mathbf{x}_{b}^{l,t+1}\\!\\!-\\!\\!\mathbf{x}_{b}^{l,*}\|_{2}^{2}$
(50)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\leq\|\mathbf{z}_{b}^{t}+\tau\eta_{b}^{l,t}\mathbf{x}_{b}^{l,t}-(\bm{\gamma}_{b}^{*}+\tau\eta_{b}^{l,t}\mathbf{x}_{b}^{l,*})\|_{2}^{2}.$
(51)

Taking the expectation of the term in (51) with respect to the random selected
neighbor AP $l$ and putting the relation
$\eta_{b}^{l,t}=\frac{c_{lb}^{t}\eta_{b}^{t}}{p_{l}}$ into it, we obtain

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\mathbb{E}\left\\{\|\mathbf{z}_{b}^{t}\\!+\\!\tau\eta_{b}^{l,t}\mathbf{x}_{b}^{l,t}\\!-\\!(\bm{\gamma}_{b}^{*}+\tau\eta_{b}^{l,t}\mathbf{x}_{b}^{l,*})\|_{2}^{2}\right\\}\\!=\\!\|\mathbf{z}_{b}^{t}-\bm{\gamma}_{b}^{*}\|_{2}^{2}+$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!(\tau\eta_{b}^{t})^{2}\sum_{l\in\mathcal{N}_{b}}\frac{(c_{lb}^{t})^{2}}{p_{l}}\|\mathbf{x}_{b}^{l,t}\\!-\\!\mathbf{x}_{b}^{l,*}\|_{2}^{2}\\!+\\!2\tau\eta_{b}^{t}\left\langle\mathbf{z}_{b}^{t}\\!-\\!\bm{\gamma}_{b}^{*},\mathbf{x}_{b}^{t}\\!-\\!\mathbf{x}_{b}^{*}\right\rangle.~{}~{}~{}~{}~{}$
(52)

According to the Lemma 1 and the step 7 of Algorithm 1, we can put
$\mathbf{x}=\mathbf{z}_{b}^{t}=\text{prox}_{\beta\eta_{b}^{t}g}(\bm{\zeta}_{b}^{t}-\tau\eta_{b}^{t}\mathbf{x}_{b}^{t})$
and
$\mathbf{y}=\bm{\gamma}_{b}^{*}=\text{prox}_{\beta\eta_{b}^{t}g}(\bm{\zeta}_{b}^{*}-\tau\eta_{b}^{t}\mathbf{x}_{b}^{*})$
into (V-B). Then the first term of LHS of (V-B) can be bounded as follows

$\displaystyle\|\mathbf{z}_{b}^{t}-\bm{\gamma}_{b}^{*}\|_{2}^{2}\leq\|\bm{\zeta}_{b}^{t}-\tau\eta_{b}^{t}\mathbf{x}_{b}^{t}-(\bm{\zeta}_{b}^{*}-\tau\eta_{b}^{t}\mathbf{x}_{b}^{*})\|_{2}^{2}$
$\displaystyle-\left(1+\frac{1}{\beta\eta_{b}^{t}\mathcal{L}_{g}}\right)\|\bm{\zeta}_{b}^{t}-\tau\eta_{b}^{t}\mathbf{x}_{b}^{t}-\mathbf{z}_{b}^{t}-(\bm{\zeta}_{b}^{*}-\tau\eta_{b}^{t}\mathbf{x}_{b}^{*}-\bm{\gamma}_{b}^{*})\|_{2}^{2}$
$\displaystyle\leq-2\tau\eta_{b}^{t}\left\langle\bm{\zeta}_{b}^{t}-\bm{\zeta}_{b}^{*},\mathbf{x}_{b}^{t}-\mathbf{x}_{b}^{*}\right\rangle-\|\bm{\zeta}_{b}^{t}-\mathbf{z}_{b}^{t}-(\bm{\zeta}_{b}^{*}-\bm{\gamma}_{b}^{*})\|_{2}^{2}$
$\displaystyle+2\tau\eta_{b}^{t}\left\langle\bm{\zeta}_{b}^{t}-\mathbf{z}_{b}^{t}-(\bm{\zeta}_{b}^{*}-\bm{\gamma}_{b}^{*}),\mathbf{x}_{b}^{t}-\mathbf{x}_{b}^{*}\right\rangle+\|\bm{\zeta}_{b}^{t}-\bm{\zeta}_{b}^{*}\|_{2}^{2}$
$\displaystyle=\|\bm{\zeta}_{b}^{t}-\bm{\zeta}_{b}^{*}\|_{2}^{2}-\|\bm{\zeta}_{b}^{t}-\mathbf{z}_{b}^{t}-(\bm{\zeta}_{b}^{*}-\bm{\gamma}_{b}^{*})\|_{2}^{2}$
$\displaystyle-2\tau\eta_{b}^{t}\left\langle\mathbf{z}_{b}^{t}-\bm{\gamma}_{b}^{*},\mathbf{x}_{b}^{t}-\mathbf{x}_{b}^{*}\right\rangle,$
(53)

where $\mathcal{L}_{g}$ denotes the smoothness constant of sparsity-promoting
function, the second inequality is obtained by omitting the non-positive term
multiplied by the factor $1/(\beta\eta_{b}^{t}\mathcal{L}_{g})$, and the RHS
of the last equation is obtained by amalgamating the second and the third
terms of its left hand. Combining the inequality in (V-B), we obtain the upper
bound for (V-B) in a single expression as follows:

$\displaystyle\\!\\!\\!\\!\\!\\!\mathbb{E}\left\\{\|\mathbf{z}_{b}^{t}+\tau\eta_{b}^{l,t}\mathbf{x}_{b}^{l,t}-(\bm{\gamma}_{b}^{*}+\tau\eta_{b}^{l,t}\mathbf{x}_{b}^{l,*})\|_{2}^{2}\right\\}\leq\|\bm{\zeta}_{b}^{t}-\bm{\zeta}_{b}^{*}\|_{2}^{2}$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!-\\!\|\bm{\zeta}_{b}^{t}\\!-\\!\mathbf{z}_{b}^{t}\\!-\\!(\bm{\zeta}_{b}^{*}\\!-\\!\bm{\gamma}_{b}^{*})\|_{2}^{2}\\!+\\!(\tau\eta_{b}^{t})^{2}\sum_{l\in\mathcal{N}_{b}}\frac{(c_{lb}^{t})^{2}\|\mathbf{x}_{b}^{l,t}\\!-\\!\mathbf{x}_{b}^{l,*}\|_{2}^{2}}{p_{l}}.~{}~{}~{}~{}~{}$
(54)

To obtain the bound of the second term of $\mathcal{W}_{b}^{t}$, we add the
following missing sum at the both sides of the inequality based on (50) and
(51).

$\displaystyle\mathbb{E}\left\\{\sum_{l\in\mathcal{N}_{b}}(\tau\eta_{b}^{l,t+1})^{2}\|\mathbf{x}_{b}^{l,t+1}-\mathbf{x}_{b}^{l,*}\|_{2}^{2}\right\\}$
(55)
$\displaystyle=\mathbb{E}\left\\{(\tau\eta_{b}^{l,t+1})^{2}\|\mathbf{x}_{b}^{l,t+1}-\mathbf{x}_{b}^{l,*}\|_{2}^{2}\right\\}$
$\displaystyle+\mathbb{E}\left\\{\sum_{k\in\mathcal{N}_{b},k\neq
l}(\tau\eta_{b}^{k,t+1})^{2}\|\mathbf{x}_{b}^{k,t+1}-\mathbf{x}_{b}^{k,*}\|_{2}^{2}\right\\}$
$\displaystyle=\mathbb{E}\left\\{(\tau\eta_{b}^{l,t+1})^{2}\|\mathbf{x}_{b}^{l,t+1}-\mathbf{x}_{b}^{l,*}\|_{2}^{2}\right\\}$
$\displaystyle+\mathbb{E}\left\\{\sum_{k\in\mathcal{N}_{b},k\neq
l}(\tau\eta_{b}^{l,t})^{2}\|\mathbf{x}_{b}^{k,t}-\mathbf{x}_{b}^{k,*}\|_{2}^{2}\right\\}$
$\displaystyle=\mathbb{E}\left\\{(\tau\eta_{b}^{l,t+1})^{2}\|\mathbf{x}_{b}^{l,t+1}-\mathbf{x}_{b}^{l,*}\|_{2}^{2}\right\\}$
$\displaystyle+\sum_{l\in\mathcal{N}_{b}}(1-p_{l})(\tau\eta_{b}^{l,t})^{2}\|\mathbf{x}_{b}^{l,t}-\mathbf{x}_{b}^{l,*}\|_{2}^{2}$
$\displaystyle=\sum_{l\in\mathcal{N}_{b}}(\tau\eta_{b}^{l,t})^{2}\|\mathbf{x}_{b}^{l,t}-\mathbf{x}_{b}^{l,*}\|_{2}^{2}-(\tau\eta_{b}^{t})^{2}\sum_{l\in\mathcal{N}_{b}}\frac{(c_{lb}^{t})^{2}}{p_{l}}\|\mathbf{x}_{b}^{l,t}-\mathbf{x}_{b}^{l,*}\|_{2}^{2}$
$\displaystyle+\mathbb{E}\left\\{(\tau\eta_{b}^{l,t+1})^{2}\|\mathbf{x}_{b}^{l,t+1}-\mathbf{x}_{b}^{l,*}\|_{2}^{2}\right\\},$
(56)

where the second equality stems from the fact that the second term of the RHS
of the first equality does not change at iteration $t$, and the last equality
is obtained by incorporating the definition of $\eta_{b}^{l,t}$. Adding (55)
and (56) to (50) and (51), respectively. Then, combining (V-B), it is found
that the last term of RHS of (V-B) vanishes because of the same in (56). We
arrive at

$\displaystyle\sum_{l\in\mathcal{N}_{b}}(\tau\eta_{b}^{l,t+1})^{2}\mathbb{E}\left\\{\|\mathbf{x}_{b}^{l,t+1}-\mathbf{x}_{b}^{l,*}\|_{2}^{2}\right\\}+\mathbb{E}\left\\{\|\bm{\gamma}_{b}^{t+1}-\bm{\gamma}_{b}^{*}\|_{2}^{2}\right\\}$
$\displaystyle\leq\sum_{l\in\mathcal{N}_{b}}(\tau\eta_{b}^{l,t})^{2}\mathbb{E}\left\\{\|\mathbf{x}_{b}^{l,t}-\mathbf{x}_{b}^{l,*}\|_{2}^{2}\right\\}-\frac{(\tau\eta_{b}^{l,t})^{2}\|\mathbf{x}_{b}^{l,t+1}-\mathbf{x}_{b}^{l,*}\|_{2}^{2}}{\tau\eta_{b}^{l,t}\mathcal{L}_{\Psi_{l}}}$
$\displaystyle+\mathbb{E}\left\\{\|\bm{\zeta}_{b}^{t}-\bm{\zeta}_{b}^{*}\|_{2}^{2}\right\\}-\mathbb{E}\left\\{\|\bm{\zeta}_{b}^{t}-\mathbf{z}_{b}^{t}-(\bm{\zeta}_{b}^{*}-\bm{\gamma}_{b}^{*})\|_{2}^{2}\right\\}$
$\displaystyle\leq\underbrace{\mathbb{E}\left\\{\left\|\bm{\gamma}_{b}^{t}-\bm{\gamma}_{b}^{*}\right\|_{2}^{2}\right\\}+\sum_{l\in\mathcal{N}_{b}}(\tau\eta_{b}^{l,t})^{2}\mathbb{E}\left\\{\|\mathbf{x}_{b}^{l,t}-\mathbf{x}_{b}^{l,*}\|_{2}^{2}\right\\}}_{{\mathcal{W}}_{b}^{t}}$
$\displaystyle-\bar{p}\varpi\eta_{b}^{t}d_{f}(\bm{\gamma}_{b}^{t},\bm{\gamma}_{b}^{*})-\frac{(\tau\eta_{b}^{l,t})^{2}\|\mathbf{x}_{b}^{l,t+1}-\mathbf{x}_{b}^{l,*})\|_{2}^{2}}{\tau\eta_{b}^{l,t}\mathcal{L}_{\Psi_{l}}}$
$\displaystyle-\mathbb{E}\left\\{\|\bm{\zeta}_{b}^{t}-\mathbf{z}_{b}^{t}-(\bm{\zeta}_{b}^{*}-\bm{\gamma}_{b}^{*})\|_{2}^{2}\right\\},$
(57)

where the last inequality follows from the fact that, as seen in the Lemma 1,
the difference between $\bm{\zeta}_{b}^{t}$ and $\bm{\zeta}_{b}^{*}$ has an
upper bound if the step size satisfies the condition in (41). Omitting the
last two terms of the expression (V-B) leads to

$\displaystyle{\mathcal{W}}_{b}^{t+1}\leq{\mathcal{W}}_{b}^{t}-\bar{p}\varpi\eta_{b}^{t}d_{f}(\bm{\gamma}_{b}^{t},\bm{\gamma}_{b}^{*}).$
(58)

By telescoping the inequality in (58) from 0 to $t$, we have

$\displaystyle\\!\\!\\!\\!\sum_{j=0}^{t}\\!d_{f}(\bm{\gamma}_{b}^{t},\bm{\gamma}_{b}^{*})\\!\\!\\!$
$\displaystyle\leq$
$\displaystyle\\!\\!\\!t\mathbb{E}\left\\{d_{f}(\bar{\bm{\gamma}}_{b}^{t},\bm{\gamma}_{b}^{*})\right\\}\\!\leq\\!\frac{1}{\bar{p}\varpi\eta_{b}^{0}}({\mathcal{W}}_{b}^{0}\\!-\\!{\mathcal{W}}_{b}^{t+1}).~{}~{}~{}~{}~{}$
(59)

Omitting the last term inside the round brackets and then dividing out both
sides of (59) by $t$, we obtain

$\displaystyle\mathbb{E}\left\\{d_{f}(\bar{\bm{\gamma}}_{b}^{t},\bm{\gamma}_{b}^{*})\right\\}\leq\frac{1}{t\bar{p}\varpi\eta_{b}^{0}}{\mathcal{W}}_{b}^{0}.$
(60)

In the sequel, collecting above Bregman divergence bound over all APs results
in (46). ∎

From (60), we note that the proposed algorithm converges at the scale of
$\mathcal{O}(1/t)$ and it clearly shows how the probability $\bar{p}$ affects
the convergence rate. When $\bar{p}=1$, i.e., the case for exact the full
gradient calculation, the Bregman divergence bound reaches the minimum value.
In addition, the sparsity parameter $\beta$ has vanished in inequality (V-B),
while the similarity parameter $\tau$ contained in ${\mathcal{W}}_{b}^{0}$ of
(60) does not equal to $0$, which means that a small value of similarity
parameter $\tau$ is required for keeping admissible Bregman divergence bound.
The optimum parameters $\beta$ and $\tau$ can be designed based on the steady
state mean-square deviation (MSD), i.e.,
$\mathbb{E}(\|\text{col}\\{\bm{\gamma}_{1}^{*}-{\bm{\gamma}}_{1}^{t},\bm{\gamma}_{2}^{*}-{\bm{\gamma}}_{2}^{t},\cdots,\bm{\gamma}_{B}^{*}-{\bm{\gamma}}_{B}^{t}\\}\|_{2}^{2})$
with $t\rightarrow\infty$, which will be discussed in the future work.

##  VI Numerical Results

In this section, we present extensive simulation results to validate the
effectiveness of the proposed cooperative sourced and unsourced random access
schemes in 6G wireless networks. We simulate the 6G cell-free wireless network
comprising $B=20$ APs geographically distributed in a vast area to serve $N$
potential devices. The AP-to-AP distance is $0.5$ km and the radius of the
network coverage is set as $1.8$ km. The positive constants $\theta$ is
selected to be $1/0.039$. The penalty parameters $\beta$ and $\tau$ are set as
$0.38$ and $0.03$, respectively. $p_{l}$ is set to
$\frac{1}{|\mathcal{N}_{b}|_{c}}$, and the constant $\rho$ is set to $500$.

###  VI-A Activity Detection for Cooperative Sourced Random Access

We first conduct simulations to validate the effectiveness of the proposed CAD
algorithm for device activity detection. As a performance measure, we adopt
the activity error rate (AER). The AER is a sum of the missed detection
probability, defined as the probability that a device is active but is
declared to be inactive, and the false-alarm probability, defined as the
probability that a device is inactive but the detector declares it to be
active. As a reference, we compare the proposed CAD algorithm with two
baseline schemes: the conventional ML-based multi-cell massive MIMO [13] and
the AMP-based multi-cell massive MIMO [22], where each AP only serves its
cell’s devices without multi-cell cooperation and treats the inter-cell
interference as noise. The average signal-to-noise ratio (SNR) of generic AP
$b$ and active device $k\in\mathcal{K}$ over $L$ sequence dimensions is given
by
SNR${}_{b,k}=\frac{\|\mathbf{s}_{k}\|_{2}^{2}g_{b,k}\mathbb{E}\left\\{\|\mathbf{h}_{b,k}\|_{2}^{2}\right\\}}{\mathbb{E}\left\\{\|\mathbf{W}_{b}\|_{F}^{2}\right\\}}=\frac{g_{b,k}}{\sigma^{2}}$.

Figure 2: The AER for different step sizes with $N=1,000$, $K=300$, $M=32$,
and SNR = $10$ dB.

We first validate the analytic convergence rate analysis by numerically
comparing the detection performance of CAD, ML-based, and AMP-based multi-cell
methods in Fig. 2. Herein, in order to study the influence of the step size
separately, we fix $\bar{p}$ in the proposed CAD algorithm to 1. The
adjustment parameter $\epsilon$ is set to $90$, the number of APs for
cooperation is set to $5$, and the step size $\eta_{b}=0.03$ for the CAD
algorithm with a fixed step size which is identical for all the APs. The
simulation result shows that the proposed CAD with the adaptive step size
converges much faster than that of the proposed CAD algorithm with the fixed
step size. This is due to the fact that the CAD algorithm with adaptive step
size can always choose a more appropriate step size, which allows the
algorithm to converge fast at the very first few iterations and quickly
attains admissible results in terms of the relatively small generalization
error.

Figure 3: The AER for different probabilities $\bar{p}$ with $N=1,000$,
$K=200$, $M=32$, $L=120$, and SNR = $10$ dB.

Next, we study the impact of $\bar{p}$ on the performance of the CAD
algorithm. The parameter $\epsilon$ is set to $90$ and the number of APs for
cooperation is set to $5$ in Fig. 3. The results indicate that for relatively
large $\bar{p}$, it slightly influences the convergence rate and the detection
accuracy of the proposed scheme. Although the effect is obvious for smaller
$\bar{p}$, the case still obtains a huge performance gain compared with the
ML-based and the AMP-based multi-cell algorithms in terms of AER. It is
observed that the ML-based and the AMP-based multi-cell algorithms all
converge more rapidly at the beginning, but reach a plateau afterward and
their performance can not be further improved with the number of iterations.
Agreeing with the discussions in Section V. B, CAD algorithm enjoys small
computational complexity for $0<\bar{p}<1$, since fewer gradients need to be
calculated in each iteration. Our results also indicate that it should be
possible to throw away some gradient calculations while keeping the almost the
same convergence rate and detection accuracies.

Figure 4: The AER for different numbers of APs for cooperation with $N=1,000$,
$K=200$, and SNR = $10$ dB.

In the rest of the simulations, the probability of gradient calculation
$\bar{p}$ is set to 1 for unveiling the full potential of CAD under various
system settings. Fig. 4 depicts how detection performance of the proposed CAD
algorithm changes with different choices of the number of APs for cooperation.
Initially, in the regime with a few numbers of cooperation APs, the AER
decreases sharply as the number of cooperation APs increases. However, when
the number of cooperation APs continues to increase, the performance
improvement diminishes. In addition, it is observed that such a performance
saturation point value depends on the number of AP antennas $M$ and the pilot
sequence length $L$. Specifically, increasing $M$ or $L$ helps lower the
saturation point value, which indicates that the AP becomes more capable of
detecting the activity with a low communication cost. The reason for this
phenomenon is that for an arbitrary AP, more cooperative connections result in
more intermediate estimate exchanges in the proposed CAD algorithm, leading to
the improved detection performance. However, the channel strengths from a
specific active device to far away APs are negligible, and the exchange of
intermediate estimates with the remote AP can not further improve the massive
detection performance. The results also indicate that only a small number of
APs are required for effective cooperation which strikes a tradeoff between
the detection performance and the communication cost.

Figure 5: The AER for different lengths of sequence $L$ with $N=500$, $K=100$,
and SNR = $10$ dB.

In the rest of the simulations, the number of APs for cooperation is set to
$5$ for illustration. Fig. 5 shows the detection performance versus the length
of the sequence $L$. From this figure, we can intuitively observe that the AER
of all the considered algorithms decrease as the sequence length or the number
of antennas of each AP increases. We can also see that the proposed CAD
algorithm reduces the required length of sequences for achieving an accurate
device activity detection compared with both the ML-based and the AMP-based
multi-cell algorithms. Such an advantage of cooperative strategies mainly
comes from the fact that the proposed algorithm exploits the structured
sparsity and similarity of the multiple APs. Besides, the closed-form
expressions for the proximal operator of similarity-promoting and sparsity-
promoting terms are derived to achieve higher efficiency. In contrast, the ML-
based and the AMP-based multi-cell approaches ignore such prior information
and only perform activity detection for the devices distributed in its own
cell.

Figure 6: Comparison of the proposed CAD algorithm and centralized algorithms
with $N=1,000$, $K=100$, $M=64$ and SNR = $10$ dB.

Then, we compare the proposed CAD algorithm with the corresponding centralized
CAD algorithm and the centralized AMP algorithm. Note that when all APs are
connected to a centralized centre and the signals received at all APs are
centrally processed, the information from all APs is shared. In this context,
the CAD algorithm reduce to the centralized CAD algorithm since
$|\mathcal{N}_{b}|_{c}=N$ is considered. Moreover, applying the AMP algorithm
to jointly process the signals received at all APs can lead to the centralized
AMP algorithm. In Fig. 6, we show the AER curves of these algorithms. It is
observed that the CAD algorithm has some inevitable performance loss compared
with the centralized CAD algorithm due to the less shared information. Note
that although the centralized AMP algorithm can process the signals received
at all APs, it needs to estimate a larger number of unknown parameters
compared with the proposed CAD algorithm, which introduces an unnecessary
challenge to the estimation. On the other hand, the centralized AMP algorithm
does not consider the similarities among the device state vectors. Thus, the
proposed CAD is more efficient for activity detection than that of the
centralized AMP algorithm.

Figure 7: The AER for different activity probabilities with $N=500$, $L=100$,
and SNR = $10$ dB.

Fig. 7 plots the AER curves of the considered algorithms against the different
activity probability. For the comparison of the performance of the proposed
scheme and the multi-cell algorithms, we see that the performance of all the
algorithms is degraded as the activity probability increases. This is because
the co-channel interference among devices increases as more devices are
active. Note that the AMP algorithm requires the knowledge of large-scale
fading coefficients and the number of active devices $K$, while the CAD
algorithm only requires the sample covariance of channel observations. Thus,
the proposed CAD algorithm is robust to the inaccurate knowledge of $K$ and
the variation in the channel statistics compared with the AMP algorithm.
Moreover, the CAD algorithm does not require the explicit knowledge of channel
strengths which only needs to estimate a smaller number of unknown parameters,
thus, it is more efficient for activity detection than that of the AMP-based
multi-cell algorithm. The proposed algorithm outperforms the two multi-cell
algorithms by a large margin even if the activity probability is higher than
$0.6$. In practice, the proposed algorithm is appealing for the spread
applications of IoT with a wide range of activity probability.

Figure 8: The AER for different SNRs with $N=1,000$, $K=100$, $L=50$, and
$M=32$ antennas at each AP.

Fig. 8 plots the detection performance with different SNRs. It is observed
that for the considered range of SNRs, the ML-based and the AMP-based multi-
cell algorithms perform worse than that of the CAD algorithm, and the
performance gap is enlarged as the SNR increases. The reason is that the
inter-cell interference is a severe limiting factor for reliable activity
detection in multi-cell based algorithms, especially when inter-cell
interference is dominated in high SNR. While there are no cell boundaries in
6G cell-free wireless networks and the inter-cell interference can be avoided
via the APs’ cooperation.

###  VI-B Data Detection for Cooperative Unsourced Random Access

For the proposed cooperative massive access scheme designed for cell-free
massive MIMO based IoT, we further compare the data detection for the
unsourced random access. Note that the probability of error adopted here is
different from the definition in sourced random access. The system performance
is expressed in terms of probability of missed detection of the per device and
per AP, defined as the average fraction of transmitted messages not contained
in the list, i.e.,
$P_{b}^{\mathrm{MD}}=\frac{1}{K}\sum_{k\in\mathcal{K}}\mathbb{P}(m_{b,k}\notin\mathcal{L}_{b})$,
and the probability of false-alarm, defined as the average fraction of decoded
messages that were indeed not sent, i.e.,
$P_{b}^{\mathrm{FA}}=\frac{|\mathcal{L}_{b}\backslash\\{m_{b,k}:k\in\mathcal{K}\\}|_{c}}{|\mathcal{L}_{b}|_{c}}$.
We use the probability of error
$P_{e}=\sum_{b=1}^{B}(P_{b}^{\mathrm{MD}}+P_{b}^{\mathrm{FA}})/B$ to measure
the system detection performance. As a reference, we compare the proposed
approach with the ML-based multi-cell non-cooperative activity detection for
unsourced random access under the same spectral efficiency [21].

Figure 9: Error probability $P_{e}$ for different $M$ with $J=12$,
$\bar{b}=96$, $\mathcal{Z}=32$, $L=100$, $u=3,200$, and transmit power = $0$
dB.

In Fig. 9, we set the energy per bit, $E_{b}$, and per component AWGN variance
$\sigma^{2}$ such that transmit power $E_{b}/\sigma^{2}=0$ dB. Fig. 9 confirms
that as the number of antennas at each AP, $M$, increases, the detection error
of the proposed CAD algorithm drops rapidly and faster than that of the ML-
based multi-cell algorithm, indicating that employing cooperative AP can
quickly drive the detection error to zero with a fewer number of antennas for
saving cost. We note that when the number of antennas at each AP exceeds a
certain point, the error probability in data detection vanishes with the
proposed CAD algorithm. In other words, the superiority of the proposed CAD
algorithm is evident in massive MIMO systems, which is a key technique for 6G
wireless networks.

Figure 10: Error probability $P_{e}$ for different transmit powers with
$J=12$, $\bar{b}=96$, $\mathcal{Z}=32$, $L=100$, $u=3,200$, and $M=300$.

Fig. 10 plots the error probability versus the transmit power. As expected,
increasing the transmit power improves the performance appreciably and the CAD
algorithm achieves much better performance than that of the ML-based multi-
cell method. It is also seen that the ML-based multi-cell method is not as
efficient as the proposed cooperative detection method in accommodating more
active devices $K$, which is a necessary requirement for practical IoT
applications. The simulation result clearly illustrates the advantage of
cooperative strategies which stems from that the proposed scheme not only
exploits the covariance of local received signal, but also explores the
jointly sparse structure and similarity of device state vectors among
different APs to enhance the detection performance.

##  VII Conclusion

This paper designed a grant-free cooperative random access framework for
sourced and unsourced random access in 6G cell-free wireless networks based on
the covariance of the received signals. A cooperative strategy was proposed
for activity detection which involves cooperation among adjacent APs in order
to exploit joint similarity and sparsity. The developed high-accuracy and low-
complexity CAD algorithm handles non-differentiable penalty terms and is shown
to converge to the optimality for adaptive step-sizes. Finally, simulation
results were presented to illustrate the benefit of cooperation and the
proposed algorithm for sourced and unsourced random access can both achieve
near-optimal activity detection performance.

##  Appendix A The Proof of Lemma 1

Let us elaborate on the second moments of
$\bm{\zeta}_{b}^{t}-\bm{\zeta}_{b}^{*}$ for any $\eta_{b}^{t}\leq\eta_{b,0}$
for some constant $\eta_{b,0}$.

$\displaystyle\|\bm{\zeta}_{b}^{t}-\bm{\zeta}_{b}^{*}\|_{2}^{2}\\!=\\!\|\bm{\gamma}_{b}^{t}\\!-\\!\bm{\gamma}_{b}^{*}\|_{2}^{2}\\!-\\!2\bar{p}\eta_{b}^{t}\left\langle\triangledown
f(\bm{\gamma}_{b}^{t})\\!-\\!\triangledown
f(\bm{\gamma}_{b}^{*}),\bm{\gamma}_{b}^{t}\\!-\\!\bm{\gamma}_{b}^{*}\right\rangle$
$\displaystyle+\bar{p}^{2}(\eta_{b}^{t})^{2}\|\triangledown
f(\bm{\gamma}_{b}^{t})-\triangledown f(\bm{\gamma}_{b}^{*})\|_{2}^{2}$
$\displaystyle\leq\\!\\!\|\bm{\gamma}_{b}^{t}\\!-\\!\bm{\gamma}_{b}^{*}\|_{2}^{2}\\!-\\!\bar{p}\eta_{b}^{t}(2\\!-\\!\eta_{b,0}\mathcal{L}_{f})\left\langle\triangledown
f(\bm{\gamma}_{b}^{t})\\!-\\!\triangledown
f(\bm{\gamma}_{b}^{*}),\bm{\gamma}_{b}^{t}-\bm{\gamma}_{b}^{*}\right\rangle$
$\displaystyle-\bar{p}\eta_{b}^{t}\eta_{b,0}\mathcal{L}_{f}\left\langle\triangledown
f(\bm{\gamma}_{b}^{t})-\triangledown
f(\bm{\gamma}_{b}^{*}),\bm{\gamma}_{b}^{t}-\bm{\gamma}_{b}^{*}\right\rangle$
$\displaystyle+\bar{p}^{2}\eta_{b}^{t}\eta_{b,0}\|\triangledown
f(\bm{\gamma}_{b}^{t})-\triangledown f(\bm{\gamma}_{b}^{*})\|_{2}^{2}.$ (61)

According to the properties of Bregman divergences $d_{f}(\cdot,\cdot)$ for a
convex and $\mathcal{L}_{f}$-smooth $f(\cdot)$, we have

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\|\triangledown
f(\bm{\gamma}_{b}^{t})\\!-\\!\triangledown
f(\bm{\gamma}_{b}^{*})\|_{2}^{2}\leq\mathcal{L}_{f}\left\langle\triangledown
f(\bm{\gamma}_{b}^{t})\\!-\\!\triangledown
f(\bm{\gamma}_{b}^{*}),\bm{\gamma}_{b}^{t}\\!-\\!\bm{\gamma}_{b}^{*}\right\rangle,~{}~{}~{}~{}~{}$
(62)

and

$\displaystyle\left\langle\triangledown f(\bm{\gamma}_{b}^{t})-\triangledown
f(\bm{\gamma}_{b}^{*}),\bm{\gamma}_{b}^{t}-\bm{\gamma}_{b}^{*}\right\rangle\geq
d_{f}(\bm{\gamma}_{b}^{t},\bm{\gamma}_{b}^{*}).$ (63)

Substituting (62) into (A) whose last two terms will vanish

$\displaystyle\\!\\!\\!\\!\\!\\!\|\bm{\zeta}_{b}^{t}-\bm{\zeta}_{b}^{*}\|_{2}^{2}\leq\|\bm{\gamma}_{b}^{t}-\bm{\gamma}_{b}^{*}\|_{2}^{2}$
$\displaystyle\\!\\!\\!\\!\\!\\!-\bar{p}\eta_{b}^{t}(2-\eta_{b,0}\mathcal{L}_{f})\left\langle\triangledown
f(\bm{\gamma}_{b}^{t})-\triangledown
f(\bm{\gamma}_{b}^{*}),\bm{\gamma}_{b}^{t}-\bm{\gamma}_{b}^{*}\right\rangle.~{}~{}~{}$
(64)

As a consequence, by substituting (63) into (A), we have

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\|\bm{\zeta}_{b}^{t}-\bm{\zeta}_{b}^{*}\|_{2}^{2}\leq\|\bm{\gamma}_{b}^{t}-\bm{\gamma}_{b}^{*}\|_{2}^{2}-\bar{p}\eta_{b}^{t}(2-\eta_{b,0}\mathcal{L}_{f})d_{f}(\bm{\gamma}_{b}^{t},\bm{\gamma}_{b}^{*}).$
(65)

Inspecting expression (65), we observe that the convergence of this recursion
requires $2-\eta_{b,0}\mathcal{L}_{f}>0$, thereby yielding
$\eta_{b,0}<\frac{2}{\mathcal{L}_{f}}$. According to the definition of
$\eta_{b}^{t}$ in (III-C) and the definition of smoothness of $f(\cdot)$, we
have

$\frac{1}{\mathcal{L}_{f}+\epsilon}<\eta_{b}^{t}<\frac{1}{\epsilon}.$ (66)

Altogether these results, we conclude Lemma 1.

##  Appendix B The Proof of Lemma 2

From the definition of the proximal operator in (11), it holds that

$\displaystyle\\!\\!\\!\\!\\!\\!\mathbf{z}_{b}^{t}=\text{prox}_{\beta\eta_{b}^{t}g}(\bm{\gamma}_{b}^{*}-\eta_{b}^{t}\widetilde{\bigtriangledown
f}(\bm{\gamma}_{b}^{*})-\tau\eta_{b}^{t}\mathbf{x}_{b}^{*})=\mathop{\text{missing}}{argmin}\limits_{\mathbf{y}}\left\\{\beta\eta_{b}^{t}g(\mathbf{y})\right.$
$\displaystyle\\!\\!\\!\\!\\!\\!\left.~{}~{}~{}~{}~{}~{}+\frac{1}{2}\|\mathbf{y}-(\bm{\gamma}_{b}^{*}-\eta_{b}^{t}\widetilde{\bigtriangledown
f}(\bm{\gamma}_{b}^{*})-\tau\eta_{b}^{t}\mathbf{x}_{b}^{*})\|_{2}^{2}\right\\}.$
(67)

Since the function $g(\bm{\gamma}_{b})$ defined in (6) is convex, we have

$\displaystyle\mathbf{0}\in\mathbf{z}_{b}^{t}-(\bm{\gamma}_{b}^{*}-\eta_{b}^{t}\widetilde{\bigtriangledown
f}(\bm{\gamma}_{b}^{*})-\tau\eta_{b}^{t}\mathbf{x}_{b}^{*})+\beta\eta_{b}^{t}\partial
g(\mathbf{z}_{b}^{t}).$ (68)

Combining (68) and (44) leads to $\bm{\gamma}_{b}^{*}=\mathbf{z}_{b}^{*}$.
Therefore, the relation (44) implies that equation (42) holds. Similarly,
equation (14) implies

$\displaystyle\bm{\gamma}_{b}^{t}=\text{prox}_{\tau\eta_{b}^{l,t}\Psi_{l}}(\mathbf{z}_{b}^{*}+\tau\eta_{b}^{l,t}\mathbf{x}_{b}^{l,*})$
$\displaystyle=\mathop{\text{missing}}{argmin}\limits_{\mathbf{y}}\left\\{\tau\eta_{b}^{l,t}\Psi_{l}(\mathbf{y})+\frac{1}{2}\|\mathbf{y}-(\mathbf{z}_{b}^{*}+\tau\eta_{b}^{l,t}\mathbf{x}_{b}^{l,*})\|_{2}^{2}\right\\}.$
(69)

There exists a unique solution $\bm{\gamma}_{b}^{*}$ for this problem due to
the convexity of each $\Psi_{l}$, i.e.,

$\mathbf{0}\in\bm{\gamma}_{b}^{t}-(\mathbf{z}_{b}^{*}+\tau\eta_{b}^{l,t}\mathbf{x}_{b}^{l,*})+\tau\eta_{b}^{l,t}\partial\Psi_{l}(\mathbf{x}_{b}^{l,t})).$
(70)

Using this bound and the relation
$\mathbf{x}_{b}^{l,*}\in\partial\Psi_{l}(\bm{\gamma}_{b}^{*})$, it can be
easily verified that $\bm{\gamma}_{b}^{t}=\bm{\gamma}_{b}^{*}$. Therefore,
there exists a vector $\bm{\gamma}_{b}^{*}$ satisfying equation (43). This
completes the proof.

## References

  * [1] X. Shao, X. Chen, D. W. K. Ng, C. Zhong, and Z. Zhang, “Covariance-based cooperative activity detection for massive grant-free random access,” _IEEE Global Commun. Conf. (GLOBECOM)_ , Taipei, Taiwan, Dec. 2020, pp. 1-6.
  * [2] W. Saad, M. Bennis, and M. Chen, “A vision of 6G wireless systems: Applications, trends, technologies, and open research problems,” _IEEE Network_ , vol. 34, no. 3, pp. 134-142, May 2020.
  * [3] Z. Zhang et al., “6G wireless networks: Vision, requirements, architecture, and key technologies,” _IEEE Veh. Techno. Mag.,_ vol. 14, no. 3, pp. 28-41, Sept. 2019,
  * [4] K. B. Letaief, W. Chen, Y. Shi, J. Zhang, and Y. A. Zhang, “The roadmap to 6G: AI empowered wireless networks,” _IEEE Commun. Mag.,_ vol. 57, no. 8, pp. 84-90, Aug. 2019.
  * [5] K. Senel and E. G. Larsson, “Grant-free massive MTC-enabled massive MIMO: A compressive sensing approach,” _IEEE Trans. Commun._ , vol. 66, pp. 6164-6175, Aug. 2018.
  * [6] X. Chen, D. W. K. Ng, W. Yu, E. G. Larsson, N. Al-Dhahir, and R. Schober, “Massive access for 5G and beyond,” _IEEE J. Sel. Areas Commun._ , vol. PP, no. 99, pp. 1-1, Aug. 2020.
  * [7] V. W. S. Wong, R. Schober, D. W. K. Ng, and L.-C. Wang, _Key Technologies for 5G Wireless Systems_. Cambridge, U.K.: Cambridge Univ. Press, 2017.
  * [8] Z. Chen, F. Sohrabi, and W. Yu, “Sparse activity detection for massive connectivity,” _IEEE Trans. Signal Process._ , vol. 66, no. 7, pp. 1890-1904, Apr. 2018.
  * [9] L. Liu and W. Yu, “Massive connectivity with massive MIMO-Part I: Device activity detection and channel estimation,” _IEEE Trans. Signal Process._ , vol. 66, no. 11, pp. 2933-2946, Jun. 2018.
  * [10] M. Ke, Z. Gao, Y. Wu, X. Gao, and R. Schober, “Compressive sensing-based adaptive active user detection and channel estimation: Massive access meets massive MIMO,” _IEEE Trans. Signal Process._ , vol. 68, pp. 764-779, Jan. 2020.
  * [11] Y. Cheng, L. Liu, and L. Ping, “Orthogonal AMP for massive access in channels with spatial and temporal correlations,” _IEEE J. Sel. Areas Commun._ , vol. PP, no. 99, pp. 1-1, Aug. 2020.
  * [12] X. Shao, X. Chen, and R. Jia, “A dimension reduction-based joint activity detection and channel estimation algorithm for massive access,” _IEEE Trans. Signal Process._ , vol. 68, pp. 420-435, Jan. 2020.
  * [13] S. Haghighatshoar, P. Jung, and G. Caire, “A new scaling law for activity detection in massive MIMO systems,” Jun. 2018. [Online]: arXiv:1803.02288.
  * [14] J. Dong, J. Zhang, Y. Shi, and J. H. Wang, “Faster activity and data detection in massive random access: A multi-armed bandit approach,” Jan. 2020. arXiv preprint arXiv:2001.10237.
  * [15] S. Haghighatshoar, P. Jung, and G. Caire, “Improved scaling law for activity detection in massive MIMO systems,” in _Proc. IEEE Int. Symp. Inf. Theory_ , Colorado, USA, Jun. 2018. pp. 381-385.
  * [16] L. Cheng, L. Liu, and S. Cui, “A covariance-based user activity detection and channel estimation approach with novel pilot design,” _IEEE Intern. Workshop Signal Process. Advances Wireless Commun. (SPAWC)_ , Atlanta, GA, USA, 2020, pp. 1-5.
  * [17] X. Shao, X. Chen, C. Zhong, J. Zhao, and Z. Zhang, “A unified design of massive access for cellular Internet of Things,” _IEEE Internet of Things J._ , vol. 6. no. 2, pp. 3934-3947, Apr. 2019.
  * [18] Y. Polyanskiy, “A perspective on massive random-access,” in _Proc. IEEE Int. Symp. Inf. Theory (ISIT)_ , Aachen, Germany, Jun. 2017, pp. 2523-2527.
  * [19] V. K. Amalladinne, A. Vem, D. K. Soma, K. R. Narayanan, and J. F. Chamberland, “A coupled compressive sensing scheme for uncoordinated multiple access,” Sept. 2018. [Online]. Available: http://arxiv.org/abs/1809.04745.
  * [20] J. Zhang, E. Björnson, M. Matthaiou, D. W. K. Ng, H. Yang and D. J. Love, “Prospective multiple antenna technologies for beyond 5G,” in _IEEE J. Sel. Areas Commun._ , vol. 38, no. 8, pp. 1637-1660, Aug. 2020.
  * [21] A. Fengler, G. Caire, P. Jung, and S. Haghighatshoar, “Massive MIMO unsourced random access,” Jan. 2019, [Online]: http://arxiv.org/abs/1901.00828.
  * [22] Z. Chen, F. Sohrabi, and W. Yu, “Multi-cell sparse activity detection for massive random access: Massive MIMO versus cooperative MIMO,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 8, pp. 1558-2248, Aug. 2019.
  * [23] P. Yang, D. Guo, and H. Yang, “Massive access in multi-cell wireless networks using Reed-Muller codes,” May. 2020, arXiv preprint arXiv:2003.11568.
  * [24] G. Interdonato, E. Bjornson, H. Q. Ngo, et al. “Ubiquitous cell-free massive MIMO communications,” _EURASIP J. on Wireless Commun. and Net._ , vol. 1, p. 197, Aug. 2019.
  * [25] H. Q. Ngo, A. Ashikhmin, H. Yang, E. G. Larsson, and T. L. Marzetta, “Cell-free massive MIMO versus small cells,” _IEEE Trans. Wirel. Commun._ , vol. 16, no. 3, pp. 1834-1850, Mar. 2017.
  * [26] M. Ke, Z. Gao, Y. Wu, X. Gao, and K. Wong, “Massive access in cell-free massive MIMO-based Internet of Things: cloud computing and edge computing paradigms,” _IEEE J. Sel. Areas Commun._ , vol. PP, no. 99, pp. 1-1, Aug. 2020.
  * [27] O. Ozdogan, E. Bjornson, and J. Zhang, “Performance of cell-free massive MIMO with Rician fading and phase shifts,” _IEEE Trans. on Wireless Commun._ , vol. 18, no. 11, pp. 5299-5315, Nov. 2019.
  * [28] C. Li, S. Huang, Y. Liu, and Z. Zhang, “Distributed jointly sparse multitask learning over networks,” _IEEE Trans. on Cybern._ , vol. 48, no. 1, pp. 151-164, Jan. 2018
  * [29] A. Jalali, P. Ravikumar, and S. Sanghavi, “A dirty model for multiple sparse regression,” _IEEE Trans. Inf. Theory,_ vol. 59, no. 12, pp. 7947-7968, 2013.
  * [30] A. Wiesel, “Geodesic convexity and covariance estimation,” _IEEE Trans. Signal Process._ , vol. 60, no. 12, pp. 6182, Dec. 2012.
  * [31] S. Khanna and C. R. Murthy, “Decentralized joint-sparse signal recovery: a sparse bayesian learning approach,” _IEEE Trans. on Signal and Inf. Process. over Netw._ , vol. 3, no. 1, pp. 29-45, Mar. 2017
  * [32] N. Parikh and S. Boyd, “Proximal algorithms,” _Found. Trends Optim._ , vol. 1, no. 3, pp. 123-231, 2013.
  * [33] J. Sherman and W. J. Morrison, “Adjustment of an inverse matrix corresponding to a change in one element of a given matrix,” _The Annals of Mathematical Statistics_ , vol. 21, no. 1, pp. 124-127, 1950.
  * [34] J. Chen, C. Richard, and A. H. Sayed, “Diffusion LMS for clustered multitask networks,” _IEEE Intern. Conf. Acous., Speech, Signal Proces. (ICASSP)_ , Florence, May. 2014, pp. 5487-5491.
  * [35] A. G. GAkritas, E. K. Akritas, and G. I. Malaschonok, “Various proofs of Sylvester’s (determinant) identity, _Mathematics and Computers in Simulation,_ vol. 42, no. 4-6, pp. 585-593, 1996.
  * [36] J. Barzilai and J. M. Borwein, “Two-point step size gradient methods,” _IMA J. Numer. Anal._ , vol. 8, pp. 141-148, 1988.
  * [37] S. Boyd and L. Vandenberghe, “Convex Optimization,” Cambridge, U.K.: Cambridge Univ. Press, 2004.
  * [38] T. Hofmann, A. Lucchi, S. L. Julien, and B. McWilliams, “Variance reduced stochastic gradient descent with neighbors,” _Advances in Neural Information Processing Systems (NIPS)_ , pp. 2296-2304, 2015.
  * [39] A. Defazio, “A simple practical accelerated method for finite sums,” _Advances in Neural Information Processing Systems (NIPS)_ , 2016.
  * [40] N. D. Vanli, M. G. zbalaban, and A. Ozdaglar, “A Simple proof for the iteration complexity of the proximal gradient algorithm,” _Neural Information Processing Systems (NIPS)_ , Barcelona, Spain, 2016.