# Reconfigurable Intelligent Surface-Aided 6G Massive Access: Coupled Tensor
Modeling and Sparse Bayesian Learning

Xiaodan Shao, , Lei Cheng, Xiaoming Chen, , Chongwen Huang, , and Derrick Wing
Kwan Ng Xiaodan Shao is with the College of Information Science and Electronic
Engineering, Zhejiang University, Hangzhou, 310027, China, and also with
Zhejiang Provincial Key Laboratory of Info. Proc., Commun. & Netw. (IPCAN),
Hangzhou 310027, China (Email: shaoxiaodan@zju.edu.cn).Lei Cheng is with the
College of Information Science and Electronic Engineering, Zhejiang
University, Hangzhou, 310027, China, and also with Zhejiang Provincial Key
Laboratory of Info. Proc., Commun. & Netw. (IPCAN), Hangzhou 310027, China
(e-mail: lei_cheng@zju.edu.cn).Xiaoming Chen is with the College of
Information Science and Electronic Engineering, Zhejiang University, Hangzhou
310027, China (Email: chen_xiaoming@zju.edu.cn).C. Huang is with College of
Information Science and Electronic Engineering, Zhejiang University, Hangzhou
310027, China, and with International Joint Innovation Center, Zhejiang
University, Haining 314400, China, and also with Zhejiang-Singapore Innovation
and AI Joint Research Lab and Zhejiang Provincial Key Laboratory of Info.
Proc., Commun. & Netw. (IPCAN), Hangzhou 310027, China. (E-mail:
chongwenhuang@zju.edu.cn).Derrick Wing Kwan Ng is with the School of
Electrical Engineering and Telecommunications, University of New South Wales,
Sydney, NSW 2052, Australia.(E-mail: w.k.ng@unsw.edu.au).

###### Abstract

This paper investigates a reconfigurable intelligent surface (RIS)-aided
unsourced random access (URA) scheme for the sixth-generation (6G) wireless
networks with massive sporadic traffic devices. First of all, this paper
proposes a novel joint active device separation (the message recovery of
active device) and channel estimation architecture for the RIS-aided URA.
Specifically, the RIS passive reflection is optimized before the successful
device separation. Then, by associating the data sequences to multiple rank-
one tensors and exploiting the angular sparsity of the RIS-BS channel, the
detection problem is cast as a high-order coupled tensor decomposition problem
without the need of exploiting pilot sequences. However, the inherent coupling
among multiple sparse device-RIS channels, together with the unknown number of
active devices make the detection problem at hand deviate from the widely-used
coupled tensor decomposition format. To overcome this challenge, this paper
judiciously devises a probabilistic model that captures both the element-wise
sparsity from the angular channel model and the low-rank property due to the
sporadic nature of URA. Then, based on such a probabilistic model, a iterative
detection algorithm is developed under the framework of sparse variational
inference, where each update iteration is obtained in a closed-form and the
number of active devices can be automatically estimated for effectively
avoiding the overfitting of noise. Extensive simulation results confirm the
excellence of the proposed URA algorithm, especially for the case of a large
number of reflecting elements for accommodating a significantly large number
of devices.

###### Index Terms:

Massive unsourced random access, higher-order coupled tensor, sparse Bayesian
learning, 6G, reconfigurable intelligent surface.

##  I Introduction

Massive access or massive machine-type communication (mMTC) has been
identified as one of the main application cases in the upcoming sixth-
generation (6G) wireless networks, where the number of potential devices is
expected to reach hundreds of billions by 2030. In particular, a large number
of devices are connected to the Internet via a base station (BS) with sporadic
communications, where only a small fraction of devices are active concurrently
[2]. In this context, applying conventional grant-based random access schemes
to this emerging scenario would lead to an exceedingly high access latency and
a prohibitive signaling overhead. To address this issue, grant-free random
access schemes have been regarded as a promising candidate technology for
enabling massive access in 6G wireless networks [3], where active devices
transmit their data without establishing any access grant protocol with the
BS. Currently, a commonly discussed category of grant-free random access is
sourced random access, where each device is preassigned with a unique non-
orthogonal pilot sequence and the BS identifies the device activity based on
the received signals by distinguishing the transmitted sequences [2]. However,
since the number of devices in 6G wireless network grows sufficiently large,
sourced random access schemes become increasingly inefficient to assign fixed
pilot sequences to all potential devices. For instance, in order to detect
$K_{a}$ active devices from total $\bar{K}$ devices, the length of pilot
sequence should grow with the scale of $K_{a}\log(\bar{K})$, even through an
effective sparse recovery algorithm implemented [4, 5].

To overcome this challenge, another category of grant-free random access,
i.e., unsourced random access (URA), was proposed in [6]. Indeed, the roll-out
of URA is motivated by practical Internet-of-Things scenarios, where millions
of low-cost devices have a common codebook at the moment of production [6].
Active devices maps their messages to codewords of the common codebook, and
the BS recovers the messages by detecting the transmitted codewords, without
determining the identities of active devices. If active devices wish to reveal
themselves, they can embed identity information in their payloads. In order to
decrease the size of codebook, active devices usually partition their messages
to several sub-messages [8, 7]. Active devices adds some parity check bits to
each sub-message which is mapped to a codeword of the common codebook, and
then the BS recovers all the sub-messages which are stitched to the original
messages by using a some decoder according to the parity check bits, e.g., the
tree-based decoder [8]. Compared to sourced random access, URA has a higher
spectral efficiency, since it does send pilot sequence for identity detection
[9]. Moreover, URA has a lower computational complexity, which is independent
of the total number of potential devices and depends only on the cardinality
of the active device set [9]. Hence, URA is appealing to 6G wireless networks
with a large number of devices with sporadic activities.

Despite the advancement of URA, the devices may locate in a service dead zone,
where the direct link from the devices to the BS is not available [10].
Consequently, the signals received from the devices are weak such that
performing active device separation is a challenging task for the BS. To
overcome this challenge, the use of millimeter-wave (mmWave)/terahertz (THz)
frequencies with massive multiple-input-multiple output (MIMO) was proposed to
enhance the detection accuracy of grant-free random access [11]. However, it
significantly increases the hardware and energy cost as well as signal
processing complexity due to the deployment of a large number of radio
frequency (RF) chains. To this end, emerging wireless technologies for
manipulating the wireless channels with low cost, low complexity, and high
energy efficiency are desired to unlock the potential of URA in 6G wireless
networks.

Recently, reconfigurable intelligent surface (RIS), as a promising technology
of 6G wireless networks, has been proposed to enhance the signal quality at
desired receivers [12, 13, 14, 15, 16, 17]. Specifically, an RIS is able to
establish favourable channel responses by customizing the wireless propagation
environment via a large number of reconfigurable passive reflecting elements.
Thus, the deployment of RIS can enhance the strength of the signals between
the BS and the devices in dead zones via bypassing the obstacle between them.
Moreover, by optimizing the phase shifts, the signal propagation between the
BS and the devices can be smartly reconfigured to improve the accuracy of
active device separation. On the other hand, the performance gain in RIS-aided
communication systems depends critically on the availability of channel state
information (CSI). Yet, CSI acquisition is quite challenging in practice due
to the passive nature of RIS. In general, RIS without equipping active radio
frequency chains can neither transmit nor receive pilot signals, thus it is
challenging to estimate the channels from RIS to the BS and devices to RIS
separately. Instead, the concatenated device-RIS-BS channels are usually
estimated based on the pilot sequences sent from the devices. For example, by
exploiting the information on the slow-varying channel components and the
hidden channel sparsity, [18] formulated the channel estimation problem for
RIS-aided multiuser MIMO systems as a matrix-calibration based matrix
factorization task. Besides, to enable massive random access, [10] proposed a
three-stage framework based on the approximate message passing (AMP) for joint
active device separation and channel estimation in RIS-aided sourced random
access.

Although various approaches have been proposed in the literature for channel
estimation in RIS-aided systems, e.g., [13, 14, 16, 18, 10], the required
pilot signaling overhead scales with the product of the number of RIS
reflecting elements and devices, which are prohibitively large in practical
scenarios, especially for the case of massive access. As a remedy, the study
of joint active device separation and channel estimation with a fewer number
of pilots for the RIS-aided URA scheme in 6G wireless networks is desired. To
achieve the goal, the considered problem can be formulated as a novel non-
standard coupled tensor decomposition problem with a sparse coupled factor,
assuming no knowledge about the tensor rank (i.e., the number of active
devices). To handle coupled tensor decomposition type problems, various
attempts have been conducted [19, 20]. For instance, [19] proposed a coupled
polyadic decomposition method, where each tensor is a sum of coupled rank-1
tensors sharing some common factors. Unfortunately, the results from [19] and
[20] cannot be directly applied to the joint active device separation and
channel estimation problem with a complicated coupled structure. Furthermore,
the number of active devices is unknown in practice, which determines the
number of unknown model parameters. Besides, it is generally non-deterministic
polynomial-time hard (NP-hard) to acquire the number of active devices from
tensor data. To address this problem, the authors in [21] introduced an
additive regularization term that penalizes complicated channel models to
mitigate the overfitting of noise. However, this method inevitably consumes
enormous computational resources and leads to a heavy computation burden. To
learn the tensor rank automatically, recent advances integrate the tensor rank
learning into its hyper-parameter inference steps, such that the Bayesian
theory can provide a natural recipe for automatic rank determination [22, 23,
24, 25]. Nevertheless, the effectiveness of applying variational Bayesian
inference for solving coupled tensor decomposition is still unclear in the
literature. Moreover, the coupling structure (from the RIS-aided URA scheme
[6]), the element-wise sparsity structure (from angular channel representation
[26]), and the low-rank property (from the sporadic nature of mMTC [3])
together pose unique challenges in designing proper probabilistic model and
the associated inference algorithm. In this paper, by devising novel priors to
incorporate problem structures, we design a novel probabilistic model to
investigate the trade-off between the expressive power of the model and the
tractability of the inference algorithm. The contributions of this paper are
as follows:

  1. 1.

This paper studies a new application of RIS for detecting active devices and
estimating the corresponding channels to support URA in 6G wireless networks.
Specifically, we propose a new RIS-aided URA architecture, where RIS passive
reflection is first optimized for enhancing the performance of device
separation. Then the transmitted data symbols are constructed as multiple
rank-1 tensors that allows a pragmatic method for active device separation and
device-RIS channel estimation without relying on pilot sequences.

  2. 2.

This paper proposes a novel element-wise sparsity and low-rank inducing
probabilistic model for the considered high-order coupled tensor decomposition
problem. By using variational inference, an intelligent joint active device
separation and channel estimation algorithm is derived, assuming no knowledge
of active device number and noise power.

  3. 3.

This paper provides a thorough theoretical analysis of the proposed algorithm
and the results show that it enjoys fast convergence and low computational
complexity. Importantly, the proposed algorithm is well designed for handling
URA in large-scale RIS regime compared with the existing ones.

The rest of this paper is organized as follows: Section II introduces the
system model and the channel models, and Section III proposes the RIS-aided
joint active device separation and channel estimation architecture. Section IV
designs a coupled tensor based automatic detection algorithm by applying
Bayesian learning and analyzes the performance of the proposed algorithm.
Section V provides extensive simulation results to validate the effectiveness
of the proposed algorithm. Finally, Section VI concludes the paper.

_Notations_ : We use $\mathbb{C}^{A\times B}$ to denote the space of complex
matrices of size $A\times B$, $|\cdot|$ to denote the absolute value of a
complex number, $(\cdot)^{H}$ and $(\cdot)^{T}$ to denote conjugate transpose
and transpose, respectively, $\otimes$ to denote the Kronecker product,
$\circ$ to denote vector outer product, $\left\|\cdot\right\|_{F}$ to denote
Frobenius norm of a matrix, $\left\|\cdot\right\|_{2}$ to denote the
$l_{2}$-norm of an input vector, $\left[\\!\left[\cdot\right]\\!\right]$ to
denote the Kruskal operator. $x\in[a,b]$ means that scalar $x$ lies in the
closed interval between $a$ and $b$. $\mathbb{E}$ denotes the expectation of
its input argument. $\|\cdot\|_{0}$ denotes the $l_{0}$-norm defined as the
number of nonzero elements of an input vector or matrix. $*$ denotes
conjugation. The Khatri-Rao product is denoted by $\diamond$. The Hadamard
product is denoted by $\odot$. $\mathbf{I}_{K}$ represents the $K\times K$
identity matrix. Vector $\mathbf{1}_{K}\in\mathbb{R}^{K}$ is the all-ones
vector with length $K$. $\mathbf{a}(n)$ denotes the $n$th element of vector
$\mathbf{a}$. $\mathbf{A}(n,k)$ denotes the $(n,k)$th element of matrix
$\mathbf{A}$. $\mathbf{A}(:,n)$ and $\mathbf{A}(n,:)$ denote the $n$th column
and the $n$th row of matrix $\mathbf{A}$, respectively. $\mathrm{Tr}(\cdot)$
denotes the trace of a matrix. $\mathbf{x}\in\mathbb{C}^{n}$ is said to follow
a vector-valued normal distribution with mean $\mathbf{u}$ and covariance
matrix $\bm{\Sigma}$ of the form
$\mathcal{CN}_{n}(\mathbf{x}|\mathbf{u},\bm{\Sigma})$. Matrix $\mathbf{G}$ is
said to follow the matrix-variate normal distribution with mean matrix
$\mathbf{M}$ and covariance matrix $\bm{\Upsilon}\otimes\bm{\Sigma}$ of the
form $\mathcal{CN}_{N\times
K}(\mathbf{G}|\mathbf{M},\bm{\Upsilon}\otimes\bm{\Sigma})=\frac{1}{(2\pi)^{NK/2}|\bm{\Upsilon}|_{d}^{N/2}|\bm{\Sigma}|_{d}^{K/2}}\exp(-\frac{1}{2}\mathrm{Tr}\bm{\Upsilon}^{-1}(\mathbf{G}-\mathbf{M})\bm{\Sigma}^{-1}(\mathbf{G}-\mathbf{M})^{H})$
with $|\cdot|_{d}$ being the determinant. $\mathcal{O}(\cdot)$ stands for the
big-O notation. $p(\cdot|\cdot)$ denotes conditional probability distribution.
$\mathrm{diag}(\mathbf{x})$ denotes a diagonal matrix with the diagonal
entries specified by vector $\mathbf{x}$. $\mathcal{D}[\mathbf{A}]$ is a
diagonal matrix taking the diagonal element from $\mathbf{A}$.

##  II System Model

We consider a RIS-aided 6G wireless network, as illustrated in Fig. 1, where a
BS equipped with $M$ antennas serves massive single-antenna devices. Herein,
the RIS consisting of $N$ passive reflecting elements is deployed to enhance
the devices’ communication performance. The passive reflecting elements of the
RIS are arranged as an $N_{1}\times N_{2}$ uniform rectangular array with
$N\triangleq N_{1}\times N_{2}$. Although the number of potential devices
$\bar{K}$ is numerous, only $K_{a}\ll\bar{K}$ devices are active concurrently
at a given slot, due to the bursty nature of the envisioned 6G services [3].
Note that the number of active devices, $K_{a}$, is random and unknown to the
BS. It is assumed that the direct links between active devices and the BS are
obstructed [10], and thus calls for a need to enhance the quality of
communications via deploying RIS [13]. The channels are assumed to be
invariant in a given time slot, but might independently fade over time slots.
In particular, $\mathbf{h}_{k}\in\mathbb{C}^{N}$ is used to denote the channel
vector between the RIS and the $k$th device; $\mathbf{U}\in\mathbb{C}^{M\times
N}$ is used to denote the channel matrix between the BS and the RIS.

Figure 1: RIS-aided massive unsourced random access in 6G wireless networks
with blocked direct links.

###  II-A Channel Models

Since the RIS is usually mounted at some tall building, it is expected that
there are only limited scatters around the RIS. This suggests the adoption of
a sparse angular channel model [26], which has been validated by real-world
measurements and widely adopted in RIS related research [18]. In specific, the
RIS-BS channel, i.e., ${\bf U}\in\mathbb{C}^{M\times N}$ can be expressed as

$\displaystyle\mathbf{U}=\sqrt{\bar{\mu}}\sum_{i=1}^{\bar{I}}\bar{\xi}_{i}\mathbf{a}_{\mathrm{B}}(\sigma_{i})\mathbf{a}_{\mathrm{R}}(\bar{\phi}_{i},\bar{\psi}_{i})^{H},$
(1)

where
$\mathbf{a}_{\mathrm{R}}(\bar{\phi}_{i},\bar{\psi}_{i})=\bm{\varphi}_{N_{2}}(\cos(\bar{\psi}_{i})\cos(\bar{\phi}_{i}))\otimes\bm{\varphi}_{N_{1}}(\cos(\bar{\psi}_{i})\sin(\bar{\phi}_{i}))$
with $\bm{\varphi}_{N_{p}}(x)\triangleq\frac{1}{\sqrt{N_{p}}}[1,\\\
e^{-j\frac{2\pi}{\varrho}dx},\cdots,e^{-j\frac{2\pi}{\varrho}d(N_{p}-1)x}]^{T},\forall
p=1,2$; $(\bar{\phi}_{i},\bar{\psi}_{i})$ being the azimuth and elevation
angle-of-departure (AoD) from the RIS, respectively;
$\mathbf{a}_{\mathrm{B}}(\sigma_{i})=\bm{\varphi}_{M}(\sin(\sigma_{i}))$ is
the receive array response with respect to AoDs from the RIS to the BS with
$\|\mathbf{a}_{\mathrm{B}}(\sigma_{i})\|^{2}=M$; $\sigma_{i}$ is the
corresponding azimuth angle-of-arrival (AoA) to the BS; $N_{1}$ and $N_{2}$
denote the length and the width of the rectangular array of RIS; $\varrho$
denotes the wavelength of carrier frequency; $\bar{I}$ is the total number of
spatial paths between the BS and the RIS; $\bar{\mu}$ is the RIS-BS distance-
dependent path loss; $\bar{\xi}_{i}$ denotes the complex amplitude associated
with the $i$-th path.

For the device-RIS channel $\mathbf{h}_{k}\in\mathbb{C}^{N}$, it can be
modeled as

$\displaystyle\mathbf{h}_{k}=\sqrt{\mu_{k}}\sum_{i=1}^{I_{k}}\epsilon_{i}^{k}\mathbf{a}_{\mathrm{R}}(\phi_{i}^{k},\psi_{i}^{k}),$
(2)

where $\epsilon_{i}^{k}$ denotes the complex amplitude associated with the
$i$-th path; $\mu_{k}$ is the device-RIS distance-dependent path loss;
$\phi_{i}^{k}$ and $\psi_{i}^{k}$ are the corresponding azimuth and elevation
AoA to the RIS, respectively; $I_{k}$ denotes the number of paths between the
$k$-th device and the RIS. Following the grid-based scheme in [18], the
representation of channel vector $\mathbf{h}_{k}$ can be further simplified.
In particular, two sampling grids
$\bm{\nu}=[\nu_{1},\cdots,\nu_{{N_{1}}^{{}^{\prime}}}]^{T}$ with length
${N_{1}}^{{}^{\prime}}\geq N_{1}$, and
$\bm{\varsigma}=[\varsigma_{1},\cdots,\varsigma_{{N_{1}}^{{}^{\prime}}}]^{T}$
with length ${N_{2}}^{{}^{\prime}}\geq N_{2}$, are employed such that
$\mathbf{h}_{k}$ can be represented as

$\displaystyle\mathbf{h}_{k}=\mathbf{A}_{\mathrm{R}}\bm{\lambda}_{k},$ (3)

where
$\mathbf{A}_{\mathrm{R}}\\!=\\!\left[\bm{\varphi}_{N_{1}}(\nu_{1}),\cdots,\bm{\varphi}_{N_{1}}(\nu_{{N_{1}}^{{}^{\prime}}})\right]\\!\otimes\\!\left[\bm{\varphi}_{N_{2}}(\varsigma_{1}),\cdots,\bm{\varphi}_{N_{2}}(\varsigma_{{N_{2}}^{{}^{\prime}}})\right]\in\mathbb{C}^{N\times
N_{1}^{\prime}N_{2}^{\prime}}$, and
$\bm{\lambda}_{k}\in\mathbb{C}^{N_{1}^{\prime}N_{2}^{\prime}}$ represents the
channel coefficients of $\mathbf{h}_{k}$ in the angular domain. Since the
number of paths is usually limited, $\bm{\lambda}_{k}$ is essentially sparse,
with each nonzero value being the corresponding composite path gain
$\sqrt{\mu_{k}}\epsilon_{i}^{k}$. The exploitation of such a sparse structure
has been shown to bring improved channel estimation performance and reduced
complexity [11, 26, 27, 28].

###  II-B Signal Models

According to the characteristics of RIS-aided URA, we propose a joint active
device separation and channel estimation frame structure, as shown in Fig. 2.
Specifically, the BS first estimates the RIS-BS channel and then recovers the
transmitted messages of the active devices and estimates the corresponding
device-RIS channels. Finally, the BS conducts downlink data transmission based
on the estimated CSI. Since the BS and the RIS are placed at fixed positions,
the BS-RIS channel $\mathbf{U}$ is quasi-static [29]. As such, we only need to
estimate $\mathbf{U}$ once over a long period of time. The quasi-static RIS-BS
channel can be efficiently estimated by the dual-link pilot transmission
scheme proposed in [30], where the BS transmits pilots to the RIS via the
downlink channel with a single-antenna. Then, the RIS reflects the transmitted
pilots back to the BS via the uplink channel with a set of predesigned
reflection coefficients. Based on the received pilots, the BS can accurately
obtain the estimate of the quasi-static RIS-BS channel, i.e.,
$\hat{\mathbf{U}}$.

Figure 2: The proposed joint active device separation and channel estimation
frame structure for RIS-aided URA. Figure 3: Bit mapping of transmitted
symbols for the $k$th device.

On the other hand, a block transmission framework is adopted for joint active
device separation and device-RIS channel estimation. Specifically, each active
device divides its $B$-bit message into $L$ sub-messages [8] and the length of
the $l$-th sub-message is $B_{l}$ bits, as illustrated in Fig. 3. Then, active
devices transmit their $L$ sub-messages sequently, which are jointly recovered
at the BS. Noting that since each subblock $l$ contains $K_{a}$ recovered sub-
messages, only an instance of the transmitted data can be obtained in each
subblock. However, the ultimate goal of the BS is to recover the whole set of
$B$-bit messages that were transmitted by all the active devices. To achieve
this goal, a decoder for decoding binary messages among different subblocks of
all the $K_{a}$ active devices is needed. For example, in [8], an outer tree-
based decoder was adopted to stitch the decoded binary sequences across
different subblocks (cf. Fig. 4). Based on the recovered message, the BS can
separate the active devices if they embed their identity (ID) information into
the transmitted messages. Note that the inner decoder of URA is quite
different from the conventional blind detection problem [37]. First, the
devices in URA share a common codebook, while the devices in blind detection
may be assigned separate codebooks. Moreover, the receiver of URA neither
knows the transmitted codewords, nor the number of transmitted codewords. In
contrast, the blind detection generally has the knowledge of the number of
transmitted codewords [37].

Following the URA paradigm [8], we assume that all the devices in arbitrary
subblock $l$ exploit the same constellation
$\wp=\\{\mathbf{c}_{1},\mathbf{c}_{2},\cdots,\mathbf{c}_{2^{B_{l}}}\\}$
containing $2^{B_{l}}$ elements. Considering the angular space channel in (3)
and the RIS to the BS reflecting channel $\hat{\mathbf{U}}$, the received
signals of the $l$th subblock at the BS can be expressed as

$\displaystyle\mathbf{Y}_{l}=\sum_{k=1}^{K_{a}}{\hat{\mathbf{U}}}\text{diag}(\mathbf{A}_{\mathrm{R}}\bm{\lambda}_{k})\mathbf{v}_{l}\mathbf{s}_{k,l}^{H}+\mathbf{W}_{l},$
(4)

where $\mathbf{s}_{k,l}\in\wp\subset\mathbb{C}^{\tau}$ is the transmitted
sequence of complex valued baseband symbols from the $k$th device over $\tau$
channel uses. $\mathbf{W}_{l}\in\mathbb{C}^{M\times\tau}$ is the additive
white Gaussian noise. Herein, the coefficient
$[\mathbf{v}_{l}]_{n}=e^{\jmath\theta_{n}},n\in\\{1,\cdots,N\\}$ denotes the
phase shift of RIS reflecting element $n$ [31]. $\theta_{n}\in[0,2\pi]$ and
phase shift vector $\mathbf{v}_{l}$ is constant during the $l$-th subblock and
varies from one subblock to another subblock. Let
$\mathbf{y}_{l}=\mathrm{vec}(\mathbf{Y}_{l})$ and
$\mathbf{w}_{l}=\mathrm{vec}(\mathbf{W}_{l})$ denote the vectorized versions
of $\mathbf{Y}_{l}\in\mathbb{C}^{M\times\tau}$ and $\mathbf{W}_{l}$,
respectively. Moreover, we rearrange
$\hat{\mathbf{U}}\text{diag}(\mathbf{A}_{\mathrm{R}}\bm{\lambda}_{k})\mathbf{v}_{l}$
as
$\hat{\mathbf{U}}\text{diag}(\mathbf{A}_{\mathrm{R}}\bm{\lambda}_{k})\mathbf{v}_{l}\\!\\!=\\!\\!\sum_{n=1}^{N}\\!\hat{\mathbf{U}}(:,n)\\{\mathbf{A}_{\mathrm{R}}\bm{\lambda}_{k}\\}(n)\mathbf{v}_{l}(n)\\!\\!=\\!\\!\hat{\mathbf{V}}_{l}\mathbf{A}_{\mathrm{R}}\bm{\lambda}_{k},$
where
$\hat{\mathbf{V}}_{l}=[\mathbf{v}_{l}(1)\hat{\mathbf{U}}(:,1),\mathbf{v}_{l}(2)\hat{\mathbf{U}}(:,2),\cdots,\mathbf{v}_{l}(N)\hat{\mathbf{U}}(:,N)]\in\mathbb{C}^{M\times
N}$. Hence, (4) can be equivalently rewritten as

$\displaystyle\mathbf{y}_{l}=\sum_{k=1}^{K_{a}}\mathbf{s}_{k,l}\otimes\mathbf{P}_{l}\bm{\lambda}_{k}+\mathbf{w}_{l},$
(5)

with
$\mathbf{P}_{l}=\hat{\mathbf{V}}_{l}\mathbf{A}_{\mathrm{R}}\in\mathbb{C}^{M\times
N}$. Based on the received signal $\mathbf{y}_{l},l\in\\{1,...,L\\}$, the BS
jointly estimates the channel coefficients $\bm{\lambda}_{k}$ and recovers the
messages $\mathbf{s}_{k,l}$ for active device separation.

##  III Design of The RIS-Aided URA Architecture

In this section, we first design the RIS passive reflection for facilitating
RIS-assisted active device separation and channel estimation. Then, we
formulate the joint active device separation and channel estimation problem as
a high-order coupled tensor decomposition problem.

###  III-A RIS Reflection Design

The phase-shift vector $\mathbf{v}_{l},l=1,\cdots,L$, in (5) controls the
reflection of the waves impinging on the RIS. Given an optimized phase-shift
design, an RIS can significantly improve the performance of the URA system. In
the proposed RIS-aided URA system, the active device separation and the
channel estimation are performed concurrently, thus we propose a specific
optimized/random combined phase-shift design mechanism for $L$ subblocks. At
the $l=1,\cdots,L-1$ stage, the phase shift is randomly generated, which can
be regarded as a training process. At the last subblock, i.e., $l=L$, we
design the optimal phase shift to enhance device separation and channel
estimation.

Since the locations of the devices are a priori unknown, the ideal phase-shift
design should provide a low misdetection probability for the entire coverage
area. The authors in [32] has verified that the probability of misdetection is
a monotonically decreasing function with respect to the power of cascaded
channel.

Under the multipath scenario [33], the resultant IRS reflection optimization
problem is non-convex, and obtaining the global optimal solution is in general
intractable. As such, we adopt a compromised approach. Considering practical
delay-power spectral density in sparse channels, the first arrival path
usually contains the most power of the channel due to short propagation
distance [26]. In other words, the channel gain is mainly determined by the
first path. Motivated by such an observation, we formulate the worse-case
optimization problem for maximizing the minimum channel gain over the first
path among all devices, which is given by [32]

$\displaystyle\mathop{\max}\limits_{\\{\mathbf{v}_{L}\\}}\mathop{\min}\limits_{k\in[1,K]}\frac{\varrho^{2}}{16\pi^{2}d_{\mathrm{DR},k}^{2}d_{\mathrm{RB}}^{2}}\left|\mathbf{a}^{H}_{\mathrm{R}}(\phi_{1}^{k},\psi_{1}^{k})\text{diag}(\mathbf{a}_{\mathrm{R}}(\bar{\phi}_{1},\bar{\psi}_{1}))\mathbf{v}_{L}\right|^{2}$
$\displaystyle\mathrm{s.t.}~{}~{}~{}[\mathbf{v}_{L}]_{n}=1,n=1,\cdots,N,$ (6)

where $d_{\mathrm{DR},k}$ denotes the distance from device $k$ to RIS and
$d_{\mathrm{RB}}$ denotes the distance from the RIS to BS.

By setting
$\bar{\mathbf{a}}^{H}(\phi_{1}^{k},\psi_{1}^{k})=\mathbf{a}^{H}_{\mathrm{R}}(\phi_{1}^{k},\psi_{1}^{k})\text{diag}(\mathbf{a}_{\mathrm{R}}(\bar{\phi}_{1},\bar{\psi}_{1}))$
and
$\bar{\mathbf{A}}(\phi_{1}^{k},\psi_{1}^{k})=\bar{\mathbf{a}}(\phi_{1}^{k},\psi_{1}^{k})\bar{\mathbf{a}}^{H}(\phi_{1}^{k},\psi_{1}^{k})$,
problem (III-A) can be reformulated as

$\displaystyle\mathop{\max}\limits_{\\{\mathbf{V}_{L}\\}}\mathop{\min}\limits_{k\in[1,K]}\left\\{\frac{\varrho^{2}}{16\pi^{2}d_{\mathrm{DR},k}^{2}d_{\mathrm{RB}}^{2}}\mathrm{Tr}\left(\bar{\mathbf{A}}(\phi_{1}^{k},\psi_{1}^{k})\mathbf{V}_{L}\right)\right\\}$
$\displaystyle\mathrm{s.t.}~{}~{}~{}\mathbf{V}_{L}(n,n)=1,n=1,\cdots,N,$ (7)

where $\mathbf{V}_{L}=\mathbf{v}_{L}\mathbf{v}_{L}^{H}$. Note that the rank-
one constraint on $\mathbf{V}_{L}$ is dropped in problem (III-A). Fortunately,
problem (III-A) is a standard convex semidefinite program, which can be solved
by some optimization tools, i.e., CVX [34]. Yet,
$\mathrm{rank}(\mathbf{V}_{L}^{\mathrm{opt}})$ cannot be always guaranteed,
where $\mathbf{V}_{L}^{\mathrm{opt}}$ is the optimal solution to problem
(III-A). In this context, it is difficult to obtain the optimal phase-shift
vector $\mathbf{v}_{L}^{\mathrm{opt}}$ from $\mathbf{V}_{L}^{\mathrm{opt}}$.
To solve this problem, the following Gaussian randomization approach is
applied to $\mathbf{V}_{L}^{\mathrm{opt}}$ [35]:

  1. 1.

For $j\in\\{0,1,\cdots,J-1\\}=\mathcal{J}$, generate $J$ random vectors
$\bm{\nu}_{j}\sim\mathcal{CN}_{N}(\bm{\nu}_{j}|\mathbf{0},\mathbf{V}_{L}^{\mathrm{opt}})$
and set $|[\bm{\nu}_{j}]_{n}|=1,\forall n,\forall j$.

  2. 2.

Search
$\hat{\mathbf{v}}_{L}^{\mathrm{opt}}=\arg\mathop{\max}\limits_{\bm{\nu}_{j},\forall
j\in\mathcal{J}}\mathop{\min}\limits_{k\in\mathcal{Q}}\left\\{\bm{\nu}_{j}^{H}\bar{\mathbf{A}}(\phi_{1}^{k},\psi_{1}^{k})\bm{\nu}_{j}\right\\}$.

Thus, we can obtain a sub-optimal phase-shift vector
$\hat{\mathbf{v}}_{L}^{\mathrm{opt}}$.

###  III-B Coupled Tensor Modeling

Based on the designed phase shift $\hat{\mathbf{v}}_{l},l=1,\cdots,L$, the
data symbol is transmitted from the active devices. In this paper, we
construct the constellation structure $\wp$ according to the tensor
decomposition format [36]. As shown in Fig. 3, it is assumed that the channel
use, i.e., $\tau$, can be factorized as $\tau=\Pi_{i}^{d}\tau_{i}$ for some
$d\geq 2$ and $\tau_{i}\geq 2,\forall i$. Subsequently, the complex symbols
transmitted by the $k$th device at the $l$th subblock, denoted by
$\mathbf{s}_{k,l}\in\mathbb{C}^{\tau}$, can be rewritten as the vector
representation of a rank-$1$ tensor
$\mathcal{S}_{k,l}\in\mathbb{C}^{\tau_{1}\times\tau_{2}\times\cdots\times\tau_{d}}$
of dimensions $\tau_{1},\tau_{2},\cdots,\tau_{d}$, that is

$\displaystyle\mathbf{s}_{k,l}=\mathrm{vec}(\mathcal{S}_{k,l})\in\mathbb{C}^{\Pi_{i}^{d}\tau_{i}}=\mathbb{C}^{\tau},$
(8)

with
$\mathcal{S}_{k,l}=\mathbf{x}_{1,k,l}\circ\mathbf{x}_{2,k,l}\circ\cdots\circ\mathbf{x}_{d,k,l},\forall
k,l$, where each $\mathbf{x}_{i,k,l}$ is generated from a sub-constellation
$\wp_{i}$, which is defined as a discrete subset of $\mathbb{C}^{\tau_{i}}$.
As shown in Fig. 3, bit information can be mapped to symbol as follows.
$B_{l}$ coded bits are split into $d$ sets of
$\\{B_{1,1}\\},...,\\{B_{l,d}\\}$ bits, respectively, corresponding to the $d$
tensor dimensions. For the $i$th set, $B_{l,i}$-bit data is mapped to an
element $\mathbf{x}_{i,k,l}$. Then, the vector symbol $\mathbf{s}_{k,l}$ is
formed by the calculation in (8). Actually, all possible combinations of
elements in the sub-constellations result in the vector constellation
$\wp\subset\mathbb{C}^{\tau}$, which is given by
$\wp=\left\\{\mathrm{vec}(\mathbf{x}_{1}\circ\mathbf{x}_{2}\circ\cdots\circ\mathbf{x}_{d}),\mathbf{x}_{1}\\!\in\\!\wp_{1},\\!\cdots,\\!\mathbf{x}_{d}\in\wp_{d}\right\\}$.

Figure 4: High-level description of the proposed RIS-based unsourced random
access scheme, where MSG $k$ denotes the message of the $k$th device.

Substituing (8) into (5), we can reformulate (5) as the following equivalent
tensor decomposition form

$\displaystyle\\!\\!\\!\\!\\!\mathcal{Y}_{l}\\!=\\!\underbrace{\sum_{k=1}^{K_{a}}\mathbf{x}_{1,k,l}\circ\mathbf{x}_{2,k,l}\circ\mathbf{x}_{3,k,l}\circ\cdots\circ\mathbf{x}_{d,k,l}\circ\mathbf{P}_{l}\bm{\lambda}_{k}}_{\mathcal{Y}_{l}^{0}}+\mathcal{W}_{l},$
(9)

where $\mathcal{W}_{l}\in\mathbb{C}^{\tau_{1}\times\cdots\tau_{d}\times M}$ is
the additive white Gaussian noise represented in the tensor space. Therefore,
in (9), each active device transmits a sequence that is associated with a
rank-one tensor of order $d$, while $\mathcal{Y}_{l}^{0}$ is the noise-free
tensor of order $d+1$ and has dimensions $\tau_{1},\cdots,\tau_{d},M$,
respectively. Its tensor rank is at most $K_{a}$, since a tensor is said to be
rank-$K_{a}$ whenever $K_{a}$ is the smallest integer such that the tensor can
be written as a sum of $K_{a}$ rank-1 tensors [21].

Consequently, we can formulate the design of joint active device separation
and channel estimation as the following non-convex optimization problem:

$\displaystyle\mathop{\arg\min}\limits_{\mathbf{x}_{i,k,l}\in\wp_{i}^{\tau_{i}},\bm{\lambda}_{k}\in\mathbb{C}^{N}}\sum_{l=1}^{L}\left\|\mathcal{Y}_{l}-\sum_{k=1}^{K_{a}}\mathbf{x}_{1,k,l}\circ\cdots\circ\mathbf{x}_{d,k,l}\circ(\mathbf{P}_{l}\bm{\lambda}_{k})\right\|_{F}^{2}$
$\displaystyle\mathrm{s.t.}~{}\|\bm{\lambda}_{k}\|_{0}\leq\zeta_{s},~{}~{}k=1,2,\cdots,K_{a},$
(10)

where $\zeta_{s}$ is a predefined parameter for imposing the channel sparsity.
By solving problem (III-B), we can obtain
$\mathbf{x}_{1,k,1},\cdots,\mathbf{x}_{d,k,L}$ and $\bm{\lambda}_{k}$ for
$k=1,2,\cdots,K$. After that, the estimation of transmitted data symbol
$\hat{\mathbf{s}}_{k,l}$ can be formed by computing the vector outer product
according to (8). The high-level description of the proposed RIS-based URA
scheme is shown in Fig. 4.

###  III-C The Challenges Ahead

For the proposed RIS-aided URA architecture, the key is to solve problem
(III-B). However, there are two major challenges that prohibit the
straightforward application of some existing tensor decomposition advances
(e.g., using either optimization methods in [19] and [20] or applying the
Bayesian framework in [22] and [24]):

  1. 1.

Problem (III-B) is a non-standard coupled tensor decomposition problem as the
sparse profile vector is common and nonlinearly coupled with other parameters
in modeling multiple tensors $\mathcal{Y}_{l},\forall l$. In contrast, in [19]
and [20], the coupled factor matrix $\mathbf{B}_{k}$ is common for $\forall
l\in\\{1,\cdots,L\\}$ and does not pose a sparse structure. However, in
problem (III-B), the coupled information is a sparse vector $\bm{\lambda}_{k}$
come from $\mathbf{B}_{k}=\mathbf{P}_{l}\bm{\lambda}_{k}$ with known
$\mathbf{P}_{l}$. The exploitation of such a sparse structure can be exploited
to improve channel estimation performance and to reduce the algorithm
complexity.

  2. 2.

In [19] and [20], the authors assumed that the tensor rank $K_{a}$ is known.
However, the tensor rank (number of active devices $K_{a}$) in problem (III-B)
is random and unknown in practice introducing NP-hardness to the estimation
based on tensor data. The accurate information about the model complexity is
essential in avoiding overfitting of noises or underfitting of signals.
Therefore, we introduce the variational inference algorithm that can
automatically learn both the transmit data, the channel state information, and
the active device number from the measured data at the BS under the Bayesian
learning framework.

##  IV Coupled Tensor-Based Automatic Detection Algorithm: Bayesian Learning

Note that solving the discrete problem in (III-B) optimally via an exhaustive
search, however, requires $2^{K_{a}B}$ number of evaluations of the objective
function. To circumvent the complexity issue, we first relax the discrete
domain
${{\\{\mathbf{x}_{i,k,l}\in\wp_{i}^{\tau_{i}}\\}_{i=1}^{d}}_{k=1}^{K_{a}}}_{l=1}^{L}$
in (III-B) to a continuous one
${{\\{\mathbf{x}_{i,k,l}\in\mathbb{C}^{\tau_{i}}\\}_{i=1}^{d}}_{k=1}^{K_{a}}}_{l=1}^{L}$.
Besides, $K_{a}$ is unknown for the modeling of all $\mathcal{Y}_{l}$,
$l=1,2,\cdots,L$, and its estimation has been shown to be NP-hard. To tackle
this challenge, an effective scheme is to introduce two regularization terms
[23] that penalize the model complexity and avoid possible overfitting of
noise as follows

$\displaystyle\mathop{\arg\min}\limits_{{\\{\mathbf{X}_{l}^{i}\in\mathbb{C}^{\tau_{i}\times
K}\\}_{i=1}^{d}}_{l=1}^{L}\mathbf{G}\in\mathbb{C}^{M\times
K}}\sum_{l=1}^{L}\left\|\mathcal{Y}_{l}-\left[\\!\\!\left[\mathbf{X}_{l}^{1},\cdots,\mathbf{X}_{l}^{d},~{}\mathbf{P}_{l}\mathbf{G}\right]\\!\\!\right]\right\|_{F}^{2}$
$\displaystyle\\!+\\!\mathop{\sum}\limits_{k=1}^{K}\gamma_{k}\\!\left(\\!\sum_{l=1}^{L}\mathop{\sum}\limits_{i=1}^{d}\mathbf{X}_{l}^{i}(:,k)^{H}\mathbf{X}_{l}^{i}(:,k)\\!\right)\\!+\\!\mathop{\sum}\limits_{k=1}^{K}\eta_{k}\mathbf{G}(:,k)^{H}\mathbf{G}(:,k)$
$\displaystyle\mathrm{s.t.}~{}\|\mathbf{G}(:,k)\|_{0}\leq\zeta_{s},~{}~{}k=1,2,\cdots,K,$
(11)

where $\mathbf{X}_{l}^{i}\in\mathbb{C}^{\tau_{i}\times K}$ with the $k$th
column being $\mathbf{x}_{d,k,l}$, $\mathbf{G}\in\mathbb{C}^{M\times K}$ is
defined similarly but with the $k$th column being $\bm{\lambda}_{k}$. Herein,
we set the column number of all factor matrices
$\\{\mathbf{X}_{l}^{1},\cdots,\mathbf{X}_{l}^{d},\mathbf{G}\\}$ as $K$, which
is the maximum possible number of active devices $K_{a}$. The idea about
regularization term is that if $\gamma_{k}>0$ and $\eta_{k}>0$ are
sufficiently large (e.g., $10^{3}$) after inference, the elements of the $k$th
columns in the optimal ${\\{\mathbf{X}_{l}^{i}\\}_{l=1}^{L}}_{i=1}^{d}$ and
$\mathbf{G}$ approach zeros. Then, the corresponding column can be pruned out
and the number of remaining columns with significant non-zero values in each
factor matrix returns an estimate of number of active devices. However, the
choice of regularization parameters $\gamma_{k}$ and $\eta_{k}$ is not
straightforward and computationally demanding, since setting $\gamma_{k}$ and
$\eta_{k}$ too large would lead to excessive residual squared error, while
setting them too small may incur overfitting of noises. Therefore, we develop
an intelligent algorithm that can automatically learn both the factor matrices
and regularization parameters from the measured data at the BS under the
Bayesian learning framework.

###  IV-A Element-Wise Sparsity and Low-Rank Inducing Probabilistic Modeling

Firstly, to apply the Bayesian learning framework, a probabilistic model,
which encodes the knowledge of problem (IV), needs to be established. To
strike the trade-off between the flexibility of knowledge encoding and the
tractability of the inference algorithm, we propose a novel probabilistic
model by interpreting each term in problem (IV) via some probability density
functions (pdfs) [23]. We start with the last regularization term in the
objective function of (IV), which can be modeled as a circularly-symmetric
complex Gaussian prior distribution111This is inspired by the $l_{2}$ norm in
the last regularization term and the assumption that device to RIS channel is
independent among the different devices. of the columns in matrix
$\mathbf{G}$, i.e.,
$\prod_{k=1}^{K}\mathcal{CN}\left(\mathbf{G}(:,k)|\mathbf{0},(\eta_{k})^{-1}\mathbf{I}\right)$.
On the other hand, note that the $l_{0}$-norm constraint in (IV) is imposed
such that the elements in each column of channel are also sparse. This
observation inspires a sparsity-aware probabilistic modeling on each element
of $\mathbf{G}$. Therefore, by taking both the column-wise sparsity (i.e.,
low-rank) and the element-wise sparsity structure into account, we propose the
following novel prior for matrix $\mathbf{G}$:

$\displaystyle
p\left(\mathbf{G}|\\{\eta_{k}\\}_{k=1}^{K},{\\{\bm{\xi}(n,k)\\}_{n=1}^{N}}_{k=1}^{K}\right)$
$\displaystyle=\prod_{k=1}^{K}\\!\left[\mathcal{CN}\left(\mathbf{G}(:,k)|\mathbf{0},\eta_{k}^{-1}\mathbf{I}\right)\right.\left.\prod_{n=1}^{N}\mathcal{CN}\left(\\!\mathbf{G}(n,k)|0,\bm{\xi}(n,k)^{-1}\right)\right],$
(12)

with

$\displaystyle p(\\{\eta_{k}\\}_{k=1}^{K}|\bm{\imath}_{\eta})$
$\displaystyle=\prod_{k=1}^{K}\\!\\!\mathrm{gamma}(\eta_{k}|\delta,\delta)\\!\\!=\\!\\!\prod_{k=1}^{K}\\!\\!\eta_{k}^{\delta-1}\exp(-\delta\eta_{k}),$
(13) $\displaystyle
p({\\{\bm{\xi}(n,k)\\}_{n=1}^{N}}_{k=1}^{K}|\bm{\imath}_{\bm{\xi}})=\prod_{n=1}^{N}\prod_{k=1}^{K}\mathrm{gamma}(\bm{\xi}(n,k)|\delta,\delta)$
$\displaystyle=\prod_{n=1}^{N}\prod_{k=1}^{K}\bm{\xi}(n,k)^{\delta-1}\exp(-\delta\bm{\xi}(n,k)),$
(14)

where $\delta>0$ is a small number that indicates the non-informativeness of
the prior model, the natural parameter
$\bm{\imath}_{\eta}=[-\delta\mathbf{1}_{K};(\delta-1)\mathbf{1}_{K}]$, and
$\bm{\imath}_{\bm{\xi}}=[-\delta\mathbf{1}_{NK};(\delta-1)\mathbf{1}_{NK}]$.

The proposed prior in (IV-A) is employed to model the sparsity pattern of
$\mathbf{G}$. Different from the previous works [22, 23] that only adopted
Gaussian-gamma pairs for modeling the column-wise sparsity, the prior in
(IV-A) utilizes the product rule of probability to further encode the element-
wise sparsity information. The proposed prior has a clear physical
interpretation. Particularly, $\bm{\xi}(n,k)^{-1}$ can be interpreted as the
power of each element $\mathbf{G}(n,k)$, while $\eta_{k}^{-1}$ has a physical
interpretation as the power of each column in $\mathbf{G}$. Therefore, if
$\eta_{k}^{-1}$ is learnt to approach zero, regardless of $\bm{\xi}(n,k)$, the
corresponding columns in the factor matrices play no role in channel modeling
and thus can be pruned out by thresholding. Similarly, for a nonzero column
($\eta_{k}^{-1}\neq 0$), $\mathbf{G}(n,k)$ would be shrunk to zero as
$\bm{\xi}(n,k)^{-1}$ goes to zero. Therefore, the proposed prior in (IV-A) can
simultaneously promote element-wise sparsity of the channel and low-rank
property of tensor data, which mimic the sparsity constraint and the last
regularization term in (IV). Note that the conventional way to promote both
element-wise and column-wise sparsity is to use hierarchical Bayesian [39].
However, the resulting penalty is less capable in enforcing a structural low-
rank to the tensor. In this context, multiplying two priors serves as a viable
approach since both hyperparameters influence $\mathbf{G}$ directly.

Similarly, for the regularization term about $\mathbf{X}_{l}^{i}$ in problem
(IV), it can be modeled as a zero-mean circularly-symmetric complex Gaussian
prior distribution over the columns of the factor matrices as follows [22].

$\displaystyle
p\left({\\{\mathbf{X}_{l}^{i}\\}_{l=1}^{L}}_{i=1}^{d}|\\{\gamma_{k}\\}_{k=1}^{K}\right)=\prod_{l=1}^{L}\prod_{i=1}^{d}\prod_{k=1}^{K}\mathcal{CN}\left(\mathbf{X}_{l}^{i}(:,k)|\mathbf{0},\gamma_{k}^{-1}\mathbf{I}\right).$
(15)

A gamma distribution is adopted for the penalty parameter $\gamma_{k}$ to
enforce the column-wise sparsity [22]:

$\displaystyle\\!\\!\\!\\!\\!p(\\{\gamma_{k}\\}_{k=1}^{K}|\bm{\imath}_{\gamma})\\!=\\!\prod_{k=1}^{K}\\!\\!\mathrm{gamma}(\gamma_{k}|\delta,\delta)\\!\\!=\\!\\!\prod_{k=1}^{K}\\!\\!\gamma_{k}^{\delta-1}\exp(-\delta\gamma_{k}),$
(16)

where $\gamma_{k}^{-1}$ can be interpreted as the power of each column in
$\mathbf{X}_{l}^{i}$ and
$\bm{\imath}_{\gamma}=[-\delta\mathbf{1}_{K};(\delta-1)\mathbf{1}_{K}]$.
Noting that after integrating the gamma hyper-prior, the marginal distribution
of model parameters (15) is a Student’s $t$ distribution, which is strongly
peaked at zero and with heavy tails, thus promoting sparsity. Furthermore, the
gamma hyper-prior (16) is conjugate to the Gaussian prior (15). This conjugacy
permits the closed-form solution of the variational inference [39].

Finally, since the elements of the additive noise $\mathcal{W}_{l}$ obeys
white Gaussian distribution, the sum of the squared error term in problem (IV)
can be interpreted as the negative log of a likelihood function given by

$\displaystyle
p(\\{\mathcal{Y}_{l}\\}_{l=1}^{L}|{\\{\mathbf{X}_{l}^{i}\\}_{i=1}^{d}}_{l=1}^{L},\mathbf{G},\beta)$
$\displaystyle\propto\exp\left(-\beta\mathop{\sum}\limits_{l=1}^{L}\left\|{\mathcal{Y}_{l}}\right.\right.\left.\left.-\left[\\!\\!\left[\mathbf{X}_{l}^{1},~{}\mathbf{X}_{l}^{2},\cdots,\mathbf{X}_{l}^{d},~{}\mathbf{P}_{l}\mathbf{G}\right]\\!\\!\right]\right\|_{F}^{2}\right),$
(17)

where the noise precision $\beta$ obeys gamma distribution, i.e.,
$p(\beta|\bm{\imath}_{\beta})\propto\beta^{\delta-1}\exp(-\delta\beta)$ with
natural parameters $\bm{\imath}_{\beta}=[-\delta;(\delta-1)]$.

Figure 5: The proposed probabilistic model for problem (IV) that incorporates
the information of both element-wise sparsity and low-rank property. Filled-in
circle and cube indicate the observable variables, hollow circles and cube
denote the unknown variables, the boxes denote pre-specified hyperparameters,
and arrows describe conditional dependencies between variables.

Graphical illustration of the proposed complete probabilistic model is
presented in Fig. 5. Let
$\bm{\Theta}=\\{{{\\{\mathbf{X}_{l}^{i}\\}_{i=1}^{d}}_{l=1}^{L}},\mathbf{G},\beta,{\\{\gamma_{k}\\}_{k=1}^{K}},{\\{\eta_{k}\\}_{k=1}^{K}},{\\{\bm{\xi}(n,k)\\}_{n=1}^{N}}_{k=1}^{K}\\}$
collect all the unknown random variables, the joint pdf of $\bm{\Theta}$ and
$\\{\mathcal{Y}_{l}\\}_{l=1}^{L}$ is given as follows:

$\displaystyle
p(\bm{\Theta},\\{\mathcal{Y}_{l}\\}_{l=1}^{L})=p(\\{\mathcal{Y}_{l}\\}_{l=1}^{L}|{{\\{\mathbf{X}_{l}^{i}\\}_{i=1}^{d}}_{l=1}^{L}},\mathbf{G},\beta)p({{\\{\mathbf{X}_{l}^{i}\\}_{i=1}^{d}}_{l=1}^{L}}|{\\{\gamma_{k}\\}_{k=1}^{K}})p\\!\left(\\!\mathbf{G}|\\{\eta_{k}\\}_{k=1}^{K},{\\{\bm{\xi}(n,k)\\}_{n=1}^{N}}_{k=1}^{K}\\!\right)\\!$
$\displaystyle\times
p\\!\left(\\!\beta|\bm{\imath}_{\beta}\\!\right)\\!p\\!\left(\\!{\\{\gamma_{k}\\}_{k=1}^{K}}|\bm{\imath}_{\gamma}\\!\right)\\!p\\!\left(\\!{\\{\eta_{k}\\}_{k=1}^{K}}|\bm{\imath}_{\eta}\\!\right)\\!p\\!\left(\\!{\\{\bm{\xi}(n,k)\\}_{n=1}^{N}}_{k=1}^{K}|\bm{\imath}_{\bm{\xi}}\\!\right)\\!$
$\displaystyle\propto\exp\left\\{-\beta\sum_{l=1}^{L}\left\|\mathcal{Y}_{l}-\left[\\!\\!\left[\mathbf{X}_{l}^{1},~{}\mathbf{X}_{l}^{2},\cdots,\mathbf{X}_{l}^{d},~{}\mathbf{P}_{l}\mathbf{G}\right]\\!\\!\right]\right\|_{F}^{2}\right.\left.\\!-\\!\sum_{l=1}^{L}\sum_{i=1}^{d}\mathrm{Tr}(\bm{\Lambda}\mathbf{X}_{l}^{i,H}\mathbf{X}_{l}^{i})\\!+\\!L\left(\Pi_{i=1}^{d+1}\tau_{i}M\right)\ln\beta\right.$
$\displaystyle\left.\\!-\\!\delta\beta\\!-\\!\sum_{k=1}^{K}\delta\gamma_{k}+(\delta\\!-\\!1)\ln\beta-\mathrm{Tr}(\bar{\bm{\Lambda}}\mathbf{G}^{H}\mathbf{G})\\!-\\!\sum_{k=1}^{K}\mathrm{Tr}(\mathbf{Z}_{k}\mathbf{G}(:,k)\mathbf{G}(:,k)^{H})\\!+\\!(N\\!+\\!\delta\\!-\\!1)\sum_{k=1}^{K}\ln{\eta_{k}}\right.$
$\displaystyle\left.\\!-\\!\sum_{k=1}^{K}\delta\eta_{k}\\!-\\!\delta\sum_{n=1}^{N}\sum_{k=1}^{K}{\bm{\xi}(n,k)}\\!+\\!(N\\!+\\!\delta\right.\left.\left.\\!-\\!1)\sum_{n=1}^{N}\sum_{k=1}^{K}\ln{\bm{\xi}(n,k)}\\!+\\!\left(\\!L\sum_{i=1}^{d}\tau_{i}+\delta\\!-\\!1\\!\right)\sum_{k=1}^{K}\ln{\gamma_{k}}\right.\\!\\!\right\\}\\!,\\!\\!\\!$
(18)

with
$\mathbf{Z}_{k}=\text{diag}(\bm{\xi}({1,k}),\bm{\xi}({2,l}),\cdots,\bm{\xi}(N,k)\in\mathbb{C}^{N\times
N}$,
$\bm{\Lambda}\\!=\\!\text{diag}(\gamma_{1},\\!\cdots,\\!\gamma_{K}\\!)\in\mathbb{C}^{K\times
K}$, and
$\bar{\bm{\Lambda}}\\!=\\!\text{diag}(\eta_{1},\\!\cdots,\\!\eta_{K}\\!)\in\mathbb{C}^{K\times
K}.$ Given the probabilistic model
$p(\bm{\Theta},\\{\mathcal{Y}_{l}\\}_{l=1}^{L})$ in (IV-A) at the top of next
page, the next goal of Bayesian inference is to learn the model parameters
$\bm{\Theta}$ from the tensor data $\\{\mathcal{Y}_{l}\\}_{l=1}^{L}$, where
the posterior distribution of unknown model parameters, i.e.,
$p(\bm{\Theta}|\\{\mathcal{Y}_{l}\\}_{l=1}^{L})$, is needed to be sought
before further proceeding. Noting that maximizing the posterior probability
$p(\bm{\Theta}|\\{\mathcal{Y}_{l}\\}_{l=1}^{L})=p\left(\bm{\Theta},\\{\mathcal{Y}_{l}\\}_{l=1}^{L}\right)/p\left(\\{\mathcal{Y}_{l}\\}_{l=1}^{L}\right)$
is similar to addressing the problem (IV). The difference is that the problem
(IV) cannot learn the regularization parameters.

###  IV-B Algorithm Development via Variational Inference

Unfortunately, the probabilistic model described by the joint pdf (IV-A) is
still sophisticated. In particular, solving the multiple integrals for
computing the posterior distribution is generally intractable. To address this
problem, we adopt the variational inference method which establishes a
variational distribution $Q(\bm{\Theta})$ to approximate the true posterior
$p(\bm{\Theta}|\\{\mathcal{Y}_{l}\\}_{l=1}^{L})$. To achieve this goal,
$Q(\bm{\Theta})$ is the solution which minimizes the Kullback-Leibler (KL)
divergence, i.e.,

$\displaystyle\mathop{\mathrm{missing}}{minimize}\limits_{Q(\bm{\Theta})}~{}\mathrm{KL}(Q(\bm{\Theta})|p(\bm{\Theta}|\\{\mathcal{Y}_{l}\\}_{l=1}^{L}))\triangleq\mathop{\mathrm{missing}}{minimize}\limits_{Q(\bm{\Theta})}~{}$
$\displaystyle-\mathbb{E}_{Q(\bm{\Theta})}\left\\{\ln\frac{p(\bm{\Theta}|\\{\mathcal{Y}_{l}\\}_{l=1}^{L})}{Q(\bm{\Theta})}\right\\}.$
(19)

To address the problem (IV-B), mean-field approximation [39] is widely used to
offer a tractable solution of (IV-B), which assumes that the variational pdf
can be represented in a fully factorized form, i.e.,
$Q(\bm{\Theta})=\mathop{\prod}\limits_{j=1}^{J}Q({\Theta}_{j})$, where
${\Theta}_{j}\in\bm{\Theta}$ is part of $\bm{\Theta}$ with
$\bigcup_{j=1}^{J}\bm{\Theta}_{j}=\bm{\Theta}$ and
$\bigcap_{j=1}^{J}\bm{\Theta}_{j}=\varnothing$, and $J$ is the number of
subsets. With this factorization, (IV-B) reduces to

$\displaystyle\mathop{\mathrm{missing}}{minimize}\limits_{Q(\\{{\Theta}_{j}\\}_{j=1}^{J})}-\mathbb{E}_{\\{Q({\Theta}_{j})\\}_{j=1}^{J}}\left\\{\ln\left(\frac{p(\bm{\Theta}|\\{\mathcal{Y}_{l}\\}_{l=1}^{L})}{\prod_{j=1}^{J}Q({\Theta}_{j})}\right)\right\\}.$
(20)

The factorized structure of $\\{Q(\Theta_{j})\\}_{j=1}^{J}$ motivates the use
of block coordinate descent to obtain a suboptimal solution of (20).
Particularly, by fixing the variational pdfs $\\{Q(\Theta_{i})\\}_{i\neq j}$,
$Q(\Theta_{j})$ is optimized as follows:

$\displaystyle\mathop{\mathrm{missing}}{minimize}\limits_{Q({\Theta}_{j})}~{}$
$\displaystyle\int Q({\Theta}_{j})\left(-\mathbb{E}_{\prod_{i\neq
j}Q({\Theta}_{i})}\ln p(\bm{\Theta},\\{\mathcal{Y}_{l}\\}_{l=1}^{L})\right.$
$\displaystyle\left.+\ln Q({\Theta}_{j})\right)d{\Theta}_{j}.$ (21)

It has been shown in [39] that the optimal solution of (IV-B) is given by

$\displaystyle Q^{\dagger}({\Theta}_{j})=\frac{\exp(\mathbb{E}_{\prod_{i\neq
j}Q({\Theta}_{i})}\ln
p(\bm{\Theta},\\{\mathcal{Y}_{l}\\}_{l=1}^{L}))}{\int\mathbb{E}_{\prod_{i\neq
j}Q({\Theta}_{i})}\ln
p(\bm{\Theta},\\{\mathcal{Y}_{l}\\}_{l=1}^{L})d\Theta_{j}},\forall j.$ (22)

Using (22), we derive the closed-form posterior update for each variational
pdf, i.e., $Q^{\dagger}({\Theta}_{j})$ in the following subsections.

###  IV-C Learning Device-RIS Channel: Inferring Coupled Factor
$Q(\mathbf{G})$

We start by deriving the explicit expression of the optimal variational pdf
for the coupled factor $Q^{\dagger}(\mathbf{G})$, whose mean matrix gives the
MMSE estimate of the interested device-RIS channel matrix [38]. It is clear
that the likelihood function in (IV-A) results in complicated couplings among
the shared factor matrix $\mathbf{G}$, which make the derivation of
$Q^{\dagger}(\mathbf{G})$ challenging. To overcome this difficulty, we first
define
$\mathcal{Y}_{l}(d+1)\in\mathbb{C}^{M\times\tau_{1}\tau_{2}\cdots\tau_{d}}$ as
the unfolding operation [23] on a $(d+1)$th-order tensor
$\mathcal{Y}_{l}\in\mathbb{C}^{\tau_{1}\times\cdots\times\tau_{d}\times M}$
along its $(d+1)$-th mode. By substituting (IV-A) into (22), using the result
$\left\|\mathbf{A}\right\|_{F}^{2}=\mathrm{Tr}(\mathbf{A}\mathbf{A}^{H})$ and
then only keeping the terms relevant to $\mathbf{G}$, we obtain

$\displaystyle
Q^{\dagger}(\mathbf{G})\propto\exp\left\\{\mathbb{E}\left[-\beta\sum_{l=1}^{L}\left\|\mathcal{Y}_{l}-\left[\\!\\!\left[\mathbf{X}_{l}^{1},\cdots,\mathbf{X}_{l}^{d},~{}\mathbf{P}_{l}\mathbf{G}\right]\\!\\!\right]\right\|_{F}^{2}-\sum_{k=1}^{K}\mathrm{Tr}(\mathbf{Z}_{k}\mathbf{G}(:,k)\mathbf{G}(:,k)^{H})\right.\right.$
$\displaystyle\left.\left.-\mathrm{Tr}(\bar{\bm{\Lambda}}\mathbf{G}^{H}\mathbf{G})\right]\right\\}\propto\\!\exp\\!\left\\{\right.\left.\\!-\mathrm{Tr}(\sum_{l=1}^{L}\mathbf{P}_{l}^{H}\mathbf{P}_{l}\mathbf{G}\underbrace{\mathbb{E}[\beta]\mathbb{E}\left[\left(\\!\mathop{\diamond}\limits_{i=1}^{d}\mathbf{X}_{l}^{i}\\!\right)^{T}\\!\left(\\!\mathop{\diamond}\limits_{i=1}^{d}\mathbf{X}_{l}^{i}\\!\right)^{*}\right]}_{\bm{\Xi}_{l}^{-1}}\mathbf{G}^{H}\\!-\\!\mathbf{P}_{l}^{H}\mathbf{P}_{l}\right.$
$\displaystyle\left.\left.\times\underbrace{\mathbb{E}[\beta](\mathbf{P}_{l}^{H}\mathbf{P}_{l})^{-1}\mathbf{P}_{l}^{H}\mathcal{Y}_{l}(d+1)\left(\mathop{\diamond}\limits_{i=1}^{d}\mathbb{E}[\mathbf{X}_{l}^{i}]\right)^{*}\bm{\Xi}_{l}}_{\bm{\Phi}_{l}}\bm{\Xi}_{l}^{-1}\mathbf{G}^{H}-\mathbf{P}_{l}^{H}\mathbf{P}_{l}\mathbf{G}\bm{\Xi}_{l}^{-1}{\bm{\Phi}_{l}^{H}})\right\\}\right.$
$\displaystyle\times\exp\left\\{\\!-\\!\mathrm{vec}(\mathbf{G}^{H})^{H}(\mathbb{E}[\bm{\Gamma}]\\!+\\!\mathbb{E}[\mathbf{Z}])\mathrm{vec}(\mathbf{G}^{H})\\!\right\\},$
(23)

with
$\mathbf{Z}\\!\\!=\\!\\!\text{diag}(\bm{\xi}(1,1),\\!\cdots\\!,\bm{\xi}(1,K),\\!\cdots\\!,\bm{\xi}(N,1),\\!\cdots\\!,\bm{\xi}(N,K))\\!\in\\!\mathbb{C}^{NK\times
NK}$, and
$\bm{\Gamma}=\text{diag}(\mathbf{1}_{N}\otimes\left[\eta_{1},\cdots,\eta_{K}\right])\in\mathbb{C}^{NK\times
NK}$, where
$\mathop{\diamond}\limits_{i=1}^{d}\mathbf{X}_{l}^{i}=\mathbf{X}_{l}^{d}\mathop{\diamond}\limits\mathbf{X}_{l}^{d-1}\mathop{\diamond}\limits\cdots\mathop{\diamond}\limits\mathbf{X}_{l}^{1}$
denotes the multiple Khatri-Rao products. Let
$\mathbf{A}=(\mathbf{a}_{1},\cdots,\mathbf{a}_{n})\in\mathbb{C}^{K\times N}$
with $\mathbf{a}_{n}\in\mathbb{C}^{K}$. The vectorization operation on
$\mathbf{A}$ is expressed by
$\mathrm{vec}(\mathbf{A})=(\mathbf{a}_{1}^{T},\cdots,\mathbf{a}_{n}^{T})^{T}\in\mathbb{C}^{NK}$.

However, there are some complicated matrix multiplications in the last three
lines of formula (23) at the top of this page arising from
$\mathbf{P}_{l},\forall l$, which make it difficult to identify the
variational pdf $Q^{\dagger}(\mathbf{G})$. To reveal its structure, we utilize
the property of the vectorization operation that
$\mathrm{Tr}(\mathbf{A}\mathbf{X}\mathbf{B}\mathbf{X}^{T})=\mathrm{vec}(\mathbf{X})^{T}(\mathbf{B}\otimes\mathbf{A}^{T})\mathrm{vec}(\mathbf{X})$
with $\mathbf{X}\in\mathbb{C}^{m\times n}$, $\mathbf{A}\in\mathbb{C}^{m\times
m}$ and $\mathbf{B}\in\mathbb{C}^{n\times n}$ [39], and then the terms inside
the second $\exp(\cdot)$ in (23) can be reformulated as

$\displaystyle-\mathrm{Tr}\left(\sum_{l=1}^{L}\mathbf{P}_{l}^{H}\mathbf{P}_{l}\mathbf{G}\bm{\Xi}_{l}^{-1}\mathbf{G}^{H}\right.\left.-\mathbf{P}_{l}^{H}\mathbf{P}_{l}\bm{\Phi}_{l}\bm{\Xi}_{l}^{-1}\mathbf{G}^{H}\right.$
$\displaystyle\left.-\mathbf{P}_{l}^{H}\mathbf{P}_{l}\mathbf{G}\bm{\Xi}_{l}^{-1}\bm{\Phi}_{l}^{H}\right)$
$\displaystyle=-\mathrm{Tr}\left(\sum_{l=1}^{L}\bm{\Xi}_{l}^{-1}\mathbf{G}^{H}\mathbf{P}_{l}^{H}\mathbf{P}_{l}\mathbf{G}\right.$
$\displaystyle\left.-\bm{\Xi}_{l}^{-1}\bm{\Phi}_{l}^{H}\mathbf{P}_{l}^{H}\mathbf{P}_{l}\mathbf{G}-\bm{\Xi}_{l}^{-1}\mathbf{G}^{H}\mathbf{P}_{l}^{H}\mathbf{P}_{l}\bm{\Phi}_{l}\right)$
$\displaystyle=\sum_{l=1}^{L}-\mathrm{vec}(\mathbf{G}^{H})^{H}\left(\mathbf{P}_{l}^{H}\mathbf{P}_{l}\otimes\bm{\Xi}_{l}^{-1}\right)\mathrm{vec}(\mathbf{G}^{H})$
$\displaystyle+\mathrm{vec}(\bm{\Phi}_{l}^{H})^{H}\left(\mathbf{P}_{l}^{H}\mathbf{P}_{l}\otimes\bm{\Xi}_{l}^{-1}\right)\mathrm{vec}(\mathbf{G}^{H})$
$\displaystyle+\mathrm{vec}(\mathbf{G}^{H})^{H}\left(\mathbf{P}_{l}^{H}\mathbf{P}_{l}\otimes\bm{\Xi}_{l}^{-1}\right)\mathrm{vec}(\bm{\Phi}_{l}^{H}).$
(24)

Substituting (IV-C) into (23), we have

$\displaystyle
Q^{\dagger}(\mathrm{vec}(\mathbf{G}^{H}))\\!\propto\\!\prod_{l=1}^{L}\mathcal{CN}_{NK}\left(\mathrm{vec}(\mathbf{G}^{H})|\mathrm{vec}(\bm{\Phi}_{l}^{H}),(\mathbf{P}_{l}^{H}\mathbf{P}_{l})^{-1}\right.$
$\displaystyle\left.\otimes\bm{\Xi}_{l}\right)\times\mathcal{CN}_{NK}(\mathrm{vec}(\mathbf{G}^{H})|\mathbf{0},(\mathbb{E}[\bm{\Gamma}]+\mathbb{E}[\mathbf{Z}])^{-1}).$
(25)

Since (IV-C) consists of the product of two Gaussian pdfs, after applying
Lemma 1 in Appendix A, we have
$Q^{\dagger}(\mathrm{vec}(\mathbf{G}^{H}))=\mathcal{CN}_{NK}(\mathrm{vec}(\mathbf{G}^{H})|\mathbf{u},\bm{\Omega})$
with

$\displaystyle\\!\\!\\!\bm{\Omega}\\!$
$\displaystyle=\\!\left(\sum_{l=1}^{L}\mathbf{P}_{l}^{H}\mathbf{P}_{l}\otimes\bm{\Xi}_{l}^{-1}+\mathbb{E}[\bm{\Gamma}]+\mathbb{E}[\mathbf{Z}]\right)^{-1}\\!\\!\\!\\!\in\\!\mathbb{C}^{NK\times
NK},$ (26) $\displaystyle\mathbf{u}$
$\displaystyle=\bm{\Omega}\left(\sum_{l=1}^{L}\mathbf{P}_{l}^{H}\mathbf{P}_{l}\otimes\bm{\Xi}_{l}^{-1}\mathrm{vec}(\bm{\Phi}_{l}^{H})\right)\in\mathbb{C}^{NK}.$
(27)

Therefore, the pdf $Q^{\dagger}(\mathrm{vec}(\mathbf{G}))$ is a circularly
symmetric complex vector-valued normal distribution with mean $\mathbf{u}$ and
covariance $\bm{\Omega}$. An intuitive interpretation of (26) and (27) are
given as follows. The covariance matrix $\bm{\Omega}$ is updated by combining
the posterior information from other factor matrices
$\mathbf{X}_{l}^{i},\forall l,i$, the prior low-rank information
$\mathbb{E}[\bm{\Gamma}]$, and prior element-wise sparsity information
$\mathbb{E}[\mathbf{Z}]$. The tradeoff among these three terms is controlled
by the expectation of the noise precision $\mathbb{E}[\beta]$ in
$\bm{\Xi}_{l}$. In other words, the larger the current noise precision leads
to more information from other factors than that from prior information. On
the other hand, the channel mean vector $\mathbf{u}$ is rotated by the
covariance matrix $\bm{\Omega}$ to obtain the property of low-rank and
element-wise sparsity.

Once the mean vector in (27) is derived, the estimation of device-RIS
reflected channel matrix $\mathbf{M}_{G}$ can be obtained. Note that the
computation of $\bm{\Xi}_{l}$ in (23) and (27) is similar to (IV-D) and (30),
which will be introduced in the next section.

###  IV-D Learning to Separate Devices: Inferring $Q(\mathbf{X}_{l}^{i})$

After obtaining the expression of $Q(\mathbf{G})$, we need to derive the
update equation for each individual variational pdf $Q(\mathbf{X}_{l}^{i})$,
whose mean matrix is the estimated data message transmitted from the BS [38].
After substituting (IV-A) into (22) and keeping the terms relevant to
$\mathbf{X}_{l}^{i}$, the variational pdf $Q(\mathbf{X}_{l}^{i})$ can be
derived as follows:

$\displaystyle
Q^{\dagger}(\mathbf{X}_{l}^{i})\propto\exp\left\\{\mathbb{E}\left[-\beta\left\|\mathcal{Y}_{l}-\left[\\!\\!\left[\mathbf{X}_{l}^{1},~{}\mathbf{X}_{l}^{2},\cdots,\mathbf{X}_{l}^{d},~{}\mathbf{P}_{l}\mathbf{G}\right]\\!\\!\right]\right\|_{F}^{2}-\mathrm{Tr}(\bm{\Lambda}\mathbf{X}_{l}^{i,H}\mathbf{X}_{l})\right]\right\\}$
$\displaystyle\propto\exp\left\\{-\mathrm{Tr}\left(\mathbf{X}_{l}^{i}\right.\right.\left.\left.\underbrace{\mathbb{E}[\beta]\mathbb{E}\left[\left(\\!\mathop{\diamond}\limits_{j=1,j\neq
i}^{d}\mathbf{X}_{l}^{j}\mathop{\diamond}\limits\mathbf{P}_{l}\mathbf{G}\\!\right)^{T}\\!\left(\\!\mathop{\diamond}\limits_{j=1,j\neq
i}^{d}\mathbf{X}_{l}^{j}\mathop{\diamond}\limits\mathbf{P}_{l}\mathbf{G}\\!\right)^{*}\right]\\!+\\!\mathbb{E}[\bm{\Lambda}]}_{(\bm{\Sigma}_{l}^{i})^{-1}}\right.\right.$
$\displaystyle\left.\left.\times\\!\mathbf{X}_{l}^{i,H}\\!\\!-\\!\mathbf{X}_{l}^{i}(\bm{\Sigma}_{l}^{i})^{-1}\\![\underbrace{\mathbb{E}[\beta]\mathcal{Y}_{l}(i)\left(\mathop{\diamond}\limits\limits_{j=1,j\neq
i}^{d}\mathbb{E}[\mathbf{X}_{l}^{i}]\diamond\mathbf{P}_{l}\mathbb{E}[\mathbf{G}]\right)^{*}\\!\\!\bm{\Sigma}_{l}^{i}}_{\mathbf{M}_{l}^{i}}]^{H}\right.\right.\left.\left.-\mathbf{M}_{l}^{i}(\bm{\Sigma}_{l}^{i})^{-1}\mathbf{X}_{l}^{i,H}\right)\right\\}.$
(28)

It can be concluded that the variational pdf
$Q^{\dagger}(\mathbf{X}_{l}^{i})(1\leq l\leq L,1\leq i\leq d)$ is a circularly
symmetric complex matrix normal distribution $\mathcal{CN}_{N\times
K}(\mathbf{X}_{l}^{i}|\mathbf{M}_{l}^{i},\mathbf{1}_{\tau_{i}}\otimes\bm{\Sigma}_{l}^{i})$
with mean matrix $\mathbf{M}_{l}^{i}$ and covariance
$\mathbf{1}_{\tau_{i}}\otimes\bm{\Sigma}_{l}^{i}$.

Once the matrix $\mathbf{M}_{l}^{i}$ in (IV-D) at the top of next page is
derived, the estimation of the transmitted message can be obtained. All items
in $\mathbf{M}_{l}^{i}$ can be calculated directly, except for
$\bm{\Sigma}_{l}^{i}$. According to the property that
$(\mathop{\diamond}\limits_{j=1}^{d}\mathbf{A}^{j})^{T}(\mathop{\diamond}\limits_{j=1}^{d}\mathbf{A}^{j})^{*}=\mathop{\odot}\limits_{j=1}^{d}\mathbf{A}^{j,T}\mathbf{A}^{j,*}$
[22], the expectation of the coupled and individual factors terms in
$(\bm{\Sigma}_{l}^{i})^{-1}$ of (IV-D) reduces to

$\displaystyle\mathbb{E}\left[\left(\\!\mathop{\diamond}\limits_{j=1,j\neq
i}^{d}\mathbf{X}_{l}^{j}\mathop{\diamond}\limits\mathbf{P}_{l}\mathbf{G}\\!\right)^{T}\\!\left(\\!\mathop{\diamond}\limits_{j=1,j\neq
i}^{d}\mathbf{X}_{l}^{j}\mathop{\diamond}\limits\mathbf{P}_{l}\mathbf{G}\\!\right)^{*}\right]$
$\displaystyle=\mathcal{D}\left[\mathop{\odot}\limits_{j=1,j\neq
i}^{d}\mathbb{E}\left[\mathbf{X}_{l}^{j,T}\mathbf{X}_{l}^{j,*}\right]\odot\mathbb{E}\left[\mathbf{G}^{T}\mathbf{P}_{l}^{T}\mathbf{P}_{l}\mathbf{G}^{*}\right]\right],$
(29)

where $\mathop{\odot}\limits_{j=1,j\neq
i}^{d}\mathbf{A}^{j}=\mathbf{A}^{d}\odot\cdots\odot\mathbf{A}^{i+1}\odot\mathbf{A}^{i-1}\odot\mathbf{A}^{1}$
is the multiple Hadamard products. The second term
$\mathbb{E}\left[\mathbf{G}^{T}\mathbf{P}_{l}^{T}\mathbf{P}_{l}\mathbf{G}^{*}\right]$
on the right hand side (RHS) of (IV-D) can be derived by Lemma 2 in Appendix A
and the first expectation is given by

$\displaystyle\mathbb{E}\left[\mathbf{X}_{l}^{j,T}\mathbf{X}_{l}^{j,*}\right]=\mathbf{M}_{l}^{j,T}\mathbf{M}_{l}^{j,*}+{\tau_{j}}\bm{\Sigma}_{l}^{j}.$
(30)

###  IV-E Learning Hyper-Parameters: Inferring $Q(\bm{\xi})$, $Q(\bm{\eta})$,
$Q(\bm{\gamma})$, and $Q(\beta)$

Then, we infer the model parameters $\bm{\xi}(n,k)$ characterizing the
element-wise sparsity of the channel. After plugging the postulated joint pdf
(IV-A) into (22) and only remaining the terms about $\bm{\xi}(n,k)$, we have
$Q^{\dagger}({\\{\bm{\xi}(n,k)\\}_{n=1}^{N}}_{k=1}^{K})=\prod_{n=1}^{N}\prod_{k=1}^{K}Q^{\dagger}(\bm{\xi}(n,k))$
with

$\displaystyle
Q^{\dagger}(\bm{\xi}(n,k))\propto\exp\\{-\bm{\xi}(n,k)(\underbrace{\mathbb{E}\left[\mathbf{G}(n,k)^{*}\mathbf{G}(n,k)\right]+\delta}_{a_{\bm{\xi}(n,k)}})$
$\displaystyle+(\underbrace{N+\delta}_{b_{\bm{\xi}}}-1)\ln\bm{\xi}(n,k)\\}.$
(31)

Since the optimal $Q^{\dagger}(\bm{\xi}(n,k))$ in (IV-E) obeys a gamma
distribution of the form
$Q^{\dagger}(\bm{\xi}(n,k))=\mathrm{gamma}(\bm{\xi}(n,k)|a_{\bm{\xi}(n,k)},b_{\bm{\xi}})$,
the mean value of $\bm{\xi}(n,k)$ can be updated as
$\mathbb{E}[\bm{\xi}(n,k)]=b_{\bm{\xi}}/a_{\bm{\xi}(n,k)}$. Using the
correlation property of the normal distribution derived in (49), the
expectation in $a_{\bm{\xi}(n,k)}$ can be calculated as

$\\!\\!\mathbb{E}\left[\mathbf{G}(n,k)^{*}\mathbf{G}(n,k)\right]\\!=\\!\mathbf{M}_{G}(n,k)^{*}\mathbf{M}_{G}(n,k)\\!+\\!\bm{\Omega}_{n,n}^{b}(k,k).\\!\\!\\!\\!\\!$
(32)

Next, we infer the model parameters $\eta_{k}$ and $\gamma_{k}$ characterizing
the low rank of the tensor data $\mathcal{Y}_{l}$. Combining (IV-A) and (22)
and removing the terms independent of $\eta_{k}$, we obtain

$\displaystyle
Q^{\dagger}(\\{\eta_{k}\\}_{k=1}^{K})\propto\exp\\{\mathbb{E}[-\mathrm{Tr}(\bm{\Lambda}\mathbf{G}^{H}\mathbf{G})-\sum_{k=1}^{K}\delta\eta_{k}$
$\displaystyle+\left(N+\delta-1\right)\sum_{k=1}^{K}\ln{\eta_{k}}]\\},$ (33)

which can be reformulated as
$Q^{\dagger}(\\{\eta_{k}\\}_{k=1}^{K})=\prod_{k=1}^{K}Q^{\dagger}(\eta_{k})$
with

$\displaystyle\\!\\!\\!Q^{\dagger}(\eta_{k})\\!\propto\\!\exp\\{-\eta_{k}\underbrace{\mathbb{E}\left[\mathbf{G}(:,k)^{H}\mathbf{G}(:,k)+\delta\right]}_{a_{\eta_{k}}}\\!+\\!(\underbrace{N\\!+\\!\delta}_{b_{\eta}}-1)\ln\eta_{k}\\}.$
(34)

From (34), we know that the optimal $Q^{\dagger}(\eta_{k})$ obeys a gamma
distribution of the form $\mathrm{gamma}(\eta_{k}|a_{\eta_{k}},b_{\eta})$.
Taking the expectation in $a_{\eta_{k}}$ and utilizing (49) in Lemma 2 in
Appendix A, we have

$\displaystyle
a_{\eta_{k}}=\mathbf{M}_{G}(:,k)^{H}\mathbf{M}_{G}(:,k)\\!+\\!\\{\mathop{\sum}\limits_{n=1}^{N}\bm{\Omega}_{n,n}^{b}\\}(k,k)+\delta.$
(35)

Thus the expectation of parameter $\eta_{k}$ is given by
$\mathbb{E}(\eta_{k})=b_{\eta}/a_{\eta_{k}}$ with $b_{\eta}=N+\eta$. In a
similar way, the optimal $Q^{\dagger}(\gamma_{k})$ is found to be a gamma
distribution of the form
$\mathrm{gamma}(\gamma_{k}|a_{\gamma_{k}},b_{\gamma})$ with
$a_{\gamma_{k}}=\mathbb{E}\left[\sum_{l=1}^{L}\sum_{i=1}^{d}\mathbf{X}_{l}^{i}(:,k)^{H}\mathbf{X}_{l}^{i}(:,k)+\delta\right]$
and $b_{\gamma}=L\sum_{i=1}^{d}\tau_{i}+\delta$. Herein,
$\mathbb{E}\left[\mathbf{X}_{l}^{i}(:,k)^{H}\mathbf{X}_{l}^{i}(:,k)\right]=\mathbf{M}_{l}(:,k)^{i,H}\mathbf{M}_{l}(:,k)^{i}+{\tau_{i}}\bm{\Sigma}_{l}^{i}(k,k)$.
Then we can calculate the expectation of parameter $\gamma_{k}$ by
$\mathbb{E}(\gamma_{k})=b_{\gamma}/a_{\gamma_{k}}$.

Finally, the posterior distribution of $\beta$ is updated by the following
equation:

$\displaystyle
Q^{\dagger}(\beta)\propto(\underbrace{(\Pi_{i=1}^{d+1}\tau_{i}M)L+\delta}_{b_{\beta}}-1)\ln\beta$
$\displaystyle-\beta\underbrace{\sum_{l=1}^{L}\mathbb{E}\left[\left\|\mathcal{Y}_{l}-\left[\\!\\!\left[\mathbf{X}_{l}^{1},~{}\mathbf{X}_{l}^{2},\cdots,\mathbf{X}_{l}^{d},~{}\mathbf{P}_{l}\mathbf{G}\right]\\!\\!\right]\right\|_{F}^{2}+\delta\right]}_{a_{\beta}}.\\!\\!\\!$
(36)

By comparing (IV-E) to the functional form of gamma distribution, we have
$Q^{\dagger}(\beta)=\mathrm{gamma}(\beta|a_{\beta},b_{\beta})$. Herein,
$b_{\beta}$ in (IV-E) is related to the number of observations and $a_{\beta}$
approximates the residual of model fitting measured by the squared Frobenius
norm on observed entries.

For calculating $a_{\beta}$ in (IV-E), we apply the tensor unfolding and then
expand the Frobenius norm as follows

$\displaystyle\mathbb{E}\left[\left\|\mathcal{Y}_{l}-\left[\\!\\!\left[\mathbf{X}_{l}^{1},~{}\mathbf{X}_{l}^{2},\cdots,\mathbf{X}_{l}^{d},~{}\mathbf{P}_{l}\mathbf{G}\right]\\!\\!\right]\right\|_{F}^{2}\right]$
$\displaystyle=\mathrm{Tr}\left(\mathop{\odot}\limits_{i=1}^{d}\left(\mathbf{M}_{l}^{i,H}\mathbf{M}_{l}^{i,*}+{\tau_{i}}\bm{\Sigma}_{l}^{i}\right)^{H}\times\left(\mathbf{M}_{G}^{H}\mathbf{P}_{l}^{H}\mathbf{P}_{l}\mathbf{M}_{G}\right.\right.$
$\displaystyle\left.\left.+\mathop{\sum}\limits_{i=1}^{N}\mathop{\sum}\limits_{j=1}^{N}\\{\mathbf{P}_{l}^{H}\mathbf{P}_{l}\\}(i,j)\bm{\Omega}_{i,j}^{b}\right)+\left\|\mathcal{Y}_{l}(d+1)\right\|_{F}^{2}\right.$
$\displaystyle\left.-\mathbf{P}_{l}\mathbf{M}_{G}\left(\mathop{\diamond}\limits_{i=1}^{d}\mathbf{M}_{l}^{i}\right)^{T}\mathcal{Y}_{l}(d+1)^{H}\right.$
$\displaystyle\left.-\mathcal{Y}_{l}(d+1)\left(\mathop{\diamond}\limits_{i=1}^{d}\mathbf{M}_{l}^{i}\right)^{*}\mathbf{M}_{G}^{H}\mathbf{P}_{l}^{H}\right),$
(37)

where the RHS of (IV-E) is obtained by utilizing equation (50) in Lemma 2 in
Appendix A and the result in (30).

###### _Remark_ 1

(Learning of Active Device Number): During the iterations, the mean of the
model complexity parameters, i.e., $\mathbb{E}(\eta_{k})$ and
$\mathbb{E}(\gamma_{k})$ are learned by equations
$\mathbb{E}(\eta_{k}^{t+1})=b_{\eta}/a_{\eta_{k}}^{t+1}$ and
$\mathbb{E}(\gamma_{k}^{t+1})=b_{\gamma}/a_{\gamma_{k}}^{t+1}$, respectively,
together with those of other parameters. Since $\mathbb{E}(\eta_{k}^{-1})$ and
$\mathbb{E}(\gamma_{k}^{-1})$ have physical interpretations as the power of
each column in $\mathbf{G}$ and $\mathbf{X}_{l}^{i}$, if some
$\mathbb{E}(\eta_{k})$ and $\mathbb{E}(\gamma_{k})$ are large enough, it
indicates that their corresponding columns in $\mathbf{G}$ and
$\mathbf{X}_{l}^{i}$ can be safely pruned out. Then, the number of non-zero
columns in each factor matrix equals to the number of active devices.

###  IV-F Algorithm Summary

Since the statistics of each variational pdf rely on other variational pdfs.
Therefore, all parameters of variational pdfs need to be updated
alternatingly. For clarity, the pseudo-code of the resulting algorithm is
outlined in Algorithm 1.

Algorithm 1 Coupled Tensor-Based Automatic Detection (CTAD) Algorithm

Input: $\\{\mathcal{Y}_{l}\\}_{l=1}^{L}$ and total iterations $T$

Initialization: $\mathbf{M}_{G}^{0}$,
${\\{\mathbf{M}_{l}^{i,0}\\}_{l=1}^{L}}_{i=1}^{d}$, $\bm{\Sigma}_{G}^{0}$,
${\\{\bm{\Sigma}_{l}^{i,t}\\}_{l=1}^{L}}_{i=1}^{d}$, $a_{\beta}^{0}$,
$\\{a_{\eta_{k}}^{0},a_{\gamma_{k}}^{0}\\}_{k=1}^{K}$, and
${\\{a_{{\bm{\xi}}_{n,k}}^{0}\\}_{n=1}^{N}}_{k=1}^{K}$

for $t=1:T$ do

Updates the parameters of $Q(\mathbf{G})^{t+1}$:
$\displaystyle\bm{\Xi}_{l}^{t+1}$
$\displaystyle=\left(\frac{a_{\beta}^{t}}{b_{\beta}}\mathcal{D}\left[\mathop{\odot}\limits_{i=1}^{d}\left((\mathbf{M}_{l}^{i,t})^{H}\mathbf{M}_{l}^{i,t}+{\tau_{i}}\bm{\Sigma}_{l}^{i,t}\right)^{*}\right]\right)^{-1},$
(38)
$\displaystyle{\bm{\Omega}}^{t+1}=\left(\sum_{l=1}^{L}\mathbf{P}_{l}^{H}\mathbf{P}_{l}\otimes\left(\bm{\Xi}_{l}^{t+1}\right)^{-1}\right.$
$\displaystyle\left.+\text{diag}\left(\mathbf{1}_{N}\otimes\left[\frac{b_{\eta}}{a_{\eta_{1}}^{t}},\cdots,\frac{b_{\eta}}{a_{\eta_{K}}^{t}}\right]\right)\right.+\text{diag}\left(\frac{b_{\bm{\xi}}}{a_{\bm{\xi}_{1,1}}^{t}},\right.$
$\displaystyle\left.\left.\cdots,\frac{b_{\bm{\xi}}}{a_{\bm{\xi}_{1,K}}^{t}},\cdots,\frac{b_{\bm{\xi}}}{a_{\bm{\xi}_{N,1}}^{t}},\cdots,\frac{b_{\bm{\xi}}}{a_{\bm{\xi}_{N,K}}^{t}}\right)\right)^{-1},$
(39)
$\displaystyle\mathbf{u}^{t+1}={\bm{\Omega}}^{t+1}\mathop{\sum}\limits_{l=1}^{L}\mathbf{P}_{l}^{H}\mathbf{P}_{l}\otimes(\bm{\Xi}_{l}^{t+1})^{-1}$
$\displaystyle\mathrm{vec}\\!\left(\\!\left(\\!\frac{a_{\beta}^{t}}{b_{\beta}}(\mathbf{P}_{l}^{H}\mathbf{P}_{l})^{-1}\\!\mathbf{P}_{l}^{H}\mathcal{Y}_{l}(d+1)\left(\mathop{\diamond}\limits_{i=1}^{d}\mathbf{M}_{l}^{i,t}\right)^{*}\bm{\Xi}_{l}^{t+1}\\!\right)^{H}\\!\right)\\!,\\!\\!\\!\\!$
(40) $\displaystyle\mathbf{M}_{G}$
$\displaystyle=\mathrm{reshape}(\mathbf{u},N\times K).$ (41)

Updates the parameters of $Q(\mathbf{X}_{l}^{i})^{t+1}$:
$\displaystyle\bm{\Sigma}_{l}^{i,t+1}\\!=\\!\\!\left(\frac{a_{\beta}^{t}}{b_{\beta}}\mathcal{D}\\!\left[\\!\mathop{\odot}\limits_{j=1,j\neq
i}^{d}\\!\left(\\!(\mathbf{M}_{l}^{j,t})^{H}\mathbf{M}_{l}^{j,t}\\!+\\!{\tau_{j}}\bm{\Sigma}_{l}^{j,t}\right)^{*}\\!\right.\right.\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!$
$\displaystyle\left.\left.\\!\odot\\!\left(\\!\mathbf{M}_{G}^{t,H}\mathbf{P}_{l}^{H}\mathbf{P}_{l}\mathbf{M}_{G}^{t}+\mathop{\sum}\limits_{i=1}^{N}\mathop{\sum}\limits_{j=1}^{N}\\{\mathbf{P}_{l}^{H}\mathbf{P}_{l}\\}(i,j)\bm{\Omega}_{i,j}^{b,t}\\!\right)^{*}\right]\right.$
$\displaystyle\left.+\text{diag}\left(\frac{b_{\gamma}}{a_{\gamma_{1}}^{t}},\cdots,\frac{b_{\gamma}}{a_{\gamma_{K}}^{t}}\right)\right)^{-1},$
(42) $\displaystyle\mathbf{M}_{l}^{i,t}$
$\displaystyle=\frac{a_{\beta}^{t}}{b_{\beta}}\mathcal{Y}_{l}(i)\left(\mathop{\diamond}\limits\limits_{j=1,j\neq
i}^{d}\mathbf{M}_{l}^{i,t}\diamond\mathbf{P}_{l}\mathbf{M}_{G}^{t}\right)^{*}\bm{\Sigma}_{l}^{i,t+1}.$
(43)

Updates the parameters of $Q(\bm{\xi}(n,k))^{t+1}$:
$\displaystyle\\!\\!\\!\\!a_{\bm{\xi}(n,k)}^{t+1}=\mathbf{M}_{G}^{t+1,*}(n,k)\mathbf{M}_{G}^{t+1}(n,k)+\bm{\Omega}_{n,n}^{b,t+1}(k,k).$
(44)

Updates the parameters of $Q(\eta_{k})^{t+1}$:
$\displaystyle\\!\\!\\!a_{\eta_{k}}^{t+1}$
$\displaystyle\\!=\\!(\mathbf{M}_{G}^{t}(:,k)^{H}\mathbf{M}_{G}(:,k)^{t}\\!+\\![\mathop{\sum}\limits_{n=1}^{N}\bm{\Omega}_{n,n}^{b,t+1}](k,k))\\!+\\!\delta.$
(45)

Updates the parameters of $Q(\gamma_{k})^{t+1}$:  $\displaystyle
a_{\gamma_{k}}^{t+1}$
$\displaystyle=\sum_{l=1}^{L}\sum_{i=1}^{d}(\mathbf{M}_{l}^{i,t}(:,k))^{H}\mathbf{M}_{l}^{i,t}(:,k)+{\tau_{i}}\bm{\Sigma}_{l}^{i,t}(k,k)+\delta.$
(46)

Updates the parameters of $Q(\bm{\beta})^{t+1}$:  $\displaystyle
a_{\beta}^{t+1}=\sum_{l=1}^{L}\pounds_{l}^{t+1}+\delta,$ (47) where
$\pounds_{l}^{t+1}$ is calculated by (IV-E).

end for

Output: $\mathbf{M}_{G}^{t+1}$ and
${\\{\mathbf{M}_{l}^{i,t+1}\\}_{l=1}^{L}}_{k=1}^{K}$

We now provide the details of Algorithm 1: Eq. (38) denotes the variance of
factor matrix $\mathbf{G}$, and Eqs. (39) and (40) denote the variance and
mean of vector $\mathrm{vec}(\mathbf{G})$, respectively, which is derived in
Section IV-C. Eq. (41) denotes the mean in the matrix form of factor matrix
$\mathbf{G}$, which is reshaped from vector $\mathbf{u}$. Eqs. (5) and (43)
denote the variance and mean of factor matrix $\mathbf{X}$, respectively,
which are derived in Section IV-D. Eqs. (44)-(47) denote the hyper-parameters
of $\bm{\xi}$, $\bm{\eta}$, $\bm{\gamma}$, and $\beta$, respectively, which
are derived in Section IV-E. After applying Algorithm 1, we can obtain
${\\{\mathbf{M}_{l}^{i,t+1}\\}_{l=1}^{L}}_{k=1}^{K}$ and
$\mathbf{M}_{G}^{t+1}$, which correspond to the estimation of data
$\mathbf{x}_{i,k,l},\forall i,k,l$, and the channel sparsity profile
$\bm{\lambda}_{k},\forall k$. Note that scalar indeterminacy can be resolved
by designing each sub-constellation $\wp_{i}$ based on a Grassmannian codebook
[40].

###  IV-G Performance Analysis

To gain further insights from the above proposed algorithm, this subsection
analyzes its performance from the perspective of computational complexity and
the convergence property.

####  IV-G1 Computational Complexity

For each iteration, the computational complexity is dominated by updating each
factor matrix and mainly arises from the matrix multiplication, which is in
the order of
$\mathcal{O}(L(\mathop{\sum}\limits_{i=1}^{d}\tau_{i}+M)K^{3}+Ld\mathop{\prod}_{i=1}^{d}M\tau_{i}K^{2}+LN^{2}M)$.
It can be seen that the complexity of the proposed CTAD algorithm scales
linearly with the subblocks but only polynomially with the potential number of
active devices. In contrast, the computational complexity of the active device
separation algorithm in [10] scales polynomially with the total number of
devices, i.e., $\mathcal{O}(\bar{K}^{4})$. Thus, the proposed CTAD algorithm
is computational efficiency.

####  IV-G2 Convergence Property

The functional minimization of the KL divergence in (IV-B) is convex with
respect to a single variational pdf $Q(\Theta_{j})$ when the others
$\\{Q(\Theta_{j})\\}_{i\neq j}$ are fixed. This guarantees monotonic decrease
of the KL divergence in (IV-B) over iterations. In addition, the proposed CTAD
algorithm is an instance of the block coordinate descent optimization strategy
over the functional space. Therefore, the CTAD algorithm is able to decrease
the objective in every iteration. Since the proposed algorithm is developed
under the mean-field variational inference framework [39], the limit point of
the proposed algorithm is guaranteed to converge to a stationary point of the
objective function in (IV-B).

####  IV-G3 Identifiability Property

Uniqueness is a key characteristic of the proposed tensor decomposition
problem. If $LM\geq 2\varpi$ and
$2K_{a}+d\leq\min\\{\kappa_{1},\kappa_{2},\cdots,\kappa_{L}\\}$ with
$\kappa_{l}=\sum_{i=1}^{d}r_{\mathbf{X}_{l}^{i}}+\min(M,r_{\mathbf{G}})$, then
$\mathbf{X}_{l}^{i}$ and $\mathbf{G}$ are almost surely identifiable from the
received data $\mathcal{Y}_{l},l=1,2,\cdots,L$. Herein,
$r_{\mathbf{X}_{l}^{i}}$ is the Kruskal-rank of each factor matrix
$\mathbf{X}_{l}^{i}$ and $\varpi$ is an upper bound on the number of nonzero
elements in each column of $\mathbf{G}$.

_Remark 2_ : In this paper, we consider a scenario of Rayleigh fading channels
without line-of-sight (LoS) paths between the devices and the RIS. In fact,
due to the massive number of devices in a cell, there may exist LoS paths
between the devices and the RIS. In such a scenario, the proposed tensor-based
URA scheme is also applicable by redesigning the prior for matrix $\mathbf{G}$
in (IV-A). Hence, the proposed URA scheme has a wide application.

##  V Numerical Results

In this section, we present extensive simulation results to validate the
effectiveness of the proposed CTAD algorithm in 6G wireless networks. For the
URA scheme, we consider the case that every active device embeds its ID in the
payload, in addition to the $B=360$ information bits [8]. It is assumed that
there are total of $\bar{K}=4096$ devices in the network, therefore,
$B_{\mathrm{ID}}=\log_{2}(K)=12$ bits are required to encode the device ID.
Hence, each device transmits a total of
$B_{\mathrm{tot}}=B_{\mathrm{ID}}+B=372$ bits. At the BS, the active devices
are identified based on the ID embedded in the decoded payload. The active
device separation error is evaluated by the packet error rate (PER) metric
defined as $\mathrm{PER}=\frac{\bar{K}_{a}}{{K_{a}}}$, where $\bar{K_{a}}$ is
the number of incorrectly decoded messages. If multiple messages with the same
device ID are decoded, they are discarded and counted as an error. We assume
that the pilot sequences obey an independent and identically distributed
(i.i.d.) Gaussian distribution with zero mean and unit variance. Every
$\epsilon_{i}^{k}$ of the BS-device channel is drawn from the complex Gaussian
distributed random variables with zero mean and unit variance. The distance
$d_{\mathrm{DR},k}$ from the $k$th device to the RIS is randomly generated
from $500$ m to $1000$ m and from the RIS to the BS is $100$ m. Note that the
models in (1) and (2) are 3D channel models. We draw the central azimuth and
elevation AoD from the RIS (AoA to the RIS) of each cluster uniformly over
$[-180^{\circ},180^{\circ}]$ and $[-90^{\circ},90^{\circ}]$, respectively;
draw the central azimuth AoA at the BS of each cluster uniformly over
$[-90^{\circ},90^{\circ}]$. The angular spread of each subpath is set as
$15^{\circ}$. The over-complete bases $\bm{\nu}$ and $\bm{\varsigma}$ in (3)
are uniform sampling grids covering $[-1,1]$. The channel estimation accuracy
is evaluated in terms of normalized mean square error (NMSE) given by
$\mathrm{NMSE}=\frac{1}{T_{r}}\mathop{\sum}\limits_{i=1}^{T_{r}}\frac{\left\|{\mathbf{G}}^{i}-\hat{\mathbf{G}}^{i}\right\|_{F}^{2}}{\left\|{\mathbf{G}}^{i}\right\|_{F}^{2}}$,
where $T_{r}$ is the number of Monte Carlo runs and $\hat{\mathbf{G}}^{i}$ is
the device-RIS channel estimated at the $i$-th run. The tree-based decoder [8]
is employed in the simulation, where the $B$-bit message of each device is
divided into $L$ subblocks of size $B_{1},\cdots,B_{L}$ satisfying the
following conditions: $\sum_{l}B_{l}=B$, $B_{1}=R$, and $B_{l}<R$ for all
$i=2,\cdots,L$. Herein, all subblocks $i=2,\cdots,L$ are augmented to size $R$
by appending the parity bits which are set to $p_{b}(0,168,270)$. Unless
stated otherwise, we set $R=270$,
$N_{1}^{\prime}/N_{1}=N_{2}^{\prime}/N_{2}=1.5$, $\tau_{p}=N$,
$\delta=10^{-6}$, $d=2$, $\tau_{1}=80$, and $\tau_{2}=80$.

For comparison, we consider a baseline scheme which employs a two-stage
estimator for addressing the optimization problem in (IV): First, the
alternating least square (ALS) algorithm [21] is applied to fit the model in
(IV) and $\bar{\mathbf{G}}_{l}=\mathbf{P}_{l}\mathbf{G}$ is regarded unknown
to estimate. As for the second stage, we compute the minimum $l_{1}$-norm
solutions for the sparse factor matrices $\mathbf{G}$ from the estimation of
${\bar{\mathbf{G}}}_{l}$ and use the existing convex optimization solver [34]
to solve this least absolute shrinkage and selection operator (LASSO) problem.
Since the number of active devices $K_{a}$ is unknown, we set the initial
upper bounds of $K_{a}$ as $K$ in the two-stage estimator.

Figure 6: The estimation accuracy of the number of active devices for
different transmit power with $K_{a}=50$, $N=500$, $\bar{K}=4096$. Figure 7:
The tensor rank estimate versus number of iterations with $K_{a}=50$, $N=400$,
transmit power $=15$ dBm, and $\bar{K}=4096$.

When measuring the rank quality of estimation, we consider the mean value of
the rank estimation and probability of the successful recovery in the form of
$\mathrm{Pr}(K=K_{a})$. Each simulation is repeated 200 times to obtain the
averaged result. Fig. 6 suggests that for $M=300$ and the data transmit power
$\geq 5$ dBm, the proposed CTAD algorithm can recover the true tensor rank
with probability 1 in the considered simulation setting. This shows the high
accuracy and robustness of the proposed algorithm when the noise power is
moderate. For $M=350$, a cut-off point for accurate rank estimation is relaxed
to transmit power $\geq-5$ dBm, which indicates that the success rate of rank
estimation is increased as the antenna number $M$ at the BS grows in the lower
transmit power region. This is due to the fact that $\mathbf{P}_{l}$ is
equivalent to the measurement matrix in compressed sensing and the increasing
number of BS antennas results in the increment of measurement length for a
better observation. Fig. 7 shows the mean value of estimated rank along with
the number of iterations. It can be seen that only a few iterations are
required for the convergence of the proposed algorithm on average. In
particular, the proposed algorithm has gradually recovered the number of
active devices, while the estimated rank remains stable. This is because, as
illustrated in Remark 1, the regularization parameters of these active devices
are kept over the next iteration, while the values of the other devices are
pruned out by thresholding.

Figure 8: PER versus different tensor sizes $d$ and $\tau_{i}$ with $M=350$,
$N=500$, $L=3$, $K_{a}=15$, and $\bar{K}=4096$. Figure 9: PER versus different
$K$ with $K_{a}=50$, $M=256$, transmit power $=10$ dBm, and $\bar{K}=4096$.
Figure 10: The convergence performance with $K_{a}=50$, $M=350$, $N=500$,
transmit power $=15$ dBm, and $\bar{K}=4096$.

Fig. 8 plots the PER against transmit power for varying tensor sizes in an
RIS-aided massive unsourced access system. It can be seen that URA scheme in
[36] cannot recover any active device without the RIS under unfavorable
propagation conditions. The reason is that we assume that the direct link
between the BS and the devices is completely blocked and unavailable, and
hence the BS cannot receive the signal sending from the devices. However, this
issue can be resolved by deploying an RIS in URA networks. Moreover, we can
observe that the optimized/random combined phase shift achieves significant
performance gains over the proposed CTAD with fully random phase shifts.
Again, this confirms the effectiveness of the optimized/random combined
beamforming design for RIS-aided URA communications. Note that the PER of the
CTAD algorithm with $d=4$ and $d=3$ are higher than that of the CTAD algorithm
with $d=2$. This can be explained through the number of DoF per active device.
Since a variable in Grassmannian of lines in dimension $\tau_{i}$ has
$\tau_{i}-1$ DoF [40], the sum-DoF of the active devices in the considered
model can be defined as
$\mathrm{DoF}(K_{a})=K_{a}\sum\limits_{i=1}^{d}(\tau_{i}-1)$. As such, for a
fixed number of active devices $K_{a}$, we can see that the available DoF of
$d=4$ and $d=3$ are lower than that of $d=2$, which makes $d=2$ more suitable
for RIS-aided URA communications.

In Fig. 9, we provide the PER performance curves of the proposed CTAD
algorithm under different initial upper bounds $k$ of $K_{a}$ available, where
the genie-aided two-stage approach with exact $K_{a}$ and the two-stage
approach with incorrect active device numbers $K$ are set as benchmarks.
First, it can be observed that the proposed CTAD algorithm can offer the best
PER results due to its superior capabilities in determining the exact tensor
rank, i.e., $K_{a}$. In contrast, the two-stage approach with $K$ overfits the
noise heavily and the performance degrades severely in channel estimation.
Second, the genie-aided two-stage approach with $K_{a}$ still performs worse
than the proposed CTAD algorithm. The reason is that the proposed CTAD
algorithm exploits the prior statistical distributions of the sparse channel,
which is ignored by the two-stage approach. Third, increasing the number of
reflecting elements causes a PER degradation, as the dimension of channel
estimates improves with an increment in the reflecting elements.

Then, we show the convergence behaviour of the proposed CTAD algorithm in Fig.
10 by setting $\bm{\xi}=\mathbf{0}$ and $\bm{\xi}\neq\mathbf{0}$,
respectively. It is observed that the NMSE of the proposed algorithm with
$\bm{\xi}\neq\mathbf{0}$ decreases quickly with the increment of the number of
iterations and the algorithm converges to stable values in less than 25
iterations on average. Also, Fig. 10 verifies that the proposed algorithm
significantly outperforms the algorithm that does not explore the element-wise
sparsity ($\bm{\xi}=\mathbf{0}$). The reason is that when
$\bm{\xi}=\mathbf{0}$, although the algorithm searches for a low-rank
structure by automatically adjusting parameters $\eta_{k}$ and $\gamma_{k}$,
it is not able to explore element sparsity as each column of the channel
matrix shares the same hyperparameter. Hence, the estimation accuracy of model
parameters is unsatisfactory, which in turn degrades its ability to decompose
tensors. In contrast, the proposed CTAD algorithm with
$\bm{\xi}\neq\mathbf{0}$ infers the dynamics of the network by using the power
of columns of channel.

Figure 11: NMSE versus $L$ with $M=256$, $N=350$, transmit power $=15$ dBm,
and $\bar{K}=4096$. Figure 12: The NMSE versus sampling resolution $\iota$
with $M=250$, $N=500$, $K_{a}=50$, and $\bar{K}=4096$. Figure 13: PER versus
transmit power under the case where the direct channel link is not completely
blocked.

Next, we focus on the NMSE performance of the proposed CTAD algorithm with
different numbers of subblocks $L$ and a fixed rate $B/\tau$ in Fig. 11. As
can be seen from the obtained results, utilizing a single subblock $L=1$
cannot accurately estimate the channel, since $\mathbf{G}$ does not admit a
unique solution from the compressed model $\mathbf{P}_{l}\mathbf{G}$. While
the proposed CTAD algorithm is quite accurate for both $L=3,4$ and the
performance degrades with decreasing $L$. From this simulation result, we can
see that the coupled tensor factorization criterion in problem (IV) is
critical to the reconstruction of the channel, since the shared parameter
$\mathbf{G}$ in the $L$ fitting terms serves as an anchor to fix the
permutation and scaling ambiguities. As the difference in PER performances of
$L=3$ and $L=4$ is relatively smaller, $L$ is set to 3 in the rest of the
simulations for unveiling the full potential of the CTAD algorithm under
various system settings.

In Fig. 12, we study the effect of the grid lengths $N_{1}^{\prime}$ and
$N_{2}^{\prime}$. We adopt $\iota=N_{1}^{\prime}/N_{1}=N_{2}^{\prime}/N_{2}$
to represent the ratio between the grid length and the number of reflecting
elements. It is seen that the NMSE decreases as $\iota$ increases, since
increasing the sampling grid length leads to a higher angular resolution and
hence results in a sparser $\mathbf{G}$. The proposed algorithm substantially
decreases the NMSE compared with the two-stage approach, and the performance
gaps among the two algorithms decrease as the ratio decreases. This is because
the NMSE performance is limited by the sparsity of the channel. Besides, the
performances of all estimation schemes increase as the transmit power grows
due to a better received signal quality at the BS.

Finally, we compare the proposed CTAD algorithm and the URA schemes designed
in [8] and [36]. In Fig. 13, we consider a scenario that the direct channel
link between the BS and devices is not completely blocked such that a weak
signal can reach the BS through the direct link. It is assumed that the direct
channel link obeys an i.i.d. Gaussian distribution with zero mean and $0.3$
variance. Without the aid of RIS, the URA schemes designed in [8] and [36]
perform active device detection only with weak direct link signals.
Accordingly, the proposed CTAD algorithm performs active device separation
only based on IRS-associated signals. It can be seen that the PER performance
of the proposed algorithm is significantly better than that of URA schemes
designed in [8] and [36]. There are two main reasons, the first is that the
proposed algorithm exploits the prior statistical distributions of the sparse
channel, which is ignored by the algorithms in [8] and [36]. On the other
hand, the RIS can reflect impinging electromagnetic waves with a controllable
phase shift via the help of a smart controller. Through the passive
beamforming designed in Section III-A, the RIS is able to provide an extra
high-quality channel link to overcome the unfavourable propagation conditions
of URA systems.

##  VI Conclusion

In this paper, we proposed a new RIS-aided URA architecture and formulated
joint active device separation and channel estimation as a coupled high-order
tensor problem. To solve the problem, a novel probabilistic modeling and an
automatic learning algorithm were proposed under the framework of Bayesian
inference. It is found that the deployment of a RIS can effectively improve
the performance of URA in 6G wireless networks.

##  Appendix A Lemma 1 and Lemma2

In this section, we provide some results on normal distribution and
$\mathbf{G}$, which will be used for the derivation of individual variables in
Section IV.

_Lemma 1_ : Given the vector-valued normal distribution
$\mathcal{CN}_{NK}(\mathbf{x}|\mathbf{a}_{i},\mathbf{A}_{i})$ with the mean
vector $\mathbf{a}_{i}$ and covariance $\mathbf{A}_{i}$ for $i=1,\cdots,n$, we
have

$\displaystyle\\!\\!\mathop{\prod}\limits_{i=1}^{n}\mathcal{CN}_{NK}(\mathbf{x}|\mathbf{a}_{i},\mathbf{A}_{i})\propto\exp\left[-\frac{1}{2}\mathbf{x}^{H}\left(\mathop{\sum}\limits_{i=1}^{n}\mathbf{A}_{i}^{-1}\right)\mathbf{x}\\!-\\!\mathbf{x}^{H}\right.$
$\displaystyle\\!\\!\left.\times\left(\\!\mathop{\sum}_{i=1}^{n}\mathbf{A}_{i}^{-1}\mathbf{a}_{i}\\!\right)\\!-\\!\left(\\!\mathop{\sum}_{i=1}^{n}\mathbf{A}_{i}^{-1}\mathbf{a}_{i}\\!\right)^{H}\mathbf{x}\right]\\!\\!=\\!\mathcal{CN}_{NK}(\mathbf{x}|\mathbf{c},\mathbf{C}),$
(48)

with
$\mathbf{C}=\left(\mathop{\sum}\limits_{i=1}^{n}\mathbf{A}_{i}^{-1}\right)^{-1}$
and
$\mathbf{c}=\mathbf{C}\left(\mathop{\sum}_{i=1}^{n}\mathbf{A}_{i}^{-1}\mathbf{a}_{i}\right)$.

_Lemma 2_ : If $\mathrm{vec}(\mathbf{G}^{H})$ obeys the vector-valued normal
distribution
$\mathcal{CN}_{NK}(\mathrm{vec}(\mathbf{G}^{H})|\mathbf{u},\bm{\Omega})$, then
$\mathbf{G}$ is said to follow the matrix-variate normal distribution
$\mathcal{CN}_{N\times K}(\mathbf{G}|\mathbf{M}_{G},\bm{\Omega})$ with mean
matrix $\mathbf{M}_{G}$ and covariance matrix $\bm{\Omega}$. Herein, the mean
matrix $\mathbf{M}_{G}\in\mathbb{C}^{N\times K}$ is obtained by rearranging
vector $\mathbf{u}$. This further implies that

$\displaystyle\mathbb{E}[\mathbf{G}^{H}\mathbf{G}]$
$\displaystyle=\mathbf{M}_{G}^{H}\mathbf{M}_{G}+\mathop{\sum}\limits_{n=1}^{N}\bm{\Omega}_{n,n}^{b},$
(49)
$\displaystyle\mathbb{E}[\mathbf{G}^{H}\mathbf{P}_{l}^{H}\mathbf{P}_{l}\mathbf{G}]$
$\displaystyle=\mathbf{M}_{G}^{H}\mathbf{P}_{l}^{H}\mathbf{P}_{l}\mathbf{M}_{G}+\mathop{\sum}\limits_{i=1}^{N}\mathop{\sum}\limits_{j=1}^{N}\\{\mathbf{P}_{l}^{H}\mathbf{P}_{l}\\}(i,j)\bm{\Omega}_{i,j}^{b},$
(50)

where $\bm{\Omega}_{i,j}^{b}\in\mathbb{C}^{K\times K}$ is the $(i,j)$th block
of $\bm{\Omega}$.

## References

  * [1] X. Shao, L. Cheng, X. Chen, C. Huang, and D. W. K. Ng, “A Bayesian tensor approach to enable RIS for 6G massive unsourced random access,” _IEEE Global Commun. Conf. (GLOBECOM)_ , Madrid, Spain, Dec. 2021.
  * [2] H. Xiao, W. Chen, J. Fang, B. Ai, and I. J. Wassell, “A grant-free method for massive machine-type communication with backward activity level estimation,” _IEEE Trans. Signal Process._ , vol. 68, pp. 6665-6680, Dec. 2020.
  * [3] X. Shao, X. Chen, D. W. K. Ng, C. Zhong, and Z. Zhang, “Cooperative active device separation: Sourced and unsourced massive random access paradigms,” _IEEE Trans. Signal Process._ , vol. 68, pp. 6578-6593, Dec. 2020.
  * [4] M. Ke, Z. Gao, Y. Wu, X. Gao, and R. Schober, “Compressive sensing-based adaptive active user detection and channel estimation: Massive access meets massive MIMO,” _IEEE Trans. Signal Process._ , vol. 68, pp. 764-779, Jan. 2020.
  * [5] X. Shao, X. Chen, C. Zhong, J. Zhao, and Z. Zhang, “A Unified design of massive access for cellular Internet of Things,” _IEEE Intern. Things J._ , vol. 6, no. 2, pp. 3934-3947, 2019.
  * [6] Y. Polyanskiy, “A perspective on massive random-access,” in _Proc. IEEE Int. Symp. Inf. Theory_ , Aachen, Germany, Jun. 2017, pp. 2523-2527.
  * [7] A. Fengler, S. Haghighatshoar, P. Jung, and G. Caire, “Non-Bayesian activity detection, large-scale fading coefficient estimation, and unsourced random access with a massive MIMO receiver,” _IEEE Trans. Inform. Theory_ , vol. 67, no. 5, pp. 2925-2951, May 2021.
  * [8] A. Fengler, G. Caire, P. Jung, and S. Haghighatshoar, “Massive MIMO unsourced random access,” Jan. 2019, [Online]. Available: http://arxiv.org/abs/1901.00828.
  * [9] X. Xie, Y. Wu, J. Gao, and W. Zhang, “Massive unsourced random access for massive MIMO correlated channels”, _IEEE Global Commun. Conf. (GLOBECOM)_ , Taipei, Taiwan, Dec. 2021.
  * [10] S. Xia and Y. Shi, “Intelligent reflecting surface for massive device connectivity: Joint active device separation and channel estimation,” _IEEE Intern. Conf. Acoustics, Speech and Signal Process. (ICASSP)_ , Barcelona, Spain, 2020, pp. 5175-5179.
  * [11] X. Shao, X. Chen, C. Zhong, and Z. Zhang, “Exploiting simultaneous low-rank and sparsity in delay-angular domain for Millimeter-Wave/Terahertz wideband massive access,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 4, pp. 2336-2351, Apr. 2022.
  * [12] X. Yu, D. Xu, Y. Sun, D. W. K. Ng, and R. Schober, “Robust and secure wireless communications via intelligent reflecting surfaces,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 11, pp. 2637-2652, Nov. 2020.
  * [13] Q. Wu, S. Zhang, B. Zheng, C. You, and R. Zhang, “Intelligent reflecting surface aided wireless communications: A tutorial,” _IEEE Trans. Commun._ , vol. 69, no. 5, pp. 3313-3351, May 2021.
  * [14] C. Huang, S. Hu, G. C. Alexandropoulos, A. Zappone, C. Yuen, R. Zhang, R. D. Renzo, and M. Debbah, “Holographic MIMO surfaces for 6G wireless networks: Opportunities, challenges, and trends,” _IEEE Wireless Commun._ , vol. 27, no. 5, pp. 118-125, Oct. 2020.
  * [15] J. Yao, Z. Zhang, X. Shao, _et al_ , “Concentrative intelligent reflecting surface aided computational imaging via fast block sparse Bayesian learning,” in _Proc. IEEE Veh. Technol. Conf. (VTC)_ , Jun. 2021, pp. 1-6.
  * [16] C. Huang, A. Zappone, G. C. Alexandropoulos, M. Debbah, and C. Yuen, “Reconfigurable intelligent surfaces for energy efficiency in wireless communication,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 8, pp. 4157-4170, Aug. 2019.
  * [17] X. Shao, C. You, W. Ma, X. Chen, and R. Zhang, “Target sensing with intelligent reflecting surface: Architecture and performance,” _IEEE J. Sel. Areas Commun._ , vol. PP, no, 99, pp. 1-1, Mar. 2022.
  * [18] H. Liu, X. Yuan, and Y. -J. A. Zhang, “Matrix-calibration-based cascaded channel estimation for reconfigurable intelligent surface aided multiuser MIMO,” _IEEE J. Sel. Areas Commun.,_ vol. 38, no. 11, pp. 2621-2636, Nov. 2020.
  * [19] M. S. Rensen and L. De Lathauwer, “Coupled canonical polyadic decompositions and (coupled) decompositions in multilinear rank-($L_{r,n},L_{r,n},1$) terms-part I: Uniqueness” _SIAM J. Matrix Analysis and Appli._ , vol. 36, no. 2, pp. 496-522, 2015.
  * [20] L. De Lathauwer and E. Kofidis, “Coupled matrix-tensor factorizations-the case of partially shared factors” _Asilomar Conference on Signals, Systems, and Computers._ 2017, pp. 711-715.
  * [21] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” _SIAM Rev._ , vol. 51, pp. 455-500, 2009.
  * [22] L. Cheng, Y.-C. Wu, and H. V. Poor, “Probabilistic tensor canonical polyadic decomposition with orthogonal factors,” _IEEE Trans. Signal Process.,_ vol. 65, no. 3, pp. 663-676, Feb., 2017.
  * [23] Q. Zhao, L. Zhang, and A. Cichocki, “Bayesian CP actorization of incomplete tensors with automatic rank determination” _IEEE Trans. Pattern Analy. Machine Intell._ , vol. 37, no. 9, pp. 1751-1763, Sept. 2015.
  * [24] L. Cheng and Q. Shi, “Towards overfitting avoidance: Tuning-free tensor-aided multi-user channel estimation for 3D massive MIMO communications,” _IEEE J. Sel. Topics Signal Process._ , vol. 15, no. 3, pp. 832-846, Apr. 2021.
  * [25] L. Cheng, X. Tong, S. Wang, Y. Wu, and H. V. Poor, “Learning nonnegative factors from Tensor data: Probabilistic modeling and inference algorithm,” _IEEE Trans. Signal Process._ , vol. 68, pp. 1792-1806, Feb. 2020.
  * [26] W. U. Bajwa, J. Haupt, A. M. Sayeed, and R. Nowak, “Compressed channel sensing: A new approach to estimating sparse multipath channels,” _Proc. IEEE_ , vol. 98, no. 6, pp. 1058-1076, Jun. 2010.
  * [27] A. L. Swindlehurst, E. Ayanoglu, P. Heydari, and F. Capolino, “Millimeter-wave massive MIMO: The next wireless revolution?” _IEEE Commun. Mag._ , vol. 52, no. 9, pp. 56-62, Sep. 2014.
  * [28] X. Shao, X. Chen, and R. Jia, “A dimension reduction-based joint activity detection and channel estimation algorithm for massive access,” _IEEE Trans. Signal Process._ , vol. 68, pp. 420-435, Jan. 2020.
  * [29] B. Zheng, C. You, W. Mei, and R. Zhang, “A survey on channel estimation and practical passive beamforming design for intelligent reflecting surface aided wireless communications,” _IEEE Commun. Surv. Tuts._ , vol. 24, no. 2, pp. 1035-1071, Jun. 2022.
  * [30] C. Hu, L. Dai, S. Han, and X. Wang, “Two-timescale channel estimation for reconfigurable intelligent surface aided wireless communications,” _IEEE Trans. Commun._ , vol. 69, no. 11, pp. 7736-7747, Nov. 2021.
  * [31] X. Yu, V. Jamali, D. Xu, D. W. K. Ng, and R. Schober, “Smart and reconfigurable wireless communications: From IRS modeling to algorithm design,” _IEEE Wireless Commun._ , vol. 28, no. 6, pp. 118-125, Dec. 2021.
  * [32] F. Laue, V. Jamali, and R. Schober, “IRS-assisted active device separation,” _International Workshop Signal Process. Advances Wireless Commun. (SPAWC)_ , July 2022.
  * [33] E. Björnson and L. Sanguinetti, “Rayleigh fading modeling and channel hardening for reconfigurable intelligent surfaces,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 4, pp. 830-834, Apr. 2021.
  * [34] S. Boyd and L. Vandenberghe, _Convex Optimization_ , New York, NY, USA: Cambridge University Press, 2004.
  * [35] Z.-Q. Luo, W.-K. Ma, A. M. So, Y. Ye, and S. Zhang, “Semidefinite relaxation of quadratic optimization problems,” _IEEE Signal Process. Mag._ , vol. 27, pp. 20-34, May 2010.
  * [36] A. Decurninge, I. Land, and M. Guillaud, “Tensor-based modulation for unsourced massive random access,” _IEEE Wireless Commun. Let._ , vol. 10, no. 3, pp. 552-556, 2020.
  * [37] L. Chen and X. Yuan, “Blind multiuser detection in massive MIMO channels with clustered sparsity,” _IEEE Wireless Commun. Lett._ , vol. 8, no. 4, pp. 1052-1055, Aug. 2019.
  * [38] S. M. Kay, _Fundamentals of Statistical Signal Srocessing_ , Prentice Hall PTR, 1993.
  * [39] M. J. Wainwright and M. I. Jordan, “Graphical models, exponential families, and variational inference,” _Found. Trends Mach. Learn.,_ vol. 1, no. 102, pp. 1-305, Jan. 2008.
  * [40] K. Ngo, A. Decurninge, M. Guillaud, and S. Yang, “Cube-Split: A structured Grassmannian constellation for non-coherent SIMO communications,” _IEEE Trans. Wireless Commun.,_ vol. 19, no. 3, pp. 1948-1964, Mar. 2020.