# Towards Overfitting Avoidance: Tuning-free Tensor-aided Multi-user Channel
Estimation for 3D Massive MIMO Communications

Lei Cheng and Qingjiang Shi L. Cheng is with the Shenzhen Research Institute
of Big Data, Shenzhen, Guangdong, P. R. China (e-mail: leicheng@sribd.cn).Q.
Shi is with the School of Software Engineering at Tongji University, Shanghai,
201804, China. He is also with the Shenzhen Research Institute of Big Data,
Shenzhen, 518172, China. Email<EMAIL_ADDRESS>
###### Abstract

Channel estimation has long been deemed as one of the most critical problems
in three-dimensional (3D) massive multiple-input multiple-output (MIMO), which
is recognized as the leading technology that enables 3D spatial signal
processing in the fifth-generation (5G) wireless communications and beyond.
Recently, by exploring the angular channel model and tensor decompositions,
the accuracy of single-user channel estimation for 3D massive MIMO
communications has been significantly improved given a limited number of pilot
signals. However, these existing approaches cannot be straightforwardly
extended to the multi-user channel estimation task, where the base station
(BS) aims at acquiring the channels of multiple users at the same time. The
difficulty is that the coupling among multiple users’ channels makes the
channel estimation deviate from widely-used tensor decompositions. It gives a
non-standard tensor decomposition format that has not been well tackled. To
overcome this challenge, besides directly fitting the new tensor model for
channel estimation to the wireless data via block coordinate descent (BCD)
method, which is prone to the overfitting of noises or requires regularization
parameter tuning, we further propose a novel tuning-free channel estimation
algorithm that can automatically control the channel model complexity and thus
effectively avoid the overfitting. Numerical results are presented to
demonstrate the excellent performance of the proposed algorithm in terms of
both estimation accuracy and overfitting avoidance.

###### Index Terms:

Joint model-and-data-driven wireless communications, 3D massive MIMO, channel
estimation, tensor methods, tuning-free.

##  I Introduction

In recent years, massive multiple-input multiple-output (MIMO) has gradually
evolved from being a theoretical concept to a leading practical technology for
the next generation wireless communications [1, 2]. To further embrace the
forthcoming era of Internet of Things (IoT), it calls for advanced three-
dimensional (3D) spatial signal processing techniques (e.g., 3D beamforming)
[5, 6] to allow high-quality communications among multiple users (including
unmanned aerial vehicles (UAVs) in the sky [3] and unmanned ground vehicles
(UGVs) on roads [4]). To achieve this, rather than mechanically tilting
conventional antenna arrays, 3D massive MIMO, in which the BS is with a 3D
antenna array, has emerged as an enabling technique to broaden the scope of BS
[7, 8, 9]. However, its promise can only be fulfilled when accurate channel
state information (CSI) of multiple users is available at the BS.

As one of the most critical problems in wireless communications, channel
estimation has been continuously studied for many decades, synergizing nearly
all the ideas of signal processing methods including optimizations [10],
statistics [11], algebras [12], and machine learning [13]. Its ultimate goal
is to estimate wireless channels as accurate as possible given a limited
number of pilot signals, while its challenges vary significantly under
different channel models and wireless systems including millimeter-wave
massive MIMO [47], IoT [48] and machine-type communications (MTC) [49]. That
is why “no free lunch theorem” [14] in machine learning also holds in channel
estimation research, in the sense that there is no panacea that suits every
wireless scenario. In particular, for 3D massive MIMO communications, it is
widely recognized that there is a unique challenge in exploiting the inherent
3D spatial structure inside the channel coefficients [15, 16, 17], and thus
needs tailored algorithm designs.

To reveal the underlying 3D spatial structure, an emerging trend is to
leverage the angular channel model, which has been validated by real-world
measurements [8, 9]. In this model, the channel coefficient is modelled as the
summation of different propagation paths, each of which is specified by angle
parameters and fading parameters. Since this channel model mimics the signal
propagations in the physical world, its parameters are with clear
interpretations. On the other hand, the mathematical form of the angular
channel model shares a lot of similarities with the models in array signal
processing, and thus has triggered tremendous research progress on massive
MIMO channel estimation from an array signal processing perspective [18]. In
particular, many fundamental ideas in array signal processing, including
discrete Fourier transformation (DFT) [19], multiple signal classification
(MUSIC) [20], and estimation of signal parameters via rotational invariance
technique (ESPRIT) [21], have tapped into the algorithm design of massive MIMO
channel estimation and brought significant performance improvement [22, 23,
24]. Furthermore, these inspiring ideas have integrated with advanced tensor
methods to achieve more accurate 3D massive MIMO channel estimation even with
limited pilot signals [27, 17, 28, 25, 26].

However, previous works on tensor-aided 3D massive MIMO communications [17,
28, 27, 25, 26] mainly investigated single-user channel estimation. By
equivalently formulating channel estimation problems as standard tensor
decompositions, a vast number of off-the-shelf tensor decomposition tools can
be utilized [29]. One might contemplate the straightforward extension of
existing works to the scenario where the BS estimates multiple users’ channels
simultaneously. Unfortunately, the coupling among different users’ channels
make the problem formulation deviate from widely-used tensor decompositions.
Instead, it gives a non-standard tensor decomposition format that has not been
well tackled. This unique challenge requires a novel tensor-aided multi-user
channel estimation algorithm design for 3D massive MIMO communications.

The most straightforward approach is to fit the new tensor decomposition model
to the observation data via solving an optimization problem. In particular,
with the widely adopted least-squares (LS) model fitting criterion, it can be
shown that the optimization problem enjoys a block multi-convex property [30],
in the sense that although the original problem is not convex, after fixing
other variables other than one variable, the remaining problem is convex. It
motivates the leverage of block coordinate descent (BCD) method [31] to solve
the model fitting problem. This approach can be interpreted as a maximum-
likelihood (ML) approach under the assumption that the signals are corrupted
by additive white Gaussian noises (AWGNs) [14]. From the viewpoint of machine
learning, it is well known that the ML solution is prone to the overfitting of
noises if the model complexity is not set correctly [14]. In the angular
channel model, the model complexity is determined by the number of independent
propagation paths, which however is unknown in practice [8, 9]. To mitigate
the overfitting, a typical method is to introduce an additive regularization
term that penalizes complicated channel models [14]. However, for the best
channel estimation performance, this approach requires tuning the
regularization parameters carefully to balance the data fitting and the model
complexity control, which inevitably consumes enormous computation resources.
Therefore, in this paper, we aim at answering the following question: _could
we develop a tuning-free channel estimation algorithm that can automatically
learn the optimal channel model complexity from the wireless data?_

This question invites a data-driven approach to the wireless research, in
order to let the wireless data tell its desired channel model complexity. This
goal just coincides with the fundamental philosophy of Bayesian methods[32].
In particular, the Bayesian Occam Razor principle states that the multiple
integrations in the Bayes rule will automatically drive the inferred model to
the simplest one that can still explain the data well. This has enabled
tuning-free algorithm designs for Bayesian neural network [33], sparse
Bayesian learning [34], and more recently Bayesian structured tensor
decompositions [35, 36]. Its great success in automatic model complexity
control inspires us to rethink the multi-user 3D massive MIMO channel
estimation problem from a Bayesian perspective. In particular, by establishing
the probabilistic model and designing the efficient inference algorithm, in
this paper, we propose a novel tuning-free tensor-aided multi-user channel
estimation algorithm for 3D massive MIMO communications. Numerical results
have corroborated its excellent performance in terms of both channel
estimation accuracy and overfitting avoidance.

The remainder of this paper is organized as follows. In Section II, after
introducing the system model, the multi-user channel estimation problem is
formulated as a non-standard tensor decomposition problem. To fit the new
tensor model to the wireless data, a BCD-based method is briefly introduced in
Section III, which however is prone to the overfitting of noises. To avoid the
overfitting via a tuning-free approach, a novel algorithm based on Bayesian
modelling and inference is proposed in Section IV. Simulation results are
presented in Section V to show the effectiveness of the proposed algorithm.
Finally, conclusions are drawn in Section VI.

Notation: Boldface lowercase and uppercase letters will be used for vectors
and matrices, respectively. Tensors are written as calligraphic letters.
$\mathbb{E}[\leavevmode\nobreak\ \cdot\leavevmode\nobreak\ ]$ denotes the
expectation of its argument and $j\triangleq\sqrt{-1}$. Superscripts $T$, $*$
and $H$ denote transpose, conjugate and Hermitian, respectively. $\bm{A}^{-1}$
denotes the inverse of a matrix $\bm{A}$. The operator
$\textrm{Tr}\left({\bm{A}}\right)$ denotes the trace of a matrix $\bm{A}$.
$\bigparallel\cdot{\bigparallel}_{F}$ represents the Frobenius norm of the
argument. $\mathcal{CN}(\bm{x}|\bm{u},\bm{R})$ stands for the probability
density function of a circularly-symmetric complex Gaussian vector $\bm{x}$
with mean $\bm{u}$ and covariance matrix $\bm{R}$. The operator
$\mathfrak{Re}\\{\cdot\\}$ represents the real part of the argument. The
symbol $\propto$ represents a linear scalar relationship between two real-
valued functions. The $N\times N$ diagonal matrix with diagonal elements
$a_{1}$ through $a_{N}$ is represented as
$\mathrm{diag}\\{a_{1},a_{2},...,a_{N}\\}$, while $\bm{I}_{M}$ represents the
$M\times M$ identity matrix. The $(i,j)^{th}$ element, the $i^{th}$ row, and
the $j^{th}$ column of a matrix $\bm{A}$ are represented by $\bm{A}_{i,j}$,
$\bm{A}_{i,:}$ and $\bm{A}_{:,j}$, respectively.

##  II System Model And Problem Formulation: When Angular Channel Model Meets
Multi-user Massive MIMO

Consider a massive MIMO system where the BS is equipped with a 3D uniform
cuboid antenna array (UCA), as shown in Figure 1, and each user is equipped
with a single antenna. Let $M$ and $N$ denote the number of antennas at the BS
and the number of users, respectively. In the BS, with the first antenna
assumed to be the origin of the coordinate system, the number of antennas in
the x-direction, y-direction and z-direction are $I_{1}$, $I_{2}$ and $I_{3}$,
respectively (i.e., $M=I_{1}I_{2}I_{3}$). Obviously, the UCA includes the
uniform rectangular array (URA) and the uniform linear array (ULA) as its
special cases by setting some of $\\{I_{1},I_{2},I_{3}\\}$ to be one.

In this paper, we consider the uplink transmission where all the users
simultaneously transmit their pilot signals to the BS through narrow-band non-
line-of-sight (NLOS) channels111The discussions on incorporating the LOS path
are presented in _Remark 1_ (at the end of Section IV. A). Each user is
assigned a unique pilot sequence $\bm{s}_{n}=[s_{n}(1),...,s_{n}(L)]^{T}$ with
length $L$, which is assumed to be smaller than the channel coherence length.
The channel state information (CSI) from the $n^{th}$ user to the $m^{th}$
antenna at the BS is modeled as a complex coefficient $h_{m}^{n}$. Then, the
received discrete-time complex baseband signal at the BS can be modeled as

$\displaystyle\bm{Y}$
$\displaystyle=\sum_{n=1}^{N}\bm{s}_{n}\bm{h}_{n}^{T}+\bm{W}=\bm{S}\bm{H}+\bm{W},$
(1)

where vector $\bm{h}_{n}=[h_{1}^{n},h_{2}^{n},...,h_{M}^{n}]^{T}$ collects
channel coefficients for the $n^{th}$ user, and each element $w_{l,m}$ in the
noise matrix $\bm{W}$ denotes the additive white Gaussian noise (AWGN) at the
BS, i.e., $w_{l,m}\sim\mathcal{CN}(0,\beta^{-1})$ is spatially and temporally
independent. Pilot matrix $\bm{S}\in\mathbb{C}^{L\times N}$ is with the
$n^{th}$ column being $\bm{s}_{n}$, and channel matrix
$\bm{H}\in\mathbb{C}^{N\times M}$ is with the $n^{th}$ row being
$\bm{h}_{n}^{T}$.

The goal of multi-user channel estimation is to estimate the channel matrix
$\bm{H}$ from the received data $\bm{Y}$ at the BS with the help of the pilot
matrix $\bm{S}$. From data model (1), a standard least-squares (LS) solution
can be obtained immediately:

$\displaystyle\hat{\bm{H}}^{\mathrm{LS}}=(\bm{S}^{H}\bm{S})^{-1}\bm{S}^{H}\bm{Y}.$
(2)

When using the LS estimator, since no prior information is incorporated, it is
well known that the estimation accuracy heavily relies on the pilot length $L$
[39]. That is, to ensure accurate channel estimation, long pilot sequences are
required to be transmitted at user sides, which however will consume
invaluable spectral resources. This is not desirable in practical massive MIMO
systems, and thus calls for alternative solutions that can significantly
improve the accuracy of channel estimation even with limited pilots.

Figure 1: A massive MIMO system where the base station (BS) is equipped with a
three dimensional (3D) uniform cuboid antenna array (UCA). Under both spatial
and frequency narrow-band assumption, the relative delays among different non-
line-of-sight (NLOS) propagation paths and the effect of different subcarriers
are assumed to be negligible [46].

To achieve this, recent research works have repeatedly shown the glimmers of
hope from channel model structure exploitation. In particular, a vast amount
of research works [18, 12, 22, 23, 24, 17, 28, 27, 25, 26] have shown the
effectiveness of the angular channel model, which has been validated by real-
world measurements [8, 9]. Not only does it depict the signal propagations in
the physical word via angle parameters and fading parameters, it also bridges
the design of massive MIMO systems and array signal processing techniques.
More specifically, it assumes that for the $n^{th}$ user, the channel model
consists of $R^{n}$ propagation paths, each of which is determined by path
gain $\xi_{r^{n}}$, elevation angle $\theta_{r^{n}}$ and azimuth angle
$\phi_{r^{n}}$, i.e.,

$\displaystyle
h_{x_{m},y_{m},z_{m}}^{n}=\sum_{r=1}^{R^{n}}\xi_{r^{n}}\mathrm{exp}\Big{\\{}j\frac{2\pi}{\lambda_{c}}\big{[}x_{m}\sin\theta_{r^{n}}\cos\phi_{r^{n}}$
$\displaystyle+y_{m}\sin\theta_{r^{n}}\sin\phi_{r^{n}}+z_{m}\cos\theta_{r^{n}}\big{]}\Big{\\}},$
(3)

where $\lambda_{c}$ is the wavelength of the carrier signal and
$(x_{m},y_{m},z_{m})$ is the coordinate of the $m^{th}$ antenna. Notice that
in (3), under both spatial and frequency narrow-band assumption, the relative
delays among different propagation paths and the effect of different
subcarriers are assumed to be negligible [46]. Although the channel model (3)
shares a lot of similarities with the array signal processing model [20],
[21], there is a slight difference. Due to the block-fading assumption, the
path gain $\xi_{r^{n}}$ is assumed to be unchanged during the channel
estimation. In contrast, in most array signal processing applications [20],
[21], the source signals are assumed to be time-varying.

Using the angular channel model (3), instead of directly estimating the
channel matrix $\bm{H}$ with $MN$ unknown parameters, one could estimate the
model parameters
$\\{\\{\xi_{r^{n}},\theta_{r^{n}},\phi_{r^{n}}\\}_{r=1}^{R^{n}}\\}_{n=1}^{N}$
and then reconstruct the channel coefficients. By this approach, only
$\sum_{n=1}^{N}3R^{n}$ unknown parameters need to be estimated. Since the path
number $R^{n}$ is usually much smaller than the antenna number $M$, the
adoption of the angular channel model significantly reduces the number of
unknown parameters, and thus allows more accurate channel estimation.

However, estimating these unknown parameters
$\\{\\{\xi_{r^{n}},\theta_{r^{n}},\phi_{r^{n}}\\}_{r=1}^{R^{n}}\\}_{n=1}^{N}$
from the observation data $\bm{Y}$ is quite challenging, since they are
nonlinearly coupled in the channel model (3). In particular, motivated by the
AWGN assumption, the following LS-based optimization problem can be
formulated:

$\displaystyle\min_{\\{\\{\xi_{r^{n}},\theta_{r^{n}},\phi_{r^{n}}\\}_{r=1}^{R^{n}}\\}_{n=1}^{N}}\sum_{l=1}^{L}\sum_{m=1}^{M}\bigparallel\bm{Y}_{l,m}-\sum_{n=1}^{N}s_{n}(l)$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\times\sum_{r=1}^{R^{n}}\xi_{r^{n}}\mathrm{exp}\Big{\\{}j\frac{2\pi}{\lambda_{c}}\big{[}x_{m}\sin\theta_{r^{n}}\cos\phi_{r^{n}}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
+y_{m}\sin\theta_{r^{n}}\sin\phi_{r^{n}}+z_{m}\cos\theta_{r^{n}}\big{]}\Big{\\}}{\bigparallel}_{F}^{2}.$
(4)

Similar optimization problems have been investigated in array signal
processing society for many decades [19], and it is widely agreed that
directly optimizing these variables
$\\{\\{\xi_{r^{n}},\theta_{r^{n}},\phi_{r^{n}}\\}_{r=1}^{R^{n}}\\}_{n=1}^{N}$
is prohibitively expensive in computations. Instead, subspace methods (e.g,
MUSIC [20] and ESPRIT [21]) have come up as the main tools to enable the
accurate estimation of these unknown parameters in a computationally efficient
manner. Its key idea is to recast the parameter estimation problem as the low
dimensional signal subspace recovery problem, for which an array of
dimensionality reduction tools (e.g., low-rank matrix decompositions [14]) are
off-the-shelf. Inspired by this idea and further exploiting the 3D structure
of the antenna array at the BS, recent studies leverage low-rank tensor
decompositions to achieve better signal subspace recovery and subsequently
more accurate channel estimation [27, 17, 28, 25, 26]. However, these works
are limited to the single-user case, while their extensions to the multi-user
scenario are not straightforward.

To see this, following the tensor modelling in previous works [27, 17, 28, 25,
26], we re-organize the channel coefficients
$\\{h_{x_{m},y_{m},z_{m}}^{n}\\}_{m=1}^{M}$ into a 3D tensor
$\mathcal{H}^{n}\in\mathbb{C}^{I_{1}\times I_{2}\times I_{3}}$. In particular,
let set $S_{x}$ collect all the antennas’ x-axis coordinates
$\\{x_{m}\\}_{m}^{M}$, with repeated values eliminated and remaining values
sorted via the ascending order, (i.e., $S_{x}(i_{1})$ is the $i_{1}^{th}$
largest number in $\\{x_{m}\\}_{m=1}^{M}$). Similarly, let $S_{y}$ and $S_{z}$
collect ordered coordinates on the y-axis and z-axis, respectively. Then, the
channel coefficient $h_{x_{m},y_{m},z_{m}}^{n}$ can be equivalently re-indexed
as $h_{i_{1},i_{2},i_{3}}^{n}$, where the index $(i_{1},i_{3},i_{3})$
satisfies $S_{x}(i_{1})=x_{m}$, $S_{y}(i_{2})=y_{m}$ and $S_{y}(i_{3})=z_{m}$.
Using the new indexing scheme, the channel model (3) can be equivalently re-
expressed as

$\displaystyle
h_{i_{1},i_{2},i_{3}}^{n}=\sum_{r=1}^{R^{n}}\xi_{r^{n}}\mathrm{exp}\Big{\\{}j\frac{2\pi}{\lambda_{c}}\big{[}S_{x}(i_{1})\sin\theta_{r^{n}}\cos\phi_{r^{n}}$
$\displaystyle+S_{y}(i_{2})\sin\theta_{r^{n}}\sin\phi_{r^{n}}+S_{3}(i_{3})\cos\theta_{r^{n}}\big{]}\Big{\\}},$
$\displaystyle=\sum_{r=1}^{R^{n}}\xi_{r^{n}}\mathrm{exp}(S_{x}(i_{1})u_{r^{n}})\mathrm{exp}(S_{y}(i_{2})v_{r^{n}})\mathrm{exp}(S_{z}(i_{3})p_{r^{n}}),$
(5)

where $u_{r^{n}}=j\frac{2\pi}{\lambda_{c}}\sin\theta_{r^{n}}\cos\phi_{r^{n}}$,
$v_{r^{n}}=j\frac{2\pi}{\lambda_{c}}\sin\theta_{r^{n}}\sin\phi_{r^{n}}$ and
$p_{r^{n}}=j\frac{2\pi}{\lambda_{c}}\cos\theta_{r^{n}}$. Comparing expression
(5) to the definition of tensor canonical polyadic decomposition (CPD) [40],
it is easy to identify that each 3D channel tensor $\mathcal{H}^{n}$ follows a
rank-$R^{n}$ tensor CPD format:

$\displaystyle\mathcal{H}^{n}$
$\displaystyle\triangleq\llbracket\bm{U}^{(n)},\bm{V}^{(n)},\left[\bm{\xi}^{n}\right]^{T}\diamond\bm{P}^{(n)}\rrbracket$
$\displaystyle=\sum_{r=1}^{R^{n}}\bm{u}^{(n)}_{r}\circ\bm{v}^{(n)}_{r}\circ\bm{p}^{(n)}_{r}\xi_{r^{n}},$
(6)

where $\bm{U}^{(n)}\in\mathbb{C}^{I_{1}\times R^{n}}$ is with its
$(i_{1},r)^{th}$ element being $\mathrm{exp}(S_{x}(i_{1})u_{r^{n}})$;
$\bm{V}^{(n)}\in\mathbb{C}^{I_{2}\times R^{n}}$ is with its $(i_{2},r)^{th}$
element being $\mathrm{exp}(S_{y}(i_{2})v_{r^{n}})$; and
$\bm{P}^{(n)}\in\mathbb{C}^{I_{3}\times R^{n}}$ is with its $(i_{3},r)^{th}$
element being $\mathrm{exp}(S_{z}(i_{3})p_{r^{n}})$. $\bm{u}^{(n)}_{r}$,
$\bm{v}^{(n)}_{r}$ and $\bm{p}^{(n)}_{r}$ are the $r^{th}$ columns in matrix
$\bm{U}^{(n)}$, $\bm{V}^{(n)}$ and $\bm{P}^{(n)}$, respectively. Symbol
$\circ$ denotes vector outer product, $\diamond$ denotes Khatri-Rao product,
and vector
$\bm{\xi}^{n}=[\xi_{1^{n}},\xi_{2^{n}},...,\xi_{R^{n}}]^{T}\in\mathbb{C}^{R^{n}\times
1}$. Then, optimization problem (4) can be equivalently formulated as:

$\displaystyle\min_{\\{\\{\xi_{r^{n}},\theta_{r^{n}},\phi_{r^{n}}\\}_{r=1}^{R^{n}}\\}_{n=1}^{N}}\sum_{l=1}^{L}\bigparallel\mathcal{Y}_{l}-\sum_{n=1}^{N}s_{n}(l)$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\times\llbracket\bm{U}^{(n)},\bm{V}^{(n)},\left[\bm{\xi}^{n}\right]^{T}\diamond\bm{P}^{(n)}\rrbracket{\bigparallel}_{F}^{2}.$
(7)

In (7), $\mathcal{Y}_{l}\in\mathbb{C}^{I_{1}\times I_{2}\times I_{3}}$ is a 3D
tensor that collects measurements $\\{\bm{Y}_{l,m}\\}_{m=1}^{M}$ according to
$\left[\mathcal{Y}_{l}\right]_{i_{1},i_{2},i_{3}}=\bm{Y}_{l,m^{*}}$ where
$x_{m^{*}}=S_{x}(i_{1})$, $y_{m^{*}}=S_{y}(i_{2})$ and
$z_{m^{*}}=S_{z}(i_{3})$. At last, inspired by the subspace methods, rather
than searching parameters
$\\{\\{\xi_{r^{n}},\theta_{r^{n}},\phi_{r^{n}}\\}_{r=1}^{R^{n}}\\}_{n=1}^{N}$
exhaustively, it is viable to firstly estimate the factor matrices
$\\{\bm{U}^{(n)},\bm{V}^{(n)},\left[\bm{\xi}^{n}\right]^{T}\diamond\bm{P}^{(n)}\\}_{n=1}^{N}$
from the data $\\{\mathcal{Y}_{l}\\}_{l=1}^{L}$ and then reconstruct each
channel tensor $\mathcal{H}^{n}$ via
$\llbracket\bm{U}^{(n)},\bm{V}^{(n)},\left[\bm{\xi}^{n}\right]^{T}\diamond\bm{P}^{(n)}\rrbracket$.
For the brevity of notations, let factor matrices
$\\{\bm{U}^{(n)},\bm{V}^{(n)},\left[\bm{\xi}^{n}\right]^{T}\diamond\bm{P}^{(n)}\\}_{n=1}^{N}$
simply be denoted by
$\\{\bm{\Xi}^{(1),n},\bm{\Xi}^{(2),n},\bm{\Xi}^{(3),n}\\}_{n=1}^{N}$. Then,
the channel estimation problem can be formulated as:

$\displaystyle\min_{\\{\\{\bm{\Xi}^{(k),n}\\}_{k=1}^{3}\\}_{n=1}^{N}}\sum_{l=1}^{L}\bigparallel\mathcal{Y}_{l}-\sum_{n=1}^{N}s_{n}(l)$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\times\llbracket\bm{\Xi}^{(1),n},\bm{\Xi}^{(2),n},\bm{\Xi}^{(3),n}\rrbracket{\bigparallel}_{F}^{2}.$
(8)

If $N=1$ (i.e. single-user case), problem (8) can be decomposed into a set of
standard tensor CPD problems that enjoy appealing uniqueness property (see
Appendix A), for which there are abundant “out of the box” algorithms [29].
However, when the BS serves multiple users simultaneously, the summand inside
the Frobenius norm prohibits the straightforward utilization of standard
tensor decomposition tools, and thus make the multi-user channel estimation
problem much more challenging than the single-user counterpart [27, 17, 28,
25, 26]. In particular, when $N>1$, the factor matrices
$\\{\\{\bm{\Xi}^{(k),n}\\}_{k=1}^{3}\\}_{n=1}^{N}$ are intricately coupled
together after expanding the Frobenius norm (as elaborated in Appendix B).
This coupling is much more complicated than those appeared in existing single-
user channel estimation works [27, 17, 28, 25, 26], making their extensions
(either using optimizations or Bayesian methods) to the multi-user scenario
not straightforward. This paper makes the first attempt to tackle this
challenge.

##  III Direct Fitting via Block Coordinate Descent: How to Avoid Overfitting?

It is not difficult to show the non-convexity of problem (8), since all the
factor matrices $\\{\\{\bm{\Xi}^{(k),n}\\}_{k=1}^{3}\\}_{n=1}^{N}$ are coupled
together via multi-linear products. However, a closer inspection could reveal
its appealing block multi-convex property [30], based on which BCD method [31]
can be leveraged. More specifically, although problem (8) is not convex with
respect to $\\{\\{\bm{\Xi}^{(k),n}\\}_{k=1}^{3}\\}_{n=1}^{N}$, after fixing
all the variables to their latest updates other than a single factor matrix
$\bm{\Xi}^{(k),n}$, in the iteration $t+1$, the remaining subproblem can be
formulated as:

$\displaystyle\min_{\bm{\Xi}^{(k),n}}\sum_{l=1}^{L}\bigparallel\\!\left[\mathcal{B}_{l,p\neq
n}\right]^{\kappa}(k)\\!\\!-\\!\\!s_{n}(l)\bm{\Xi}^{(k),n}\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\left[\bm{\Xi}^{(j),n}\right]^{\kappa}\right)^{T}\\!\\!{\bigparallel}_{F}^{2},$
(9)

where

$\displaystyle\left[\mathcal{B}_{l,p\neq n}\right]^{\kappa}$
$\displaystyle\triangleq\mathcal{Y}_{l}\\!\\!-\\!\\!\sum_{p=1,p\neq
n}^{N}s_{p}(l)\Big{\llbracket}\\!\left[\bm{\Xi}^{(1),p}\right]^{\kappa},\left[\bm{\Xi}^{(2),p}\right]^{\kappa},\left[\bm{\Xi}^{(3),p}\right]^{\kappa}\\!\Big{\rrbracket},$
(10)

and $\kappa$ denotes the most recent update index, i.e., $\kappa=t+1$ when
$j<k$ or $p<n$, and $\kappa=t$ otherwise. $\left[\mathcal{B}_{l,n\neq
k}\right]^{\kappa}(k)$ is a matrix obtained by unfolding the tensor
$\left[\mathcal{B}_{l,n\neq k}\right]^{\kappa}$ along its $k^{th}$ dimension,
and the multiple Khatri-Rao products $\mathop{\diamond}\limits_{n=1,n\neq
k}^{N}{\bm{A}}^{(n)}={\bm{A}}^{(N)}\diamond{\bm{A}}^{(N-1)}\diamond\cdots\diamond{\bm{A}}^{(k+1)}\diamond{\bm{A}}^{(k-1)}\diamond\cdots\diamond{\bm{A}}^{(1)}$.
After checking the positive semi-definiteness of the Hessian matrix,
subproblem (9) can be shown to be convex. Then, by setting the derivative of
the objective function in (9) to be zero, the closed-form optimal solution can
be obtained as follows:

$\displaystyle\left[\bm{\Xi}^{(k),n}\right]^{t+1}\\!\\!\\!=\\!\left[\sum_{l=1}^{L}\left[\mathcal{B}_{l,n\neq
k}\right]^{\kappa}(k)s_{n}(l)^{*}\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\left[\bm{\Xi}^{(j),n}\right]^{\kappa}\right)^{*}\right]$
$\displaystyle\\!\\!\\!\times\\!\\!\\!\left[\sum_{l=1}^{L}|s_{n}(l)|^{2}\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\left[\bm{\Xi}^{(j),n}\right]^{\kappa}\right)^{T}\\!\\!\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\left[\bm{\Xi}^{(j),n}\right]^{\kappa}\right)^{*}\\!\right]^{-1}.$ (11)

Since each subproblem is convex, after iteratively updating each
$\left[\bm{\Xi}^{(k),n}\right]^{t+1}$ via (11), the resultant BCD algorithm,
which is summarized in Algorithm 1 at the top of this page, is guaranteed to
converge to a critical point of the objective function of (8) [31].

Algorithm 1: BCD Based Multi-user Channel Estimation

Initializations: Choose path number estimates $\\{\hat{R}^{n}\\}_{n=1}^{N}$
and initial values
$\\{\\{\left[\bm{\Xi}^{(k),n}\right]^{0}\\}_{k=1}^{3}\\}_{n=1}^{N}$.

Iterations: For the iteration $t+1$ ($t\geq 0$),

Update factor matrix $\left[\bm{\Xi}^{(k),n}\right]^{t+1}$

$\displaystyle\left[\bm{\Xi}^{(k),n}\right]^{t+1}\\!=\\!\left[\sum_{l=1}^{L}\left[\mathcal{B}_{l,n\neq
k}\right]^{\kappa}(k)s_{n}(l)^{*}\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\left[\bm{\Xi}^{(j),n}\right]^{\kappa}\right)^{*}\right]$
$\displaystyle\\!\\!\\!\times\\!\\!\\!\left[\sum_{l=1}^{L}|s_{n}(l)|^{2}\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\left[\bm{\Xi}^{(j),n}\right]^{\kappa}\right)^{T}\\!\\!\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\left[\bm{\Xi}^{(j),n}\right]^{\kappa}\right)^{*}\\!\right]^{-1},$

where $\left[\mathcal{B}_{l,n\neq k}\right]^{\kappa}$ is computed using (10);
$\kappa$ denotes the most recent update index, i.e., $\kappa=t+1$ when $j<k$
or $p<n$, and $\kappa=t$ otherwise.

Until Convergence

Channel Estimation:

$\hat{\mathcal{H}}^{n}=\Big{\llbracket}\left[\bm{\Xi}^{(1),n}\right]^{t+1},\left[\bm{\Xi}^{(2),n}\right]^{t+1},\left[\bm{\Xi}^{(3),n}\right]^{t+1}\Big{\rrbracket},\forall
n.$

However, to implement Algorithm 1, prior knowledge about the path numbers
$\\{R^{n}\\}_{n=1}^{N}$ are required, which however is difficult to acquire in
practice. On the other hand, as seen in (6), path number $R^{n}$ controls the
number of rank-1 component in the CPD model, and thus controls the channel
model complexity. In [40], it has been shown that generally $R^{n}$ is non-
deterministic polynomial-time hard (NP-hard) to estimate. With over-estimated
path numbers $\\{\hat{R}^{n}\\}_{n=1}^{N}$, (or equivalently too complicated
channel models), directly fitting the tensor channel models
$\\{\mathcal{H}\\}_{n=1}^{N}$ to the observation data
$\\{\mathcal{Y}_{l}\\}_{l=1}^{N}$ via Algorithm 1 will be prone to the
overfitting of noises, and thus will cause performance deterioration in
channel estimation. To avoid the overfitting, a widely-adopted approach is to
introduce an additional regularization term that penalizes the model
complexity as follows [14]:

$\displaystyle\min_{\\{\\{\bm{\Xi}^{(k),n}\\}_{k=1}^{3}\\}_{n=1}^{N}}\sum_{l=1}^{L}\bigparallel\mathcal{Y}_{l}-\sum_{n=1}^{N}s_{n}(l)$
$\displaystyle\times\llbracket\bm{\Xi}^{(1),n},\bm{\Xi}^{(2),n},\bm{\Xi}^{(3),n}\rrbracket{\bigparallel}_{F}^{2}+\sum_{n=1}^{N}\sum_{k=1}^{3}\gamma_{k}^{n}g(\bm{\Xi}^{(k),n}),$
(12)

where the regularization function $g(\cdot)$ (e.g., $l_{1}$ norm and $l_{2}$
norm) is pre-selected. For the best channel estimation performance, the
regularization parameters $\\{\\{\gamma_{k}^{n}\\}_{k=1}^{3}\\}_{n=1}^{N}$
need to be finely tuned, which is however computationally demanding. Then, an
immediate question is: could we develop a tuning-free approach such that the
model complexity can be optimally learnt from the data? This is fundamentally
important in achieving overfitting avoidance for channel estimation.

##  IV Towards A tuning-free Approach: A Bayesian Perspective

This question has been partially answered in the research of Bayesian
modelling and inference. In the early pioneering works of Mackay [33] and
Tipping [34] on Bayesian neural network and relevance vector machine,
sparsity-enhancing priors were employed to encode an over-parameterized model.
Together with the Bayesian Occam Razor principle, which indicates that
Bayesian inference will automatically seek the simplest model that can still
explain the data adequately, the inference algorithm will drive redundant
model parameters to be all zeros and thus effectively control the model
complexity. This idea has triggered flourishing research on Bayesian
compressive sensing [37], sparse Bayesian learning [38], and more recently
Bayesian structured tensor decompositions [35, 36]. However, for the tensor-
aided multi-user channel estimation problem in (8), since it does not follow a
standard tensor decomposition format, there is no existing Bayesian solution.
Therefore, in this paper, we develop such an algorithm from the first
principle of Bayesian methods.

###  IV-A Sparsity-promoting Probabilistic Modelling

Firstly, the probabilistic model, which encodes the knowledge of problem (8),
needs to be established. Motivated by the LS cost function in (8) (and
equivalently the AWGN assumption in data model (1)), a Gaussian likelihood
function is adopted as follows:

$\displaystyle
p(\\{\mathcal{Y}_{l}\\}_{l=1}^{L}|\\{\\{\bm{\Xi}^{(k),n}\\}_{k=1}^{3}\\}_{n=1}^{N})$
$\displaystyle\propto\exp\Bigg{\\{}-\beta\sum_{l=1}^{L}\bigparallel\mathcal{Y}_{l}-\sum_{n=1}^{N}s_{n}(l)\llbracket\bm{\Xi}^{(1),n},\bm{\Xi}^{(2),n},\bm{\Xi}^{(3),n}\rrbracket{\bigparallel}_{F}^{2}\Bigg{\\}},$
(13)

where $\beta^{-1}$ is the noise power. To reflect the non-informativeness of
the noise power, gamma distribution
$p(\beta)=\mathrm{gamma}(\beta|\epsilon,\epsilon)$ with $\epsilon$ being very
small (e.g., $10^{-6}$) is employed as its prior.

For factor matrices $\\{\\{\bm{\Xi}^{(k),n}\\}_{k=1}^{3}\\}_{n=1}^{N}$, since
$\\{\bm{\Xi}^{(k),n}\\}_{k=1}^{3}$ determines the $n^{th}$ user’s channel
tensor $\mathcal{H}^{n}$ and the channels of different users are assumed to be
statistically independent, we have
$p(\\{\\{\bm{\Xi}^{(k),n}\\}_{k=1}^{3}\\}_{n=1}^{N})=\prod_{n=1}^{N}p(\\{\bm{\Xi}^{(k),n}\\}_{k=1}^{3})$.
In the channel model (6), it can be observed that the channel tensor
$\mathcal{H}^{n}$ is the summation of $R^{n}$ rank-1 tensors, each of which is
determined by the $r^{th}$ columns in the three factor matrices, i.e.,
$\\{\bm{\Xi}^{(k),n}_{:,r}\\}_{k=1}^{3}$. By treating each
$\\{\bm{\Xi}^{(k),n}_{:,r}\\}_{k=1}^{3}$ as an independent building block of
the channel model, we have
$p(\\{\bm{\Xi}^{(k),n}\\}_{k=1}^{3})=\prod_{k=1}^{3}\prod_{r=1}^{R^{n}}p(\bm{\Xi}^{(k),n}_{:,r})$.
Since the exact path number $R^{n}$ is unknown, an upper bound on its value
$\bar{R}^{n}$ is assumed to give an over-parameterized model. Then, inspired
by previous Bayesian sparsity modelling works [33, 34], a sparsity-promoting
Gaussian-gamma prior distribution is utilized to model
$\\{\bm{\Xi}^{(k),n}_{:,r}\\}_{r=1}^{\bar{R}^{n}}$. Finally, the sparsity-
promoting prior for all the factor matrices is:

$\displaystyle
p(\\{\\{\bm{\Xi}^{(k),n}\\}_{k=1}^{3}\\}_{n=1}^{N}|\\{\\{\gamma_{r}^{n}\\}_{r=1}^{\bar{R}^{n}}\\}_{n=1}^{N})$
$\displaystyle=\prod_{n=1}^{N}p(\\{\bm{\Xi}^{(k),n}\\}_{k=1}^{3}|\\{\gamma_{r}^{n}\\}_{r=1}^{\bar{R}^{n}})$
$\displaystyle=\prod_{n=1}^{N}\prod_{k=1}^{3}\prod_{r=1}^{\bar{R}^{n}}p(\bm{\Xi}^{(k),n}_{:,r}|\gamma_{r}^{n})$
$\displaystyle=\prod_{n=1}^{N}\prod_{k=1}^{3}\prod_{r=1}^{\bar{R}^{n}}\mathcal{CN}\left(\bm{\Xi}^{(k),n}_{:,r}|\bm{0}_{I_{k}},(\gamma_{r}^{n})^{-1}\bm{I}_{I_{k}}\right),$
(14) $\displaystyle
p(\\{\\{\gamma_{r}^{n}\\}_{r=1}^{\bar{R}^{n}}\\}_{n=1}^{N})=\prod_{n=1}^{N}\prod_{r=1}^{\bar{R}^{n}}\mathrm{gamma}(\gamma_{r}^{n}|\epsilon,\epsilon),$
(15)

where $\epsilon$ is a very small number (e.g., $10^{-6}$) that indicates the
non-informativeness of the prior model. Consequently, the proposed
probabilistic model is a three-layer Bayes network, as illustrated in Figure
2.

_Remark 1:_ In (14), due to the NLOS assumption adopted in this paper, there
is no need to explicitly model the significant power differences of paths. If
the LOS path is considered (with its power much larger than other paths),
order statistics [51] might be exploited to model this structural information,
which is an interesting future direction to investigate.

Figure 2: Probabilistic model for 3D massive MIMO multi-user channel
estimation.

###  IV-B Variational Inference: Block Coordinate Descent Over Functional
Space

Table I: Optimal variational probability density functions Variational pdfs | Remarks
---|---
$Q^{\dagger}\left(\bm{\Xi}^{(k),n}\right)=\mathcal{CMN}(\bm{\Xi}^{(k),n}|\bm{M}^{(k),n},\bm{I}_{I_{k}},\bm{\Sigma}^{(k),n}),\forall n,k$ | Circularly-symmetric complex matrix normal distribution
| with mean $\bm{M}^{(k),n}$ and covariance matrix $\bm{\Sigma}^{(k),n}$
$Q^{\dagger}\left(\gamma_{r}^{n}\right)=\mathrm{gamma}(\gamma_{r}^{n}|a_{r}^{n},b_{r}^{n}),\forall n,r$ | Gamma distribution with shape $a_{r}^{n}$ and rate $b_{r}^{n}$
$Q^{\dagger}\left(\beta\right)=\mathrm{gamma}(\beta|c,d)$ | Gamma distribution with shape $c$ and rate $d$

Let $\bm{\Theta}$ be the set that collects all the unknown variables in the
probabilistic model, i.e,
$\bm{\Theta}=\\{\\{\\{\bm{\Xi}^{(k),n}\\}_{k=1}^{3}\\}_{n=1}^{N},\\{\\{\gamma_{r}^{n}\\}_{r=1}^{\bar{R}^{n}}\\}_{n=1}^{N},\beta\\}$.
The goal of Bayesian inference is to infer the posterior distribution of each
unknown variable. Following the Bayes rule, the posterior distribution
$p(\bm{\Theta}_{i}|\\{\mathcal{Y}_{l}\\}_{l=1}^{L})=\int\frac{p(\bm{\Theta},\\{\mathcal{Y}_{l}\\}_{l=1}^{L})}{\int
p(\bm{\Theta},\\{\mathcal{Y}_{l}\\}_{l=1}^{L})d\bm{\Theta}}d\bm{\Theta}_{j\neq
i}$, where $\bm{\Theta}_{i}$ is part of $\bm{\Theta}$ with
$\cup_{i=1}^{I}\bm{\Theta}_{i}=\bm{\Theta}$ and
$\cap_{i=1}^{I}\bm{\Theta}_{i}=\O$. However, the intricacy of the
probabilistic model does not allow the tractable solution of the multiple
integrations involved [14]. Fortunately, this challenge is not totally new,
and is common for modern Bayesian inference tasks such as Bayesian deep
learning [41] and Bayesian tensor methods [35, 36]. In these works,
variational inference (VI) is advocated since it scales well to complicated
models with a large number of parameters [43, 44]. In essence, VI recasts the
intractable multiple integration problem into a functional optimization
problem. In particular, it solves the following problem [44]:

$\displaystyle\min_{Q(\bm{\Theta})}\mathrm{KL}\big{(}Q\left(\bm{\Theta}\right)\parallel
p\left(\bm{\Theta}\mid\\{\mathcal{Y}_{l}\\}_{l=1}^{L}\right)\big{)}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\triangleq-\mathbb{E}_{Q\left(\bm{\Theta}\right)}\left\\{\ln\frac{p\left(\bm{\Theta}\mid\\{\mathcal{Y}_{l}\\}_{l=1}^{L}\right)}{Q\left(\bm{\Theta}\right)}\right\\}$
$\displaystyle\mathrm{s.t.}\leavevmode\nobreak\ \leavevmode\nobreak\
Q(\bm{\Theta})\in\mathcal{F},$ (16)

where $\mathrm{KL}(\cdot||\cdot)$ denotes the Kullback-Leibler (KL) divergence
between two arguments, and $\mathcal{F}$ is a pre-selected family of
probability density functions (pdfs). The rationale behind problem (16) is:
although the exact posterior distribution
$p\left(\bm{\Theta}\mid\\{\mathcal{Y}_{l}\\}_{l=1}^{L}\right)$ has no closed-
form, we can still seek a tractable variational probability distribution
$Q(\bm{\Theta})$ in one family $\mathcal{F}$ that is the closest to the true
posterior distribution
$p\left(\bm{\Theta}\mid\\{\mathcal{Y}_{l}\\}_{l=1}^{L}\right)$ in terms of the
KL divergence.

The choice of probability distribution family $\mathcal{F}$ is an art, since
it needs to be flexible enough to ensure the freedoms of variational pdfs,
while simple enough to enable efficient functional optimization. Mean-field
family is such a good choice, as evidenced by a lot of Bayesian inference
works [42, 14]. It assumes that
$Q(\bm{\Theta})=\prod_{i=1}^{I}Q(\bm{\Theta}_{i})$. Then, inspired by its
factorized structure, the idea of BCD could be migrated to the functional
space. More specifically, after fixing other variational pdfs
$\\{Q(\bm{\Theta}_{j})\\}_{j\neq i}$ to their latest update results, the pdf
$Q(\bm{\Theta}_{i})$ can be updated via solving the following problem:

$\displaystyle\min_{Q(\bm{\Theta}_{i})}\int
Q(\bm{\Theta}_{i})\Big{(}-\mathbb{E}_{\prod_{j\neq
i}Q(\bm{\Theta}_{j})}\left[\ln
p(\bm{\Theta},\\{\mathcal{Y}_{l}\\}_{l=1}^{L})\right]$ $\displaystyle+\ln
Q(\bm{\Theta}_{i})\Big{)}d{\bm{\Theta}_{i}}.$ (17)

Using variational calculus, the optimal solution of subproblem (17) can be
shown to be [43, 44]:

$\displaystyle
Q^{\dagger}\left(\bm{\Theta}_{i}\right)=\frac{\exp\left(\mathbb{E}_{\prod_{j\neq
i}Q\left(\bm{\Theta}_{j}\right)}\left[\ln{p\left(\bm{\Theta},\\{\mathcal{Y}_{l}\\}_{l=1}^{L}\right)}\right]\right)}{\int\exp\left(\mathbb{E}_{\prod_{j\neq
i}Q\left(\bm{\Theta}_{j}\right)}\left[\ln
p\left(\bm{\Theta},\\{\mathcal{Y}_{l}\\}_{l=1}^{L}\right)\right]\right)d\bm{\Theta}_{i}},$
(18)

where

$\displaystyle\ln
p\left(\bm{\Theta},\\{\mathcal{Y}_{l}\\}_{l=1}^{L}\right)=\left(\prod_{k=1}^{3}I_{k}L\right)\ln\beta-\beta\sum_{l=1}^{L}\bigparallel\mathcal{Y}_{l}-\sum_{n=1}^{N}s_{n}(l)$
$\displaystyle\times\llbracket\bm{\Xi}^{(1),n},\bm{\Xi}^{(2),n},\bm{\Xi}^{(3),n}\rrbracket{\bigparallel}_{F}^{2}+(\epsilon-1)\ln\beta-\epsilon\beta$
$\displaystyle-\sum_{n=1}^{N}\sum_{k=1}^{3}\left[\mathrm{Tr}\left(\bm{\Xi}^{(k),n}\bm{\Gamma}^{n}\left[\bm{\Xi}^{(k),n}\right]^{H}\right)+I_{k}\sum_{r=1}^{\bar{R}^{n}}\ln\gamma_{r}^{n}\right]$
$\displaystyle+\sum_{n=1}^{N}\sum_{r=1}^{\bar{R}^{n}}\left[(\epsilon-1)\ln\gamma_{r}^{n}-\epsilon\gamma_{r}^{n}\right],$
(19)

and
$\bm{\Gamma}^{n}=\mathrm{diag}\\{\gamma_{1}^{n},...,\gamma_{\bar{R}^{n}}^{n}\\}$.

###  IV-C Tuning-free Algorithm Derivation

After substituting (19) into (18), the optimal solutions
$\\{Q^{\dagger}\left(\bm{\Theta}_{i}\right)\\}_{i=1}^{I}$ can be obtained.
Although straightforward as it may seem, multiple integrations involved in
(18) and complicated tensor algebras in (19) both make the derivations
technically challenging and tedious. On the other hand, the coupling among
different users’ channel parameters deviates the algorithm derivations from
those developed in related works on single-user channel estimation [27, 17,
28, 25, 26]. Consequently, it needs much effort to derive the optimal
variational pdfs for the multi-user channel estimation problem. To keep the
brevity of the main body, we move the lengthy derivations to Appendix C and
only present the final optimal solutions
$\\{Q^{\dagger}\left(\bm{\Theta}_{i}\right)\\}_{i=1}^{I}$ in Table I at the
top of this page.

In Table I, the optimal variational pdf for each factor matrix
$Q(\bm{\Xi}^{(k),n})$ is a circularly-symmetric complex matrix normal
distribution [45], where the covariance matrix

$\displaystyle\bm{\Sigma}^{(k),n}=$
$\displaystyle\Bigg{[}\sum_{l=1}^{L}|s_{n}(l)|^{2}\mathbb{E}\left[\beta\right]\mathbb{E}\Big{[}\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\bm{\Xi}^{(j),n}\right)^{T}$
$\displaystyle\times\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\bm{\Xi}^{(j),n}\right)^{*}\Big{]}+\mathbb{E}\left[\bm{\Gamma}^{n}\right]\Bigg{]}^{-1},$
(20)

and mean matrix

$\displaystyle\bm{M}^{(k),n}$
$\displaystyle=\sum_{l=1}^{L}\mathfrak{B}_{l,p\neq
n}(k)s_{n}(l)^{*}\mathbb{E}\left[\beta\right]\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\mathbb{E}\left[\bm{\Xi}^{(j),n}\right]\right)^{*}\bm{\Sigma}^{(k),n},$
(21)

with

$\displaystyle\mathfrak{B}_{l,p\neq n}$
$\displaystyle\triangleq\mathcal{Y}_{l}\\!-\\!\sum_{p=1,p\neq
n}^{N}s_{p}(l)\Big{\llbracket}\mathbb{E}\left[\bm{\Xi}^{(1),p}\right],\mathbb{E}\left[\bm{\Xi}^{(2),p}\right],\mathbb{E}\left[\bm{\Xi}^{(3),p}\right]\Big{\rrbracket}.$
(22)

On the other hand, the optimal variational distributions for each
$\gamma_{r}^{n}$ and $\beta$ are gamma distributions, with parameters

$\displaystyle a_{r}^{n}=\epsilon+\sum_{k=1}^{3}I_{k},$ (23) $\displaystyle
b_{r}^{n}=\epsilon+\sum_{k=1}^{3}\mathbb{E}\left[\left[\bm{\Xi}^{(k),n}_{:,r}\right]^{H}\bm{\Xi}^{(k),n}_{:,r}\right],$
(24) $\displaystyle c=\epsilon+\prod_{k=1}^{3}I_{k}L,$ (25) $\displaystyle
d=\epsilon+\sum_{l=1}^{L}\mathbb{E}\left[\bigparallel\mathcal{Y}_{l}-\sum_{n=1}^{N}s_{n}(l)\llbracket\bm{\Xi}^{(1),n},\bm{\Xi}^{(2),n},\bm{\Xi}^{(3),n}\rrbracket{\bigparallel}_{F}^{2}\right].$
(26)

In (20)-(26), there are several expectations that need to be computed. Some of
them can be directly obtained from their parameters. In particular,
$\mathbb{E}[\bm{\Xi}^{(k),n}]=\bm{M}^{(k),n}$,
$\mathbb{E}\left[\gamma_{r}^{n}\right]=\frac{a_{r}^{n}}{b_{r}^{n}}$ and
$\mathbb{E}\left[\beta\right]=\frac{c}{d}$. Some of them have already computed
in previous works [35, 17]:

$\displaystyle\mathbb{E}\left[\left[\bm{\Xi}^{(k),n}_{:,r}\right]^{H}\bm{\Xi}^{(k),n}_{:,r}\right]=\left[\bm{M}^{(k),n}_{:,r}\right]^{H}\bm{M}^{(k),n}_{:,r}+I_{k}\bm{\Sigma}^{(k),n}_{r,r},$
$\displaystyle\mathbb{E}\Big{[}\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\bm{\Xi}^{(j),n}\right)^{T}\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\bm{\Xi}^{(j),n}\right)^{*}\Big{]}$ (27)
$\displaystyle=\mathop{\odot}\limits_{j=1,j\neq
k}^{3}\left[\left[\bm{M}^{(j),n}\right]^{T}\left[\bm{M}^{(j),n}\right]^{*}+I_{j}\left[\bm{\Sigma}^{(j),n}\right]^{*}\right],$
(28)

where the multiple Hadamard products $\mathop{\odot}\limits_{n=1,n\neq
k}^{N}{\bm{A}}^{(n)}={\bm{A}}^{(N)}\odot{\bm{A}}^{(N-1)}\odot\cdots\diamond{\bm{A}}^{(k+1)}\odot{\bm{A}}^{(k-1)}\odot\cdots\odot{\bm{A}}^{(1)}$.
However, due to the coupling among different users’ channel parameters, there
is one complicated expectation in (26) that has not been tackled so far. In
Appendix D, we show that

$\displaystyle\mathbb{E}\Bigg{[}\bigparallel\mathcal{Y}_{l}-\sum_{n=1}^{N}s_{n}(l)\llbracket\bm{\Xi}^{(1),n},\bm{\Xi}^{(2),n},\bm{\Xi}^{(3),n}\rrbracket{\bigparallel}_{F}^{2}\Bigg{]},$
$\displaystyle=$
$\displaystyle\bigparallel\mathcal{Y}_{l}{\bigparallel}_{F}^{2}-\mathrm{Tr}\Bigg{(}2\mathfrak{Re}\Big{(}\mathcal{Y}_{l}(1)\sum_{n=1}^{N}s_{n}(l)^{*}\left(\mathop{\diamond}\limits_{j=2}^{3}\bm{M}^{(j),n}\right)^{*}$
$\displaystyle\times\left[\bm{M}^{(1),n}\right]^{H}\Big{)}-\sum_{n=1}^{N}\sum_{p=1,p\neq
n}^{N}s_{n}(l)s_{p}(l)^{*}\bm{M}^{(1),n}$
$\displaystyle\times\left(\mathop{\diamond}\limits_{j=2}^{3}\bm{M}^{(j),n}\right)^{T}\left(\mathop{\diamond}\limits_{j=2}^{3}\bm{M}^{(j),p}\right)^{*}\left[\bm{M}^{(1),p}\right]^{H}\Bigg{)}$
$\displaystyle+\mathrm{Tr}\Bigg{(}\sum_{n=1}^{N}|s_{n}(l)|^{2}\left[\left[\bm{M}^{(1),n}\right]^{H}\bm{M}^{(1),n}+I_{1}\bm{\Sigma}^{(1),n}\right]$
$\displaystyle\mathop{\odot}\limits_{k=2}^{3}\left[\left[\bm{M}^{(k),n}\right]^{H}\bm{M}^{(k),n}+I_{k}\bm{\Sigma}^{(k),n}\right]^{*}\Bigg{)}.$
(29)

From (20)-(29), it is easy to see that the parameters of each optimal
variational pdf $Q^{\dagger}\left(\bm{\Theta}_{i}\right)$ rely on the
statistics of other variational pdfs
$\\{Q^{\dagger}\left(\bm{\Theta}_{i}\right)\\}_{j\neq i}$. By alternatively
updating these variational pdfs, a tuning-free iterative channel estimation
algorithm can be summarized in Algorithm 2 at the top of this page.

Algorithm 2: VI Based Multi-user Channel Estimation

Initializations: Choose $\bar{R}^{n}>R^{n},\forall n$, initial values
$\\{\\{\left[\bm{M}^{(k),n}\right]^{0},\left[\bm{\Sigma}^{(k),n}\right]^{0}\\}_{k=1}^{3}\\}_{n=1}^{N}$
and $\epsilon$. Let
$\\{[a_{r}^{n}]^{0},[b_{r}^{n}]^{0},c^{0},d^{0}\\}=\epsilon,\forall r,n$.

Iterations:

For the iteration $t+1$ ($t\geq 0$),

Update the parameters of $Q(\bm{\Xi}^{(k),n})^{t+1}$:

$\displaystyle\left[\bm{\Sigma}^{(k),n}\right]^{t+1}\\!\\!=\\!\\!\Bigg{[}\\!\sum_{l=1}^{L}|s_{n}(l)|^{2}\frac{c^{t}}{d^{t}}\\!\\!\mathop{\odot}\limits_{j=1,j\neq
k}^{3}\\!\\!\Bigg{[}\left[\bm{M}^{(j),n}\right]^{\kappa,T}\\!\\!\left[\bm{M}^{(j),n}\right]^{\kappa,*}$
$\displaystyle+I_{j}\left[\bm{\Sigma}^{(j),n}\right]^{\kappa,*}\Bigg{]}+\mathrm{diag}\left\\{\frac{[a^{n}_{1}]^{t}}{[b^{n}_{1}]^{t}},...,\frac{[a^{n}_{\bar{R}^{n}}]^{t}}{[b^{n}_{\bar{R}^{n}}]^{t}}\right\\}\Bigg{]}^{-1},$
(30)

$\displaystyle\left[\bm{M}^{(k),n}\right]^{t+1}\\!\\!=\\!\\!\sum_{l=1}^{L}\Bigg{(}\mathcal{Y}_{l}\\!-\\!\\!\\!\sum_{p=1,p\neq
n}^{N}s_{p}(l)\Big{\llbracket}\left[\bm{M}^{(1),p}\right]^{\kappa},\left[\bm{M}^{(2),p}\right]^{\kappa}\\!\\!\\!\\!,$
$\displaystyle\left[\bm{M}^{(3),p}\right]^{\kappa}\Big{\rrbracket}\Bigg{)}(k)s_{n}(l)^{*}\frac{c^{t}}{d^{t}}\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\left[\bm{M}^{(j),n}\right]^{\kappa}\right)^{*}\left[\bm{\Sigma}^{(k),n}\right]^{t+1}\\!\\!\\!\\!\\!\\!,$
(31)

where $\kappa$ denotes the most recent update index, i.e., $\kappa=t+1$ when
$j<k$ or $p<n$, and $\kappa=t$ otherwise.

Update the parameters of $Q(\gamma_{r}^{n})^{t+1}$:

$\displaystyle\left[a_{r}^{n}\right]^{t+1}=\epsilon+\sum_{k=1}^{3}I_{k},$ (32)
$\displaystyle\left[b_{r}^{n}\right]^{t+1}=\epsilon+\sum_{k=1}^{3}\left[\bm{M}^{(k),n}_{:,r}\right]^{t+1,H}\left[\bm{M}^{(k),n}_{:,r}\right]^{t+1}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ +I_{k}\bm{\left
missing}[\bm{\Sigma}^{(k),n}_{r,r}\right]^{t+1}.$ (33)

Update the parameters of $Q(\beta)^{t+1}$:

$\displaystyle c^{t+1}=\epsilon+\prod_{k=1}^{3}I_{k}L,$ (34) $\displaystyle
d^{t+1}=\epsilon+\sum_{l=1}^{L}\mathfrak{f}_{l}^{t+1},$ (35)

where $\mathfrak{f}_{l}^{t+1}$ is computed using (29) with
$\\{\bm{M}^{(k),n},\bm{\Sigma}^{(k),n}\\}$ being replaced by
$\\{\left[\bm{M}^{(k),n}\right]^{t+1},\left[\bm{\Sigma}^{(k),n}\right]^{t+1}\\},\forall
n,k$.

Until Convergence

Channel Estimation:

$\hat{\mathcal{H}}^{n}=\Big{\llbracket}\left[\bm{M}^{(1),n}\right]^{t+1},\left[\bm{M}^{(2),n}\right]^{t+1},\left[\bm{M}^{(3),n}\right]^{t+1}\Big{\rrbracket},\forall
n.$

###  IV-D Intuitive Interpretation of Updating Equations

####  IV-D1 Intuitive interpretation of (20) and (21)

The covariance matrix $\bm{\Sigma}^{(k),n}$ of the approximate posterior
distribution $Q(\bm{\Xi}^{(k),n})$ computed in (20) combines the prior
information from $\mathbb{E}\left[\bm{\Gamma}^{n}\right]$ and the information
from other factor matrices. It is then used as the rotation matrix in the
estimation of the factor matrix mean $\bm{M}^{(k),n}$ in (21), which takes the
linear combination of the observation data and other factor matrices. If there
is no prior information $\mathbb{E}\left[\bm{\Gamma}^{n}\right]$ and no noise
precision estimate $\mathbb{E}\left[\beta\right]$, the update equation (21) is
very similar to the BCD update in (11), since VI essentially performs BCD
steps over the functional space.

####  IV-D2 Intuitive interpretation of (23)-(26)

From (23) and (24), it can be seen that
$\mathbb{E}\left[\gamma_{r}^{n}\right]=\frac{a_{r}^{n}}{b_{r}^{n}}$ is
proportional to the inverse of the sum of the $r^{th}$ column powers in all
three factor matrices. Therefore, if the $r^{th}$ columns are learnt to be
nearly zero, it will give a very large
$\mathbb{E}\left[\gamma_{r}^{n}\right]=\frac{a_{r}^{n}}{b_{r}^{n}}$, which
will further encourage the sparsity of the $r^{th}$ columns in (21). On the
other hand, it is straightforward to see that (25) is related to the number of
observations and (26) approximates the model fitting error.

###  IV-E Further Discussions and Insights

To gain more insights from the proposed algorithm, discussions on its
automatic model complexity control, convergence property, and computational
complexity are presented in this subsection.

####  IV-E1 Automatic model complexity control

In Algorithm 2, although the initial channel model is over-parameterized,
there is no need to manually tune any parameter to control the model
complexity for overfitting avoidance, since the parameters
$\\{\\{\frac{c_{r}^{n}}{d_{r}^{n}}\\}_{r=1}^{\bar{R}^{n}}\\}_{n=1}^{N}$, which
are the expectations of
$\\{\\{\gamma_{r}^{n}\\}_{r=1}^{\bar{R}^{n}}\\}_{n=1}^{N}$, effectively shrink
the values of redundant columns in the factor matrices. In particular, if
$\frac{c_{r}^{n}}{d_{r}^{n}}$ is learnt to be very large, they would
contribute to the covariance matrix of the factor matrix (as seen in (30)) and
then rescale the $r^{th}$ column of the factor matrix to approach zero values
(as seen in (31)). On the other hand, parameters
$\\{\\{\frac{c_{r}^{n}}{d_{r}^{n}}\\}_{r=1}^{\bar{R}^{n}}\\}_{n=1}^{N}$ will
be updated together with other parameters in the algorithm, following the
principle of the employed Bayesian framework.

####  IV-E2 Convergence property

The algorithm is developed under the framework of mean-field VI, which
inherently performs BCD steps over the functional space. Its convergence
result has been established in [43]. In particular, it has been shown that
when the variational pdf is optimized using (18) in each iteration (just as
what we have done in this paper), the limit point generated by the BCD steps
over the functional space of variational pdfs is guaranteed to be at least a
stationary point of the KL divergence in (16) under the assumption of mean-
field family [43] .

####  IV-E3 Computational complexity

In each iteration, the computational complexity of Algorithm 2 is dominated by
the steps of updating the factor matrices, costing
$O(\sum_{n=1}^{N}\prod_{k=1}^{3}3I_{k}(\bar{R}^{n})^{2}+\sum_{n=1}^{N}\sum_{k=1}^{3}(\bar{R}^{n})^{3})$.
The overall complexity is about
$O(q(\sum_{n=1}^{N}\prod_{k=1}^{3}3I_{k}(\bar{R}^{n})^{2}+\sum_{n=1}^{N}\sum_{k=1}^{3}(\bar{R}^{n})^{3}))$
where $q$ is the number of iterations required for convergence. Thus it can be
seen that the complexity is comparable to that of Algorithm 1, in which the
computational complexity is
$O(q^{\prime}(\sum_{n=1}^{N}\prod_{k=1}^{3}3I_{k}(\hat{R}^{n})^{2}+\sum_{n=1}^{N}\sum_{k=1}^{3}(\hat{R}^{n})^{3}))$
where $q^{\prime}$ is the number of iterations at convergence.

##  V Numerical Results and Discussions

In this section, numerical results are presented to assess the channel
estimation performance of the proposed tuning-free algorithm (i.e., Algorithm
2). Consider a UCA with $M=512$ antenna elements, which are deployed in a 3D
grid with dimensions $I_{1}=8$, $I_{2}=8$, $I_{3}=8$ and the inter-grid
spacing $d_{x}=d_{y}=d_{z}=\lambda_{c}/2$. Assume that there are $N=5$ users
simultaneously transmitting signals to the BS. For each user, there are
$R^{n}=3$ propagation paths with elevation angles randomly selected from
$[-\pi/2,\pi/2]$ and azimuth angles randomly selected from $[-\pi,\pi]$. The
pilot length is $L=10$, and each pilot symbol is sampled from a zero-mean
circularly-symmetric complex Gaussian distribution with unit variance. The
path gains $\\{\xi_{r^{n}}\\}_{r,n}$ are drawn from a zero-mean circularly-
symmetric complex Gaussian distribution with unit variance, and without any
correlation across $r$ and $n$. The signal-to-noise ratio (SNR) is defined as
$10\log_{10}\left(\frac{\sum_{l=1}^{L}\sum_{n=1}^{N}s_{n}(l)\bigparallel\llbracket\bm{U}^{(n)},\bm{V}^{(n)},\left[\bm{\xi}^{n}\right]^{T}\diamond\bm{P}^{(n)}\rrbracket{\bigparallel}_{F}^{2}}{\parallel\mathcal{W}\parallel_{F}^{2}}\right)$
where $\mathcal{W}\in\mathbb{C}^{I_{1}\times I_{2}\times I_{3}\times L}$ is a
tensor collecting all the noise samples. For the proposed tuning-free
algorithm, initial mean $\left[\bm{M}^{(k),n}\right]^{0}$ for each matrix
$\bm{\Xi}^{(k),n}$ is drawn from a zero-mean circularly-symmetric complex
matrix normal distribution with an identity covariance matrix, and the initial
covariance matrix is set as
$\left[\bm{\Sigma}^{(k),n}\right]^{0}=\bm{I}_{\bar{R}^{n}\times\bar{R}^{n}}$.
The upper bound for channel path $\bar{R}^{n}=\min\\{I_{1},I_{2},I_{3}\\}=8$
unless stated otherwise, which is a common practice in Bayesian tensor
decompositions [17],[27],[35],[36]. Each point in the following figures is an
average of 100 Monte-Carlo trials.

###  V-A Convergence Property and Automatic Channel Model Complexity Learning

The convergence behavior of the proposed tuning-free algorithm is shown in
Figure 3 under two different SNRs, where the mean-square-error (MSE) of
channel estimation
$\frac{1}{MN}\sum_{n=1}^{N}||\hat{\mathcal{H}}^{n}-\mathcal{H}^{n}||_{F}^{2}$
is adopted as the measure. From Figure 3, it can be seen that the MSEs of the
proposed algorithm decrease significantly in the first tens of iterations and
then gradually converge to stable values.

Figure 3: The convergence behavior of the proposed algorithm under SNR = 10 dB
and SNR = 20 dB ($\bar{R}^{n}$ = 8, $R^{n}=3,L=10$). Figure 4: Performance of
channel estimation versus different channel path upper bound values
$\bar{R}^{n}$ (SNR = 20 dB, $R^{n}=3,L=10$).

To see whether the proposed algorithm is sensitive to the initial upper bound
value $\bar{R}^{n}$, under SNR = 20 dB, the MSEs of the proposed tuning-free
algorithm (labeled as VI-$\bar{R}^{n}$) are presented in Figure 4, in which
the MSEs of the LS method (labeled as LS), the BCD method (i.e., Algorithm 1)
with incorrect path numbers $\\{\bar{R}^{n}\\}_{n=1}^{N}$ (labeled as
BCD-$\bar{R}^{n}$) and the genie-aided BCD method with exact path numbers
$\\{{R}^{n}\\}_{n=1}^{N}$ (labeled as BCD-$R^{n}$) are served as benchmarks.
From Figure 4, it can be seen that the proposed algorithm with different
values of $\bar{R}^{n}\in\\{6,8,10,30,50\\}$ shows indistinguishable channel
estimation performances to those of the genie-aided BCD-$R^{n}$ method. Notice
that $\bar{R}^{n}\in\\{30,50\\}$ is much larger than the true path number
(tensor rank) $R^{n}=3$. This shows that with different values of upper bound
$\bar{R}^{n}$, the proposed tuning-free algorithm still can learn the model
complexity well and then give accurate channel estimation results. On the
other hand, with much larger $\bar{R}^{n}\in\\{30,50\\}$, the
BCD-$\bar{R}^{n}$ algorithm overfits the noises heavily, and even cannot
outperform the LS method in channel estimations.

Figure 5: The estimates of $\\{\gamma_{r}^{n}\\}_{r=1}^{\bar{R}^{n}}$ for user
1 and user 3 in different Monte-Carlo trials ($\bar{R}^{n}=8$, $R^{n}=3$,
$L=10$, SNR = 20 dB). Figure 6: The estimates of
$\\{\gamma_{r}^{n}\\}_{r=1}^{\bar{R}^{n}}$ for user 1 and user 3 in different
Monte-Carlo trials ($\bar{R}^{n}=30$, $R^{n}=3$, $L=10$, SNR = 20 dB).

As discussed in Section IV. E, the estimation results of
$\\{\gamma_{r}^{n}\\}_{r=1}^{\bar{R}^{n}}$ under different ${\bar{R}}^{n}$s
determine the channel model complexity learning performance of the proposed
method. Since $\\{\gamma_{r}^{n}\\}_{r=1}^{\bar{R}^{n}}$ in different Monte-
Carlo trials are possibly with different sparsity patterns (i.e., the very
small values might appear in different subscripts $r$), averaging them over
Monte-Carlo trials is not informative. Therefore, in Figure 5 and Figure 6, we
present the estimation results of $\\{\gamma_{r}^{n}\\}_{r=1}^{\bar{R}^{n}}$
for user 1 and user 3 in three independent Monte-Carlo trials under
${\bar{R}}^{n}=8$ and ${\bar{R}}^{n}=30$ respectively. From these two figures,
it can be seen that under both ${\bar{R}}^{n}=8$ and ${\bar{R}}^{n}=30$, only
three $\gamma_{r}^{n}$s were estimated to be very small, indicating that there
are three significant channel paths for each user. Since the exact path number
$R_{n}=3$, it shows that the proposed method can accurately recover the
channel model complexity and thus avoid the overfitting.

###  V-B Channel Estimation Performance

To assess the channel estimation performance at different SNRs, the MSEs of
different algorithm are shown in Figure 7. From Figure 7, it is obvious that
the three tensor-aided methods (VI-$\bar{R}^{n}$, BCD-$\bar{R}^{n}$, and
BCD-$R^{n}$) achieve much more accurate channel estimation than the LS method,
due to the exploitation of tensor structures in the adopted angular channel
model. It can be also observed that the genie-aided BCD-$R^{n}$ method
achieves the lowest MSE for a wide range of SNRs, since it fits the channel
coefficients into the observation data assuming the accurate channel model
complexity, which however is not available in practice. With a wrong guess of
the path numbers, the MSEs of the BCD-$\bar{R}^{n}$ method are much higher
than those of the BCD-$R^{n}$ algorithm, due to the overfitting of noises. In
contrast, although the proposed VI-$\bar{R}^{n}$ algorithm is also with a
wrong guess of the path numbers, its MSEs are nearly the same as those of the
genie-aided BCD-$R^{n}$ method. This shows the effectiveness of the Bayesian
method in automatic model complexity control and overfitting avoidance.

Figure 7: Performance of channel estimation versus SNRs ($\bar{R}^{n}$ = 8,
$R^{n}=3,L=10$).

On the other hand, we present the running time of the three iterative tensor-
aided channel estimation algorithms (VI-$\bar{R}^{n}$, BCD-$\bar{R}^{n}$, and
BCD-$R^{n}$) in Table II. From Table II, it can be observed that the proposed
algorithm is with comparable running time to that of the BCD-$\bar{R}^{n}$
approach, which validates the complexity analysis in Section IV. E. Notice
that these two approaches cost much more time than the genie-aided BCD-$R^{n}$
algorithm, since they need to update over-determined model parameters.

Table II: Running time (second) of different channel estimation algorithms ($\bar{R}^{n}=8$, ${R}^{n}=3$, $L=10$) SNR | BCD-$R^{n}$ | BCD-$\bar{R}^{n}$ | VI-$\bar{R}^{n}$
---|---|---|---
10 dB | 1.7942 | 6.5550 | 6.6890
20 dB | 1.5422 | 4.1581 | 4.9155

To see how the model complexity of channels affects different algorithms,
under SNR = 20 dB, the MSEs of channel estimations versus different path
numbers in the channel model are presented in Figure 8. In previous simulation
studies, $R^{n}=3,\forall n$ is considered. Here we further consider different
path number values $R^{n}=\\{2,3,4,5,6\\}$, each of which indicates different
channel model complexities. With a higher $R^{n}$, there are more unknown
channel coefficients. Then, it is expected that the channel estimation
performance would degrade given the same amount of the observation data. This
conjecture has been validated by Figure 8, in which the MSEs indeed increase
as $R^{n}$ increases. On the other hand, it can be seen that the proposed
VI-$\bar{R}^{n}$ algorithm achieves indistinguishable performances as those of
the genie-aided BCD-$R^{n}$ method. This shows that the proposed algorithm can
learn a wide range of model complexities and then effectively shrink redundant
channel model parameters for overfitting avoidance.

Figure 8: Performance of channel estimation versus different path numbers (SNR
= 20 dB, $\bar{R}^{n}=8,L=10$). Figure 9: Performance of channel estimation
versus different pilot lengths (SNR = 20 dB, $\bar{R}^{n}=8,R^{n}=3$).

Finally, in Figure 9, we show how the proposed VI-$\bar{R}^{n}$ algorithm
saves the pilot resources for channel estimation, compared to the standard LS
method. From Figure 9, it is clear that given the MSE $10^{-3}$, the proposed
VI-$\bar{R}^{n}$ algorithm needs only around $20$ pilot signals while the LS
method requires about $50$ pilot signals. The gain comes from both the tensor
structure exploitation of channel model and the Bayesian philosophy in
overfitting avoidance.

##  VI Conclusions and Future Research

In this paper, the multi-user channel estimation problem for 3D massive MIMO
communications was investigated through the lens of Bayesian tensor methods.
The channel estimation problem was firstly recasted as a factor matrix
learning problem for a non-standard tensor decomposition model, which requires
a novel learning algorithm design with an integrated feature of overfitting
avoidance. To achieve this goal, a tuning-free channel estimation algorithm
was proposed in this paper under the framework of Bayesian modelling and
inference. Numerical studies have shown the excellent channel estimation
performance of the proposed method in terms of both accuracy and overfitting
avoidance.

In future research, the exploitation of the shift-invariance property [21],
[50] in cubic antenna array at the base station might give a new tensor model
for massive MIMO communications, which also needs novel tuning-free channel
estimation algorithm design. _We believe that the integration of tensor model,
array signal processing and Bayesian method will bring us closer to the era of
“Joint Model-and-Data-Driven Wireless Communications”._

##  Appendix A Uniqueness Property of Tensor CPD

In [40], a sufficient condition for the uniqueness of tensor CPD is stated as
follows.

Uniqueness condition for CPD [40]. _If
$\llbracket\bm{A}^{(1)},\bm{A}^{(2)},...,\bm{A}^{(N)}\rrbracket$ =
$\llbracket\bm{\Xi}^{(1)},\bm{\Xi}^{(2)},...,\bm{\Xi}^{(N)}\rrbracket$, and
$\sum_{n=1}^{N}k_{n}\geq 2L+(N-1)$ where $k_{i}$ denotes the k-rank of matrix
$\bm{A}^{(i)}$ and $L$ is the tensor rank. Then the following equations hold:
$\bm{\Xi}^{(1)}=\bm{A}^{(1)}\bm{\Delta}\bm{\Lambda}^{(1)}$,
$\bm{\Xi}^{(2)}=\bm{A}^{(2)}\bm{\Delta}\bm{\Lambda}^{(2)}$, …,
$\bm{\Xi}^{(N)}=\bm{A}^{(N)}\bm{\Delta}\bm{\Lambda}^{(N)}$ where $\bm{\Delta}$
is a permutation matrix and diagonal matrix $\bm{\Lambda}^{(n)}$ satisfies
$\prod_{n=1}^{N}\bm{\Lambda}^{(n)}=\bm{I}_{L}$ ._

Thus tensor CPD is unique up to trivial scaling and permutation ambiguities.

##  Appendix B Complicated Coupling after Expanding the Frobenius norm in (8)

Since $||\mathcal{X}||_{F}^{2}=||\mathcal{X}(k)||_{F}^{2}$ where
$\mathcal{X}(k)$ is the $k^{th}$ unfolding matrix of the tensor $\mathcal{X}$,
after using the unfolding property of tensor CPD [40], problem (8) is
equivalent to

$\displaystyle\min_{\\{\\{\bm{\Xi}^{(k),n}\\}_{k=1}^{3}\\}_{n=1}^{N}}\sum_{l=1}^{L}\bigparallel\mathcal{Y}_{l}(k)-$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\sum_{n=1}^{N}s_{n}(l)\bm{\Xi}^{(k),n}\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\left[\bm{\Xi}^{(j),n}\right]\right)^{T}{\bigparallel}_{F}^{2}.$ (36)

Further using the fact
$||\mathbf{X}||_{F}^{2}=\mathrm{Tr}(\mathbf{X}\mathbf{X}^{H})$ to expand the
Frobenius norm in (36), we have the following problem:

$\displaystyle\min_{\\{\\{\bm{\Xi}^{(k),n}\\}_{k=1}^{3}\\}_{n=1}^{N}}\sum_{l=1}^{L}\mathrm{Tr}\Bigg{(}\Big{[}\mathcal{Y}_{l}(k)-\sum_{n=1}^{N}s_{n}(l)\bm{\Xi}^{(k),n}$
$\displaystyle\times\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\left[\bm{\Xi}^{(j),n}\right]\right)^{T}\Big{]}\Big{[}\mathcal{Y}_{l}(k)-\sum_{n=1}^{N}s_{n}(l)\bm{\Xi}^{(k),n}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\times\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\left[\bm{\Xi}^{(j),n}\right]\right)^{T}\Big{]}^{H}\Bigg{)}.$ (37)

$\displaystyle\min_{\\{\\{\bm{\Xi}^{(k),n}\\}_{k=1}^{3}\\}_{n=1}^{N}}\sum_{l=1}^{L}\mathrm{Tr}\Bigg{(}\mathcal{Y}_{l}(k)\mathcal{Y}_{l}(k)^{H}-2\mathfrak{Re}\Big{(}\mathcal{Y}_{l}(k)\sum_{n=1}^{N}s_{n}^{*}(l)\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\left[\bm{\Xi}^{(j),n}\right]\right)^{*}\big{[}\bm{\Xi}^{(k),n}\big{]}^{H}\Big{)}$
$\displaystyle+\underbrace{\Big{[}\sum_{n=1}^{N}s_{n}(l)\bm{\Xi}^{(k),n}\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\left[\bm{\Xi}^{(j),n}\right]\right)^{T}\Big{]}\Big{[}\sum_{n=1}^{N}s_{n}(l)\bm{\Xi}^{(k),n}\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\left[\bm{\Xi}^{(j),n}\right]\right)^{T}\Big{]}^{H}}_{\mathfrak{t}}\Bigg{)}.$
(38)

After manipulating algebras, problem (37) becomes (38) at the top of the next
page. In the term $\mathfrak{t}$, it is clear that the product of two
summation terms, i.e.,

$\displaystyle\Big{[}\sum_{n=1}^{N}s_{n}(l)\bm{\Xi}^{(k),n}\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\left[\bm{\Xi}^{(j),n}\right]\right)^{T}\Big{]}$
$\displaystyle\times\Big{[}\sum_{n=1}^{N}s_{n}(l)\bm{\Xi}^{(k),n}\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\left[\bm{\Xi}^{(j),n}\right]\right)^{T}\Big{]}^{H}$ (39)

will result in complicated coupling among the factor matrices
$\\{\\{\bm{\Xi}^{(k),n}\\}_{k=1}^{3}\\}_{n=1}^{N}$. Although problem (38)
seems complicated, if we only optimize a single factor matrix
$\bm{\Xi}^{(k),n}$ while fixing other variables, problem (38) will become
problem (9) in Section III, which is a convex problem and can be easily
solved.

##  Appendix C The Derivations of The Optimal Variational Pdfs in Table I

After substituting (19) into (18) and only keep terms relevant to
$\bm{\Xi}^{(k),n}$, we have

$\displaystyle Q^{\dagger}\left(\bm{\Xi}^{(k),n}\right)$
$\displaystyle\propto\exp\Bigg{\\{}\mathbb{E}\Bigg{[}-\beta\sum_{l=1}^{L}\bigparallel\mathcal{Y}_{l}-s_{n}(l)\llbracket\bm{\Xi}^{(1),n},\bm{\Xi}^{(2),n},$
$\displaystyle\bm{\Xi}^{(3),n}\rrbracket-\sum_{p=1,p\neq
n}^{N}s_{p}(l)\llbracket\bm{\Xi}^{(1),p},\bm{\Xi}^{(2),p},\bm{\Xi}^{(3),p}\rrbracket{\bigparallel}_{F}^{2}$
$\displaystyle-\mathrm{Tr}\left(\bm{\Xi}^{(k),n}\bm{\Gamma}^{n}\left[\bm{\Xi}^{(k),n}\right]^{H}\right)\Bigg{]}\Bigg{\\}}.$
(40)

Then, we utilize the result
$\parallel\bm{A}\parallel_{F}^{2}=\mathrm{Tr}(\bm{A}\bm{A}^{H})$ to expand the
Frobenius norm. After a series of algebra manipulations,
$Q^{\dagger}\left(\bm{\Xi}^{(k),n}\right)$ can be organized to be

$\displaystyle Q^{\dagger}\left(\bm{\Xi}^{(k),n}\right)$
$\displaystyle\propto\exp\Bigg{\\{}\mathbb{E}\Bigg{[}\mathrm{Tr}\Bigg{(}-\bm{\Xi}^{(k),n}\Bigg{(}\sum_{l=1}^{L}|s_{n}(l)|^{2}\beta\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\bm{\Xi}^{(j),n}\right)^{T}$
$\displaystyle\times\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\bm{\Xi}^{(j),n}\right)^{*}+\bm{\Gamma}^{n}\Bigg{)}\left[\bm{\Xi}^{(k),n}\right]^{H}+2\mathfrak{Re}\Bigg{(}\bm{\Xi}^{(k),n}$
$\displaystyle\times\sum_{l=1}^{L}s_{n}(l)\beta\left(\mathop{\diamond}\limits_{j=1,j\neq
k}^{3}\bm{\Xi}^{(j),n}\right)^{T}\Bigg{(}\mathcal{Y}_{l}(k)\\!-\\!\sum_{p=1,p\neq
n}^{N}s_{p}(l)$
$\displaystyle\times\Big{\llbracket}\bm{\Xi}^{(1),p},\bm{\Xi}^{(2),p},\bm{\Xi}^{(3),p}\Big{\rrbracket}(k)\Bigg{)}\Bigg{)}\Bigg{)}\Bigg{]}\Bigg{\\}}.$
(41)

After distributing the expectations and comparing the functional form of (41)
to that of circularly-symmetric complex matrix normal distribution [45], it
can be concluded that
$Q^{\dagger}\left(\bm{\Xi}^{(k),n}\right)=\mathcal{CMN}(\bm{\Xi}^{(k),n}|\bm{M}^{(k),n},\bm{I}_{I_{k}},\bm{\Sigma}^{(k),n})$
with its mean $\bm{M}^{(k),n}$ and covariance matrix $\bm{\Sigma}^{(k),n}$
being defined in (20) and (21).

Similarly, after substituting (19) and (18), and fixing all the variables
other than $\\{\\{\gamma_{r}^{n}\\}_{r=1}^{\bar{R}^{n}}\\}_{n=1}^{N}$, we have

$\displaystyle
Q^{\dagger}\left(\\{\\{\gamma_{r}^{n}\\}_{r=1}^{\bar{R}^{n}}\\}_{n=1}^{N}\right)\propto\exp\Bigg{\\{}\mathbb{E}\Bigg{[}\sum_{n=1}^{N}\sum_{k=1}^{3}-\mathrm{Tr}\Big{(}\bm{\Xi}^{(k),n}\bm{\Gamma}^{n}$
$\displaystyle\left[\bm{\Xi}^{(k),n}\right]^{H}\Big{)}+I_{k}\sum_{r=1}^{\bar{R}^{n}}\ln\gamma_{r}^{n}+(\epsilon-1)\ln\gamma_{r}^{n}-\epsilon\gamma_{r}^{n}\Bigg{]}\Bigg{\\}}.$
(42)

Using the fact that
$\mathrm{Tr}\Big{(}\bm{\Xi}^{(k),n}\bm{\Gamma}^{n}\left[\bm{\Xi}^{(k),n}\right]^{H}\Big{)}=\sum_{r=1}^{\bar{R}^{n}}\gamma_{r}^{n}\left[\bm{\Xi}^{(k),n}\right]_{:,r}^{H}\left[\bm{\Xi}^{(k),n}\right]_{:,r}$,
it can be shown that

$\displaystyle
Q^{\dagger}\left(\\{\\{\gamma_{r}^{n}\\}_{r=1}^{\bar{R}^{n}}\\}_{n=1}^{N}\right)\propto\prod_{n=1}^{N}\prod_{r=1}^{\bar{R}^{n}}\exp\Bigg{\\{}\mathbb{E}\Bigg{[}-\gamma_{r}^{n}\sum_{k=1}^{3}\Big{(}\left[\bm{\Xi}^{(k),n}\right]_{:,r}^{H}$
$\displaystyle\left[\bm{\Xi}^{(k),n}\right]_{:,r}+I_{k}\ln\gamma_{r}^{n}\Big{)}+(\epsilon-1)\ln\gamma_{r}^{n}-\epsilon\gamma_{r}^{n}\Bigg{]}\Bigg{\\}}.$
(43)

It is easy to conclude that
$Q^{\dagger}\left(\\{\\{\gamma_{r}^{n}\\}_{r=1}^{\bar{R}^{n}}\\}_{n=1}^{N}\right)=\prod_{n=1}^{N}\prod_{r=1}^{\bar{R}^{n}}Q^{\dagger}(\gamma_{r}^{n})$,
where

$\displaystyle
Q^{\dagger}(\gamma_{r}^{n})\propto\exp\Bigg{\\{}\left(\sum_{k=1}^{3}I_{k}+\epsilon-1\right)\ln\gamma_{r}^{n}$
$\displaystyle-\gamma_{r}^{n}\left(\epsilon+\sum_{k=1}^{3}\mathbb{E}\left[\left[\bm{\Xi}^{(k),n}\right]_{:,r}^{H}\left[\bm{\Xi}^{(k),n}\right]_{:,r}\right]\right)\Bigg{\\}}.$
(44)

By comparing (44) to the functional form of gamma distribution, we have
$Q^{\dagger}\left(\gamma_{r}^{n}\right)=\mathrm{gamma}(\gamma_{r}^{n}|a_{r}^{n},b_{r}^{n})$,
where $a_{r}^{n},b_{r}^{n}$ is defined by (23) and (24) respectively.

Finally, we use (18) and (19) again to derive the optimal variational pdf
$Q^{\dagger}\left(\beta\right)$. It can be shown that

$\displaystyle
Q^{\dagger}\left(\beta\right)\propto\exp\Bigg{\\{}\left(\prod_{k=1}^{3}I_{k}L+\epsilon-1\right)\ln\beta$
$\displaystyle-\beta\Bigg{(}\epsilon+\sum_{l=1}^{L}\mathbb{E}\Bigg{[}\bigparallel\mathcal{Y}_{l}-\sum_{n=1}^{N}s_{n}(l)\llbracket\bm{\Xi}^{(1),n},\bm{\Xi}^{(2),n},\bm{\Xi}^{(3),n}\rrbracket{\bigparallel}_{F}^{2}\Bigg{]}\Bigg{)}\Bigg{\\}}.$
(45)

After comparing (45) to the functional form of gamma distribution, it is easy
to identify $Q^{\dagger}\left(\beta\right)=\mathrm{gamma}(\beta|c,d)$, where
$c$ and $d$ are expressed in (25) and (26) respectively.

##  Appendix D Expectation Computation for (26)

In (26), computing the expectation
$\mathbb{E}\left[\bigparallel\mathcal{Y}_{l}-\sum_{n=1}^{N}s_{n}(l)\llbracket\bm{\Xi}^{(1),n},\bm{\Xi}^{(2),n},\bm{\Xi}^{(3),n}\rrbracket{\bigparallel}_{F}^{2}\right]$
is quite complicated. We use the result
$\parallel\bm{A}\parallel_{F}^{2}=\mathrm{Tr}(\bm{A}\bm{A}^{H})$ and the
tensor unfolding property [40] to expand the Frobenius norm:

$\displaystyle\mathbb{E}\left[\bigparallel\mathcal{Y}_{l}-\sum_{n=1}^{N}s_{n}(l)\llbracket\bm{\Xi}^{(1),n},\bm{\Xi}^{(2),n},\bm{\Xi}^{(3),n}\rrbracket{\bigparallel}_{F}^{2}\right]$
$\displaystyle=\mathbb{E}\Bigg{[}\mathrm{Tr}\Bigg{(}\mathcal{Y}_{l}(1)\mathcal{Y}_{l}(1)^{H}-2\mathfrak{Re}\Bigg{(}\mathcal{Y}_{l}(1)\sum_{n=1}^{N}s_{n}(l)^{*}$
$\displaystyle\left(\mathop{\diamond}\limits_{j=2}^{3}\bm{\Xi}^{(j),n}\right)^{*}\left[\bm{\Xi}^{(1),n}\right]^{H}\Bigg{)}+\mathcal{G}_{l}(1)\mathcal{G}_{l}(1)^{H}\Bigg{)}\Bigg{]},$
(46)

where

$\displaystyle\mathcal{G}_{l}=\sum_{n=1}^{N}s_{n}(l)\llbracket\bm{\Xi}^{(1),n},\bm{\Xi}^{(2),n},\bm{\Xi}^{(3),n}\rrbracket.$
(47)

After distributing the expectations, the most complicated term is
$\mathbb{E}\left[\mathcal{G}_{l}(1)\mathcal{G}_{1}(1)^{H}\right]$. Using the
tensor unfolding property [40] again, we have

$\displaystyle\mathbb{E}\left[\mathcal{G}_{l}(1)\mathcal{G}_{1}(1)^{H}\right]$
$\displaystyle=\mathbb{E}\Bigg{[}\mathrm{Tr}\Bigg{(}\Bigg{[}\sum_{n=1}^{N}s_{n}(l)\bm{\Xi}^{(1),n}\left(\mathop{\diamond}\limits_{j=2}^{3}\bm{\Xi}^{(j),n}\right)^{T}\Bigg{]}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\times\Bigg{[}\sum_{n=1}^{N}s_{n}(l)\bm{\Xi}^{(1),n}\left(\mathop{\diamond}\limits_{j=2}^{3}\bm{\Xi}^{(j),n}\right)^{T}\Bigg{]}^{H}\Bigg{)}\Bigg{]}$
$\displaystyle=\mathrm{Tr}\Bigg{(}\sum_{n=1}^{N}\sum_{p=1}^{N}s_{n}(l)s_{p}(l)^{*}\mathbb{E}\Bigg{[}\bm{\Xi}^{(1),n}\left(\mathop{\diamond}\limits_{j=2}^{3}\bm{\Xi}^{(j),n}\right)^{T}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\times\left(\mathop{\diamond}\limits_{j=2}^{3}\bm{\Xi}^{(j),p}\right)^{*}\left[\bm{\Xi}^{(1),p}\right]^{H}\Bigg{]}\Bigg{)}.$
(48)

Further using the results in (27) and (28), we have

$\displaystyle\mathbb{E}\left[\mathcal{G}_{l}(1)\mathcal{G}_{1}(1)^{H}\right]$
$\displaystyle=\mathrm{Tr}\Bigg{(}\sum_{n=1}^{N}\sum_{p=1,p\neq
n}^{N}s_{n}(l)s_{p}(l)^{*}\bm{M}^{(1),n}\left(\mathop{\diamond}\limits_{j=2}^{3}\bm{M}^{(j),n}\right)^{T}$
$\displaystyle\times\left(\mathop{\diamond}\limits_{j=2}^{3}\bm{M}^{(j),p}\right)^{*}\left[\bm{M}^{(1),p}\right]^{H}\Bigg{)}$
$\displaystyle+\mathrm{Tr}\Bigg{(}\sum_{n=1}^{N}|s_{n}(l)|^{2}\left[\left[\bm{M}^{(1),n}\right]^{H}\bm{M}^{(1),n}+I_{1}\bm{\Sigma}^{(1),n}\right]$
$\displaystyle\mathop{\odot}\limits_{k=2}^{3}\left[\left[\bm{M}^{(k),n}\right]^{H}\bm{M}^{(k),n}+I_{k}\bm{\Sigma}^{(k),n}\right]^{*}\Bigg{)}.$
(49)

After putting (49) into (46), the result of (29) can be obtained.

## References

  * [1] E. G. Larsson, O. Edfors, F. Tufvesson and T. L. Marzetta, “Massive MIMO for next generation wireless systems,” _IEEE Communications Magazine_ , vol. 52, no. 2, pp. 186-195, Feb. 2014.
  * [2] E. Bjornson, E. G. Larsson and T. Marzetta, “Massive MIMO: Ten myths and one critical question,” _IEEE Communications Magazine_ , vol. 54, no. 10, pp. 114-123, Feb. 2016.
  * [3] L. Li, T.-H. Chang and S. Cai, “UAV positioning and power control for two-way wireless relaying,” _IEEE Trans. on Wireless Communications_ , vol. 19, no. 2, pp.1008-1024, Feb. 2020.
  * [4] S. Wang, M. Xia, and Y-C. Wu, “Backscatter data collection with unmanned ground vehicle: mobility management and power allocation,” _IEEE Trans. on Wireless Communications_ , vol. 18, no. 4, pp. 2314-2328, Apr. 2019.
  * [5] X. Li, Z. Liu, N. Qin, and S. Jin, “FFR based joint 3D beamforming interference coordination for multi-cell FD-MIMO downlink transmission systems, _IEEE Trans. on Vehicular Technology_ , vol. 69, no. 3, pp. 3105-3118, Mar. 2020.
  * [6] Y. Huang, Q. Wu, T. Wang, G. Zhou, and R. Zhang, “3D beam tracking for cellular-connected UAV,” _IEEE Wireless Communications Letters_ , vol. 9, no. 5, pp. 736-740, May. 2020.
  * [7] S. M. Razavizadeh, M. Ahn, and I. Lee, “Three-dimensional beamforming: A new enabling technology for 5G wireless networks,” _IEEE Signal Processing Magazine_ , vol. 31, no. 6, pp. 94-101, Nov. 2014.
  * [8] Y. H. Nam, B. L. Ng, K. Sayana, Y. Li, J. Zhang, Y. Kim, and J. Lee, “Full-dimension MIMO (FD-MIMO) for next generation cellular technology,” _IEEE Communications Magazine_ , vol. 51, no. 6, pp. 172-179, 2013.
  * [9] Y. Kim, H. Ji, J. Lee, Y. H. Nam, B. L. Ng, I. Tzanidis, and J. Zhang “Full dimension MIMO (FD-MIMO): The next evolution of MIMO in LTE systems,” _IEEE Wireless Communications_ , vol. 21, no. 2, pp. 26-33, 2014.
  * [10] T.-H. Chang, W. -C. Chiang, Y. -W. Peter Hong, and C. -Y. Chi, “Training sequence design for discriminatory channel estimation in wireless MIMO systems,” _IEEE Trans. on Signal Processing_ , vol. 58, no. 12, pp. 6223-6237, Dec. 2010.
  * [11] C. K. Wen, S. Jin, K. K. Wong, J. C. Chen, and P. Ting, “Channel estimation for massive MIMO using Gaussian-mixture Bayesian learning,” _IEEE Trans. on Wireless Communications_ , vol. 14, no. 3, pp.1356-1368, 2014.
  * [12] C. Qian, X. Fu, and N. D. Sidiropoulos, “Algebraic channel estimation algorithms for FDD massive MIMO systems,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 13, no. 5, pp. 961-973, Jun. 2019.
  * [13] Y. Yang, F. Gao, Z. Zhong, B. Ai, and A. Alkhateeb, “Deep transfer learning based downlink channel prediction for FDD massive MIMO systems,” arXiv preprint arXiv:1912.12265, 2019.
  * [14] K. P. Murphy, _Machine learning: a probabilistic perspective_ , MIT press, 2012.
  * [15] R. Shafin, L. Liu, Y. Li, A. Wang, and J. Zhang, “Joint angle and delay estimation for 3D massive MIMO systems based on parametric channel modelling”, _IEEE Trans. on Wireless Communications_ , vol. 16, no. 8, pp. 5370-5383, Aug. 2017.
  * [16] J. Kaleva, N. J. Myers, A. Tölli, R. W. Heath, and U. Madhow, “Short range 3D MIMO mmWave channel reconstruction via geometry-aided AoA estimation,” _in 2019 IEEE Asilomar Conference on Signals, Systems, and Computers_ , pp. 427-431, 2019.
  * [17] L. Cheng, C. Xing, and Y-C. Wu, “Irregular array manifold aided channel estimation in massive MIMO communications,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 13, no. 5, pp. 974-988, Sep. 2019.
  * [18] F. Gao, Z. Tian, E. G. Larsson, M. Pesavento, and S. Jin, “Introduction to the special issue on array signal processing for angular models in massive MIMO communications,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 13, no. 5, pp. 882-885, Sep. 2019.
  * [19] V. Trees, _Detection, estimation, and modulation theory, optimum array processing_ , John Wiley & Sons, 2004.
  * [20] M. Pesavento, A. B. Gershman, and M. Haardt, “Unitary root-MUSIC with a real-valued eigen-decomposition: A theoretical and experimental performance study,” _IEEE Trans. on Signal Processing_ , vol. 48, no.5, pp. 1306-1314, 2000.
  * [21] R. Richard and T. Kailath, “ESPRIT-estimation of signal parameters via rotational invariance techniques, ” _IEEE Trans. on Acoustics, Speech, and Signal Processing_ , vol. 37, no. 7, pp. 984-995, 1989.
  * [22] D. Fan, F. Gao, G. Wang, Z. Zhong, and A. Nallanathan, “Angle domain signal processing aided channel estimation for indoor 60GHz TDD/FDD massive MIMO systems,” _IEEE Journal on Selected Areas in Communications_ , vol. 35, no. 9, pp. 1948-1961, 2017.
  * [23] Z. Guo, X. Wang and W. Heng, “Millimeter-wave channel estimation based on 2-D beamspace MUSIC Method,” _IEEE Trans. on Wireless Communications_ , vol. 16, no. 8 pp. 5384-5394, 2017.
  * [24] R. Shafin, L. Liu, and J. Zhang, “DoA Estimation and RMSE characterization for 3D massive-MIMO/FD-MIMO OFDM system,” in _2015 IEEE Global Communications Conference (GLOBECOM)_ , pp. 1-6, Dec. 2015.
  * [25] D. C. Araújo, A. L. De Almeida, J. P. Da Costa, and R. T. de Sousa,“Tensor-based channel estimation for massive MIMO-OFDM systems,” _IEEE Access_ , vol. 7, pp. 42133-42147, 2019.
  * [26] F. Wen, and C. Liang, “Improved tensor-MODE based direction-of-arrival estimation for massive MIMO systems,” _IEEE Communications Letters_ , vol. 9, no.12, pp. 2182-2185, 2015.
  * [27] L. Cheng, Y-C. Wu, J. Zhang, and L. Liu, “Subspace identification for DOA estimation in massive / full-dimension MIMO system: bad data mitigation and automatic source enumeration,” _IEEE Trans. on Signal Processing_ , vol. 63, no. 22, pp. 5897-5909, Nov 2015.
  * [28] L. Cheng, Y-C. Wu, S. Ma, J. Zhang and L. Liu, “Channel estimation in full-dimensional massive MIMO system using one training symbol,” _in Proceedings of the IEEE 18th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)_ , Hokkaido, Japan, July 2017.
  * [29] N. D. Sidiropoulos, L. D. Lathauwer, X. Fu, K. Huang, E. E. Papalexakis and C. Faloutsos, “Tensor decomposition for signal processing and machine learning,” _IEEE Trans. on Signal Processing_ , vol. 65, no. 13, pp. 3551-3582, 2017.
  * [30] Y. Xu and W. Yin, “A block coordinate descent method for regularized multi-convex optimization with applications to nonnegative tensor factorization and completion,” _SIAM Journal on Imaging Sciences_ , vol. 6, no. 3, pp. 1758-1789, 2013.
  * [31] H. Shi, S. Tu, Y. Xu and W. Yin, “A primer on coordinate descent algorithms,” 2016, arXiv preprint arXiv:1610.00040.
  * [32] M. J. Beal, _Variational algorithms for approximate Bayesian inference_ , London, University of London, 2003.
  * [33] D. J. MacKay, “Probable networks and plausible predictions-a review of practical Bayesian methods for supervised neural networks,” _Computation in Neural Systems_ , vol. 6, no. 3, pp. 469-505, 1995.
  * [34] M. E. Tipping, “Sparse Bayesian learning and the relevance vector machine,” _Journal of Machine Learning Research_ , vol. 1, pp. 211-244, Jun. 2001.
  * [35] L. Cheng, Y-C. Wu, and H. V. Poor, “Probabilistic tensor canonical polyadic decomposition with orthogonal factors,” _IEEE Trans. on Signal Processing_ , vol. 65, no. 3, pp. 663-676, Feb. 2017.
  * [36] L. Cheng, X. Tong, S. Wang, Y-C. Wu, and H. V. Poor, “Learning nonnegative factors from tensor data: probabilistic modelling and inference algorithm,” _IEEE Trans. on Signal Processing_ , accepted, Feb. 2020.
  * [37] S. Ji, Y. Xue, and L. Carin, “Bayesian compressive sensing,” _IEEE Trans. on Signal Processing_ , vol. 16, no, 56, pp. 2346-56, May. 2008.
  * [38] D. Wipf, and B. Rao, “Sparse Bayesian learning for basis selection,” _IEEE Trans. on Signal processing_ , vol. 52, no. 8, pp. 2153-64, Jun. 2004.
  * [39] S. M. Kay, “Fundamentals of statistical signal processing, volume i: Estimation theory”, PTR Prentice-Hall, Englewood Cliffs, 1993.
  * [40] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” _SIAM Review_ , vol. 51, no. 3, pp. 455-500, Aug. 2009.
  * [41] Q. Liu and D. Wang, “Stein variational gradient descent: A general purpose bayesian inference algorithm,” _in Advances in Neural Information Processing Systems (NeuIPS)_ , pp. 2378-2386, 2016.
  * [42] M. Hoffman, D. Blei, J. Paisley, and C. Wang, “Stochastic variational inference,” _Journal of Machine Learning Research_ , vol. 14, pp. 1303-1347, 2013.
  * [43] M. J. Wainwright and M. I. Jordan, “Graphical models, exponential families, and variational inference,” _Foundations and Trends in Machine Learning_ , vol. 1, no. 102, pp. 1-305, Jan. 2008.
  * [44] C. Zhang, J. Butepage, H. Kjellstrom and S. Mandt, “Advances in variational inference,” _IEEE Trans. on Pattern Analysis and Machine Intelligence_ , vol. 41, no. 8. pp. 2008-2026, Aug. 2019.
  * [45] A. K. Gupta and D. K. Nagar, _Matrix Variate Distributions_ , CRC Press,1999.
  * [46] M. Wang, F. Gao, S. Jin, and H. Lin, “An overview of enhanced massive MIMO with array signal processing techniques,” _IEEE J. Sel. Topics Signal Process._ , vol. 13, no. 5, pp. 886-901, Sep. 2019.
  * [47] X. Gao, L. Dai, S. Zhou, A. M. Sayeed, and L. Hanzo, “Wideband beamspace channel estimation for millimeter-wave MIMO systems relying on lens antenna arrays,” _IEEE Trans. Signal Process._ , vol. 67, no. 18, pp. 4809-4824, Sep. 2019.
  * [48] Z. Ding, L. Dai, and H. V. Poor, “MIMO-NOMA design for small packet transmission in the Internet of things,” _IEEE Access_ , vol. 4, pp. 1393-1405, Apr. 2016
  * [49] L. Liu, E. G. Larsson, W. Yu, P. Popovski, C. Stefanovic, and E. de Carvalho, “Sparse signal processing for grant-free massive connectivity: A future paradigm for random access protocols in the Internet of Things,” _IEEE Signal Process. Mag._ , vol. 35, no. 5, pp. 88-99, Sep. 2018.
  * [50] S. Sahnoun and P. Comon, “Joint source estimation and localization,” _IEEE Trans. Signal Process._ , vol. 63, no. 10, pp. 2485-2495, May 2015.
  * [51] B. C. Arnold, N. Balakrishnan, and N. H. Nagaraja, _A First Course in Order Statistics_ , Society for Industrial and Applied Mathematics, 2008.