# On the Impact of Channel Estimation on the Design and Analysis of IRSA based
Systems

Chirag Ramesh Srivatsa, and Chandra R. Murthy Chirag Ramesh Srivatsa is with
the Robert Bosch Centre for Cyber-Physical Systems, Indian Institute of
Science, Bangalore, India. Chandra R. Murthy is with the Dept. of Electrical
Communication Engineering, Indian Institute of Science, Bangalore, India
(e-mail: {chiragramesh, cmurthy}@iisc.ac.in).This work has been presented in
part in [1]. 0000-0002-3732-4733 0000-0003-4901-9434

###### Abstract

Irregular repetition slotted aloha (IRSA) is a distributed grant-free random
access protocol where users transmit multiple replicas of their packets to a
base station (BS). The BS recovers the packets using successive interference
cancellation. In this paper, we first derive channel estimates for IRSA,
exploiting the sparsity structure of IRSA transmissions, when non-orthogonal
pilots are employed across users to facilitate channel estimation at the BS.
Allowing for the use of non-orthogonal pilots is important, as the length of
orthogonal pilots scales linearly with the total number of devices, leading to
prohibitive overhead as the number of devices increases. Next, we present a
novel analysis of the throughput of IRSA under practical channel estimation
errors, with the use of multiple antennas at the BS. Finally, we theoretically
characterize the asymptotic throughput performance of IRSA using a density
evolution based analysis. Simulation results underline the importance of
accounting for channel estimation errors in analyzing IRSA, which can even
lead to 70% loss in performance in severely interference-limited regimes. We
also provide novel insights on the effect of parameters such as pilot length,
SNR, number of antennas at the BS, etc, on the system throughput.

###### Index Terms:

Irregular repetition slotted aloha, pilot contamination, density evolution,
channel estimation

##  I Introduction

Massive machine-type communications (mMTC) is an evolving 5G use-case,
expected to serve around $10^{6}$ devices per square kilometer [2]. The users
in mMTC applications are sporadically active and transmit short packets to a
central base station (BS) [3]. Grant-free random access (GFRA) protocols are
appropriate in mMTC applications since they incur a low control and signaling
overhead [4, 5]. Typically, in these protocols, users transmit packets
(consisting of a header containing pilot symbols followed by the data payload)
by randomly accessing resource elements (REs). Since the length of orthogonal
pilots scales linearly with the number of users, the overhead of assigning
orthogonal pilots becomes prohibitively expensive. Thus, pilot contamination
is inevitable due to the use of non-orthogonal pilots, and has to be accounted
for while analyzing the performance of GFRA protocols for mMTC.

One popular GFRA protocol is irregular repetition slotted aloha (IRSA) [6, 7],
which is the focus of this paper. Users in IRSA transmit replicas of their
packets on a randomly selected subset of the available REs. The indices of the
REs in which they transmit make up the access pattern matrix (APM). Existing
works in IRSA assume availability of perfect channel state information (CSI)
at the BS, which is difficult to achieve, especially when non-orthogonal
pilots are employed. Channel estimation errors and pilot contamination due to
non-orthogonal pilots can erase much of the gains promised by IRSA protocols.
Thus, one of the main goals of this paper is to understand the impact of
estimated CSI on the performance of IRSA when non-orthogonal pilots are used.

The decoding in IRSA is an iterative process involving successive interference
cancellation (SIC) [8], where the users are decoded via a combination of
inter-RE and intra-RE SIC [9]. _Inter-RE SIC_ refers to the removal of packet
replicas from a different RE than the one the packet was decoded in, while
_intra-RE SIC_ refers to the removal of a packet from the same RE in which the
packet was decoded, in order to facilitate decoding additional packets that
may have been transmitted in that RE. The conventional version of IRSA used
only inter-RE SIC to decode users and assumed a collision model, wherein only
singleton REs can be decoded [6]. Here, a singleton RE refers to an RE where a
single user’s packet is received without collision. Since no packets can be
decoded in REs where collisions occur, the maximum possible throughput is one
packet per RE, the same as the throughput with perfectly coordinated multiple
access. This maximum can be achieved asymptotically as the number of users and
REs go to infinity, when the soliton distribution is used to generate the
repetition factors of the users [10].

When the BS is equipped with multiple antennas, it can potentially decode
multiple packets in a single RE, i.e., if the signal to interference plus
noise ratios (SINRs) of the packets are sufficiently high. Thus, using an
_SINR-threshold_ model has also been considered for IRSA, where users can be
decoded if and only if their SINR is higher than a predetermined threshold
[11]. After decoding users with sufficiently high SINRs, with a combination of
intra-RE and inter-RE SIC, the packet replicas of the decoded users can be
removed from all the REs in which they have transmitted packets. Then, all the
REs can be revisited to see if further users can be decoded from the residual
signal. This procedure is continued iteratively until no further packets can
be decoded. This yields a higher throughput compared to the collision model,
and can potentially achieve a throughput greater than one packet per RE. Thus,
a second goal of this paper is to characterize the performance of IRSA under
estimated CSI as a function of system parameters such as the number of
antennas at the BS, the pilot length, the SINR threshold, etc.

The throughput of the IRSA family of multiple access protocols is analyzed
using the density evolution (DE) approach, wherein two probability densities
are obtained as functions of each other [6]. This iterative recipe provides
the asymptotic performance of the system. The asymptotic throughput has been
obtained for IRSA via DE for the MAC [10], accounting for path loss [9], for
the scalar Rayleigh fading channel [11], with multiuser detectors [12], for
the polarized MIMO channel in satellite networks [13], and other enhanced
variants of IRSA [14, 15]. The authors in [16] have proposed an algorithm to
detect the subset of active users in IRSA. In contrast, this work focuses on
the effect of channel estimation errors on the performance of IRSA. A closely
related protocol is pattern division multiple access (PDMA) [17], where users
transmit their packets across a subset of REs governed by a binary APM. A
difference with PDMA is that the APM is designed in a centralized manner to
maximize the so-called _constellation-constrained capacity,_ which is not
scalable to a massive number of users in mMTC. Thus, a theoretical analysis of
the throughput of the IRSA protocol under pilot contamination, accounting for
the effect of channel estimation errors, path loss, fading, and multiple
antennas at the BS, is not yet available in the literature, to the best of our
knowledge.

###  I-A Contributions

Our main contributions in this paper are as follows:

  1. 1.

We derive channel estimates for IRSA under three schemes: the first one
exploits the sparsity in the APM to estimate the channels of the users, and
the other two assume knowledge of the APM and output minimum mean square error
(MMSE) estimates. (See Theorem 1 in Sec. III.)

  2. 2.

We present a novel analysis of the SINR in IRSA accounting for channel
estimation errors, where estimates are acquired via non-orthogonal pilots
under the three estimation schemes. (See Theorem 2 in Sec. IV.)

  3. 3.

We theoretically analyze the throughput of IRSA via DE, when users perform
path loss inversion based power control. The analysis reveals the asymptotic
performance of the protocols as the number of users and REs get large. (See
Theorem 3 in Sec. V-C and also Sec. V-D.)

Through extensive simulations, we show that channel estimation errors lead to
a significant loss of throughput compared to the ideal scenario with perfect
CSI at the BS, even resulting in up to 70% loss in severely interference-
limited regimes. In particular, in mMTC applications, since it is not possible
to assign orthogonal pilots to all users, the resulting pilot contamination
can significantly degrade the SINR, leading to poor performance. On the
positive side, this loss in performance can be recuperated by optimizing
system parameters such as pilot length, number of antennas, frame length,
signal to noise ratio, and SINR threshold. In particular, we show that the
pilot length required to obtain near-optimal performance is orders of
magnitude lower than the pilot length needed to assign orthogonal pilots to
all users. For example, a pilot length of $\tau=12$ is sufficient to obtain
optimal performance with $M=150$ users, whereas the use of orthogonal pilot
sequences requires $\tau=150$ pilot symbols (See Fig. 1). This is possible
because only a small fraction of users transmit in a given RE in IRSA;
exploiting this sparsity in user access allows one to obtain accurate channel
estimates even when the pilots are non-orthogonal (See Algorithm 1.)

Our analysis also reveals an inflection load, beyond which the system becomes
interference-limited, resulting in a dramatic reduction of the throughput. The
asymptotic throughput obtained via DE serves as an upper bound for the
achievable throughput, and facilitates numerical optimization of the
throughput with respect to the system parameters.

Notation: The symbols $a$, $\mathbf{a}$, $\mathbf{A}$,
$[\mathbf{A}]_{i,\mathrel{\mathop{\mathchar 58\relax}}}$,
$[\mathbf{A}]_{\mathrel{\mathop{\mathchar 58\relax}},j}$, $\mathbf{0}_{N}$,
$\mathbf{1}_{N},$ and $\mathbf{I}_{N}$ denote a scalar, a vector, a matrix,
the $i$th row of $\mathbf{A}$, the $j$th column of $\mathbf{A}$, all-zero
vector of length $N$, all ones vector of length $N$, and an identity matrix of
size $N\times N$, respectively. $[\mathbf{a}]_{\mathcal{S}}$ and
$[\mathbf{A}]_{\mathrel{\mathop{\mathchar 58\relax}},\mathcal{S}}$ denote the
elements of $\mathbf{a}$ and the columns of $\mathbf{A}$ indexed by the set
$\mathcal{S}$ respectively. $\text{diag}(\mathbf{a})$ is a diagonal matrix
with diagonal entries given by $\mathbf{a}$. The set of real and complex
matrices of size $N\times M$ are denoted as $\mathop{{}\mathbb{R}}^{N\times
M}$ and $\mathop{{}\mathbb{C}}^{N\times M}$.
$\mathcal{N}(\mathbf{a},\mathbf{A})$ and $\mathcal{CN}(\mathbf{a},\mathbf{A})$
denote the real and complex Gaussian distribution, respectively, with mean
$\mathbf{a}$ and covariance $\mathbf{A}$. $[N]$ denotes the set
$\\{1,2,\ldots,N\\}$. $|\cdot|$, $\|\cdot\|$, $[\cdot]^{T}$, $[\cdot]^{*}$,
$[\cdot]^{H}$, $\mathop{{}\mathbb{E}}[\cdot]$, and
$\mathop{{}\mathbb{E}_{\mathbf{a}}}\left[\cdot\right]$ denote the magnitude
(or cardinality of a set), $\ell_{2}$ norm, transpose, conjugate, hermitian,
expectation, and the expectation conditioned on $\mathbf{a}$, respectively.

##  II System Model

An IRSA system is considered with $M$ single-antenna users communicating with
a central BS equipped with $N$ antennas. The users are assumed to be
arbitrarily located within a cell, with the BS located at the cell center. The
fading is modeled as block-fading, quasi-static and Rayleigh distributed. The
time-frequency resource is divided into REs, and $T$ REs together constitute a
_frame_. The REs can be slots, subcarriers or both. In each frame, the users
contend for the channel by randomly selecting a subset of REs, and they
transmit replicas of their packets in the selected REs. Each packet replica
comprises of a header containing pilot symbols and payload containing data and
error correction symbols.

The access of REs in a given frame by all the users can be represented by a
binary access pattern matrix (APM) $\mathbf{G}\in\\{0,1\\}^{T\times M}$. The
entries of $\mathbf{G}$ are denoted by $g_{tm}=[\mathbf{G}]_{tm}$, and
$g_{tm}=1$ if the $m$th user transmits its packet in the $t$th RE, and
$g_{tm}=0$ otherwise. The $m$th user samples their repetition factor $d_{m}$
from a preset probability distribution. They then choose $d_{m}$ REs from a
total of $T$ REs uniformly at random for transmission. We note that, due to
the distributed nature of the protocol, the $M$ columns of $\mathbf{G}$ are
i.i.d., and $\mathbf{G}$ is independently generated from one frame to the
next.

At the BS, the received signal in the $t$th RE is a superposition of the
packets transmitted by the users that are scheduled to transmit in the same
RE. In the pilot phase, if $g_{tm}=1$, the $m$th user transmits a
$\tau$-length pilot $\mathbf{p}_{m}\in\mathbb{C}^{\tau}$ in the $t$th RE, with
each pilot symbol transmitted at an average power $P_{\tau}$, and thus,
$\mathbb{E}[\|\mathbf{p}_{m}\|^{2}]=\tau P_{\tau}$. The pilot signal
${\mathbf{Y}_{t}^{p}}\in\mathbb{C}^{N\times\tau}$ received at the BS using its
$N$ antennas and in the $t$th RE is given by

$\displaystyle{\mathbf{Y}_{t}^{p}}$
$\displaystyle=\sum\nolimits_{m=1}^{M}g_{tm}\mathbf{h}_{tm}\mathbf{p}_{m}^{H}+{\mathbf{N}_{t}^{p}},$
(1)

where ${\mathbf{N}_{t}^{p}}\in\mathbb{C}^{N\times\tau}$ is the complex
additive Gaussian noise at the BS with
$[{\mathbf{N}_{t}^{p}}]_{nj}\mathrel{\overset{\text{i.i.d.}}{\scalebox{1.5}[1.0]{$\sim$}}}\mathcal{CN}(0,N_{0})$
$\forall\ n\in[N],\ j\in[\tau]$ and $t\in[T]$, and $N_{0}$ is the noise
variance. Here $\mathbf{h}_{tm}=[h_{tm1},\ldots,h_{tmN}]^{T}$ is the uplink
channel vector of the $m$th user in the $t$th RE, with
$h_{tmn}\mathrel{\overset{\text{i.i.d.}}{\scalebox{1.5}[1.0]{$\sim$}}}\mathcal{CN}(0,\beta_{m}\sigma_{h}^{2}),\
\forall\ t\in[T],\ m\in[M]$ and $n\in[N]$, where $\beta_{m}$ is the path loss
coefficient and $\sigma_{h}^{2}$ is the fading variance.

In the data phase, users transmit their data symbols. Considering one of the
data symbols, the $m$th user transmits a data symbol $x_{m}$ with
$\mathbb{E}[x_{m}]=0$ and $\mathbb{E}[|x_{m}|^{2}]=P$, i.e., with transmit
power $P$. The corresponding received data signal
$\mathbf{y}_{t}\in\mathbb{C}^{N}$ at the BS in the $t$th RE is

$\displaystyle\mathbf{y}_{t}$
$\displaystyle=\sum\nolimits_{m=1}^{M}g_{tm}\mathbf{h}_{tm}x_{m}+\mathbf{n}_{t},$
(2)

where $\mathbf{n}_{t}\in\mathbb{C}^{N}$ is the complex additive white Gaussian
noise at the BS with
$[\mathbf{n}_{t}]_{n}\mathrel{\overset{\text{i.i.d.}}{\scalebox{1.5}[1.0]{$\sim$}}}\mathcal{CN}(0,N_{0}),$
$\forall\ n\in[N]$ and $t\in[T]$.

###  II-A SIC-based Decoding

The received data is processed iteratively at the BS. The BS computes channel
estimates for all users in all REs using the pilot symbols. It uses these
channel estimates to combine the received data signal across the BS antennas
and attempts to decode the user’s data packet, treating interference from
other users as noise. If it successfully decodes any user, which can be
verified via a cyclic redundancy check, it performs SIC in all REs in which
that user has transmitted, with both inter-RE and intra-RE SIC. The BS
proceeds with the next iteration, where the channels are re-estimated for the
remaining users, and this decoding process proceeds iteratively.

In this work, we abstract the decoding of a user’s packet using an SINR-
threshold model. That is, if the SINR of a packet in a given RE in any
decoding iteration exceeds a threshold $\gamma_{\text{th}}$, then the packet
can be decoded correctly [9, 11]. _Packet capture_ occurs when a packet can be
decoded correctly as per the SINR-threshold model, even though it collides
with another packet, and is thus considered a good abstraction of the decoding
in the physical layer.

We now describe the performance evaluation of IRSA via the SINR-threshold
model. We first compute channel estimates and SINR achieved by all users in
all REs. If we find a user with SINR $\geq\gamma_{\text{th}}$ in some RE, we
mark the data packet as having been decoded successfully and remove the
contribution of the user’s packet from all REs that contain a replica of that
packet. In the next iteration, the channels are re-estimated from the residual
pilot symbols after SIC, the SINRs are recomputed in all REs, and the decoding
of users’ packets continues. The decoding process proceeds in iterations and
stops when no additional users are decoded in two successive iterations. The
system throughput $\mathcal{T}$ is calculated as the number of correctly
decoded unique packets divided by the number of REs.

The crucial step in evaluating the performance of the above decoding procedure
is the calculation of the SINR at the receiver, which depends on the CSI
available at the BS. Sec. III describes the channel estimation process. The
detailed derivation of the SINR is presented later in Sec. IV.

##  III Channel Estimation

In this section, the channel estimates for all users are derived under three
schemes. The first scheme, termed the sparsity-based estimation scheme,
estimates both the APM and the channels of the users. In contrast with this,
the other schemes exploit the knowledge of $\mathbf{G}$ and output MMSE
estimates. This is not a strong assumption and can be made possible by using
pseudo-random pattern matrices generated from a seed that is available at the
BS and the users.

Channel estimation is performed based on the received pilot signal, which
contains the pilots transmitted by all the users who have transmitted in that
RE. The estimates are recomputed in every iteration, and hence the signals and
channel estimates are indexed by the decoding iteration $k$. Let the set of
users who have not yet been decoded in the first $k-1$ iterations be denoted
by $\mathcal{S}_{k}$, and for some $m\in\mathcal{S}_{k}$, let
$\mathcal{S}_{k}^{m}\triangleq\mathcal{S}_{k}\setminus\\{m\\}$, with
$\mathcal{S}_{1}=[M]$. The received pilot signal at the BS, in the $t$th RE,
and during the $k$th decoding iteration, is given by

$\displaystyle{\mathbf{Y}_{t}^{pk}}$
$\displaystyle=\sum\nolimits_{i\in\mathcal{S}_{k}}g_{ti}\mathbf{h}_{ti}\mathbf{p}_{i}^{H}+{\mathbf{N}_{t}^{p}}.$
(3)

We now discuss three channel estimation schemes for IRSA.

###  III-A Sparsity-based APM and Channel Estimation

The first scheme is the sparsity-based estimation scheme in which we estimate
the APM and the channels in each decoding iteration. We consider the conjugate
transpose of the received pilot signal in the $t$th RE from (3) as
$\overline{\mathbf{Y}}_{t}^{pk}\triangleq{\mathbf{Y}_{t}^{pkH}}$, with
$\overline{\mathbf{N}}_{t}\triangleq{\mathbf{N}_{t}^{pH}}$. Let
$\mathbf{P}\in\mathbb{C}^{\tau\times M}$ contain the known pilots of the $M$
users as its columns and
$\mathbf{P}^{k}=[\mathbf{P}]_{\mathrel{\mathop{\mathchar
58\relax}},\mathcal{S}_{k}}$. The signal $\overline{\mathbf{Y}}_{t}^{pk}$ can
be factorized into the product of two matrices as follows:

$\displaystyle\underbrace{\overline{\mathbf{Y}}_{t}^{pk}}_{\tau\times N}$
$\displaystyle=\underbrace{\left[\mathbf{p}_{{i_{1}}},\ldots,\mathbf{p}_{{i_{M^{k}}}}\right]}_{\mathbf{P}^{k}}\underbrace{\begin{bmatrix}{g_{t{i_{1}}}}\mathbf{h}_{t{i_{1}}}^{H}\\\
\vdots\\\
{g_{t{i_{M^{k}}}}}\mathbf{h}_{t{i_{M^{k}}}}^{H}\end{bmatrix}}_{\mathbf{Z}_{t}^{k}}+\underbrace{\overline{\mathbf{N}}_{t}}_{\tau\times
N},$ (4)

where $\mathcal{S}_{k}=\\{i_{1},i_{2},\ldots,i_{M^{k}}\\}$, with
$M^{k}=|\mathcal{S}_{k}|$. Here, $\mathbf{Z}_{t}^{k}\in\mathbb{C}^{M^{k}\times
N}$ contains the $t$th row of the unknown APM $\mathbf{G}$, and the unknown
channels. The rows of $\mathbf{Z}_{t}^{k}$ are either all-zero or all-nonzero
depending on whether the corresponding $g_{ti}=0$ or $1$. This results in an
under-determined system of equations, where the columns of
$\mathbf{Z}_{t}^{k}$ share the same support. This structure is called as a
multiple measurement vector (MMV) recovery problem in compressed sensing. The
estimation of $\mathbf{Z}_{t}^{k}$ from (4) can be performed using well known
MMV recovery algorithms from compressed sensing literature to recover
$\\{g_{ti}\\}$ in the each of the $T$ REs.

Multiple sparse Bayesian learning111Any MMV algorithm can be used to recover
joint-sparse columns of $\mathbf{Z}_{t}^{k}$, but we use MSBL due to its high
performance. MSBL also outputs a ”plug-in” MMSE channel estimate which can
then be used to find a meaningful SINR expression since the estimate is
uncorrelated with the estimation error [18]. (MSBL) [19] is an empirical
Bayesian algorithm that can recover $\mathbf{Z}_{t}^{k}$ from linear under-
determined observations $\overline{\mathbf{Y}}_{t}^{pk}$. In MSBL, a Gaussian
prior is imposed on the columns of $\mathbf{Z}_{t}^{k}$ as

$\displaystyle
p(\mathbf{Z}_{t}^{k};\\!\bm{\gamma}_{kt})\\!=\\!\textstyle{\prod\nolimits_{n=1}^{N}}p([\mathbf{Z}_{t}^{k}]_{\mathrel{\mathop{\mathchar
58\relax}},n};\\!\bm{\gamma}_{kt})\\!=\\!\textstyle{\prod\nolimits_{n=1}^{N}}\mathcal{CN}(\mathbf{0}_{M_{t}},\bm{\Gamma}_{kt}),\\!$
(5)

where $\bm{\Gamma}_{kt}=$ diag$(\bm{\gamma}_{kt})$ and the columns of
$\mathbf{Z}_{t}^{k}$ are i.i.d. The elements of
$\bm{\gamma}_{kt}\in\mathbb{R}_{+}^{M^{k}}$ are unknown hyperparameters for
the undecoded users. Recovering the hyperparameters would yield $g_{tm}$ since
$[\bm{\gamma}_{kt}]_{m}$ models the variance of the $m$th user’s channel in
the $t$th RE. The hyperparameters are estimated by iteratively maximizing the
log-likelihood $\log p(\overline{\mathbf{Y}}_{t}^{pk};\bm{\gamma}_{kt})$, with
$p(\overline{\mathbf{Y}}_{t}^{pk};\bm{\gamma}_{kt})=\prod_{n=1}^{N}p([\overline{\mathbf{Y}}_{t}^{pk}]_{\mathrel{\mathop{\mathchar
58\relax}},n};\bm{\gamma}_{kt})$.

Input :
$\tau,N,N_{0},\mathcal{S}_{k},\mathbf{P},\overline{\mathbf{Y}}_{t}^{pk},\gamma_{\text{pr}},j_{\max}$

1 Compute: $M^{k}=|\mathcal{S}_{k}|,$
$\mathbf{P}^{k}=[\mathbf{P}]_{\mathrel{\mathop{\mathchar
58\relax}},\mathcal{S}_{k}}$

2 Initialize: $\bm{\gamma}_{kt}^{0}=\bm{1}_{M^{k}}$

3 for  _$j=0,1,2,\ldots,j_{\max}$_ do

4 Compute ${\bm{\Gamma}}_{kt}^{j}=\text{diag}(\bm{\gamma}_{kt}^{j})$

5
$\bm{\Sigma}_{kt}^{j+1}\\!\\!=\\!\bm{\Gamma}_{kt}^{j}-\bm{\Gamma}_{kt}^{j}\mathbf{P}^{kH}\\!(N_{0}\mathbf{I}_{\tau}\\!+\\!\mathbf{P}^{k}\bm{\Gamma}_{kt}^{j}\mathbf{P}^{kH})^{-1}\mathbf{P}^{k}\bm{\Gamma}_{kt}^{j}$

6
$\bm{\mu}_{ktn}^{j+1}=N_{0}^{-1}\bm{\Sigma}_{kt}^{j+1}\mathbf{P}^{kH}[\overline{\mathbf{Y}}_{t}^{pk}]_{\mathrel{\mathop{\mathchar
58\relax}},n},\ 1\leq n\leq N$

7
$\textstyle{[{\bm{\gamma}}_{kt}^{j+1}]_{i}=\dfrac{1}{N}\sum\limits_{n=1}^{N}([\bm{\Sigma}_{kt}^{j+1}]_{i,i}+|[\bm{\mu}_{ktn}^{j+1}]_{i}|^{2}),\
\\!\forall\ i\in[M^{k}]}$

8 end for

$\textbf{Output:}\quad\hat{{g}}_{tm}^{k}=\begin{cases}1,\
[{\bm{\gamma}}_{kt}^{j_{\max}}]_{m}\geq\gamma_{\text{pr}}\\\ 0,\
[{\bm{\gamma}}_{kt}^{j_{\max}}]_{m}<\gamma_{\text{pr}}\end{cases}\\!\\!,\
\forall\ m\in[M^{k}]$,
$\hat{\mathbf{Z}}_{t}^{k}=[\bm{\mu}_{kt1}^{j_{\max}}\bm{\mu}_{kt2}^{j_{\max}}\ldots\bm{\mu}_{ktN}^{j_{\max}}]$

Algorithm 1 APM and Channel Estimation in $t$th RE

Let $j$ denote the iteration index in MSBL. Stated in our notation, the
overall estimation procedure is summarized in Algorithm 1. The MSBL algorithm
converges to a saddle point or a local maximizer of the overall log-likelihood
[19]. Further, the MSBL algorithm has been empirically shown to correctly
recover the support of $\mathbf{Z}_{t}^{k}$, provided $\tau$ and $N$ are large
enough [19], if the signal to noise ratio is good enough. The algorithm is run
for $j_{\max}$ iterations in each of the $T$ REs. As the iterations proceed,
the hyperparameters corresponding to users with $g_{ti}=0$ converge to zero,
resulting in sparse estimates. At the end of the iterations, the estimated
coefficient $\hat{g}_{tm}^{k}$ for the $m$th user in the $t$th RE in the $k$th
decoding iteration is obtained by thresholding
$[{\bm{\gamma}}_{kt}^{j_{\max}}]_{m}$ at a value $\gamma_{\text{pr}}$. This
can result in errors in estimating $g_{ti}$, and the errors in APM estimation
can be described by

$\displaystyle\mathcal{F}_{t}^{k}$ $\displaystyle=\\{i\in[M^{k}]\ |\
\hat{g}_{ti}^{k}(1-{g}_{ti})=1\\},$ (6a) $\displaystyle\mathcal{M}_{t}^{k}$
$\displaystyle=\\{i\in[M^{k}]\ |\ (1-\hat{g}_{ti}^{k}){g}_{ti}=1\\},$ (6b)

where $\mathcal{F}_{t}^{k}$ is the set of false positive users, and
$\mathcal{M}_{t}^{k}$ is the set of false negative users. These errors affect
decoding of other users in two ways: both kinds of users contaminate the
channel estimates of other users, and users in $\mathcal{M}_{t}^{k}$ interfere
with the data decoding of other users as well. The effect of errors in
detection of users is described in detail in [16].

The algorithm also outputs the maximum aposteriori probability estimates of
the channels $\hat{\mathbf{Z}}_{t}^{k}$ in each of the $T$ REs. The estimate
$\hat{\mathbf{H}}_{t}^{k}=\hat{\mathbf{Z}}_{t}^{kH}\in\mathbb{C}^{N\times
M^{k}}$ of the channels of the $M^{k}$ users is described in Theorem 1 and can
be calculated as
$\hat{\mathbf{H}}_{t}^{k}={{\mathbf{Y}}^{pk}_{t}}\mathbf{P}^{k}{\hat{\bm{\Gamma}}_{kt}}(\mathbf{P}^{kH}\mathbf{P}^{k}{\hat{\bm{\Gamma}}_{kt}}+N_{0}\mathbf{I}_{M^{k}})^{-1},$
where ${\hat{\bm{\Gamma}}_{kt}}=$ diag$({\bm{\gamma}}_{kt}^{j_{\max}})$. This
estimate is a "plug-in" MMSE estimate and it contains estimates for
erroneously detected users as well. An added advantage of MSBL is that the
path loss coefficient can be calculated by averaging the estimated
hyperparameters across REs as
$\hat{\beta}_{i}^{k}=(\sum\nolimits_{t=1}^{T}\hat{g}_{ti}^{k}[{\bm{\gamma}}_{kt}^{j_{\max}}]_{i})/(\sigma_{h}^{2}\sum\nolimits_{t=1}^{T}\hat{g}_{ti}^{k})$.
Thus, Algorithm 1 does not require any prior information about the APM or
$\\{\beta_{i}\\}_{i=1}^{M}$ to estimate the channels.

###  III-B MMSE Channel Estimation with Known APM

We now derive the MMSE channel estimates for all users in each RE, exploiting
the knowledge of the APM $\mathbf{G}$ and $\\{\beta_{i}\\}_{i=1}^{M}$. By
using a common seed at the BS and the users, the APM can be generated at the
BS and thus, we can assume that the BS has knowledge of $\mathbf{G}$. Let
$\mathcal{G}_{t}=\\{i\in[M]|g_{ti}=1\\}$ be the set of users who have
transmitted in the $t$th RE. Let
$M_{t}^{k}=|\mathcal{G}_{t}\cap\mathcal{S}_{k}|$ be the number of users who
have transmitted in the $t$th RE and have not been decoded in the first $k-1$
iterations, $\mathbf{H}_{t}^{k}\in\mathbb{C}^{N\times M_{t}^{k}}$ denote the
channel matrix which contains the channels of the $M_{t}^{k}$ users,
$\mathbf{P}_{t}^{k}\in\mathbb{C}^{\tau\times M_{t}^{k}}$ denote a matrix that
contains the pilot sequences of the $M_{t}^{k}$ users and
$\mathbf{B}_{t}^{k}\triangleq\sigma_{h}^{2}\text{diag}(\beta_{i_{1}},\beta_{i_{2}},\ldots,\beta_{i_{M_{t}^{k}}})$
be a diagonal matrix containing the path loss coefficients of the $M_{t}^{k}$
users, with
$\mathcal{G}_{t}\cap\mathcal{S}_{k}=\\{i_{1},i_{2},\ldots,i_{M_{t}^{k}}\\}$.
Thus, the received signal from (3) can be written as
$\mathbf{Y}_{t}^{pk}=\mathbf{H}_{t}^{k}\mathbf{P}_{t}^{kH}+{\mathbf{N}}_{t}^{p}$,
where $\mathbf{P}_{t}^{k}=[\mathbf{P}]_{\mathrel{\mathop{\mathchar
58\relax}},\mathcal{G}_{t}\cap\mathcal{S}_{k}}$. The MMSE estimate
$\hat{\mathbf{H}}_{t}^{k}$ of $\mathbf{H}_{t}^{k}$ is presented in Theorem 1,
and can be written as

$\displaystyle\hat{\mathbf{H}}_{t}^{k}$
$\displaystyle={\mathbf{Y}}_{t}^{pk}(\mathbf{P}_{t}^{k}\mathbf{B}_{t}^{k}\mathbf{P}_{t}^{kH}+N_{0}\mathbf{I}_{\tau})^{-1}\mathbf{P}_{t}^{k}\mathbf{B}_{t}^{k},$
(7a)
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}{\mathbf{Y}}_{t}^{pk}\mathbf{P}_{t}^{k}\mathbf{B}_{t}^{k}(\mathbf{P}_{t}^{kH}\mathbf{P}_{t}^{k}\mathbf{B}_{t}^{k}+N_{0}\mathbf{I}_{M_{t}^{k}})^{-1},$
(7b)

where $(a)$ follows from $(\mathbf{AB}+\mathbf{I})^{-1}\mathbf{A}$ $=$
$\mathbf{A}(\mathbf{BA}+\mathbf{I})^{-1}$. Here, the estimate can be
calculated via an inverse of either a $\tau\times\tau$ matrix or an
$M_{t}^{k}\times M_{t}^{k}$ matrix as required. The MSBL estimate converges to
the MMSE estimate when the hyperparameters are estimated well enough, as will
be seen in Sec. VI.

###  III-C Low Complexity MMSE with Known APM

We now describe a low complexity MMSE (LCMMSE) estimate that does not require
a matrix inversion computation. For this purpose, the received signal in (3)
is right-multiplied by the pilot $\mathbf{p}_{m}$ to obtain

$\displaystyle\mathbf{y}_{tm}^{pk}$
$\displaystyle={\mathbf{Y}_{t}^{pk}}\mathbf{p}_{m}=\sum\nolimits_{i\in\mathcal{S}_{k}}g_{ti}\mathbf{h}_{ti}(\mathbf{p}_{i}^{H}\mathbf{p}_{m})+{\mathbf{N}_{t}^{p}}\mathbf{p}_{m},$
(8)

which is used to find an MMSE estimate of the channel ${\mathbf{h}}_{tm}$ of
the $m$th user in the $t$th RE. The LCMMSE channel estimate
$\hat{\mathbf{h}}_{tm}^{k}$ is described in Theorem 1 and is calculated as

$\displaystyle\hat{\mathbf{h}}_{tm}^{k}\\!$
$\displaystyle=\\!\dfrac{g_{tm}\beta_{m}\|\mathbf{p}_{m}\|^{2}\sigma_{h}^{2}}{N_{0}\|\mathbf{p}_{m}\|^{2}+\sum\nolimits_{i\in\mathcal{S}_{k}}|\mathbf{p}_{i}^{H}\mathbf{p}_{m}|^{2}g_{ti}\beta_{i}\sigma_{h}^{2}}\mathbf{y}_{tm}^{pk}\triangleq\eta_{tm}^{k}\mathbf{y}_{tm}^{pk}.$

Similar to the MMSE estimate, the LCMMSE estimate uses the knowledge of the
APM and $\\{\beta_{i}\\}_{i=1}^{M}$. While the MMSE estimator uses the signal
${\mathbf{Y}_{t}^{pk}}$ to compute the estimates, and thus exploits all the
information available at the BS, the LCMMSE estimator uses only
$\mathbf{y}_{tm}^{pk}$, i.e., the projection of ${\mathbf{Y}_{t}^{pk}}$ onto
${\mathbf{p}}_{m}$, to estimate ${\mathbf{h}}_{tm}$.

The channel estimates under the three schemes and their error variances are
given by the following theorem.

Table I: Channel estimates and error variances under three estimation schemes.

| Sparsity-based estimation with MSBL | MMSE | LCMMSE
---|---|---|---
$\hat{\mathbf{H}}^{k}_{t}$ | ${{\mathbf{Y}}^{pk}_{t}}\mathbf{P}^{k}{\hat{\bm{\Gamma}}_{kt}}(\mathbf{P}^{kH}\mathbf{P}^{k}{\hat{\bm{\Gamma}}_{kt}}+N_{0}\mathbf{I}_{M^{k}})^{-1}$ | ${{\mathbf{Y}}^{pk}_{t}}\mathbf{P}_{t}^{k}\mathbf{B}_{t}^{k}(\mathbf{P}_{t}^{kH}\mathbf{P}_{t}^{k}\mathbf{B}_{t}^{k}+N_{0}\mathbf{I}_{M_{t}^{k}})^{-1}$ | ${{\mathbf{Y}}^{pk}_{t}}\mathbf{P}_{t}^{k}\text{ diag}(\eta_{ti_{1}}^{k},\ldots,\eta_{ti_{M_{t}^{k}}}^{k})$
$\delta_{ti}^{k}$ | $\beta_{i}\sigma_{h}^{2}\left(\frac{N_{0}\|\mathbf{c}_{ti}^{k}\|^{2}+\sum\nolimits_{j\in\mathcal{S}_{k}^{i}}|r_{jti}^{k}|^{2}\hat{g}_{tj}^{k}g_{tj}\beta_{j}\sigma_{h}^{2}}{N_{0}\|\mathbf{c}_{ti}^{k}\|^{2}+\sum\nolimits_{j\in\mathcal{S}_{k}}|r_{jti}^{k}|^{2}\hat{g}_{tj}^{k}g_{tj}\beta_{j}\sigma_{h}^{2}}\right)$ | $\beta_{i}\sigma_{h}^{2}\left(\frac{N_{0}\|\mathbf{c}_{ti}^{k}\|^{2}+\sum\nolimits_{j\in\mathcal{S}_{k}^{i}}|r_{jti}^{k}|^{2}g_{tj}\beta_{j}\sigma_{h}^{2}}{N_{0}\|\mathbf{c}_{ti}^{k}\|^{2}+\sum\nolimits_{j\in\mathcal{S}_{k}}|r_{jti}^{k}|^{2}g_{tj}\beta_{j}\sigma_{h}^{2}}\right)$ | $\beta_{i}\sigma_{h}^{2}\left(\frac{N_{0}\|\mathbf{p}_{i}\|^{2}+\sum\nolimits_{j\in\mathcal{S}_{k}^{i}}|\mathbf{p}_{j}^{H}\mathbf{p}_{i}|^{2}g_{tj}\beta_{j}\sigma_{h}^{2}}{N_{0}\|\mathbf{p}_{i}\|^{2}+\sum\nolimits_{j\in\mathcal{S}_{k}}|\mathbf{p}_{j}^{H}\mathbf{p}_{i}|^{2}g_{tj}\beta_{j}\sigma_{h}^{2}}\right)$

######  Theorem 1.

The channel estimate $\hat{\mathbf{H}}^{k}_{t}$ of ${\mathbf{H}}^{k}_{t}$ in
the $t$th RE in the $k$th decoding iteration, under the three estimation
schemes, namely MSBL, MMSE, and LCMMSE, is given in Table I. Specifically, the
estimate of the channel ${\mathbf{h}}_{ti}$ of the $i$th user is calculated as
$\hat{\mathbf{h}}_{ti}^{k}=[\hat{\mathbf{H}}_{t}^{k}]_{\mathrel{\mathop{\mathchar
58\relax}},i}$. Further, the covariance of the estimation error
$\tilde{\mathbf{h}}_{ti}^{k}\triangleq\hat{\mathbf{h}}_{ti}^{k}-{\mathbf{h}}_{ti}$
is $\delta_{ti}^{k}\mathbf{I}_{N}$, where $\delta_{ti}^{k}$ is listed in Table
I, with $\mathbf{c}_{ti}^{k}=[\mathbf{C}_{t}^{k}]_{\mathrel{\mathop{\mathchar
58\relax}},i}$ and
$r_{jti}^{k}\triangleq\mathbf{p}_{j}^{H}\mathbf{c}_{ti}^{k}$. For MSBL, we
have
$\mathbf{C}_{t}^{k}\triangleq\mathbf{P}^{k}\mathbf{D}_{t}^{k}(\mathbf{P}^{kH}\mathbf{P}^{k}\mathbf{D}_{t}^{k}+N_{0}\mathbf{I}_{M^{k}})^{-1}$,
where $\mathbf{D}_{t}^{k}\triangleq$
diag$(d_{ti_{1}}^{k},d_{ti_{2}}^{k},\ldots,d_{ti_{M^{k}}}^{k})$ with
$d_{ti}^{k}=\hat{g}_{ti}^{k}g_{ti}\beta_{i}\sigma_{h}^{2}$. For MMSE, we have
$\mathbf{C}_{t}^{k}\triangleq\mathbf{P}_{t}^{k}\mathbf{B}_{t}^{k}(\mathbf{P}_{t}^{kH}\mathbf{P}_{t}^{k}\mathbf{B}_{t}^{k}+N_{0}\mathbf{I}_{M_{t}^{k}})^{-1}$.

###### Proof.

See Appendix A. ∎

_Remark:_ The LCMMSE estimate is composed of two components: a scaling
coefficient $\eta_{tm}^{k}$ and the post-combined received pilot signal
$\mathbf{y}_{tm}^{pk}$. From (8), we see that the received pilot signal
$\mathbf{y}_{tm}^{pk}$ contains pilots of other users, if pilot sequences are
not orthogonal. With orthogonal pilots,
$\mathbf{p}_{i}^{H}\mathbf{p}_{m}=0,\forall i\neq m$, the LCMMSE estimate is

$\hat{\mathbf{h}}_{tm}^{k}=\dfrac{g_{tm}\beta_{m}\sigma_{h}^{2}}{N_{0}+\|\mathbf{p}_{m}\|^{2}g_{tm}\beta_{m}\sigma_{h}^{2}}\left(g_{tm}\mathbf{h}_{tm}\|\mathbf{p}_{m}\|^{2}+{\mathbf{N}_{t}^{p}}\mathbf{p}_{m}\right),$

and
$\delta_{tm}^{k}=g_{tm}\beta_{m}\sigma_{h}^{2}N_{0}/(N_{0}+g_{tm}\beta_{m}\sigma_{h}^{2}\|\mathbf{p}_{m}\|^{2})$,
i.e., there is no pilot contamination, and the LCMMSE estimate coincides with
the MMSE estimate.

_Complexity:_ The MMSE scheme has a complexity of
$\mathcal{O}(\tau^{2}M_{t}^{k})$ floating point operations (flops) since it
involves inverting a $\tau\times\tau$ matrix. The MSBL scheme, with $s$
iterations, has a complexity of $\mathcal{O}(s\tau^{2}M^{k})$ flops [19]. The
LCMMSE scheme has the lowest complexity of $\mathcal{O}(M_{t}^{k})$ flops
since it does not need any matrix inversion.

##  IV SINR Analysis

In this section, the SINR of each user in all the REs where it transmits data
is derived, accounting for pilot contamination and channel estimation errors.
Let $\rho_{tm}^{k}$ denote the SINR of the $m$th user in the $t$th RE in the
$k$th decoding iteration. Similar to (2), the received data signal in the
$t$th RE and $k$th decoding iteration can be written as

$\displaystyle{\mathbf{y}_{t}^{k}}$
$\displaystyle=\sum\nolimits_{i\in\mathcal{S}_{k}}g_{ti}\mathbf{h}_{ti}x_{i}+{\mathbf{n}_{t}}.$
(9)

A combining vector $\mathbf{a}_{tm}^{k}$ is used to decode the $m$th user in
the $t$th RE and $k$th decoding iteration, and thus we obtain

$\displaystyle\tilde{y}_{tm}^{k}=\mathbf{a}_{tm}^{kH}{\mathbf{y}_{t}^{k}}$
$\displaystyle={\mathbf{a}}_{tm}^{kH}\hat{\mathbf{h}}_{tm}^{k}g_{tm}x_{m}-{\mathbf{a}}_{tm}^{kH}\tilde{\mathbf{h}}_{tm}^{k}g_{tm}x_{m}$
$\displaystyle\qquad+{\mathbf{a}}_{tm}^{kH}\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}g_{ti}\mathbf{h}_{ti}x_{i}+{\mathbf{a}}_{tm}^{kH}\mathbf{n}_{t},$
(10)

where $\tilde{\mathbf{h}}_{tm}^{k}$ is as defined in Theorem 1. From the
above, we see that the signal used to decode the $m$th user’s data is composed
of four terms. The term
$T_{1}\triangleq{\mathbf{a}}_{tm}^{kH}\hat{\mathbf{h}}_{tm}^{k}g_{tm}x_{m}$ is
the useful signal component of the $m$th user; the term
$T_{2}\triangleq{\mathbf{a}}_{tm}^{kH}\tilde{\mathbf{h}}_{tm}^{k}g_{tm}x_{m}$
is contributed by the channel estimation error $\tilde{\mathbf{h}}_{tm}^{k}$
of the $m$th user; the term
$T_{3}\triangleq\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}{\mathbf{a}}_{tm}^{kH}\mathbf{h}_{ti}g_{ti}x_{i}$
captures the inter-user interference from the users who have also transmitted
in the $t$th RE and have not yet been decoded upto the $k$th decoding
iteration; and the last term
$T_{4}\triangleq{\mathbf{a}}_{tm}^{kH}\mathbf{n}_{t}$ is the additive noise
component.

In order to compute the SINR, the power in the received signal is calculated
conditioned on the knowledge of the estimates [20]. Since MMSE estimates are
employed, all three estimates are uncorrelated with the channel estimation
error, and thus $T_{2}$ is uncorrelated with $T_{1}$. The additive noise is
uncorrelated with the signal, and since the users’ data signals are
independent, $T_{3}$ is uncorrelated with the other terms. Thus, all four
components in the received signal are uncorrelated and the total power is the
sum of the powers of the individual components

$\displaystyle\mathbb{E}_{\mathbf{z}}[|\tilde{y}_{tm}^{k}|^{2}]$
$\displaystyle=\sum\nolimits_{i=1}^{4}\mathbb{E}_{\mathbf{z}}[|T_{i}|^{2}],$
(11)

where $\mathbf{z}$ contains the channel estimates of the users. The SINR for
all the users is now presented.

######  Theorem 2.

The signal to interference plus noise ratio (SINR) achieved by the $m$th user
in the $t$th RE in the $k$th decoding iteration can be written as

$\displaystyle\rho_{tm}^{k}$
$\displaystyle=\dfrac{{\tt{Gain}}_{tm}^{k}}{N_{0}/P+{\tt{MUI}}_{tm}^{k}+{\tt{Est}}_{tm}^{k}},\
\forall m\in\mathcal{S}_{k},$ (12)

where ${\tt{Gain}}_{tm}^{k}$ represents the useful signal power of the $m$th
user, ${\tt{MUI}}_{tm}^{k}$ represents the multi-user interference power of
other users, and ${\tt{Est}}_{tm}^{k}$ represents the interference power
caused due to the channel estimation errors. Under MMSE and LCMMSE channel
estimation, these can be expressed as

$\displaystyle{\tt{Gain}}_{tm}^{k}$
$\displaystyle=g_{tm}\dfrac{|{\mathbf{a}}_{tm}^{kH}\hat{\mathbf{h}}_{tm}^{k}|^{2}}{\|{\mathbf{a}}_{tm}^{k}\|^{2}},\
{\tt{MUI}}_{tm}^{k}=\textstyle{\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}}g_{ti}\dfrac{|{\mathbf{a}}_{tm}^{kH}\hat{\mathbf{h}}_{ti}^{k}|^{2}}{\|{\mathbf{a}}_{tm}^{k}\|^{2}},$
$\displaystyle{\tt{Est}}_{tm}^{k}$
$\displaystyle=\textstyle{\sum\nolimits_{i\in\mathcal{S}_{k}}}g_{ti}\delta_{ti}^{k}.$

With the sparsity-based scheme, the SINR denominator contains an additional
term, ${\tt{FNU}}_{tm}^{k}$, which represents the interference power caused
due to false negative users. The corresponding terms with MSBL can be
expressed as

$\displaystyle{\tt{Gain}}_{tm}^{k}\\!\\!$
$\displaystyle=\hat{g}_{tm}^{k}g_{tm}\frac{|{\mathbf{a}}_{tm}^{kH}\hat{\mathbf{h}}_{tm}^{k}|^{2}}{\|{\mathbf{a}}_{tm}^{k}\|^{2}},\
{\tt{MUI}}_{tm}^{k}\\!=\\!\\!\\!\textstyle{\sum\limits_{i\in\mathcal{S}_{k}^{m}}}\hat{g}_{ti}^{k}g_{ti}\dfrac{|{\mathbf{a}}_{tm}^{kH}\hat{\mathbf{h}}_{ti}^{k}|^{2}}{\|{\mathbf{a}}_{tm}^{k}\|^{2}},$
$\displaystyle{\tt{Est}}_{tm}^{k}\\!\\!$
$\displaystyle=\\!\textstyle{\sum\nolimits_{i\in\mathcal{S}_{k}}}\hat{g}_{ti}^{k}g_{ti}\delta_{ti}^{k},\
{\tt{FNU}}_{tm}^{k}\\!=\\!\textstyle{\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}}(1-\hat{g}_{ti}^{k})g_{ti}\beta_{i}\sigma_{h}^{2}.$

Here, the estimates
$\hat{\mathbf{h}}_{ti}^{k}=[\hat{\mathbf{H}}_{t}^{k}]_{\mathrel{\mathop{\mathchar
58\relax}},i}$ and the error variances $\delta_{ti}^{k}$ are obtained from
Theorem 1 for all the three schemes.

###### Proof.

See Appendix B. ∎

The SINR expression derived in Theorem 2 is applicable to any arbitrary
receive combining scheme given by the matrix $\mathbf{A}_{t}^{k}$, with
$\mathbf{a}_{tm}^{k}=[\mathbf{A}_{t}^{k}]_{\mathrel{\mathop{\mathchar
58\relax}},m}$. When regularized zero forcing (RZF) combining is used, the
combining matrix is

$\displaystyle{\mathbf{A}_{t}^{k}}={\hat{\mathbf{H}}^{k}_{t}}({\hat{\mathbf{H}}^{kH}_{t}}{\hat{\mathbf{H}}^{k}_{t}}+\lambda\mathbf{I}_{M_{t}^{k}})^{-1},$
(13)

where $\lambda$ is the regularization parameter. The SINR with RZF can be
computed by substituting the columns of the above matrix into (12). We now
describe two popular combining schemes, which are special cases of RZF, in
which simpler expressions for the SINR can be computed. The expressions are
written for MMSE/LCMMSE, and can be extended to MSBL as detailed in Theorem 2.

Input :
$\tau,N,T,M,N_{0},\mathbf{G},\mathbf{P},\\{{\mathbf{Y}}_{t}^{p}\\}_{t=1}^{T},k_{\max}$

1 Initialize: $\mathcal{S}_{1}=[M]$, $\mathcal{G}_{t}=\\{i\in[M]|g_{ti}=1\\}$

2 for  _$k=1,2,3,\ldots,k_{\max}$_ do

3 for  _$t=1,2,\ldots,T$_ do

4 Find $M_{t}^{k}=|\mathcal{G}_{t}\cap\mathcal{S}_{k}|,$
$\mathbf{P}_{t}^{k}=[\mathbf{P}]_{\mathrel{\mathop{\mathchar
58\relax}},\mathcal{G}_{t}\cap\mathcal{S}_{k}},\mathbf{Y}_{t}^{pk}$

5 Compute $\hat{\mathbf{h}}_{ti}^{k},\ \forall i\in\mathcal{S}_{k}$ via
Theorem 1

6 Evaluate the SINR $\rho_{ti}^{k}$ via Theorem 2

7 If $\rho_{ti}^{k}\geq\gamma_{\text{th}}$, remove user $i$ from
$\mathcal{S}_{k}$ and perform IC in all REs where $g_{ti}=1$

8

9 end for

10

11 end for

Output: $\mathcal{T}=(M-|\mathcal{S}_{k_{\max}}|)/T,\
\mathsf{PLR}=|\mathcal{S}_{k_{\max}}|/M$

Algorithm 2 Performance Evaluation of IRSA

####  IV-1 Maximal Ratio Combining (MRC)

MRC is obtained from RZF as $\lambda\rightarrow\infty$ and the combining
matrix becomes ${\mathbf{A}_{t}^{k}}={\hat{\mathbf{H}}^{k}_{t}}$. Thus
${\mathbf{a}}_{tm}^{k}=\hat{\mathbf{h}}_{tm}^{k}$, and SINR can be computed as

$\rho_{tm}^{k}=\dfrac{Pg_{tm}\|\hat{\mathbf{h}}_{tm}^{k}\|^{2}}{N_{0}+\sum\nolimits_{i\in\mathcal{S}_{k}}Pg_{ti}\delta_{ti}^{k}+\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}Pg_{ti}\frac{|\hat{\mathbf{h}}_{tm}^{kH}\hat{\mathbf{h}}_{ti}^{k}|^{2}}{\|\hat{\mathbf{h}}_{tm}^{k}\|^{2}}}.$
(14)

####  IV-2 Zero Forcing (ZF)

The RZF combiner reduces to the ZF combiner as $\lambda\rightarrow 0$. The
inverse of the gram-matrix of the channel estimates exists with probability
one when $N\geq M_{t}^{k}$ and ${\hat{\mathbf{H}}^{k}_{t}}$ has full column
rank.222 We note that the condition $N\geq M_{t}^{k}$ is not hard to satisfy
in IRSA. For example, with $L=2,3,4$, each RE will be occupied by $6,9,12$
users on an average, respectively, if the average repetition factor is
$\bar{d}=3$. Thus any $N$ greater than, say, $16$ would be sufficient to
decode the users in most REs. Hence, we can compute the combining matrix as
${\mathbf{A}_{t}^{k}}={\hat{\mathbf{H}}^{k}_{t}}({\hat{\mathbf{H}}^{kH}_{t}}{\hat{\mathbf{H}}^{k}_{t}})^{-1}.$
Using the above, it is easy to show that the SINR expression simplifies as
[20]

$\displaystyle\rho_{tm}^{k}=\dfrac{Pg_{tm}}{(N_{0}+\sum\nolimits_{i\in\mathcal{S}_{k}}Pg_{ti}\delta_{ti}^{k})[({\hat{\mathbf{H}}^{kH}_{t}}{\hat{\mathbf{H}}^{k}_{t}})^{-1}]_{mm}}.$
(15)

Note that the third term in the denominator of (14) has been suppressed with
ZF combining. However, due to pilot contamination, the term
$[({\hat{\mathbf{H}}^{kH}_{t}}{\hat{\mathbf{H}}^{k}_{t}})^{-1}]_{mm}$ may
contain contributions from the channels of all users. As a consequence, the
gram matrix could be ill-conditioned, and the denominator term could be large.
Thus, the pilot length, which determines the pilot contamination incurred, is
crucial in comparing the performance obtained by the combining schemes. The
system throughput can now be calculated from the above SINR expressions via
the decoding model described in Sec. II-A, and is described in Algorithm 2 for
MMSE/LCMMSE. For MSBL, the initial step in each RE instead consists of finding
$M^{k}=|\mathcal{S}_{k}|,$ and
$\mathbf{P}^{k}=[\mathbf{P}]_{\mathrel{\mathop{\mathchar
58\relax}},\mathcal{S}_{k}}$. We also estimate $\\{g_{ti}\\}$ and
$\\{\mathbf{h}_{ti}\\}$ via Algorithm 1 before finding the SINR.

Before proceeding with the analysis of the throughput, we briefly discuss the
SINR in the massive MIMO regime, which helps us in interpreting the SINR
expressions. We note that the results presented in Sec. VI hold true for any
$N$. However, when $N$ is large, a simpler expression for SINR with MRC can be
obtained as follows.

######  Lemma 1.

As the number of antennas $N$ gets large, the SINR with MRC converges almost
surely to

$\displaystyle\overline{\rho}_{tm}^{k}$ $\displaystyle=\dfrac{N{\tt
Sig}_{tm}^{k}}{\epsilon_{tm}^{k}\left(N_{0}/P+{\tt NCoh}_{tm}^{k}\right)+{\tt
Coh}_{tm}^{k}},$ (16)

where ${\tt Sig}_{tm}^{k}$ is the desired signal, ${\tt NCoh}_{tm}^{k}$
represents the non-coherent interference, and ${\tt Coh}_{tm}^{k}$ represents
the coherent interference. Each of these can be found in Table II. Here,
$\delta_{tm}^{k}$ and ${\mathbf{c}}_{tm}^{k}$ are obtained from Theorems 1 and
2, respectively, for the three estimation schemes.

###### Proof.

See Appendix C. ∎

_Remark:_ ${\tt NCoh}_{tm}^{k}$ arises due to channel estimation errors and is
independent of $N$, while ${\tt Coh}_{tm}^{k}$ is due to pilot contamination
and increases linearly with $N$. Further, since $\overline{\rho}_{tm}^{k}$ is
independent of the fading states of each user, it assures successful recovery
of packets with high probability if
$\overline{\rho}_{tm}^{k}\gg\gamma_{\text{th}}$. Similarly, the packet will
not be decodable with probability close to $1$ if
$\overline{\rho}_{tm}^{k}\ll\gamma_{\text{th}}$. However, it turns out that in
order to characterize the throughput of IRSA, it is necessary to capture the
statistics of the SINR when
$\overline{\rho}_{tm}^{k}\approx\gamma_{\text{th}}$. The small fluctuations in
$\rho_{tm}^{k}$ around $\overline{\rho}_{tm}^{k}$ due to fading, and the
resulting probability of packet decoding error, need to be calculated
accurately. Hence, the calculation of the statistics of the SINR using (12) is
vital to find the throughput of IRSA. We address this in the next section.

##  V Theoretical Analysis of Throughput

Density Evolution (DE) analysis has been applied to characterize the
asymptotic performance of message passing-based decoding on graphs for low
density parity check codes [21] and IRSA [6]. In this section, the
representation of IRSA decoding as a bipartite graph is discussed first. Then
the graph perspective distributions are defined, the failure probabilities are
derived, and finally, the asymptotic throughput of IRSA is characterized. It
is assumed that users perform path loss inversion-based power control. We note
that a closed form expression for the throughput cannot be derived even for
the most basic variant of IRSA due to the underlying graph structure [6].
Hence, we need to resort to DE, which provides an iterative recipe to compute
the throughput.

SIC-based decoding can be viewed as message passing on a bipartite graph [6],
and thus IRSA, which uses SIC decoding, can be decoded on graphs. A typical
IRSA frame can be represented as a bipartite graph, which is made up of $M$
user nodes (one node for each user), $T$ RE nodes (one node for each RE), and
the edges between them. An edge connects a user node to an RE node if and only
if that user has transmitted a packet in that corresponding RE. During
decoding, edges that connect to users whose SINR is above a threshold are
removed from each RE. Each decoding iteration consists of several intra-RE SIC
and inter-RE SIC steps. Once an SIC step is performed, the corresponding edge
in the bipartite graph is removed. Decoding is successful if, at the end of
the SIC process, all edges in the graph get removed. A decoding failure is
declared if not all edges have been removed or no new edge is removed from the
graph in two consecutive iterations.

Table II: Deterministic equivalent approximation to the SINR.

| Sparsity-based estimation with MSBL | MMSE | LCMMSE
---|---|---|---
$\epsilon_{tm}^{k}$ | $N_{0}\|{\mathbf{c}}_{tm}^{k}\|^{2}+\sum\nolimits_{i\in\mathcal{S}_{k}}g_{ti}\beta_{i}\sigma_{h}^{2}|\mathbf{c}_{tm}^{kH}\mathbf{p}_{i}|^{2}$ | $N_{0}\|{\mathbf{c}}_{tm}^{k}\|^{2}+\sum\nolimits_{i\in\mathcal{S}_{k}}g_{ti}\beta_{i}\sigma_{h}^{2}|\mathbf{c}_{tm}^{kH}\mathbf{p}_{i}|^{2}$ | $N_{0}\|{\mathbf{p}}_{m}\|^{2}+\sum\nolimits_{i\in\mathcal{S}_{k}}g_{ti}\beta_{i}\sigma_{h}^{2}|\mathbf{p}_{m}^{H}\mathbf{p}_{i}|^{2}$
${\tt Sig}_{tm}^{k}$ | $\hat{g}_{tm}^{k}g_{tm}(\epsilon_{tm}^{k})^{2}$ | $g_{tm}(\epsilon_{tm}^{k})^{2}$ | $g_{tm}\beta_{m}^{2}\sigma_{h}^{4}\|{\mathbf{p}}_{m}\|^{4}$
${\tt NCoh}_{tm}^{k}$ | $\hat{g}_{tm}^{k}g_{tm}\delta_{tm}^{k}+\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}g_{ti}\beta_{i}\sigma_{h}^{2}$ | $g_{tm}\delta_{tm}^{k}+\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}g_{ti}\beta_{i}\sigma_{h}^{2}$ | $g_{tm}\delta_{tm}^{k}+\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}g_{ti}\beta_{i}\sigma_{h}^{2}$
${\tt Coh}_{tm}^{k}$ | $N\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}g_{ti}\beta_{i}^{2}\sigma_{h}^{4}|\mathbf{c}_{tm}^{kH}\mathbf{p}_{i}|^{2}$ | $N\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}g_{ti}\beta_{i}^{2}\sigma_{h}^{4}|\mathbf{c}_{tm}^{kH}\mathbf{p}_{i}|^{2}$ | $N\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}g_{ti}\beta_{i}^{2}\sigma_{h}^{4}|\mathbf{p}_{m}^{H}\mathbf{p}_{i}|^{2}$

###  V-A Graph Perspective Degree Distributions

The total number of packets transmitted by a user in a given frame is referred
to as the repetition factor of that user. It is equal to the degree of the
user node at the start of decoding, and is the same as the number of edges
connected to that user node in the bipartite graph representation of SIC
decoding. The _node-perspective user degree distribution_ is defined as the
set of probabilities $\\{\phi_{d}\\}_{d=2}^{d_{\max}}$, where $\phi_{d}$
represents the probability that a user has a repetition factor $d$ with
$d_{\max}$ being the maximum number of REs in which any user is allowed to
transmit. Here, $\phi_{d}$ is nonzero for $d\geq 2$ since each user transmits
at least $2$ packets in IRSA.

The total number of packets received in an RE is referred to as the collision
factor of that RE. It is equal to the degree of the RE node at the start of
decoding, which is the number of edges connected to that RE node. The _node-
perspective RE degree distribution_ is defined as the set of probabilities
$\\{\psi_{c}\\}_{c=0}^{M}$, where $\psi_{c}$ represents the probability that
an RE has a collision factor $c$. The polynomial representations of the node-
perspective user and RE degree distributions are

$\displaystyle\phi(x)=\sum\nolimits_{d=2}^{d_{\max}}\phi_{d}x^{d},\ \
\psi(x)=\sum\nolimits_{c=0}^{M}\psi_{c}x^{c},$ (17)

respectively. The corresponding _edge-perspective user and RE degree
distributions_ are defined as
$\lambda(x)=\sum_{d=2}^{d_{\max}}\lambda_{d}x^{d-1}$
$=\phi^{\prime}(x)/\phi^{\prime}(1);$ $\xi(x)=\sum_{c=1}^{M}\xi_{c}x^{c-1}$
$=\psi^{\prime}(x)/\psi^{\prime}(1)$, respectively, where
$\lambda_{d}=d\phi_{d}/\phi^{\prime}(1)$ represents the probability that an
edge is connected to a user with repetition factor $d$ and
$\xi_{c}=c\psi_{c}/\psi^{\prime}(1)$ represents the probability that an edge
is connected to an RE with collision factor $c$.

The _input load_ $L$ of the system is defined as the ratio of the number of
users to the number of REs, $L\triangleq M/T$. The average repetition factor
is $\bar{d}=\phi^{\prime}(1)=\sum_{d}d\phi_{d}$ and the average collision
factor is $\bar{c}=\psi^{\prime}(1)=\sum_{c}c\psi_{c}$, making the load
$L={M}/T=\bar{c}/\bar{d}$. Since $\bar{c}=L\bar{d}$, fixing the load and the
node-perspective user degree distribution fixes the other three degree
distributions as well. The probability that a generic user, from a total of
${M}$ users, transmits within an RE is $\bar{c}/{M}$. Since the users transmit
their packets independently of each other, $\psi_{c}$ follows a binomial
distribution. Thus, the coefficients of the polynomials representing the node
and edge-perspective RE degree distributions are respectively given by

$\displaystyle\psi_{c}$
$\displaystyle=\binom{M}{c}\left(\dfrac{\bar{c}}{{M}}\right)^{c}\left(1-\dfrac{\bar{c}}{{M}}\right)^{{M}-c},$
(18a) $\displaystyle\text{and \ }\xi_{c}$
$\displaystyle=\binom{{M}-1}{c-1}\left(\dfrac{\bar{c}}{{M}}\right)^{c-1}\left(1-\dfrac{\bar{c}}{{M}}\right)^{{M}-c}.$
(18b)

For a fixed $L={M}/T$, as ${M},T\rightarrow\infty$, the node-perspective and
edge-perspective RE degree distributions, which are binomial, become Poisson
distributed [22]:

$\psi_{c}=\frac{\left(\bar{c}\right)^{c}\exp\left(-\bar{c}\right)}{c!}\text{
and
}\xi_{c}=\frac{\left(\bar{c}\right)^{c-1}\exp\left(-\bar{c}\right)}{\left(c-1\right)!}.$
(19)

###  V-B Failure Probabilities

In the case of a decoding failure, failure messages are exchanged along the
edges between the user and the RE nodes. The probability that an edge carries
a failure message from an RE node to a user node in the $i$th iteration is
denoted by $p_{i}$. The probability that an edge carries a failure message
from a user node to an RE node in the $i$th iteration is denoted by $q_{i}$.

The failure probability $q_{i}$ is calculated using the edge-perspective user
degree distribution as

$\displaystyle
q_{i}=\sum\nolimits_{d=2}^{d_{\max}}\lambda_{d}q_{i}^{(d)}=\sum\nolimits_{d=2}^{d_{\max}}\lambda_{d}p_{i-1}^{d-1}=\lambda(p_{i-1}).$
(20)

Here, $q_{i}^{(d)}$ is the probability that an edge carries a failure message
in the $i$th iteration given that it is connected to a user node with
repetition factor $d$. The edges carry a failure message from a user if and
only if all the other $d-1$ incoming edges to that user carry failure messages
in the previous iteration, i.e., $q_{i}^{(d)}=p_{i-1}^{d-1}$.

The failure probability $p_{i}$ is calculated using the edge-perspective RE
degree distribution as

$\displaystyle p_{i}$
$\displaystyle=\sum\nolimits_{c=1}^{M}\xi_{c}p_{i}^{(c)}\
{\xrightarrow{{M}\rightarrow\infty}}\ p_{i}=\
\sum\nolimits_{c=1}^{\infty}\xi_{c}p_{i}^{(c)},$ (21)

where $p_{i}^{(c)}$ is the probability that an edge carries a failure message
in the $i$th iteration given that it is connected to an RE node with collision
factor $c$. DE is applicable as ${M}$ and $T\rightarrow\infty$ with $L={M}/T$
kept fixed [6]. Hence the above probability is computed as an infinite
summation.

In the SINR-threshold model, decoding failure happens at an RE node if the
SINR of all users who have transmitted in that RE and have not yet been
decoded is below the SINR threshold. This constitutes a failure message from
the RE node [11]. In order to determine $p_{i}^{(c)}$, any one of the $c$
packets is considered to be a reference packet, which can get decoded with a
combination of intra-RE and inter-RE SIC. Separating the intra-RE and inter-RE
SIC, $p_{i}^{(c)}$ can be evaluated as

$\displaystyle p_{i}^{(c)}$
$\displaystyle=1-\sum\limits_{r=1}^{c}\theta_{r}\binom{c-1}{r-1}q_{i}^{r-1}(1-q_{i})^{c-r}.$
(22)

Here, $\theta_{r}$ denotes the probability that the reference packet gets
decoded in the current decoding iteration starting from degree $r$ using only
intra-RE SIC, and $\binom{c-1}{r-1}q_{i}^{r-1}(1-q_{i})^{c-r}$ denotes the
probability that the collision factor of the RE node reduces from $c$ to $r$
using only inter-RE SIC [9]. Substituting for $p_{i}^{(c)}$ from (22), we
obtain $p_{i}$ as a function of $q_{i}$:

$\displaystyle p_{i}$
$\displaystyle=1-\sum_{c=1}^{\infty}\sum_{r=1}^{c}\xi_{c}\theta_{r}\binom{c-1}{r-1}q_{i}^{r-1}(1-q_{i})^{c-r}.$
(23)

Thus, we need to evaluate $\theta_{r}$, which we describe in the next
subsection. We compute $p_{i}$ and $q_{i}$ recursively from each other as
observed in (23) and (20). Finally, we will compute the packet loss rate and
the throughput using the value of $p_{\infty}$ obtained from the above
recursive equations.

###  V-C Characterization of $\theta_{r}$

We now describe a procedure to evaluate the success probability $\theta_{r}$,
which is the probability of decoding the reference packet in an RE with degree
$r$ via intra-RE SIC only. There are $r$ users whose packets have not yet been
decoded in the given RE. The reference packet can get decoded in any of the
intra-RE SIC steps. The packets with SINR higher than that of the reference
packet get decoded first. Further, the reference packet can only be decoded if
decoding has been successful for higher SINR packets, i.e., if they satisfied
SINR $\geq\gamma_{\text{th}}$ as well. Thus, $\theta_{r}$ is the joint
probability that the reference packet and the packets with higher SINRs all
get decoded.

Clazzer et al. [11] evaluate $\theta_{r}$ as the probability "$D(r)$" under a
Rayleigh fading SISO channel setup with a perfect CSI assumption. The same
method cannot be applied here, since we consider MIMO Rayleigh fading and
account for imperfect CSI due to pilot contamination and channel estimation
errors. In particular, in a MIMO setup, it is possible that multiple users’
SINR simultaneously exceed the decoding threshold. Further, their work is
limited to the case where the decoding threshold $\gamma_{\text{th}}$ is such
that only one user can be decoded in any decoding iteration, while we make no
such assumptions.

Since $\theta_{r}$ is evaluated based on the SINR of multiple users in a
_single RE_ , we consider only one RE wherein $r$ users have transmitted their
packets. These users are decoded via only intra-RE iterations since there is
only a single RE under consideration. Let the set of users who have not yet
been decoded in the first $k-1$ intra-RE decoding iterations be denoted by
$\mathcal{S}_{k}$, and
$\mathcal{S}_{k}^{m}\triangleq\mathcal{S}_{k}\setminus\\{m\\}$, with
$\mathcal{S}_{1}=[r]$.333The set $\mathcal{S}_{k}$ as defined here is a slight
abuse of notation. In Sec. III, the set $\mathcal{S}_{k}$ consisted of users
being decoded via both intra-RE and inter-RE iterations, whereas here,
$\mathcal{S}_{k}$ consists of users being decoded via only intra-RE
iterations. In each intra-RE decoding iteration, a single user with the
highest SINR is decoded if their SINR $\geq\gamma_{\text{th}}$.

The SINR of the $m$th user in the $k$th intra-RE decoding iteration,
$\rho_{m}^{k}$, is calculated as seen before in Theorem 2. Specifically, when
users are only decoded via intra-RE SIC within one RE, we obtain the SINR as

$\displaystyle\rho_{m}^{k}=\frac{|{\mathbf{a}}_{m}^{kH}\hat{\mathbf{h}}_{m}^{k}|^{2}}{\|{\mathbf{a}}_{m}^{k}\|^{2}(N_{0}/P+\textstyle{\sum\nolimits_{i\in\mathcal{S}_{k}}}{\delta_{i}^{k}})+\textstyle{\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}}|{\mathbf{a}}_{m}^{kH}\hat{\mathbf{h}}_{i}^{k}|^{2}}.$
(24)

Here, $\delta_{i}^{k}$ is the error variance of the $i$th user in the $k$th
intra-RE decoding iteration, $\hat{\mathbf{h}}_{m}^{k}$ is the channel
estimate of the $m$th user, both obtained from Theorem 1, and
$\mathbf{a}_{m}^{k}$ is the combining vector for the $m$th user.444Since the
decoding process with intra-RE SIC involves only the RE in consideration, the
RE index and the APM are dropped in this section. Let $\rho_{\max}^{k}$ denote
the SINR of the user with the highest SINR in the $k$th intra-RE decoding
iteration, calculated as $\rho_{\max}^{k}=\max_{m\in\mathcal{S}_{k}}\
\rho_{m}^{k}$. Let $s$ be the index of the intra-RE decoding iteration in
which the reference packet is decoded, with $1\leq s\leq r$. Thus,
$\theta_{r}$ is calculated as

$\displaystyle\theta_{r}$
$\displaystyle=\text{Pr}(\rho_{\max}^{1}\geq\gamma_{\text{th}},\rho_{\max}^{2}\geq\gamma_{\text{th}},\ldots,\rho_{\max}^{s}\geq\gamma_{\text{th}}).$
(25)

Recall that the reference packet is tagged uniformly at random from the users.
With path loss inversion based power control, users have identical channel
statistics, and thus, $\theta_{r}$ is independent of which packet is tagged as
the reference packet.

The computation of the success probability $\theta_{r}$ is involved because
there is no clear relation between the peak SINRs across decoding iterations.
Also, the channel estimates of different users are correlated, across both the
user index and the decoding iteration index, making it difficult to use order
statistics. Further, $\theta_{r}$ is dependent on a large number of random
channel vectors, the order statistics of the peak SINRs, and the pilot
sequences of all the users. As a consequence, $\theta_{r}$ cannot be found in
closed form, and needs to be empirically evaluated. However, we present three
approximations to $\theta_{r}$, which are valid when perfect CSI is available
at the BS, i.e., there is no pilot contamination or estimation errors. These
lead to interpretable expressions for the SINR and $\theta_{r}$, and provide
upper bounds on the throughput with estimated CSI.

######  Theorem 3.

When perfect CSI is available at the BS, and MRC is used for decoding,
$\theta_{1}$ is given by

$\displaystyle\theta_{1}$
$\displaystyle=\Gamma_{\text{inc}}(N,\rho_{0}^{-1}\gamma_{\text{th}})/\Gamma(N),$
(26)

where $\rho_{0}\triangleq P\sigma_{h}^{2}/N_{0}$,
$\Gamma_{\text{inc}}(s,x)=\int_{x}^{\infty}t^{s-1}\,e^{-t}\,{\rm d}t$ is the
upper incomplete gamma function, and $\Gamma(s)$ is the ordinary gamma
function. For $r\geq 2$, the SINR with MRC and large $N$ can be computed as
${\rho}_{m}^{k}=N(\rho_{0}^{-1}+N\textstyle{\sum_{i\in\mathcal{S}_{k}^{m}}}t_{mi})^{-1},$
where
$t_{mi}\triangleq|\mathbf{h}_{m}^{H}\mathbf{h}_{i}|^{2}/(\|\mathbf{h}_{m}\|^{2}\|\mathbf{h}_{i}\|^{2})$.
With $t_{0}\triangleq\gamma_{\text{th}}^{-1}-N^{-1}\rho_{0}^{-1}$,
$\theta_{2}$ can be calculated as

$\displaystyle\theta_{2}$ $\displaystyle=\mathbbm{1}\\{t_{0}\geq
1\\}+(1-(1-t_{0})^{N})\mathbbm{1}\\{0\leq t_{0}\leq 1\\}.$ (27)

Three approximations to $\theta_{r}$ for $r\geq 3$ and large $N$ are described
below. Approximating $\rho_{\max}^{1}$ as $\rho_{1}^{1}$, and assuming $u_{m}$
as i.i.d. Gamma distributed with shape $r-1$ and rate $N$, we obtain the
_Gamma_ approximation:

$\displaystyle\emph{Gamma:}\ \ \ \theta_{r}$
$\displaystyle=1-\Gamma_{\text{inc}}(r-1,Nt_{0})/\Gamma(r-1).$ (28)

Approximating $\rho_{\max}^{1}=\rho_{1}^{1}$ and
$u_{m}\mathrel{\overset{\text{i.i.d.}}{\scalebox{1.5}[1.0]{$\sim$}}}\mathcal{N}((r-1)\mu_{T},(r-1)\sigma_{T}^{2})$,
where $\mu_{T}\triangleq(N+1)^{-1}$, and $\sigma_{T}^{2}\triangleq
N(N+1)^{-2}(N+2)^{-1}$, we obtain the _Normal_ approximation:

$\displaystyle\emph{Normal:}\ \ \ \theta_{r}$
$\displaystyle=1-\mathcal{Q}\left(\frac{t_{0}-(r-1)\mu_{T}}{\sqrt{r-1}\sigma_{T}}\right),$
(29)

where $\mathcal{Q}(\cdot)$ is the standard Normal Q-function. Finally, in the
_Deterministic_ approximation, the SINR becomes
${\rho}^{k}_{m}=N/(\rho_{0}^{-1}+r-k)$, and $\theta_{r}$ becomes

$\displaystyle\emph{Deterministic:}\ \ \
\theta_{r}=\mathbbm{1}\textstyle{\\{r\leq\lfloor
N/{\gamma_{\text{th}}}-\rho_{0}^{-1}+1\rfloor\\}}.$ (30)

###### Proof.

See supplementary material. ∎

_Remark:_ The above approximations provide closed form expressions for
$\theta_{r}$ and are valid when $N$ is large [23]. The first two
approximations have SINRs that are obtained by applying the theory of
deterministic equivalents to only the norms of the channels, and yields an
SINR that is affected only by the randomness in the multi-user interference
components. This is supported by the fact that the interference components
converge to their deterministic equivalents slower than the norms converge to
their deterministic equivalents [23]. The deterministic approximation follows
directly from Lemma 1, where the SINR is a deterministic quantity, and hence
$\theta_{r}$ is a binary function of $r$. With finite number of antennas, due
to small scale fading, the SINR of the users vary around this approximate
SINR. These variations affect the value of $\theta_{r}$, and are not captured
by the deterministic approximation, even though we obtain simple closed form
expressions for it. As a consequence, the throughput computed using the
deterministic approximation can be far from the actual throughput in certain
regimes and close to the actual throughput in other regimes, as will be seen
in Sec. VI-A.

###  V-D Final Evaluation of Throughput

After evaluating $\theta_{r}$ from Sec. V-C and substituting for $\xi_{c}$
from (19), we can simplify (23) to

$\displaystyle p_{i}$
$\displaystyle=1-e^{-\bar{c}q_{i}}\sum_{r=1}^{\infty}\theta_{r}\dfrac{({\bar{c}q_{i}})^{r-1}}{(r-1)!}\triangleq
f(q_{i}).$ (31)

Thus, $q_{i}=\lambda(p_{i-1})$ and $p_{i}=f(q_{i})$ are calculated alternately
as functions of each other as seen in (20) and (31). The procedure can be
initialized with either $q_{0}=1$ or $p_{0}=f(1)$.

The failure probability at the end of decoding is
$p_{\infty}=\lim_{i\rightarrow\infty}p_{i}$ and $(p_{\infty})^{d}$ is the
probability that a packet transmitted from a user with repetition factor $d$
does not get decoded at the receiver. Therefore, the asymptotic packet loss
rate (PLR), which is the fraction of packets that are not decoded at the BS,
is calculated as

PLR
$\displaystyle=\phi(p_{\infty})=\sum\nolimits_{d=2}^{{d_{\max}}}\phi_{d}(p_{\infty})^{d}.$
(32)

The asymptotic throughput of the system can now be obtained from the
asymptotic PLR as

$\displaystyle\mathcal{T}=L(1-{\textsf{PLR}}).$ (33)

The iterations $p_{i}=f(\lambda(p_{i-1}))$ converge asymptotically to
$p_{\infty}=0$ if the system load $L<L^{*}$ [6]. Here, $L^{*}$ is called the
_inflection load_ of the system: for any $L\geq L^{*}$, the system becomes
interference limited and the PLR does not converge to $0$ as $L$ increases.
Thus, for $L<L^{*}$, $p_{\infty}=0$ and therefore the asymptotic
$\textsf{PLR}=0$, and the throughput equals $L$. For $L\geq L^{*}$, the
throughput decreases monotonically with $L$.

##  VI Numerical Results

In this section, the previously derived SINR analysis is used to evaluate the
throughput of IRSA with estimated channels via Monte Carlo simulations, and
provide insights into the dependence of the system performance on the various
system parameters. In each simulation, independent realizations of the user
locations, the APM, and the fades experienced by the users are generated. The
throughput for each simulation is calculated as described in Sec. II-A, and
the effective system throughput $\mathcal{T}$ is calculated by averaging over
the simulations.

The results in this section are for $T=50$ REs, $N_{s}=10^{3}$ Monte Carlo
runs, $\lambda=10^{-2}$, $\alpha=3.76$, $\sigma_{h}^{2}=1$, SINR threshold
$\gamma_{\text{th}}=10$, MSBL threshold $\gamma_{\text{pr}}=10^{-6}$, cell
radius $r_{\max}=1000$ m, and reference distance $r_{0}=100$ m [20]. The
number of users contending for the $T$ REs is computed based on the load $L$
as $M=\lfloor LT\rceil$. The soliton distribution [10] with $d_{\max}=27$
maximum repetitions is used to generate the repetition factor $d_{m}$ for the
$m$th user, whose access pattern is formed by uniformly randomly choosing
$d_{m}$ REs from $T$ REs [6]. The APM is formed by stacking the pattern
vectors of all the users. The location of each user is uniformly sampled from
within a cell of radius $r_{\max}$ centered at the BS. The path loss
coefficient is calculated as $\beta_{m}=(r_{m}/r_{0})^{-\alpha}$ where $r_{m}$
is the radial distance of the $m$th user from the BS. The signal to noise
ratio (SNR) for the $m$th user is calculated as
$P\sigma_{h}^{2}\beta_{m}/N_{0}$. The received SNR of a user at the edge of
the cell at the BS is termed as the _cell edge SNR_ , and is denoted by
$\text{SNR}_{\text{edge}}$. The power levels of all users is chosen such that
the signal from a user at a distance $r_{\max}$ from the BS is received at
$\text{SNR}_{\text{edge}}$. This ensures that all users’ signals are received
at an SINR that at least $\text{SNR}_{\text{edge}}$ on average, in _singleton_
REs. If $\text{SNR}_{\text{edge}}\geq\gamma_{\text{th}}$, i.e., it is such
that the cell edge user’s signal is decodable, then all users’ signals are
decodable with high probability in singleton REs. The power levels of users is
set to $P=P_{\tau}=20$ dBm [20] and $N_{0}$ is chosen such that the cell edge
SNR is $10$ dB, unless otherwise stated. The pilot sequence for each user is
generated as
$\mathbf{p}_{m}\sim\mathcal{CN}(\mathbf{0}_{\tau},P_{\tau}\mathbf{I}_{\tau})$.

Figure 1: Comparison between MMSE, MSBL, and LCMMSE schemes. Figure 2: Impact
of pilot length $\tau$ on rate with MMSE.

Fig. 1 shows the effect of pilot length on the system throughput at different
$L$ under the three estimation schemes, with $N=16$. MMSE scheme performs the
best and reaches the optimal throughputs of $\mathcal{T}=L$ for very low pilot
lengths. MSBL scheme achieves the optimal throughputs for $L=1,2,3$ at
$\tau=4,8,12$ respectively, and beyond that, the performance is the same as
that of MMSE. This shows that with a few additional pilot symbols, we can do
away with the assumption of knowing the APM and path loss coefficients. LCMMSE
scheme matches MMSE for $L=1$ and for higher $L$, it needs a lot more pilot
symbols. This is because of both pilot contamination and low quality channel
estimates. Also, we note that the use of orthogonal pilots would require
$\tau\geq 50,100,150$ for $L=1,2,3$, respectively. The optimal throughput of
$\mathcal{T}=L$ is achieved with far fewer pilot symbols under all the three
estimation schemes. This is because a small subset of users transmit in any RE
in IRSA. Finally, under all the three schemes, we can achieve $\mathcal{T}\geq
1$, which is the maximum throughput achievable under perfectly coordinated
orthogonal access, i.e., grant-based orthogonal access. This shows the utility
of using IRSA as a GFRA protocol for MMTC, especially due to it’s high
performance at medium to high $L$. To summarize, the pilot length has a
significant impact on the performance of IRSA and yields near-optimal
throughputs at significantly lower pilot lengths than that required for
orthogonal pilot transmission. The drop in $\mathcal{T}$ at low pilot lengths
under estimated channels underscores the importance of accounting for the
effect of imperfect CSI in analyzing the performance of IRSA.

We focus on MMSE/LCMMSE hereafter in order to avoid clutter in the plots,
since MSBL matches the performance of MMSE with slightly higher $\tau$. In
Fig. 2, we investigate the effect of $L$, $\tau_{c}$ and $\gamma_{\text{th}}$
on the achievable rate $\mathcal{R}$ of the system with MMSE, with $N=16$.
Here, the rate is obtained as
$\mathcal{R}=(1-\tau/\tau_{c})\mathcal{T}\log_{2}(1+\gamma_{\text{th}})$
(bps/Hz), where $\tau_{c}$ is the total length of any user’s packet. Firstly,
we look at the effect of changing $\gamma_{\text{th}}$ by fixing
$\tau_{c}=100$. For $L=2$, $\gamma_{\text{th}}=20$ offers a higher rate than
$\gamma_{\text{th}}=10$, provided $\tau\geq 3$. Thus, at low loads, increasing
$\gamma_{\text{th}}$ (correspondingly, selecting a higher order modulation and
coding scheme) leads to better achievable rates. In contrast, when $L=4$,
$\gamma_{\text{th}}=10$ outperforms $\gamma_{\text{th}}=20$, because the
system is highly interference limited. Next, comparing $L=2,3,4$ for
$\tau_{c}=100$ and $\gamma_{\text{th}}=10$, we see that the rate improves with
$L$, provided the pilot length is large enough. Finally, decreasing $\tau_{c}$
reduces the achievable rate, as the relative overhead due to pilots increases.
Thus, at high loads, the throughput $\mathcal{T}$ limits the achievable rate,
while at low loads, the SINR threshold $\gamma_{\text{th}}$ is the primary
factor in determining the achievable rate.

Figure 3: Effect of number of antennas $N$ with MMSE. Figure 4: Effect of cell
edge SNR with MMSE.

In Fig. 3, we investigate the effect of the number of antennas at the BS, by
plotting the throughput with MMSE channel estimation for different $L$ and
$\tau$. Intuitively, we expect that, to achieve the optimal throughput of $L$,
we would require slightly more than $L\bar{d}$ antennas at the BS, since
$L\bar{d}$ users transmit packets per RE on average. The orthogonal pilots
curve is obtained by allocating $\tau=M=\lfloor LT\rceil$ for each $L$. Under
all configurations, it is observed that increasing $N$ has a significant
impact, and the peak throughput achieved reaches its maximum of
$\mathcal{T}=L$. Further, $\tau=10$ achieves a very similar performance as
that of orthogonal pilots, and $\tau=5$ performs poorly at low $N$ and high
$L$. For $L=2$, the throughput reaches the peak $\mathcal{T}=2$ for $N\geq 8$
for all three values of $\tau$. Similarly, for a high load of $L=3$, the
throughput reaches the peak, $\mathcal{T}=3$, for $N\geq 16$. For $L=2,3$,
since the average repetition factor $\bar{d}=3$, each RE is occupied by $6,9$
users, respectively. Thus, a slightly higher number of antennas is sufficient
to recover all the packets, provided accurate channel estimates are available
(i.e., $\tau$ is large enough). It is observed that at $L=2,\ N=4$ and $L=3,\
N=8$, improving $\tau$ greatly improves the throughput. Increasing the number
of antennas increases the array gain and the decoding capability of the
regularized zero forcing decoder at the BS, which in turn leads to more users
getting decoded. This shows the effectiveness of the number of antennas in
improving the throughput. Also, when $N=12$, the dramatic drop in the
throughput of $\mathcal{T}=3.8$ for $\tau=200$ (orthogonal pilots) to
$\mathcal{T}=1.2$ for $\tau=5$, which is around $70\%$ loss in performance,
shows that it is crucial to account for estimated CSI while analyzing the
performance of IRSA systems.

Fig. 4 shows the impact of cell edge SNR on the packet loss rate $\sf PLR$
with MMSE, with $N=16$. For SNR $<-5$ dB, the $\sf PLR$ is high, and in the
noise-limited regime ($-5<$ SNR $<0$ dB), an increase in cell edge SNR sharply
decreases the $\sf PLR$. For $L=4$, $\tau=5$, the system becomes interference-
limited, and thus the performance saturates at high SNR. This is because, at
low $\tau$, both signal and interference powers get scaled equally, and the
SINR remains roughly constant. Increasing $\tau$ from $5$ to $10$ and then to
orthogonal pilots, we observe that the $\sf PLR$ falls from $0.5$ to
$10^{-2.5}$ to $10^{-5}$. The higher $\tau$ and SNR result in accurate channel
estimates, and thus very low $\sf PLR$ is observed. Similarly, at $L=2$, the
drop of $\sf PLR$ from $10^{-1.7}$ to $10^{-2.8}$ to $10^{-3.9}$ for
$\tau=5,10$ and orthogonal pilots emphasizes the need to account for estimated
CSI when analyzing the performance of IRSA. In summary, the overall
performance can be improved by increasing the pilot length, number of
antennas, or cell edge SNR, but these need to be increased judiciously,
keeping the other parameters in mind.

Figure 5: Effect of regularization parameter and $\tau$ with MMSE.

Fig. 5 shows the effect of the regularization parameter, $\lambda$, on the
throughput of the system when MMSE estimation is employed, with $L=4$. As
$\lambda$ is varied from $10^{-6}$ to $1$, the curves go from ZF on the left
to RZF in the middle and finally to MRC on the right. For $N=4$, increasing
$\tau$ from $5$ to $10$ to $30$ only marginally improves the throughput. This
is because the system is highly interference limited, and hence channel
inversion does not work well at low $N$. For $\tau=5$, increasing $N$ from $4$
to $16$ to $32$ improves the performance due to the interference suppression
capability of RZF. Similar observations can be made for $\tau=10$ as well. MRC
does not have the interference suppression capability of RZF, and thus the
performance saturates at a low value for all $\tau$. We note that the optimal
throughput of $\mathcal{T}=4$ is obtained over a wide range of $\lambda$, and
thus precise optimization of $\lambda$ is not necessary to obtain near-optimal
throughputs.

Figure 6: Impact of load on $\mathsf{PLR}$ with LCMMSE.

Fig. 6 studies the impact of $L$ and $\tau$ on the system packet loss rate,
$\mathsf{PLR}$, evaluated with $N=16$, $\gamma_{\text{th}}=16$, and
$\lambda=1$. As the pilot length $\tau$ increases, better quality channel
estimates are obtained, and the corresponding SINR increases. In particular,
the system requires higher pilot lengths due to the use of LCMMSE estimates.
The loss rates reduce with increase in $\tau$, and gets closer to the
orthogonal loss rate. The $\sf PLR$ of perfectly coordinated orthogonal access
is the lowest. Similar to existing works, there is an error floor region where
the $\sf PLR$ is very low (upto $L=2$ for orthogonal pilots) after which the
$\sf PLR$ increases rapidly and is called the waterfall region. Here $L=2$
marks the inflection load, where the system transitions from the error floor
to the waterfall region.

Figure 7: Impact of power control on throughput.

In Fig. 7, the impact of power control on the throughput with LCMMSE is
characterized. For this plot, users transmit at powers that are dependent on
their distances from the BS. Specifically, the $m$th user, who is located at a
distance $r_{m}$ from the BS transmits at a power
$P(r_{m}/r_{0})^{\alpha-\zeta}$, making $\zeta$ the effective path loss
exponent. The cell edge SNR is fixed to 10 dB, and the throughputs are
obtained by varying $\zeta$ and $P$. When $\zeta=0$, the signals of the users
undergo pure fading, and the system achieves a peak throughput of
$\mathcal{T}=1.52$ at $L=1.6$. Further, as $L$ is increased, the throughput
drops to 0. The throughput of the system increases as $\zeta$ increases, until
$\zeta=2/3$. The exact $\zeta$ that yields the highest throughput is dependent
on other system parameters such as SNR, $\gamma_{\text{th}}$, and $N$. As
$\zeta$ is increased, the channel coefficients of the users become more
disparate, and thus offer a higher degree of capture effect. Beyond $\zeta=3$,
the throughput decreases as the exponent is so high that the received signal
power becomes comparable to the noise. For higher $\zeta$, the throughput
saturates as $L$ is increased since a few users are always decoded due to path
loss disparity. The channel fades and the path loss coefficients contribute to
the disparity amongst the channel coefficients of the users, and thus such a
system has higher throughputs than a system with only path loss [9] or only
fading [11]. Thus, it is useful to consider the combined effects of fading and
path loss in optimizing the performance.

###  VI-A Theoretical Validation of Throughput

Figure 8: Effect of $T$ on the throughput.

The results in this subsection are presented for $\tau=10$, cell edge SNR = 10
dB, $N=16$, $\lambda=10^{-2}$, $\gamma_{\text{th}}=16$, $d_{\max}=8$ maximum
repetitions and $N_{s}=10^{3}$ Monte Carlo runs. To reduce clutter in the
plots, we present the results for the lowest complexity (LCMMSE) channel
estimation scheme.

Fig. 8 investigates the effect of increasing the number of REs on the
throughput. The peak throughput increases from $\mathcal{T}=1.52$ at $L=1.6$
for $T=50$ to $\mathcal{T}=1.85$ at $L=1.85$ for $T=500$. Since $\bar{d}$ is
fixed, each user has a larger number of REs to choose from as $T$ is
increased. Thus, the interference reduces, and the throughput increases until
it reaches a peak and then drops off. The success probability $\theta_{r}$ is
evaluated empirically via $10^{4}$ Monte Carlo runs, and this in turn yields
the asymptotic theoretical throughput, which is marked as "DE". This can be
achieved as $M,T\rightarrow\infty$ with a fixed $L$. It is seen that this
asymptotic throughput increases linearly with the load until it hits a maximum
at the inflection load of the system, which occurs at $L^{*}=2$ in this case.
The throughput drops sharply beyond this load. The asymptotic throughput
provides an upper bound on the throughput achievable with finitely many REs
for low to moderate loads. At very low and high loads, the throughput achieved
with finitely many REs exactly matches with the DE asymptotic throughput. A
convenient operating point would be to set the system load to, say, 90% of the
inflection load, as, in this case, only finitely many REs would be sufficient
to achieve the asymptotic throughput. Finally, it can be observed that the
throughput of the system can be increased by increasing $T$, but only when the
system is operated at a load that is lower than the inflection load. Beyond
the inflection load, the system is always interference-limited and increasing
$T$ does not help.

Figure 9: Rate for different SINR thresholds.

In Fig. 9, the asymptotic rate of the system is plotted versus the system load
for different SINR thresholds with $\tau_{c}=100$. For a fixed
$\gamma_{\text{th}}$, $\mathcal{R}$ increases until the inflection load and
then drops off to zero. It is observed that a high $\mathcal{R}$ can be
achieved at lower loads by choosing a high $\gamma_{\text{th}}$, whereas, at
high loads, in order to serve more users, $\gamma_{\text{th}}$ must be kept
low. The choice of the threshold $\gamma_{\text{th}}$ decides the rate of
transmission, which in turn is related to the modulation and coding scheme to
be used. In summary, the SINR threshold $\gamma_{\text{th}}$, which depends on
the modulation and coding scheme employed and determines the data rate, can be
chosen based on the system parameters such as the number of antennas, training
duration, number of users/REs, and the transmit power.

Figure 10: Validation of theoretical approximations. Figure 11: Comparison of
approximations with simulation.

We now validate the approximations derived in Theorem 3 with the simulations
obtained with MRC, $d_{\max}=27$ maximum repetitions, and
$\gamma_{\textrm{th}}=10$. Fig. 10, reveals an inflection SNR∗ of $0$ dB and
$-7$ dB for $L=1,N=16$ and $L=2,N=64$ respectively, which behaves similar to
the inflection load $L^{*}$. Both the normal and the gamma approximations
match well with the asymptotic throughput obtained from the DE process. This
is because the deterministic approximation results in an SINR that is
completely deterministic and $\theta_{r}$ that is a binary function of $r$,
and consequently does not capture the statistics of the SINRs very well.
Further, the deterministic approximation results in a throughput that acts as
a step function since $\theta_{r}$ depends binarily on $N$,
$\gamma_{\text{th}}$, and SNR. As we go from $L=1,N=16$ to $L=2,N=64$, the
approximations become closer, and both the normal and the gamma approximations
match perfectly with the asymptotic throughput. In summary, the theoretical
curves with the approximations match the simulations when $N$ is increased, as
expected.

Fig. 11 examines the effect of $T$ on the approximations with $L=2$ and SNR
$=10$ dB. With finitely many REs, such as $T=50,100,300$, the throughput
achieves the optimal throughput $\mathcal{T}=2$ for $N=24,18,16$. The
asymptotic throughput obtained with DE provide an inflection $N^{*}=12$, which
matches perfectly with the normal approximation. The gamma approximation does
not match as well as the normal approximation. Here, the curves are with MRC
and perfect CSI, and the presented curves are valid upper bounds to the
throughputs with estimated CSI. These can be achieved with high enough $\tau$
as observed in Fig. 1, and thus the derived results provide very good
approximations to the asymptotic throughput achievable with estimated CSI.

##  VII Conclusions

This paper studied the effect of estimated CSI on the throughput of IRSA,
which is a distributed medium access protocol for mMTC involving repetition of
packets across different randomly selected REs. Decoding the users’ packets at
the BS involves successive interference cancellation. First, the channel
estimates were derived under three schemes: a sparsity-based scheme with MSBL,
MMSE, and LCMMSE. The corresponding SINR of all the users were obtained under
all three schemes accounting for pilot contamination, channel estimation
errors, path loss as well as multiple antennas at the BS. It was seen that
these errors significantly reduce the peak achievable throughput, even
resulting in up to 70% loss in certain regimes. Further, a density evolution
based analysis was presented to characterize the asymptotic performance of the
protocol when users perform path loss inversion based power control. Here,
several approximations to the success probability $\theta_{r}$ were derived
and it was seen that these approximations match well as the number of antennas
at the BS becomes large. Finally, several new insights into the design of
IRSA-based systems was discussed, namely, the improvement of the system
throughput, the evaluation of the operating load beyond which the system
becomes interference limited, and the choice of the decoding threshold
$\gamma_{\text{th}}$. The results underscored the importance of accounting for
practical channel estimation in studying the throughput offered by the IRSA
protocol. Future work could involve using differential evolution techniques
[24] to obtain the optimal repetition distribution that maximizes the
throughput in the finite frame length regime.

## References

  * [1] C. R. Srivatsa and C. R. Murthy, “Throughput analysis of PDMA/IRSA under practical channel estimation,” in _Proc. SPAWC_ , July 2019.
  * [2] H. Shariatmadari, R. Ratasuk, S. Iraji, A. Laya, T. Taleb, R. Jäntti, and A. Ghosh, “Machine-type communications: current status and future perspectives toward 5G systems,” _IEEE Commun. Mag._ , vol. 53, 2015.
  * [3] T. Xia, M. M. Wang, C. Jiang, J. Zhang, L. Wang, and X. You, “Practical machine-type communication for energy internet of things: An introduction,” _IEEE Commun. Mag._ , vol. 3, no. 1, pp. 48–59, 2019.
  * [4] L. Liu, E. G. Larsson, W. Yu, P. Popovski, C. Stefanovic, and E. de Carvalho, “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, 2018.
  * [5] 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, 2018.
  * [6] G. Liva, “Graph-based analysis and optimization of contention resolution diversity slotted ALOHA,” _IEEE Trans. Commun._ , vol. 59, no. 2, pp. 477–487, February 2011.
  * [7] E. Paolini, G. Liva, and M. Chiani, “Coded slotted ALOHA: A graph-based method for uncoordinated multiple access,” _IEEE Trans. Inf. Theory_ , vol. 61, no. 12, pp. 6815–6832, Dec 2015.
  * [8] S. Verdu _et al._ , _Multiuser detection_. Cambridge university press, 1998.
  * [9] E. E. Khaleghi, C. Adjih, A. Alloum, and P. Muhlethaler, “Near-far effect on coded slotted ALOHA,” in _Proc. PIMRC_ , Oct 2017.
  * [10] K. R. Narayanan and H. D. Pfister, “Iterative collision resolution for slotted ALOHA: An optimal uncoordinated transmission policy,” in _Proc. ISTC_ , Aug 2012, pp. 136–139.
  * [11] F. Clazzer, E. Paolini, I. Mambelli, and C. Stefanovic, “Irregular repetition slotted ALOHA over the Rayleigh block fading channel with capture,” in _Proc. ICC_ , May 2017.
  * [12] M. Ghanbarinejad and C. Schlegel, “Irregular repetition slotted ALOHA with multiuser detection,” in _Proc. WONS_ , 2013, pp. 201–205.
  * [13] J. Su, G. Ren, B. Zhao, and J. Ding, “Enhancing irregular repetition slotted ALOHA with polarization diversity in LEO satellite networks,” _KSII TIIS_ , vol. 14, no. 9, pp. 3907–3923, 2020.
  * [14] F. Clazzer, C. Kissling, and M. Marchese, “Enhancing contention resolution ALOHA using combining techniques,” _IEEE Trans. Commun._ , vol. 66, no. 6, pp. 2576–2587, 2017.
  * [15] F. Clazzer, F. Lázaro, G. Liva, and M. Marchese, “Detection and combining techniques for asynchronous random access with time diversity,” in _SCC 2017_ , 2017, pp. 1–6.
  * [16] C. R. Srivatsa and C. R. Murthy, “User activity detection for irregular repetition slotted ALOHA based MMTC,” _arXiv:2111.06140_ , 2021.
  * [17] S. Chen, B. Ren, Q. Gao, S. Kang, S. Sun, and K. Niu, “Pattern division multiple access-a novel nonorthogonal multiple access for fifth-generation radio networks,” _IEEE Trans. Veh. Technol._ , vol. 66, no. 4, pp. 3185–3196, April 2017.
  * [18] B. Hassibi and B. Hochwald, “How much training is needed in multiple-antenna wireless links?” _IEEE Trans. Inf. Theory_ , vol. 49, no. 4, pp. 951–963, 2003.
  * [19] D. P. Wipf and B. D. Rao, “An empirical Bayesian strategy for solving the simultaneous sparse approximation problem,” _IEEE Trans. Signal Process._ , vol. 55, no. 7, pp. 3704–3716, 2007.
  * [20] E. Björnson, J. Hoydis, and L. Sanguinetti, _Massive MIMO Networks: Spectral, Energy, and Hardware Efficiency_ , 2017.
  * [21] Wang Lin, Xiao Juan, and Guanrong Chen, “Density evolution method and threshold decision for irregular LDPC codes,” in _Proc. ICCCAS_ , vol. 1, 2004, pp. 25–28 Vol.1.
  * [22] A. Papoulis and H. Saunders, _Probability, random variables and stochastic processes_ , 1989.
  * [23] R. Couillet and M. Debbah, _Random matrix methods for wireless communications_. Cambridge University Press, 2011.
  * [24] R. Storn and K. Price, “Differential evolution – a simple and efficient heuristic for global optimization over continuous spaces,” _Journal of Global Optimization_ , 1997.

##  Appendix A Proof of Theorem 1

####  A-1 MMSE

We first vectorize the signal as

$\displaystyle\overline{\mathbf{y}}_{t}^{k}$
$\displaystyle\triangleq\text{vec}({\mathbf{Y}}_{t}^{pk})=(\mathbf{P}_{t}^{k*}\otimes\mathbf{I}_{N}){\mathbf{h}}_{t}^{k}+\overline{\mathbf{n}}_{t},$
(34)

where ${\mathbf{h}}_{t}^{k}\triangleq\text{vec}(\mathbf{H}_{t}^{k})$,
$\overline{\mathbf{n}}_{t}\triangleq\text{vec}(\mathbf{N}_{t}^{p})$, and
$\otimes$ is the Kronecker product. The MMSE estimator is
$\hat{{\mathbf{h}}}_{t}^{k}\triangleq\mathop{{}\mathbb{E}_{\mathbf{z}}}\left[{\mathbf{h}}_{t}^{k}\right]$,
where $\mathbf{z}=\overline{\mathbf{y}}_{t}^{k}$. The error
$\tilde{\mathbf{h}}_{t}^{k}$ $\triangleq$
$\hat{\mathbf{h}}_{t}^{k}-{\mathbf{h}}_{t}^{k}$ is uncorrelated with
$\mathbf{z}$ and the estimate. The conditional statistics of a Gaussian random
vector $\mathbf{x}$ are

$\displaystyle\mathop{{}\mathbb{E}_{\mathbf{z}}}\left[{\mathbf{x}}\right]$
$\displaystyle=\mathop{{}\mathbb{E}}\left[{\mathbf{x}}\right]+\mathbf{K}_{\mathbf{x}\mathbf{z}}\mathbf{K}_{\mathbf{z}\mathbf{z}}^{-1}\left(\mathbf{z}-\mathop{{}\mathbb{E}}\left[{\mathbf{z}}\right]\right),$
(35) $\displaystyle\mathbf{K}_{\mathbf{xx}|\mathbf{z}}$
$\displaystyle=\mathbf{K}_{\mathbf{x}\mathbf{x}}-\mathbf{K}_{\mathbf{x}\mathbf{z}}\mathbf{K}_{\mathbf{z}\mathbf{z}}^{-1}\mathbf{K}_{\mathbf{z}\mathbf{x}}.$
(36)

Here, $\mathbf{K}_{\mathbf{x}\mathbf{x}},$
$\mathbf{K}_{\mathbf{xx}|\mathbf{z}},$ and $\mathbf{K}_{\mathbf{x}\mathbf{z}}$
are the unconditional covariance of $\mathbf{x}$, the conditional covariance
of $\mathbf{x}$ conditioned on $\mathbf{z}$, and the cross-covariance of
$\mathbf{x}\ \&\ \mathbf{z}$ respectively. From (35), the MMSE channel
estimate $\hat{\mathbf{h}}_{t}^{k}$ can be calculated as

$\displaystyle\hat{\mathbf{h}}_{t}^{k}=\mathop{{}\mathbb{E}}{[\mathbf{h}_{t}^{k}]}+\mathop{{}\mathbb{E}}{[\mathbf{h}_{t}^{k}{\overline{\mathbf{y}}_{t}^{kH}}]}\mathop{{}\mathbb{E}}[\overline{\mathbf{y}}_{t}^{k}{\overline{\mathbf{y}}_{t}^{kH}}]^{-1}(\overline{\mathbf{y}}_{t}^{k}-\mathop{{}\mathbb{E}}{[\overline{\mathbf{y}}_{t}^{k}]}).$
(37)

The terms in the above expression can be evaluated as

$\displaystyle\mathop{{}\mathbb{E}}{[\mathbf{h}_{t}^{k}{\overline{\mathbf{y}}_{t}^{kH}}]}=\mathbf{B}_{t}^{k}\mathbf{P}_{t}^{kT}\otimes\mathbf{I}_{N},$
$\displaystyle\mathop{{}\mathbb{E}}[\overline{\mathbf{y}}_{t}^{k}{\overline{\mathbf{y}}_{t}^{kH}}]=(\mathbf{P}_{t}^{k*}\mathbf{B}_{t}^{k}\mathbf{P}_{t}^{kT}+N_{0}\mathbf{I}_{\tau})\otimes\mathbf{I}_{N},$
$\displaystyle\qquad\hat{\mathbf{h}}_{t}^{k}=(\mathbf{B}_{t}^{k}\mathbf{P}_{t}^{kT}(\mathbf{P}_{t}^{k*}\mathbf{B}_{t}^{k}\mathbf{P}_{t}^{kT}+N_{0}\mathbf{I}_{\tau})^{-1}\otimes\mathbf{I}_{N})\overline{\mathbf{y}}_{t}^{k},$

and thus, the MMSE estimate $\hat{\mathbf{H}}_{t}^{k}$ of $\mathbf{H}_{t}^{k}$
is

$\displaystyle\hat{\mathbf{H}}_{t}^{k}$
$\displaystyle={\mathbf{Y}}_{t}^{pk}(\mathbf{P}_{t}^{k}\mathbf{B}_{t}^{k}\mathbf{P}_{t}^{kH}+N_{0}\mathbf{I}_{\tau})^{-1}\mathbf{P}_{t}^{k}\mathbf{B}_{t}^{k},$
(38)
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}{\mathbf{Y}}_{t}^{pk}\mathbf{P}_{t}^{k}\mathbf{B}_{t}^{k}(\mathbf{P}_{t}^{kH}\mathbf{P}_{t}^{k}\mathbf{B}_{t}^{k}+N_{0}\mathbf{I}_{M_{t}^{k}})^{-1},$
(39)

where $(a)$ follows from $(\mathbf{AB}+\mathbf{I})^{-1}\mathbf{A}$ $=$
$\mathbf{A}(\mathbf{BA}+\mathbf{I})^{-1}$.

####  A-2 LCMMSE

The LCMMSE estimator is
$\hat{\mathbf{h}}_{tm}^{k}\triangleq\mathop{{}\mathbb{E}_{\mathbf{z}}}\left[{\mathbf{h}}_{tm}\right]$,
where $\mathbf{z}=\mathbf{y}_{tm}^{pk}$ is the received pilot signal. The
error $\tilde{\mathbf{h}}_{tm}^{k}$ $\triangleq$
$\hat{\mathbf{h}}_{tm}^{k}-{\mathbf{h}}_{tm}$ is uncorrelated with
$\mathbf{y}_{tm}^{pk}$ and the estimate. From (35), the LCMMSE channel
estimate $\hat{\mathbf{h}}_{tm}^{k}$ can be calculated

$\displaystyle\hat{\mathbf{h}}_{tm}^{k}=\mathop{{}\mathbb{E}}{[\mathbf{h}_{tm}{\mathbf{y}_{tm}^{pkH}}]}\mathop{{}\mathbb{E}}[\mathbf{y}_{tm}^{pk}{\mathbf{y}_{tm}^{pkH}}]^{-1}\mathbf{y}_{tm}^{pk}$
$\displaystyle=\dfrac{g_{tm}\beta_{m}\|\mathbf{p}_{m}\|^{2}\sigma_{h}^{2}}{N_{0}\|\mathbf{p}_{m}\|^{2}+\sum\nolimits_{i\in\mathcal{S}_{k}}|\mathbf{p}_{i}^{H}\mathbf{p}_{m}|^{2}g_{ti}\beta_{i}\sigma_{h}^{2}}\mathbf{y}_{tm}^{pk}\triangleq\eta_{tm}^{k}\mathbf{y}_{tm}^{pk}.$

####  A-3 MSBL

In each iteration of MSBL, two steps are performed. The first step, termed the
E-step, updates the covariance $\bm{\Sigma}_{kt}^{j+1}$ and mean
$\bm{\mu}_{ktn}^{j+1}$ of the posterior
$p([\mathbf{Z}_{t}^{k}]_{\mathrel{\mathop{\mathchar
58\relax}},n}|[\mathbf{Y}_{t}]_{\mathrel{\mathop{\mathchar
58\relax}},n},\bm{\gamma}_{kt}^{j})$

$\displaystyle\bm{\Sigma}_{kt}^{j+1}$
$\displaystyle\\!=\\!\bm{\Gamma}_{kt}^{j}-\bm{\Gamma}_{kt}^{j}\mathbf{P}^{kH}\\!(N_{0}\mathbf{I}_{\tau}\\!+\\!\mathbf{P}^{k}\bm{\Gamma}_{kt}^{j}\mathbf{P}^{kH})^{-1}\mathbf{P}^{k}\bm{\Gamma}_{kt}^{j},$
(40) $\displaystyle\bm{\mu}_{ktn}^{j+1}$
$\displaystyle=N_{0}^{-1}\bm{\Sigma}_{kt}^{j+1}\mathbf{P}^{kH}[\overline{\mathbf{Y}}_{t}^{pk}]_{\mathrel{\mathop{\mathchar
58\relax}},n},\ n\in[N].$ (41)

The second step, termed the M-step, updates the hyperparameter for the $i$th
user in the $t$th RE as

$\displaystyle\textstyle{[{\bm{\gamma}}_{kt}^{j+1}]_{i}=\dfrac{1}{N}\sum\limits_{n=1}^{N}([\bm{\Sigma}_{kt}^{j+1}]_{i,i}+|[\bm{\mu}_{ktn}^{j+1}]_{i}|^{2}),\
i\in[M^{k}].}$ (42)

This step estimates the variance of the channel of the $i$th user in the $t$th
RE. Based on the estimate $\hat{g}_{ti}^{k}$ and the true $g_{ti}$, the set of
users $[M^{k}]$ can be divided into four disjoint subsets

$\displaystyle\mathcal{A}_{t}^{k}$ $\displaystyle=\\{i\in[M^{k}]\ |\
\hat{g}_{ti}^{k}{g}_{ti}=1\\},$ (43a) $\displaystyle\mathcal{F}_{t}^{k}$
$\displaystyle=\\{i\in[M^{k}]\ |\ \hat{g}_{ti}^{k}(1-{g}_{ti})=1\\},$ (43b)
$\displaystyle\mathcal{M}_{t}^{k}$ $\displaystyle=\\{i\in[M^{k}]\ |\
(1-\hat{g}_{ti}^{k}){g}_{ti}=1\\},$ (43c) $\displaystyle\mathcal{I}_{t}^{k}$
$\displaystyle=\\{i\in[M^{k}]\ |\ (1-\hat{g}_{ti}^{k})(1-{g}_{ti})=1\\}.$
(43d)

$\mathcal{A}_{t}^{k}$ is the set of true positive users, $\mathcal{F}_{t}^{k}$
is the set of false positive users, $\mathcal{M}_{t}^{k}$ is the set of false
negative users, and $\mathcal{I}_{t}^{k}$ is the set of true negative users.
False positive and false negative users form the errors in APM estimation. As
the decoding iterations proceed, more users get decoded, and the errors in APM
estimation decrease. The MSBL channel estimate
$\hat{\mathbf{H}}_{t}^{k}={{\mathbf{Y}}^{pk}_{t}}\mathbf{P}^{k}{\hat{\bm{\Gamma}}_{kt}}(\mathbf{P}^{kH}\mathbf{P}^{k}{\hat{\bm{\Gamma}}_{kt}}+N_{0}\mathbf{I}_{M^{k}})^{-1}$
is output in the E-step from Algorithm 1, where ${\hat{\bm{\Gamma}}_{kt}}=$
diag$({\bm{\gamma}}_{kt}^{j_{\max}})$. The false negative users’ channels do
not get estimated even though they contribute towards
${{\mathbf{Y}}^{pk}_{t}}$. The false positive users’ channels get estimated
even though they haven’t transmitted, and thus, an erroneous channel estimate
is output for those users. Since $[\bm{\gamma}_{kt}]_{i}$ models the variance
of the $i$th users signal in the $t$th RE, it models
$g_{ti}\beta_{i}\sigma_{h}^{2}$. Thus, the estimated hyperparameter
$[{\bm{\gamma}}_{kt}^{j_{\max}}]_{i}$ would recover both $\hat{g}_{ti}^{k}$
and $\hat{\beta}_{i}^{k}$. Since the path loss is same across REs, a higher
quality estimate for the path loss can be estimated by averaging across REs,
and thus we obtain
$\hat{\beta}_{i}^{k}=(\sum\nolimits_{t=1}^{T}\hat{g}_{ti}^{k}[{\bm{\gamma}}_{kt}^{j_{\max}}]_{i})/(\sigma_{h}^{2}\sum\nolimits_{t=1}^{T}\hat{g}_{ti}^{k})$.

####  A-4 Error variances

The conditional covariance of ${\mathbf{h}}_{ti}$ is calculated conditioned on
$\mathbf{z}=\hat{\mathbf{h}}_{ti}^{k}$. In MMSE, with
$\mathbf{c}_{ti}^{k}=[\mathbf{C}_{t}^{k}]_{\mathrel{\mathop{\mathchar
58\relax}},i}$ and
$\mathbf{C}_{t}^{k}\triangleq\mathbf{P}_{t}^{k}\mathbf{B}_{t}^{k}(\mathbf{P}_{t}^{kH}\mathbf{P}_{t}^{k}\mathbf{B}_{t}^{k}+N_{0}\mathbf{I}_{M_{t}^{k}})^{-1}$,
we have

$\displaystyle\mathbf{K}_{{\mathbf{h}}_{ti}{\mathbf{h}}_{ti}}$
$\displaystyle=\mathbb{E}[{\mathbf{h}}_{ti}{\mathbf{h}}_{ti}^{H}]=\beta_{i}\sigma_{h}^{2}\mathbf{I}_{N},$
$\displaystyle\mathbf{K}_{{\mathbf{h}}_{ti}\mathbf{z}}$
$\displaystyle=\mathbb{E}[{\mathbf{h}}_{ti}\hat{\mathbf{h}}_{ti}^{kH}]=\mathbf{p}_{i}^{H}\mathbf{c}_{ti}^{k}g_{ti}\beta_{i}\sigma_{h}^{2}\mathbf{I}_{N},$
$\displaystyle\mathbf{K}_{{\mathbf{z}}\mathbf{z}}$
$\displaystyle=(N_{0}\|\mathbf{c}_{ti}\|^{2}+\textstyle{\sum\nolimits_{j\in\mathcal{S}_{k}}}|\mathbf{p}_{j}^{H}\mathbf{c}_{ti}^{k}|^{2}g_{tj}\beta_{j}\sigma_{h}^{2})\mathbf{I}_{N}.$

Thus, the conditional covariance is

$\displaystyle\mathbf{K}_{{\mathbf{h}}_{ti}{\mathbf{h}}_{ti}|\mathbf{z}}=\mathbf{K}_{{\mathbf{h}}_{ti}{\mathbf{h}}_{ti}}-\mathbf{K}_{{\mathbf{h}}_{ti}\mathbf{z}}\mathbf{K}_{\mathbf{z}\mathbf{z}}^{-1}\mathbf{K}_{\mathbf{z}{\mathbf{h}}_{ti}}$
$\displaystyle=\beta_{i}\sigma_{h}^{2}\left(\frac{N_{0}\|\mathbf{c}_{ti}^{k}\|^{2}+\sum\nolimits_{j\in\mathcal{S}_{k}^{i}}|r_{jti}^{k}|^{2}g_{tj}\beta_{j}\sigma_{h}^{2}}{N_{0}\|\mathbf{c}_{ti}^{k}\|^{2}+\sum\nolimits_{j\in\mathcal{S}_{k}}|r_{jti}^{k}|^{2}g_{tj}\beta_{j}\sigma_{h}^{2}}\right)\mathbf{I}_{N}\triangleq\delta_{ti}^{k}\mathbf{I}_{N},$

where $r_{jti}^{k}\triangleq\mathbf{p}_{j}^{H}\mathbf{c}_{ti}^{k}$ and
$\delta_{ti}^{k}$ accounts for pilot contamination. The conditional
autocorrelation follows as

$\displaystyle\mathbb{E}_{\mathbf{z}}[{\mathbf{h}}_{tm}{\mathbf{h}}_{tm}^{H}]$
$\displaystyle=\mathbf{K}_{{\mathbf{h}}_{tm}{\mathbf{h}}_{tm}|\mathbf{z}}+\mathbb{E}_{\mathbf{z}}[{\mathbf{h}}_{tm}]\mathbb{E}_{\mathbf{z}}[{\mathbf{h}}_{tm}]^{H}$
$\displaystyle=\delta_{tm}^{k}\mathbf{I}_{N}+\hat{\mathbf{h}}_{tm}^{k}\hat{\mathbf{h}}_{tm}^{kH}.$

The unconditional and conditional means of the estimation error are
$\mathbb{E}[\tilde{\mathbf{h}}_{tm}^{k}]=\mathbb{E}[\hat{\mathbf{h}}_{tm}^{k}-{\mathbf{h}}_{tm}]=0$
and
$\mathbb{E}_{\mathbf{z}}[\tilde{\mathbf{h}}_{tm}^{k}]=\mathbb{E}_{\mathbf{z}}[\hat{\mathbf{h}}_{tm}^{k}-{\mathbf{h}}_{tm}]=\hat{\mathbf{h}}_{tm}^{k}-\hat{\mathbf{h}}_{tm}^{k}=0.$
The conditional autocovariance of the error therefore simplifies as

$\displaystyle\mathbf{K}_{\tilde{\mathbf{h}}_{tm}^{k}\tilde{\mathbf{h}}_{tm}^{k}|\mathbf{z}}=\mathbb{E}_{\mathbf{z}}[\tilde{\mathbf{h}}_{tm}^{k}\tilde{\mathbf{h}}_{tm}^{kH}]$
$\displaystyle\ \ \
=\mathbb{E}_{\mathbf{z}}[{\mathbf{h}}_{tm}{\mathbf{h}}_{tm}^{H}]-\hat{\mathbf{h}}_{tm}^{k}\hat{\mathbf{h}}_{tm}^{kH}=\delta_{tm}^{k}\mathbf{I}_{N},$

and thus, $\delta_{tm}^{k}$ is also the variance of the estimation error.
Substituting $\mathbf{C}_{t}^{k}=\mathbf{P}_{t}^{k}\text{
diag}(\eta_{ti_{1}}^{k},\ldots,\eta_{ti_{M_{t}^{k}}}^{k}),$ we get the error
variance for LCMMSE as described in Table I.

The MSBL estimate error is also uncorrelated with the estimate and the error
variance can be derived similar to the MMSE scheme since the MSBL estimate is
a “plug-in" MMSE estimate. Since only true positive users’ channels are
estimated, the error variance is calculated only for the subset of true
positive users (users with $\hat{g}_{ti}^{k}g_{ti}=1$), and thus, each
$g_{ti}$ is accompanied by $\hat{g}_{ti}^{k}$ similar to [16]. Further, since
the error variance models the true interference from other true positive
users, the true path loss coefficient accompanies $\hat{g}_{ti}^{k}g_{ti}$.
Hence we define
$\mathbf{C}_{t}^{k}\triangleq\mathbf{P}^{k}\mathbf{D}_{t}^{k}(\mathbf{P}^{kH}\mathbf{P}^{k}\mathbf{D}_{t}^{k}+N_{0}\mathbf{I}_{M^{k}})^{-1}$
and $\mathbf{D}_{t}^{k}\triangleq$
diag$(d_{ti_{1}}^{k},d_{ti_{2}}^{k},\ldots,d_{ti_{M^{k}}}^{k})$, with
$d_{ti}^{k}=\hat{g}_{ti}^{k}g_{ti}\beta_{i}\sigma_{h}^{2}$. Substituting for
$\mathbf{C}_{t}^{k}$, we get the error variance for MSBL, listed in Table I.

##  Appendix B Proof of Theorem 2

In order to evaluate the SINR, we first calculate the power of the received
signal, which is calculated conditioned on the knowledge of the estimates
$\mathbf{z}\triangleq\text{vec}(\hat{\mathbf{H}}_{t}^{k})$ as
$\mathbb{E}_{\mathbf{z}}[|\tilde{y}_{tm}^{k}|^{2}]=\mathbb{E}_{\mathbf{z}}[|\sum_{i=1}^{4}T_{i}|^{2}]$.
Since noise is uncorrelated with data,
$\mathbb{E}_{\mathbf{z}}[T_{1}T_{4}^{H}]=\mathbb{E}_{\mathbf{z}}[T_{2}T_{4}^{H}]=\mathbb{E}_{\mathbf{z}}[T_{3}T_{4}^{H}]={0}$.
Since MMSE channel estimates are uncorrelated with their errors [20],
$\mathbb{E}_{\mathbf{z}}[T_{1}T_{2}^{H}]={0}$. Computing the remaining power
components requires the evaluation of $\mathbb{E}_{\mathbf{z}}[x_{i}x_{j}]$
for $i\neq j$ which can be calculated as
$\mathbb{E}_{\mathbf{z}}[x_{i}x_{j}]=\mathbb{E}_{\mathbf{z}}[x_{i}]\mathbb{E}_{\mathbf{z}}[x_{j}]=0$.
Thus, all the four terms are uncorrelated and the power in the received signal
is just a sum of the powers of the individual components
$\mathbb{E}_{\mathbf{z}}[|\tilde{y}_{tm}^{k}|^{2}]=\sum_{i=1}^{4}\mathbb{E}_{\mathbf{z}}[|T_{i}|^{2}]$.
We now compute the powers of each of the components. The useful signal power
is

$\displaystyle\mathbb{E}_{\mathbf{z}}[|T_{1}|^{2}]$
$\displaystyle=\mathbb{E}_{\mathbf{z}}[|{\mathbf{a}}_{tm}^{kH}\hat{\mathbf{h}}_{tm}^{k}g_{tm}x_{m}|^{2}]=Pg_{tm}^{2}|{\mathbf{a}}_{tm}^{kH}\hat{\mathbf{h}}_{tm}^{k}|^{2}.$

The desired gain is written as

$\displaystyle{\tt{Gain}}_{tm}^{k}$
$\displaystyle\triangleq\frac{\mathbb{E}_{\mathbf{z}}[|T_{1}|^{2}]}{P\|{\mathbf{a}}_{tm}^{k}\|^{2}}=g_{tm}\frac{|{\mathbf{a}}_{tm}^{kH}\hat{\mathbf{h}}_{tm}^{k}|^{2}}{\|{\mathbf{a}}_{tm}^{k}\|^{2}}.$
(44)

The power of the estimation error is expressed as

$\displaystyle\mathbb{E}_{\mathbf{z}}[|T_{2}|^{2}]$
$\displaystyle=\mathbb{E}_{\mathbf{z}}[|{\mathbf{a}}_{tm}^{kH}\tilde{\mathbf{h}}_{tm}^{k}g_{tm}x_{m}|^{2}]=Pg_{tm}^{2}\delta_{tm}^{k}\|{\mathbf{a}}_{tm}^{k}\|^{2}.$

Next, the power of the inter-user interference term $T_{3}$ is

$\displaystyle\mathbb{E}_{\mathbf{z}}[|T_{3}|^{2}]=\mathbb{E}_{\mathbf{z}}\left[\left|{\mathbf{a}}_{tm}^{kH}\textstyle{\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}}g_{ti}\mathbf{h}_{ti}x_{i}\right|^{2}\right]$
$\displaystyle\ \ \ \
=P\textstyle{\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}}g_{ti}^{2}{\mathbf{a}}_{tm}^{kH}\mathbb{E}_{\mathbf{z}}[\mathbf{h}_{ti}\mathbf{h}_{ti}^{H}]{\mathbf{a}}_{tm}^{k}$
$\displaystyle\ \ \ \
=P\textstyle{\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}}g_{ti}^{2}{\mathbf{a}}_{tm}^{kH}(\delta_{ti}^{k}\mathbf{I}_{N}+\hat{\mathbf{h}}_{ti}^{k}\hat{\mathbf{h}}_{ti}^{kH}){\mathbf{a}}_{tm}^{k}$
$\displaystyle\ \ \ \
=P\textstyle{\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}}g_{ti}^{2}(\|{\mathbf{a}}_{tm}^{k}\|^{2}\delta_{ti}^{k}+|{\mathbf{a}}_{tm}^{kH}\hat{\mathbf{h}}_{ti}^{k}|^{2}).$

Here,
$\mathbb{E}_{\mathbf{z}}[|T_{2}|^{2}]+\mathbb{E}_{\mathbf{z}}[|T_{3}|^{2}]$
represents the contribution of estimation errors and multi-user interference
components of the other users. Since $g_{ti}$ is binary, its powers are
dropped. We now split the normalized version of the above into the sum of the
error component ${\tt{Est}}_{tm}^{k}$ and the multi-user interference
${\tt{MUI}}_{tm}^{k}$ as follows

$\displaystyle{\tt{Est}}_{tm}^{k}$
$\displaystyle\triangleq\textstyle{\sum\nolimits_{i\in\mathcal{S}_{k}}}g_{ti}\delta_{ti}^{k},\
{\tt{MUI}}_{tm}^{k}\triangleq\textstyle{\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}}g_{ti}\frac{|{\mathbf{a}}_{tm}^{kH}\hat{\mathbf{h}}_{ti}^{k}|^{2}}{\|{\mathbf{a}}_{tm}^{k}\|^{2}}.$
(45)

The noise power is calculated as

$\displaystyle\mathbb{E}_{\mathbf{z}}[|T_{4}|^{2}]$
$\displaystyle=\mathbb{E}_{\mathbf{z}}[|{\mathbf{a}}_{tm}^{kH}\mathbf{n}_{t}|^{2}]=N_{0}\|{\mathbf{a}}_{tm}^{k}\|^{2}.$
(46)

A meaningful SINR expression can be written out by dividing the useful signal
power from (44) by the sum of the interference and the noise powers (from
(45), and (46)) [20]. Note that the interference component is comprised of the
estimation error term and the signal powers of other users who have also
transmitted in the same RE. SINR can thus be calculated as in (12) for all
users. For MMSE/LCMMSE, the corresponding SINR can be calculated by plugging
in the channel estimates as detailed in Theorem 2.

In MSBL, each of $T_{1},T_{2},$ and $T_{3}$ is calculated among the subset of
true positive users in the $t$th RE, i.e., users in
$\mathcal{A}_{t}^{k}=\\{i\in[M^{k}]|\hat{g}_{ti}^{k}g_{ti}=1\\}$. Hence, each
of the powers previously derived for MMSE is accompanied by
$\hat{g}_{ti}^{k}g_{ti}$. We need to account for false negative users, i.e.,
users in $\mathcal{M}_{t}^{k}=\\{i\in[M^{k}]|(1-\hat{g}_{ti}^{k})g_{ti}=1\\}$.
These users interfere with the decoding of other users and the SINR for such
users is $0$ since they will never get decoded. Such users’ signals are
uncorrelated with the other terms, and thus, their power is

$\displaystyle\mathbb{E}_{\mathbf{z}}[|T_{5}|^{2}]=\mathbb{E}_{\mathbf{z}}[|\textstyle{\sum\nolimits_{i\in\mathcal{S}_{k}^{m}\cap\mathcal{M}_{t}^{k}}}{\mathbf{a}}_{tm}^{kH}\mathbf{h}_{ti}g_{ti}x_{i}|^{2}]$
$\displaystyle\ \ \ \
\stackrel{{\scriptstyle(b)}}{{=}}P\textstyle{\sum\nolimits_{i\in\mathcal{S}_{k}^{m}\cap\mathcal{M}_{t}^{k}}}g_{ti}^{2}{\mathbf{a}}_{tm}^{kH}\mathbb{E}[\mathbf{h}_{ti}\mathbf{h}_{ti}^{H}]{\mathbf{a}}_{tm}^{k}$
$\displaystyle\ \ \ \
=P\textstyle{\sum\nolimits_{i\in\mathcal{S}_{k}^{m}\cap\mathcal{M}_{t}^{k}}}g_{ti}^{2}{\mathbf{a}}_{tm}^{kH}(\beta_{i}\sigma_{h}^{2}\mathbf{I}_{N}){\mathbf{a}}_{tm}^{k}$
$\displaystyle\ \ \ \
=P\textstyle{\sum\nolimits_{i\in\mathcal{S}_{k}^{m}\cap\mathcal{M}_{t}^{k}}}g_{ti}^{2}\beta_{i}\sigma_{h}^{2}\|{\mathbf{a}}_{tm}^{k}\|^{2},$
(47)

where the conditional expectation is dropped in $(b)$ since the BS does not
have the knowledge of the channel estimates of false negative users. The
normalised power of the false positive users is
${\tt{FNU}}_{tm}^{k}\triangleq\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}(1-\hat{g}_{ti}^{k})g_{ti}\beta_{i}\sigma_{h}^{2}.$

##  Appendix C Proof of Lemma 1

It is known that, as the number of antennas gets large, both
$\|\hat{\mathbf{h}}_{tm}^{k}\|^{2}$ and
$|\hat{\mathbf{h}}_{tm}^{kH}\hat{\mathbf{h}}_{ti}^{k}|^{2}$ converge almost
surely (a.s.) to their deterministic equivalents [23]. Evaluating the
deterministic equivalents as in [23] and plugging into the SINR expression
instead of the original terms, we can find an approximation to the SINR in the
high antenna regime. As $N$ gets large, the SINR with MRC converges almost
surely
($\rho_{tm}^{k}\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}\overline{\rho}_{tm}^{k}$)
to

$\displaystyle\overline{\rho}_{tm}^{k}$ $\displaystyle=\dfrac{N{\tt
Sig}_{tm}^{k}}{\epsilon_{tm}^{k}\left(N_{0}/P+{\tt NCoh}_{tm}^{k}\right)+{\tt
Coh}_{tm}^{k}},$ (48)

where ${\tt Sig}_{tm}^{k}$ is the desired gain, ${\tt NCoh}_{tm}^{k}$ is the
non-coherent interference, and ${\tt Coh}_{tm}^{k}$ is the coherent
interference. For LCMMSE, ${\tt NCoh}_{tm}^{k}\triangleq
g_{tm}\delta_{tm}^{k}+\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}g_{ti}\beta_{i}\sigma_{h}^{2}$,
${\tt Sig}_{tm}^{k}\triangleq
g_{tm}\beta_{m}^{2}\sigma_{h}^{4}\|{\mathbf{p}}_{m}\|^{4}$, ${\tt
Coh}_{tm}^{k}\triangleq
N\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}g_{ti}\beta_{i}^{2}\sigma_{h}^{4}|\mathbf{p}_{m}^{H}\mathbf{p}_{i}|^{2}$,
and $\epsilon_{tm}^{k}\triangleq
N_{0}\|{\mathbf{p}}_{m}\|^{2}+\sum\nolimits_{i\in\mathcal{S}_{k}}g_{ti}\beta_{i}\sigma_{h}^{2}|\mathbf{p}_{m}^{H}\mathbf{p}_{i}|^{2}$.
For MMSE, $\epsilon_{tm}^{k}\triangleq
N_{0}\|{\mathbf{c}}_{tm}^{k}\|^{2}+\sum\nolimits_{i\in\mathcal{S}_{k}}g_{ti}\beta_{i}\sigma_{h}^{2}|\mathbf{c}_{tm}^{kH}\mathbf{p}_{i}|^{2}$,
${\tt Sig}_{tm}^{k}\triangleq g_{tm}(\epsilon_{tm}^{k})^{2}$, ${\tt
Coh}_{tm}^{k}\triangleq
N\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}g_{ti}\beta_{i}^{2}\sigma_{h}^{4}|\mathbf{c}_{tm}^{kH}\mathbf{p}_{i}|^{2}$,
${\tt NCoh}_{tm}^{k}$ $\triangleq
g_{tm}\delta_{tm}^{k}+\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}g_{ti}\beta_{i}\sigma_{h}^{2}$.
For MSBL, $\epsilon_{tm}^{k}\triangleq
N_{0}\|{\mathbf{c}}_{tm}^{k}\|^{2}+\sum\nolimits_{i\in\mathcal{S}_{k}}g_{ti}\beta_{i}\sigma_{h}^{2}|\mathbf{c}_{tm}^{kH}\mathbf{p}_{i}|^{2}$,
${\tt
NCoh}_{tm}^{k}\triangleq\hat{g}_{tm}^{k}g_{tm}\delta_{tm}^{k}+\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}g_{ti}\beta_{i}\sigma_{h}^{2}$,
${\tt Sig}_{tm}^{k}\triangleq\hat{g}_{tm}^{k}g_{tm}(\epsilon_{tm}^{k})^{2}$,
and ${\tt Coh}_{tm}^{k}\triangleq
N\sum\nolimits_{i\in\mathcal{S}_{k}^{m}}g_{ti}\beta_{i}^{2}\sigma_{h}^{4}|\mathbf{c}_{tm}^{kH}\mathbf{p}_{i}|^{2}$.
Here, $\delta_{tm}^{k}$ and ${\mathbf{c}}_{tm}^{k}$ are obtained from Theorems
1 and 2, respectively, for the three estimation schemes. Here,
$\delta_{tm}^{k}$ and ${\mathbf{c}}_{tm}^{k}$ are obtained from Theorems 1 and
2, respectively, for the three estimation schemes. The above expressions are
obtained by replacing each of the terms involving $\hat{\mathbf{h}}_{tm}^{k}$
in (14) with their respective deterministic equivalents.

##  Appendix D Proof of Theorem 3

Let $k$ denote the intra-RE decoding iteration. When perfect CSI is available
at the BS and the users perform path loss inversion, the SINR of the $m$th
user in an RE is computed as

$\displaystyle{\rho}_{m}^{k}=\dfrac{P\|\mathbf{h}_{m}\|^{4}}{N_{0}\|\mathbf{h}_{m}\|^{2}+P\sum_{i\in\mathcal{S}_{k}^{m}}|\mathbf{h}_{m}^{H}\mathbf{h}_{i}|^{2}}.$
(49)

For $r=1$, ${\rho}_{1}^{1}=P\|\mathbf{h}_{m}\|^{2}/N_{0}$, and $\theta_{1}$
reduces to

$\displaystyle\theta_{1}$
$\displaystyle=\text{Pr}({\rho}_{1}^{1}\geq\gamma_{\text{th}})=\Gamma_{\text{inc}}(N,\rho_{0}^{-1}\gamma_{\text{th}})/\Gamma(N),$
(50)

where $\rho_{0}\triangleq P\sigma_{h}^{2}/N_{0}$,
$\Gamma_{\text{inc}}(s,x)=\int_{x}^{\infty}t^{s-1}\,e^{-t}\,{\rm d}t$ is the
upper incomplete gamma function and $\Gamma(s)$ is the ordinary gamma
function. The interference is written as
$t_{mi}=|\mathbf{h}_{m}^{H}\mathbf{h}_{i}|^{2}/(\|\mathbf{h}_{m}\|^{2}\|\mathbf{h}_{i}\|^{2})$,
where $t_{mi}\sim\text{Beta}(\alpha=1,\beta=N)$. We have
$\|\mathbf{h}_{i}\|^{2}/N\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}\sigma_{h}^{2}$
and
$\|\mathbf{h}_{i}\|^{4}/N^{2}\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}\sigma_{h}^{4}$
as $N\rightarrow\infty$ [23]. Thus, we can approximate the SINR as

$\displaystyle{\rho}_{m}^{k}$ $\displaystyle\approx
N(\rho_{0}^{-1}+N\textstyle{\sum_{i\in\mathcal{S}_{k}^{m}}}t_{mi})^{-1}.$ (51)

Here, we have applied the theory of deterministic equivalents to only the
channel norms and not to the interference. This is supported by the fact that
the interference components converge to their deterministic equivalents slower
than the norms converge to their deterministic equivalents [23].

For $r=2$, since $t_{12}=t_{21}$,
${\rho}_{1}^{1}={\rho}_{2}^{1}=N/(\rho_{0}^{-1}+Nt_{12})$. Thus,
${\rho}_{\max}^{1}=N/(\rho_{0}^{-1}+Nt_{12})$ and
${\rho}_{\max}^{2}=N\rho_{0}$ with ${\rho}_{\max}^{1}\leq{\rho}_{\max}^{2}$.
Thus, the success probability reduces to $\theta_{r}=$
Pr$({\rho}_{\max}^{1}\geq\gamma_{\rm th})$. Let
$t_{0}\triangleq\gamma_{\text{th}}^{-1}-N^{-1}\rho_{0}^{-1}$. Hence,
$\theta_{2}$ is calculated as

$\displaystyle\theta_{2}$
$\displaystyle\approx\text{Pr}({\rho}_{\max}^{1}\geq\gamma_{\text{th}})=\text{Pr}(t_{12}\leq
t_{0})$ $\displaystyle=\mathbbm{1}\\{t_{0}\geq
1\\}+(1-(1-t_{0})^{N})\mathbbm{1}\\{0\leq t_{0}\leq 1\\}.$ (52)

For $r\geq 3$, ${\rho}_{m}^{k}$ need not be a monotonically increasing
function of $k$ as seen in (51), and thus we cannot order the SINRs
meaningfully to compute a closed form expression for $\theta_{r}$. With
$u_{m}=\sum_{i\in[r]\setminus m}t_{mi}$, the maximum SINR in the first intra-
RE iteration is calculated as
$\rho_{\max}^{1}=\max_{m\in[r]}N(\rho_{0}^{-1}+Nu_{m})^{-1}$. Here, $u_{m}$ is
not independent across $m$ and it is not clear which $u_{m}$ is the minimum.
Thus, we approximate $\rho_{\max}^{1}$ as $\rho_{1}^{1}$, and upon dropping
the other SINR terms from (25), $\theta_{r}$ becomes

$\displaystyle\theta_{r}$
$\displaystyle\approx\text{Pr}(\rho_{1}^{1}\geq\gamma_{\text{th}})=\text{Pr}(u_{1}\leq
t_{0}).$ (53)

We now discuss two approximations to $u_{m}$ to evaluate $\theta_{r}$, with
the assumption that $u_{m}$ is independent across $m$.

Since $\lim_{N\rightarrow\infty}$ Beta$(\alpha=1,\beta=N)=\exp(\lambda=N),$ we
approximate $t_{mi}\sim\exp(N)$, which is a good approximation at high $N$
[22]. Even with this approximation, $u_{m}$ is identically Gamma distributed
across users but not independent. Thus, with the independence assumption,
$u_{m}$ is i.i.d. Gamma distributed with shape parameter $r-1$ and rate
parameter $N$, i.e.,
$u_{m}\mathrel{\overset{\text{i.i.d.}}{\scalebox{1.5}[1.0]{$\sim$}}}$
Gamma$(r-1,N)$. Thus, we obtain the _Gamma_ approximation:

$\displaystyle\theta_{r}$ $\displaystyle\approx
1-\Gamma_{\text{inc}}(r-1,Nt_{0})/\Gamma(r-1).$ (54)

Similarly, when we assume $t_{mi}$ is Normal distributed, $u_{m}$ is
identically Normal distributed across users but not independent. Let
$\mu_{T}=(N+1)^{-1}$ and $\sigma_{T}^{2}=N(N+1)^{-2}(N+2)^{-1}$. If we
approximate $t_{mi}\sim\mathcal{N}(\mu_{T},\sigma_{T}^{2})$ and $u_{m}$ is
independent across $m$, then
$u_{m}\mathrel{\overset{\text{i.i.d.}}{\scalebox{1.5}[1.0]{$\sim$}}}\mathcal{N}((r-1)\mu_{T},(r-1)\sigma_{T}^{2})$.
Thus, we obtain the _Normal_ approximation:

$\displaystyle\theta_{r}$ $\displaystyle\approx
1-\mathcal{Q}\left(\frac{t_{0}-(r-1)\mu_{T}}{\sqrt{r-1}\sigma_{T}}\right),$
(55)

where $\mathcal{Q}(\cdot)$ is the standard Normal Q-function.

A simpler expression can be obtained for $\theta_{r}$ by applying the theory
of deterministic equivalents to not just the channel norms but also to the
interference. Thus,
$|\mathbf{h}_{i}^{H}\mathbf{h}_{m}|^{2}/N\stackrel{{\scriptstyle\text{a.s.}}}{{\longrightarrow}}\sigma_{h}^{4}$,
as $N\rightarrow\infty$ [23]. Thus, the SINR becomes

$\displaystyle{\rho}^{k}_{m}$ $\displaystyle=N/(\rho_{0}^{-1}+r-k),$

which is not random and is a deterministic function of $N$ and $\rho_{0}$.
This expression for SINR follows from Lemma 1. Thus, we obtain the
_deterministic_ approximation:

$\displaystyle\theta_{r}$
$\displaystyle=\text{Pr}({\rho}^{1}_{1}\geq\gamma_{\text{th}})=\mathbbm{1}\textstyle{\\{r\leq\lfloor
N/{\gamma_{\text{th}}}-\rho_{0}^{-1}+1\rfloor\\}}.$ (56)