# Performance Analysis of Joint Active User Detection and Channel Estimation
for Massive Connectivity ††thanks: Jia-Cheng Jiang and Hui-Ming Wang are with
the School of Information and Communications Engineering, Xi'an Jiaotong
University, Xi'an, 710049, Shaanxi, P. R. China, and also with the Ministry of
Education Key Lab for Intelligent Networks and Network Security, Xi'an,
710049, Shaanxi, P. R. China. (e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>††thanks: H. Vincent Poor is with the Department of
Electrical Engineering, Princeton University, Princeton, NJ 08544 USA (e-mail:
<EMAIL_ADDRESS>
Jia-Cheng Jiang, Hui-Ming Wang, _Senior Member_ , _IEEE_ and H. Vincent Poor,
_Life Fellow_ , _IEEE_

###### Abstract

This paper considers joint active user detection (AUD) and channel estimation
(CE) for massive connectivity scenarios with sporadic traffic. The state-of-
art method under a Bayesian framework to perform joint AUD and CE in such
scenarios is approximate message passing (AMP). However, the existing
theoretical analysis of AMP-based joint AUD and CE can only be performed with
a given fixed point of the AMP state evolution function, lacking the analysis
of AMP phase transition and Bayes-optimality. In this paper, we propose a
novel theoretical framework to analyze the performance of the joint AUD and CE
problem by adopting the replica method in the Bayes-optimal condition.
Specifically, our analysis is based on a general channel model, which reduces
to particular channel models in multiple typical MIMO communication scenarios.
Our theoretical framework allows ones to measure the optimality and phase
transition of AMP-based joint AUD and CE as well as to predict the
corresponding performance metrics under our model. To reify our proposed
theoretical framework, we analyze two typical scenarios from the massive
random access literature, i.e., the isotropic channel scenario and the
spatially correlated channel scenario. Accordingly, our performance analysis
produces some novel results for both the isotropic Raleigh channel and
spatially correlated channel case.

##  I Introduction

Recently, under the basic characteristics of mMTC, i.e, a large number of user
devices and sporadic user traffic, the grant-free access strategy has been
considered to allow the active devices to access the wireless network without
a grant [1, 2], which reduces both the access latency and signal processing
overhead. Accordingly, the base station (BS) should simultaneously identify
all the active user devices under grant-free random access. Furthermore, the
BS is required to accurately acquire the channel state information for
decoding uplink signals and executing downlink precoding after user activity
detection. Hence, both active user detection (AUD) and channel estimation (CE)
are required at the BS based on pilot sequences sent by the user devices.
Typically, the number of user devices in mMTC is very large, and thus
assigning orthogonal pilot sequences to all user devices is prohibitive in
practice, motivating the application of non-orthogonal pilot sequences in such
situations. A central problem in the mMTC scenario is to jointly perform AUD
and CE based on non-orthogonal pilot sequences.

Different from the sparse detection problem [3, 4], the joint AUD and CE
requires both the support set and the corresponding amplitudes of the sparse
vector, and can be formulated as a compressed sensing (CS) problem. Several
CS-based solutions have been reported [5, 6, 7, 8, 9]. In [5], a dimension
reduction method to reduce the pilot sequence length and computational
complexity for joint AUD and CE has been proposed, which projects the original
device state matrix onto a low-dimensional space by exploiting its sparse and
low-rank structure. To fully utilize the statistical information of wireless
channels, some Bayesian-based CS methods have been designed to achieve joint
AUD and CE with higher accuracy. Specifically, the authors in [6] proposed an
approximate message passing (AMP)-based joint AUD and CE strategy based on the
randomly generated non-orthogonal pilot sequences for a massive MIMO system,
and the statistical knowledge of both the channel and user sparsity was
modeled using a prior distribution to facilitate the performance of AMP. Along
this line, the authors in [7] design an AMP-based algorithm to adaptively
detect the active devices by exploiting the virtual angular domain sparsity of
the channels in an orthogonal frequency division multiplexing (OFDM) broadband
system. Furthermore, an expectation propagation (EP)-based joint AUD and CE
algorithm was proposed in [8] for massive access with a single-antenna BS,
where the computationally intractable posterior probability of the involved
sparse signals was approximated by a multivariate Gaussian distribution,
enhancing the performance of joint AUD and CE. In addition, we have claimed in
[9] that the inherent temporal correlations of both the user channels and the
active indicators between adjacent time slots can be used to enhance both the
detection performance and channel estimation performance.

One critical issue for joint AUD and CE for massive connectivity is to analyze
the performance of user detection and channel estimation. The AMP algorithm
provides a corresponding theoretical framework called state evolution to
accurately track the performance of AMP in each iteration. By using state
evolution analysis, the authors in [10, 6] provide analytical
characterizations of the missed detection and false alarm probabilities for
device detection and channel error covariance for channel estimation under a
massive MIMO Raleigh channel assumption. However, although state evolution
provides a theoretical framework for analyzing the performance metrics of
joint AUD and CE, there are some critical issues that are left untouched.
First, the performance analysis is based on the fixed point of the state
evolution function, but the connection between the performance and the
specific choices of system parameters such as pilot length, transmit power,
fraction of number of active users and number of antennas was not established.
Second, for specific choices of system parameters, the AMP iterations are
blocked in a sub-optimal fixed point, so that the Bayes-optimal AUD and CE
performance cannot be achieved via the AMP framework [11]. The region of the
system parameters where the AMP framework is sub-optimal cannot be measured
via the state evolution analysis in [10, 6]. Finally, some related works [12,
13] have shown that the performance of AMP exhibits a _phase transition_
phenomenon where the AMP algorithm will exhibit disconnect performance
variations with variations of the system parameters, which has not been
analyzed for massive connectivity.

To address the above problems, in this paper, we propose a theoretical
framework to analyze the performance of joint AUD and CE based on the replica
method [14, 15], which is a classical tool for analyzing large systems that
comes primarily from physics. The proposed framework established the
connection between the system parameters and the performance of joint AUD and
CE. Based on that, both the Bayes-optimal and the AMP achievable mean square
error (MSE) can be predicted, and the region of system parameters where the
AMP algorithm is sub-optimal can be determined. In addition, the phase
transition phenomenon of our system can be analyzed, which guides the system
design for massively connected networks.

Note that, different from [6], where the performance analysis is performed
based on the isotropic Raleigh channel model, the analysis in this paper is
considered based on a more general channel model following the approach of
[16]. Such a channel model has strong flexibility and reduces to particular
channel models in multiple typical MIMO communication scenarios, such as the
isotropic Raleigh channel model [10], and the spatially correlated channel
model [7, 17, 18], etc. In particular, our performance analysis produces some
novel results for both the isotropic Raleigh channel and spatially correlated
channel case, and provides verifications for the analytical results in [10,
6]. Our main contributions can be summarized as follow.

  * •

First, we consider a general grouping channel model for the joint AUD and CE
scenario. We design a theoretical framework carefully tailored to our
considered scenario based on the general idea of replica method. Based on
that, we establish relations between the joint AUD and CE performance metrics
and the system parameters: pilot length, transmit power, fraction of number of
active users and number of BS antennas, under our considered channel model.

  * •

Second, we analyze the isotropic channel scenario based on our theoretical
framework. Concretely, we prove that the Bayes-optimal/AMP-achievable joint
AUD and CE performance can be evaluated by a scalar valued function, which
also provides a phase transition diagram. We further provide the analysis in
the asymptotic MIMO regime and prove that the phase transition phenomenon
disappears and the AMP can always achieve Bayes-optimal performance when the
number of BS antennas is very large.

  * •

Third, we analyze the performance of joint AUD and CE in the spatially
correlated channel case. We prove that the performance of each user group can
be separately evaluated by a particular free entropy function when the user
groups are in mutually orthogonal subspaces, implying the per-group processing
(PGP) will never bring the performance loss. In addition, we show that the
spatially correlated channel is more likely to promote a phase transition
compared with the isotropic channel, due to the small number of subspace
dimensions and corresponding antenna gains, leading that there will be a
performance gap between AMP and Bayes-optimal performance. However, the pilot
length required to step over the phase transition can be significantly reduced
in that case.

Figure 1: Our model of the massive device communication network.

##  II System Model

###  II-A Massive Random Access Scenario

We consider a uplink massive random access scenario in a single-cell cellular
network with $N$ user devices. Each user is equipped with a single antenna,
and the BS is equipped with $M$ antennas. The user traffic is assumed to be
sporadic, i.e., most of user devices are idle at any given time.

Following the approach of [16], in this paper, we consider the following
channel model. We suppose that $N$ users are divided into $G$ groups based on
the similarity of their covariance matrixes, where each group has $K_{g}$
users and the number of users is assumed same in each group for concise, i.e.,
$K_{g}=K=N/G$. Our model of the massive device communication network is shown
in Fig. 1. We denote $\bm{h}_{gk}\in\mathbb{C}^{M\times 1}$ as the channel
response of $k$th user in the group $g$. We further assume that users in the
same group have an identical channel probability density function (PDF)
denoted as
$Q_{g}(\bm{h}_{kg})=\mathcal{CN}(\bm{h}_{kg};\bm{0},\mathcal{C}_{g})$ with
covariance matrix $\mathcal{C}_{g}$. For the $k$th user in the group $g$, the
user state, i.e., active or idle, is characterized by an activity indicator,
denoted as $a_{gk}$, i.e., $a_{gk}=1$ if it is active and $a_{gk}=0$
otherwise. We denote $\bm{f}_{kg}\in\mathbb{C}^{T\times 1}$ as the non-
orthogonal pilot sequence that is assigned to the $k$th user in the group $g$
with independent and identically distributed (i.i.d.) components generated
from a complex Gaussian distribution with zero mean and variance $1/T$, so
that the pilot sequence has a unit norm111In the asymptotic regime that we
consider, the pilot sequences have an unit power.. Accordingly, the overall
channel input-output is

$\bm{Y}=\bm{F}\bm{S}+\bm{W}=\sum\limits_{g,k}a_{gk}\bm{f}_{gk}\bm{h}_{gk}^{T}+\bm{W},$
(1)

where $\bm{F}\in\mathbb{C}^{T\times
N}\triangleq[\bm{f}_{1},\dots,\bm{f}_{GK}]$ is the pilot matrix,
$\bm{S}\in\mathbb{C}^{N\times M}\triangleq[\bm{s}_{11},\dots,\bm{s}_{GK}]^{T}$
with $\bm{s}_{gk}\in\mathbb{C}^{N\times 1}\triangleq a_{gk}\bm{h}_{gk}$, and
$\bm{W}$ is the additive white Gaussian noise (AWGN) matrix with i.i.d.
elements distributed as $\mathcal{CN}(0,\sigma_{w}^{2})$. Each user device has
a constant transmit power $P_{t}$, which is absorbed into the large-scale
fading of the channel coefficients in (1) for concise.

We assume that the user devices are synchronized and each user decides whether
or not to access the channel with a probability $\rho$ in an i.i.d. manner
[10, 6, 7, 19], and we further assume the fraction of active users $\rho$ and
the PDF $Q_{g}(\bm{h}_{kg})$ of the $k$th user channel in the group $g$ are
known for the BS222Note that the parameters learning strategy under mMTC
scenario has been considered in the related work [7, 20] based on the
expectation-maximization method, which can be easily extended under our model.
Since we focus on the theoretical analysis under our framework, the parameters
learning issue is beyond the scope of this paper. As a consequence, the
corresponding PDF of the matrix $\bm{S}$ can be formulated as

$p(\bm{S})=\prod_{g=1}^{G}\prod_{k=1}^{K}[(1-\rho)\delta(\bm{s}_{gk})+\rho
Q_{g}(\bm{s}_{gk})].$ (2)

Based on our model (1) and the prior distribution of the transmitted signals
(2)333Note that we assume joint AUD and CE is performed within one coherence
time so that the channel coefficients remain unchanged., our task is to
identify the active users, i.e., determine the active indicator $a_{gk}$ for
each user, as well as estimate the corresponding channel coefficients
$\bm{h}_{gk}$ for the active users.

_Remark_ : Note that our channel model assumption has a strong flexibility and
reduces to particular channel models in multiple typical MIMO communication
scenarios, such as the isotropic Rayleigh fading channels [10, 6], where the
channel coefficients of each user device are i.i.d., and the spatially
correlated channels [7, 17, 18], where different user groups are sufficiently
well separated in the angular domain, and the channel covariance matrixes
exhibit a low rank structure. We also note that although we consider the case
where each group has the same number of users and the fraction of active users
for each group is also the same, our following proposed theoretical framework
can be easily extended to the case where the active ratio and user number are
different for each group.

Based on our general channel model (2), in the following, we propose a novel
theoretical framework by adopting replica method [14, 15] to provide phase
transition analysis, Bayes-optimality analysis of our scenario, and our
theoretical framework can also provide the performance prediction of the AMP-
based joint AUD and CE, which is the state-of-art algorithm under the Bayesian
framework.

##  III Replica Analysis on Joint AUD and CE

Our theoretical framework is based on statistical physics. In the statistical
physics literature, the _free entropy_ of a system reflects the macro
performance that characterizes some thermodynamic properties of the system
[15]. Some related works in the signal processing literature have also shown
that evaluating the fixed point of free entropy function provides the minimum
MSE (MMSE) prediction for signal recovery and the achievable MSE for the AMP
framework [12, 13]. For the massive connectivity scenario, evaluating the free
entropy function of our system (1) provides an analytic tool to measure the
performances metrics of joint AUD and CE.

Towards this end, we adopt the replica method in this paper to evaluate the
free entropy function of our system. The replica method [14, 15], which is a
classical tool in the statistical physics literature, provides an efficiency
way to reduce the calculation of a certain free entropy into an optimization
problem over specific covariance matrices. Such a reduction is based on a set
of typical replica assumptions that include the self-averaging property, the
validity of ``replica trick'', the ability to exchange certain limits and the
replica symmetry [11]. It has been shown in many related works that the
replica method provides correct predictions on many typical applications, such
as signal processing [21] and physical layer communications [22, 23].

Our analysis based on the replica method framework [11] is considered under a
certain _asymptotic regime_ with $K\rightarrow\infty$, $T\rightarrow\infty$
and a fixed ratio $T/K\rightarrow\alpha$. Next, we provide the background of
the statistical physics literature and obtain the free entropy function with
respect to the system parameters under our scenario via replica method.
Although these system parameters can not be infinity in the practice, we still
utilize the asymptotic regime to facilitate analysis, since the system
parameters are typically large in the massive connectivity. Authors in [6]
have shown that the recovery performances in the practical settings match the
theoretical results in the large system limit.

$\begin{split}\Phi(\mathcal{E})=&-{\rm
Tr}\left(\left(\bm{\Delta}+\frac{\mathcal{E}G}{\alpha}\right)^{-1}\left(\sum\limits_{g}\rho\mathcal{C}_{g}-\mathcal{E}G\right)\right)-\alpha
M-\alpha\log\left|\bm{\Delta}+\frac{\mathcal{E}G}{\alpha}\right|\\\
&+\sum\limits_{g}\int{\rm d}\bm{s}_{g}[(1-\rho)\delta(\bm{s}_{g})+\rho
Q_{g}(\bm{s}_{g})]\int{\rm D}\bm{z}\log\left\\{\int{\rm
d}\bm{x}_{g}[(1-\rho)\delta(\bm{x}_{g})+\rho Q_{g}(\bm{x}_{g})]\right.\\\
&\left.\exp\left\\{-\bm{x}_{g}^{H}(\bm{\Delta}+\frac{\mathcal{E}G}{\alpha})^{-1}\bm{x}_{g}+2\mathfrak{R}\left(\bm{x}_{g}^{H}[(\bm{\Delta}+\frac{\mathcal{E}G}{\alpha})^{-1}\bm{s}_{g}+(\bm{\Delta}+\frac{\mathcal{E}G}{\alpha})^{-\frac{1}{2}}\bm{z}]\right)\right\\}\right\\}.\end{split}$
(7)

###  III-A Free Entropy Function Derivation via Replica Method

Under the statistical physics literature, a probabilistic inference approach
to reconstruct the transmitted signal $\bm{S}$ aims to sample $\bm{X}$ from
the following posterior distribution based on (1), given by

$p(\bm{X}|\bm{Y})=\frac{p(\bm{X})}{Z}\prod_{t=1}^{T}\frac{\exp\left(-(\bm{y}_{t}-\bm{f}_{t}\bm{X})\bm{\Delta}^{-1}(\bm{y}_{t}-\bm{f}_{t}\bm{X})^{H}\right)}{\pi^{M}|\bm{\Delta}|},$
(3)

where $\bm{\Delta}\triangleq\sigma_{w}^{2}\bm{I}$, $Z$ is the partition
function of this distribution, and
$p(\bm{X})=\prod_{g,k}[(1-\rho)\delta(\bm{x}_{gk})+\rho Q_{g}(\bm{x}_{gk})]$,
which can be regarded as the prior distribution of $\bm{X}$. Under the Bayes-
optimal assumption, the prior distribution over $\bm{X}$ matches the PDF of
transmitted signal $p(\bm{S})$. As a consequence, $Z$ can be formulated as

$\displaystyle Z=\int{\rm d}\bm{X}$
$\displaystyle\prod_{t=1}^{T}\frac{1}{\pi^{M}|\bm{\Delta}|}e^{-(\bm{y}_{t}-\bm{f}_{t}\bm{X})\bm{\Delta}^{-1}(\bm{y}_{t}-\bm{f}_{t}\bm{X})^{H}}$
$\displaystyle\times\prod_{g=1}^{G}\prod_{k=1}^{K}[(1-\rho)\delta(\bm{x}_{gk})+\rho
Q_{g}(\bm{x}_{gk})],$ (4)

where $\bm{f}_{t}$ is the $t$th row of the pilot matrix $\bm{F}$. We note that
above integration is performed over each element of matrix $\bm{X}$ and such a
integral expression will be used throughout this paper. According to the
statistical physics literature, the free entropy is defined as
$\Phi\triangleq\log Z$. Utilizing the self-averaging assumption [14, 15]:
$\lim\limits_{N\to\infty}\Pr\left[\left|\frac{Z}{N}-\frac{\mathbb{E}(Z)}{N}\right|\geq\theta\right]=0$,
for any tolerance $\theta>0$, the free entropy can be reformulated as
$\Phi=\frac{1}{N}\mathbb{E}_{\bm{F},\bm{S},\bm{W}}(\log Z)$. Hence, to
determine the free entropy function, one requires to compute the average of
the log-partition function $\mathbb{E}_{\bm{F},\bm{S},\bm{W}}(\log Z)$. Based
on the standard procedure of the replica method [11], we can alternatively
determine the average log-partition function by the following relation
$\mathbb{E}_{\bm{F},\bm{S},\bm{W}}(\log Z)=\lim\limits_{n\to
0}\frac{1}{n}\log\left(\mathbb{E}_{\bm{F},\bm{S},\bm{W}}Z^{n}\right)$.
According to the _replica trick_ in [11, 14, 15], the quantity
$\mathbb{E}_{\bm{F},\bm{S},\bm{W}}Z^{n}$ is carried out as if $n$ were an
integer, and we take the fact that $n$ is a real number into consideration
after obtaining a manageable enough expression. As a consequence, the free
entropy can be formulated as

$\Phi=\lim\limits_{N\to\infty}\lim\limits_{n\to
0}\frac{1}{Nn}\log\left(\mathbb{E}_{\bm{F},\bm{S},\bm{W}}Z^{n}\right).$ (5)

Based on this, we have the following theorem.

###### Theorem 1

Define the reconstruction MSE matrix over the system (1) as
$\mathcal{E}\triangleq\frac{1}{G}\sum_{g}\mathcal{E}_{g}$, where

$\mathcal{E}_{g}\triangleq\frac{1}{K}\sum\limits_{k=1}^{K}\left(\hat{\bm{x}}_{gk}(\bm{Y})-\bm{s}_{gk}\right)\left(\hat{\bm{x}}_{gk}(\bm{Y})-\bm{s}_{gk}\right)^{H},$
(6)

where
$[\hat{\bm{x}}_{11}(\bm{Y}),\dots,\hat{\bm{x}}_{GK}(\bm{Y})]=\hat{\bm{X}}$ is
the sample of the posterior distribution $p(\bm{X}|\bm{Y})$ with given
$\bm{Y}$ defined in (3). Under our certain asymptotic regime defined in Sec.
III, and using the standard assumptions of the replica method, i.e., the self-
averaging assumption, the replica-symmetry assumption, and the replica trick
in (5), the free entropy (5) can be calculated by optimizing a function with
respect to $\mathcal{E}$, which is given in the equation (7). We note that the
term ${\rm D}\bm{z}$ is a complex Gaussian integration measure.

###### Proof 1

Please see the appendix -A.

We call the function (7) as _free entropy function_ 444Since the free entropy
function is closely related to the free entropy, we also use $\Phi$ to
represent it. . According to [12, 13], optimizing $\mathcal{E}$ that maximizes
the free entropy function (7) corresponds to the minimum MSE (MMSE) of
specific choices of system parameters, i.e., $\rho$, $\alpha$,
$\mathcal{C}_{g}$, $\forall g$, in the Bayes-optimal condition. Setting the
first order derivative of $\Phi(\mathcal{E})$ with respect to the matrix
$\mathcal{E}$ into zero, we get the following equation.

$\displaystyle\mathcal{E}=\frac{1}{G}\sum\limits_{g}$
$\displaystyle\mathbb{E}_{\bm{s}_{g},\bm{z}}\left[\left(\eta_{g}\left(\bm{s}_{g}+(\bm{\Delta}+\frac{\mathcal{E}G}{\alpha})^{\frac{1}{2}}\bm{z}\right)-\bm{s}_{g}\right.\right)$
$\displaystyle\times\left(\left.\eta_{g}\left(\bm{s}_{g}+(\bm{\Delta}+\frac{\mathcal{E}G}{\alpha})^{\frac{1}{2}}\bm{z}\right)-\bm{s}_{g}\right)^{H}\right],$
(8)

where $\bm{z}$ obeys $\mathcal{CN}(\bm{z};\bm{0},\bm{I})$ and
$\eta_{g}(\cdot)$ is the Bayes-optimal denoiser of the noisy measurement
$\hat{\bm{s}}_{g}\triangleq\bm{s}_{g}+(\bm{\Delta}+\frac{\mathcal{E}G}{\alpha})^{\frac{1}{2}}\bm{z}$
with the latent signal $\bm{s}_{g}$ distributed as
$(1-\rho)\delta(\bm{s}_{g})+\rho Q_{g}(\bm{s}_{g})$. We note that
$\eta_{g}(\cdot)$ is the corresponding MMSE denoiser of noisy measurement
$\hat{\bm{s}}_{g}$.

We can also observe from (8) that there is a close match between the fixed
point of the AMP state evolution function in [24] and the stationary point of
the free energy function (7). In addition, according to [12, 13], the AMP
algorithm will be blocked by a particular local maximum point of (7). This
means that seeking the stationary points of the free entropy function (7) also
provides accurate prediction of the performance of the AMP algorithm. We
remark that the recovery MSE of AMP is sometimes sub-optimal since only the
global maximum point that maximizes the free entropy function (7) corresponds
to the MMSE, and we call it AMP-achievable MSE.

In the following, we establish connections between the fixed point of the AMP
state evolution function predicted by the free entropy function (7), and the
performance metrics of joint AUD and CE. Typically, we adopt the likelihood
ratio test (LRT) detection for AUD, and adopt MMSE criterion for CE. We note
that the author in [25] claims that, in the case of large i.i.d. zero-mean
Gaussian sensing matrix, the AMP methods exhibit fast convergence. Since in
the massive access literature, the system is large enough and the i.i.d. zero-
mean Gaussian sensing matrix is adopted, guaranteeing the convergence of the
AMP algorithm in the following sections.

$\begin{split}\Phi(\tau)=&-\alpha
M\left(\frac{\sigma_{w}^{2}}{\frac{1}{\alpha}\tau+\sigma_{w}^{2}}+\log(\sigma_{w}^{2}+\frac{1}{\alpha}\tau)\right)+M\sum\limits_{g}\frac{(1-\rho)\sigma_{g}^{2}}{\sigma_{g}^{2}+\sigma_{w}^{2}+\frac{1}{\alpha}\tau}\\\
&+\sum\limits_{g}\int{\rm
D}\bm{z}\rho\log\left[(1-\rho)\exp\left\\{-\frac{||\bm{z}||^{2}\sigma_{g}^{2}}{\sigma_{w}^{2}+\frac{1}{\alpha}\tau}\right\\}+\rho\left(\frac{\frac{1}{\alpha}\tau+\sigma_{w}^{2}}{\frac{1}{\alpha}\tau+\sigma_{w}^{2}+\sigma_{g}^{2}}\right)^{M}\right]\\\
&+\sum\limits_{g}\int{\rm
D}\bm{z}(1-\rho)\log\left[(1-\rho)\exp\left\\{-\frac{||\bm{z}||^{2}\sigma_{g}^{2}}{\sigma_{g}^{2}+\sigma_{w}^{2}+\frac{1}{\alpha}\tau}\right\\}+\rho\left(\frac{\frac{1}{\alpha}\tau+\sigma_{w}^{2}}{\frac{1}{\alpha}\tau+\sigma_{w}^{2}+\sigma_{g}^{2}}\right)^{M}\right].\end{split}$
(14)

###  III-B Prediction of AMP-based Joint AUD and CE

After executing the AMP algorithm on the transmitted signal $\bm{Y}$ in (1),
the overall estimation problem is decoupled as a sequence of vector-valued
estimation problems. For the users in the group $g$, the decoupled signal
model is
$\hat{\bm{s}}_{g}\triangleq\bm{s}_{g}+\bm{\Sigma}^{\frac{1}{2}}\bm{z}$, where
we define
$\bm{\Sigma}\triangleq\bm{\Delta}+\frac{\mathcal{E}^{\star}G}{\alpha}$ as the
equivalent noise covariance matrix. The $\mathcal{E}^{\star}$ denotes the
fixed point of AMP state evolution function that fulfills equation (8), so
that it is also a stationary point of the free entropy function (7). Note that
we have dropped the subscript about the user indicator $k$, much as the
following expressions, since all the users within a group share a common
probabilistic model.

As a consequence, for the user devices in the group $g$, the AUD rule based on
log-likelihood ratio (LLR) with a threshold $l_{g}$ can be formulated as

$\displaystyle\text{LLR}(\hat{\bm{s}}_{g})=$
$\displaystyle\log\left(\frac{p(\hat{\bm{s}}_{g}|a_{g}=1)}{p(\hat{\bm{s}}_{g}|a_{g}=0)}\right)$
$\displaystyle=$
$\displaystyle\hat{\bm{s}}_{g}^{H}\tilde{\bm{\Sigma}}\hat{\bm{s}}_{g}+\log\frac{\left|\bm{\Sigma}\right|}{\left|\mathcal{C}_{g}+\bm{\Sigma}\right|}>l_{g},$
(9)

where we define
$\tilde{\bm{\Sigma}}\triangleq\bm{\Sigma}^{-1}-\left(\mathcal{C}_{g}+\bm{\Sigma}\right)^{-1}$.
Accordingly, the missed detection and false alarm probabilities $P_{M}$ and
$P_{F}$ can be obtained via

$\displaystyle
P_{F}=\Pr\left\\{\hat{\bm{s}}_{g}^{H}\tilde{\bm{\Sigma}}\hat{\bm{s}}_{g}>l_{g}^{\prime}|a_{g}=0\right\\},$
(10) $\displaystyle
P_{M}=\Pr\left\\{\hat{\bm{s}}_{g}^{H}\tilde{\bm{\Sigma}}\hat{\bm{s}}_{g}<l_{g}^{\prime}|a_{g}=1\right\\},$
(11)

where $l_{g}^{\prime}\triangleq
l_{g}-\log\left|\bm{\Sigma}\right|+\log\left|\mathcal{C}_{g}+\bm{\Sigma}\right|$.

We then consider the performance analysis for CE. Since the CE for the active
users are performed after the AUD, we thus consider the mean of conditional
posterior distribution $p(\bm{s}_{g}|\hat{\bm{s}}_{g},a_{g}=1)$ as the channel
estimator for active user devices. Based on the definition of noisy
measurement $\hat{\bm{s}}_{g}$, after performing the AMP algorithm, the CE for
the active user devices in the group $g$ is

$\hat{\bm{h}}_{g}=\left(\bm{\Sigma}^{-1}+\mathcal{C}_{g}^{-1}\right)^{-1}\bm{\Sigma}^{-1}\hat{\bm{s}}_{g}.$
(12)

For the ease of analysis, we only consider the CE error in the case where the
AUD is executed perfectly. Under such a circumstance, the CE error matrix can
be formulated as

$\displaystyle\mathbb{E}\left[\left(\bm{s}_{g}-\hat{\bm{h}}_{g}\right)\left(\bm{s}_{g}-\hat{\bm{h}}_{g}\right)^{H}\right]=\left(\mathcal{C}_{g}^{-1}+\bm{\Sigma}^{-1}\right)^{-1},$
(13)

where the expectation is taken over $p(\bm{s}_{g},\hat{\bm{s}}_{g}|a_{g}=1)$.
As a consequence, the performance of joint AUD and CE can be analyzed after
obtaining the stationary point $\mathcal{E}^{\star}$ of the free entropy (7).
We note that the equation (7) and the corresponding performance of joint AUD
and CE in this section is based on the general channel model (2). In the
following sections, we will show our theoretical framework can reduce to two
typical scenarios in the massive random access literature, i.e., the isotropic
channel scenario and the spatially correlated channel scenario. Concretely, in
Sec IV, we provide the performance analysis of the isotropic channel case.
While in the Sec. V, we provide the performance analysis of the spatially
correlated channel case.

##  IV Isotropic Channel

In this section, we consider a typical scenario in wireless communications,
where there are many small reflectors around the BS, and the Rayleigh fading
MIMO channel with isotropic channel is assumed [10, 6]. In this case, the
Rayleigh channel components of each user group are assumed to be i.i.d.
complex Gaussian with zero means and unit variances across all the BS
antennas, and function $Q_{g}(\bm{h}_{g})$ in the prior distribution over the
channel in the $g$th user group is
$Q_{g}(\bm{h}_{g})=\mathcal{CN}(\bm{h}_{g};0,\sigma_{g}^{2}\bm{I})$, where
$\sigma_{g}^{2}$ is the product of the large scale coefficient and the user
transmitted power. In the following, we will analyze the performance of joint
AUD and CE under this case by first deriving the corresponding free entropy
function as well as its expression under an asymptotic massive MIMO regime.
Then, the performance analysis of joint AUD and CE is given.

###  IV-A Free Entropy Function

Before deriving the concrete free entropy function, we note that it is proved
in [6] that the equivalent noise covariance matrix in the AMP state evolution
iterations always in the form of a diagonal matrix with identical diagonal
entries in the isotropic Rayleigh channel case, so that the fixed point of the
state evolution function (8) remains the same form. We therefore restrict the
domain of $\mathcal{E}$ in the form of $\mathcal{E}=\tau G^{-1}\mathbf{I}$,
and the newly derived free entropy function with respect to $\tau$ is in the
following theorem.

###### Theorem 2

With the assumption of the isotropic Rayleigh channel, i.e.,
$Q_{g}(\bm{s}_{g})=\mathcal{CN}(\bm{s}_{g};0,\sigma_{g}^{2}\mathbf{I})$ and
the matrix $\mathcal{E}$ is diagonal with identical entries, i.e.,
$\mathcal{E}=\tau G^{-1}\mathbf{I}$, the free entropy function can be
reformulated as a scalar valued function, given by (14).

###### Proof 2

Please see the appendix -B.

Note that the MMSE under the isotropic channel case can be derived by
calculating the global maximum point of the free entropy function (14), and
the fixed point of the AMP state evolution can be predicted by the largest
$\tau$ that associated with a local maximum of (14). As we shall see in the
following Sec. IV-C, there exists a region of system parameters, i.e., pilot
length $L$, transmit power $P_{t}$, fraction of number of active user devices
$\rho$ and number of BS antennas $M$, where the MSE of AMP is blocked by a
local maximum point, so that the AMP algorithm is suboptimal. However, the
following proposition indicates that when considering an asymptotic MIMO
regime with the number of BS antennas go to infinity, the region where the AMP
algorithm is suboptimal vanishes.

###### Proposition 1

In an asymptotic massive MIMO regime, i.e., $M\to\infty$, by neglecting some
constants and the irrelevant factor, the free entropy function (14) is reduced
to

$\displaystyle\Phi(\tau)=$
$\displaystyle-\frac{\alpha\sigma_{w}^{2}}{\frac{1}{\alpha}\tau+\sigma_{w}^{2}}-\alpha\log(\sigma_{w}^{2}+\frac{1}{\alpha}\tau)$
$\displaystyle-\rho\sum\limits_{g}\log(1+\frac{\sigma_{g}^{2}}{\sigma_{w}^{2}+\frac{1}{\alpha}\tau}),$
(15)

which has and only has one local maximum point, and can be derived by
calculating the equation

$\tau=\rho\sum\limits_{g}\left(\sigma_{g}^{-2}+\left(\frac{1}{\alpha}\tau+\sigma_{w}^{2}\right)^{-1}\right)^{-1}.$
(16)

###### Proof 3

Please see the appendix -C.

Therefore, with an asymptotic number of BS antennas, since there is only one
local maximum point, the AMP algorithm always achieves Bayes-optimal
performance in the isotropic Rayleigh channel case, and the fixed point of AMP
state evolution can be derived by solving (16).

###  IV-B Prediction of AMP-based AUD and CE

(a) $P_{t}=18$dBm.

(b) $P_{t}=23$dBm.

(c) $P_{t}=33$dBm.

Figure 2: Free entropy as a function of MSE with $M=2$ BS antennas under
different settings of $\alpha$ and $P_{t}$. (The ``red cross'' denotes the
MMSE point, the ``red circle'' denotes the AMP-achievable MSE point, the ``red
square'' denotes the unachievable local maxima MSE point, and such marks will
also be used in the following figures.)

(a) $M=4$.

(b) $M=8$.

(c) $M\to\infty$.

Figure 3: Free entropy as a function of MSE with $P_{t}=33$dBm transmit power
under different settings of $\alpha$ and $M$.

Define $\tau^{\star}$ as the largest stationary point of (14) that associated
with the fixed point of the AMP state evolution. Together with (9), the LLR in
the isotropic channel case is

$\displaystyle\text{LLR}(\hat{\bm{s}_{g}})=M\log\frac{\sigma_{w}^{2}+\frac{1}{\alpha}}{\sigma_{g}^{2}+\sigma_{w}^{2}+\frac{1}{\alpha}\tau^{\star}}$
$\displaystyle+\hat{\bm{s}_{g}}^{H}\left((\sigma_{w}^{2}+\frac{1}{\alpha}\tau^{\star})^{-1}-(\sigma_{g}^{2}+\sigma_{w}^{2}+\frac{1}{\alpha}\tau^{\star})^{-1}\right)\hat{\bm{s}_{g}}.$
(17)

Obviously, the corresponding detection sufficient statistic is
$||\hat{\bm{s}_{g}}||^{2}$. As a consequence, we have

$\displaystyle P_{D}(M)$
$\displaystyle=\gamma\left(M,\left(\sigma_{g}^{2}+\sigma_{w}^{2}+\frac{1}{\alpha\tau^{\star}}\right)^{-1}l_{g}^{\prime}\right)\Gamma^{-1}(M),$
(18) $\displaystyle P_{F}(M)$
$\displaystyle=\gamma\left(M,\left(\sigma_{w}^{2}+\frac{1}{\alpha\tau^{\star}}\right)^{-1}l_{g}^{\prime}\right)\Gamma^{-1}(M).$
(19)

where $\gamma(M,\cdot)\Gamma(M)^{-1}$ is the CDF of the Chi-square
distribution with $2M$ degrees of freedom, and

$\displaystyle l_{g}^{\prime}\triangleq\sigma_{g}^{-2}$
$\displaystyle\left(\sigma_{w}^{2}+\frac{1}{\alpha\tau^{\star}}\right)\left(\sigma_{g}^{2}+\sigma_{w}^{2}+\frac{1}{\alpha\tau^{\star}}\right)$
$\displaystyle\left(l_{g}-M\log\frac{\sigma_{w}^{2}+\frac{1}{\alpha}\tau^{\star}}{\sigma_{g}^{2}+\sigma_{w}^{2}+\frac{1}{\alpha}\tau^{\star}}\right),$

with arbitrary threshold $l_{g}$. Further, in the asymptotic massive MIMO
regime, as $M$ go to infinity, we have
$\lim\nolimits_{M\to\infty}P_{F}(M)=P_{M}(M)=0$.

As for the CE performance, according to (12) and $\mathcal{E}=\tau
G^{-1}\mathbf{I}$, after performing the AMP algorithm, the CE error for an
active user in the group $g$ is

$\mathbb{E}\left[\left(\bm{s}_{g}-\hat{\bm{h}}_{g}\right)\left(\bm{s}_{g}-\hat{\bm{h}}_{g}\right)^{H}\right]=\epsilon\mathbf{I},$
(20)

where
$\epsilon=\left(\left(\sigma_{w}^{2}+\frac{1}{\alpha}\tau^{\star}\right)^{-1}+\sigma_{g}^{-2}\right)^{-1}$.

Therefore, in the isotropic Rayleigh channel, perfect AUD and the Bayes-
optimal channel estimation can be achieved as long as the number of antennas
is large enough. We also note the above results validate and extend the
analytical results in [6]. Besides, studying the free entropy function (14)
provides a phase transition diagram and an optimality analysis of the
isotropic channel case with a finite number of BS antennas, as shown in the
following section.

###  IV-C Verification in the Isotropic Channel Scenario

In this section, we provide numerical examples to verify our results based on
the replica method in the isotropic channel scenario. Specifically, we
consider that each user accesses the channel with a probability $\rho=0.1$ and
we assume the served user devices have been divided into $5$ user groups and
the distance $d_{g}$ between each user group is randomly distributed in the
regime $[0.1,1]$km. The path loss model of the wireless channel for each user
group $g$ is given as ${\rm PL}_{g}=-128.1-36.7\log_{10}(d_{g})$ in dB. Since,
we assume no power adaption, we denote the transmit power for each user as
$P_{t}$, and we have $\sigma_{g}^{2}=P_{t}\times{\rm PL}_{g}$. We assume the
bandwidth of the wireless channel are $1$MHz and the power spectral density of
the AWGN at the BS is $-169$dBm/Hz. The

####  IV-C1 Phase Transition of the Free Entropy Function

First, we examine the phase transition of the free entropy function and
demonstrating the optimal recovery MSE and the AMP-achievable MSE under our
scenario, showing that for different settings of the system parameters, i.e.,
$P_{t}$ and $\alpha$, the MSE performance can be divided into some performance
regions.

Specifically, we can notice that in the Fig. 2(c), there exists phase
transitions and the MSE performance can be divided into several regions. In
the first region, with $\alpha=0.525$, the free entropy function has only one
local maximum, indicating the optimal MSE performance, and the AMP-achievable
MSE coincide. In the second region, with $\alpha=0.550$, the second local
maximum point with a lower function value than the first local maximum of the
free entropy function appears, indicating the unachievable MSE point. In the
third region, with $\alpha=0.575$, the smaller local maximum point leads to a
larger value of the free entropy function, which indicates the MMSE under such
a choice of parameters. However, a larger local maximum point blocks the MSE
performance of AMP, since it always converges to the local maximum point
associated with the largest MSE. Hence, in this region, there exists a gap
between the AMP-achievable MSE and the MMSE. Finally, in the region with
$\alpha=0.600\sim 0.675$, one local maximum point disappears, and the AMP-
achievable MSE and the MMSE coincide again. We can intuitively infer that
there exists a hard threshold between the latter two regions, i.e, when the
pilot length exceeds the threshold, the local maximum point associated with
the largest MSE $\tau$ discards, and the AMP algorithm converges to the local
maximum associated with a lower MSE, leading that there is a phase transition
of the MSE performance of AMP.

We remark that in order to reduce the gap between the AMP-achievable MSE and
the MMSE in the case that the AMP algorithm is sub-optimal, one can consider
the idea of _seeding matrix_ introduced in [11] in the design of pilot matrix
$\bm{F}$, which is beyond the scope of this manuscript.

In addition, Figs. 2(a)-2(c) demonstrate the free entropy as a function of MSE
with $M=2$ BS antennas, and Figs. 2(a)-2(c) demonstrate the free entropy
function with $P_{t}=33$dBm transmit power. We can notice that the region
where AMP is suboptimal persists, but becomes smaller and eventually
disappears as the number of BS antennas increases and transmit power
decreases. Interestingly, we find that in the isotropic channel scenario,
$M=8$ antennas is enough to avoid phase transition in the case with the
transmit power is lower than $P_{t}=33$dBm.

####  IV-C2 Performance Analysis of AUD

Figure 4: Detection performance prediction of AMP with different settings of
pilot length.

We then provide numerical results to analyze the missed detection and false
alarm probabilities in the isotropic channel scenario. The number of user
devices is set as $N=10000$. Fig. 4 demonstrates the prediction of AMP-based
AUD by equation (19) and (18) versus pilot length, with $\tau^{\star}$
predicted by the free entropy function. The simulation curves are depicted by
the empirical AMP algorithm [26, 13]. Different from the standard AMP
algorithm, the empirical AMP algorithm adopts an empirical state evolution
function in the iterations, substituting the standard state evolution function
in the AMP algorithm, in order to avoid the complex expectation calculations.
The numerical results shows that our predicted performances are consistent
with that of the empirical AMP algorithm for most of the settings. The
empirical AMP and the proposed theoretical prediction have the similar
tendency and there exists phase transition phenomenons in both of them. We
notice that when phase transition occurs, the required length of pilot
sequences for our prediction is slightly shorter than that for the empirical
algorithm. Since the prediction error is about $0.04$ over the length of pilot
sequence, we think the theoretical results can provide a prediction for the
detection performance in the practice. We also note that our predicted
performance corresponding to the AMP achievable MSE provides a bound of the
performance of the empirical AMP algorithm. The performance loss of the
empirical AMP algorithm is because that the state evolution function that
fulfills the fixed point condition (8) of the free entropy function is
substituted by an empirical one.

$\displaystyle\Phi_{g}(\mathcal{E}_{g})=-\alpha{\rm
Tr}\left(\left(\bm{\Delta}+\frac{\mathcal{E}_{g}}{\alpha}\right)^{-1}\left(\frac{1}{\alpha}\rho\mathcal{C}_{g}+\bm{\Delta}\right)\right)-\alpha\log\left|\bm{\Delta}+\frac{\mathcal{E}_{g}}{\alpha}\right|+\int{\rm
d}\bm{s}_{g}[(1-\rho)\delta(\bm{s}_{g})+\rho Q_{g}(\bm{s}_{g})]\int{\rm
D}\bm{z}$ $\displaystyle\log\left(\int{\rm
d}\bm{x}_{g}[(1-\rho)\delta(\bm{x}_{g})+\rho
Q_{g}(\bm{x}_{g})]\right.\exp\left(-\bm{x}_{g}^{H}(\bm{\Delta}+\frac{\mathcal{E}_{g}}{\alpha})^{-1}\bm{x}_{g}+2\mathfrak{R}\left(\bm{x}_{g}^{H}[(\bm{\Delta}+\frac{\mathcal{E}_{g}}{\alpha})^{-1}\bm{s}_{g}+(\bm{\Delta}+\frac{\mathcal{E}_{g}}{\alpha})^{-\frac{1}{2}}\bm{z}]\right)\right).$
(21)

$\displaystyle{\Phi}_{g}(\bm{\Xi}_{g})=$
$\displaystyle-\alpha\sum\limits_{m=1}^{r_{g}}\left(\frac{\sigma_{w}^{2}}{\frac{1}{\alpha}\xi_{g,m}+\sigma_{w}^{2}}+\log\left(\frac{1}{\alpha}\xi_{g,m}+\sigma_{w}^{2}\right)\right)+\sum\limits_{m=1}^{r_{g}}\frac{(1-\rho)\lambda_{g,m}}{\lambda_{g,m}+\sigma_{w}^{2}+\frac{1}{\alpha}\xi_{g,m}}$
$\displaystyle+\int{\rm
D}\bm{z}\rho\log\left[(1-\rho)\prod_{m=1}^{r_{g}}\exp\left\\{-\frac{|z_{m}|^{2}\lambda_{g,m}}{\sigma_{w}^{2}+\frac{1}{\alpha}\xi_{g,m}}\right\\}+\rho\prod_{m=1}^{r_{g}}\frac{\frac{1}{\alpha}\xi_{g,m}+\sigma_{w}^{2}}{\frac{1}{\alpha}\xi_{g,m}+\sigma_{w}^{2}+\lambda_{g,m}}\right]$
$\displaystyle+\int{\rm
D}\bm{z}(1-\rho)\log\left[(1-\rho)\prod_{m=1}^{r_{g}}\exp\left\\{-\frac{|z_{m}|^{2}\lambda_{g,m}}{\lambda_{g,m}+\sigma_{w}^{2}+\frac{1}{\alpha}\xi_{g,m}}\right\\}+\rho\prod_{m=1}^{r_{g}}\left(\frac{\frac{1}{\alpha}\xi_{g,m}+\sigma_{w}^{2}}{\frac{1}{\alpha}\xi_{g,m}+\sigma_{w}^{2}+\lambda_{g,m}^{2}}\right)\right].$
(22)

####  IV-C3 Performance Analysis of CE

Figure 5: Channel estimation performance prediction of AMP with different
settings of pilot length.

The results of the CE error prediction is demonstrated in Fig. 5. We can
observe that the CE error prediction performance is consistent compared with
the detection error prediction. Our numerical results show that when the phase
transition occurs, the performances of the joint AUD and CE are highly
improved. The predicted performance corresponding to the AMP achievable MSE
reveals the minimum length of pilot sequences to make the AMP algorithm
accomplish the phase transition. We think this is a guide to how long pilot
sequences are needed when designing the system. For the prediction error that
may occur when using the empirical AMP algorithm in the practice, we can
appropriately increase the number of pilots to ensure that the phase
transition occurs.

##  V Spatially Correlated Channel

In this section, we concern another typical channel scenario in the wireless
communications literature, where the scattering is localized around the user
devices and the BS is elevated and thus has no scatterers in its near filed
[27], resulting in the spatially correlated channel. Assuming no line-of-sight
propagation, the channel of each user $k$ in the group $g$ is distributed as
$Q_{g}(\mathbf{h}_{g})=\mathcal{CN}(\bm{0},\mathcal{C}_{g})$, where
$\mathcal{C}_{g}=\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}$ with a rank
$r_{g}\ll M$. Considering the number of BS antennas is sufficiently large, the
$\bm{\Lambda}_{g}$ is approximated diagonal [17, 18].

We further suppose that different user groups are sufficiently well separated
in the angle of arrival (AoA) domain and the angular spread (AS) of each group
is sufficiently small. Accordingly, we assume that channels coefficients in
all the user groups are in different mutually orthogonal subspaces such that
$\bm{U}_{g}^{H}\bm{U}_{j}=\bm{0}$, for $j\neq g$. We note that although
directly achieving the mutually orthogonal subspaces is too restrictive, some
user scheduling strategies can be adopted to guarantee the user groups in
mutually orthogonal subspaces are served simultaneously [16]. In the following
of this section, we will see that the above spatially correlated channel
assumption reduces the expression of the free entropy function (7) and
provides some novel propositions for joint AUD and CE.

###  V-A Free Entropy Function

Before deriving the free entropy function in the spatially correlated channel
scenario, we recall that the matrix $\mathcal{E}$ is defined as
$\mathcal{E}G=\sum_{g}\mathcal{E}_{g}$, where $\mathcal{E}_{g}$ is the
corresponding recovery MSE matrix of each user group $g$. As a consequence, we
have the following lemma.

###### Lemma 1

Under the spatially correlated channel assumption, i.e.,
$\bm{U}_{g}^{H}\bm{U}_{j}=\bm{0}$, for $j\neq g$, all the stationary points of
the free entropy function (7) fulfill the equation:
$\mathcal{E}^{(\text{s})}=\sum_{g}\bm{U}_{g}\bm{\Xi}_{g}^{(\text{s})}\bm{U}_{g}^{H}$
for some symmetric positive semidefinite matrix $\bm{\Xi}_{g}^{(\text{s})}$.

###### Proof 4

Please see the appendix -D.

Since the statistical characteristics of our considered scenario are all
reflected in the stationary points of the free entropy function, we restrict
the form of the matrix $\mathcal{E}$ as
$\mathcal{E}=\sum_{g}\bm{U}_{g}\bm{\Xi}_{g}\bm{U}_{g}^{H}$. As a consequence,
the following theorem holds.

###### Theorem 3

Under the spatially correlated channel assumption, the free entropy function
(7) can be decoupled into a summation of $G$ independent free entropy
functions, i.e., $\Phi(\mathcal{E})=\sum_{g=0}^{G}\Phi_{g}(\mathcal{E}_{g})$,
where $\Phi_{0}(\mathcal{E}_{0})$ is a constant, and the specific expression
of each $\Phi_{g}(\mathcal{E}_{g})$ is formulated in the equation (21).

###### Proof 5

Please see the appendix -E.

Accordingly, the recovery MSE of each user group can be evaluated separately
via the equation (21), indicating that the transmitted signals from users
outside one group will not affect the reconstruction within such a group. In
the following, the performance analysis of joint AUD and CE in the spatially
correlated channel scenario will be provided based on the _Theorem 3_.

(a) $\alpha=0.11$.

(b) $\alpha=0.13$.

(c) $\alpha=0.15$.

Figure 6: Free entropy as a function of MSE with $P_{t}=18$dBm transmit power
under different settings of $\alpha$ in spatially correlated channel scenario.
(The ``red cross'' denotes the MMSE point, the ``red circle'' denotes the
predicted AMP-achievable MSE point, the ``red square'' denotes the
unachievable local maxima MSE point, and the ``blue triangle'' denotes the
simulated AMP MSE point.)

###  V-B Prediction of AMP-Based joint AUD and CE

According to the _Lemma 1_ , we define $\bm{\Xi}_{g}^{\star},\forall g$ as the
matrixs that associated with the fixed point $\mathcal{E}^{\star}$ of AMP
state evolution. We further restrict the structure of the matrix
$\bm{\Xi}_{g}^{\star}$ to be diagonal555Restricting the matrix
$\bm{\Xi}_{g}^{\star}$ into diagonal is reasonable. As proved in [13], the
state of state evolution function always maintains the same structure with
iterations. The initial state of the state evolution in the AMP framework is
always set to be $\mathcal{E}_{g}^{(0)}=\mathbb{E}(\bm{s}_{g}\bm{s}_{g}^{H})$,
and in our situation, we have
$\mathcal{E}_{g}^{(0)}=\rho\mathcal{C}_{g}=\rho\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}$.
Hence, $\bm{\Xi}_{g}^{\star}$ will also be diagonal since the matrix
$\bm{\Lambda}_{g}$ is approximately diagonal in our scenario with a large
number of antennas., with the $m$th diagonal element denoted as
$\xi_{g,m}^{\star}$. By combining _Theorem 3_ with _Lemma 1_ , the matrix
$\bm{\Xi}_{g}^{\star}$ can be obtained by seeking the local maximum point of
the equation (22), where the quantities $\xi_{g,m}$, $\lambda_{g,m}$ are
defined as the $m$th element of the diagonal of $\bm{\Xi}_{g}$ and
$\bm{\Lambda}_{g}$.

Then, we consider the performance analysis of joint AUD and CE. For AUD,
defining $\underline{\bm{s}_{g}}\triangleq\bm{U}_{g}^{H}\hat{\bm{s}}_{g}$ and
together with (9), the detection sufficient statistic in such a scenario can
be formulated as

$\text{T}(\underline{\bm{s}_{g}})=\sum\limits_{m=1}^{r_{g}}\frac{\lambda_{g,m}|\underline{s_{g,m}}|^{2}}{\left(\sigma_{w}^{2}+\frac{1}{\alpha}\xi_{g,m}^{\star}\right)\left(\sigma_{w}^{2}+\frac{1}{\alpha}\xi_{g,m}^{\star}+\lambda_{g,m}\right)}.$
(23)

and the involved conditional probabilities of $\underline{\bm{s}_{g}}$ are
given as

$\displaystyle p(\underline{\bm{s}_{g}}|a_{g}=1)$
$\displaystyle=\mathcal{CN}(\underline{\bm{s}_{g}};\bm{0},\bm{\Lambda}_{g}+\bm{\Delta}+\frac{1}{\alpha}\bm{\Xi}_{g}^{\star}),$
$\displaystyle p(\underline{\bm{s}_{g}}|a_{g}=0)$
$\displaystyle=\mathcal{CN}(\underline{\bm{s}_{g}};\bm{0},\bm{\Delta}+\frac{1}{\alpha}\bm{\Xi}_{g}^{\star}).$
(24)

As a result, we can then derive the specific functional forms of the detection
performance metrics in the following proposition.

###### Proposition 2

The missed detection probability with the definition $P_{D}={\rm
Pr}\\{\text{T}(\underline{\bm{s}_{g}})>l_{g}|a_{g}=1\\}$ and the false alarm
probability with the definition $P_{F}={\rm
Pr}\\{\text{T}(\underline{\bm{s}_{g}})>l_{g}|a_{g}=0\\}$ based on the
conditional probabilities (24) and the detection sufficient statistic (23) are
given as

$\displaystyle P_{D}=\sum\limits_{m=1}^{r_{g}}\prod_{j\neq
m}\frac{\omega_{g,j}}{\omega_{g,j}-\omega_{g,m}}\exp\\{-\omega_{g,m}l_{g}\\},$
$\displaystyle P_{F}=\sum\limits_{m=1}^{r_{g}}\prod_{j\neq
m}\frac{\tilde{\omega}_{g,j}}{{\tilde{\omega}}_{g,j}-{\tilde{\omega}}_{g,m}}\exp\\{-{\tilde{\omega}}_{g,m}l_{g}\\},$
(25)

where

$\omega_{g,m}\triangleq(\sigma_{w}^{2}+\frac{1}{\alpha}\xi_{g,m}^{\star})\lambda_{g,m}^{-1},$

and

$\tilde{\omega}_{g,m}\triangleq(\sigma_{w}^{2}+\frac{1}{\alpha}\xi_{g,m}^{\star}+\lambda_{g,m})\lambda_{g,m}^{-1}.$

###### Proof 6

Please see the appendix -F.

As for the analysis for CE, according to (13), the CE error of the user group
$g$ associated with the $m$th eigenvalue $\lambda_{g,m}$ can be formulated as
$(\sigma_{w}^{2}+\lambda_{g,m}+\frac{1}{\alpha}\xi_{g,m}^{\star})^{-1}$.

After that, we can see that the performance of joint AUD and CE can be
predicted by the local maximum point of the free entropy function (22). The
verification of such a result, a phase transition diagram and an optimality
analysis of joint AUD and CE in the spatially correlated channel scenario will
be shown in the following Sec. V-D.

###  V-C Discussion

Before proceeding, we provide some discussions about our theoretical results.
Multiplying $\bm{U}_{g}^{*}$ for each group by the right of the received
signal $\bm{Y}$ in the original model (1), we can obtain $G$ a sub-model
expressed as

$\underline{\bm{Y}}_{g}=\underline{\bm{F}}_{g}\underline{\bm{X}}_{g}+\underline{\bm{W}}_{g},$
(26)

where $\underline{\bm{Y}}_{g}\triangleq\bm{Y}\bm{U}_{g}^{*}$,
$\underline{\bm{X}}_{g}\triangleq\bm{X}\bm{U}_{g}^{*}$,
$\underline{\bm{W}}_{g}\triangleq\bm{W}\bm{U}_{g}^{*}$ and
$\underline{\bm{F}}_{g}$ is the corresponding pilot matrix for the user device
for the group $g$. We note that the free entropy function of (26) reflecting
its statistical performances corresponds to (22), so that separately using
model (26) for each group is equivalent to jointly processing all the groups
with (1), if the mutually orthogonal subspaces condition holds. Such a per-
group processing (PGP) idea has been successfully adopted in the case of
downlink transmission [16] and our theoretical results provide the theoretical
foundations for the PGP strategy in the uplink joint AUD and CE problem.

Complementally, when the mutually orthogonal condition does not perfectly
hold, by carefully designing the processing matrix for each group, the PGP
strategy can also be applied by permitting the affordable inter-group
inference. According to the approach in [16], when the BS equipped with a
uniform linear array, the inter-group inference can be absorbed into the
additive noise matrix with independent elements. Hence, by slightly modifying
the free entropy function (22), our framework can also be adopted to evaluate
performance of joint AUD and CE, where the mutually orthogonal condition does
not perfectly hold.

###  V-D Verification in the Spatially Correlated Channel Scenario

In this section, we provide numerical results to verify the analysis based on
the replica method in spatially correlated channel scenario. We consider the
channel coefficients of all user groups are in orthogonal subspaces. Without
loss of generality, we investigate the performance analysis results for one
certain user group. We consider the typical spatially correlated channel model
in [27]. For visualization, we assume the number of BS antennas is $M=64$. The
path loss model of the wireless channel, the distribution of the distance
between the BS station and the user devices, the bandwidth the power spectral
density and the total number of user devices are considered same with that in
the settings of the isotropic channel case.

####  V-D1 Phase Transition of the Free Entropy Function

We first consider the center angle of user angle domain is set as $45^{\circ}$
distributed as a normal distribution with a $1^{\circ}$ AS, where most of the
power of channel coefficients is concentrated on a two-dimensional subspace.
Figs. 6(a)-6(c) shows the different MSE performance regions of the free
entropy function (22) in different settings. We can notice that the
performance regions in the spatially correlated channel case is similar with
that in the isotropic channel case. Specifically, Fig. 6(b) demonstrates the
performance gap between the AMP-achievable MSE and the MMSE. The phase
transitions appear in both our prediction and the simulated AMP algorithm.
Differently from the isotropic channel scenario, it is observed that although
there exists a large number of antennas in the BS, the phase transition
phenomenon still exists. Fig. 6 also shows that our analytical MSE results
predicted in the large system limit and AMP recovery MSE in practical
settings, i.e., $T=220,260,300$ and $K=2000$ are consistent. Such results are
evidences that our analytical results are valid when the length of pilot
sequences and the number of user devices are large but finite.

(a) $P_{t}=13$dBm.

(b) $P_{t}=23$dBm.

(c) $P_{t}=33$dBm.

Figure 7: NMSE performance of AMP versus pilot length in different settings of
transmit power and AS.

Then, we demonstrate the case where the AS is not very small, which is more in
line with the practical scenario. The phase transition phenomenons in all the
considering settings can be observed in Fig. 7. With the transmit power
increasing, the phase transition phenomenon is becoming more and more obvious.
In addition, with the AS increases, the number of the efficient dimensions in
the subspace of the user channel becomes larger and the total energy will be
allocated to these effective dimensions, leading that the phase transition
phenomenon is weakened. Such an observation matches our previous results that
the spatially correlated channel is more likely to promote a phase transition
compared with the isotropic channel. Another observation is that smaller AS is
beneficial to NMSE performance when pilot resources are abundant, while the
pilot resources are scarce, the opposite is true.

####  V-D2 Performance Analysis of AUD

Figure 8: Detection performance prediction.

Fig. 8 depicts the prediction of AMP-based AUD as a function of pilot length
with $1^{\circ}$ AS. It is observed that the phase transition of the detection
error probabilities appears in both $P_{t}=13$dBm and $P_{t}=18$dBm settings,
which is different from the results in the isotropic channel case. This is
because the energy of a large number of antennas is concentrated in a small
amount of subspace dimensions, increasing the signal-to-noise ratio in
considered dimensions.

In summary, we note that the spatially correlated channel promotes the
appearance of phase transition compared with the isotropic channel. This is
because the energy of a large number of antennas is concentrated in a small
amount of subspace dimensions. This reduces the effective dimensions of the BS
antennas, while increasing the power on each dimension. Despite that, the
pilot length required to step over the phase transition is highly reduced in
the spatially correlated channel scenario, i.e., about $\alpha=0.14$ for the
spatially correlated channel case and about $\alpha=0.625$ for the isotropic
channel case with $5$ user groups.

##  VI Conclusion

In this paper, we have provided an analysis of joint AUD and CE problem under
massive connectivity based on the replica method. Particularly, we have
established a novel theoretical framework under a general channel model that
reduces to multiple typical MIMO channel models. Based on the general
framework, we have analyzed two typical scenarios in massive connectivity,
i.e., the isotropic channel case and the spatially correlated channel case. We
have provided the analysis of the Bayes-optimility, phase transition, and the
predictions of the performance of joint AUD and CE in both the two cases. In
addition, we have shown that the spatially correlated channel is more likely
to promote a phase transition compared with the isotropic channel, due to the
small number of subspace dimensions and corresponding antenna gains. However,
thanks to the spatially correlated channel, where all the user group can be
perfectly partitioned, the pilot length required to step over the phase
transition can be significantly reduced.

Some future directions of research are also implied by this paper. 1) Our
theoretical framework was established in the Bayes-optimal condition, if the
condition is not met, what changes will exist in the analytical results. 2)
Exploiting the potential user grouping information may bring some advantages,
and how to design an algorithm to achieve joint user grouping, AUD and CE. 3)
How to design the processing matrix in the PGP-based strategy for joint AUD
and CE.

###  -A Proof of Theorem 1

According to the expression of partition function $Z$ in (4), we notice that
$Z^{n}$ can be written as

$\displaystyle Z^{n}=\int$ $\displaystyle\prod_{a}{\rm
d}\bm{X}^{(a)}\prod_{a,g,k}\left((1-\rho)\delta(\bm{x}_{gk}^{(a)})+\rho
Q_{g}(\bm{x}_{gk}^{(a)})\right)$
$\displaystyle\times\prod_{t}\frac{1}{\pi^{nM}|\bm{\Delta}|^{n}}e^{-\sum\limits_{a}(\bm{y}_{t}-\bm{f}_{t}\bm{X}^{(a)})\bm{\Delta}^{-1}(\bm{y}_{t}-\bm{f}_{t}\bm{X}^{(a)})^{H}},$

$\displaystyle 1$ $\displaystyle=\int{\rm
d}\bm{\Theta}\exp\left(\sum\limits_{g}{\rm
Tr}\left(\sum\limits_{a}\hat{\mathcal{Q}}_{g}^{(a)H}\left(\frac{K}{2}\mathcal{Q}_{g}^{(a)}-\frac{1}{2}\sum\limits_{k}\mathbf{x}_{kg}^{(a)}\mathbf{x}_{kg}^{(a)^{H}}\right)+\sum\limits_{a}\hat{\mathcal{Q}}_{g}^{(a)T}\left(\frac{K}{2}\mathcal{Q}_{g}^{(a)*}-\frac{1}{2}\sum\limits_{k}\mathbf{x}_{kg}^{(a)*}\mathbf{x}_{kg}^{(a)^{T}}\right)\right.\right.$
$\displaystyle\left.-\sum\limits_{a\neq
b}\hat{\mathcal{T}}_{g}^{(ab)H}\left(\frac{K}{2}\mathcal{T}_{g}^{(ab)}-\frac{1}{2}\sum\limits_{k}\mathbf{x}_{kg}^{(a)}\mathbf{x}_{kg}^{(b)^{H}}\right)\right.-\left.\sum\limits_{a\neq
b}\hat{\mathcal{T}}_{g}^{(ab)T}\left(\frac{K}{2}\mathcal{T}_{g}^{(ab)*}-\frac{1}{2}\sum\limits_{k}\mathbf{x}_{kg}^{(a)*}\mathbf{x}_{kg}^{(b)^{T}}\right)\right.$
$\displaystyle\left.\left.-\sum\limits_{a}\hat{\mathcal{M}}_{g}^{(a)H}\left(K\mathcal{M}_{g}^{(a)}-\sum\limits_{k}\mathbf{x}_{kg}^{(a)}\mathbf{s}_{kg}^{H}\right)-\sum\limits_{a}\hat{\mathcal{M}}_{g}^{(a)T}\left(K\mathcal{M}_{g}^{(a)*}-\sum\limits_{k}\mathbf{x}_{kg}^{(a)*}\mathbf{s}_{kg}^{T}\right)\right)\right),$
(30)

$\displaystyle\mathbb{E}Z^{n}=\int{\rm
d}\bm{\Theta}\exp\bigg{(}Kn\sum\limits_{g}{\rm
Tr}\left(\frac{1}{2}\left(\hat{\mathcal{Q}}_{g}^{H}\mathcal{Q}_{g}+\hat{\mathcal{Q}}_{g}^{T}\mathcal{Q}_{g}^{*}\right)-\frac{1}{2}(n-1)\left(\hat{\mathcal{T}}_{g}^{H}\mathcal{T}_{g}+\hat{\mathcal{T}}_{g}^{T}\mathcal{T}_{g}^{*}\right)\right.-\left.\left(\hat{\mathcal{M}}_{g}^{H}\mathcal{M}_{g}+\hat{\mathcal{M}}_{g}^{T}\mathcal{M}_{g}^{*}\right)\right)\bigg{)}$
$\displaystyle\times\prod_{t}\frac{\beta_{t}}{\pi^{Mn}|\bm{\Delta}|^{n}}\prod_{g}\bigg{(}\int{\rm
d}\bm{s}_{g}((1-\rho)\delta(\bm{s}_{g})+\rho Q_{g}(\bm{s}_{g}))\int{\rm
D}\bm{z}\bigg{(}\int\left({\rm d}\bm{x}_{g}((1-\rho)\delta(\bm{x}_{g})+\rho
Q_{g}(\bm{x}_{g}))\right)$
$\displaystyle\times\exp\left(-\frac{1}{2}\bm{x}_{g}^{H}(\hat{\mathcal{Q}}_{g}^{H}+\hat{\mathcal{Q}}_{g}+\hat{\mathcal{T}}_{g}^{H}+\hat{\mathcal{T}}_{g})\bm{x}_{g}\right.+\bm{x}_{g}^{H}\big{(}\big{(}\frac{\hat{\mathcal{T}}_{g}^{H}+\hat{\mathcal{T}}_{g}}{2}\big{)}^{\frac{1}{2}}\bm{z}+\hat{\mathcal{M}}_{g}\bm{s}_{g}\big{)}+\big{(}\bm{z}^{H}\big{(}\frac{\hat{\mathcal{T}}_{g}^{H}+\hat{\mathcal{T}}_{g}}{2}\big{)}^{\frac{1}{2}}+\bm{s}_{g}^{H}\hat{\mathcal{M}}_{g}^{H}\big{)}\bm{x}_{g}\bigg{)}\bigg{)}^{n}\bigg{)}^{K}.$
(32)

$\displaystyle\beta_{t}=$
$\displaystyle\left|\mathbf{I}_{Mn}+\sigma_{w}^{-2}\mathcal{G}\right|^{-1}$
$\displaystyle=$
$\displaystyle\left|n\left(\bm{\Delta}+\frac{K}{T}\sum\limits_{g}\left(\mathcal{Q}_{g}-\mathcal{T}_{g}\right)\right)^{-1}\left(\frac{K}{T}\left(\sum\limits_{g}\rho\mathcal{C}_{g}-\mathcal{M}_{g}-\mathcal{M}_{g}^{H}+\mathcal{T}_{g}\right)+\bm{\Delta}\right)+\mathbf{I}_{M}\right|^{-1}\times\left|\mathbf{I}_{M}+\sigma_{w}^{-2}\frac{K}{T}\sum\limits_{g}(\mathcal{Q}_{g}-\mathcal{T}_{g})\right|^{-n}$
$\displaystyle\approx$ $\displaystyle\exp\bigg{(}-n{\rm
Tr}\bigg{(}\bigg{(}\bm{\Delta}+\frac{K}{T}\sum\limits_{g}\bigg{(}\mathcal{Q}_{g}-\mathcal{T}_{g}\bigg{)}\bigg{)}^{-1}\bigg{(}\frac{K}{T}\bigg{(}\sum\limits_{g}\rho\mathcal{C}_{g}-\mathcal{M}_{g}-\mathcal{M}_{g}^{H}+\mathcal{T}_{g}\bigg{)}+\bm{\Delta}\bigg{)}\bigg{)}$
$\displaystyle\quad\quad\quad-n\log\left|\mathbf{I}_{M}+\sigma_{w}^{-2}\frac{K}{T}\sum\limits_{g}(\mathcal{Q}_{g}-\mathcal{T}_{g})\right|\bigg{)}.$
(33)

$\displaystyle\tilde{\Phi}(\bm{\mathcal{Q}},\bm{\mathcal{T}},\bm{\mathcal{M}},\bm{\hat{\mathcal{Q}}},\bm{\hat{\mathcal{T}}},\bm{\hat{\mathcal{M}}})=\sum\limits_{g}{\rm
Tr}\left(\frac{1}{2}\left(\hat{\mathcal{Q}}_{g}^{H}\mathcal{Q}_{g}+\hat{\mathcal{Q}}_{g}^{T}\mathcal{Q}_{g}^{*}\right)+\frac{1}{2}\left(\hat{\mathcal{T}}_{g}^{H}\mathcal{T}_{g}+\hat{\mathcal{T}}_{g}^{T}\mathcal{T}_{g}^{*}\right)-\left(\hat{\mathcal{M}}_{g}^{H}\mathcal{M}_{g}+\hat{\mathcal{M}}_{g}^{T}\mathcal{M}_{g}^{*}\right)\right)$
$\displaystyle-\alpha{\rm
Tr}\bigg{(}\bigg{(}\bm{\Delta}+\frac{1}{\alpha}\sum\limits_{g}\bigg{(}\mathcal{Q}_{g}-\mathcal{T}_{g}\bigg{)}\bigg{)}^{-1}\bigg{(}\frac{1}{\alpha}\bigg{(}\sum\limits_{g}\rho\mathcal{C}_{g}-\mathcal{M}_{g}-\mathcal{M}_{g}^{H}+\mathcal{T}_{g}\bigg{)}+\bm{\Delta}\bigg{)}\bigg{)}-\alpha\log\left|\bm{\Delta}+\frac{1}{\alpha}\sum\limits_{g}(\mathcal{Q}_{g}-\mathcal{T}_{g})\right|$
$\displaystyle+\sum\limits_{g}\bigg{(}\int{\rm
d}\bm{s}_{g}((1-\rho)\delta(\bm{s}_{g})+\rho Q_{g}(\bm{s}_{g}))\int{\rm
D}\bm{z}\log\bigg{(}\int\left({\rm
d}\bm{x}_{g}((1-\rho)\delta(\bm{x}_{g})+\rho
Q_{g}(\bm{x}_{g}))\right)\exp\left(-\frac{1}{2}\bm{x}_{g}^{H}(\hat{\mathcal{Q}}_{g}^{H}+\hat{\mathcal{Q}}_{g}\right.$
$\displaystyle+\hat{\mathcal{T}}_{g}^{H}+\hat{\mathcal{T}}_{g})\bm{x}_{g}+\bm{x}_{g}^{H}\big{(}\big{(}\frac{\hat{\mathcal{T}}_{g}^{H}+\hat{\mathcal{T}}_{g}}{2}\big{)}^{\frac{1}{2}}\bm{z}+\hat{\mathcal{M}}_{g}\bm{s}_{g}\big{)}+\big{(}\bm{z}^{H}\big{(}\frac{\hat{\mathcal{T}}_{g}^{H}+\hat{\mathcal{T}}_{g}}{2}\big{)}^{\frac{1}{2}}+\bm{s}_{g}^{H}\hat{\mathcal{M}}_{g}^{H}\big{)}\bm{x}_{g}\bigg{)}\bigg{)}\bigg{)}-\alpha\log\pi^{M}.$
(35)

As a consequence, by considering the equation (1), the average free entropy
can be derived as

$\displaystyle\mathbb{E}Z^{n}=\int\prod_{a}{\rm
d}\bm{X}^{(a)}\prod_{a,g,k}\left((1-\rho)\delta(\bm{x}_{gk}^{(a)})+\rho
Q_{g}(\bm{x}_{gk}^{(a)})\right)$ $\displaystyle\times\int\prod_{g,k}{\rm
d}\bm{s}_{gk}\left((1-\rho)\delta(\bm{s}_{gk})+\rho
Q_{g}(\bm{s}_{gk})\right)\prod_{t}\frac{1}{\pi^{nM}|\bm{\Delta}|^{n}}\beta_{t}.$
(27)

The expectation is taken over $\bm{F}$, $\bm{S}$, and $\bm{W}$, and the
quantity $\beta_{t}$ is defined as

$\displaystyle\beta_{t}\triangleq\mathbb{E}_{\bm{F},\bm{W}}\big{(}\exp(-\sigma_{w}^{-2}\sum\limits_{a}||\bm{v}_{t}^{(a)}||^{2})\big{)},$

where we define

$\bm{v}_{t}^{(a)}\triangleq\sum_{g,k}f_{tgk}\bm{s}_{gk}-\sum_{g,k}f_{tgk}\bm{x}_{gk}^{(a)}+\bm{w}_{t},a=1,\dots,n.$

We note that $f_{tgk}$ is the corresponding element in the vector
$\bm{f}_{t}$, and the vector $\bm{w}_{t}$ is the corresponding noise vector.
Since we consider a Gaussian random pilot matrix with zero mean,
$\bm{v}_{t}^{(a)}$ is still a Gaussian vector with zero mean. The covariance
matrices are calculated as follow.

$\begin{split}\mathbb{E}_{\bm{F},\bm{W}}(\bm{v}_{t}^{(a)}\bm{v}_{t}^{(a)H})=&\frac{1}{T}\sum_{g,k}(\bm{s}_{kg}-\bm{x}_{kg}^{(a)})(\bm{s}_{kg}-\bm{x}_{kg}^{(a)})^{H}+\bm{\Delta},\\\
\mathbb{E}_{\bm{F},\bm{W}}(\bm{v}_{t}^{(a)}\bm{v}_{t}^{(b)H})=&\frac{1}{T}\sum_{g,k}(\bm{s}_{kg}-\bm{x}_{kg}^{(a)})(\bm{s}_{kg}-\bm{x}_{kg}^{(b)})^{H}+\bm{\Delta}.\end{split}$

We further define a vector
$\tilde{\bm{v}}_{t}\triangleq[\bm{v}_{t}^{(1)},\dots,\bm{v}_{t}^{(n)}]^{T}$,
which is a stitched vector of $\bm{v}_{t}^{(a)}$. It is obvious that the
vector $\tilde{\bm{v}}_{t}$ is a Gaussian vector with zero mean, and we define
$\mathcal{G}$ as its covariance matrix. Then, the term $\beta_{t}$ can be
rewritten as

$\displaystyle\mathbb{E}_{\bm{F},\bm{W}}\big{(}\exp(-\sigma_{w}^{-2}\sum\limits_{a}||\bm{v}_{t}^{(a)}||^{2})\big{)}=|\mathbf{I}_{Mn}+\sigma_{w}^{-2}\mathcal{G}|^{-1}.$

We notice that the matrix $\mathcal{G}$ can be separated into $n\times n$
blocks of size $M\times M$, and the block in the $a$th row and $b$th column
satisfies
$\mathcal{G}_{ab}=\mathbb{E}_{\bm{F},\bm{W}}(\bm{v}_{t}^{(a)}\bm{v}_{t}^{(b)H})$.
Giving the overlaps

$\displaystyle\rho\mathcal{C}_{g}$
$\displaystyle=\frac{1}{K}\sum\limits_{k=1}^{K}\mathbf{s}_{kg}\mathbf{s}_{kg}^{H},\
\mathcal{M}_{g}^{(a)}=\frac{1}{K}\sum\limits_{k=1}^{K}\mathbf{x}_{kg}^{(a)}\mathbf{s}_{kg}^{H},$
$\displaystyle\mathcal{Q}_{g}^{(a)}$
$\displaystyle=\frac{1}{K}\sum\limits_{k=1}^{K}\mathbf{x}_{kg}^{(a)}\mathbf{x}_{kg}^{(a)^{H}},\
\mathcal{T}_{g}^{(ab)}=\frac{1}{K}\sum\limits_{k=1}^{K}\mathbf{x}_{kg}^{(a)}\mathbf{x}_{kg}^{(b)H},$
(28)

we have

$\displaystyle\mathcal{G}_{aa}$
$\displaystyle=\frac{K}{T}\left(\sum\limits_{g=1}^{G}\rho\mathcal{C}_{g}-\mathcal{M}_{g}^{(a)}-\mathcal{M}_{g}^{(a)H}+\mathcal{Q}_{g}^{(a)}\right)+\bm{\Delta},$
$\displaystyle\mathcal{G}_{ab}$
$\displaystyle=\frac{K}{T}\left(\sum\limits_{g=1}^{G}\rho\mathcal{C}_{g}-\mathcal{M}_{g}^{(a)}-\mathcal{M}_{g}^{(b)H}+\mathcal{T}_{g}^{(ab)}\right)+\bm{\Delta}.$
(29)

We note that the term $\beta_{t}$ is completely determined by the overlaps in
(28). For further simplifying the equation (27), we introduce an identity seen
in (30), which is similarly provided in [11, 28]. For concise, we use the
notation

$\displaystyle{\rm d}$ $\displaystyle\bm{\Theta}$ $\displaystyle=$
$\displaystyle\prod_{g,a}{\rm d}\hat{\mathcal{Q}}_{g}^{(a)}{\rm
d}\hat{\mathcal{Q}}_{g}^{(a)*}{\rm d}\hat{\mathcal{M}}_{g}^{(a)}{\rm
d}\hat{\mathcal{M}}_{g}^{(a)*}{\rm d}\mathcal{Q}_{g}^{(a)}{\rm
d}\mathcal{Q}_{g}^{(a)*}{\rm d}\mathcal{M}_{g}^{(a)}{\rm
d}\mathcal{M}_{g}^{(a)*}$ $\displaystyle\times\prod_{g,a,b}{\rm
d}\hat{\mathcal{T}}_{g}^{(ab)}{\rm d}\hat{\mathcal{T}}_{g}^{(ab)*}{\rm
d}\mathcal{T}_{g}^{(ab)}{\rm d}\mathcal{T}_{g}^{(ab)*}.$

Note that we have introduced the conjugated parameters
$\\{\hat{\mathcal{Q}}_{g}^{(a)},\hat{\mathcal{M}}_{g}^{(a)},\hat{\mathcal{T}}_{g}^{(ab)}\\}$
to enforce the consistency conditions (28). We then consider the replica-
symmetry assumption, which is a general consideration in the replica analysis.
Under this, we have the following equations.

$\displaystyle\hat{\mathcal{M}}_{g}^{(a)}=\hat{\mathcal{M}}_{g},\quad\hat{\mathcal{Q}}_{g}^{(a)}=\hat{\mathcal{Q}}_{g},\quad\hat{\mathcal{T}}_{g}^{(ab)}=\hat{\mathcal{T}}_{g},$
$\displaystyle\mathcal{M}_{g}^{(a)}=\mathcal{M}_{g},\quad\mathcal{Q}_{g}^{(a)}=\mathcal{Q}_{g},\quad\mathcal{T}_{g}^{(ab)}=\mathcal{T}_{g}.$
(31)

Further combining the Hubbard-Stratonovich transform [14]

$\begin{split}&e^{\frac{1}{2}\sum\limits_{a\neq
b}\bm{x}_{g}^{(a)H}(\hat{\mathcal{T}}_{g}+\hat{\mathcal{T}}_{g}^{H})\bm{x}_{g}^{(b)}+\frac{1}{2}\sum\limits_{a}\bm{x}_{g}^{(a)H}(\hat{\mathcal{T}}_{g}+\hat{\mathcal{T}}_{g}^{H})\bm{x}_{g}^{(a)}}\\\
=&\int{\rm
D}\bm{z}e^{\sum\limits_{a}\bm{x}_{g}^{(a)H}(\frac{\hat{\mathcal{T}}_{g}+\hat{\mathcal{T}}_{g}^{H}}{2})^{1/2}\bm{z}+\bm{z}^{H}\sum\limits_{a}(\frac{\hat{\mathcal{T}}_{g}+\hat{\mathcal{T}}_{g}^{H}}{2})^{1/2}\bm{x}_{g}^{(a)}},\end{split}$

the equation (27) can be re-formulated as (32).

We now turn to consider the $n\to 0$ limit. Recall the expression of
$\beta_{t}$ in (32). By combining the equations in (29), the matrix
$\mathcal{G}$ can be rewritten as

$\displaystyle\mathcal{G}=\bm{\amalg}_{n}$
$\displaystyle\otimes\left(\frac{K}{T}(\sum_{g}\rho\mathcal{C}_{g}-\mathcal{M}_{g}-\mathcal{M}_{g}^{H}+\mathcal{T}_{g})+\bm{\Delta}\right)$
$\displaystyle+\mathbf{I}_{n}\otimes\left(\frac{K}{T}\sum_{g}(\mathcal{Q}_{g}-\mathcal{T}_{g})\right),$

where $\bm{\amalg}_{n}$ stands for the $n\times n$ matrix with elements all
equal to one. The eigenvalue set of $\mathcal{G}$ consists of the $M$
eigenvalues of matrix
$\frac{K}{T}(\sum_{g}\rho\mathcal{C}_{g}-\mathcal{M}_{g}-\mathcal{M}_{g}^{H}+\mathcal{T}_{g})+\bm{\Delta}$,
and $(n-1)$ groups of $M$ eigenvalues with each group consisting the
eigenvalues of $\frac{K}{T}\sum_{g}(\mathcal{Q}_{g}-\mathcal{T}_{g})$. As a
consequence, $\beta_{t}$ can be rewritten and further approximated in (33).

Combining the final expression in (33) and with the approximation

$\int{\rm D}\bm{z}f(\bm{z})^{n}=1+n\int{\rm D}\bm{z}\log
f(\bm{z})\approx\exp\left(n\int{\rm D}\bm{z}\log f(\bm{z})\right)$

in the $n\to 0$ limit, we then obtain

$\mathbb{E}Z^{n}=\int{\rm
d}\bm{\Theta}\exp\left(Kn{\tilde{\Phi}}(\bm{\mathcal{Q}},\bm{\mathcal{T}},\bm{\mathcal{M}},\bm{\hat{\mathcal{Q}}},\bm{\hat{\mathcal{T}}},\bm{\hat{\mathcal{M}}})\right),$
(34)

where we have denoted
$\bm{\mathcal{Q}}\triangleq\\{\\{\mathcal{Q}_{g}\\}_{g=1}^{G},\\{\mathcal{Q}^{*}_{g}\\}_{g=1}^{G}\\}$,
$\bm{\mathcal{T}}\triangleq\\{\\{\mathcal{T}_{g}\\}_{g=1}^{G},\\{\mathcal{T}^{*}_{g}\\}_{g=1}^{G}\\}$,
$\bm{\mathcal{M}}\triangleq\\{\\{\mathcal{M}_{g}\\}_{g=1}^{G},\\{\mathcal{M}^{*}_{g}\\}_{g=1}^{G}\\}$,
$\bm{\hat{\mathcal{Q}}}\triangleq\\{\\{\hat{\mathcal{Q}}_{g}\\}_{g=1}^{G},\\{\hat{\mathcal{Q}}^{*}_{g}\\}_{g=1}^{G}\\}$,
$\bm{\hat{\mathcal{T}}}\triangleq\\{\\{\hat{\mathcal{T}}_{g}\\}_{g=1}^{G},\\{\hat{\mathcal{T}}^{*}_{g}\\}_{g=1}^{G}\\}$,
$\bm{\hat{\mathcal{M}}}\triangleq\\{\\{\hat{\mathcal{M}}_{g}\\}_{g=1}^{G},\\{\hat{\mathcal{M}}^{*}_{g}\\}_{g=1}^{G}\\}$
as the collections of the involved variables, and the function $\tilde{\Phi}$
is specified in the equation (35).

The saddle point method [28] is performed by taking the extremum of the (35)
with respect to the free parameters. Accordingly, we have

$\Phi={\rm
extr}\left({\tilde{\Phi}}(\bm{\mathcal{Q}},\bm{\mathcal{T}},\bm{\mathcal{M}},\bm{\hat{\mathcal{Q}}},\bm{\hat{\mathcal{T}}},\bm{\hat{\mathcal{M}}})\right),$

where ${\rm extr}(\cdot)$ returns the extremum of the involved function with
respect to its arguments. Setting derivatives of $\tilde{\Phi}$ with respect
to $\mathcal{M}_{g}$, $\mathcal{T}_{g}$, $\mathcal{Q}_{g}-\mathcal{T}_{g}$,
$\forall g$ and all their conjugations into zero, then we can obtain the self-
consistent equations as follow.

$\displaystyle\forall g,\hat{\mathcal{M}}_{g}^{H}=$
$\displaystyle\alpha\mathcal{W}_{g},\quad\hat{\mathcal{Q}}_{g}^{H}+\hat{\mathcal{T}}_{g}^{H}=\alpha\mathcal{W}_{g}.$
$\displaystyle\forall g,\hat{\mathcal{T}}_{g}^{H}=$
$\displaystyle\alpha\mathcal{W}_{g}(\frac{T}{K}\bm{\Delta}+\sum\limits_{g}(\mathcal{\rho}\mathcal{C}_{g}-\mathcal{M}_{g}-\mathcal{M}_{g}^{H}+\mathcal{T}_{g}))\mathcal{W}_{g},$

where we denote

$\mathcal{W}_{g}\triangleq(\frac{T}{K}\bm{\Delta}+\sum\nolimits_{g}(\mathcal{Q}_{g}-\mathcal{T}_{g}))^{-1}.$

Under Bayes-optimal condition, we have $\mathcal{Q}_{g}=\rho\mathcal{C}_{g}$
and $\mathcal{M}_{g}=\mathcal{T}_{g}=\mathcal{M}_{g}^{H}$, resulting in
$\hat{\mathcal{M}}_{g}^{H}=\hat{\mathcal{T}}_{g}^{H}=\alpha\mathcal{W}_{g}$
and $\hat{\mathcal{Q}}_{g}^{H}=\bm{0}$, $\forall g$. By further considering
the Nishimori identity introduced in the Sec. II-B of [11], we have the
relation
$\mathcal{E}G=\sum_{g}\mathcal{E}_{g}=\sum_{g}\mathcal{Q}_{g}-\mathcal{T}_{g}$.
As a consequence, after ignoring the irrelevant constant, the expression of
the free entropy function with respect to $\mathcal{E}$ can be obtained as
(7).

###  -B Proof of Theorem 2

$\displaystyle\Phi(\tau)$ $\displaystyle=-{\rm
Tr}\left(\frac{\sum_{g}\rho\sigma_{g}^{2}-\tau}{\sigma_{w}^{2}+\frac{\tau}{\alpha}}\right)-\alpha
M-\alpha
M\log\left(\sigma_{w}^{2}+\frac{\tau}{\alpha}\right)+\sum\limits_{g}\int{\rm
d}\bm{s}_{g}((1-\rho)\delta(\bm{s}_{g})+\rho Q_{g}(\bm{s}_{g}))\int{\rm
D}\bm{z}$ $\displaystyle\times\log\left(\int{\rm
d}\bm{x}_{g}((1-\rho)\delta(\bm{x}_{g})+\rho
Q_{g}(\bm{x}_{g}))\right.\left.\exp\left(-\bm{x}_{g}^{H}(\sigma_{w}^{2}+\frac{\tau}{\alpha})^{-1}\bm{x}_{g}+2\mathfrak{R}\left(\bm{x}_{g}^{H}\left((\sigma_{w}^{2}+\frac{\tau}{\alpha})^{-1}\bm{s}_{g}+(\sigma_{w}^{2}+\frac{\tau}{\alpha})^{-\frac{1}{2}}\bm{z}\right)\right)\right)\right).$
(36)

$\displaystyle\tilde{\mathcal{I}}_{g}^{(1)}$ $\displaystyle=(1-\rho)\int{\rm
D}\bm{z}\log\left((1-\rho)+\rho\int{\rm
d}\bm{x}_{g}Q_{g}(\bm{x}_{g})\exp\left(-\bm{x}_{g}^{H}(\sigma_{w}^{2}+\frac{\tau}{\alpha})^{-1}\bm{x}_{g}+2\mathfrak{R}\left(\bm{x}_{g}^{H}(\sigma_{w}^{2}+\frac{\tau}{\alpha})^{-\frac{1}{2}}\bm{z}\right)\right)\right)$
$\displaystyle=(1-\rho)\int{\rm
D}\bm{z}\log\left((1-\rho)+\rho\frac{\exp\left(||\bm{z}||^{2}\left(\sigma_{w}^{2}+\frac{\tau}{\alpha}\right)^{-1}\left(\sigma_{g}^{-2}+\left(\sigma_{w}^{2}+\frac{\tau}{\alpha}\right)^{-1}\right)^{-1}\right)}{\sigma_{g}^{-2M}\left(\sigma_{g}^{-2}+\left(\sigma_{w}^{2}+\frac{\tau}{\alpha}\right)^{-1}\right)^{M}}\right)$
$\displaystyle=(1-\rho)\int{\rm
D}\bm{z}\log\left((1-\rho)+\rho\left(\frac{\frac{1}{\alpha}\tau+\sigma_{w}^{2}}{\frac{1}{\alpha}\tau+\sigma_{w}^{2}+\sigma_{g}^{2}}\right)^{M}\exp\left(\frac{||\bm{z}||^{2}\sigma_{g}^{2}}{\sigma_{g}^{2}+\sigma_{w}^{2}+\frac{1}{\alpha}\tau}\right)\right).$
(37) $\displaystyle\tilde{\mathcal{I}}_{g}^{(2)}=$ $\displaystyle\rho\int{\rm
D}\bm{z}{\rm d}\bm{s}_{g}Q_{g}(\bm{s}_{g})\log\left((1-\rho)+\rho\int{\rm
d}\bm{x}_{g}Q_{g}(\bm{x}_{g})\exp\left(-\bm{x}_{g}^{H}(\sigma_{w}^{2}+\frac{\tau}{\alpha})^{-1}\bm{x}_{g}\right.\right.\left.\left.+2\mathfrak{R}\left(\bm{x}_{g}^{H}\left((\sigma_{w}^{2}+\frac{\tau}{\alpha})^{-1}\bm{s}_{g}+(\sigma_{w}^{2}+\frac{\tau}{\alpha})^{-\frac{1}{2}}\bm{z}\right)\right)\right)\right)$
$\displaystyle=$ $\displaystyle\rho\int{\rm
D}\bm{z}\log\left((1-\rho)+\rho\int{\rm
d}\bm{x}_{g}Q_{g}(\bm{x}_{g})\exp\left(\bm{x}_{g}^{H}(\sigma_{w}^{2}+\frac{\tau}{\alpha})^{-1}\bm{x}_{g}+2\mathfrak{R}\left(\sqrt{\frac{\sigma_{g}^{2}+\sigma_{w}^{2}+\frac{\tau}{\alpha}}{\left(\sigma_{w}^{2}+\frac{\tau}{\alpha}\right)^{2}}}\bm{x}_{g}^{H}\bm{z}\right)\right)\right)$
$\displaystyle=$ $\displaystyle\rho\int{\rm
D}\bm{z}\log\left((1-\rho)+\rho\left(\frac{\frac{1}{\alpha}\tau+\sigma_{w}^{2}}{\frac{1}{\alpha}\tau+\sigma_{w}^{2}+\sigma_{g}^{2}}\right)^{M}\exp\left(\frac{\sigma_{g}^{2}}{\sigma_{w}^{2}+\frac{\tau}{\alpha}}||\bm{z}||^{2}\right)\right).$
(38)

$\displaystyle\Phi(\mathcal{E})=$ $\displaystyle-{\rm
Tr}\left(\left(\bm{\Delta}+\frac{\mathcal{E}G}{\alpha}\right)^{-1}\left(\sum\limits_{g}\rho\mathcal{C}_{g}-\mathcal{E}G\right)\right)-\alpha
M-\alpha\log\left|\bm{\Delta}+\frac{\mathcal{E}G}{\alpha}\right|$
$\displaystyle+\sum\limits_{g}\int{\rm
d}\bm{s}_{g}((1-\rho)\delta(\bm{s}_{g})+\rho Q_{g}(\bm{s}_{g}))\int{\rm
D}\bm{z}\log\left((1-\rho)+\rho\frac{|\bm{V}_{g}|}{|\mathcal{C}_{g}|}\exp\\{\bm{\mu}_{g}^{H}\bm{V}_{g}^{-1}\bm{\mu}_{g}\\}\right).$
(39)

$\displaystyle\bm{V}_{g}=$
$\displaystyle\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}-\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}(\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}+\bm{\Delta}+\frac{1}{\alpha}\tilde{\bm{U}}\tilde{\bm{\Xi}}\tilde{\bm{U}}^{H})^{-1}\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}$
$\displaystyle=$
$\displaystyle\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}-\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}(\bm{U}_{g}(\bm{\Lambda}_{g}+\sigma_{w}^{2}\mathbf{I}_{g}+\frac{1}{\alpha}\bm{\Xi}_{g})\bm{U}_{g}^{H}+\sum\limits_{j\neq
g}\bm{U}_{j}(\sigma_{w}^{2}\mathbf{I}_{j}+\frac{1}{\alpha}\bm{\Xi}_{j})\bm{U}_{j}^{H})^{-1}\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}$
$\displaystyle=$
$\displaystyle\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}-\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}(\bm{U}_{g}(\bm{\Lambda}_{g}+\sigma_{w}^{2}\mathbf{I}_{g}+\frac{1}{\alpha}\bm{\Xi}_{g})^{-1}\bm{U}_{g}^{H}+\sum\limits_{j\neq
g}\bm{U}_{j}(\sigma_{w}^{2}\mathbf{I}_{j}+\frac{1}{\alpha}\bm{\Xi}_{j})^{-1}\bm{U}_{j}^{H})\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}$
$\displaystyle=$
$\displaystyle\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}-\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}(\bm{U}_{g}(\bm{\Lambda}_{g}+\sigma_{w}^{2}\mathbf{I}_{g}+\frac{1}{\alpha}\bm{\Xi}_{g})^{-1}\bm{U}_{g}^{H}+\sum\limits_{j\neq
g}\bm{U}_{j}(\sigma_{w}^{2}\mathbf{I}_{j})^{-1}\bm{U}_{j}^{H})\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}$
$\displaystyle=$
$\displaystyle\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}-\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}(\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}+\bm{\Delta}+\frac{1}{\alpha}\bm{U}_{g}\bm{\Xi}_{g}\bm{U}_{g}^{H})^{-1}\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}=\mathcal{C}_{g}-\mathcal{C}_{g}(\mathcal{C}_{g}+\bm{\Delta}+\frac{1}{\alpha}\mathcal{E}_{g})^{-1}\mathcal{C}_{g}.$
(40)

Substituting $\mathcal{E}=\tau G^{-1}\mathbf{I}$ into the general free entropy
equation (7) in Sec. III-A, we can get the equation (36). Do some math, we
further get $\Phi(\tau)=-\alpha
M\left(\frac{\sigma_{w}^{2}}{\frac{1}{\alpha}\tau+\sigma_{w}^{2}}+\log(\sigma_{w}^{2}+\frac{1}{\alpha}\tau)\right)+\sum\nolimits_{g}\tilde{\mathcal{I}}_{g}^{(1)}+\sum\nolimits_{g}\tilde{\mathcal{I}}_{g}^{(2)}$,
where the quantities $\tilde{\mathcal{I}}_{g}^{(1)}$ and
$\tilde{\mathcal{I}}_{g}^{(2)}$ are shown in (37)-(38). We note that the
derivation of $\tilde{\mathcal{I}}_{g}^{(2)}$ utilizes the fact that
$\bm{s}_{g}$ and $\bm{z}$ is independent. As a consequence, we get (14) in the
Sec. IV-A of the paper.

###  -C Proof of Proposition 1

By adopting the fact that $\exp{x}\leq x+1$ for $x>0$, we get following
inequalities

$\begin{split}&\exp\left(-\frac{\sigma_{g}^{2}}{\sigma_{w}^{2}+\frac{1}{\alpha}\tau}\right)<\frac{\frac{1}{\alpha}\tau+\sigma_{w}^{2}}{\frac{1}{\alpha}\tau+\sigma_{w}^{2}+\sigma_{g}^{2}},\\\
&\exp\left(-\frac{\sigma_{g}^{2}}{\sigma_{g}^{2}+\sigma_{w}^{2}+\frac{1}{\alpha}\tau}\right)>\frac{\frac{1}{\alpha}\tau+\sigma_{w}^{2}}{\frac{1}{\alpha}\tau+\sigma_{w}^{2}+\sigma_{g}^{2}}.\end{split}$

Using the law of large numbers, the term $||\bm{z}||^{2}\to M$. As a
consequence, as $M\to\infty$, we have

$\begin{split}(1-\rho)&\exp\left(-\frac{||\bm{z}||^{2}\sigma_{g}^{2}}{\sigma_{w}^{2}+\frac{1}{\alpha}\tau}\right)+\rho\left(\frac{\frac{1}{\alpha}\tau+\sigma_{w}^{2}}{\frac{1}{\alpha}\tau+\sigma_{w}^{2}+\sigma_{g}^{2}}\right)^{M}\\\
\approx&\rho\left(\frac{\frac{1}{\alpha}\tau+\sigma_{w}^{2}}{\frac{1}{\alpha}\tau+\sigma_{w}^{2}+\sigma_{g}^{2}}\right)^{M},\\\
(1-\rho)&\exp\left(-\frac{||\bm{z}||^{2}\sigma_{g}^{2}}{\sigma_{g}^{2}+\sigma_{w}^{2}+\frac{1}{\alpha}\tau}\right)+\rho\left(\frac{\frac{1}{\alpha}\tau+\sigma_{w}^{2}}{\frac{1}{\alpha}\tau+\sigma_{w}^{2}+\sigma_{g}^{2}}\right)^{M}\\\
\approx&(1-\rho)\exp\left(-\frac{\sigma_{g}^{2}}{\sigma_{g}^{2}+\sigma_{w}^{2}+\frac{1}{\alpha}\tau}\right)^{M}.\end{split}$

As a consequence, the free entropy function (14) is reduced to (15), as we can
see in the Sec. IV-A of the paper. The first order derivative of function (15)
is given by
$\frac{\partial\Phi}{\partial\tau}=\frac{1}{\tau+\alpha\sigma_{w}^{2}}(\rho\sum_{g}\frac{\sigma_{g}^{2}}{\frac{1}{\alpha}\tau+\sigma_{w}^{2}+\sigma_{g}^{2}}-\frac{\tau}{\frac{1}{\alpha}\tau+\sigma_{w}^{2}})$,
and we can prove the term
$\rho\sum_{g}\frac{\sigma_{g}^{2}}{\frac{1}{\alpha}\tau+\sigma_{w}^{2}+\sigma_{g}^{2}}$
is a monotonically decreasing with respect to $\tau$, and the term
$\frac{\tau}{\frac{1}{\alpha}\tau+\sigma_{w}^{2}}$ is monotonically increasing
with respect to $\tau$. Since there exists interaction in the regions of such
two terms, there has and only has one point that fulfills the equation
$\rho\sum_{g}\frac{\sigma_{g}^{2}}{\frac{1}{\alpha}\tau+\sigma_{w}^{2}+\sigma_{g}^{2}}=\frac{\tau}{\frac{1}{\alpha}\tau+\sigma_{w}^{2}}$,
and the point corresponds to the global maximum of the function (15).

###  -D Proof of Lemma 1

We note that the stationary point of the free energy function (7) in the Sec.
III-A with respect to every $\mathcal{E}_{g}^{(\text{s})}$ for each user group
fulfills the following equation

$\displaystyle\mathcal{E}_{g}^{(\text{s})}=$
$\displaystyle\mathbb{E}_{\bm{s}_{g},\bm{z}}\left[\left(\eta_{g}\left(\bm{s}_{g}+\bm{\Sigma}^{(\text{s})\frac{1}{2}}\bm{z}\right)-\bm{s}_{g}\right.\right)$
$\displaystyle\times\left(\left.\eta_{g}\left(\bm{s}_{g}+\bm{\Sigma}^{(\text{s})\frac{1}{2}}\bm{z}\right)-\bm{s}_{g}\right)^{H}\right],$

where we define
$\bm{\Sigma}^{(\text{s})}\triangleq\bm{\Delta}+\frac{\mathcal{E}^{(\text{s})}G}{\alpha}$.
Further denoting
$\hat{\bm{s}}_{g}^{(\text{s})}\triangleq\bm{s}_{g}+\bm{\Sigma}^{(\text{s})\frac{1}{2}}\bm{z}$,
we have

$\displaystyle\mathcal{E}_{g}^{(\text{s})}=$
$\displaystyle\mathbb{E}_{\bm{s}_{g},\hat{\bm{s}}_{g}^{(\text{s})}}\left(\left(\eta_{g}\left(\hat{\bm{s}}_{g}^{(\text{s})}\right)-\bm{s}_{g}\right)\left(\eta_{g}\left(\hat{\bm{s}}_{g}^{(\text{s})}\right)-\bm{s}_{g}\right)^{H}\right)$
$\displaystyle=$
$\displaystyle\mathbb{E}_{\hat{\bm{s}}_{g}^{(\text{s})}}\left(\mathbb{E}_{\bm{s}_{g}|\hat{\bm{s}}_{g}^{(\text{s})}}\left(\left(\eta_{g}\left(\hat{\bm{s}}_{g}^{(\text{s})}\right)-\bm{s}_{g}\right)\left(\eta_{g}\left(\hat{\bm{s}}_{g}^{(\text{s})}\right)-\bm{s}_{g}\right)^{H}\right)\right).$

Combining the channel model in Sec. V, we have

$\displaystyle
p(\bm{s}_{g}|\hat{\bm{s}}_{g}^{(\text{s})})=(1-{\hat{\rho}})\delta(\bm{s}_{g})$
$\displaystyle+$
$\displaystyle{\hat{\rho}}\mathcal{CN}\big{(}\bm{s}_{g};(\bm{\Sigma}^{(\text{s})-1}+\mathcal{C}_{g}^{-1})^{-1}\bm{\Sigma}^{(\text{s})-1}\hat{\bm{s}}_{g}^{(\text{s})},(\bm{\Sigma}^{(\text{s})-1}+\mathcal{C}_{g}^{-1})^{-1}\big{)},$

where

$\hat{\rho}\triangleq\frac{\rho\mathcal{CN}(\hat{\bm{s}}_{g}^{(\text{s})};\bm{0},\mathcal{C}_{g}+\bm{\Sigma}^{(\text{s})})}{\rho\mathcal{CN}(\hat{\bm{s}}_{g}^{(\text{s})};\bm{0},\mathcal{C}_{g}+\bm{\Sigma}^{(\text{s})})+(1-\rho)\mathcal{CN}(\hat{\bm{s}}_{g};\bm{0},\bm{\Sigma}^{(\text{s})})}.$

We further use the matrix inversion lemma, and we obtain
$(\mathcal{C}_{g}^{-1}+\bm{\Sigma}^{(\text{s})-1})^{-1}=\mathcal{C}_{g}-\mathcal{C}_{g}(\mathcal{C}_{g}+\bm{\Sigma}^{(\text{s})})^{-1}\mathcal{C}_{g}$.
Together with $\mathcal{C}_{g}=\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H}$, we
finally have
$\mathcal{E}_{g}^{(\text{s})}=\bm{U}_{g}\bm{\Xi}_{g}^{(\text{s})}\bm{U}_{g}^{H}$,
with

$\begin{split}\bm{\Xi}_{g}^{(\text{s})}&=\mathbb{E}_{\hat{\bm{s}}_{g}}\left[\hat{\rho}(1-\hat{\rho})\left(\bm{\Lambda}_{g}\bm{U}_{g}^{H}-\bm{\Lambda}_{g}\bm{U}_{g}^{H}(\mathcal{C}_{g}+\bm{\Sigma}^{(\text{s})})^{-1}\mathcal{C}_{g}\right)\right.\\\
&\times\bm{\Sigma}^{(\text{s})-1}\hat{\bm{s}}_{g}\hat{\bm{s}}_{g}^{H}\bm{\Sigma}^{(\text{s})-1}\left(\bm{U}_{g}\bm{\Lambda}_{g}-\mathcal{C}_{g}(\mathcal{C}_{g}+\bm{\Sigma}^{(\text{s})})^{-1}\bm{U}_{g}\bm{\Lambda}_{g}\right)\\\
&\left.+\hat{\rho}\left(\bm{\Lambda}_{g}-\bm{\Lambda}_{g}\bm{U}_{g}^{H}(\mathcal{C}_{g}+\bm{\Sigma}^{(\text{s})})^{-1}\bm{U}_{g}\bm{\Lambda}_{g}\right)\right].\end{split}$

###  -E Proof of Theorem 3

Note that by considering the specific functional form of
$Q_{g}(\bm{s}_{g})=\rho\mathcal{CN}(\bm{s}_{g};\bm{0},\bm{U}_{g}\bm{\Lambda}_{g}\bm{U}_{g}^{H})$,
the free entropy function (7) can be rewritten as shown in the equation (39),
where we define

$\displaystyle\bm{V}_{g}=$
$\displaystyle\mathcal{C}_{g}-\mathcal{C}_{g}(\mathcal{C}_{g}+\bm{\Delta}+\frac{1}{\alpha}G\mathcal{E})^{-1}\mathcal{C}_{g},$
$\displaystyle\bm{\mu}_{g}=$
$\displaystyle\bm{V}_{g}\left(\left(\bm{\Delta}+\frac{G\mathcal{E}}{\alpha}\right)^{-1}\bm{s}_{g}+\left(\bm{\Delta}+\frac{G\mathcal{E}}{\alpha}\right)^{-\frac{1}{2}}\bm{z}\right).$

Using the matrix inversion lemma together with the condition
$\mathcal{E}=\sum_{g}\bm{U}_{g}\bm{\Xi}_{g}\bm{U}_{g}^{H}$, we obtain the
derivations in (40), where the quantity
$\tilde{\bm{U}}\triangleq[\bm{U}_{1},\dots,\bm{U}_{G}]$ is a _tall unitary_
matrix and $\tilde{\bm{\Xi}}$ is a $\sum_{g}r_{g}$ square block diagonal
matrix with $\bm{\Xi}_{g}$ as the $g$th block. Similarly, we have
$\bm{\mu}_{g}=\bm{V}_{g}((\bm{\Delta}+\frac{\mathcal{E}_{g}}{\alpha})^{-1}\bm{s}_{g}+(\bm{\Delta}+\frac{\mathcal{E}_{g}}{\alpha})^{-\frac{1}{2}}\bm{z})$.
As a result, we get
$\Phi(\mathcal{E})=\sum_{g=0}^{G}\Phi_{g}(\mathcal{E}_{g})$, with
$\Phi_{g}(\mathcal{E}_{g})$ in the form of (21) in the Sec. V-A of the paper.

###  -F Proof of Proposition 2

The definition of the detection probability is $P_{D}={\rm
Pr}\\{\text{T}(\underline{\bm{s}_{g}})>l_{g}|a_{g}=1\\}$ with threshold
$l_{g}$. To specify the expression of $P_{D}$, we define a random variable
$C_{g,m}$ with one realization

$c_{g,m}=2|\underline{s_{g,m}}|^{2}\left(\sigma_{w}^{2}+\frac{1}{\alpha}\xi_{g,m}^{\star}+\lambda_{g,m}\right)^{-1},$

and as a result, we have $C_{g,m}\sim\chi^{2}(2)$, i.e.,
$p(c_{g,m})=\frac{1}{2}\exp(-\frac{1}{2}c_{g,m})$. Similarly, we further
define a random variable $D_{g,m}$ with one realization

$d_{g,m}=\frac{1}{2}\lambda_{g,m}c_{g,m}\left(\sigma_{w}^{2}+\frac{1}{\alpha}\xi_{g,m}^{\star}\right)^{-1}.$

Then, we get the PDF of random variable $D_{g,m}$ as

$\displaystyle p_{D_{g,m}}(d_{g,m})=$
$\displaystyle(\sigma_{w}^{2}+\frac{1}{\alpha}\xi_{g,m}^{\star})\lambda_{g,m}^{-1}$
$\displaystyle\exp(-(\sigma_{w}^{2}+\frac{1}{\alpha}\xi_{g,m}^{\star})\lambda_{g,m}^{-1}d_{g,m})\mathbb{I}(d_{g,m}>0),$

where $\mathbb{I}(d_{g,m}>0)$ is the indicator function. Note that
$p_{D_{g,m}}(d_{g,m};\omega_{g,m})$ is an exponential distribution with
parameter
$\omega_{g,m}\triangleq(\sigma_{w}^{2}+\frac{1}{\alpha}\xi_{g,m}^{\star})\lambda_{g,m}^{-1}$.
We further assume $\omega_{g,m}$ are all distinct for each $m$ which is common
in spatially correlated channel. Based on the theorem of sum of independent
exponentially distributed random variables in [29], the PDF of the sufficient
statistic with condition $a_{g}=1$ can be written as

$p(\text{T}(\underline{\bm{s}_{g}})|a_{g}=1)=\sum\limits_{m=1}^{r_{g}}\prod_{j\neq
m}\frac{\omega_{g,j}p_{D_{g,m}}(\text{T}(\underline{\bm{s}_{g}});\omega_{g,m})}{\omega_{g,j}-\omega_{g,m}},$
(41)

Combining the definition, the missed detection probability $P_{D}$ can be then
obtained via the cumulative distribution function (CDF) of (41). Similarly,
the sufficient statistic with condition $a_{g}=0$ can be formulated as

$p(\text{T}(\underline{\bm{s}_{g}})|a_{g}=0)=\sum_{m=1}^{r_{g}}\prod_{j\neq
m}\frac{\tilde{\omega}_{g,j}p_{D_{g,m}}(\text{T}(\underline{\bm{s}_{g}});\tilde{\omega}_{g,m})}{\tilde{\omega}_{g,j}-\tilde{\omega}_{g,j}},$
(42)

where
$\tilde{\omega}_{g,m}\triangleq(\sigma_{w}^{2}+\frac{1}{\alpha}\xi_{g,m}^{\star}+\lambda_{g,m})\lambda_{g,m}^{-1}$.
As a consequence, based on the CDF of (41), the expression of false alarm
probability with definition $P_{F}={\rm
Pr}\\{\text{T}(\underline{\bm{s}_{g}})>l_{g}|a_{g}=0\\}$ can be then
specified, as shown in (25) in the Sec.V-B of the paper.

## References

  * [1] A. T. Abebe and C. G. Kang, ``Comprehensive grant-free random access for massive low latency communication,'' in _Proc. IEEE Int. Conf. Commun. (ICC)_ , 2017, pp. 1–6.
  * [2] K. Senel and E. G. Larsson, ``Grant-free massive MTC-enabled massive MIMO: A compressive sensing approach,'' _IEEE Trans. Commun._ , vol. 66, no. 12, pp. 6164–6175, Dec. 2018.
  * [3] P. Salvo Rossi, D. Ciuonzo, K. Kansanen, and T. Ekman, ``Performance analysis of energy detection for MIMO decision fusion in wireless sensor networks over arbitrary fading channels,'' _IEEE Trans. Wireless Commun._ , vol. 15, no. 11, pp. 7794–7806, 2016.
  * [4] D. Ciuonzo, A. Aubry, and V. Carotenuto, ``Rician MIMO channel- and jamming-aware decision fusion,'' _IEEE Trans. Signal Process._ , vol. 65, no. 15, pp. 3866–3880, 2017.
  * [5] 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, 2020.
  * [6] 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.
  * [7] 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, 2020.
  * [8] J. Ahn, B. Shim, and K. B. Lee, ``EP-based joint active user detection and channel estimation for massive machine-type communications,'' _IEEE Trans. Commun._ , vol. 67, no. 7, pp. 5178–5189, 2019.
  * [9] J.-C. Jiang and H.-M. Wang, ``Massive random access with sporadic short packets: Joint active user detection and channel estimation via sequential message passing,'' _IEEE Trans. Wireless Commun._ , vol. 20, no. 7, pp. 4541–4555, 2021.
  * [10] L. Liu, E. G. Larsson, W. Yu, P. Petar _et al._ , ``Sparse signal processing for grant-free massive connectivity: A future paradigm for random access protocols in the internet of things,'' _IEEE Signal Process. Mag._ , vol. 35, no. 5, pp. 88–99, Sep. 2018.
  * [11] F. Krzakala, M. M zard, F. Sausset, Y. Sun, and L. Zdeborov , ``Probabilistic reconstruction in compressed sensing: Algorithms, phase diagrams, and threshold achieving matrices,'' _Journal of Statistical Mechanics: Theory and Experiment_ , 2012.
  * [12] J. Zhu and D. Baron, ``Performance limits with additive error metrics in noisy multimeasurement vector problems,'' _IEEE Trans. Signal Process._ , vol. 66, no. 20, pp. 5338–5348, 2018.
  * [13] G. Hannak, A. Perelli, N. Goertz, G. Matz, and M. E. Davies, ``Performance analysis of approximate message passing for distributed compressed sensing,'' _IEEE J. Sel. Topics Signal Process._ , vol. 12, no. 5, pp. 857–870, 2018.
  * [14] T. Tanaka, ``A statistical-mechanics approach to large-system analysis of CDMA multiuser detectors,'' _IEEE Trans. Inf. Theory_ , vol. 48, no. 11, pp. 2888–2910, 2002.
  * [15] M. Mézard and A. Montanari, _Information, Physics, and Computation_. Oxford, 2009.
  * [16] A. Adhikary, J. Nam, J. Ahn, and G. Caire, ``Joint spatial division and multiplexing-the large-scale array regime,'' _IEEE Trans. Inf. Theory_ , vol. 59, no. 10, pp. 6441–6463, 2013.
  * [17] C. Wen, S. Jin, K. Wong, J. Chen, and P. Ting, ``Channel estimation for massive MIMO using Gaussian-mixture bayesian learning,'' _IEEE Trans. Wireless Commun._ , vol. 14, no. 3, pp. 1356–1368, 2015.
  * [18] H. Yin, D. Gesbert, M. Filippou, and Y. Liu, ``A coordinated approach to channel estimation in large-scale multiple-antenna systems,'' _IEEE J. Sel. Areas in Commun._ , vol. 31, no. 2, pp. 264–273, 2013.
  * [19] G. Hannak, M. Mayer, A. Jung, G. Matz, and N. Goertz, ``Joint channel estimation and activity detection for multiuser communication systems,'' in _Proc. IEEE Int. Conf. Commun. Workshop (ICCW)_ , Jun. 2015, pp. 2086–2091.
  * [20] M. Ke, Z. Gao, Y. Wu, X. Gao, and K.-K. Wong, ``Massive access in cell-free massive MIMO-based internet of things: Cloud computing and edge computing paradigms,'' 2020, [Online]. Available: https://arxiv.org/abs/2007.09617, preprint.
  * [21] S. Rangan, A. K. Fletcher, and V. K. Goyal, ``Asymptotic analysis of MAP estimation via the replica method and applications to compressed sensing,'' _IEEE Trans. Inf. Theory_ , vol. 58, no. 3, pp. 1902–1923, 2012\.
  * [22] H. He, C. Wen, and S. Jin, ``Bayesian optimal data detector for hybrid mmWave MIMO-OFDM systems with low-resolution ADCs,'' _IEEE J. Sel. Topics Signal Process._ , vol. 12, no. 3, pp. 469–483, 2018.
  * [23] H. Wang, C. Wen, and S. Jin, ``Bayesian optimal data detector for mmWave OFDM system with low-resolution ADC,'' _IEEE J. Sel. Areas Commun._ , vol. 35, no. 9, pp. 1962–1979, 2017.
  * [24] M. Bayati and A. Montanari, ``The dynamics of message passing on dense graphs, with applications to compressed sensing,'' _IEEE Trans. Inf. Theory_ , vol. 57, no. 2, pp. 764–785, Feb. 2011.
  * [25] S. Rangan, P. Schniter, A. K. Fletcher, and S. Sarkar, ``On the convergence of approximate message passing with arbitrary matrices,'' _IEEE Trans. Inf. Theory_ , vol. 65, no. 9, pp. 5339–5351, 2019.
  * [26] 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.
  * [27] E. Björnson, J. Hoydis, and L. Sanguinetti, ``Massive MIMO networks: Spectral, energy, and hardware efficiency,'' _Foundations and Trends in Signal Processing_ , vol. 11, no. 3-4, pp. 154–655, 2017.
  * [28] J. Barbier and F. Krzakala, ``Approximate message-passing decoder and capacity achieving sparse superposition codes,'' _IEEE Trans. Inf. Theory_ , vol. 63, no. 8, pp. 4894–4927, 2017.
  * [29] M. Bibinger, ``Notes on the sum and maximum of independent exponentially distributed random variables with different scale parameters,'' 2013, [Online]. Available: https://arxiv.org/abs/1307.3945, preprint.