# Matrix Factorization Based Blind Bayesian Receiver for Grant-Free Random
Access in mmWave MIMO mMTC

Zhengdao Yuan, Fei Liu, Qinghua Guo, , Xiaojun Yuan, , Zhongyong Wang, and
Yonghui Li Z. Yuan is with the Artificial Intelligence Technology Engineering
Research Center, Open University of Henan, Zhengzhou 450002, China. He was
with the School of Electrical, Computer and Telecommunications Engineering,
University of Wollongong, Wollongong, NSW 2522, Australia (e-mail:
yuan_zhengdao@foxmail.com).F. Liu and Z. Wang are with the School of
Geoscience and Techonology and the School of Electrical and Information
Engineering, Zhengzhou University, Zhengzhou 450000, China, (e-mail:
<EMAIL_ADDRESS>zywangzzu@gmail.com)Q. Guo is with the School of
Electrical, Computer and Telecommunications Engineering, University of
Wollongong, Wollongong, NSW 2522, Australia (e-mail: qguo@uow.edu.au). X. Yuan
is with the National Key Laboratory on Wireless Communications, University of
Electronic Science and Technology of China, Chengdu 611731, China (e-mail:
xjyuan@uestc.edu.cn).Yonghui Li is with the School of Electrical and
Information Engineering, University of Sydney, Sydney, NSW 2006, Australia
(e-mail: yonghui.li@ sydney.edu.au)

###### Abstract

Grant-free random access is promising for massive connectivity with sporadic
transmissions in massive machine type communications (mMTC), where the hand-
shaking between the access point (AP) and users is skipped, leading to high
access efficiency. In grant-free random access, the AP needs to identify the
active users and perform channel estimation and signal detection.
Conventionally, pilot signals are required for the AP to achieve user activity
detection and channel estimation before active user signal detection, which
may still result in substantial overhead and latency. In this paper, to
further reduce the overhead and latency, we explore the problem of grant-free
random access without the use of pilot signals in a millimeter wave (mmWave)
multiple input and multiple output (MIMO) system, where the AP performs blind
joint user activity detection, channel estimation and signal detection
(UACESD). We show that the blind joint UACESD can be formulated as a
constrained composite matrix factorization problem, which can be solved by
exploiting the structures of the channel matrix and signal matrix. Leveraging
our recently developed unitary approximate message passing based matrix
factorization (UAMP-MF) algorithm, we design a message passing based Bayesian
algorithm to solve the blind joint UACESD problem. Extensive simulation
results demonstrate the effectiveness of the blind grant-free random access
scheme.

###### Index Terms:

Grant-free random access, user activity detection, massive machine type
communications, millimeter wave communications, signal detection, approximate
message passing, matrix factorization.

##  I Introduction

Massive machine type communications (mMTC) is one of the major application
scenarios in the fifth generation (5G) wireless communications, which can be
used to accommodate massive connections with sporadic transmissions in
internet of things (IoT) systems [1, 2]. The conventional grant-based access
technique needs a handshaking procedure to exchange the control signaling
between the access point (AP) and active users to establish the communication
links between them, which may result in excessive overhead, long and uncertain
latency. Due to the sporadic short-burst machine-type IoT data traffic, the
high overhead and long latency can be unacceptable as the communication
becomes inefficient due to the small amount of payload data [3]. Recently,
grant-free random access has attracted much attention, where the handshaking
procedure is skipped and users can transmit data straight away at any time
slot. In the 5G and future generation of wireless communications, the
millimeter-wave (mmWave) frequencies between 30-300 GHz will be employed to
alleviate the current spectrum shortage in sub-6GHz bands. In addition, the
multiple input and multiple output (MIMO) technology is used to further
improve the spectrum efficiency and communication reliability. Hence, in this
work, we consider the issue of grant-free random access in a mmWave MIMO
systems.

To achieve grant-free random access, the AP normally needs to perform user
activity detection to identify active users, channel estimation to acquire
channel state information of active users and multi-user signal detection to
detect non-orthogonal signals of the active users. User activity detection can
be coupled with channel estimation and/or active user signal detection, e.g.,
pilot assisted joint user activity detection and channel estimation is
followed by multi-user signal detection in[4, 5, 6, 7, 8, 9, 10]; channel
estimation is followed by performing user activity detection and signal
detection jointly in [11] and [12]; and joint user activity detection, channel
estimation and active user signal detection are performed with the aid of
pilot signals in [13] and [14]. Exploiting the fact that only a small fraction
of users in the network are active at a time, joint user activity detection
and active user signal detection or joint user activity detection and channel
estimation can be formulated as compressive sensing (CS) problems [15], [16],
which can be solved using sparse signal recovery algorithms. In [11] and [17],
with the assumption that the channel state information is available at the AP,
various methods such as those based on orthogonal matching pursuit (OMP) [18],
approximate message passing (AMP) [11] and prior-information aided adaptive
subspace pursuit [19] were developed for joint user activity detection and
active-user signal detection. However, in many scenarios, wireless channels
vary over time and have to be estimated frequently. In the existing works,
pilot signals are often employed for either user activity detection or joint
user activity detection and active user signal detection. The use of pilot
signals can result in a substantial overhead, especially in a MIMO grant-free
random access system where the number of channel coefficients to be estimated
is relatively large. In addition, the transmission of pilot signals (followed
by data) can also lead to considerable latency. To make the grant-free random
access more efficient, it is significant to further reduce the overhead and
latency, e.g., by reducing the number of pilot signals or even avoiding the
use of pilot signals [7, 20, 21].

In this work, we consider the problem of grant-free random access in a mmWave
MIMO system as illustrated in Fig. 1, where the AP is equipped with an antenna
array and provides service to a number of IoT users in the area. In
particular, we focus on uplink transmission. The users do not transmit pilot
signals, thereby the overhead and latency due to pilot signals can be avoided,
leading to significant improvement in access efficiency. However, this leads
to a challenging task for the AP, i.e., the AP needs to perform blind joint
user activity detection, channel estimation and signal detection (UACESD). In
this paper, we show that the blind joint UACESD can be formulated as a
constrained composite matrix factorization problem, which can be solved by
exploiting the strong constraints (structures) of the channel matrix and
signal matrix, i.e., the mmWave MIMO channel matrix exhibits sparsity in the
beam domain [22], [23] and the elements of the signal matrix are discrete-
valued, which are randomly drawn from the symbol alphabet. With these
constraints, we recover the channel matrix and signal matrix, and determine
the rank of the matrices simultaneously, thereby achieving joint blind UACESD.
It is worth mentioning that the blind joint UACESD problem is solved by
exploiting the constraints of the relevant matrices rather than relying on
pilot signals. Then we design a Bayesian method for the formulated constrained
composite matrix factorization problem. In particular, leveraging a recently
proposed unitary approximate message passing for matrix factorization (UAMP-
MF) algorithm [24], a message passing algorithm is designed in this work.
Extensive simulations are carried out to demonstrate the effectiveness of the
proposed grant-free random access scheme without using pilot signals.

Figure 1: Illustration of grant-free random access.

The rest of the paper is organized as follows. An brief introduction to
(unitary) approximate message passing and its application to matrix
factorization are presented in Section II. The system model for grant free
random access in mmWave MIMO and the problem formulation of joint blind UACESD
are described in Section III. The message passing based Bayesian blind joint
UACESD algorithm is developed in Section IV. Numerical simulation results are
provided in Section V, followed by conclusions in Section VI.

Notations \- The notations used in the paper are as follows. Boldface lower-
case and upper-case letters denote vectors and matrices, respectively, and
superscript $(\cdot)^{T}$ represents the transpose operation. A Gaussian
distribution of $x$ with mean $\hat{x}$ and variance $\nu_{x}$ is represented
by $\mathcal{N}(x;{\hat{x}},\nu_{x})$. We also simply use $\mathcal{N}(m,v)$
to represent a Gaussian distribution with mean $m$ and variance $v$. Notation
$\otimes$ represents the Kronecker product. The relation $f(x)=cg(x)$ for some
positive constant $c$ is written as $f(x)\propto g(x)$. We use
$\boldsymbol{a}\cdot\boldsymbol{b}$ and $\boldsymbol{a}\cdot/\boldsymbol{b}$
to represent the element-wise product and division between vectors
$\boldsymbol{a}$ and $\boldsymbol{b}$, respectively. The notation
$\boldsymbol{a}^{.-1}$ denotes the element-wise inverse operation to vector
$\boldsymbol{a}$. We use $|\boldsymbol{A}|]{{}^{.2}}$ to denote element-wise
magnitude squared operation for $\boldsymbol{A}$, and use
$||\boldsymbol{a}||^{2}$ to denote the squared $l_{2}$ norm of
$\boldsymbol{a}$. The notation $<\boldsymbol{a}>$ denotes the average
operation for $\boldsymbol{a}$, i.e., the sum of the elements of
$\boldsymbol{a}$ divided by the number of its elements.The probability density
function for the normal random matrix $\boldsymbol{X}$ is denoted as
$\boldsymbol{X}\sim\mathcal{MN}\big{(}\boldsymbol{X};\hat{\boldsymbol{X}},\boldsymbol{U}_{X},\boldsymbol{V}_{X}\big{)}$,
where $\hat{\boldsymbol{X}}$ denotes the mean of the random matrix
$\boldsymbol{X}$, $\boldsymbol{U}_{X}$ and $\boldsymbol{V}_{X}$ determine the
covariance among rows and columns of $\boldsymbol{X}$, respectively. We use 1
and 0 to denote an all-one vector and an all-zero vector with proper length,
respectively. Sometimes, we use a subscript $n$ for 1, i.e., $\textbf{1}_{n}$
to indicate its length $n$. The superscript of $\boldsymbol{a}^{t}$ denotes
the $t$-th iteration for $\boldsymbol{a}$. We use $[\boldsymbol{a}]_{n}$ to
denoted the $n$-the element of $\boldsymbol{a}$.

##  II Preliminary: (Unitary) Approximate Message Passing for Matrix
Factorization

The AMP algorithm was original developed based on the loopy belief propagation
with Gaussian and Taylor-series approximations [25] for compressive sensing
with the following model

$\mathbb{y}=\mathbb{A}\boldsymbol{x}+\mathbb{w},$ (1)

where $\mathbb{y}$ is an observation vector, $\mathbb{A}$ is a known
measurement matrix, $\mathbb{x}$ is a (sparse) vector to be recovered and
$\mathbb{w}$ is a noise vector with zero mean and covariance matrix
$\beta^{-1}\mathbf{I}$. AMP was later extended to solve estimation problems
with a generalized linear observation model in [26]. AMP enjoys low complexity
and it works well in the case of large independent and identically distributed
(sub)Gaussian $\mathbf{A}$, but it can easily diverge in the case of a generic
matrix $\mathbf{A}$ [27].

Inspired by the work in [28], the work in [29] shows that the robustness of
AMP can be improved remarkably through simple pre-processing, i.e., performing
a unitary transformation to the original linear model [29]. As any matrix
$\mathbf{A}$ has an SVD $\mathbf{A=U\Lambda V}$ with $\mathbf{U}$ and
$\mathbf{V}$ being two unitary matrices, performing a unitary transformation
with $\mathbf{U}^{H}$ leads to the following transformed model

$\boldsymbol{r}=\mathbf{\Phi}\boldsymbol{x}+\omega,$ (2)

where $\mathbf{r=U}^{H}\mathbf{y}$,
$\mathbf{\Phi}=\mathbf{U}^{H}\mathbf{A}=\mathbf{\Lambda}\mathbf{V}$,
$\mathbf{\Lambda}$ is a rectangular diagonal matrix, and
${\omega}=\mathbf{U}^{H}\mathbf{w}$ remains a zero-mean Gaussian noise vector
with the same covariance matrix $\beta^{-1}\mathbf{I}$. We can apply the
vector step size AMP [26] to model (2), leading to the first version of UAMP
(called UAMPv1) shown in Algorithm 1 (by replacing $\mathbf{r}$ and
$\mathbf{\Phi}$ with $\mathbf{y}$ and $\mathbf{A}$ in Algorithm 1
respectively, the original AMP algorithm is recovered).

Algorithm 1 UAMP (UAMPv2 executes operations in [ ])

Initialize $\boldsymbol{\tau}_{x}^{(0)}(\mathrm{or}~{}{\tau}_{x}^{(0)})>0$ and
${{\boldsymbol{x}}^{(0)}}$. Set $\boldsymbol{s}^{(-1)}=\mathbf{0}$ and $t=0$.
Define vector ${\boldsymbol{\lambda}=\mathbf{\Lambda\Lambda}^{H}\textbf{1}}$.
Repeat

1:$\boldsymbol{\tau}_{p}$ = $\mathbf{|\Phi|}^{.2}\boldsymbol{\tau}^{t}_{x}$
$\left[\mathrm{or}~{}\boldsymbol{\tau}_{p}=\tau^{t}_{x}\boldsymbol{\lambda}\right]$

2:$\boldsymbol{p}=\mathbf{\Phi}{{\boldsymbol{x}}^{t}}-\boldsymbol{\tau}_{p}\cdot\boldsymbol{s}^{t-1}$

3:$\boldsymbol{\tau}_{s}=\mathbf{1}./(\boldsymbol{\tau}_{p}+\beta^{-1}\mathbf{1})$

4:$\boldsymbol{s}^{t}=\boldsymbol{\tau}_{s}\cdot(\boldsymbol{r}-\boldsymbol{p})$

5:$\mathbf{1}./\boldsymbol{\tau}_{q}$ =
$|\mathbf{\Phi}^{H}|^{.2}\boldsymbol{\tau}_{s}$
$\left[\mathrm{or}~{}\mathbf{1}./\boldsymbol{\tau}_{q}=(\frac{1}{N}\boldsymbol{\lambda}^{H}\boldsymbol{\tau}_{s})\mathbf{1}\right]$

6:$\boldsymbol{q}={{\boldsymbol{x}}^{t}}+\boldsymbol{\tau}_{q}\cdot(\mathbf{\Phi}^{H}\boldsymbol{s}^{t})$

7:$\boldsymbol{\tau}_{x}^{t+1}$ = $\boldsymbol{\tau}_{q}\cdot
g_{x}^{\prime}(\boldsymbol{q},\boldsymbol{\tau}_{q})$
$\left[\mathrm{or}~{}\tau_{x}^{t+1}\\!=\\!\frac{1}{N}\mathbf{1}^{H}\left(\boldsymbol{\tau}_{q}\cdot
g_{x}^{\prime}(\boldsymbol{q},\tau_{q})\right)\right]$

8:${{\mathbf{x}}^{t+1}}=g_{x}(\boldsymbol{q},\boldsymbol{\tau}_{q})$

9:$t=t+1$

Until terminated

Applying an average operation to two vectors: $\boldsymbol{\tau}_{x}$ in Line
7 and $|\mathbf{\Phi}^{H}|^{.2}\boldsymbol{\tau}_{s}$ in Line 5 of UAMPv1 in
Algorithm 1 leads to the second version of UAMP with lower complexity, where
the number of matrix-vector products is reduced from 4 in UAMPv1 (or the
original AMP) to 2 per iteration (refer to [30] for detailed derivation).

In the (U)AMP algorithms, the function
$g_{x}(\mathbf{q},\boldsymbol{\tau}_{q})$ is related to the prior of
$\mathbf{x}$, which returns a column vector with the $n$th element
$[g_{x}(\mathbf{q},\boldsymbol{\tau}_{q})]_{n}$ given as

$[g_{x}(\mathbf{q},\boldsymbol{\tau}_{q})]_{n}=\frac{\int
x_{n}p(x_{n})\mathcal{N}(x_{n};q_{n},\tau_{q_{n}})dx_{n}}{\int
p(x_{n})\mathcal{N}(x_{n};q_{n},\tau_{q_{n}})dx_{n}},$ (3)

where $p(x_{n})$ represents a prior for $x_{n}$. The function
$g_{x}^{\prime}(\mathbf{q},\boldsymbol{\tau}_{q})$ returns a column vector and
the $n$th element is denoted by
$[g_{x}^{\prime}(\mathbf{q},\boldsymbol{\tau}_{q})]_{n}$, where the derivative
is taken with respect to $q_{n}$.

In [24], leveraging the variational inference [31] and UAMP, an algorithm
called UAMP-MF is developed to solve the matrix factorization problem with
model

$\boldsymbol{R}=\boldsymbol{A}\boldsymbol{X}+\boldsymbol{W}$ (4)

where $\boldsymbol{R}$ is an observation matrix, $\boldsymbol{W}$ is a noise
matrix, and matrices $\boldsymbol{A}$ and $\boldsymbol{X}$ are two factor
matrices to be recovered. To achieve high performance, instead of using the
mean field approximation with full factorization, the two matrices are
decoupled and are treated as two high dimensional latent variables, leading to
the updates of two distributions on the two matrices. By exploiting the
structure of involved variational messages and through a covariance matrix
whitening process, UAMPv1 is incorporated to variational inference to
efficiently deal with the distribution updates of the two matrices. The method
can be implemented using message passing with UAMP as its key component. UAMP-
MF inherits the low complexity and robustness of UAMP, and can be used to deal
with many matrix factorization problems [24]. The UAMP-MF algorithm is shown
in Algorithm 2, where the quantity $C$ in Line 25 is computed as

$\displaystyle
C=\text{Tr}\Big{(}\big{(}\boldsymbol{Y}-\hat{\boldsymbol{H}}\hat{\boldsymbol{X}}\big{)}^{\textrm{T}}\big{(}\boldsymbol{Y}-\hat{\boldsymbol{H}}\hat{\boldsymbol{X}}\big{)}\Big{)}+\text{Tr}\Big{(}\hat{\boldsymbol{X}}\hat{\boldsymbol{X}}^{\textrm{T}}\text{Tr}(\boldsymbol{U}_{H})\boldsymbol{V}_{H}$
$\displaystyle+\text{Tr}(\boldsymbol{V}_{X})\boldsymbol{U}_{X}\hat{\boldsymbol{H}}^{\textrm{T}}\hat{\boldsymbol{H}}+\text{Tr}(\boldsymbol{V}_{X})\boldsymbol{U}_{X}\text{Tr}(\boldsymbol{U}_{H})\boldsymbol{V}_{H}\Big{)}.$
(5)

In this work, we will formulate the blind joint UACESD as a constrained
composite matrix factorization problem, and leveraging the UAMP-MF algorithm,
we will develop an efficient message passing algorithm to solve the
constrained composite matrix factorization problem.

Algorithm 2 UAMP-MF

Initialization: $\boldsymbol{U}_{H}=\boldsymbol{I}_{M}$,
$\boldsymbol{V}_{H}=\boldsymbol{I}_{N}$,
$\hat{\boldsymbol{H}}=\textbf{1}_{MN}$,
$\boldsymbol{V}_{X}=\boldsymbol{I}_{L}$.
$\boldsymbol{\Xi}_{X}=\textbf{1}_{NL}$, $\boldsymbol{S}_{X}=\textbf{0}_{NL}$,
$\boldsymbol{\Xi}_{H}=\textbf{1}_{MN}$, $\boldsymbol{S}_{H}=\textbf{0}_{MN}$.
Repeat

1:$\overline{\boldsymbol{W}}_{X}=\hat{\boldsymbol{H}}^{\textrm{H}}\hat{\boldsymbol{H}}+M\boldsymbol{V}_{H}$

2:$[\boldsymbol{C}_{X},\boldsymbol{D}_{X}]=\text{eig}(\overline{\boldsymbol{W}}_{X})$

3:$\boldsymbol{R}_{X}=\boldsymbol{D}_{X}^{-\frac{1}{2}}\boldsymbol{C}_{X}^{\textrm{H}}\hat{\boldsymbol{H}}^{\textrm{H}}\boldsymbol{Y}$,
$\boldsymbol{\Phi}_{X}=\boldsymbol{D}_{X}^{-\frac{1}{2}}\boldsymbol{C}_{X}^{\textrm{H}}$

4:$\boldsymbol{V}_{P_{X}}=|\boldsymbol{\Phi}_{X}|^{.2}\boldsymbol{\Xi}_{X}$

5:$\boldsymbol{P}_{X}=\boldsymbol{\Phi}_{X}\hat{\boldsymbol{X}}-\boldsymbol{V}_{\boldsymbol{P}_{X}}\cdot\boldsymbol{S}_{X}$

6:$\boldsymbol{V}_{S_{X}}=\textbf{1}./(\boldsymbol{V}_{P_{X}}+\hat{\lambda}^{-1}\mathbf{1})$

7:$\boldsymbol{S}_{X}=\boldsymbol{V}_{S_{X}}\cdot(\boldsymbol{R}_{X}-\boldsymbol{P}_{X})$

8:$\boldsymbol{V}_{Q_{X}}=\textbf{1}./(|\boldsymbol{\Phi}_{X}^{\textrm{H}}|^{.2}\boldsymbol{V}_{S_{X}})$

9:$\boldsymbol{Q}_{X}=\hat{\boldsymbol{X}}+\boldsymbol{V}_{Q_{X}}\cdot(\boldsymbol{\Phi}_{X}^{\textrm{H}}\boldsymbol{S}_{X})$

10:$\boldsymbol{\Xi}_{X}=\boldsymbol{V}_{Q_{X}}\cdot\boldsymbol{G}_{X}^{\prime}(\boldsymbol{Q}_{X},\boldsymbol{V}_{Q_{X}})$

11:$\hat{\boldsymbol{X}}=\boldsymbol{G}_{X}(\boldsymbol{Q}_{X},\boldsymbol{V}_{Q_{X}})$

12:$\boldsymbol{U}_{X}=\text{diag}(\text{mean}(\boldsymbol{\Xi}_{X},2))$

13:$\overline{\boldsymbol{W}}_{H}=\hat{\boldsymbol{X}}\hat{\boldsymbol{X}}^{\textrm{H}}+L\boldsymbol{U}_{X}$

14:$[\boldsymbol{C}_{H},\boldsymbol{D}_{H}]=\text{eig}(\overline{\boldsymbol{W}}_{H})$

15:$\boldsymbol{R}_{H}=\boldsymbol{D}_{H}^{-\frac{1}{2}}\boldsymbol{C}_{H}^{\textrm{H}}\hat{\boldsymbol{X}}\boldsymbol{Y}^{\textrm{H}}$,
$\boldsymbol{\Phi}_{H}=\boldsymbol{D}_{H}^{-\frac{1}{2}}\boldsymbol{C}_{H}^{\textrm{H}}$

16:$\boldsymbol{V}_{P_{H}}=|\boldsymbol{\Phi}_{H}|^{.2}\boldsymbol{\Xi}_{H}^{\textrm{H}}$

17:$\boldsymbol{P}_{H}=\boldsymbol{\Phi}_{H}\hat{\boldsymbol{H}}^{\textrm{H}}-\boldsymbol{V}_{\boldsymbol{P}_{H}}\cdot\boldsymbol{S}_{H}$

18:$\boldsymbol{V}_{S_{H}}=\mathbf{1}./(\boldsymbol{V}_{P_{H}}+\hat{\lambda}^{-1}\mathbf{1})$

19:$\boldsymbol{S}_{H}=\boldsymbol{V}_{S_{H}}\cdot(\boldsymbol{R}_{H}-\boldsymbol{P}_{H})$

20:$\boldsymbol{V}_{Q_{H}}=\mathbf{1}./(|\boldsymbol{\Phi}_{H}^{\textrm{H}}|^{.2}\boldsymbol{V}_{S_{H}})$

21:$\boldsymbol{Q}_{H}=\hat{\boldsymbol{H}}^{\textrm{H}}+\boldsymbol{V}_{Q_{H}}\cdot(\boldsymbol{\Phi}_{H}^{\textrm{H}}\boldsymbol{S}_{H})$

22:$\boldsymbol{\Xi}_{H}=\boldsymbol{V}_{Q_{H}}\cdot\boldsymbol{G}_{H}^{\prime}(\boldsymbol{Q}_{H}^{\textrm{H}},\boldsymbol{V}_{Q_{H}}^{\textrm{T}})$

23:$\hat{\boldsymbol{H}}=\boldsymbol{G}_{H}(\boldsymbol{Q}_{H}^{\textrm{H}},\boldsymbol{V}_{Q_{H}}^{\textrm{T}})$

24:$\boldsymbol{U}_{H}=\text{diag}(\text{mean}(\boldsymbol{\Xi}_{H},1))$

25:$\hat{\lambda}=ML/C$ with $C$ given in (5)

Until terminated

##  III System Model and Problem Formulation

###  III-A System Model

Consider a mmWave MIMO system, where an AP equipped with $M$ antennas is used
to serves $U$ user devices, each of which has a single antenna due to the
consideration of the size and the cost of user devices. However, we note that
the extension of this work to the case of user devices with multiple antennas
is straightforward. The number of active devices in a time slot is denoted by
$N$, which can be much smaller than $U$. We consider grant-free uplink
transmission. We assume that the MIMO channel is statistic over $L$
consecutive symbol intervals, which is called a frame. The information bit
sequence $\boldsymbol{b}_{n}$ of active user $n$ is mapped to a symbol
sequence $\boldsymbol{x}_{n}^{T}\in\mathbb{C}^{1\times L}$. Then the signal
received by the AP over a frame can be represented as

$\boldsymbol{Y}=\boldsymbol{H}\boldsymbol{X}+\boldsymbol{W},$ (6)

where $\boldsymbol{Y}=[Y_{ml}]\in\mathbb{C}^{M\times L}$ represents the
received signal matrix,
$\boldsymbol{X}=[\boldsymbol{x}_{1},...,\boldsymbol{x}_{N}]^{T}\in\mathbb{C}^{N\times
L}$ denotes the transmitted signal matrix,
$\boldsymbol{H}=[H_{mn}]\in\mathbb{C}^{M\times N}$ is the MIMO channel matrix
and $\boldsymbol{W}=[W_{ml}]\in\mathbb{C}^{M\times L}$ denotes the additive
temporally and spatially white Gaussian noise with zero mean and variance
$\sigma_{\omega}^{2}$. We assume that no pilot signals are transmitted by the
active users, and their IDs are included in the symbol sequences.

Figure 2: Illustration of the system in base-band.

The beamspace channel model is widely used due to the highly directional
nature of propagation at mmWave frequencies. For simplicity, we assume that
the antenna array employed by the AP is a uniform linear one with $M$
antennas. Then the steering vector of the antenna array can be represented as

$\boldsymbol{a}(\theta)=\left[1,e^{-j2\pi\vartheta},...,e^{-j2\pi\vartheta(M-1)}\right]^{T},$
(7)

where the normalized spatial angle
$\vartheta=\frac{d}{\lambda^{\prime}}\sin(\theta)$ with $\theta\in[0,\pi]$
being the physical angle, $d$ denotes the antenna spacing and
$\lambda^{\prime}$ is the wavelength. The beamspace representation using
uniformly spaced spatial angles $\vartheta_{k}=k/K$ with $k=0:K-1$ can be
represented as

$\boldsymbol{H}=\boldsymbol{F}\boldsymbol{G}$ (8)

where $\boldsymbol{F}\in\mathbb{C}^{M\times K}(K>M)$ is a partial DFT matrix
[32], i.e., the first $M$ rows of the $K$-point DFT matrix. Due to the limited
number of scatterers in a typical mmWave environment,
$\boldsymbol{G}\in\mathbb{C}^{K\times N}$ is a sparse matrix. Substituting (8)
into (6) we have

$\displaystyle\boldsymbol{Y}$ $\displaystyle=$
$\displaystyle\boldsymbol{H}\boldsymbol{X}+\boldsymbol{W}$ (9)
$\displaystyle=$
$\displaystyle\boldsymbol{F}\boldsymbol{G}\boldsymbol{X}+\boldsymbol{W}.$

In this paper, to further reduce the latency of the grant-free random access
system and avoid the overhead, we assume that the users do not transmit pilot
symbols. The task of the receiver at the AP is to directly detect the active
users and decode their signals blindly. Specifically,

  * •

The matrix $\boldsymbol{F}$ is known, but neither $\boldsymbol{G}$ nor
$\boldsymbol{X}$ in (9) is known, which need to be estimated jointly based on
$\boldsymbol{Y}$. This is a composite matrix factorization problem. The
factorization is possible if we can fully exploit the structures of the
matrices, i.e., $\boldsymbol{G}$ is a sparse matrix and the entries of
$\boldsymbol{X}$ are discrete-valued symbols, which are strong constraints for
the composite factorization problem.

  * •

As the number of active users is unknown, the number of columns of
$\boldsymbol{H}$ (or the number of rows of $\boldsymbol{X}$) (9) is unknown,
which needs to be determined.

In this work, we will design a Bayesian method to recover $\boldsymbol{G}$ and
$\boldsymbol{X}$. Leveraging the UAMP-MF algorithm, an efficient message
passing algorithm will be developed.

###  III-B Probabilistic Formulation

We will develop a Bayesian method to recover $\boldsymbol{G}$ and
$\boldsymbol{X}$ simultaneously, where we specify proper priors for them to
reflect the constraints on them. As the entries of matrix $\boldsymbol{G}$ are
sparse and independent, we use the sparsity-inducing Bernoulli-Gaussian
distribution [33] as the prior of $\boldsymbol{G}$, i.e.,

$p(\boldsymbol{G})=\prod_{n,k}p(g_{kn})=\prod_{k,n}\left((1-\epsilon)\delta(g_{kn})+\epsilon\mathcal{N}(g_{kn};0,\nu)\right),$
(10)

where $g_{kn}$ is the $(k,n)$th element of matrix $\boldsymbol{G}$,
$\delta(\cdot)$ denotes the Dirac delta function, and the parameter $\epsilon$
represents the sparsity rate. Here we assume that the number of active users
$N$ is known, and its estimation will be discussed later in Section IV.D. The
$n$-th row of $\boldsymbol{X}$ stand for the transmitted symbols of the $n$-th
active users, which are discrete valued and randomly drawn from the alphabet
$\mathcal{A}$. So the prior distribution of $\boldsymbol{X}$ can be expressed
as

$\displaystyle
p(\boldsymbol{X})=\prod_{n,l}p(x_{n,l})=\frac{1}{|\mathcal{A}|}\prod_{n,l}\sum_{a=1}^{|\mathcal{A}|}\delta(x_{n,l}-\alpha_{a}),$
(11)

where $|\mathcal{A}|$ denotes the size of the alphabet, and
$\mathcal{A}=\left\\{\alpha_{1},\dots,\alpha_{\left|\mathcal{A}\right|}\right\\}$.

As the noise power is normally unknown, so its estimation is also taken into
account. We define $\lambda$ as the precision of the noise, i.e.,
$\lambda=1/\sigma_{w}^{2}$. It is treated as a random variable with an
improper prior $f_{\lambda}(\lambda)\propto 1/\lambda$ [30]. With model (9),
we have the following joint conditional distribution and its factorization

$\displaystyle p(\boldsymbol{X},\boldsymbol{H},\lambda|\boldsymbol{Y})\propto
p(\boldsymbol{Y}|\boldsymbol{X},\boldsymbol{H},\lambda)p(\boldsymbol{X})p(\boldsymbol{H}|\boldsymbol{G})p(\boldsymbol{G})p(\lambda),$
(12)

where

$\displaystyle p(\boldsymbol{Y}|\boldsymbol{X},\boldsymbol{H},\lambda)$
$\displaystyle=$
$\displaystyle\mathcal{MN}\left(\boldsymbol{Y};\boldsymbol{H}\boldsymbol{X},\boldsymbol{I}_{M},\lambda^{-1}\boldsymbol{I}_{L}\right)$
(13) $\displaystyle p(\boldsymbol{H}|\boldsymbol{G})$ $\displaystyle=$
$\displaystyle\delta\left(\boldsymbol{H}-\boldsymbol{F}\boldsymbol{G}\right).$
(14)

We aim to compute the (proximate) marginals
$p(\boldsymbol{x}_{n}|\boldsymbol{Y})$, based on which the transmitted symbols
can be detected (at the same time, the channel matrix is also estimated).
However the exact inference is intractable as the involved variables are of
high dimension. In this work, we resort to approximate inference techniques,
and in particular, with a factor graph representation of the problem, a
message passing algorithm will be developed, leveraging the UAMP-MF algorithm.

##  IV  UAMP-MF Based Message Passing Algorithm Design

###  IV-A Factor Graph Representation

To facilitate the factor graph representation, we define the relevant local
functions based on the factorization in (12):
$f_{\boldsymbol{Y}}(\boldsymbol{Y},\boldsymbol{X},\boldsymbol{H},\lambda)=p(\boldsymbol{Y}|\boldsymbol{X},\boldsymbol{H},\lambda)$,
$f_{\mathcal{D}}(x_{n,l})=p(x_{n,l})$,
$f_{\boldsymbol{H}}(\boldsymbol{H},\boldsymbol{G})=p(\boldsymbol{H}|\boldsymbol{G})$,
$f_{\lambda}(\lambda)=p(\lambda)$, and
$f_{\boldsymbol{G}}(\boldsymbol{G})=p(\boldsymbol{G})$. The local factors and
their corresponding probability functions are summarized in Table 1. The
factor graph representation of (12) is shown in Fig. 3, where we divide the
factor graph into two sub-graphs and develop a message passing algorithm in
which UAMP-MF is incorporated. Specifically, the UAMP-MF algorithm is mainly
used to handle the message passing in Sub-graph I.

TABLE I: Local functions and distributions in (12) Local Function | Distribution | Function
---|---|---
$f_{\boldsymbol{Y}}(\boldsymbol{Y},\boldsymbol{X},\boldsymbol{H},\lambda)$ | $p(\boldsymbol{Y}|\boldsymbol{X},\boldsymbol{H},\lambda)$ | $\mathcal{MN}\left(\boldsymbol{Y};\boldsymbol{H}\boldsymbol{X},\boldsymbol{I}_{M},\lambda\boldsymbol{I}_{L}\right)$
$f_{\boldsymbol{H}}(\boldsymbol{H},\boldsymbol{G})$ | $p(\boldsymbol{H}|\boldsymbol{G})$ | $\delta\left(\boldsymbol{H}-\boldsymbol{F}\boldsymbol{G}\right)$
$f_{\mathcal{D}}(x_{n,l})$ | $p(x_{n,l})$ | $\frac{1}{|\mathcal{A}|}\sum_{a=1}^{|\mathcal{A}|}\delta(x_{n,l}-\alpha_{a})$
$f_{\boldsymbol{G}}(\boldsymbol{G})$ | $p(\boldsymbol{G})$ | $\prod_{n,k}\left((1-\epsilon)\delta(g_{kn})+\epsilon\mathcal{N}(g_{n,l};0,\nu)\right)$
$f_{\lambda}(\lambda)$ | $p(\lambda)$ | $1/\lambda$
Figure 3: Factor graph representation of equation (12).

###  IV-B Message Passing Algorithm Design

In the derivation of the message passing algorithm, we use the notation
$m_{f_{A\to B}}(x)$ to represent the message from (factor or variable) node
$A$ to node $B$, which is a function of $x$. It is noted that the message
passing algorithm is an iterative one, where each iteration involves a forward
recursion and a backward recursion, and the computation of a message in an
iteration may require some messages computed in the last iteration. The UAMP-
MF based algorithm developed in this work is summarized in Algorithm 3.

####  IV-B1 Message Update in Sub-Graph I

According to the derivation of UAMP-MF in [24], we can obtain the backward
message about $\boldsymbol{X}$ by running Lines 1-9 of Algorithm 2, which
produces two matrices $\boldsymbol{Q}_{X}$ and $\boldsymbol{V}_{Q_{X}}$. Due
to the decoupling of (U)AMP, we have the message about each entry of
$\boldsymbol{X}$, i.e.,

$m_{x_{nl}\to
f_{\mathcal{D}}}=\mathcal{CN}(x_{nl};q_{x_{nl}},\nu_{q^{x}_{nl}})$ (15)

where $q_{x_{nl}}$ and $\nu_{q^{x}_{nl}}$ are respectively the $(n,l)$-th
element of $\boldsymbol{Q}_{X}$ and $\boldsymbol{V}_{Q_{X}}$. This facilitates
the computation of the belief (marginal) of $x_{n,l}$ with the separable prior
(11), i.e.,

$b(x_{n,l})\propto\mathcal{CN}(x_{nl};q_{x_{nl}},\nu_{q^{x}_{nl}})\sum_{a=1}^{|\mathcal{A}|}\delta(x_{n,l}-\alpha_{a}).$
(16)

To make the message passing tractable, we then project $b(x_{n,l})$ to be
Gaussian with mean $\hat{x}_{nl}$ and variance $\nu_{x_{nl}}$. This is
equivalent to perform the minimum mean squared error estimation (MMSE), i.e.,
compute the a posteriori mean and variance of $x_{nl}$, with the following
pseudo scalar observation models

$q_{x_{nl}}=x_{nl}+w_{x_{nl}},n=1,...,N,l=1,...,L$ (17)

where $w_{x_{nl}}$ represents a Gaussian model noise with mean zero and
variance $\nu_{q^{x}_{nl}}$. It is not hard to show that the _a posteriori_
mean $\hat{x}_{nl}$ and variance $\nu_{x_{nl}}$ of $x_{nl}$ are given by

$\displaystyle\hat{x}_{nl}$ $\displaystyle=$
$\displaystyle\sum\nolimits_{a=1}^{|\mathcal{A}|}\alpha_{a}\beta_{nl}^{a}$
(18) $\displaystyle\nu_{x_{nl}}$ $\displaystyle=$
$\displaystyle\sum\nolimits_{a=1}^{|\mathcal{A}|}\beta_{nl}^{a}|\alpha_{a}-\hat{x}_{nl}|^{2},$
(19)

where

$\beta_{nl}^{a}={\xi_{nl}^{a}}/{\sum\nolimits_{a=1}^{|\mathcal{A}|}\xi_{nl}^{a}}$
(20)

with

$\xi_{nl}^{a}=\exp\left({-\nu_{q^{x}_{nl}}^{-1}|\alpha_{a}-q_{x_{nl}}|^{2}}\right).$
(21)

The above lead to Lines 1-5 of Algorithm 3.

Then we can stack $\hat{x}_{nl}$ and $\nu_{x_{nl}}$ to form matrices
$\hat{\boldsymbol{X}}=[\hat{x}_{nl}]$ and
$\boldsymbol{\Xi}_{X}=[\nu_{x_{nl}}]$. We can see that each element $x_{n,l}$
has its own variance. To facilitate subsequent processing, we make an
approximation by assuming that the elements of each row in
$\hat{\boldsymbol{X}}$ share a same variance, which is the average of their
variances. Then they can be collectively characterized by a matrix normal
distribution, i.e.,
$b(\boldsymbol{X})=\mathcal{MN}(\boldsymbol{X};\hat{\boldsymbol{X}},\boldsymbol{U}_{X},\boldsymbol{V}_{X})$
with $\boldsymbol{U}_{X}=\text{diag}(\text{mean}(\boldsymbol{\Xi}_{X},2))$ and
$\boldsymbol{V}_{X}=\boldsymbol{I}_{L}$, where
$\text{mean}(\boldsymbol{\Xi}_{X},2)$ represents the average operation on the
rows of $\boldsymbol{\Xi}_{X}$. These lead to Lines 6 and 7 of Algorithm 3.

###  IV-C Message Update in Sub-Graph II

Algorithm 3 UAMP-MF Based Blind Joint UACESD

Initialization: $\hat{N}=U_{r}$, $\boldsymbol{U}_{H}=\boldsymbol{I}_{M}$,
$\hat{\boldsymbol{H}}=\mathcal{R}_{c}(M,\hat{N})$,
$\boldsymbol{V}_{X}=\boldsymbol{I}_{L}$,
$[\boldsymbol{U}_{F},\boldsymbol{\Lambda}_{F},\boldsymbol{V}_{F}]=\text{SVD}(\boldsymbol{F})$,
$\boldsymbol{\Phi}_{G}=\boldsymbol{U}^{H}_{F}\boldsymbol{F}$, and
$\boldsymbol{\Lambda}_{G}=(\boldsymbol{\Lambda}_{F}\boldsymbol{\Lambda}_{F}^{\textrm{H}})\boldsymbol{1}_{K}\boldsymbol{1}_{N}^{\textrm{T}}$.
Repeat

1:Compute $\boldsymbol{V}_{Q_{X}}$ and $\boldsymbol{Q}_{X}$ with Lines 1-9 in
Algorithm 2

2:$\forall
n,l,a:\xi_{nl}^{a}=\exp\left({-\nu_{q^{x}_{nl}}^{-1}|\alpha_{a}-q_{x_{nl}}|^{2}}\right)$

3:$\forall
n,l,a:\beta_{nl}^{a}={\xi_{n,l}^{a}}/{\sum\nolimits_{a=1}^{|\mathcal{A}|}\xi_{nl}^{a}}$

4:$\forall
n,l:\hat{x}_{nl}=\sum\nolimits_{a=1}^{|\mathcal{A}|}\alpha_{a}\beta_{nl}^{a}$

5:$\forall
n,l:\nu_{x_{nl}}=\sum\nolimits_{a=1}^{|\mathcal{A}|}\beta_{nl}^{a}|\alpha_{a}-\hat{x}_{nl}|^{2}$

6:$\hat{\boldsymbol{X}}=\big{[}\hat{x}_{nl}\big{]}$,
$\boldsymbol{\Xi}_{X}=\big{[}\nu_{x_{nl}}\big{]}$

7:$\boldsymbol{U}_{X}=\text{diag}(\text{mean}(\boldsymbol{\Xi}_{X},2))$

8:Compute $\boldsymbol{V}_{Q_{H}}$ and $\boldsymbol{Q}_{H}$ with Lines 13 \-
21 in Algorithm 2

9:$\boldsymbol{R}_{G}=\boldsymbol{U}^{H}_{F}\boldsymbol{Q}_{\boldsymbol{H}}^{\textrm{H}}$,

10:$\tau=\left<1./\boldsymbol{V}_{Q_{H}}^{\textrm{T}}\right>$

11:$\boldsymbol{V}_{P_{G}}=\boldsymbol{\Lambda}_{G}\cdot\boldsymbol{V}_{G}$

12:$\boldsymbol{P}_{G}=\boldsymbol{\Phi}_{G}\hat{\boldsymbol{G}}-\boldsymbol{V}_{P_{G}}\cdot\boldsymbol{S}_{G}$

13:$\boldsymbol{V}_{S_{G}}=1./(\boldsymbol{V}_{P_{G}}+\tau^{-1})$

14:$\boldsymbol{S}_{G}=\boldsymbol{V}_{S_{G}}\cdot(\boldsymbol{R}_{G}-\boldsymbol{P}_{G})$

15:$1./\boldsymbol{V}_{Q_{G}}=\boldsymbol{1}_{K}\left(\frac{1}{K}\boldsymbol{1}_{M}^{T}(\boldsymbol{\Lambda}_{G}\cdot\boldsymbol{V}_{S_{G}})\right)$

16:$\boldsymbol{Q}_{G}=\hat{\boldsymbol{G}}+\boldsymbol{V}_{Q_{G}}\cdot(\mathbf{\Phi}_{G}^{H}\boldsymbol{S}_{G})$

17:$\forall n,k$ update $\hat{g}_{kn}$ and $\nu_{g_{kn}}$ by (36)

18:stack $\hat{g}_{kn}$ and $\nu_{g_{kn}}$ into matrices
$\hat{\boldsymbol{G}}$ and $\boldsymbol{V}_{G}$

19:update sparsity rate $\epsilon$ by (45)

20:$\boldsymbol{\Xi}_{H}=\left|\boldsymbol{U}_{F}\right|^{2}\left(1./(\tau+1./\boldsymbol{V}_{P_{G}})\right)$

21:$\hat{\boldsymbol{H}}=\boldsymbol{U}_{F}\left(\boldsymbol{V}_{H}\cdot(\boldsymbol{R}_{G}\cdot\tau+\hat{\boldsymbol{P}}_{G}./\boldsymbol{V}_{P_{G}})\right)$

22:$\boldsymbol{V}_{H}=\text{diag}(\text{mean}(\boldsymbol{\Xi}_{H},1))$

23:$\hat{\lambda}=ML/C$ with $C$ given in (5)

Until terminated

According to the derivations of the UAMP-MF algorithm, the outgoing message of
Sub-graph I is about matrix $\boldsymbol{H}$, which is also decoupled Gaussian
distributions with mean matrix
$\boldsymbol{Q}_{\boldsymbol{H}}\in\mathbb{C}^{M\times N}$ and variance matrix
$\boldsymbol{V}_{Q_{H}}\in\mathbb{C}^{M\times N}$ (the entries in the matrix
are the variances of the corresponding entries in $\boldsymbol{H}$), which can
be computed using Lines 13-21 of Algorithm 2.

Next, we need to compute the message passed from the factor node $f_{H}$ to
variable node $\boldsymbol{G}$. Considering the hard constraint
$f_{H}(\boldsymbol{H},\boldsymbol{G})=\delta(\boldsymbol{H}-\boldsymbol{F}\boldsymbol{G})$,
matrix $\boldsymbol{Q}_{\boldsymbol{H}}$ can be regarded as an observation
matrix with the following pseudo model

$\displaystyle\boldsymbol{Q}_{\boldsymbol{H}}=\boldsymbol{F}\boldsymbol{G}+\boldsymbol{W}_{H}$
(22)

where $\boldsymbol{W}_{H}$ is a model noise matrix, whose entries are
independent with means zero and variances given by the corresponding entries
in $\boldsymbol{V}_{Q_{H}}$. To enable low complexity inference, we make an
approximation that the model noise is white, and the entries have a same
precision $\tau$ (or variance $\tau^{-1}$), which is the average of the
individual precisions in the matrix $1./\boldsymbol{V}_{Q_{H}}$, i.e.,
$\tau=\left<1./\boldsymbol{V}_{Q_{H}}^{\textrm{T}}\right>$. After doing this,
(22) becomes a standard linear problem with known system transfer matrix
$\boldsymbol{F}$, hence UAMP in Algorithm 1 can be incorporated to solve it.

We choose UAMPv2 due to its lower complexity. It is noted that UAMPv2 in
Algorithm 1 is developed for a single vector problem. Next, we extend it to
the case for a matrix problem in (22). With the SVD
$\boldsymbol{F}=\boldsymbol{U}_{F}\boldsymbol{\Lambda}_{F}\boldsymbol{V}_{F}$,
we carry out a unitary transformation with $\boldsymbol{U}^{H}_{F}$ to (22),
yielding

$\boldsymbol{R}_{G}=\boldsymbol{\Phi}_{G}\boldsymbol{G}+\boldsymbol{W}_{H}^{\prime},$
(23)

where

$\displaystyle\boldsymbol{R}_{G}=\boldsymbol{U}^{H}_{F}\boldsymbol{Q}_{\boldsymbol{H}}\in\mathbb{C}^{M\times
N},$ (24)
$\displaystyle\boldsymbol{\Phi}_{G}=\boldsymbol{U}^{H}_{F}\boldsymbol{F}\in\mathbb{C}^{M\times
K},$ (25)

and $\boldsymbol{W}_{H}^{\prime}=\boldsymbol{U}^{H}_{F}\boldsymbol{W}_{H}$.
Since $\boldsymbol{F}$ is a unitary matrix, $\boldsymbol{W}_{H}^{\prime}$ is
still a white Gaussian noise matrix and has the same variance as
$\boldsymbol{W}_{H}$. Matrix $\boldsymbol{\Lambda}_{F}$ is an $M\times K$
rectangular diagonal matrix. In Algorithm 1, a vector
$\boldsymbol{\lambda}=\boldsymbol{\Lambda}\boldsymbol{\Lambda}^{H}\boldsymbol{1}$
is defined. To make the subsequent expression concise, we define a matrix
$\boldsymbol{\Lambda}_{G}\in\mathbb{C}^{M\times N}$ as

$\displaystyle\boldsymbol{\Lambda}_{G}=\boldsymbol{\lambda}_{F}\boldsymbol{1}_{N}^{\textrm{T}}$
(26)

where
$\boldsymbol{\lambda}_{F}=(\boldsymbol{\Lambda}_{F}\boldsymbol{\Lambda}_{F}^{\textrm{H}})\boldsymbol{1}_{K}$.

Following UAMPv2 in Algorithm 1, we can compute the mean and variances about
the entries in matrix $\boldsymbol{G}$ in the following. Firstly, we compute
matrices $\boldsymbol{V}_{P_{G}}$ and $\boldsymbol{P}_{G}$ as

$\displaystyle\boldsymbol{V}_{P_{G}}$ $\displaystyle=$
$\displaystyle\boldsymbol{\Lambda}_{G}\cdot\boldsymbol{V}_{G}$ (27)
$\displaystyle\boldsymbol{P}_{G}$ $\displaystyle=$
$\displaystyle\boldsymbol{\Phi}_{G}\hat{\boldsymbol{G}}-\boldsymbol{V}_{P_{G}}\cdot\boldsymbol{S}_{G}$
(28)

where $\hat{\boldsymbol{G}}$ and $\boldsymbol{V}_{G}$ represent the mean
matrix and variance matrix. The matrices are obtained based on the a
posteriori distribution $b(\boldsymbol{G})$, which are computed in (39). It is
noted that, before the first iteration, we use the initialization
$\hat{\boldsymbol{G}}=\boldsymbol{0}_{K\times N}$ and
$\boldsymbol{V}_{G}=\boldsymbol{1}_{K\times N}$. In (28) matrix
$\boldsymbol{S}_{G}$ is initiated as 0 for the first iteration, and is updated
in $\eqref{eq:hatSG}$.

Then, we update intermediate matrices $\boldsymbol{V}_{S_{G}}$ and
$\boldsymbol{S}_{G}$ as

$\displaystyle\boldsymbol{V}_{S_{G}}=1./(\boldsymbol{V}_{P_{G}}+\tau^{-1}{\boldsymbol{1}_{M\times
N}})$ (29)
$\displaystyle\boldsymbol{S}_{G}=\boldsymbol{V}_{S_{G}}\cdot(\boldsymbol{R}_{G}-\boldsymbol{P}_{G}),$
(30)

where $\tau$ represents the model noise precision. According to Line 5 of
Algorithm 1, we can get $\boldsymbol{V}_{Q_{G}}$ and $\boldsymbol{Q}_{G}$ as

$\displaystyle 1./\boldsymbol{V}_{Q_{G}}$ $\displaystyle=$
$\displaystyle\boldsymbol{1}_{K}\left(\frac{1}{K}\boldsymbol{1}_{M}^{T}(\boldsymbol{\Lambda}_{G}\cdot\boldsymbol{V}_{S_{G}})\right)$
(31) $\displaystyle\boldsymbol{Q}_{G}$ $\displaystyle=$
$\displaystyle\hat{\boldsymbol{G}}+\boldsymbol{V}_{Q_{G}}\cdot(\mathbf{\Phi}_{G}^{H}\boldsymbol{S}_{G}).$
(32)

The above lead to Lines 9-16 of Algorithm 3.

Due to the decoupling of (U)AMP, we have the following scalar pseudo
observation models

$g_{kn}=q_{g_{kn}}+w_{kn}$ (33)

where $q_{g_{kn}}=\left[\boldsymbol{Q}_{G}\right]_{kn}$,
$\nu^{q}_{g_{kn}}=\left[\boldsymbol{V}_{Q_{G}}\right]_{kn}$, and $w_{kn}$
represents the Gaussian noise with variance $\nu^{q}_{g_{kn}}$. Then, we can
compute the a posteriori distributions of $\\{g_{kn}\\}$ with the pseudo
models and the priors $\\{p(g_{kn})\\}$. As discussed previously, matrix
$\boldsymbol{G}\in\mathbb{C}^{K\times N}$ is sparse and we employ the sparsity
promoting Bernoulli-Gaussian prior in (10). The a posteriori distribution
$b(g_{kn})$ can be obtained as

$\displaystyle
b(g_{kn})=\frac{p(g_{kn})\mathcal{CN}(g_{kn};q_{g_{kn}},\nu^{q}_{g_{kn}})}{\int_{g_{kn}}p(g_{kn})\mathcal{CN}(g_{kn};q_{g_{kn}},\nu^{q}_{g_{kn}})}$
$\displaystyle=\frac{(1-\epsilon)\mathcal{CN}(q_{g_{kn}};0,\nu^{q}_{g_{kn}})\delta(g_{kn})+\alpha_{kn}\mathcal{CN}(g_{kn};\gamma_{kn},\nu_{\gamma_{kn}})}{(1-\epsilon)\mathcal{CN}(\hat{q}_{g_{kn}};0,\nu^{q}_{g_{kn}})+\alpha_{kn}}$
$\displaystyle=(1-\pi_{kn})\delta(g_{kn})+\pi_{kn}\mathcal{CN}(g_{kn};\gamma_{kn},\nu_{\gamma_{kn}})$
(34)

where

$\displaystyle\gamma_{kn}=\frac{q_{g_{kn}}\nu}{\nu^{q}_{g_{kn}}+\nu},\nu_{\gamma_{kn}}=\frac{\nu^{q}_{g_{kn}}\nu}{\nu^{q}_{g_{kn}}+\nu}$
$\displaystyle\alpha_{kn}=\epsilon\mathcal{CN}(q_{g_{kn}};0,\nu^{q}_{g_{kn}}+\nu)$
$\displaystyle\pi_{kn}=\frac{\alpha_{kn}}{(1-\epsilon)\mathcal{CN}(q_{g_{kn}};0,\nu^{q}_{g_{kn}})+\alpha_{kn}}.$
(35)

To avoid intractable message passing due to the complex form of the belief
$b(g_{kn})$, we project $b(g_{kn})$ to a Gaussian form, i.e.,

$\displaystyle\text{Proj}\left(b(g_{kn})\right)=\mathcal{CN}(g_{kn};\hat{g}_{kn},\nu_{g_{kn}})$
(36)

where

$\displaystyle\hat{g}_{kn}$ $\displaystyle=$
$\displaystyle\pi_{nk}\gamma_{kn},$ (37) $\displaystyle\nu_{g_{kn}}$
$\displaystyle=$ $\displaystyle\pi_{nk}^{2}\nu_{\gamma_{kn}}.$ (38)

The above lead to Line 17 of Algorithm 3.

We stack $\hat{g}_{kn},\nu_{g_{kn}}\forall n,k$ into matrices
$\hat{\boldsymbol{G}}$ and $\boldsymbol{V}^{\prime}_{G}$,

$\displaystyle\hat{\boldsymbol{G}}=\\{\hat{g}_{kn}\\},\boldsymbol{V}^{\prime}_{G}=\\{\nu_{g_{kn}}\\}.$
(39)

According to Line 7 of Algorithm 1, we need to average the variance matrix
$\boldsymbol{V}^{\prime}_{G}$ column wise, i.e.,

$\boldsymbol{\nu}_{g}^{T}=\frac{1}{K}\boldsymbol{1}_{K}^{\textrm{T}}\boldsymbol{V}^{\prime}_{G}$
(40)

then we get variance matrix $\boldsymbol{V}_{G}$ by

$\boldsymbol{V}_{G}=\boldsymbol{1}_{K}\boldsymbol{\nu}_{g}^{T}.$ (41)

As the knowledge about the sparsity rate $\epsilon$ is often unknown, we can
also incorporate its estimation into the message passing algorithm, with the
expectation maximization framework as in [34], i.e.,

$\epsilon=\mathop{\text{arg
max}}\limits_{\epsilon\in(0,1)}\sum_{k,n}\int_{g_{kn}}\Big{(}\ln
p(g_{kn}|\epsilon)b(g_{kn})\Big{)},$ (42)

and the value of $\epsilon$ can be obtained by letting the derivative be zero,
i.e.,

$\displaystyle\frac{\partial}{\partial\epsilon}\sum_{k,n}\int_{g_{kn}}\Big{(}\ln
p(g_{kn}|\epsilon)b(g_{kn})\Big{)}$ $\displaystyle\ \ \ \
=\sum_{k,n}\int_{g_{kn}}\Big{(}\frac{\partial}{\partial\epsilon}\ln
p(g_{kn}|\epsilon)b(g_{kn})\Big{)}=0$ (43)

With the Bernoulli-Gaussian distribution,

$\displaystyle\frac{\partial}{\partial\epsilon}\ln p(g_{kn}|\epsilon)$
$\displaystyle=$
$\displaystyle\frac{-\delta(g_{kn})+\mathcal{N}(g_{kn};0,\nu)}{(1-\epsilon)\delta(g_{kn})+\epsilon\mathcal{N}(g_{kn};0,\nu)}$
(44) $\displaystyle=$ $\displaystyle\begin{cases}1/\epsilon,&g_{kn}\neq 0\\\
-1/(1-\epsilon),&g_{kn}=0\end{cases}$

Substituting (36) and (44) into (43), we can obtain that [34]

$\displaystyle\epsilon=\frac{1}{KN}\sum\nolimits_{k,n}\pi_{kn}.$ (45)

The above correspond to Lines 17-19 of Algorithm 3.

Finally, we can compute the a posteriori distribution $b(\boldsymbol{H})$,
whose entries are approximated to be independent and Gaussian with variances
and means given by the corresponding entries of the matrices

$\displaystyle\boldsymbol{\Xi}_{H}$ $\displaystyle=$
$\displaystyle\left|\boldsymbol{U}_{F}\right|^{2}\left(1./(\tau+1./\boldsymbol{V}_{P_{G}})\right),$
(46) $\displaystyle\hat{\boldsymbol{H}}$ $\displaystyle=$
$\displaystyle\boldsymbol{U}_{F}\left(\boldsymbol{V}_{H}\cdot(\boldsymbol{R}_{G}\cdot\tau+\hat{\boldsymbol{P}}_{G}./\boldsymbol{V}_{P_{G}})\right).$
(47)

Similar to the case for updating $b(\boldsymbol{X})$, to accommodate
$\\{b(h_{mn})\\}$ with a matrix normal distribution, we make an approximation
by assuming the entries in each column of $\boldsymbol{H}$ share a common
variance, which is the average of their variances. Then
$b(\boldsymbol{H})=\mathcal{MN}(\boldsymbol{H};\hat{\boldsymbol{H}},\boldsymbol{U}_{H},\boldsymbol{V}_{H})$
with $\boldsymbol{V}_{H}=\text{diag}(\text{mean}(\boldsymbol{\Xi}_{H},1))$ and
$\boldsymbol{U}_{H}=\boldsymbol{I}_{M}$, where
$\text{mean}(\boldsymbol{\Xi}_{H},1)$ represents the average operation on the
columns of $\boldsymbol{\Xi}_{H}$. After that the noise precision is updated.
The above lead to Lines 20-23 of Algorithm 3.

###  IV-D Estimation of the Number of Active Users

In the above, we assume the number of active users $N$ is known, which needs
to be estimated in a practical scenario. In the blind UACESD algorithm, the
number of active users, $N$, corresponds to the number of rows in
$\boldsymbol{X}$ or the number of columns in $\boldsymbol{H}$. We follow the
method in [35] to estimate $N$ in an iterative manner based on the singular
values of the estimated matrix $\hat{\boldsymbol{H}}$. Initially, we set it to
the maximum number of possible active users $\bar{N}$, and update its value in
an iterative manner. With the sorted singular values $\sigma_{n}$ of the
estimated matrix $\hat{\boldsymbol{H}}$, we compute the pairwise ratios
$R_{n}=\sigma_{n}/\sigma_{n+1}$. The rank $N$ is estimated based on the
largest ratio, i.e.,

$\displaystyle\hat{N}=\text{arg max}_{n}R_{n}$ (48)

which is the largest gap in successive singular values. However, this estimate
candidate is accepted only if this maximizing ratio exceeds the average ratio,

$\displaystyle\bar{R}=\frac{1}{\hat{N}-2}\sum\nolimits_{i\neq\hat{N}}R_{i}.$
(49)

If the candidate $\hat{N}$ is accepted, the matrices $\boldsymbol{H}$ and
$\boldsymbol{X}$ are pruned and then the blind UACESD algorithm is run. If the
candidate is not accepted, the blind UACESD algorithm is run for one more
iteration, and the estimation of $N$ is performed until a candidate $\hat{N}$
is accepted.

###  IV-E Handling Ambiguity

An inherent problem with the matrix factorization problem in (9) is the
ambiguity. Due to the strong constraints on matrices $\boldsymbol{H}$ and
$\boldsymbol{X}$, there is only a phase ambiguity. To solve the problem, we
have two approaches. One approach is to use differential modulations for all
users, so the phase ambiguity does not have impact on demodulation. The other
approach is that all users transmit few pilot symbols, which are known to the
AP, so that the ambiguity can be mitigated. As the number of pilot symbols can
be very small, the overhead can be ignored. In this paper, we use the first
approach in the simulations.

##  V Simulation Results

In this section, we evaluate the performance of the UAMP-MF based blind joint
UASD algorithm. We assume that the access point is equipped with $M=100$
antennas and serves $U=300$ users in an area. In most of the simulations, we
assume that the rate of active users is $10\%$ and $16.7\%$, i.e., 30 and 50
active users during a transmission block. The size of matrix $F$ is $100\times
150$, the number of paths per user is 10, and the signals arrive at the AP
with random angles. Differential QPSK modulation is employed. Both coded
system and uncoded system are considered. For the coded system, we use a rate
$1/2$ convolutional code with generator polynomial $[5,7]_{8}$. The simulation
results are obtained by averaging over $10^{5}$ Monte Carlo trials. The signal
to noise power ration (SNR) is defined as

$\text{SNR}=\frac{||\boldsymbol{H}\boldsymbol{X}||^{2}/MNL}{\sigma^{2}},$ (50)

where the numerator is the power of the signal per antenna per user, and
$\sigma^{2}$ is the power of noise. To examine the performance of user
activity detection, we define active user error rate (AER) as

$\displaystyle\text{AER}=\frac{\text{\\#of active users}-\text{\\#of active
users identified successfully}}{\text{\\#of active users}}.$

To the best of our knowledge, there are no existing algorithms in the
literature that can be used to solve the blind joint UACESD problem formulated
in this paper. Hence, in this section, we compare the joint blind UACESD
algorithm (named ’Blind-UACESD’) with some performance bounds. The first bound
is the performance of the system where the number of active users $N$ is
known, and joint channel estimation and signal detection is performed, which
is denoted by ’Blind-CESD’. The second one is the performance of the system
where the active user number $N$ is known and a frame of transmitted signal
$\boldsymbol{X}$ is used as pilot to estimate the channel, which is named as
’CESD’. In the last one, we assume the perfect channel state information,
i.e., the channel matrix $\boldsymbol{H}$ is known, so the bilinear problem is
degenerated to a linear problem. The system in this case is denoted by ’SD’.
The bit error rate (BER), frame error rate (FRE), the normalized mean squared
error (NMSE) of channel estimation, and AER are used to evaluate the
performance.

We first consider an uncoded system. Figures 4 and 5 show the BER and FER
performance of the proposed algorithm and the relevant performance bounds. We
set the number of active users to $N=30$ and 50, and the length of the frame
$L=200$. From Fig. 4 we can see that the BER performance of Blind-UACESD and
Blind-CESD are almost the same, and Blind-UACESD performs very well and there
is about 1dB performance gap compared with CESD and SD. We can also see that,
when the number of active users is increased to 50, Blind-UACESD exhibits
certain performance loss in the relatively lower SNR range. This is because it
is difficult to identify the number of active users in the low SNR range.
Figure 5 shows the FER performance of the system versus SNR, from which the
same observations as Fig. 4 are made. Figure 5 also shows the AER performance
of Blind-UACESD when the number of active users is 50. We see that the user
activity detection performs well, and AER was not observed in the simulations
when SNR is above 3dB.

Figure 4: BER performance of the system, where $L=200$ and $N=30$ and $50$.
Figure 5: FER performance in an uncoded system, where $L=200$ and $N=30$ and
50.

In Fig.6, we vary the number of active users from 10 to 60 and examine the
performance of the system when SNR =3dB and 4dB. It can be observed that, the
BER performance of the system increases slightly with the number of active
users. When SNR=3dB and the number of active users is 60, the algorithm
experiences difficulties in active user detection, thereby leading to higher
BER. In general, the system with the proposed Blind-UACESD performs well, and
its performance is close to those performance bounds. To examine the channel
estimation performance, the performance of NMSE versus SNR is shown in Fig. 7,
where the ambiguity is removed. We can be seen from the figure that, the NMSE
of channel estimation of Blind-UACESD is almost the same as that obtained
using pilot signals (i.e., CESD). Again, when the number of active users is
50, there is also performance loss in the low SNR range, compared with the
case when the number of active users is 30.

Figure 6: BER performance versus the number of active users, where SNR = 3dB
and 4dB. Figure 7: MSE performance of channel estimation in an uncoded system,
where $L=200$ and $N=30$ and 50.

Then we consider a coded system where all the users employ a rate-$1/2$
convolutional code with generator polynomial $[5,7]_{8}$. Figures 8 and 9 show
the BER, FER, and AER performance of the system, respectively. The frame
length $L=300$. We can see that, compared to the results in the uncoded case,
the performance of the coded system is improved considerably.

In the design of the Blind-UACESD algorithm, the variance (reciprocal of the
precision) of the noise is assumed unknown, and its estimation is incorporated
in the algorithm. Lastly, we examine the performance of noise variance
estimation and the results are shown in Fig. 10 where we vary the SNR. We can
see that the estimated noise variance matches the real noise variance very
well at different SNRs.

Figure 8: BER performance with a coded system, where $L=300$ and $N=30$ and
50. Figure 9: FER performance with a coded system, where $L=300$ and $N=30$.
Figure 10: Noise variance estimation versus SNR.

##  VI Conclusions

In this paper, we have investigated the issue of joint user activity
detection, channel estimation and signal detection without the use of pilot
signals in mmWave MIMO-based grant-free random access to further reduce the
overhead and latency. The blind joint user activity detection, channel
estimation and signal detection is formulated as a contained composite matrix
factorization problem, and we have shown that it can be solved by exploiting
the structures of the factoring matrices, i.e., the channel matrix is sparse
in the beam space domain and the entries of the data matrix are discrete
valued. Leveraging UAMP-MF, we have developed a message passing Bayesian
algorithm to solve the problem. Extensive simulation results are provided to
demonstrate the effectiveness of the grant-free random access scheme without
the use of pilot signals.

## References

  * [1] F. Boccardi, R. W. Heath, A. Lozano, T. L. Marzetta, and P. Popovski, “Five disruptive technology directions for 5G,” _IEEE Communications Magazine_ , vol. 52, no. 2, pp. 74–80, 2014.
  * [2] H. Tullberg, P. Popovski, Z. Li, M. A. Uusitalo, A. Hoglund, O. Bulakci, M. Fallgren, and J. F. Monserrat, “The metis 5G system concept: Meeting the 5G requirements,” _IEEE Communications Magazine_ , vol. 54, no. 12, pp. 132–139, 2016.
  * [3] M. T. Islam, A.-e. M. Taha, and S. Akl, “A survey of access management techniques in machine type communications,” _IEEE Communications Magazine_ , vol. 52, no. 4, pp. 74–81, 2014.
  * [4] Z. Chen, F. Sohrabi, and W. Yu, “Sparse activity detection for massive connectivity,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 7, pp. 1890–1904, 2018.
  * [5] J. Fu, G. Wu, Y. Zhang, L. Deng, and S. Fang, “Active user identification based on asynchronous sparse bayesian learning with SVM,” _IEEE Access_ , vol. 7, pp. 108 116–108 124, 2019.
  * [6] T. Hara and K. Ishibashi, “Blind multiple measurement vector amp based on expectation maximization for grant-free NOMA,” _IEEE Wireless Communications Letters_ , vol. 11, no. 6, pp. 1201–1205, 2022.
  * [7] X. Zhang, F. Labeau, L. Hao, and J. Liu, “Joint active user detection and channel estimation via bayesian learning approaches in MTC communications,” _IEEE Transactions on Vehicular Technology_ , vol. 70, no. 6, pp. 6222–6226, 2021.
  * [8] X. Zhang, P. Fan, J. Liu, and L. Hao, “Bayesian learning-based multiuser detection for grant-free NOMA systems,” _IEEE Transactions on Wireless Communications_ , vol. 21, no. 8, pp. 6317–6328, 2022.
  * [9] W. Chen, H. Xiao, L. Sun, and B. Ai, “Joint activity detection and channel estimation in massive MIMO systems with angular domain enhancement,” _IEEE Transactions on Wireless Communications_ , vol. 21, no. 5, pp. 2999–3011, 2022.
  * [10] Z. Zhang, Y. Li, C. Huang, Q. Guo, L. Liu, C. Yuen, and Y. L. Guan, “User activity detection and channel estimation for grant-free random access in leo satellite-enabled internet of things,” _IEEE Internet of Things Journal_ , vol. 7, no. 9, pp. 8811–8825, 2020.
  * [11] C. Wei, H. Liu, Z. Zhang, J. Dang, and L. Wu, “Approximate message passing-based joint user activity and data detection for NOMA,” _IEEE Communications Letters_ , vol. 21, no. 3, pp. 640–643, 2017.
  * [12] Y. Chi, L. Liu, G. Song, C. Yuen, Y. L. Guan, and Y. Li, “Message passing in c-ran: Joint user activity and signal detection,” in _GLOBECOM 2017 - 2017 IEEE Global Communications Conference_ , 2017, pp. 1–6.
  * [13] F. Wei, W. Chen, Y. Wu, J. Ma, and T. A. Tsiftsis, “Message-passing receiver design for joint channel estimation and data decoding in uplink grant-free scma systems,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 1, pp. 167–181, 2019.
  * [14] Y. Du, B. Dong, W. Zhu, P. Gao, Z. Chen, X. Wang, and J. Fang, “Joint channel estimation and multiuser detection for uplink grant-free NOMA,” _IEEE Wireless Communications Letters_ , vol. 7, no. 4, pp. 682–685, 2018.
  * [15] 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 Processing Magazine_ , vol. 35, no. 5, pp. 88–99, 2018.
  * [16] Y. Zhang, Q. Guo, Z. Wang, J. Xi, and N. Wu, “Block sparse bayesian learning based joint user activity detection and channel estimation for grant-free NOMA systems,” _IEEE Transactions on Vehicular Technology_ , vol. 67, no. 10, pp. 9631–9640, 2018.
  * [17] B. Wang, L. Dai, Y. Yuan, and Z. Wang, “Compressive sensing based multi-user detection for uplink grant-free non-orthogonal multiple access,” pp. 1–5, 2015\.
  * [18] B. Wang, L. Dai, Y. Zhang, T. Mir, and J. Li, “Dynamic compressive sensing-based multi-user detection for uplink grant-free NOMA,” _IEEE Communications Letters_ , vol. 20, no. 11, pp. 2320–2323, 2016.
  * [19] Y. Du, B. Dong, Z. Chen, X. Wang, Z. Liu, P. Gao, and S. Li, “Efficient multi-user detection for uplink grant-free NOMA: Prior-information aided adaptive compressive sensing perspective,” _IEEE Journal on Selected Areas in Communications_ , vol. 35, no. 12, pp. 2812–2828, 2017.
  * [20] T. Ding, X. Yuan, and S. C. Liew, “Sparsity learning-based multiuser detection in grant-free massive-device multiple access,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 7, pp. 3569–3582, 2019.
  * [21] J. Zhang, X. Yuan, and Y.-J. A. Zhang, “Blind signal detection in massive MIMO: Exploiting the channel sparsity,” _IEEE Transactions on Communications_ , vol. 66, no. 2, pp. 700–712, 2018.
  * [22] Z. Guo, X. Wang, and W. Heng, “Millimeter-wave channel estimation based on 2-d beamspace music method,” _IEEE Transactions on Wireless Communications_ , vol. 16, no. 8, pp. 5384–5394, 2017.
  * [23] M. Xiao, S. Mumtaz, Y. Huang, L. Dai, Y. Li, M. Matthaiou, G. K. Karagiannidis, E. Bjornson, K. Yang, C.-L. I, and A. Ghosh, “Millimeter wave communications for future mobile networks,” _IEEE Journal on Selected Areas in Communications_ , vol. 35, no. 9, pp. 1909–1935, 2017.
  * [24] Z. Yuan, Q. Guo, Y. C. Eldar, and Y. Li, “Unitary approximate message passing for matrix factorization,” _ArXiv_ , vol. abs/2208.00422, 2022.
  * [25] D. L. Donoho, A. Maleki, and A. Montanari, “Message passing algorithms for compressed sensing: I. motivation and construction,” pp. 1–5, 2010.
  * [26] S. Rangan, “Generalized approximate message passing for estimation with random linear mixing,” pp. 2168–2172, 2011.
  * [27] S. Rangan, P. Schniter, A. K. Fletcher, and S. Sarkar, “On the convergence of approximate message passing with arbitrary matrices,” _IEEE Transactions on Information Theory_ , vol. 65, no. 9, pp. 5339–5351, 2019.
  * [28] Q. Guo, D. D. Huang, S. Nordholm, J. Xi, and Y. Yu, “Iterative frequency domain equalization with generalized approximate message passing,” _IEEE Signal Process. Lett._ , vol. 20, no. 6, pp. 559–562, June. 2013.
  * [29] Q. Guo and J. Xi, “Approximate message passing with unitary transformation,” _arXiv preprint arXiv:1504.04799_ , Apr. 2015.
  * [30] Z. Yuan, Q. Guo, and M. Luo, “Approximate message passing with unitary transformation for robust bilinear recovery,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 617–630, 2021.
  * [31] M. I. Jordan, Z. Ghahramani, T. S. Jaakkola, and L. K. Saul, “An introduction to variational methods for graphical models,” _Mach. Learn._ , vol. 37, no. 2, pp. 183–233, Nov. 1999.
  * [32] D. Tse and P. Viswanath, _Fundamentals of Wireless Communication_. Cambridge University Press, 2005.
  * [33] D. Baron, S. Sarvotham, and R. G. Baraniuk, “Bayesian compressive sensing via belief propagation,” _IEEE Transactions on Signal Processing_ , vol. 58, no. 1, pp. 269–280, 2010.
  * [34] J. P. Vila and P. Schniter, “Expectation-maximization gaussian-mixture approximate message passing,” _IEEE Transactions on Signal Processing_ , vol. 61, no. 19, pp. 4658–4672, 2013.
  * [35] J. T. Parker, P. Schniter, and V. Cevher, “Bilinear generalized approximate message passing-part I: Derivation,” _IEEE Transactions on Signal Processing_ , vol. 62, no. 22, pp. 5839–5853, 2014.