# Device Activity Detection and Channel Estimation for Millimeter-Wave Massive
MIMO

Yinchuan Li, _Member_ , _IEEE_ , Yuancheng Zhan, Le Zheng, _Senior Member_ ,
_IEEE_ , Xiaodong Wang, _Fellow_ , _IEEE_ Y. Li is with the Electrical
Engineering Department, Columbia University, New York, NY 10027, USA (e-mail:
yinchuan.li.cn@gmail.com). Y. Zhan is with the School of Electrical and
Electronic Engineering, Nanyang Technological University, Singapore 639798
(e-mail: zhan0530@e.ntu.edu.sg). Le Zheng is with the Radar Research
Laboratory, School of Information and Electronics, Beijing Institute of
Technology, Beijing 100081, China. He is also with the Chongqing Innovation
Center, Beijing Institute of Technology, Chongqing 401120, China
(e-mail:le.zheng.cn@gmail.com). X. Wang is with the Electrical Engineering
Department, Columbia University, New York, NY 10027, USA (e-mail:
wangx@ee.columbia.edu).

###### Abstract

Millimeter-Wave Massive MIMO is important for beyond 5G or 6G wireless
communication networks. The goal of this paper is to establish successful
communication between the cellular base stations and devices, focusing on the
problem of joint user activity detection and channel estimation. Different
from traditional compressed sensing (CS) methods that only use the sparsity of
user activities, we develop several Approximate Message Passing (AMP) based CS
algorithms by exploiting the sparsity of user activities and mmWave channels.
First, a group soft-thresholding AMP is presented to utilize only the user
activity sparsity. Second, a hard-thresholding AMP is proposed based on the
on-grid CS approach. Third, a super-resolution AMP algorithm is proposed based
on atomic norm, in which a greedy method is proposed as a super-resolution
denoiser. And we smooth the denoiser based on Monte Carlo sampling to have
Lipschitz continuity and present state evolution results. Extensive simulation
results show that the proposed method outperforms the previous state-of-the-
art methods.

###### Index Terms:

Compressed sensing, atomic norm, approximate message passing, millimeter-Wave,
massive connectivity, state evolution, machine-type communications, massive
multiple-input multiple-output (MIMO).

##  I Introduction

###  I-A Motivation

Massive connectivity is common in beyond 5G (B5G) or 6G wireless communication
networks, with profound implications for Internet of Things (IoT) and Machine
Type Communications (MTC). In this case, each cellular base station (BS) is
connected to a large number of devices (on the order of $10^{4}$ to $10^{6}$).
However, only a few of these large numbers of devices are usually active.
Since to save power, the device is in a sleep state most of the time, and will
not be activated when it is not triggered [1, 2, 3]. For such networks,
accurate user activity detection and channel estimation are crucial to
establish successful communication between BS and devices.

###  I-B Prior Works

At beginning, the traditional cellular networks is based on arranging active
devices into time or frequency domains. The drawback is that the overhead is
huge, because a large number of occasionally active users need to be scheduled
on a separate control channel. To address this problem, a method based on
random access protocol [4, 5] is proposed, the main idea of this kinds of
contention-based methods is that the orthogonal signature sequence sent to the
BS is randomly selects by each active user, if the selected preamble is not
used by another user, then the connection is established. However, for a large
number of devices, collisions are unavoidable, so contention resolution is
required. Therefore, unlicensed non-orthogonal user access schemes are very
popular. The main function of the non-orthogonal pilot sequence is to improve
the detection performance under the condition of a large number of devices and
a limited coherence time of the channel. Afterwards, the base station
completes user activity detection and channel estimation based on the pilot
sequences sent by active users to the BS simultaneously.

Furthermore, the sparsity in user activity pattern in these systems can be
exploited to formulate compressed sensing (CS) problems [6, 7, 8, 9, 10, 11,
12, 13]. In [6, 7], under the condition that the channel state information
(CSI) is well performed at the BS, a user activity and data detection
algorithm is proposed based on sparse structures.

For the case that CSI is not available at the BS, a modified Bayesian
compressed sensing method is proposed in [8] for jointly user activity
detection and channel estimation. Aiming at similar problems, [9] based on
compressed sensing techniques [10, 11, 12], adopts methods such as basis
pursuit to denoise and jointly decode information in orthogonal frequency
division multiplexing (OFDM) systems. An algorithm [14] based on Bayesian
learning is developed, utilizing proposed hyper-priors to capture structural
signal characteristics and appropriate approximations to facilitate algorithm
derivation.

Unfortunately, none of the above works has conducted strict performance
analysis for large-scale connectivity non-orthogonal multiple access schemes.
With the development of the approximate message passing (AMP) algorithm [15],
the state evolution (SE) has shown its benefit for performance analysis in
compressed sensing problem. Therefore, the AMP algorithm is recently used for
activity detection and channel estimation problems [13, 16]. The advantage of
the AMP algorithm is that it can perform state evolution analysis, and
relevant conclusions can be used to model missed detection and false positive
probability. However, the analysis in [13, 16] is quite involved. To improve
this, [1, 2] simplifies the related characterization under the assumption of a
large number of BS antennas. Furthermore, an important conclusion is drawn
that the channel estimation error is much more restrictive for MIMO mechanisms
than the detection error. [17, 18] propose model-based neural network
architectures for sporadic user detection and channel estimation in massive
machine-type communications. The deep unfolding framework is applied to design
customized neural network layers through the unrolling of two iterative
optimization algorithms. However, the ability of the networks to adapt to the
potential mismatch between training or test data and signal models on which
they are based is yet to be examined. Additionally, an algorithm is necessary
to minimize the number of learnable parameters.

Millimeter wave massive MIMO is very important for B5G/6G wireless
communication, that is, the number of paths is sparse compared to the number
of base station antennas, this sparsity in mmWave channels can actually be
better than the sparsity in user activity. In this paper, we propose several
AMP algorithms to exploit the sparsity of mmWave channels for device activity
detection and channel estimation.

###  I-C Main Contributions

In this paper, we consider the problem of user activity detection and channel
estimation in mmWave communication systems with massive MIMO, where the
channel is a superposition of multiple complex sinusoids.

We setup the compressed sensing problem by exploiting two kinds of sparsities:
1) the user activity is sparse since only a fraction of devices are active; 2)
the mmWave channel is sparse since the number of paths is usually much less
than the number of BS antennas. To solve this compressed sensing problem, four
methods are proposed based on AMP, on-grid approximation, soft-thresholding,
hard-thresholding, atomic norm (AN) [19] and greedy techniques. Specifically,
we extend the approach in [1] to the scenario of mmWave channel estimation.
The mmWave channel exhibits two types of sparsity: channel sparsity and user
activity sparsity. While [1] considers only one type of sparsity, our work
models both types of sparsity simultaneously, providing a more comprehensive
and innovative approach. In addition, the design of the denoisers in our
framework is quite original, except for the first group-soft denoiser.
Besides, the combination of the atomic norm with AMP is an innovative aspect
of our work, which leads to a fast solver and it leads to significant
performance improvements. Our main contributions are listed as follows.

1) A group soft-thresholding AMP is presented to exploit only the sparsity in
user activity. Since in mmWave channel, the channel is not Rayleigh fading and
the prior distribution may not be the Bernoulli Gaussian distribution [1]. The
minimum mean-squared error (MMSE) denoiser does not apply. We hence take the
group soft-thresholding function as the denoiser, which also has lower
computationally load than the MMSE denoiser; 2) To further utilize the
sparsity in the mmWave channel, we approximately define an overcomplete
dictionary to construct the sparsity in mmWave channel. And then a hard-
thresholding AMP is proposed based on the on-grid CS approach [20, 21].
Although the sparsity in mmWave channel is exploited, the frequencies in
channel may not fall onto discrete grids [22, 23, 24, 25] and hence hard-
thresholding AMP suffers from performance loss; 3)

To overcome this difficulty, we exploit the atomic norm technique to address
the off-grid problem to obtain super-resolution performance. Then, a direct
SDP solver is proposed, which, however, has high computational complexity; 4)
As a result, we finally combine the AMP and atomic norm techniques [19, 26,
27, 28] for user detection and channel estimation, which has moderate
computational complexity. A greedy method is proposed as the super-resolution
denoiser in the AMP, which includes the initialization, selection, local-
optimization, least-squares and residual update steps.

In addition, the state evolution of the above proposed methods is analyzed. To
ensure that the state evolution is accurate, we present methods for smoothing
the denoisers to enable they are Lipschitz continuous [29]. In particular, the
Monte Carlo sampling [30] is used to smooth the greedy denoiser and to
calculate the derivative of the denoiser, which is required in the AMP
iterations. And we prove that the derivative obtained by the Monte Carlo
sampling is accurate as long as the number of Monte Carlo simulation goes to
infinity. At last, simulation results show the advantage of the proposed
methods for user activity detection and channel estimation in massive MIMO
mmWave communication systems.

###  I-D Organization

The remainder of the paper is organized as follows. We present the background
in Section II. In Section III, we propose the group soft-thresholding AMP and
the hard-thresholding AMP algorithms, and present the state evolution analysis
results. In Section IV, we propose the SDP direct solver and the combined
approach based on the AMP and atomic norm, which can overcome the off-grid
problem. Besides, the smoothed denoiser is also presented. Simulation results
are presented in Section V. Finally, Section VI concludes the paper.

##  II Background

###  II-A Problem Formulation

Consider the uplink of a cellular network in which the base station is
equipped with $M$ antennas and contains $N$ users each equipped with an
antenna. Let $\bm{h}_{n}\in{\mathbb{C}}^{M\times 1}$ be the $n$-th user
complex uplink channel vector to the BS. We adopt a block-fading model where
$\bm{h}_{n}$ remains constant within a coherence time block, and varies
independently from block to block. For the mm-Wave channel [31, 32],
$\bm{h}_{n}$ can be modeled as

$\displaystyle\bm{h}_{n}=\sum_{\ell=1}^{L_{n}}\bar{c}_{n,\ell}\bm{a}(f_{n,\ell}),$
(1)

where $\bar{c}_{n,\ell}\in\mathbb{C}$ denotes the $\ell$-th multipath gain;
$L_{n}$ denotes the number of paths; the vector
$\bm{a}(f)\in\mathbb{C}^{M\times 1}$ denotes the steering responses of the
receive array, given by

$\displaystyle\bm{a}(f)=\frac{1}{\sqrt{M}}[1,e^{j2\pi
f},...,,e^{j2\pi(M-1)f}]^{T},$ (2)

where $f=\frac{d\sin\theta}{2\lambda_{w}}$ with $\theta$, $d$ and
$\lambda_{w}$ denoting the angle of arrival, the separation of adjacent
antenna elements and the wavelength of the transmitted signal, respectively.

Within each coherence block each user is active with probability $\epsilon$.
Denote

$\displaystyle{\alpha_{n}}=\left\\{\begin{array}[]{l}1,{\text{ if user
}}n{\text{ is active,}}\\\ 0,{\text{ otherwise}},\end{array}\right.$ (5)

so that ${\rm Pr}(\alpha_{n}=1)=\epsilon$, ${\rm
Pr}(\alpha_{n}=0)=1-\epsilon$.

Denote

$\displaystyle\bm{u}_{n}=[u_{n,1},u_{n,2},...,u_{n,Q}]^{T}\in\mathbb{C}^{Q\times
1}$ (6)

as the unique data symbols that is known by the BS.

Since we consider massive connection scenarios, the device number is usually
greater than the length of the pilot sequence, i.e., $Q<N$. Therefore, it is
impossible to assign mutually orthogonal sequences to all users.

Denote $\bm{Y}\in\mathbb{C}^{Q\times M}$ as the matrix of the received signal
at the BS across $M$ antennas over $Q$ pilot symbols, given by

$\displaystyle\bm{Y}=\sqrt{\rho}\sum_{n=1}^{N}\alpha_{n}\bm{u}_{n}\bm{h}_{n}^{T}+\bm{Z},$
(7)

where $\rho$ is the pilot transmit power and the noise matrix
$\bm{Z}\in\mathbb{C}^{Q\times M}$ contains i.i.d.
${\mathcal{C}N}(\bm{0},\sigma_{w}^{2})$ samples. Define

$\displaystyle\bm{U}=$
$\displaystyle~{}[\bm{u}_{1},\bm{u}_{2},...,\bm{u}_{N}]\in\mathbb{C}^{Q\times
N},$ (8) $\displaystyle\bm{X}=$
$\displaystyle~{}[\bm{x}_{1},\bm{x}_{2},...,\bm{x}_{N}]^{T}\in\mathbb{C}^{N\times
M},$ (9)

where

$\displaystyle\bm{x}_{n}=$
$\displaystyle~{}\alpha_{n}\sqrt{\rho}\bm{h}_{n}=\alpha_{n}\sqrt{\rho}\sum_{\ell=1}^{L_{n}}\bar{c}_{n,\ell}\bm{a}(f_{n,\ell})$
$\displaystyle=$
$\displaystyle~{}\alpha_{n}\sum_{\ell=1}^{L_{n}}c_{n,\ell}\bm{a}(f_{n,\ell})$
(10)

with $c_{n,\ell}=\sqrt{\rho}\bar{c}_{n,\ell}$. To simplify the notation, we
let $\bm{h}_{n}=\sum_{\ell=1}^{L_{n}}c_{n,\ell}\bm{a}(f_{n,\ell})$ in the
following. Then, (7) can be rewritten as

$\displaystyle\bm{Y}=\bm{U}\bm{X}+\bm{Z}.$ (11)

Our goal is to estimate the user activities (i.e., to obtain
$\\{{\alpha_{n}\\}}$ ) and also to estimate the corresponding channels (i.e.,
to determine $\\{\bm{h}_{n}:\alpha_{n}=1\\}$).

###  II-B The AMP Solution and the Traditional Denoiser

Denote

$\displaystyle\bm{X}^{t}=$
$\displaystyle~{}[\bm{x}_{1}^{t},...,\bm{x}_{N}^{t}]^{T}\in\mathbb{C}^{N\times
M},$ (12) $\displaystyle\bm{R}^{t}=$
$\displaystyle~{}[\bm{r}_{1}^{t},...,\bm{r}_{Q}^{t}]^{T}\in\mathbb{C}^{Q\times
M},$ (13)

as the estimate of $\bm{X}$ and the corresponding residual at iteration $t$,
respectively. By initializing to $\bm{X}^{0}=\bm{0}$ and $\bm{R}^{0}=\bm{Y}$,
the vector AMP [33] algorithm is updated according to iterations:

$\displaystyle\bm{x}_{n}^{t+1}=$
$\displaystyle~{}\eta_{t}((\bm{R}^{t})^{H}\bm{u}_{n}+\bm{x}_{n}^{t}),~{}n=1,2,...,N,$
(14) $\displaystyle\bm{R}^{t+1}=$
$\displaystyle~{}\bm{Y}-\bm{U}\bm{X}^{t+1}+\frac{1}{Q}\bm{R}^{t}\sum_{n=1}^{N}{\eta_{t}^{\prime}}\left((\bm{R}^{t})^{H}\bm{u}_{n}+\bm{x}_{n}^{t}\right),$
(15)

where $\eta_{t}(\cdot):\mathbb{C}^{M\times 1}\to\mathbb{C}^{M\times 1}$ is an
appropriately designed non-linear function known as denoiser [34, 35, 36, 37],
and $\eta_{t}^{\prime}(\cdot)$ is the first order derivative of
$\eta_{t}(\cdot)$ with respect to the input argument.

In [1], the channel vector is modeled as an independent quasi-static flat
fading component, i.e., $\bm{h}_{n}$ in (1) is replaced by
$\bm{h}_{n}\sim\mathcal{CN}(0,\beta_{n}\bm{I}_{n})$ with $\beta_{n}$ being the
path-loss and shadowing component. Then, based on the conditional expectation,
a MMSE denoiser is given by

$\displaystyle\eta_{t}(\bm{\tilde{x}}_{n}^{t})=\mathbb{E}\left(\bm{x}_{n}|\bm{\tilde{x}}_{n}^{t}\right)$
(16)

with $\bm{\tilde{x}}_{n}^{t}=(\bm{R}^{t})^{H}\bm{u}_{n}+\bm{x}_{n}^{t}$, which
is specifically designed for the Rayleigh fading channel. For MMSE denoiser,
once $\hat{\bm{X}}$ is obtained, the activity user can be detected and the
channel can be estimated according to Definition 1 in [1].

However, in our case the channel is not Rayleigh fading and the prior
distribution of ${\bm{x}}_{n}$ may not be the Bernoulli Gaussian distribution
[1]. Hence, the MMSE denoiser does not apply. Moreover, the MMSE denoiser can
be computationally complex, so in practice some simplified form [36] or
approximated denoiser [37] can be used instead.

##  III Proposed On-grid Solutions Based on AMP

An effective approach for user detection and channel estimation is to exploit
the sparsity of $\bm{X}$. Since $\bm{x}_{n}$’s are zero for those users are
not active, $\bm{X}$ is row sparse. Moreover, each $\bm{x}_{n}$ is also sparse
since there are sparse nature in $\bm{h}_{n}$, i.e., for the channel of mmWave
communication, the number of paths is usually much smaller than the number of
BS antennas, i.e., $L_{n}\ll M$ for $n=1,2,...,N$. In this section, we develop
on-grid CS approaches based on the AMP algorithm by utilizing the user and
channel sparsities.

###  III-A Group Soft-thresholding AMP Based on User Activity Sparsity

Following [1], we use the vector AMP algorithm (14)-(15) for mmWave
communication by exploiting the user sparsity. Once $\hat{\bm{X}}$ is
obtained, we compare $\|\hat{\bm{x}}_{n}\|_{1}$ with a threshold $\varsigma$
to detect the activity users, i.e., a user is determined as activity if
$\|\hat{\bm{x}}_{n}\|_{1}>\varsigma$. The threshold $\varsigma$ is chosen to
have a desired probabilities of miss detection and false alarm performance.
Furthermore, when a device $n$ is detected to be active, its channel is
estimated to be $\hat{\bm{h}}_{n}=\hat{\bm{x}}_{n}$. Later, we will propose
several denoisers that take advantage of both user sparsity and channel
sparsity. Before that, we introduce another denoiser based on group soft-
thresholding that only uses user sparsity, which is faster than MMSE denoiser
in [1].

####  III-A1 Group soft-thresholding denoiser

Since MMSE in [1] involves inversion operation, its computational complexity
is $\mathcal{O}(M^{3})$. To be more efficient, we estimate $\bm{X}$ by
formulating a group LASSO problem, i.e., we aim to find a group sparse
$\bm{\hat{X}}$ that can minimize $\|\bm{Y}-\bm{U}\bm{\hat{X}}\|_{F}^{2}$ by
solving

$\displaystyle\bm{\hat{X}}=\arg\min_{\bm{X}\in{\mathbb{C}}^{N\times
M}}\frac{1}{2}\|\bm{Y}-\bm{U}\bm{X}\|_{F}^{2}+\gamma\|\bm{X}\|_{2,1},$ (17)

where $\|\bm{X}\|_{2,1}\triangleq\sum_{n=1}^{N}\|\bm{x}_{n}\|_{2}$, $\gamma$
is the regularization parameter related to $N$, $Q$ and $\bm{X}$. We use the
vector AMP (14)-(15) to solve (17), where it is straightfoward to use the
group soft-thresholding function as the denoiser. In particular, for any input
$\bm{\tilde{x}}$ we have [38]

$\displaystyle\eta_{t}^{\rm soft}(\bm{\tilde{x}})=$
$\displaystyle~{}\arg\min_{\bm{x}\in{\mathbb{C}}^{M\times
1}}\frac{1}{2}\|\bm{\tilde{x}}-\bm{x}\|_{2}^{2}+\lambda_{t}\|\bm{x}\|_{2}$
$\displaystyle=$
$\displaystyle~{}\frac{\bm{\tilde{x}}}{\|\bm{\tilde{x}}\|_{2}}\left(\|\bm{\tilde{x}}\|_{2}-\lambda_{t}\right)\mathbb{I}(\|\bm{\tilde{x}}\|_{2}\geq\lambda_{t}),$
(18)

where $\mathbb{I}(\cdot)$ denotes the indicator function. Obviously, the
computational complexity of the group soft-thresholding denoiser is
$\mathcal{O}(M)$, which is much smaller than that of the MMSE denoiser.

Note that for complex value function, the first order derivative
$\eta_{t}^{\prime}(\bm{\tilde{x}})=\frac{\partial\eta_{t}(\bm{\tilde{x}})}{\partial\bm{\tilde{x}}}$
is calculated by

$\displaystyle\frac{\partial\eta_{t}(\bm{\tilde{x}})}{\partial\bm{\tilde{x}}}=$
$\displaystyle~{}\frac{1}{2}\left(\frac{\partial\eta_{t}(\bm{\tilde{x}})}{\partial\bm{\tilde{x}}^{\Re}}-i\frac{\partial\eta_{t}(\bm{\tilde{x}})}{\partial\bm{\tilde{x}}^{\Im}}\right)$
$\displaystyle=$
$\displaystyle~{}\frac{1}{2}\left(\frac{\partial\eta_{t}^{\Re}(\bm{\tilde{x}})}{\partial\bm{\tilde{x}}^{\Re}}+\frac{\partial\eta_{t}^{\Im}(\bm{\tilde{x}})}{\partial\bm{\tilde{x}}^{\Im}}\right)$
$\displaystyle+\frac{i}{2}\left(\frac{\partial\eta_{t}^{\Re}(\bm{\tilde{x}})}{\partial\bm{\tilde{x}}^{\Im}}-\frac{\partial\eta_{t}^{\Im}(\bm{\tilde{x}})}{\partial\bm{\tilde{x}}^{\Re}}\right)\in{\mathbb{C}}^{M\times
M},$ (19)

where $\eta_{t}^{\Re}(\cdot)$ and $\eta_{t}^{\Im}(\cdot)$ denote the real and
imaginary parts, respectively; and $\bm{\tilde{x}}^{\Re}$ and
$\bm{\tilde{x}}^{\Im}$ are the real and imaginary parts of the input
$\bm{\tilde{x}}\in{\mathbb{C}}^{M\times 1}$, respectively. Hence, we have

$\displaystyle~{}\frac{\partial\eta_{t}^{\rm
soft}(\bm{\tilde{x}})}{\partial\bm{\tilde{x}}}=\frac{\partial(\bm{\tilde{x}}-\lambda_{t}{\bm{\tilde{x}}}/{\|\bm{\tilde{x}}\|_{2}})}{\partial\bm{\tilde{x}}}\mathbb{I}(\|\bm{\tilde{x}}\|_{2}\geq\lambda_{t})$
$\displaystyle=$
$\displaystyle~{}\left[{\bm{I}}_{M}-\frac{\lambda_{t}}{2}\left(\frac{\partial\frac{{\bm{\tilde{x}}}^{\Re}+i{\bm{\tilde{x}}}^{\Im}}{{\|\bm{\tilde{x}}\|_{2}}}}{\partial\bm{\tilde{x}}^{\Re}}-i\frac{\partial\frac{{\bm{\tilde{x}}}^{\Re}+i{\bm{\tilde{x}}}^{\Im}}{{\|\bm{\tilde{x}}\|_{2}}}}{\partial\bm{\tilde{x}}^{\Im}}\right)\right]\mathbb{I}(\|\bm{\tilde{x}}\|_{2}\geq\lambda_{t})$
$\displaystyle=$
$\displaystyle~{}\left[{\bm{I}}_{M}-\frac{\lambda_{t}}{\|\bm{\tilde{x}}\|_{2}}{\bm{I}}_{M}+{\lambda_{t}}\frac{{\bm{\tilde{x}}}{\bm{\tilde{x}}}^{H}}{2\|\bm{\tilde{x}}\|_{2}^{3}}\right]\mathbb{I}(\|\bm{\tilde{x}}\|_{2}\geq\lambda_{t}).$
(20)

In practice, it is convenient to set $\lambda_{t}=\tau\sigma_{t}$ in which
$\tau$ is a constant and $\sigma_{t}=\frac{1}{\sqrt{MQ}}\|\bm{R}^{t}\|_{F}$
[38]. In this paper, we name the AMP with the group soft-thresholding as GST-
AMP.

####  III-A2 State evolution

Note that since the ground truth ${\bm{X}}^{\star}$ is unknown, we cannot
directly calculate $\|\hat{\bm{X}}-{\bm{X}}^{\star}\|_{2}^{2}$ to obtain the
estimation performance of (14)-(15) against the iteration number.
Alternatively, we calculate the state evolution to predict the estimation
performance in certain asymptotic regime. A remarkable property of the AMP
algorithm is that when applied to the compressed sensing problem with the
entries of the sensing matrix $\bm{U}$ generated from i.i.d. Gaussian
distribution, its estimation performance in certain asymptotic regime can be
accurately predicted by the so-called state evolution. In particular, when
$Q,K,N\to\infty$, ($K$ is the active users) the ratios converge to positive
constants, i.e., $K/N\to\epsilon$ and $Q/N\to\omega$.

Denote $\bm{\Theta}_{0}=\frac{1}{N}\sum_{n=1}^{N}\bm{x}_{n}\bm{x}_{n}^{H}$,
then a series of state evolution matrices are generated by the following
iterations:

$\displaystyle\bm{\Theta}_{t+1}(\bm{X},\sigma_{w}^{2})=$
$\displaystyle~{}\frac{1}{N}\sum_{n=1}^{N}\mathbb{E}\Bigg{[}\Big{(}\eta_{t+1}\left(\bm{x}_{n}+\bm{\Sigma}_{t}^{1/2}\bm{v}_{n}\right)-\bm{x}_{n}\Big{)}$
$\displaystyle~{}\Big{(}\eta_{t+1}\left(\bm{x}_{n}+\bm{\Sigma}_{t}^{1/2}\bm{v}_{n}\right)-\bm{x}_{n}\Big{)}^{H}\Bigg{]}.$
(21)

where the expectation is with respect to $\bm{v}_{n}\sim{\cal
CN}(\bm{0},\bm{I}_{M})$ for $n=1,2,...,N$; recall that $\sigma_{w}^{2}$ is the
variance of the noise; $\eta_{t}(\cdot)$ is the same as (III-A1) with
$\lambda_{t}$ replaced by $\bar{\lambda}_{t}=\tau\sqrt{{\rm
Tr}(\bm{\Sigma}_{t})}$, where $\rm{Tr}(\cdot)$ is the trace of the input
matrix; and

$\displaystyle\bm{\Sigma}_{t}=\frac{1}{\omega}\bm{\Theta}_{t}(\bm{X},\sigma_{w}^{2})+\sigma_{w}^{2}\bm{I}_{M}$
(22)

is referred to as the state. Based on the standard state evolution analysis
(e.g. Theorem 1 in [34] and Theorem 1 in [35]) when $Q$ and $N$ are large
enough, the state evolution can predict the MSE of the algorithm, i.e.,

$\displaystyle\bm{\Theta}_{t+1}(\bm{X},\sigma_{w}^{2})\approx\frac{1}{N}\sum_{n=1}^{N}(\bm{x}_{n}^{t+1}-\bm{x}_{n})(\bm{x}_{n}^{t+1}-\bm{x}_{n})^{H}.$
(23)

The above equations give a general framework to analyze the performance of the
AMP. Obviously, the smaller
${\text{Tr}}(\bm{\Theta}_{t+1}(\bm{X},\sigma_{w}^{2}))$ is, the smaller the
reconstruction error will be.

###  III-B Hard-thresholding AMP Based on Channel Sparsity

The GST-AMP exploits the sparsity in the user activity pattern, but does not
use the sparsity of the mmWave channel. We next propose a new method Hard-
thresholding AMP that can exploit these two sparsity. Recall that the channel
is sparse in spatial domain, i.e., $L_{n}$ is small compared to $M$. For
better user detection and channel estimation, the channel sparsity is
investigated and an improved AMP is proposed with a different
$\eta_{t}(\cdot)$. Define an overcomplete dictionary

$\displaystyle\bm{A}=[\bm{a}(f_{1}^{\prime}),\bm{a}(f_{2}^{\prime}),...,\bm{a}(f_{\tilde{M}}^{\prime})]\in{\mathbb{C}}^{M\times\tilde{M}},$
(24)

where $\\{f_{m}^{\prime}\\}$ is sets of uniformly spaced frequency points and
$\tilde{M}>M$ is the number of grid points. If the frequencies
$\\{f_{n,\ell}\\}_{1\leq\ell\leq L_{n}}$ in (1) fall on the grids for
$n=1,2,...,N$, let $\bm{\tilde{c}}_{n}$ be the corresponding $L_{n}$-sparse
coefficient vector, i.e., $\bm{h}_{n}=\bm{A}\bm{\tilde{c}}_{n}$ and
$\|\bm{\tilde{c}}_{n}\|_{0}=L_{n}$, then we have
$\bm{A}^{H}{\bm{h}}_{n}=\bm{A}^{H}\bm{A}\bm{\tilde{c}}_{n}=\bm{\tilde{c}}_{n}$
is $L_{n}$-sparse. Further, we have $\bm{A}^{H}\bm{x}$ is $L_{n}$-sparse or
$0$-sparse when $\alpha_{n}=0$. Hence, to enhance the performance of channel
estimation, we denoise $\bm{\tilde{x}}$ by solving

$\displaystyle\arg\min_{\bm{x}}\|\bm{A}^{H}\bm{x}\|_{1},~{}\text{s.t.}~{}\|\bm{x}-\bm{\tilde{x}}\|_{2}^{2}\leq\gamma,$
(25)

where $\gamma$ is a threshold that can be adaptively changed in each
iteration. (25) introduces a way to exploit the sparsity of $\bm{x}$. It is
known that we can use a hard-thresholding function as the denoiser to solve
(25), for any input $\bm{\tilde{x}}$ we have

$\displaystyle\eta_{t}^{\rm
hard}(\bm{\tilde{x}})={\bm{A}}{\mathbb{V}}({\bm{A}}^{H}\bm{\tilde{x}}),$ (26)

where
${\mathbb{V}}(\bm{v})=[{\mathbb{V}}(v_{1}),{\mathbb{V}}(v_{2}),...,{\mathbb{V}}(v_{\tilde{M}})]^{T}$
with ${\bm{v}}\in{\mathbb{C}}^{\tilde{M}\times 1}$ and
${\mathbb{V}}(v_{m})=v_{m}{\mathbb{I}}(|v_{m}|>\lambda_{t})$. In addition, we
have

$\displaystyle\frac{\partial\eta_{t}^{\rm
hard}(\bm{\tilde{x}})}{\partial\bm{\tilde{x}}}=$
$\displaystyle~{}\frac{\partial({\bm{A}}{\mathbb{V}}({\bm{A}}^{H}\bm{\tilde{x}}))}{\partial\bm{\tilde{x}}}={\bm{A}}{\rm
diag}(\bar{\mathbb{V}}({\bm{A}}^{H}\bm{\tilde{x}})){\bm{A}}^{H},$ (27)

where ${\rm diag}(\cdot)$ denotes the diagonal operator, whose output is an
diagonal matrix with diagonal entries are the input vector;
$\bar{\mathbb{V}}(\bm{v})=[\bar{\mathbb{V}}(v_{1}),\bar{\mathbb{V}}(v_{2}),...,\bar{\mathbb{V}}(v_{\tilde{M}})]^{T}$
with ${\bm{v}}\in{\mathbb{C}}^{\tilde{M}\times 1}$ and
$\bar{\mathbb{V}}(v_{m})={\mathbb{I}}(|v_{m}|>\lambda_{t})$. In this paper, we
name the AMP with the hard-thresholding as HT-AMP. Based on (25) and (11), the
algorithm provides the solution to

$\displaystyle\\{\bm{\hat{x}}_{n}\\}_{n=1}^{N}=$
$\displaystyle~{}\mathop{\arg\min}\limits_{\\{\bm{x}_{n}\\}_{n=1}^{N}}\sum_{n=1}^{N}\|\bm{A}^{H}\bm{x}_{n}\|_{1},$
(28) s.t.
$\displaystyle~{}{\left\|\bm{Y}-\sum_{n=1}^{N}\bm{u}_{n}\bm{x}_{n}^{T}\right\|_{F}}\leq\zeta,$

where $\zeta>0$ is a weight factor. The objective in (28) can be understood as
we first utilize the $\ell_{1}$-norm
$\|\bm{A}^{H}\bm{x}_{n}\|_{1}=\|\alpha_{n}\bm{\tilde{c}}_{n}\|_{1}$ to exploit
the channel sparsity, then use $\sum_{n=1}^{N}$ to exploit the user sparsity
since
$\sum_{n=1}^{N}\|\alpha_{n}\bm{\tilde{c}}_{n}\|_{1}=\left\|\bm{[}\|\alpha_{1}\bm{\tilde{c}}_{1}\|_{1},...,\|\alpha_{N}\bm{\tilde{c}}_{N}\|_{1}]^{T}\right\|_{1}$.

The denoiser used within AMP can take on almost any form, and the SE of the
HT-AMP can be obtained via (III-A2)-(23) with $\eta(\cdot)$ function in
(III-A2) replaced by (26). However, the denoiser needs to be Lipschitz
continuous [29], otherwise the SE prediction is not accurate enough. One
simple idea is to “smooth” the denoisers. In particular, we can decompose
${\mathbb{V}}(v)$ to [34]

$\displaystyle{\mathbb{V}}(v)=S_{t}(v)+D_{t}(v),$ (29)

where $S_{t}(v)$ is a weakly differentiable function and can be set as the
Lipschitz continuous function of $v$, while $D_{t}(v)$ is a piece-wise
constant function and is not continuous. Then, we introduce the Gaussian
kernel, denoted by $G_{\varepsilon}(v)$, where the standard deviation is
$\varepsilon>0$. In this way, we have the following smoothed denoiser [34]

$\displaystyle\tilde{\mathbb{V}}(v)=S_{t}(v)+D_{t}(v)*G_{\varepsilon}(v),$
(30)

where $*$ denotes the convolution operator. By using the smoothed denoiser, we
can obtain the accurate state evolution of the HT-AMP. We refer the reader to
[34] for more details of the smoothed denoiser. Alternatively, one can smooth
the denoiser according to Monte Carlo sampling [30], which is similar to the
method introduced in the next section, hence the details are omitted here for
simplicity.

Note that the frequencies in the channel signal $\bm{h}_{n}$ are in continuous
space. Hence the on-grid CS approach HT-AMP used here yields performance loss.
In the nest section, we introduce the off-grid CS approach based on atomic
norm and super resolution AMP to solve this problem.

##  IV Proposed Off-grid Solutions based on AMP & Atomic Norm

In this section, we introduce AN techniques to model the sparse representation
to achieve the super-resolution performance. Then we propose super resolution
AMP based on greedy denoiser by exploiting the sparsity in mmWave channel.

###  IV-A Towards Super Resolution Solver Based on Atomic Norm

We investigate better denoiser for the channel estimation. Notice that for the
channel of mmWave communication, the number of paths is usually much smaller
than the number of BS antennas, i.e., $L_{n}\ll M$ for $n=1,2,...,N$. Hence,
there are also sparse nature in $\bm{h}_{n}$. However, the channel signal
$\bm{h}_{n}$ is a linear combination of complex sinusoids, the phase of which
is arbitrary, i.e. in continuous space. Therefore, using $\ell_{1}$-norm to
construct a sparse representation must place the frequencies onto discrete
grids. Obviously, this step will introduce errors. To get rid of the grid
operation, we use the atomic norm to achieve off-grid property.

Define ${\bm{a}}(f)$ as an atom and ${\mathcal{A}}=\\{\bm{a}(f):f\in[0,1)\\}$
as the atom set.
$\bm{x}_{n}=\alpha_{n}\sum_{\ell=1}^{L_{n}}c_{n,\ell}\bm{a}(f_{n,\ell})$ in
(II-A) can be regarded as
$\bm{x}_{n}=\sum_{\ell=1}^{L_{n}}c_{n,\ell}\bm{a}(f_{n,\ell})$ with
$c_{n,\ell}=\sqrt{\rho}\bar{c}_{n,\ell}$ when user is active and
$c_{n,\ell}=0$ otherwise. Then, the $\ell_{p}$-atomic norm related to (II-A)
is defined as follows.

######  Definition 1.

The $\ell_{p}$-atomic norm for $\bm{x}\in\mathbb{C}^{M\times 1}$ is

$\displaystyle\|\bm{x}\|_{{\mathcal{A}},p}=\mathop{\inf}\limits_{f_{\ell}\in[0,1)\hfill\atop
c_{\ell}\in\mathbb{C}\hfill}\left\\{\left.\left({\sum\limits_{\ell}{|c_{\ell}|^{p}}}\right)^{1/p}\right|\bm{x}=\sum\limits_{\ell}{c_{\ell}{\bm{a}}(f_{\ell})}\right\\}.$
(31)

The $\ell_{p}$-atomic norm can enforce sparsity in the atom set $\cal A$.
Then, we can formulate the following optimization problem according to (11)

$\displaystyle\\{\bm{\hat{x}}_{n}\\}_{n=1}^{N}=$
$\displaystyle~{}\mathop{\arg\min}\limits_{\\{\bm{x}_{n}\\}_{n=1}^{N}}\sum_{n=1}^{N}\|\bm{x}_{n}\|_{{\cal
A},p},$ (32) s.t.
$\displaystyle~{}{\left\|\bm{Y}-\sum_{n=1}^{N}\bm{u}_{n}\bm{x}_{n}^{T}\right\|_{F}}\leq\zeta,$

where $\zeta>0$ is a weight factor. However, (32) cannot be directly solved
because the $\ell_{p}$-atomic norm is essentially a semi-infinite program. In
the next section, we propose the super resolution AMP for solving (32). Before
that, we introduce a way to solve this problem when $p=1$ so that we can
better understand it.

In the case of $p=1$, we can use the equivalent form of the atomic norm for
the atom set $\cal A$ [19] to solve (32). Specifically, we have

$\displaystyle\|\bm{x}\|_{{\cal
A},1}=\mathop{\inf}\limits_{\bm{v}\in\mathbb{C}^{M\times
1},t\in\mathbb{R}}\left\\{\begin{array}[]{l}\frac{1}{2M}{\rm{Tr}}({\rm
Toep}(\bm{v}))+\frac{t}{2},\\\ {\rm s.t.}\left[{\begin{array}[]{*{20}{c}}{{\rm
Toep}(\bm{v})}&\bm{x}\\\ {{\bm{x}^{H}}}&t\end{array}}\right]\succeq
0\end{array}\right\\},$ (37)

where $\text{Toep}(\cdot)$ denotes the Hermitian Toeplitz matrix whose first
column is the input vector, i.e.,

$\displaystyle\text{Toep}(\bm{v})=\left[{\begin{array}[]{*{20}{c}}{v_{1}}&{{v_{2}^{*}}}&{\cdots}&{{v_{M}^{*}}}\\\
{v_{2}}&{v_{1}}&{\cdots}&{v_{M-1}^{*}}\\\
{\vdots}&{\vdots}&{\ddots}&{\vdots}\\\
{v_{M}}&{v_{M-1}}&{\cdots}&{v_{1}}\end{array}}\right].$ (42)

The relationship of (37) and (42) is given by

$\displaystyle{\rm Toep}(\bm{v})=$
$\displaystyle~{}\sum_{\ell}|c_{\ell}|\bm{a}(f_{\ell})\bm{a}(f_{\ell})^{H},$
(43) $\displaystyle t=$ $\displaystyle~{}\sum_{\ell}|c_{\ell}|.$ (44)

Then, we can apply (37) and transform (32) to the following semidefinite
program (SDP):

$\displaystyle\\{\bm{\hat{x}}_{n}\\}_{n=1}^{N}=$
$\displaystyle~{}\mathop{\arg\min}\limits_{\\{\bm{x}_{n},\bm{v}_{n},t_{n}\\}_{n=1}^{N}}\frac{\zeta}{2M}\sum_{n=1}^{N}{\text{Tr}}\left({\rm
Toep}(\bm{v}_{n})\right)+\sum_{n=1}^{N}\frac{\zeta t_{n}}{2},$ s.t.
$\displaystyle~{}\left\|\bm{Y}-\sum_{n=1}^{N}\bm{u}_{n}\bm{x}_{n}^{T}\right\|_{F}\leq\zeta,$
(45)
$\displaystyle~{}\left[{\begin{array}[]{*{20}{c}}{\text{Toep}(\bm{v}_{n})}&\bm{x}_{n}\\\
{{\bm{x}_{n}^{H}}}&t_{n}\end{array}}\right]\succeq 0,~{}n=1,2,...,N.$ (48)

The above SDP problem can be solved by CVX [39], SeDuMi [40], SDPT3 [41] and
any other convex solvers. However, these solver tend to be slow as the size of
the problem becomes large. In addition, in the case of $p=0$, the problem (32)
is non-convex, so cannot be solved with a convex solver. Hence, we next
propose the super-resolution AMP with faster speed.

###  IV-B Super-Resolution AMP Based on Greedy Denoiser

####  IV-B1 Greedy denoiser and its smooth form

We combine the atomic norm optimization with the framework of AMP to solve
(32), i.e., the AMP algorithm uses the atomic norm-based denoiser in each
iteration. We call such AMP algorithm that employs a super-resolution denoiser
as S-AMP. For clarity, we come back to the 1D-array case and denote the super-
resolution denoiser as $\eta_{t}(\cdot)$ which can be applied to
$\bm{\tilde{x}}=\bm{x}+\bm{\Sigma}_{t}^{1/2}\bm{v}$ with $\bm{v}\sim{\cal
CN}(\bm{0},\bm{I}_{M})$ and will return an estimate of $\bm{x}$ that is
hopefully closer to $\bm{x}$ than $\bm{x}+\bm{\Sigma}_{t}^{1/2}\bm{v}$. S-AMP
assumes that $\bm{x}_{n}$ is sparse in continuous frequency domain. Suppose we
use the denoiser

$\displaystyle\eta_{t}(\bm{\tilde{x}})=$
$\displaystyle~{}\arg\min_{\bm{x}\in{\mathbb{C}}^{M\times
1}}\|\bm{x}\|_{\mathcal{A},p},$ (49) s.t.
$\displaystyle~{}\|\bm{\tilde{x}}-\bm{x}\|_{2}\leq\gamma,$

where $\gamma$ is a threshold that can be adaptively changed in each
iteration. The algorithm can obtain the estimate of $\bm{X}$ via iterating
(14) and (15). Note that (49) is non-convex and its global optimum may not be
obtained in polynomial time. Specifically, inspired by recent work on atomic
norm minimization based on the conditional-gradient method [42, 43], we adopt
the following procedure as the denoiser in the case of $p=0$, which includes
the initialization, selection, local-optimization, least-squares and residual
update steps. The iterative process is below and the detail is in VI-C Super-
resolution AMP (S-AMP)

Algorithm 1 Super-resolution AMP (S-AMP)

0: Observation vector $\bm{\tilde{x}}$, Maximum iterations $I$, $I^{\prime}$,
Step size $\mu_{i}$, Stopping criteria $\epsilon$, $\delta$, Thresholds
$\gamma$, $\tilde{\delta}$

0: Estimated spectral lines $\mathcal{T}$

1: Initialization:

2: $\mathcal{T}\leftarrow\emptyset$

3: $\bm{r}_{\text{res}}\leftarrow\bm{\tilde{x}}$

4: repeat

5: Selection:

6: Compute highest correlated atom $f^{\diamond}$

7: Update $\mathcal{T}:\mathcal{T}\leftarrow\mathcal{T}\cup\\{f^{\diamond}\\}$

8: Local optimization:

9: Define $\bm{c}$ and $\bm{f}$ using $\mathcal{T}$

10: Refine $\bm{f}$ by solving the optimization problem

11: Calculate gradient and Hessian matrix

12: Update $\bm{f}$ using Newton’s method:
$\bm{f}^{i+1}=\bm{f}^{i}-\mu_{i}\bm{K}(\bm{f}^{i})^{-1}\bm{p}(\bm{f}^{i})$

13: Least-squares:

14: Estimate $\bm{c}$ by solving the least-squares problem

15: Remove any atoms in $\mathcal{T}$ with coefficients smaller than
$\tilde{\delta}$

16: Residual update:

17: $\bm{r}_{\text{res}}=\bm{\tilde{x}}-\bm{\Phi}(\bm{\hat{f}})\bm{\hat{c}}$

18: until $\|\bm{r}_{\text{res}}\|_{2}^{2}\leq\epsilon$ or maximum iteration
number $I^{\prime}$ is reached

The greedy denoiser in steps 1) to 5) can be directly used in AMP iterations
(14) and (15). However, as mentioned above, we also need to smooth the greedy
denoiser, i.e., find a smooth version that without discontinuities and behave
almost the same as the original noise reducer, so as to satisfy the state
evolution equations. Suppose $\eta_{t}(\cdot)$ is discontinuous. The smoothed
denoiser can be obtained via [29]

$\displaystyle\tilde{\eta}_{t}(\bm{\tilde{x}})=\int_{{\bm{\zeta}}\in\mathbb{C}^{M\times
1}}\eta_{t}(\bm{\tilde{x}}-\bm{\zeta})\frac{1}{r^{M}(2\pi)^{M/2}}e^{-\|{\bm{\zeta}}\|_{2}^{2}/2r^{2}}d{\bm{\zeta}},$
(50)

where $d{\bm{\zeta}}=d\zeta_{1}d\zeta_{2}...d\zeta_{M}$.

Since (50) integrates over $\mathbb{C}^{M\times 1}$ based on the unknown form
$\eta_{t}(\cdot)$, we calculate this via Monte Carlo sampling [30]. First we
randomly generated some complex Gaussian vectors
$\bm{b}^{1},\bm{b}^{2},...,\bm{b}^{J_{1}}$ with standard deviation $r$. Then,
we smooth the denoiser by

$\displaystyle\hat{\tilde{\eta}}_{t}(\bm{\tilde{x}})=\frac{1}{J_{1}}\sum_{j=1}^{J_{1}}\eta_{t}(\bm{\tilde{x}}+\bm{b}^{j}),$
(51)

where $\bm{b}^{j}\sim{\cal CN}(\bm{0},r^{2}\bm{I}_{M})$ for all $j$. As a
result, $\bm{x}_{n}^{t+1}$ can be updated via (14) with $\eta_{t}(\cdot)$
replaced by the smooth greedy denoiser $\hat{\tilde{\eta}}_{t}(\cdot)$. In
addition, the SE of S-AMP can be obtained via (III-A2)-(23) with the eta
function in (III-A2) replaced by (51).

####  IV-B2 Derivative

To perform the AMP iterations in (14) and (15),
$\sum_{n=1}^{N}\frac{\partial\hat{\tilde{\eta}}_{t}}{\partial\bm{\tilde{x}}}\left(\bm{\tilde{x}}_{n}^{t}\right)$
should also be calculated. However, this is difficult since we do not have
access to the explicit form of the greedy denoiser $\eta_{t}(\cdot)$. To solve
this problem, we also calculate
$\sum_{n=1}^{N}\frac{\partial\hat{\tilde{\eta}}_{t}}{\partial\bm{\tilde{x}}}\left(\bm{\tilde{x}}_{n}^{t}\right)$
by using the Monte Carlo simulation [30]. That is, for $1\leq m\leq M$ and
$1\leq n\leq N$, we first generate $J_{2}$ vectors
$\bm{d}_{m,n}^{1},\bm{d}_{m,n}^{2},...,\bm{d}_{m,n}^{J_{2}}$ whose $m$-th
element obeys i.i.d. complex Gaussian distribution ${\cal CN}(0,1)$ and the
rest of the elements are zero.

Next, define

$\displaystyle\bm{\upsilon}_{m}\triangleq\sum_{n=1}^{N}\frac{\partial\hat{\tilde{\eta}}_{t}(\bm{\tilde{x}_{n}})}{\partial\tilde{x}_{n}(m)}\in{\mathbb{C}}^{M\times
1}.$ (52)

Then, from the following Theorem 1, (the proof process is in Section VI-B), we
can find that for large $J_{1}$, $J_{2}$ and for very small $\varepsilon$,
$\bm{\upsilon}_{m}$ can be approximately obtained by averaging

$\displaystyle\bm{\upsilon}_{m}\approx\frac{1}{J_{2}}\sum_{k=1}^{J_{2}}\bm{\upsilon}_{m,k},$
(53)

where

$\displaystyle\bm{\upsilon}_{m,k}=$
$\displaystyle~{}\sum_{n=1}^{N}\frac{d_{m,n}^{k}(m)}{\varepsilon}\left(\hat{\tilde{\eta}}_{t}(\bm{\tilde{x}}_{n}+\varepsilon\bm{d}_{m,n}^{k})-\hat{\tilde{\eta}}_{t}(\bm{\tilde{x}}_{n})\right)$
$\displaystyle=$
$\displaystyle~{}\sum_{n=1}^{N}\frac{d_{m,n}^{k}(m)}{\varepsilon
J_{1}}\sum_{j=1}^{J_{1}}\left(\eta_{t}(\bm{\tilde{x}}_{n}+{\bm{b}}^{j}+\varepsilon\bm{d}_{m}^{k})-\eta_{t}(\bm{\tilde{x}}_{n}+{\bm{b}}^{j})\right),$
(54)

where $d_{m,n}^{k}(m)$ is the $m$-th element of $\bm{d}_{m,n}^{k}$.

###  IV-C Computational Complexity Analysis

(a)  Figure 1: Convergence behaviors of the proposed detectors and their state
evolutions.

In this section, we analyze the computational complexity of the proposed
methods in the Table I. In addition to the denoiser, the complexity of the AMP
iterations in (14) and (15) mainly comes from $\bm{U}\bm{X}^{t+1}$, which is
$\mathcal{O}(QNM)$ per iteration. For the denoisers, the computational
complexity of the group soft-thresholding denoiser in (III-A1) is
$\mathcal{O}(M)$, which is smaller than that of the MMSE denoiser
($\mathcal{O}(M^{3})$) in [1]. In addition, the computational complexities of
the hard-thresholding denoiser in (26) and (27) are respectively
$\mathcal{O}(M\tilde{M})$ and $\mathcal{O}(M^{2}\tilde{M})$. As for the SDP
problem, if it is solved by the interior point method, the complexity is
${\cal O}(M^{6}N)$ per iteration [44, 24]. The main complexity of the greedy
denoiser comes from (64) and (65), with respectively complexities ${\cal
O}(M^{3})$ and ${\cal O}(M^{2}|\mathcal{T}|)$ per inner iteration (we named
the iterations of steps 1)-5) in the greedy denoiser as the inner iterations
to distinguish them from the AMP iterations). Since $M<N$ and $M<Q$, the
computational complexities of GST-AMP, HT-AMP, SDP and S-AMP methods are
$\mathcal{O}(QNM)$, $\mathcal{O}(QNM)$, ${\cal O}(M^{6}N)$ and
$\mathcal{O}(QNM)$, respectively.

TABLE I: Comparison of Computational Complexities Method | Complexity
---|---
AMP iterations | $\mathcal{O}(QNM)$
Group soft-thresholding denoiser | $\mathcal{O}(M)$
MMSE denoiser [1] | $\mathcal{O}(M^{3})$
Hard-thresholding denoiser | $\mathcal{O}(M\tilde{M})$
SDP (interior point method) [44, 24] | $\mathcal{O}(M^{6}N)$
Greedy denoiser (gradient) | $\mathcal{O}(M^{3})$
Greedy denoiser (Hessian) | $\mathcal{O}(M^{2}|\mathcal{T}|)$
GST-AMP | $\mathcal{O}(QNM)$
HT-AMP | $\mathcal{O}(QNM)$
S-AMP | $\mathcal{O}(QNM)$

##  V Simulations

###  V-A Simulation Setup

In this section, we present the numerical simulation examples to demonstrate
the performance of the proposed algorithms. In addition to the parameters
specified in the following simulations, the basic parameters are set as
follows, the number of devices in the cell is set as $N=2000$. We set the
probability $\epsilon=0.05$ and hence the active user number (pilot sequence)
is $K=100$. A quadrature phase-shift keying (QPSK) modulation is used for
pilot sequence[25]. Note that in the following simulations, we assume that
each user has a unique pilot sequence to be consistent with the method in [1,
17] for a fair comparison. The complex multi-path gains $c_{n,\ell}$ in (1)
are set as $c_{n,\ell}\sim{\cal CN}(0,1)$, and the antenna number $M=32$. The
communication system uses a signal with $Q=1000$ data symbols, and a total
bandwidth of 100 MHz, i.e., the frequency spacing between adjacent subcarriers
is 100 kHz. The transmit power of the active users is $\rho^{\rm
pilot}=30~{}\text{dBm}$. The power spectral density of the AWGN at the BS is
set as $-174~{}\text{dBm/Hz}$ with a -94dB pathloss.

To evaluate the performance achieved by the proposed detectors, we follow the
same way in [1] to calculate the missed detection and false alarm to evaluate
the proposed detectors. We invite readers to refer to [1] for more details.
Besides, we calculate the normalized mean-squared-error (NMSE)
$\|\bm{\hat{X}}-\bm{X}\|_{F}/\|\bm{X}\|_{F}$ to evaluate the proposed
detectors. In the following simulations, we compare the performance of the
proposed GST-AMP, HT-AMP and S-AMP detectors with the MMSE-AMP detector in [1]
and MMSE-NET detector in [17]. However, the proposed SDP detector is not
included in the comparison because its computational complexity is 5 orders of
magnitude higher than other proposed detectors, which is analyzed at the end
of Section IV. Moreover, we calculate the state evolution of the proposed GST-
AMP, HT-AMP and S-AMP detectors, which are respectively named GST-AMP-SE, HT-
AMP-SE and S-AMP-SE. All performance curves of GST-AMP, HT-AMP, MMSE-AMP and
MMSE-NET detectors are averaged results after at least 20000 runs, while those
of S-AMP detector are averaged results after at least 2000 runs.

###  V-B Simulation Results

####  V-B1 Convergence behavior

The convergence behavior is analyzed in Fig. 1. The NMSE in each iteration is
calculated, we can see that HT-AMP and S-AMP detectors can converge within 10
iterations, while GST-AMP detector can converge within 15 iterations.
Therefore, in subsequent simulations, we let the detectors iterate until the
difference between the two updates is less than $\epsilon=10^{-6}$, or iterate
until reach the maximum number of iterations, which is set to 10 for HT-AMP
and S-AMP detectors, and 15 for GST-AMP detector. In addition, we can see that
the state evolution curves GST-AMP-SE, HT-AMP-SE and S-AMP-SE coincide with
the NMSE curves of GST-AMP, HT-AMP and S-AMP detectors, respectively.

(a)  (b)  (c)

Figure 2: Numerical simulation results against the number of antennas ($M$) at
the BS. (a) False alarm; (b) Missed detection; (c) NMSE and state evolution.

####  V-B2 Simulation results under different $M$

Fig. 2(a) and (b) respectively show the false alarm and missed detection
versus $M$, i.e., the number of antennas at the BS. Fig. 2(c) shows the NMSE
of $\bm{X}$ in (9) versus $M$. We can see that when $M$ increases, the
performance of the proposed HT-AMP and S-AMP detectors will be significantly
better. This is because when $M$ increases, the sparsity of the millimeter
wave channel becomes better. In contrast, the MMSE-AMP detector in [1] and the
proposed GST-AMP detector do not exploit the sparsity of the mmWave channel,
so the performance does not change much, as the number of antennas $M$
increases, the number of groups also increases, and each group becomes smaller
in size. This results in a decrease in the sparsity of each group, which in
turn affects the performance of the GST-AMP. And the GST-AMP performance is
worse than that of HT-AMP and S-AMP detectors for large $M$.

Moreover, compared to S-AMP detector, the on-grid based HT-AMP detector
suffers performance loss, since the frequencies do not fall onto discrete
grids and S-AMP detector can obtain more accurate channel frequencies. Note
that although the NMSE of the MMSE-AMP detector is relatively small, its false
alarm and missed detection rate are higher than that of the HT-AMP detector
for large $M$. This is because the HT-AMP detector is less likely to fail, but
when it fails, it may produce a larger error NMSE, which will result in a
higher average value. Furthermore, we can see that when $M$ is large enough,
the state evolution curves GST-AMP-SE, HT-AMP-SE and S-AMP-SE coincide with
the NMSE curves of GST-AMP, HT-AMP and S-AMP detectors, respectively, while
there is a gap between them when $M$ is small since the detectors may have
error under non-sparse conditions.

(a)  (b)  (c)

Figure 3: Numerical simulation results as functions of the number of active
users ($K$). (a) False alarm; (b) Missed detection; (c) NMSE and SE.

####  V-B3 Simulation results under different $K$

Figs. 3(a) and (b) show the false alarm and missed detection versus the
numbers of active users $K$, respectively. Fig. 3(c) shows the NMSE of
$\bm{X}$ versus $K$. We can see that the proposed HT-AMP and S-AMP detectors
can achieve better performance than the MMSE-AMP detector and the proposed
GST-AMP detector, since they further exploit the sparsity of the mmWave
channel. In addition, we can see that when $K$ becomes larger, the sparsity in
the user activity pattern will be worse. Hence, all detectors have performance
losses in terms of NMSE, false alarm and missed detection rates, except for
GST-AMP detector where the false alarm rate becomes better. This may be due to
the lack of stability of the group soft-thresholding denoiser, i.e., there
will be some erroneous estimates which lead to deviations in the average
value, but when $K$ becomes larger, the wrong estimate is more likely to be
included in the correct result. Another possible reason is due to an
insufficient number of samples used in our experiments. Consequently, the
initial random seeds used for generating values may have a relatively
significant impact on the false alarm rate.

(a)  (b)  (c)

Figure 4: Numerical simulation results as functions of the length of data
symbols ($Q$). (a) False alarm; (b) Missed detection; (c) NMSE and state
evolution.

####  V-B4 Simulation results under different $Q$

Figs. 4(a) and (b) show the false alarm and missed detection versus the length
of data symbols $Q$, respectively. Fig. 4(c) shows the NMSE of $\bm{X}$ versus
$Q$. We can see that the proposed HT-AMP and S-AMP detectors can achieve
better performance than the MMSE-AMP detector and the proposed GST-AMP
detector. In addition, when $Q$ becomes larger, the performance of the
proposed detectors become better, since the sparsity condition becomes better.
Moreover, we can see that when $Q$ is large enough, the state evolution curves
GST-AMP-SE, HT-AMP-SE and S-AMP-SE coincide with the NMSE curves of GST-AMP,
HT-AMP and S-AMP detectors, respectively, while there is a gap between them
when $Q$ is small since the detectors may have error under non-sparse
conditions.

####  V-B5 Simulation results under different SNR

Figs. 5(a) and (b) show the false alarm and missed detection versus source
SNR, respectively. Fig. 5(c) shows the NMSE of $\bm{X}$ versus SNR. With the
SNR increases, all detectors have performance improvement in terms of NMSE,
false alarm and missed detection rates, except for GST-AMP and MMSE detectors
where the false alarm rate may increase. This phenomenon also appears in Fig.
3, which may be caused by a trade-off between the false alarm and missed
detection rates. When the missed detection rate drops a lot, it shows that we
are more inclined to regard noise as the target, and the false alarm rate
increases instead. Another possible reason is due to the limited number of
samples used in our experiments, the false alarm rate was impacted by false
signals at specific conditions, particularly in the SNR=20dB and 30dB cases.
These false signals likely contributed to the atypical trend.

(a)  (b)  (c)

Figure 5: Numerical simulation results as functions of the length of SNR. (a)
False alarm; (b) Missed detection; (c) NMSE and state evolution.

####  V-B6 Performance on Estimator

In order to analyze the performance of our AMP algorithm in user detection and
channel estimation compared with the state evolution result, we estimate
$\bm{X}$ from the received signal matrix $\bm{Y}$ in (11). In this setting, we
fix $X$ with antennas $M=32$ and change the noise power.

Fig. 6 shows the estimation results of each antenna value using weak noise
power (SNR = 50dB) and strong noise power (SNR = 0dB) for GST-AMP, GST-AMP-SE
algorithm, HT-AMP, HT-AMP-SE algorithm, and S-AMP, S-AMP-SE algorithm,
respectively. This result demonstrates that both numerical simulations can
estimate the true value with high precision when the noise power is low, and
vice versa.

(a)  (b)  (c)  (d)

(e)  (f)

Figure 6: Estimation of $X$ value with proposed methods. (a)-(c) SNR=50dB.
(d)-(e) SNR=0dB.

(a)  (b)  (c)

Figure 7: Statistics NMSE distribution between the $X$ and $\hat{X}$ with
SNR=30dB in (a) GST-AMP, (b) HT-AMP and (c) S-AMP, respectively.

To evaluate the performance of our estimator under normal settings, we
generate an estimate $\hat{X}$ from $X$ by the GST-AMP, HT-AMP and S-AMP
algorithm, respectively. We add different noises randomly 10000 times with SNR
= 30dB in $X$ and summary the statistics of the NMSE between $X$ and $\hat{X}$
in Fig. 7. It’s seen that most of the generated estimates follow the three-
sigma rule of thumb law [45], where the results fit the norm distribution with
$\mu=0.046$ and $\sigma=0.0055$ in GST-AMP, the norm distribution with
$\mu=0.022$ and $\sigma=0.0028$ in HT-AMP, and the half-norm distribution with
$\mu=0$ and $\sigma=0.0039$ in S-AMP.

####  V-B7 Running Time

We have measured the average running time for GST-AMP, HT-AMP, S-AMP, and
MMSE-AMP algorithms. We have considered the different numbers of channels
($M$) when $K=100$, and for varying levels of user activity ($K$) when $M=32$.

TABLE II: Running times of the proposed algorithms for the different channel (M) and user activity (K). Method | channels (K=100) | channels (M=32)
---|---|---
| M=8 | M=16 | M=32 | M=64 | K=50 | K=100 | K=200 | K=400 | K=1000
GST-AMP | 0.29 s | 0.48 s | 0.73 s | 1.88 s | 0.74 s | 0.73 s | 0.81 s | 0.842 s | 0.86 s
HT-AMP | 0.24 s | 0.29 s | 0.47 s | 1.36 s | 0.48 s | 0.47 s | 0.49 s | 0.54 s | 0.55 s
S-AMP | 0.74 s | 0.70s | 0.96 s | 2.14 s | 0.84 s | 0.96 s | 2.48 s | 4.46 s | 10.39 s
MMSE-AMP | 0.086 s | 0.096 s | 0.17 s | 0.41 s | 0.13 s | 0.17 s | 0.25 s | 0.44 s | 0.55 s

These results demonstrate the performance of the proposed algorithms in
practical scenarios and provide insights into their suitability for varying
channel and user activity conditions.

##  VI Conclusions

In this paper, we solved the problem of jointly user activity detection and
channel estimation for establishing successful communication between the
cellular base stations and the devices. We first proposed a group soft-
thresholding AMP detector to exploit only the sparsity in user activity. Then,
we proposed a hard-thresholding AMP detector based on the on-grid CS approach
by using the user and channel sparsities. Moreover, we proposed a direct SDP
solver and a super-resolution AMP detector to solve the off-grid problem,
where a greedy method based the atomic norm technique is used as the super-
resolution denoiser. We presented the state evolution results of the proposed
methods. A large number of experimental results validate the proposed method,
which is expected to achieve promising performance for user activity detection
and channel estimation in B5G or 6G wireless communication networks.

## Appendix

###  VI-A Table of Symbols

We have arranged and listed the symbols that appear in Table III.

TABLE III: Table of Symbols Symbol | Description
---|---
$N$ | Number of users
$M$ | Number of BS antennas
$K$ | Number of active devices
$L_{n}$ |  Number of paths of the $n$-th user
$\bar{c}_{n,\ell}$ |  Complex gain of the $\ell$-th path of the $n$-th user
$\bm{a}(f)$ | Steering vector of the receive array
$f$ | Angle of arrival
$d$ | Separation of adjacent antenna elements
$\lambda_{w}$ | Wavelength of the transmitted signal
$\epsilon$ | User activity probability
$\alpha_{n}$ | User activity indicator
$\bm{u}_{n}$ |  Data symbols of the $n$-th user
$Q$ | Length of pilot sequence
$\rho$ | Pilot transmit power
$\bm{h}_{n}$ |  Complex uplink channel vector of the $n$-th user
$\bm{X}$ | Matrix of channel vectors
$\bm{Y}$ | Matrix of received signal
$\bm{Z}$ | Noise matrix
$\bm{R}^{t}$ |  Residual matrix at iteration $t$
$\bm{x}_{n}^{t}$ |  Channel vector of the $n$-th user at iteration $t$
$\eta_{t}(\cdot)$ | Non-linear function known as denoiser
$\beta_{n}$ |  Path-loss and shadowing component of the $n$-th user
$\bm{\tilde{x}}_{n}^{t}$ |  Input argument of the denoiser at iteration $t$ for user $n$
$\gamma$ | Weighting parameter
$\bm{w}_{n}$ |  Weight vector of the $n$-th user
$\bm{\hat{x}}_{n}$ |  Estimate of the $n$-th user’s channel vector
$\bm{\hat{X}}$ | Estimate of the matrix of channel vectors
$\bm{\hat{Y}}$ | Estimate of the matrix of received signal
$\bm{\hat{Z}}$ | Estimate of the noise matrix
$\bm{\hat{R}}^{t}$ |  Estimate of the residual matrix at iteration $t$
$\bm{\hat{u}}_{n}$ |  Estimate of the data symbols of the $n$-th user

###  VI-B Theorem 1

######  Theorem 1.

Denote $\bm{d}_{m,1},\bm{d}_{m,2},...,\bm{d}_{m,N}\in{\mathbb{C}}^{M\times 1}$
whose $m$-th element obeys i.i.d. complex Gaussian distribution ${\cal
CN}(0,1)$ and the rest of the elements are zero. Then we have

$\displaystyle\bm{\upsilon}_{m}=\lim_{\varepsilon\to
0}{\mathbb{E}}\left\\{\frac{1}{\varepsilon}\sum_{n=1}^{N}d_{m,n}(m)\left({\eta_{t}}(\bm{\tilde{x}}_{n}+\varepsilon\bm{d}_{m,n})-{\eta_{t}}(\bm{\tilde{x}}_{n})\right)\right\\}$

provided that ${\eta_{t}}(\cdot)$ admits a well-defined second-order Taylor
expansion.

###### Proof:

We write the second-order Taylor expansion of
${\eta_{t}}(\bm{\tilde{x}}+\varepsilon\bm{d}_{m,n})$ as

$\displaystyle{\eta_{t}}(\bm{\tilde{x}}_{n}+\varepsilon\bm{d}_{m,n})={\eta_{t}}(\bm{\tilde{x}}_{n})$
$\displaystyle\quad+\varepsilon\frac{\partial{\eta}_{t}(\bm{\tilde{x}}_{n})}{\partial\tilde{x}_{n}(m)}d_{m,n}(m)+\varepsilon^{2}{\bm{r}}_{\eta}(\bm{\tilde{x}}_{n},d_{m,n}(m)).$
(55)

where ${\bm{r}}_{\eta}\in\mathbb{C}^{M\times 1}$ is the Lagrange remainder
vector, whose each element corresponding to each component of
$\eta_{t}(\cdot)$ and we hence have
${\mathbb{E}}\left|r_{\eta}(m)\right|<+\infty$ for $m=1,2,...,M$.

Then, subtracting ${\eta_{t}}(\bm{\tilde{x}}_{n})$ from (VI-B) and we obtain

$\displaystyle{\mathbb{E}}\left\\{\sum_{n=1}^{N}d_{m,n}(n)\left[{\eta_{t}}(\bm{\tilde{x}}_{n}+\varepsilon\bm{d}_{m,n})-{\eta_{t}}(\bm{\tilde{x}}_{n})\right]\right\\}$
$\displaystyle=$
$\displaystyle~{}\varepsilon{\mathbb{E}}\left\\{\sum_{n=1}^{N}d_{m,n}^{2}(n)\frac{\partial{\eta}_{t}(\bm{\tilde{x}}_{n})}{\partial\tilde{x}_{n}(m)}\right\\}$
$\displaystyle~{}+\varepsilon^{2}{\mathbb{E}}\left\\{\sum_{n=1}^{N}d_{m,n}(n){\bm{r}}_{\eta}(\bm{\tilde{x}}_{n},d_{m,n}(m))\right\\}$
$\displaystyle=$
$\displaystyle~{}\varepsilon\sum_{n=1}^{N}\frac{\partial{\eta}_{t}(\bm{\tilde{x}}_{n})}{\partial\tilde{x}_{n}(m)}+C_{2}\varepsilon^{2},$
(56)

where
${\mathbb{E}}\left\\{\sum_{n=1}^{N}d_{m,n}(n){\bm{r}}_{\eta}(\bm{\tilde{x}}_{n},d_{m,n}(m))\right\\}=C_{2}$
and $|C_{2}|<+\infty$ because
${\mathbb{E}}\left\\{{\bm{r}}_{\eta}(\bm{\tilde{x}}_{n},d_{m,n}(m))\right\\}<+\infty$
for $k=1,2,...,N$ and $d_{m,n}(n)$ has bounded higher-order moments since
$d_{m,n}(n)\sim{\cal CN}(0,1)$. When $\varepsilon\to 0$, we immediately see
that

$\displaystyle~{}\lim_{\varepsilon\to
0}\frac{1}{\varepsilon}{\mathbb{E}}\left\\{\sum_{n=1}^{N}d_{m,n}(n)\left[{\eta_{t}}(\bm{\tilde{x}}_{n}+\varepsilon\bm{d}_{m,n})-{\eta_{t}}(\bm{\tilde{x}}_{n})\right]\right\\}$
$\displaystyle=$
$\displaystyle~{}\sum_{n=1}^{N}\frac{\partial{\eta}_{t}(\bm{\tilde{x}}_{n})}{\partial\tilde{x}_{n}(m)}=\bm{\upsilon}_{m}.$
(57)

Then we complete the proof. ∎

###  VI-C Super-resolution AMP (S-AMP)

  1. 1.

Initialization: Define $\cal{T}$ as the set of estimated spectral lines, and
denote $\bm{r}_{\text{res}}\in\mathbb{C}^{M\times 1}$ as a residual. Then,
initialize $\cal{T}\leftarrow\emptyset$ and
$\bm{r}_{\text{res}}\leftarrow\bm{\tilde{x}}$.

  2. 2.

Selection: Compute the atom in $\cal{A}$ that has the highest correlation with
the current residual $\bm{r}_{\text{res}}$ and update $\cal{T}$ i.e, compute

$\displaystyle f^{\diamond}=\arg\max_{f\in[0,1)}\emph{corr}(f),$ (58)

where
$\emph{corr}(f):=\left|\langle{\bm{a}(f),\bm{r}_{\text{res}}}\rangle\right|$
and let ${\mathcal{T}}\leftarrow\mathcal{T}\cup\\{f^{\diamond}\\}$. For atoms,
we calculate the highest correlation by approximately finding the position of
the maximum of $\emph{corr}(f)$ on the fine grid $f_{\text{grid}}$. To be more
efficient, we can construct the oversampled fast Fourier transform on
$\bm{r}_{\text{res}}$ and directly find the largest element of
$\emph{corr}(f)$ for $f\in f_{\text{grid}}$.

  3. 3.

Local optimization: After $\cal{T}$ is updated, define
$\bm{c}=\left[{c_{1}},{c_{2}},...,{c_{|\cal{T}|}}\right]^{T}\in\mathbb{C}^{|\cal{T}|}$
and
$\bm{f}=\left[{f_{1}},{f_{2}},...,{f_{|\cal{T}|}}\right]^{T}\in[0,1)^{|\cal{T}|}$
in $\bm{x}=\bm{\Phi}(\bm{f})\bm{c}$, where

$\displaystyle\bm{\Phi}(\bm{f})=[\bm{a}(f_{1}),\bm{a}(f_{2}),...,\bm{a}(f_{|\cal{T}|})]\in\mathbb{C}^{M\times|\cal{T}|}.$
(59)

Then, we refine $\bm{f}$ in $\cal{T}$ by solving

$\displaystyle\min_{\bm{f}\in[0,1)^{|\cal{T}|},\bm{c}\in\mathbb{C}^{|\cal{T}|}}{\left\|\bm{\tilde{x}}-\bm{\Phi}(\bm{f})\bm{c}\right\|_{2}^{2}}.$
(60)

Substituting the solution
$\bm{\hat{c}}=\bm{\Phi}(\bm{f})^{\dagger}\bm{\tilde{x}}$ back to (60) yields
[46, 47]:

$\displaystyle\bm{\hat{f}}=\arg\min_{\bm{f}\in[0,1)^{|\cal{T}|}}\text{Tr}\\{\bm{P}^{\perp}(\bm{f})\bm{R}\\},$
(61)

where

$\displaystyle\bm{R}=$
$\displaystyle~{}\bm{\tilde{x}}\bm{\tilde{x}}^{H}\in\mathbb{C}^{M\times M},$
(62) $\displaystyle\bm{P}^{\perp}(\bm{f})=$
$\displaystyle~{}\bm{I}_{M}-\bm{\Phi}(\bm{f})\bm{\Phi}(\bm{f})^{\dagger}\in\mathbb{C}^{M\times
M},$ (63)

and $(\cdot)^{\dagger}$ denotes the pseudo-inverse, i.e.,
$\bm{X}^{\dagger}=(\bm{X}^{H}\bm{X})^{-1}\bm{X}^{H}$. Next, we calculate the
gradient and Hessian matrix as follows [46, 47]

$\displaystyle{\bm{p}(\bm{f})}=$
$\displaystyle~{}\nabla_{f}\left[\text{Tr}\\{\bm{P}^{\perp}(\bm{f})\bm{R}\\}\right]\in\mathbb{R}^{{|\cal{T}|}\times
1}$ (64) $\displaystyle=$ $\displaystyle~{}-2{\Re}\left\\{{\text{vec-
diag}}\left[\bm{\Phi}^{\dagger}(\bm{f})\bm{R}\bm{P}^{\perp}(\bm{f})\bm{T}(\bm{f})\right]\right\\},$
$\displaystyle\bm{K}(\bm{f})=$
$\displaystyle~{}\nabla_{f}^{2}\left[\text{Tr}\\{\bm{P}^{\perp}(\bm{f})\bm{R}\\}\right]\in\mathbb{R}^{{|\cal{T}|}\times{|\cal{T}|}}$
(65) $\displaystyle\approx 2{\Re}$
$\displaystyle\left\\{(\bm{T}(\bm{f})^{H}\bm{P}^{\perp}(\bm{f})\bm{T}(\bm{f}))\odot({\bm{\Phi}(\bm{f})^{\dagger}\bm{R}\bm{\Phi}(\bm{f})^{\dagger
H}})^{T}\right\\},$

where ${\text{vec-diag}}[\cdot]$ denotes the operator that outputs a vector
generated by the diagonal elements of an input square matrix, and
$\bm{T}(\bm{f})$ is given by (3)

$\displaystyle\bm{T}(\bm{f})=$
$\displaystyle\left[\frac{\partial\bm{a}(f)}{\partial
f}{\bigg{|}}_{f=f_{1}},\frac{\partial\bm{a}(f)}{\partial
f}{\bigg{|}}_{f=f_{2}},...,\frac{\partial\bm{a}(f)}{\partial
f}{\bigg{|}}_{f=f_{|\cal{T}|}}\right]\in\mathbb{C}^{M\times{|\cal{T}|}}$
$\displaystyle=$ $\displaystyle\left[{\begin{array}[]{*{20}{c}}1&1&\cdots&1\\\
i2\pi e^{i2\pi f_{1}}&i2\pi e^{i2\pi f_{2}}&\cdots&i2\pi e^{i2\pi
f_{|\cal{T}|}}\\\ \vdots&\vdots&\ddots&\vdots\\\
i2\pi(M-1)e^{i2\pi(M-1)f_{1}}&i2\pi(M-1)e^{i2\pi(M-1)f_{2}}&\cdots&i2\pi(M-1)e^{i2\pi(M-1)f_{|\cal{T}|}}\end{array}}\right].$
(70)

Finally, we solve (61) based on the gradient and Hessian following the
Newton’s method as

$\displaystyle\bm{f}^{i+1}=\bm{f}^{i}-\mu_{i}\bm{K}(\bm{f}^{i})^{-1}{\bm{p}(\bm{f}^{i})},~{}~{}i=0,1,...$
(71)

where $\mu_{i}$ is a step size, and the iteration goes until the maximum
iteration number $I$ is reached or the condition
$\|\bm{K}(\bm{f}^{i})^{-1}\bm{p}(\bm{f}^{i})\|_{2}<\delta$ is satisfied. The
elements in $\cal{T}$ are used as the initialization point $\bm{f}^{0}$.

  4. 4.

Least-squares: Once $\mathcal{T}$ is obtained, we can estimate $\bm{c}$ by
solving the following least-squares problem:

$\displaystyle\bm{\hat{c}}=\arg\min_{\bm{c}\in\mathbb{C}^{{|\cal{T}|}}}{\left\|\bm{\tilde{x}}-\bm{\Phi}(\bm{\hat{f}})\bm{c}\right\|_{2}^{2}}.$
(72)

Based on the estimated coefficients $\bm{\hat{c}}$, we remove any atoms in
$\mathcal{T}$ whose corresponding coefficients are smaller than a predefined
threshold $\tilde{\delta}$.

  5. 5.

Residual update:

$\displaystyle\bm{r}_{\text{res}}=\bm{\tilde{x}}-\bm{\Phi}(\bm{\hat{f}})\bm{\hat{c}}$
(73)

and repeat steps 2) to 5) until $\|\bm{r}_{\text{res}}\|_{2}^{2}\leq\epsilon$,
or the maximum iteration number $I^{\prime}$ is reached.

##  VII Acknowledgement

This work was supported by the National Key R&D Program of China (Grant No.
2018YFE0202101, 2018YFE0202103). Le Zheng is the corresponding author.

## References

  * [1] L. Liu and W. Yu, “Massive connectivity with massive MIMO-Part I: Device activity detection and channel estimation,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 11, pp. 2933–2946, 2018.
  * [2] ——, “Massive connectivity with massive MIMO-Part II: Achievable rate characterization,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 11, pp. 2947–2959, 2018.
  * [3] ——, “Massive device connectivity with massive MIMO,” in _2017 IEEE International Symposium on Information Theory (ISIT)_. IEEE, 2017, pp. 1072–1076.
  * [4] M. Hasan, E. Hossain, and D. Niyato, “Random access for machine-to-machine communication in LTE-advanced networks: Issues and approaches,” _IEEE communications Magazine_ , vol. 51, no. 6, pp. 86–93, 2013.
  * [5] E. Björnson, E. De Carvalho, J. H. Sørensen, E. G. Larsson, and P. Popovski, “A random access protocol for pilot allocation in crowded massive MIMO systems,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 4, pp. 2220–2234, 2017.
  * [6] H. Zhu and G. B. Giannakis, “Exploiting sparse user activity in multiuser detection,” _IEEE Transactions on Communications_ , vol. 59, no. 2, pp. 454–465, 2010.
  * [7] H. F. Schepker and A. Dekorsy, “Compressive sensing multi-user detection with block-wise orthogonal least squares,” in _2012 IEEE 75th Vehicular Technology Conference (VTC Spring)_. IEEE, 2012, pp. 1–5.
  * [8] X. Xu, X. Rao, and V. K. Lau, “Active user detection and channel estimation in uplink CRAN systems,” in _2015 IEEE International Conference on Communications (ICC)_. IEEE, 2015, pp. 2727–2732.
  * [9] G. Wunder, P. Jung, and M. Ramadan, “Compressive random access using a common overloaded control channel,” in _2015 IEEE Globecom Workshops (GC Wkshps)_. IEEE, 2015, pp. 1–6.
  * [10] H. F. Schepker, C. Bockelmann, and A. Dekorsy, “Exploiting sparsity in channel and data estimation for sporadic multi-user communication,” in _ISWCS 2013; The Tenth International Symposium on Wireless Communication Systems_. VDE, 2013, pp. 1–5.
  * [11] G. Wunder, H. Boche, T. Strohmer, and P. Jung, “Sparse signal processing concepts for efficient 5G system design,” _IEEE Access_ , vol. 3, pp. 195–208, 2015.
  * [12] G. Hannak, M. Mayer, A. Jung, G. Matz, and N. Goertz, “Joint channel estimation and activity detection for multiuser communication systems,” in _2015 IEEE International Conference on Communication Workshop (ICCW)_. IEEE, 2015, pp. 2086–2091.
  * [13] Z. Chen, F. Sohrabi, and W. Yu, “Sparse activity detection for massive connectivity,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 7, pp. 1890–1904, 2018.
  * [14] W. Chen, H. Xiao, L. Sun, and B. Ai, “Joint activity detection and channel estimation in massive mimo systems with angular domain enhancement,” _IEEE Transactions on Wireless Communications_ , vol. 21, no. 5, pp. 2999–3011, 2021.
  * [15] M. Bayati and A. Montanari, “The dynamics of message passing on dense graphs, with applications to compressed sensing,” _IEEE Transactions on Information Theory_ , vol. 57, no. 2, pp. 764–785, 2011.
  * [16] Z. Chen and W. Yu, “Massive device activity detection by approximate message passing,” in _2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2017, pp. 3514–3518.
  * [17] J. Johnston and X. Wang, “Model-based deep learning for joint activity detection and channel estimation in massive and sporadic connectivity,” _IEEE Transactions on Wireless Communications_ , vol. 21, no. 11, pp. 9806–9817, 2022.
  * [18] X. Ma, Z. Gao, F. Gao, and M. Di Renzo, “Model-driven deep learning based channel estimation and feedback for millimeter-wave massive hybrid mimo systems,” _IEEE Journal on Selected Areas in Communications_ , vol. 39, no. 8, pp. 2388–2406, 2021.
  * [19] G. Tang, B. N. Bhaskar, P. Shah, and B. Recht, “Compressed sensing off the grid,” _IEEE transactions on information theory_ , vol. 59, no. 11, pp. 7465–7490, 2013.
  * [20] D. L. Donoho, “Compressed sensing,” _IEEE Trans. Inf. Theory_ , vol. 52, no. 4, pp. 1289–1306, 2006.
  * [21] E. J. Candes, Y. C. Eldar, D. Needell, and P. Randall, “Compressed sensing with coherent and redundant dictionaries,” _Appl. Comput. Harmon. Anal._ , vol. 31, no. 1, pp. 59–73, 2011.
  * [22] L. Stanković, I. Orović, S. Stanković, and M. Amin, “Compressive sensing based separation of nonstationary and stationary signals overlapping in time-frequency,” _IEEE Transactions on Signal Processing_ , vol. 61, no. 18, pp. 4562–4572, 2013.
  * [23] C. Studer, P. Kuppinger, G. Pope, and H. Bolcskei, “Recovery of sparsely corrupted signals,” _IEEE Transactions on Information Theory_ , vol. 58, no. 5, pp. 3115–3130, 2012.
  * [24] Y. Li, X. Wang, and Z. Ding, “Multi-target position and velocity estimation using ofdm communication signals,” _IEEE Transactions on Communications_ , vol. 68, no. 2, pp. 1160–1174, 2019.
  * [25] Y. Li, L. Zheng, M. Lops, and X. Wang, “Interference removal for radar/communication co-existence: The random scattering case,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 10, pp. 4831–4845, 2019.
  * [26] Z. Yang and L. Xie, “Exact joint sparse frequency recovery via optimization methods,” _IEEE Trans. Signal Process._ , vol. 64, no. 19, pp. 5145–5157, 2014.
  * [27] B. N. Bhaskar, G. Tang, and B. Recht, “Atomic norm denoising with applications to line spectral estimation,” _IEEE Trans. Signal Process._ , vol. 61, no. 23, pp. 5987–5999, Dec. 2013.
  * [28] Y. Li, X. Wang, and Z. Ding, “Multi-dimensional spectral super-resolution with prior knowledge with application to high mobility channel estimation,” _IEEE Journal on Selected Areas in Communications_ , 2020.
  * [29] C. A. Metzler, A. Maleki, and R. G. Baraniuk, “From denoising to compressed sensing,” _IEEE Transactions on Information Theory_ , vol. 62, no. 9, pp. 5117–5144, 2016.
  * [30] S. Ramani, T. Blu, and M. Unser, “Monte-Carlo SURE: A black-box optimization of regularization parameters for general denoising algorithms,” _IEEE Transactions on Image Processing_ , vol. 17, no. 9, pp. 1540–1554, 2008.
  * [31] A. Alkhateeb, O. El Ayach, G. Leus, and R. W. Heath, “Channel estimation and hybrid precoding for millimeter wave cellular systems,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 8, no. 5, pp. 831–846, 2014.
  * [32] Z. Xiao, P. Xia, and X.-G. Xia, “Channel estimation and hybrid precoding for millimeter-wave mimo systems: A low-complexity overall solution,” _IEEE Access_ , vol. 5, pp. 16 100–16 110, 2017.
  * [33] J. Kim, W. Chang, B. Jung, D. Baron, and J. C. Ye, “Belief propagation for joint sparse recovery,” _arXiv preprint arXiv:1102.3289_ , 2011.
  * [34] L. Zheng, A. Maleki, H. Weng, X. Wang, and T. Long, “Does $\ell_{p}$-minimization outperform $\ell_{1}$-minimization?” _IEEE Transactions on Information Theory_ , vol. 63, no. 11, pp. 6896–6935, 2017.
  * [35] A. Javanmard and A. Montanari, “State evolution for general approximate message passing algorithms, with applications to spatial coupling,” _Information and Inference: A Journal of the IMA_ , vol. 2, no. 2, pp. 115–144, 2013.
  * [36] T. L. Marzetta, “Noncooperative cellular wireless with unlimited numbers of base station antennas,” _IEEE Transactions on Wireless Communications_ , vol. 9, no. 11, pp. 3590–3600, 2010.
  * [37] C. Guo and M. E. Davies, “Near optimal compressed sensing without priors: Parametric sure approximate message passing.” _IEEE Trans. Signal Processing_ , vol. 63, no. 8, pp. 2130–2141, 2015.
  * [38] L. Zheng, Z. Wu, M. Seok, X. Wang, and Q. Liu, “High-accuracy compressed sensing decoder based on adaptive ($\ell_{0}$,$\ell_{1}$) complex approximate message passing: Cross-layer design,” _IEEE Transactions on Circuits and Systems I: Regular Papers_ , vol. 63, no. 10, pp. 1726–1736, 2016.
  * [39] S. Boyd and L. Vandenberghe, _Convex optimization_. Cambridge, U.K.: Cambridge Univ. Press, 2004.
  * [40] J. F. Sturm, “Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones,” _Optimization Methods and Software_ , vol. 11, no. 1-4, pp. 625–653, 1999.
  * [41] K.-C. Toh, M. J. Todd, and R. H. Tütüncü, “SDPT3-—a MATLAB software package for semidefinite programming, version 1.3,” _Optimization Methods and Software_ , vol. 11, no. 1-4, pp. 545–581, 1999.
  * [42] A. Eftekhari and M. B. Wakin, “Greed is super: A fast algorithm for super-resolution,” _arXiv preprint arXiv:1511.03385_ , 2015.
  * [43] A. Fannjiang and W. Liao, “Coherence pattern–guided compressive sensing with unresolved grids,” _SIAM Journal on Imaging Sciences_ , vol. 5, no. 1, pp. 179–202, 2012.
  * [44] L. Zheng, M. Lops, and X. Wang, “Adaptive interference removal for uncoordinated radar/communication coexistence,” _IEEE J. Sel. Top. Signal Proces._ , vol. 12, no. 1, pp. 45–60, 2018.
  * [45] F. Pukelsheim, “The three sigma rule,” _The American Statistician_ , vol. 48, no. 2, pp. 88–91, 1994.
  * [46] M. Viberg, B. Ottersten, and T. Kailath, “Detection and estimation in sensor arrays using weighted subspace fitting,” _IEEE Trans. Signal Process._ , vol. 39, no. 11, pp. 2436–2449, 1991.
  * [47] M. Viberg and B. Ottersten, “Sensor array processing based on subspace fitting,” _IEEE Trans. Signal Process._ , vol. 39, no. 5, pp. 1110–1121, 1991.