# Unsupervised Massive MIMO Channel Estimation with Dual-Path Knowledge-Aware
Auto-Encoders

Zhiheng Guo, Yuanzhang Xiao, , Xiang Chen This work was supported by the
National Key R&D Program of China under Grant 2019YFE0196400. (Corresponding
author: Xiang Chen.)
Zhiheng Guo and Xiang Chen are with the School of Electronics and Information
Technology, Sun Yat-sen University, Guangzhou 510275, China (e-mail:
<EMAIL_ADDRESS>chenxiang@mail.sysu.edu.cn).
Yuanzhang Xiao is with the Hawaii Advanced Wireless Technologies Institute,
University of Hawaii, Honolulu, HI(e-mail: yxiao8@hawaii.edu).

###### Abstract

In this paper, an unsupervised deep learning framework based on dual-path
model-driven variational auto-encoders (VAE) is proposed for angle-of-arrivals
(AoAs) and channel estimation in massive MIMO systems. Specifically designed
for channel estimation, the proposed VAE differs from the original VAE in two
aspects. First, the encoder is a dual-path neural network, where one path uses
the received signal to estimate the path gains and path angles, and another
uses the correlation matrix of the received signal to estimate AoAs. Second,
the decoder has fixed weights that implement the signal propagation model,
instead of learnable parameters. This knowledge-aware decoder forces the
encoder to output meaningful physical parameters of interests (i.e., path
gains, path angles, and AoAs), which cannot be achieved by original VAE.
Rigorous analysis is carried out to characterize the multiple global optima
and local optima of the estimation problem, which motivates the design of the
dual-path encoder. By alternating between the estimation of path gains, path
angles and the estimation of AoAs, the encoder is proved to converge. To
further improve the convergence performance, a low-complexity procedure is
proposed to find good initial points. Numerical results validate theoretical
analysis and demonstrate the performance improvements of our proposed
framework.

###### Index Terms:

Massive MIMO, Angle of Arrival, Channel Estimation, Variational Auto-Encoder,
Unsupervised Learning

##  I Introduction

Massive multiple-input-multiple-output (MIMO) is one of the vital technologies
to address the challenges of explosive data traffic and high quality of
service requirement in fifth generation (5G) and beyond wireless communication
systems. Theoretically, massive MIMO can enhance the capacity of a
communication system with the increment of the antennas [1]. In practice,
accurate channel estimation is key to realize the potential gain of massive
MIMO. However, there are challenges in accurate channel estimation, arising
from sophisticated channel modeling [2], costly channel state information
(CSI) [3], and high computational complexity due to a large number of
antennas.

There are several strands of works for angle-of-arrivals (AoAs) and/or channel
estimation in massive MIMO. The first strand of works are subspace-based
estimation methods, dating back to classic algorithms such as multiple signal
classification (MUSIC) [4] and estimation of signal parameters via rotational
invariant techniques (ESPRIT) [5]. Later development of subspace-based methods
include algorithms with lower complexity [6], methods for expanding the
degrees of freedom [7], methods for special antenna array configurations [8],
and methods for millimeter wave massive MIMO system [9]. However, when applied
to massive MIMO, subspace-based methods may have high computational complexity
due to eigendecomposition of large covariance matrices of the received
signals.

The second strand of works are based on compressed sensing [10, 11]. Under
certain assumptions on the sparsity and/or structures of the received signals,
compressed sensing based methods usually pose the estimation problems as
optimization problems, obtain equivalent semidefinite programs, and solve them
by interior-point methods [12]. Recent works have proposed to use the
alternating direction method of multipliers to lower the computational
complexity [13, 14] and massive MIMO [15, 16, 17]. While compressed sensing
based methods have higher estimation accuracy than classic subspace-based
methods, they also have higher computational complexity since the estimation
involves solving an optimization problem (e.g., semidefinite programming).

The third strand of works use variational Bayesian inference on the received
signals for AoA and channel estimation [18, 19]. Variational Bayesian methods
calculate the posterior distribution of the AoAs and channel gains given the
received signals, and choose the estimates that maximize the posterior
distribution. While this approach yields highly accurate estimation [18] and
can be computationally efficient [19], it may be hard to apply to more general
channel models when the posterior distribution is intractable.

Finally, recent development of machine learning has spurred the applications
of deep learning (DL) to wireless communications [20, 21, 22]. DL based
methods can be categorized into “data-driven” and “model-driven”, a
terminology coined by [23, 24]. Data-driven methods follow the end-to-end
principle of DL, where the estimation is achieved by training neural networks
on a large number of training data (e.g., received signals, covariance
matrices) and training labels [25, 26, 27]. Data-driven DL based channel
estimation methods have utilized various neural network architectures (e.g.,
convolutional neural networks [28, 29], complex-valued networks [30]) and have
been applied to various scenarios (e.g., correlated sources [31], unknown
number of users [32], 2-dimensional estimation [33]) and applications (e.g.,
vehicular networks [34], unmanned aerial vehicle surveillance systems [35]).
In contrast, model-driven methods exploit the established physical models and
properties of signal propagation in wireless channels [36], and combine DL
with traditional signal processing schemes such as MIMO detection [37, 38],
subspace-based methods [39, 40, 41, 42] and compressed sensing based methods
[43, 44]. Compared with other methods, DL based methods shift the majority of
the computation to the training stage and has lower complexity during
deployment.

###  I-A Motivation

With few exceptions [43, 44], most DL based methods are supervised, requiring
accurately labeled data samples for training. In practice, it may be
unrealistic and costly to obtain a sufficiently large amount of labeled data
samples needed by supervised estimation [25, 26, 27]. In this sense,
unsupervised methods are more desirable. One unsupervised learning method that
is particularly suitable for channel estimation is the variational auto-
encoder (VAE). The VAE consists of an encoder and a decoder, where the encoder
outputs the posterior distribution of the latent variables given the
observations [45]. For channel estimation, we can define the latent variables
as the channel parameters of interest and the observations as the received
signals. Then we can get maximum likelihood channel estimation from the
encoder output [46, 47].

However, there are two major hurdles in using standard VAEs for AoA and
channel estimation. The first hurdle is the
uninterpretability/unidentifiability of the latent variables [48].
Specifically, the low-dimensional latent variables learned by unsupervised
methods usually has no physical meaning. Therefore, if we feed the received
signals to an original VAE, the encoder output may not be the parameters we
want to estimate (i.e., AoAs and channel gains). This is because without the
labels, the network can only take the difference between the input (i.e., the
received signals) and its reconstruction based on the latent variables as the
loss function. Such a loss function also results in the second hurdle of
multiple local optima. In other words, there may be multiple sets of channel
parameters that are local or global optima of the loss function [44]. Without
labels, it is hard to know which set of channel parameters is correct.

Our proposed framework in this paper overcomes the above hurdles by
restructuring the canonical VAEs, inspired by the insights obtained from our
rigorous analysis. To break the first hurdle of uninterpretable/unidentifiable
latent variables, we “hardwire” the decoder to implement the signal
propagation model [46, 47], instead of learning the decoder from data as in
canonical VAEs. This knowledge-aware decoder enforces the output of the
encoder to be the parameters to estimate. To break the second hurdle of
multiple local and global optima, we rigorously analyze the loss landscape and
characterize the global and local optima. Our analysis leads us to implement
the encoder as a dual-path neural network [49], where one path estimates the
AoAs from the correlation matrix of the received signals and the other path
estimates the path gains and path angles from the received signals. We update
the weights of both paths in an alternating fashion, which ensures the
convergence of the training. To avoid local optima, a low-cost method is
proposed to find good initial points prior to the training. Ablation study is
performed to demonstrate the necessity of the knowledge-aware decoder, the
proposed alternating update manner for the dual-path neural network, and our
choices of initial points.

###  I-B Contribution

In this paper, we propose an unsupervised DL based AoA and channel estimation
framework in massive MIMO. Compared with existing works, our work falls in the
small category of unsupervised DL based methods [43, 44]. While existing
unsupervised methods numerically illustrate the multiplicity of local optima,
our work provides theoretical analysis and characterization of the global and
local optima. Our work also bears some similarity to variational Bayesian
inference [18, 19]. By adopting the DL approach, our work applies to more
general scenarios in which the posterior distribution of the latent variables
is not tractable. Like other DL based methods, our framework has low
computational complexity during deployment.

As an unsupervised method, our framework achieves similar performance as
supervised learning methods, and outperforms traditional methods (e.g., MUSIC
[4] and its advanced variations [50, 51]) and state-of-the-art unsupervised
learning methods [43].

We summarize the main contributions of this paper below.

  * •

To the best of our knowledge, an unsupervised AoA and channel estimation
framework with rigorous analysis is first proposed for massive MIMO in this
paper. Our rigorous analysis of global optima (Propositions 3 and 4) and in
the regime of large numbers of antennas, analysis of local optima (Proposition
5) characterize the landscape of the loss function, which guides the design of
the proposed network and the training method.

  * •

By restructuring the standard VAE, the proposed framework overcomes two
challenges in unsupervised AoA and channel estimation. It breaks the hurdle of
uninterpretable/unidentifiable latent variables by the knowledge-aware
decoder, and alleviates the problem of multiple (bad) local optima by the
dual-path encoder. It also shows an effective way of incorporating knowledge
about wireless channels into deep learning.

  * •

Informed by our analysis, a two-phase training process is proposed. The first
phase initializes the network with good initial points, and the second phase
trains the dual-path encoder to learn the global optima.

  * •

Numerical simulations are performed to validate our theoretical analysis.
Ablation study is performed to evaluate the contributions of the key
components of the proposed framework in overcoming the aforementioned
challenges of unsupervised learning. Numerical results demonstrate that our
proposed framework achieves almost the same performance as the supervised
learning method and outperforms state-of-the-art unsupervised methods.

The remainder of this paper is organized as follows. Section II introduces the
signal model. Section III describes the proposed AoA and channel estimation
framework. Section IV provides in-depth performance analysis of the proposed
framework. Section V presents numerical results to validate our analysis and
demonstrate the performance improvement of our proposed framework. Finally,
Section VI concludes the paper.

Notations: In this paper, scalars, vectors and matrices are denoted by lower-
case letters, bold lower-case letters and bold capital letters, respectively.
The conjugate, transpose, conjugate transpose, and the Frobenius norm of a
matrix are denoted by $\left(\cdot\right)^{*}$, $\left(\cdot\right)^{T}$,
$\left(\cdot\right)^{H}$, and $\left\|\cdot\right\|_{F}$, respectively. The
Kronecker product is denoted by $\odot$. The imaginary unit is $j=\sqrt{-1}$.
Superscripts $\left(\cdot\right)^{R}$ and $\left(\cdot\right)^{I}$ denote the
real and imaginary parts of a variable. The sets of $m$-by-$n$ real and
complex matrices are $\mathbb{R}^{m\times n}$ and $\mathbb{C}^{m\times n}$.
Finally, the symbol $\mathbb{E}\left[\cdot\right]$ is the expectation
operator.

##  II Signal Propagation Model

Consider a typical massive MIMO uplink system with one base station (BS) and
$K$ users. The base station is equipped with a uniform linear array of $N$
antennas, and each user equipment has one antenna. †† Note that the proposed
framework can be applied to general 2-dimensional array configuration (e.g.,
circular, rectangular); we only use uniform linear arrays as a simple case for
theoretical analysis in this paper.  The users are assumed to be stationary
for $M$ time slots, during which the BS collects $M$ snapshots of the received
signals to estimate the users’ locations. Each user $k$’s location is
specified by the distance to the BS and the AoA
$\theta_{k}\in[-\frac{\pi}{2},\frac{\pi}{2}]$. User $k$’s AoA $\theta_{k}$ is
the impinging direction of its signal to the BS, relative to the broadside of
the antenna array (i.e., the line perpendicular to the antenna array). Given
the AoA $\theta_{k}$, user $k$’s signal experiences different phase shifts
across the antenna array, as characterized by the array response vector

$\displaystyle{\bm{a}}{(\theta_{k})}=\left[1,{e^{-j{\frac{2{\pi}{d}}{\lambda}{\sin{\theta_{k}}}}}},\ldots,{e^{-j{\frac{2{\pi}{d}}{\lambda}}{({N}-1)}{\sin{\theta_{k}}}}}\right]^{T},$
(1)

where $\lambda$ is the wavelength and $d\geq\frac{\lambda}{2}$ is the distance
between adjacent antenna elements. User $k$’s signal is also attenuated due to
channel fading. In each time slot $m$, the vector of channel gains from the
$K$ users to the base station is denoted by $\bm{h}_{m}=\bm{\beta}_{m}\odot
e^{j{\bm{\psi}_{m}}}$, where $\bm{\beta}_{m}\in\mathbb{R}^{K}$ is the vector
of path gains and $\bm{\psi}_{m}\in\mathbb{R}^{K}$ is the vector of path
angles. Following the literature [52, 12, 53], we assume block fading
channels, namely the channel gains of different time slots are independent. We
also assume that the vector of channel gains $\bm{h}_{m}$ is jointly Gaussian
distributed, namely
$\bm{h}_{m}\sim\mathcal{CN}(\bm{\mu}_{\bm{h}_{m}},\bm{\Sigma}_{\bm{h}_{m}})$.
Moreover, the AoAs and channel gains are assumed to be independent of each
other [18, 19].

Based on the propagation model, the received signals during the $M$ snapshots,
denoted by $\bm{Y}\in\mathbb{C}^{N\times M}$, can be written as [52, 12, 18,
19, 44] †† Another widely-used model is
$\bm{Y}=\bm{A}(\bm{\theta})\cdot\left[\bm{h}_{1}\odot\bm{s}_{1},\ldots,\bm{h}_{M}\odot\bm{s}_{M}\right]+\bm{N}$,
where $\bm{s}_{m}\in\mathbb{C}^{K}$ is the transmit signal in time slot $m$
[4, 5, 40]. This model is equivalent to ours when the signals are known (e.g.,
pilot and reference signals).

$\bm{Y}=\bm{A}(\bm{\theta})\cdot\left[\bm{h}_{1},\ldots,\bm{h}_{M}\right]+\bm{N},$
(2)

where $\bm{\theta}=[\theta_{1},\ldots,\theta_{K}]^{T}$ is the collection of
users’ AoAs,
$\bm{A}(\bm{\theta})=\left[\bm{a}\left(\theta_{1}\right),\bm{a}\left(\theta_{2}\right),\ldots,\bm{a}\left(\theta_{K}\right)\right]\in\mathbb{C}^{N\times
K}$ is the array response matrix, and
$\bm{N}=\left[\bm{n}_{1},\ldots,\bm{n}_{M}\right]\in\mathbb{C}^{{N}{\times}{M}}$
is the complex circular symmetric white Gaussian noise with
$\bm{n}_{m}\sim\mathcal{CN}(\bm{0},\sigma^{2}\bm{I}_{N})$ for $m=1,\ldots,M$.

The aim of this paper is to estimate the AoAs, $\theta_{k}$ for
$k=1,\ldots,K$, and the path gains $\bm{\beta}_{m}$ and path angles
$\bm{\psi}_{m}$ for $m=1,\ldots,M$, based on the received signals $\bm{Y}$.

Note that we allow the channel gains of different users to be correlated
(i.e., the covariance $\bm{\Sigma}_{\bm{h}_{m}}$ may not be a diagonal
matrix). This is an extension to the independent assumption made by existing
works [4, 5, 52, 12, 40, 44].

##  III Proposed Solution

This section describes our proposed unsupervised AoA and channel estimation
framework for massive MIMO. An overview of our framework and high-level design
principles are provided first, which is followed by detailed description.

###  III-A Overview and Design Principles

A canonical approach of unsupervised learning is to use variational auto-
encoders. The VAE aims to learn the latent variables $\bm{z}$ from the data
(e.g., the received signal $\bm{Y}$ in our case). The original VAE consists of
an encoder that learns the latent variables $\bm{z}$ from the data $\bm{Y}$
and an decoder that reconstructs the data. However, the original VAE has no
control over the physical meanings of the latent variables $\bm{z}$, and
therefore cannot be directly applied to our problem.

One way to control the encoder output is to incorporate the knowledge of
wireless channels into the decoder [46, 47, 54]. Inspired by this idea, we
propose a redesigned VAE architecture for our AoA and channel estimation
problem, where the decoder is fixed and implements the signal propagation
model in (2), instead of having learnable parameters. Specifically, the latent
variable $\bm{z}$ is $K(1+M)\times 1$ column vector, namely

$\bm{z}=\left[\bm{z}_{0}^{T},\bm{z}_{1}^{T},\ldots,\bm{z}_{M}^{T}\right]^{T},$
(3)

where $\bm{z}_{0}\in\mathbb{R}^{K}$ and $\bm{z}_{m}\in\mathbb{C}^{K}$ for
$m=1,\ldots,M$. Then the decoder reconstructs the received signals based on
the latent variable using the following equation:

$\bm{\widehat{Y}}=\bm{A}(\bm{z_{0}})\cdot\left[\bm{z}_{1},\ldots,\bm{z}_{M}\right].$
(4)

Comparing (4) with the signal propagation model (2), it can be seen that
$\bm{z}_{1},\ldots,\bm{z}_{M}$ and $\bm{z}_{0}$ play the roles of the channel
gains and the AoAs. In this way, the knowledge-aware decoder forces the latent
variables to be the parameters to estimate, namely the channel gains and the
AoAs. For this reason, we also write $\bm{z}_{m}=\bm{\widehat{h}}_{m}$ and
$\bm{z}_{0}=\bm{\widehat{\theta}}$, where $\bm{\widehat{h}}_{m}$ and
$\bm{\widehat{\theta}}$ explicitly denote the estimated channel gains and
AoAs.

Similar to standard VAEs, the encoder of our restructured VAE outputs the
posterior distribution of the latent variable given the received signal,
denoted by $p(\bm{z}|\bm{Y})$. In some special cases, the posterior
distribution $p(\bm{z}|\bm{Y})$ can be expressed analytically [18, 19].
However, in the more general scenarios considered in this paper (e.g.,
correlated channel gains), it is hard to derive the analytical expression of
the posterior distribution $p(\bm{z}|\bm{Y})$. Therefore, the encoder outputs
an approximate posterior distribution, denoted by $q(\bm{z}|\bm{Y})$. The goal
is then to minimize the difference between the approximate posterior
distribution $q(\bm{z}|\bm{Y})$ and the true posterior $p(\bm{z}|\bm{Y})$,
measured by the Kullback–Leibler (KL) divergence $\mathcal{D}_{KL}(q||p)$. As
in existing works [18, 19, 55], we restrict to approximate posterior
distributions that can be factorized as

$q\left(\bm{z}|\bm{Y}\right)=q_{0}\left(\bm{\widehat{\theta}}|\bm{Y}\right)\cdot\prod_{m=1}^{M}q_{m}\left(\bm{\widehat{h}}_{m}|\bm{Y}\right).$
(5)

The following proposition characterizes the optimal approximate posterior
distribution that minimizes the KL divergence.

Figure 1: Illustration of the proposed VAE with a dual-path encoder and a
knowledge-aware decoder.

###### Proposition 1

For each time slot $m=1,\ldots,M$, the optimal posterior distribution
$q_{m}^{*}(\bm{\widehat{h}}_{m}|\bm{Y})$ of the channel gain is a circularly
symmetric Gaussian distribution, namely

$q_{m}^{*}\left(\bm{\widehat{h}}_{m}|\bm{Y}\right)=\frac{1}{\pi^{K}|\bm{\Sigma}_{\bm{\widehat{h}}_{m}}|}e^{-(\bm{\widehat{h}}_{m}-\bm{\mu}_{\bm{\widehat{h}}_{m}})^{H}\bm{\Sigma}_{\bm{\widehat{h}}_{m}}^{-1}(\bm{\widehat{h}}_{m}-\bm{\mu}_{\bm{\widehat{h}}_{m}})}.$
(6)

###### Proof:

See Appendix A. ∎

Proposition 1 shows that the optimal posterior distribution of channel gains
is the circularly symmetric Gaussian distribution, which can be completely
determined by the mean vector and the covariance matrix. Therefore, the
encoder only needs to output the estimated mean
$\bm{\mu}_{\bm{\widehat{h}}_{m}}$ and covariance
$\bm{\Sigma}_{\bm{\widehat{h}}_{m}}$ of the channel gains. For the AoAs, a
commonly-used prior distribution of the AoAs is the uniform distribution,
under which it is difficult to calculate the true posterior or the optimal
approximate posterior distribution. Thus, we let the encoder output a point
estimate $\bm{\widehat{\theta}}$ of the AoAs. This is equivalent to using the
Dirac distribution $\delta(\bm{\theta}-\bm{\widehat{\theta}})$ as the
approximate posterior distribution of the AoAs.

Based on the optimal posterior distribution in Proposition 1 and the
knowledge-aware decoder in (4), we can derive the loss function of the
proposed framework. Although the loss function is initially defined as the KL
divergence, it is well known that minimizing the KL divergence is equivalent
to maximizing the evidence lower bound (ELBO) [45]. Hence, we define the loss
function as the negative of the ELBO and calculated in the following
proposition.

###### Proposition 2

In the proposed framework with the knowledge-aware decoder, the loss function,
defined as the negative of the ELBO, can be written analytically as in (7).

$\displaystyle\mathcal{L}^{(train)}\\!$ $\displaystyle=$
$\displaystyle\\!\underbrace{\sum_{m=1}^{M}\left\\{-\log|\bm{\Sigma}_{\bm{\widehat{h}}_{m}}|+\mathrm{tr}\left(\bm{\Sigma}_{\bm{\widehat{h}}_{m}}^{-1}\bm{\Sigma}_{\bm{h}_{m}}\right)+(\bm{\mu}_{\bm{\widehat{h}}_{m}}-\bm{\mu}_{\bm{h}_{m}})^{H}\bm{\Sigma}_{\bm{h}_{m}}^{-1}(\bm{\mu}_{\bm{\widehat{h}}_{m}}-\bm{\mu}_{\bm{h}_{m}})\right\\}}_{\text{KL
divergence between the approximate posterior distribution and the prior
distribution}}$ (7)
$\displaystyle\\!+\underbrace{\mathbb{E}_{q\left(\bm{\hat{\theta}},\bm{\widehat{h}}_{1},\ldots,\bm{\widehat{h}}_{M}|\bm{Y}\right)}\left[\frac{1}{\sigma^{2}}\left\|\bm{Y}-\widehat{\bm{Y}}\right\|_{F}^{2}\right]}_{\text{mean
squared reconstruction error}},$

###### Proof:

See Appendix B. ∎

The ELBO (7) consists of two parts. The first part is the difference of the
entropy of the prior distribution and that of the approximate posterior
distribution. The second part is the mean square error (MSE) between the
received signals and the reconstructed signals. When using the ELBO as the
loss function during training, the restructured VAE aims to reconstruct the
received signal as accurately as possible while taking into account the prior
distribution of the channel gains and AoAs. This is different from the prior
work [44] that use the MSE as the loss function.

In summary, the encoder of the proposed VAE takes the received signal $\bm{Y}$
as input and outputs the estimated parameters
$\bm{\mu}_{\bm{\widehat{h}}_{m}}$ and $\bm{\Sigma}_{\bm{\widehat{h}}_{m}}$ of
the optimal posterior distribution of the channel gains and the estimated AoAs
$\bm{\widehat{\theta}}$. Given these estimates, random samples of the channel
gains are drawn and fed into the decoder, who reconstructs the received signal
$\bm{\widehat{Y}}$ by (4).

###  III-B Implementation Details of the Proposed framework

Fig. 1 shows a detailed diagram of the proposed VAE. Since neural networks
generally use real numbers, all the variables in the diagram are real-valued.
We denote the real and imaginary parts of a complex variable by superscripts
$R$ and $I$, respectively. The proposed framework consists of a data
preprocessing module, a dual-path encoder, a reparameterization and recovery
module, and a decoder, which will be described in details.

####  III-B1 Data Preprocessing

The data preprocessing module has two functionalities. First, it calculates
the empirical covariance matrix from the received signals as

$\bm{R}_{\bm{Y}\bm{Y}}=\frac{1}{M}\bm{Y}\bm{Y}^{H},$ (8)

which will be used for AoA estimation.

Second, it calculates pseudo-labels $\widetilde{\bm{\theta}}$ of the AoAs,
which will be used in the initial stage of the training. The pseudo-labels
take values from a predefined set

$\mathcal{S}_{\theta}=\left\\{\underline{\theta},\underline{\theta}+\Delta_{\theta},\ldots,\bar{\theta}\right\\},$
(9)

where $\underline{\theta}$ is the minimum AoA, $\bar{\theta}$ is the maximum
AoA, and $\Delta_{\theta}$ is the precision of the pseudo-labels. Then, the
empirical correlation between the received signal $\bm{Y}$ and the array
response vectors $\bm{a}(\theta)$ is calculated, for
$\forall\theta\in\mathcal{S}_{\theta}$.

$r(\theta,\bm{Y})=\frac{1}{M}\left|\bm{1}_{M}^{T}\bm{Y}^{H}\bm{a}(\theta)\right|,$
(10)

where $\bm{1}_{M}$ is a $M$-dimensional column vector of all ones. Finally,
the set $\widetilde{\Theta}$ of $K$ angles corresponding to the $K$ largest
correlations are selected, namely

$\widetilde{\Theta}=\arg\max_{\widetilde{\Theta}\subset\mathcal{S}_{\theta},|\widetilde{\Theta}|=K}\sum_{\theta\in\widetilde{\Theta}}r(\theta,\bm{Y}).$
(11)

The pseudo-labels $\bm{\widetilde{\theta}}$ are obtained from the set
$\widetilde{\Theta}$.

####  III-B2 The Dual-Path Encoder

As shown in Fig. 1, we implement the encoder as a dual-path neural network,
which is divided into the AoA Residual Neural Network (ResNet) [56] and the
channel ResNet. The channel ResNet takes the real and imaginary parts of the
received signal matrices $\bm{Y}_{R}$ and $\bm{Y}_{I}$, and outputs the
estimated means and covariance of the channel gains,
$\bm{\mu}_{\bm{\widehat{h}}_{m}}$ and $\bm{\Sigma}_{\bm{\widehat{h}}_{m}}$ for
$m=1,\ldots,M$. The AoA ResNet takes the empirical correlation matrix
$\bm{R}_{\bm{Y}\bm{Y}}$ calculated in (8), and outputs the estimated AoAs
$\bm{\widehat{\theta}}$.

The proposed dual-path encoder separates the estimation of channel gains and
AoAs, and enables the training process that alternates between updating the
weights of the AoA ResNet and the channel ResNet. This kind of networks were
also used in [49]. As we will show in Proposition 3, the alternating training
process facilitates the convergence to the optimal estimates.

####  III-B3 Reparameterization and Recovery

As seen from the loss function (7), the expectation is taken over the
estimated parameters. Standard Monte Carlo simulations to calculate the
empirical expectation would result in a high variance of the gradient [45]. To
avoid such higher variance, a reparameterization trick has been introduced in
the literature [45].

Since it has been proved that the optimal posterior distribution of the
channel gains is Gaussian, the reparameterization trick generates standard
Gaussian random vectors
$\bm{\epsilon}_{\bm{\widehat{h}}_{m}^{R}},\bm{\epsilon}_{\bm{\widehat{h}}_{m}^{I}}\sim\mathcal{N}(\bm{0},\bm{I}_{K})$
and samples of the estimated parameters by

$\bm{\widehat{h}}_{m}^{R}=\bm{\mu}_{\bm{\widehat{h}}_{m}^{R}}+\bm{\Sigma}_{\bm{\widehat{h}}_{m}^{R}}\cdot\bm{\epsilon}_{\bm{\widehat{h}}_{m}^{R}}~{}\text{and}~{}\bm{\widehat{h}}_{m}^{I}=\bm{\mu}_{\bm{\widehat{h}}_{m}^{I}}+\bm{\Sigma}_{\bm{\widehat{h}}_{m}^{I}}\cdot\bm{\epsilon}_{\bm{\widehat{h}}_{m}^{I}}.$
(12)

Finally, the path gain and the path angle can be recovered from the estimated
parameters by:

$\bm{\widehat{\beta}}_{m}=\mathrm{abs}\left(\bm{\widehat{h}}_{m}\right)~{}\text{and}~{}\bm{\widehat{\psi}}_{m}=\arctan\left[\mathrm{diag}\left(\bm{\widehat{h}}_{m}^{R}\right)^{-1}\cdot\bm{\widehat{h}}_{m}^{I}\right],$
(13)

where $\mathrm{abs}(\bm{\widehat{h}}_{m})$ is the magnitudes of each element
in $\bm{\widehat{h}}_{m}$, and $\mathrm{diag}(\bm{\widehat{h}}_{m}^{R})$ is
the diagonal matrix with $\bm{\widehat{h}}_{m}^{R}$ as its diagonal.

####  III-B4 The Decoder

As discussed before, the knowledge-aware decoder takes the samples of the AoAs
$\bm{\widehat{\theta}}$ and path gains $\bm{\widehat{\beta}}_{m}$ and path
angles $\bm{\widehat{\psi}}_{m}$, and outputs the reconstructed signals based
on the propagation model (2). Our redesigned VAE “hardwires” the decoder to
implement the signal propagation model, instead of learning the decoder from
data as in original VAEs. This knowledge-aware decoder enforces the output of
the encoder to be the parameters to estimate, and thus break the first hurdle
of uninterpretable/unidentifiable latent variables.

###  III-C Training Process

As will be proved in Proposition 4 and Proposition 5 in Section IV, the loss
function (7) has multiple global and local optima. Hence, it is important to
start the training process from good initial points that are far away from
local minima. To obtain good initial points, the training process is divided
into two phases – the initialization phase and the normal training phase.

In the initialization phase, the network uses the pseudo-labels
$\bm{\widetilde{\theta}}$ calculated by the preprocessing module, and updates
the weights to minimize the modified loss function as follows:

$\mathcal{L}^{\left(init\right)}=\widetilde{\mathcal{L}}^{(train)}+\gamma\cdot\mathbb{E}\left[\left\|\bm{\widetilde{\theta}}-\bm{\widehat{\theta}}\right\|_{2}^{2}\right],$
(14)

where $\widetilde{\mathcal{L}}^{(train)}$ replaces the estimates
$\bm{\widehat{\theta}}$ with the pseudo-labels $\bm{\widetilde{\theta}}$ in
the loss function (7), and $\gamma$ is a scaling factor to ensure that the two
terms in (14) are on the same order. After the initialization, the encoder
neural network is trained in the training phase with the loss function
$\mathcal{L}^{(train)}$ in (7).

As seen from (14), at the end of the initialization phase, the AoA ResNet
learns to predict the pseudo-labels $\bm{\widetilde{\theta}}$. Then during the
training phase, the AoA ResNet is more likely to stay in the neighborhood of
the pseudo-labels. Our presumption is that the pseudo-labels are in close
vicinity of the true AoAs. Under this presumption, the initialization phase
with the modified loss function (14) can prevent the model from converging to
bad local optima.

###  III-D Computational Complexity During Deployment

We also analyze the computational complexity of our proposed method during the
deployment. The complexity is measured by the number of multiplications. The
input to the channel ResNet is the $N\times N$ covariance matrix with the real
and imaginary parts fed into the network by $C_{in}=2$ channels. The first
layer is a convolutional layer with a $W_{\text{kernel}}^{(1)}\times
W_{\text{kernel}}^{(1)}$ kernel, a stride $S$, and zero padding. The output of
the first layer is a $N_{\text{feature}}^{(1)}\times
N_{\text{feature}}^{(1)}\times C_{out}$ feature map, where
$N_{\text{feature}}^{(1)}=\lfloor\frac{N-W_{\text{kernel}}^{(1)}}{S}\rfloor+1\approx\frac{N}{S}$,
and $C_{out}$ is the number of output channels. The number of multiplications
of the first layer is
$(N_{\text{feature}}^{(1)})^{2}(W_{\text{kernel}}^{(1)})^{2}C_{in}C_{out}$.
Since the first layer is the largest in ResNet, the total complexity is upper
bounded by
$L(N_{\text{feature}}^{(1)})^{2}(W_{\text{kernel}}^{(1)})^{2}C_{in}C_{out}$,
where $L$ is the number of layers. In our networks, we choose the parameters
such that
$(W_{\text{kernel}}^{(1)})^{2}C_{in}C_{out}<N_{\text{feature}}^{(1)}$.
Therefore, the total complexity is upper bounded by
$L(N_{\text{feature}}^{(1)})^{3}\approx L\left(\frac{N}{S}\right)^{3}$. The
AoA ResNet is smaller than the channel ResNet because the input is smaller. In
summary, the total complexity of our method during deployment is upper bounded
by $\mathcal{O}\left(L\left(\frac{N}{S}\right)^{3}\right)$.

For the subspace-based methods such as MUSIC, the main computation complexity
comes from eigenvalue decomposition of the covariance matrix, whose complexity
is $\mathcal{O}(N^{3})$. Hence, our proposed method has similar, at least not
higher, computational complexity compared to MUSIC, and has lower complexity
compared to compressed sensing based methods.

##  IV Performance Analysis

This section analyzes the global optima and local optima of the loss function
(7), which sheds light on the design of the proposed framework and the
training process.

The multiplicity of global optima occurs when there are multiple sets of AoAs,
path gains and path angles that result in the same received signal. These
multiple sets of parameters all minimize the loss function (7). Thus, it might
seem impossible to get accurate channel estimation without supervision of true
labels. However, the following proposition shows that we can alleviate the
problem of multiple global optima by separating the AoA estimation and the
channel estimation.

###### Proposition 3

When the AoA estimation is accurate (i.e.,
$\bm{\widehat{\theta}}=\bm{\theta}$), the optimal estimates of the channel
gains that minimize the loss function in (7) are

$\bm{\mu}_{\bm{\widehat{h}}_{m}}\\!\\!=\\!\left(\bm{A}^{H}\bm{A}+\sigma^{2}\bm{\Sigma}_{\bm{h}_{m}}^{-1}\right)^{-1}\\!\left(\bm{A}^{H}\bm{A}\bm{h}_{m}+\sigma^{2}\bm{\Sigma}_{\bm{h}_{m}}^{-1}\bm{\mu}_{\bm{h}_{m}}\right).$
(15)

###### Proof:

See Appendix C. ∎

Proposition 3 indicates that once the accurate AoA estimation is obtained,
there is a unique set of channel gains that minimize the loss function. This
motivates us to design the encoder as a dual-path neural network, where the
AoA estimation is separated from the estimation of path gains and path angles.
It allows us to update the weights of the two networks in an alternating
fashion. In this way, if the AoA estimation network converges to the optimal
solution, the channel estimation network consequently has a unique global
optimum.

Note that the optimal estimates of channel gains in (15) maximize the
posterior distribution of the channel gains given the signals, and may not
equal to the true means of the channel gains. They converge to the true means
when the signal-to-noise ratio (SNR) goes to infinity.

Now that we have characterized the optimal estimation of the channel gains, we
turn to the analysis of AoA estimation. In contrast to channel estimation, the
optimal AoA estimates that minimize the loss function are not unique. This is
because the phases of the received signals and those of the reconstructed
signals are the same if for all $m=1,\ldots,M$ and $n=0,\ldots,N-1$, there
exists an integer $l_{m,n}\in\mathbb{Z}$ such that

$\left(\phi_{k,m}-\widehat{\phi}_{k,m}\right)-\frac{2\pi
nd}{\lambda}\left(\sin{\theta_{k}}-\sin{\widehat{\theta}_{k}}\right)=2\pi
l_{m,n}.$

The following proposition characterizes the multiplicity of globally optimal
AoA estimates.

(a) $N=32$ with $d/\lambda=0.5$

(b) $N=64$ with $d/\lambda=0.5$

(c) $N=32$ with $d/\lambda=2$

(d) $N=64$ with $d/\lambda=2$

Figure 2: Illustration of the loss landscape with respect to the estimates of
AoAs and path angles (estimation of path gains is assumed to be accurate).

(a) True AoAs $=(-90.0^{\circ},11.0^{\circ})$ with $d/\lambda=0.5$

(b) True AoAs $=(-90.0^{\circ},11.0^{\circ})$ with $d/\lambda=1$

Figure 3: Illustration of the loss landscape with respect to AoA estimates
(channel estimation is assumed to be accurate).

###### Proposition 4

Suppose that the channel estimation is accurate, namely
$\widehat{\bm{h}}_{m}=\bm{h}_{m}$ for $m=1,\ldots,M$. Then the AoA estimation
minimizes the loss function if each user $k$’s AoA estimation takes one of the
values:

$\widehat{\theta}_{k}=\arcsin\left(\sin\theta_{k}-l\cdot\frac{\lambda}{d}\right),$
(16)

where $l$ is an integer that satisfies

$\left\lceil\frac{d}{\lambda}\left(\sin{\theta_{k}}-1\right)\right\rceil\leq
l\leq\left\lfloor\frac{d}{\lambda}\left(\sin{\theta_{k}}+1\right)\right\rfloor,$
(17)

where $\left\lceil x\right\rceil$ is the minimum integer no smaller than $x$
and $\lfloor x\rfloor$ is the maximum integer no larger than $x$.

###### Proof:

See Appendix D. ∎

Proposition 4 indicates that even if the estimation of channel gains is
accurate, there are multiple sets of AoAs that minimize the loss function. It
also shows that the number of global optima increases with the carrier
frequency, suggesting that the problem is more severe in millimeter wave
systems due to higher band frequencies. Another important observation is that
the number of global optima does not depend on the number of antenna elements.
Therefore, simply increasing the number of antenna elements may not solve this
problem. A straightforward way to alleviate the problem of multiple global
optima is to apply a spatial filter [57] or to divide the cell into several
proper sectors [58].

Fig. 2 and Fig. 3 illustrate the multiple global optima in AoA estimation
under different numbers of antennas and antenna spacing. Fig. 2 shows the loss
landscape when we vary AoA and path angle estimates of one user with other
users’ AoA and channel estimates being accurate. Fig. 3 shows the loss
landscape when we vary the AoA estimates of two users with all other estimates
being accurate. Both figures validate the observations from Proposition 4 that
the number of global optima increases with the carrier frequency (i.e.,
$d/\lambda$) and does not decrease with the number of antennas.

The loss landscapes in Fig. 2 and Fig. 3 exhibit a plethora of stationary
points, which is a more serious issue compared to the multiplicity of global
optima. This is because the estimation at stationary points can be highly
inaccurate, and the networks, trained by first-order methods, are susceptible
to stationary points. Hence, it is important to analyze such stationary
points.

###### Proposition 5

Suppose that the channel estimation is accurate, namely
$\bm{\widehat{h}}_{m}=\bm{h_{m}}$ for all $m$, and that the AoA estimation is
accurate except for user $k$, namely $\widehat{\theta}_{i}=\theta_{i}$ for all
$i\neq k$. Given an arbitrary small number $\epsilon>0$, there exists a
threshold $\underline{N}(\epsilon)$ such that for any antenna array that has
more than $\underline{N}(\epsilon)$ antennas (i.e.,
$N>\underline{N}(\epsilon)$), the stationary points of the loss function are
within $\epsilon$ of the solutions to the following equation:

$\cos{\widehat{\theta}_{k}}\cdot\left[\frac{\cos{\left(\zeta_{\widehat{\theta}_{k}}(N-1)\right)}}{\zeta_{\widehat{\theta}_{k}}}-\frac{\cos{\left(\eta_{\widehat{\theta}_{k}}(N-1)\right)}}{\eta_{\widehat{\theta}_{k}}}\right]=0,$
(18)

where $\eta_{\widehat{\theta}_{k}}=\frac{2\pi
d}{\lambda}(\sin{\theta_{k}}+\sin{\widehat{\theta}_{k}})$ and
$\zeta_{\widehat{\theta}_{k}}=\frac{2\pi
d}{\lambda}2\sin{\widehat{\theta}_{k}}$.

###### Proof:

See Appendix E. ∎

Figure 4: Illustration of the loss landscape with respect to AoA estimates.
The true AoA is $11.0^{\circ}$ and the number of antennas is set to $32$. The
red dots are the stationary points calculated from Proposition 5.

Proposition 5 characterizes the stationary points in the regime of large
numbers of antennas, which is especially relevant for massive MIMO. It
stresses the importance of good initial points for AoA estimation.

From (18), we can see that there is always a stationary point around
$\hat{\theta}_{k}=\pm 90^{\circ}$, namely when the impinging signal is
parallel to the antenna array. This stationary point can usually be avoided
due to sectorization. Other stationary points $\hat{\theta}_{k}$ roughly
satisfy

$\frac{\cos{\left(\zeta_{\widehat{\theta}_{k}}(N-1)\right)}}{\zeta_{\widehat{\theta}_{k}}}=\frac{\cos{\left(\eta_{\widehat{\theta}_{k}}(N-1)\right)}}{\eta_{\widehat{\theta}_{k}}},$

which depends on the true AoA $\theta_{k}$, the spacing between antennas
$d/\lambda$, and the number $N$ of antennas.

Fig. 4 illustrates the loss landscape and the stationary points calculated
from Proposition 5. From the figure, there is always a stationary point around
$\pm 90^{\circ}$, and the characterization of stationary points from
Proposition 5 is close to simulation results. These observations validate
Proposition 5. Moreover, Fig. 4 shows that the density of stationary points is
higher around the true AoA, which makes it almost impossible to make an
accurate estimation with randomly initialized network parameters. This
observation underscores the importance of good initialization.

##  V Numerical Simulation

In this section, we first evaluate the performance of the proposed framework
against representative methods. Then we conduct ablation study to demonstrate
how the proposed framework overcomes the challenges of interpretability and
multiple local optima in unsupervised learning.

We set the number of users $K=3$, the number of snapshots $M=40$, the number
of antennas $N=32$, and the antenna spacing $d/\lambda=0.5$. The cell is
divided into three 120-degree sectors and the precision of the pseudo-labels
$\Delta_{\theta}=0.01^{\circ}$. The neural network is trained using the
adaptive moment estimation (Adam) algorithm with a fixed learning rate of
$10^{-4}$. Both AoAs and channel ResNets use the ResNet architecture with 18
layers as shown in [56] (called 18-layer ResNet in the paper). The batch size
of the input data is $128$. The training data set has $40000$ samples and the
test set has $10000$ samples, all generated from the signal propagation model
(2). The maximum number of epochs of the initialization phase is $200$, and
that of the training phase is $800$. Training is performed separately for each
SNR.

(a) Estimation of AoAs

(b) Estimation of path gains (PG) and path angles (PA)

Figure 5: Comparison of estimation accuracy against the canonical and an
advanced MUSIC algorithms, a state-of-the-art unsupervised learning method,
and the supervised learning method under independent and correlated channels.

###  V-A Performance Comparison

We compare the estimation accuracy of the proposed framework with the
following four methods: (1) the canonical MUSIC algorithm [4], (2) an advanced
MUSIC algorithm called KR–MUSIC [50, 51], (3) a state-of-the-art unsupervised
learning method [43] (“Unsupervised–CNN” in Fig. 5), and (4) the supervised
learning method.

Our proposed method and the supervised learning method output path gains and
path angles directly. In contrast, MUSIC algorithms and Unsupervised–CNN do
not provide estimates of path gains and path angles. Therefore, we provide a
baseline estimation of path gains and path angles via the least square (LS)
algorithm based on the AoA estimation from MUSIC under independent channels.

The performance comparison is done in two cases where the users have
independent and correlated channels, respectively. The case of correlated
channels demonstrates that the proposed framework can work beyond the typical
independent assumption made in some existing works. For the supervised
learning method, only the case of independent channels is shown because the
MSE under correlated channels is almost the same.

Fig. 5 shows that the proposed framework is superior to state-of-the-art
unsupervised learning methods proposed in [43] and achieves almost identical
performance as the supervised learning method. This demonstrates the advantage
of the proposed method: the removal of labels in our method comes at almost no
cost. This achievement is not trivial, because the classic unsupervised MUSIC
algorithm has much higher AoA estimation errors.

###  V-B Ablation Study

(a) Estimation of AoAs

(b) Estimation of path gains (PG) and path angles (PA)

Figure 6: Ablation study to show how the proposed framework overcomes the
challenge of uninterpretable latent variables.

####  V-B1 Dealing with uninterpretable latent variables

The proposed framework uses the knowledge-aware decoder, which implements the
signal propagation model, to induce the encoder to output the parameters of
interest. We perform ablation study to demonstrate the effectiveness of the
knowledge-aware decoder. When the knowledge-aware decoder with fixed
parameters is replaced with the standard decoder with learnable parameters,
the method is equivalent to the original VAE. The original VAE minimizes the
difference between the received signals and the reconstructions. Its encoder
output is used as the estimates of AoAs, path gains, and path angles. Fig. 6
shows that the AoA estimation of the original VAE has extremely higher MSEs.
Due to the inaccurate AoA estimation, the estimation of path gains and path
angles is much more inaccurate.

To further illustrate the role of the proposed knowledge-aware decoder, we
provide the original VAE with labels of true AoAs. The “VAE with AoA labels”
adds the MSE between the true AoAs and the estimated AoAs to the loss
function. In this way, some components of the encoder output are enforced to
be the AoAs, while the other components have no physical meaning. Fig. 6 shows
that due to the access to true AoAs, the VAE with AoA labels has slightly
lower MSEs than the proposed framework in AoA estimation and much higher MSE
in channel estimation.

In summary, this ablation study shows that the original VAE cannot enforce
physical meaning of the encoder output, unless labels are provided. In other
words, the proposed knowledge-aware decoder solves the issue of
uninterpretable latent variables in unsupervised learning.

####  V-B2 Dealing with multiple local optima

(a) Estimation of AoAs

(b) Estimation of path gains (PG) and path angles (PA)

Figure 7: Ablation study to show how the proposed framework overcomes the
challenge of multiple local optima.

Our theoretical analysis in Proposition 4 identifies multiple global optima
and indicates that the issue can be resolved by dividing the cell into several
proper sectors (i.e., sectorization). Our analysis in Proposition 5 also
stresses the importance of finding good initial points due to local optima
around the true AoAs.

To demonstrate how the proposed framework deals with multiple local optima, we
evaluate the performance of the proposed framework without sectorization and
initialization phase.

We also show the MSEs of the pseudo-labels obtained in the initialization
phase, in order to illustrate the starting point of our training phase. As
shown in Fig. 7, the proposed framework would have much higher MSEs without
sectorization. It also shows that the proposed initialization phase greatly
helps to improve the estimation accuracy over the case where the initial
points are randomly chosen. The curve of the proposed method without
initialization also illustrates the MSEs of local optima. The curve of the
proposed method without sectorization illustrates the MSEs of global optima
that lie in other sectors than the one the users are located at.

##  VI Conclusions

In this paper, an unsupervised learning framework based on a redesigned VAE is
proposed for joint AoA and channel estimation for massive MIMO. The proposed
framework solves two challenges of unsupervised learning in the context of AoA
and channel estimation. It solves the first challenge of
uninterpretable/unidentifiable latent variables by using a knowledge-aware
decoder. The knowledge-aware decoder implements the known signal propagation
model and enforces the encoder output to be the parameters to estimate. The
proposed framework solves the second challenge of multiple local optima by
using a dual-path encoder and adopting a two-phase training process. We
analytically derive the optimal approximate posterior distributions of the
estimated parameters and present the ELBO of the proposed framework as the
loss function. Theoretical analysis of the loss landscape is carried out to
inform our design of the framework and is validated by numerical simulations.
Numerical experiments demonstrate that our proposed framework achieves almost
the same performance as the supervised learning method and outperforms the
state-of-the-art unsupervised methods. Ablation study is also performed to
evaluate the contributions of the key components of the proposed framework in
overcoming the aforementioned challenges of unsupervised learning. Our
research exploits a new way for channel and AoA estimation, which has proved
the effectiveness of employing the unsupervised learning method in the massive
MIMO system.

Important future works include extending the convergence and performance
analysis to the multi-path scenarios, and using transfer learning to reduce
the sample complexity and training cost (e.g., train at one SNR and transfer
to other SNRs).

## Appendix A

According to [55, Lemma 4.1], the optimal approximate posterior distribution
satisfies

$\displaystyle
q_{m}^{*}\left(\bm{z}_{m}|\bm{Y}\right)=\frac{\exp{\mathbb{E}_{q_{-m}^{*}(\bm{z}_{-m}|\bm{Y})}\left[\log{p\left(\bm{z}_{m}|\bm{z}_{-m},\bm{Y}\right)}\right]}}{\int\exp{\mathbb{E}_{q_{-m}^{*}(\bm{z}_{-m}|\bm{Y})}\left[\log{p\left(\bm{z}_{m}|\bm{z}_{-m},\bm{Y}\right)}\right]}~{}d\bm{z}_{m}},$

where $\bm{z}_{-m}$ is the collection of all the latent variables except
$\bm{z}_{m}$, and $q_{-m}^{*}(\bm{z}_{-m}|\bm{Y})=\prod_{i=1,i\neq
m}^{M}q_{i}^{*}(\bm{z}_{i}|\bm{Y})$. For $m=1,\ldots,M$, the latent variable
$\bm{z}_{m}$ is the channel gain $\bm{h}_{m}$. Hence, we have

$\displaystyle\log q_{m}^{*}\left(\bm{z}_{m}|\bm{Y}\right)$ $\displaystyle=$
$\displaystyle\mathbb{E}_{q_{-m}^{*}(\bm{z}_{-m}|\bm{Y})}\left[\log{p\left(\bm{z}_{m}|\bm{z}_{-m},\bm{Y}\right)}\right]+C_{1}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{q_{-m}^{*}(\bm{z}_{-m}|\bm{Y})}\left[\log{p\left(\bm{z}_{m},\bm{z}_{-m},\bm{Y}\right)}\right]$
$\displaystyle-~{}\mathbb{E}_{q_{-m}^{*}(\bm{z}_{-m}|\bm{Y})}\left[\log{p\left(\bm{z}_{-m},\bm{Y}\right)}\right]+C_{1}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{q_{-m}^{*}(\bm{z}_{-m}|\bm{Y})}\left[\log{p\left(\bm{z}_{m},\bm{z}_{-m},\bm{Y}\right)}\right]+C_{2}$
$\displaystyle=$
$\displaystyle\mathbb{E}_{q_{-m}^{*}(\bm{z}_{-m}|\bm{Y})}\left[\textstyle\sum_{i=1}^{M}\log{p\left(\bm{y}_{i}|\bm{z}_{0},\bm{z}_{i}\right)}\right]$
$\displaystyle+~{}\mathbb{E}_{q_{-m}^{*}(\bm{z}_{-m}|\bm{Y})}\left[\textstyle\sum_{i=0}^{M}\log{p\left(\bm{z}_{i}\right)}\right]+C_{2}$
$\displaystyle=$
$\displaystyle\log{p\left(\bm{y}_{m}|\bm{z}_{0},\bm{z}_{m}\right)}+\log{p\left(\bm{z}_{m}\right)}+C_{3}$
$\displaystyle=$
$\displaystyle-\frac{1}{\sigma^{2}}\left[\bm{y}_{m}-\bm{A}(\bm{z}_{0})\cdot\bm{z}_{m}\right]^{H}\left[\bm{y}_{m}-\bm{A}(\bm{z}_{0})\cdot\bm{z}_{m}\right]$
$\displaystyle-\left(\bm{z}_{m}-\bm{\mu}_{\bm{h}_{m}}\right)^{H}\bm{\Sigma}_{\bm{h}_{m}}^{-1}\left(\bm{z}_{m}-\bm{\mu}_{\bm{h}_{m}}\right)+C_{4},$

where $C_{1}$, $C_{2}$, $C_{3}$, $C_{4}$ are constant with respect to
$q_{m}^{*}\left(\bm{z}_{m}\right)$. The last term in the above equation is a
quadratic function of $\bm{z}_{m}$. Therefore, the optimal approximate
posterior distribution $q_{m}^{*}\left(\bm{z}_{m}|\bm{Y}\right)$ is circularly
symmetric Gaussian.

## Appendix B

According to [45], the ELBO is defined as

$\displaystyle-\mathcal{D}_{KL}\left(q\left(\bm{z}|\bm{Y}\right)||p\left(\bm{z}\right)\right)+\mathbb{E}_{q\left(\bm{z}|\bm{Y}\right)}\left[\log{p\left(\bm{Y}|\bm{z}\right)}\right].$
(20)

Therefore, the loss function, defined as the negative of the ELBO, is

$\mathcal{D}_{KL}\left(q\left(\bm{z}|\bm{Y}\right)||p\left(\bm{z}\right)\right)-\mathbb{E}_{q\left(\bm{z}|\bm{Y}\right)}\left[\log{p\left(\bm{Y}|\bm{z}\right)}\right]=\mathbb{E}_{q\left(\bm{z}|\bm{Y}\right)}\left[\log\frac{q\left(\bm{z}|\bm{Y}\right)}{p\left(\bm{z}\right)}\right]-\mathbb{E}_{q\left(\bm{z}|\bm{Y}\right)}\left[\log{p\left(\bm{Y}|\bm{z}\right)}\right].$
(21)

Next, we look at each term in the above expression.

The term $\mathbb{E}_{q\left(\bm{z}|\bm{Y}\right)}\left[\log
q\left(\bm{z}|\bm{Y}\right)\right]$ is the negative of the differential
entropy of the approximate posterior distribution
$q\left(\bm{z}|\bm{Y}\right)$. Since the posterior distribution
$q_{0}\left(\bm{z}_{0}|\bm{Y}\right)$ of the AoA $\bm{z}_{0}=\bm{\theta}$ is
assumed to be a Dirac distribution, its differential entropy is zero. Since
the optimal approximate posterior distribution
$q_{m}\left(\bm{z}_{m}|\bm{Y}\right)$ of the channel gains
$\bm{z}_{m}=\bm{h}_{m}$ for $m=1,\ldots,M$ is proved to be circularly
symmetric Gaussian, its differential entropy is $\log\left[(\pi
e)^{K}|\bm{\Sigma}_{\widehat{\bm{h}}_{m}}|\right]$. Hence, we have

$\mathbb{E}_{q\left(\bm{z}|\bm{Y}\right)}\left[\log
q\left(\bm{z}|\bm{Y}\right)\right]=\textstyle\sum_{m=0}^{M}\mathbb{E}_{q_{m}\left(\bm{z}_{m}|\bm{Y}\right)}\left[\log
q_{m}\left(\bm{z}_{m}|\bm{Y}\right)\right]=-\textstyle\sum_{m=1}^{M}\log|\bm{\Sigma}_{\widehat{\bm{h}}_{m}}|+C_{1},$
(22)

where $C_{1}$ is a constant. The term
$\mathbb{E}_{q\left(\bm{z}|\bm{Y}\right)}\left[\log
p\left(\bm{z}\right)\right]$ can be calculated as

$\displaystyle\mathbb{E}_{q\left(\bm{z}|\bm{Y}\right)}\left[\log
p\left(\bm{z}\right)\right]$ $\displaystyle=$
$\displaystyle\mathbb{E}_{q_{0}\left(\bm{z}_{0}|\bm{Y}\right)}\left[\log
p\left(\bm{z}_{0}\right)\right]+\sum_{m=1}^{M}\mathbb{E}_{q_{m}\left(\bm{z}_{m}|\bm{Y}\right)}\left[\log
p\left(\bm{z}_{m}\right)\right]$ $\displaystyle=$
$\displaystyle\sum_{m=1}^{M}\mathbb{E}_{q_{m}\left(\bm{z}_{m}|\bm{Y}\right)}\left[-\left(\bm{z}_{m}-\bm{\mu}_{\bm{h}_{m}}\right)^{H}\bm{\Sigma}_{\bm{h}_{m}}^{-1}\left(\bm{z}_{m}-\bm{\mu}_{\bm{h}_{m}}\right)\right]$
$\displaystyle+C_{2},$ $\displaystyle=$
$\displaystyle-\sum_{m=1}^{M}\left(\bm{\mu}_{\widehat{\bm{h}}_{m}}-\bm{\mu}_{\bm{h}_{m}}\right)^{H}\bm{\Sigma}_{\bm{h}_{m}}^{-1}\left(\bm{\mu}_{\widehat{\bm{h}}_{m}}-\bm{\mu}_{\bm{h}_{m}}\right)$
$\displaystyle-\sum_{m=1}^{M}\mathrm{tr}\left(\bm{\Sigma}_{\bm{h}_{m}}^{-1}\bm{\Sigma}_{\widehat{\bm{h}}_{m}}\right)+C_{3},$

where $C_{2}$ and $C_{3}$ are constants.

The term $\mathbb{E}_{q\left(\bm{z}|\bm{Y}\right)}\left[\log
p\left(\bm{Y}|\bm{z}\right)\right]$ can be calculated as

$\displaystyle\mathbb{E}_{q\left(\bm{z}|\bm{Y}\right)}\left[\log
p\left(\bm{Y}|\bm{z}\right)\right]$ $\displaystyle=$
$\displaystyle\\!\\!\sum_{m=1}^{M}\mathbb{E}_{q\left(\bm{z}|\bm{Y}\right)}\left[\log
p\left(\bm{y}_{m}|\bm{z}_{0},\bm{z}_{m}\right)\right]$ $\displaystyle=$
$\displaystyle\\!\\!\sum_{m=1}^{M}\mathbb{E}_{q\left(\bm{z}|\bm{Y}\right)}\left[-\frac{\left[\bm{y}_{m}-\bm{A}(\bm{z}_{0})\bm{z}_{m}\right]^{H}\left[\bm{y}_{m}-\bm{A}(\bm{z}_{0})\bm{z}_{m}\right]}{\sigma^{2}}\right]$
$\displaystyle\\!\\!+C_{4},$ $\displaystyle=$
$\displaystyle\\!\\!\mathbb{E}_{q\left(\bm{z}|\bm{Y}\right)}\left[-\frac{\left\|\bm{Y}-\bm{A}(\bm{z_{0}})\cdot\left[\bm{z}_{1},\ldots,\bm{z}_{M}\right]\right\|_{F}^{2}}{\sigma^{2}}\right]+C_{4}$
$\displaystyle=$
$\displaystyle\\!\\!\mathbb{E}_{q\left(\bm{z}|\bm{Y}\right)}\left[-\frac{1}{\sigma^{2}}\left\|\bm{Y}-\bm{\widehat{Y}}\right\|_{F}^{2}\right]+C_{4},$

where the second equality comes from the channel model in (2), $C_{4}$ is a
constant, and the last equality comes from the definition of the decoder in
(4).

Combining (22)–(B), we obtain the loss function in (7).

## Appendix C

To get the optimal estimates of channel gains, we calculate the partial
derivative of the loss function (7) with respect to the estimates
$\bm{\mu}_{\widehat{\bm{h}}_{m}}$ of channel gains for $m=1,\ldots,M$.

We consider only the terms affected by the estimates
$\bm{\mu}_{\widehat{\bm{h}}_{m}}$, namely
$(\bm{\mu}_{\widehat{\bm{h}}_{m}}-\bm{\mu}_{\bm{h}_{m}})^{H}\bm{\Sigma}_{\bm{h}_{m}}^{-1}(\bm{\mu}_{\widehat{\bm{h}}_{m}}-\bm{\mu}_{\bm{h}_{m}})$
and
$\mathbb{E}_{q\left(\bm{z}|\bm{Y}\right)}\left[\frac{1}{\sigma^{2}}\left\|\bm{y}_{m}-\bm{\widehat{y}}_{m}\right\|_{2}^{2}\right]$.
We look at the MSE term first.

The received signal over $m$th snapshot can be written as

$\displaystyle\bm{y}_{m}=\bm{A}(\bm{\theta})\bm{h}_{m}+\bm{n}_{m}\sim\mathcal{CN}\left(\bm{A}(\bm{\theta})\bm{h}_{m},\sigma^{2}\bm{I}_{N}\right).$
(25)

The reconstructed signal can be written as
$\bm{\widehat{y}}_{m}=\bm{A}(\widehat{\bm{\theta}})\bm{\widehat{h}}_{m}$ with
$\bm{\widehat{h}}_{m}\sim\mathcal{CN}\left(\bm{\mu}_{\bm{\widehat{h}}_{m}},\bm{\Sigma}_{\bm{\widehat{h}}_{m}}\right)$.

When the AoA estimate is accurate, (i.e.,
$\bm{\widehat{\theta}}=\bm{\theta}$), we write
$\bm{A}=\bm{A}(\bm{\theta})=\bm{A}(\widehat{\bm{\theta}})$ for notational
simplicity, and have

$\displaystyle\bm{\widehat{y}}_{m}\sim\mathcal{CN}\left(\bm{A}\bm{\mu}_{\bm{\widehat{h}}_{m}},\bm{A}\bm{\Sigma}_{\bm{\widehat{h}}_{m}}\bm{A}^{H}\right).$
(26)

Since the noise $\mathbf{n}_{m}$ and the channel gains $\bm{\widehat{h}}_{m}$
are independent, $\bm{y}_{m}$ and $\bm{\widehat{y}}_{m}$ are independent.
Therefore, we have

$\displaystyle\bm{y}_{m}-\bm{\widehat{y}}_{m}\sim\mathcal{CN}\left(\bm{A}(\bm{h}_{m}-\bm{\mu}_{\bm{\widehat{h}}_{m}}),\sigma^{2}\bm{I}_{N}+\bm{A}\bm{\Sigma}_{\bm{\widehat{h}}_{m}}\bm{A}^{H}\right).$

As a result, the MSE term can be calculated as

$\mathbb{E}\left[\left\|\bm{y}_{m}-\bm{\widehat{y}}_{m}\right\|_{2}^{2}\right]=(\bm{h}_{m}-\bm{\mu}_{\bm{\widehat{h}}_{m}})^{H}\bm{A}^{H}\bm{A}(\bm{h}_{m}-\bm{\mu}_{\bm{\widehat{h}}_{m}})+~{}\sigma^{2}N+\mathrm{tr}\left(\bm{A}\bm{\Sigma}_{\bm{\widehat{h}}_{m}}\bm{A}^{H}\right).$
(27)

We collect the terms that depend on the estimate
$\bm{\mu}_{\bm{\widehat{h}}_{m}}$ into

$f(\bm{\mu}_{\bm{\widehat{h}}_{m}})=(\bm{\mu}_{\widehat{\bm{h}}_{m}}-\bm{\mu}_{\bm{h}_{m}})^{H}\bm{\Sigma}_{\bm{h}_{m}}^{-1}(\bm{\mu}_{\widehat{\bm{h}}_{m}}-\bm{\mu}_{\bm{h}_{m}})+~{}\frac{1}{\sigma^{2}}(\bm{h}_{m}-\bm{\mu}_{\bm{\widehat{h}}_{m}})^{H}\bm{A}^{H}\bm{A}(\bm{h}_{m}-\bm{\mu}_{\bm{\widehat{h}}_{m}}),$
(28)

which is a concave quadratic function of the estimate
$\bm{\mu}_{\bm{\widehat{h}}_{m}}$. Therefore, the estimate
$\bm{\mu}_{\bm{\widehat{h}}_{m}}$ minimizes
$f(\bm{\mu}_{\bm{\widehat{h}}_{m}})$ if and only if the derivatives with
respect to the real part $\bm{\mu}_{\bm{\widehat{h}}_{m}}^{R}$ and the
imaginary part $\bm{\mu}_{\bm{\widehat{h}}_{m}}^{I}$ are zero. According to
the $\mathbb{CR}$-calculus [59], the derivatives can be calculated as

$\displaystyle\frac{\partial
f(\bm{\mu}_{\bm{\widehat{h}}_{m}})}{\partial\bm{\mu}_{\bm{\widehat{h}}_{m}}^{R}}\\!\\!\\!\\!$
$\displaystyle=$ $\displaystyle\\!\\!\\!\\!\frac{\partial
f(\bm{\mu}_{\bm{\widehat{h}}_{m}},\bm{\mu}_{\bm{\widehat{h}}_{m}}^{*})}{\partial\bm{\mu}_{\bm{\widehat{h}}_{m}}}+\frac{\partial
f(\bm{\mu}_{\bm{\widehat{h}}_{m}},\bm{\mu}_{\bm{\widehat{h}}_{m}}^{*})}{\partial\bm{\mu}_{\bm{\widehat{h}}_{m}}^{*}}$
$\displaystyle=$
$\displaystyle\\!\\!\\!\\!2\cdot\bm{\mathfrak{Re}}\left\\{\bm{\Sigma}_{\bm{h}_{m}}^{-1}(\bm{\mu}_{\widehat{\bm{h}}_{m}}-\bm{\mu}_{\bm{h}_{m}})-\frac{1}{\sigma^{2}}\bm{A}^{H}\bm{A}(\bm{h}_{m}-\bm{\mu}_{\bm{\widehat{h}}_{m}})\right\\}$

and

$\displaystyle\frac{\partial
f(\bm{\mu}_{\bm{\widehat{h}}_{m}})}{\partial\bm{\mu}_{\bm{\widehat{h}}_{m}}^{I}}\\!\\!\\!\\!$
$\displaystyle=$ $\displaystyle\\!\\!\\!\\!j\\!\cdot\\!\left[\frac{\partial
f(\bm{\mu}_{\bm{\widehat{h}}_{m}},\bm{\mu}_{\bm{\widehat{h}}_{m}}^{*})}{\partial\bm{\mu}_{\bm{\widehat{h}}_{m}}}-\frac{\partial
f(\bm{\mu}_{\bm{\widehat{h}}_{m}},\bm{\mu}_{\bm{\widehat{h}}_{m}}^{*})}{\partial\bm{\mu}_{\bm{\widehat{h}}_{m}}^{*}}\right]$
$\displaystyle=$
$\displaystyle\\!\\!\\!\\!2\cdot\bm{\mathfrak{Im}}\left\\{\bm{\Sigma}_{\bm{h}_{m}}^{-1}(\bm{\mu}_{\widehat{\bm{h}}_{m}}-\bm{\mu}_{\bm{h}_{m}})-\frac{1}{\sigma^{2}}\bm{A}^{H}\bm{A}(\bm{h}_{m}-\bm{\mu}_{\bm{\widehat{h}}_{m}})\right\\},$

where $\bm{\mathfrak{Re}}\left\\{\cdot\right\\}$ and
$\bm{\mathfrak{Im}}\left\\{\cdot\right\\}$ are the real and imaginary parts.

Combining (C) and (C), we know that the optima estimate should satisfy

$\displaystyle\left(\bm{A}^{H}\bm{A}+\sigma^{2}\bm{\Sigma}_{\bm{h}_{m}}^{-1}\right)\bm{\mu}_{\bm{\widehat{h}}_{m}}=\bm{A}^{H}\bm{A}\bm{h}_{m}+\sigma^{2}\bm{\Sigma}_{\bm{h}_{m}}^{-1}\bm{\mu}_{\bm{h}_{m}}.$
(31)

Since the matrix
$\left(\bm{A}^{H}\bm{A}+\sigma^{2}\bm{\Sigma}_{\bm{h}_{m}}^{-1}\right)$ is
positive definitive and hence invertible, we have

$\displaystyle\bm{\mu}_{\bm{\widehat{h}}_{m}}\\!\\!=\\!\left(\bm{A}^{H}\bm{A}+\sigma^{2}\bm{\Sigma}_{\bm{h}_{m}}^{-1}\right)^{-1}\\!\left(\bm{A}^{H}\bm{A}\bm{h}_{m}+\sigma^{2}\bm{\Sigma}_{\bm{h}_{m}}^{-1}\bm{\mu}_{\bm{h}_{m}}\right).$
(32)

## Appendix D

Under the assumption that the channel estimation is accurate (i.e.,
$\bm{\widehat{h}}_{m}=\bm{h}_{m}$ for $m=1,\ldots,M$, one type of globally
optimal AoA estimates $\bm{\widehat{\theta}}$ are the ones that yield the same
phase shifts of the received signals as the true AoAs $\bm{\theta}$ do.
Mathematically, this condition is that for all $n=1,\ldots,N$, there exists an
integer $l_{n}\in\mathbb{Z}$ such that

$\frac{2\pi(n-1)d}{\lambda}(\sin{\theta_{k}}-\sin{\widehat{\theta}_{k}})=2\pi
l_{n}.$ (33)

From (33), we have $l_{1}=0$,
$l_{2}=\frac{d}{\lambda}(\sin{\theta_{k}}-\sin{\widehat{\theta}_{k}})$, and
$l_{n}=(n-1)l_{2}$ for $n=3,4,\ldots,N$. So it boils down to finding
$\widehat{\theta}_{k}$ such that
$\frac{d}{\lambda}(\sin{\theta_{k}}-\sin{\widehat{\theta}_{k}})$ is an
integer.

Since
$\sin{\theta_{k}}-\sin{\widehat{\theta}_{k}}\in\left[\sin{\theta_{k}}-1,\sin{\theta_{k}}+1\right]$,
we have

$\left\lceil\frac{d}{\lambda}\left(\sin{\theta_{k}}-1\right)\right\rceil\leq
l_{2}\leq\left\lfloor\frac{d}{\lambda}\left(\sin{\theta_{k}}+1\right)\right\rfloor,$
(34)

where $\left\lceil x\right\rceil$ is the minimum integer no smaller than $x$
and $\lfloor x\rfloor$ is the maximum integer no larger than $x$.

Therefore, the solutions to (33) are

${\widehat{\theta}_{k}}=\arcsin\left(\sin\theta_{k}-l_{2}\cdot\frac{\lambda}{d}\right)$
(35)

for integer $l_{2}$ that satisfies (34).

## Appendix E

###  E-A Calculating the Gradient

To study the stationary points, we first calculate the gradient of the loss
function with respect to the AoA estimates $\bm{\widehat{\theta}}$. Note that
the AoA estimates only affect the MSE term
$\mathbb{E}\|\bm{Y}-\bm{\widehat{Y}}\|_{F}^{2}$ in the loss function (7).

Since the received signal in the $m$th snapshot is

$\displaystyle\bm{y}_{m}=\bm{A}(\bm{\theta})\bm{h}_{m}+\bm{n}_{m}\sim\mathcal{CN}\left(\bm{A}(\bm{\theta})\bm{h}_{m},\sigma^{2}\bm{I}_{N}\right),$
(36)

and since the reconstructed signal is

$\displaystyle\bm{\widehat{y}}_{m}=\bm{A}(\widehat{\bm{\theta}})\bm{\widehat{h}}_{m}\sim\mathcal{CN}\left(\bm{A}(\widehat{\bm{\theta}})\bm{\mu}_{\bm{\widehat{h}}_{m}},\bm{A}(\widehat{\bm{\theta}})\bm{\Sigma}_{\bm{\widehat{h}}_{m}}\bm{A}(\widehat{\bm{\theta}})^{H}\right),$
(37)

we have

$\displaystyle\bm{y}_{m}-\bm{\widehat{y}}_{m}\sim\mathcal{CN}\left(\bm{A}\bm{h}_{m}-\bm{\widehat{A}}\bm{\mu}_{\bm{\widehat{h}}_{m}},\sigma^{2}\bm{I}_{N}+\bm{\widehat{A}}\bm{\Sigma}_{\bm{\widehat{h}}_{m}}\bm{\widehat{A}}^{H}\right),$

where we define $\bm{A}=\bm{A}(\bm{\theta})$ and
$\bm{\widehat{A}}=\bm{A}(\widehat{\bm{\theta}})$ for simplicity.

Thus, the MSE can be calculated as

$\mathbb{E}\left[\left\|\bm{y}_{m}-\bm{\widehat{y}}_{m}\right\|_{2}^{2}\right]=(\bm{A}\bm{h}_{m}-\bm{\widehat{A}}\bm{\mu}_{\bm{\widehat{h}}_{m}})^{H}(\bm{A}\bm{h}_{m}-\bm{\widehat{A}}\bm{\mu}_{\bm{\widehat{h}}_{m}})\\\
+\sigma^{2}N+\mathrm{tr}\left(\bm{\widehat{A}}\bm{\Sigma}_{\bm{\widehat{h}}_{m}}\bm{\widehat{A}}^{H}\right).$

Since the AoA estimates $\bm{\widehat{\theta}}$ affect $\bm{\widehat{A}}$
only, we collect the terms in (E-A) that depend on $\bm{\widehat{A}}$ as
follows:

$f_{m}(\bm{\widehat{A}})\triangleq\underbrace{\left(-\bm{h}_{m}^{H}\bm{A}^{H}\bm{\widehat{A}}\bm{\mu}_{\bm{\widehat{h}}_{m}}-\bm{\mu}_{\bm{\widehat{h}}_{m}}^{H}\bm{\widehat{A}}^{H}\bm{A}\bm{h}_{m}\right)}_{\triangleq
f_{m,1}(\bm{\widehat{A}})}+\underbrace{\bm{\mu}_{\bm{\widehat{h}}_{m}}^{H}\bm{\widehat{A}}^{H}\bm{\widehat{A}}\bm{\mu}_{\bm{\widehat{h}}_{m}}}_{\triangleq
f_{m,2}(\bm{\widehat{A}})}+\underbrace{\mathrm{tr}\left(\bm{\widehat{A}}\bm{\Sigma}_{\bm{\widehat{h}}_{m}}\bm{\widehat{A}}^{H}\right)}_{\triangleq
f_{m,3}(\bm{\widehat{A}})}$ (38)

Note that the matrix
$\bm{\widehat{A}}=\left[\bm{\widehat{a}}_{1},\ldots,\bm{\widehat{a}}_{K}\right]$
has $K$ columns, where the $k$-th column $\bm{\widehat{a}}_{k}$ is a function
of user $k$’s AoA estimate $\widehat{\theta}_{k}$, namely
$\bm{\widehat{a}}_{k}=e^{-j\bm{\Psi}(\widehat{\theta}_{k})}$.

Next, we calculate the partial derivative of the function
$f(\bm{\widehat{A}})$ in (38) with respect to user $k$’s AoA estimate
$\widehat{\theta}_{k}$ using the $\mathbb{CR}$-calculus [59]:

$\displaystyle\frac{\partial
f(\bm{\widehat{A}})}{\partial\widehat{\theta}_{k}}=\left[\frac{\partial
f(\bm{\widehat{A}})}{\partial\bm{\widehat{a}}_{k}}\right]^{T}\cdot\frac{\partial\bm{\widehat{a}}_{k}}{\partial\widehat{\theta}_{k}}+\left[\frac{\partial
f(\bm{\widehat{A}})}{\partial\bm{\widehat{a}}_{k}^{*}}\right]^{T}\cdot\frac{\partial\bm{\widehat{a}}_{k}^{*}}{\partial\widehat{\theta}_{k}}.$
(39)

To calculate $\frac{\partial
f(\bm{\widehat{A}})}{\partial\bm{\widehat{a}}_{k}}$ and $\frac{\partial
f(\bm{\widehat{A}})}{\partial\bm{\widehat{a}}_{k}^{*}}$, it is useful to note
that

$\displaystyle\bm{\widehat{A}}\bm{\mu}_{\mathbf{\widehat{h}}_{m}}=\sum_{k=1}^{K}\mu_{\mathbf{\widehat{h}}_{m},k}\bm{\widehat{a}}_{k}~{}\text{and}~{}\bm{\widehat{A}}^{*}\bm{\mu}_{\mathbf{\widehat{h}}_{m}}^{*}=\sum_{k=1}^{K}\mu_{\mathbf{\widehat{h}}_{m},k}^{*}\bm{\widehat{a}}_{k}^{*}.$
(40)

For the first term in (38), we have

$\displaystyle\frac{\partial
f_{m,1}(\bm{\widehat{A}})}{\partial\bm{\widehat{a}}_{k}}=\frac{\partial\left(-\bm{h}_{m}^{H}\bm{A}^{H}\bm{\widehat{A}}\bm{\mu}_{\mathbf{\widehat{h}}_{m}}\right)}{\partial\bm{\widehat{a}}_{k}}=-\mu_{\mathbf{\widehat{h}}_{m},k}\bm{A}^{*}\bm{h}_{m}^{*},$
(41)

and

$\displaystyle\frac{\partial
f_{m,1}(\bm{\widehat{A}})}{\partial\bm{\widehat{a}}_{k}^{*}}=\frac{\partial\left(-\bm{\mu}_{\mathbf{\widehat{h}}_{m}}^{H}\bm{\widehat{A}}^{H}\bm{A}\bm{h}_{m}\right)}{\partial\bm{\widehat{a}}_{k}^{*}}=-\mu_{\mathbf{\widehat{h}}_{m},k}^{*}\bm{A}\bm{h}_{m}.$
(42)

For the second term in (38), we have

$\displaystyle\frac{\partial
f_{m,2}(\bm{\widehat{A}})}{\partial\bm{\widehat{a}}_{k}}=\frac{\partial\left(\bm{\mu}_{\mathbf{\widehat{h}}_{m}}^{H}\bm{\widehat{A}}^{H}\bm{\widehat{A}}\bm{\mu}_{\mathbf{\widehat{h}}_{m}}\right)}{\partial\bm{\widehat{a}}_{k}}=\mu_{\mathbf{\widehat{h}}_{m},k}\bm{\widehat{A}}^{*}\bm{\mu}_{\mathbf{\widehat{h}}_{m}}^{*},$
(43)

and

$\displaystyle\frac{\partial
f_{m,2}(\bm{\widehat{A}})}{\partial\bm{\widehat{a}}_{k}^{*}}=\frac{\partial\left(\bm{\mu}_{\mathbf{\widehat{h}}_{m}}^{H}\bm{\widehat{A}}^{H}\bm{\widehat{A}}\bm{\mu}_{\mathbf{\widehat{h}}_{m}}\right)}{\partial\bm{\widehat{a}}_{k}^{*}}=\mu_{\mathbf{\widehat{h}}_{m},k}^{*}\bm{\widehat{A}}\bm{\mu}_{\mathbf{\widehat{h}}_{m}}.$
(44)

For the third term in (38), we have

$\displaystyle\frac{\partial
f_{m,3}(\bm{\widehat{A}})}{\partial\bm{\widehat{a}}_{k}}=\frac{\partial~{}\mathrm{tr}\left(\bm{\widehat{A}}\bm{\Sigma}_{\mathbf{\widehat{h}}_{m}}\bm{\widehat{A}}^{H}\right)}{\partial\bm{\widehat{a}}_{k}}=\left(\bm{\widehat{A}}^{*}\bm{\Sigma}_{\mathbf{\widehat{h}}_{m}}^{*}\right)_{:,k},$
(45)

and

$\displaystyle\frac{\partial
f_{m,3}(\bm{\widehat{A}})}{\partial\bm{\widehat{a}}_{k}^{*}}=\frac{\partial~{}\mathrm{tr}\left(\bm{\widehat{A}}\bm{\Sigma}_{\mathbf{\widehat{h}}_{m}}\bm{\widehat{A}}^{H}\right)}{\partial\bm{\widehat{a}}_{k}^{*}}=\left(\bm{\widehat{A}}\bm{\Sigma}_{\mathbf{\widehat{h}}_{m}}\right)_{:,k},$
(46)

where $\left(\cdot\right)_{:,k}$ is the $k$-th column of a matrix.

In addition, we have

$\displaystyle\frac{\partial\bm{\widehat{a}}_{k}}{\partial\widehat{\theta}_{k}}=-j\cdot\bm{\widehat{a}}_{k}\odot\bm{\Phi}(\widehat{\theta}_{k})~{}\text{and}~{}\frac{\partial\bm{\widehat{a}}_{k}^{*}}{\partial\widehat{\theta}_{k}}=j\cdot\bm{\widehat{a}}_{k}^{*}\odot\bm{\Phi}(\widehat{\theta}_{k}),$
(47)

where $\bm{\Phi}(\widehat{\theta}_{k})=\frac{2\pi
d\cos(\widehat{\theta}_{k})}{\lambda}[0,1,\ldots,N-1]^{T}$.

Combining (41)–(47), we have

$\displaystyle\frac{\partial
f_{m}(\bm{\widehat{A}})}{\partial\widehat{\theta}_{k}}$ $\displaystyle=$
$\displaystyle\left[-\mu_{\mathbf{\widehat{h}}_{m},k}\bm{A}^{*}\bm{h}_{m}^{*}+\mu_{\mathbf{\widehat{h}}_{m},k}\bm{\widehat{A}}^{*}\bm{\mu}_{\mathbf{\widehat{h}}_{m}}^{*}+\left(\bm{\widehat{A}}^{*}\bm{\Sigma}_{\mathbf{\widehat{h}}_{m}}^{*}\right)_{:,k}\right]^{T}\cdot\left[-j\cdot\bm{\widehat{a}}_{k}\odot\bm{\Phi}(\widehat{\theta}_{k})\right]$
$\displaystyle+\left[-\mu_{\mathbf{\widehat{h}}_{m},k}^{*}\bm{A}\bm{h}_{m}+\mu_{\mathbf{\widehat{h}}_{m},k}^{*}\bm{\widehat{A}}\bm{\mu}_{\mathbf{\widehat{h}}_{m}}+\left(\bm{\widehat{A}}\bm{\Sigma}_{\mathbf{\widehat{h}}_{m}}\right)_{:,k}\right]^{T}\cdot\left[j\cdot\bm{\widehat{a}}_{k}^{*}\odot\bm{\Phi}(\widehat{\theta}_{k})\right]$
$\displaystyle=$
$\displaystyle\bm{\mathfrak{Im}}\left\\{\left[\bm{J}_{m}(\bm{\widehat{A}})\right]^{T}\cdot\left[\bm{\widehat{a}}_{k}\odot\bm{\Phi}(\widehat{\theta}_{k})\right]\right\\},$

where

$\bm{J}_{m}(\bm{\widehat{A}})\triangleq\mu_{\mathbf{\widehat{h}}_{m},k}\left(\bm{\widehat{A}}^{*}\bm{\mu}_{\mathbf{\widehat{h}}_{m}}^{*}-\bm{A}^{*}\bm{h}_{m}^{*}\right)+\left(\bm{\widehat{A}}^{*}\bm{\Sigma}_{\mathbf{\widehat{h}}_{m}}^{*}\right)_{:,k}.$
(49)

Under the assumption that the channel estimation is accurate (i.e.,
$\bm{\mu}_{\bm{\widehat{h}}_{m}}=\bm{h}_{m}$ for $m=1,\ldots,M$), the
derivative can be finally written as

$\displaystyle\\!\\!\frac{\partial\mathcal{L}^{(train)}}{\partial\widehat{\theta}_{k}}=\sum_{m=1}^{M}\frac{\partial
f_{m}(\bm{\widehat{A}})}{\partial\widehat{\theta}_{k}}$ $\displaystyle=$
$\displaystyle\\!\\!\sum_{m=1}^{M}\bm{\mathfrak{Im}}\left\\{\bm{J}_{m}(\bm{\widehat{A}})^{T}\cdot\left[\bm{\widehat{a}}_{k}\odot\bm{\Phi}(\widehat{\theta}_{k})\right]\right\\}$
$\displaystyle=$
$\displaystyle\\!\\!\sum_{m=1}^{M}|h_{m,k}|^{2}\sum_{n=0}^{N-1}\frac{2\pi
dn}{\lambda}\cos{\widehat{\theta}_{k}}\left[\sin{(\eta_{\widehat{\theta}_{k}}n)}-\sin{(\zeta_{\widehat{\theta}_{k}}n)}\right],$

where $\eta_{\widehat{\theta}_{k}}=\frac{2\pi
d}{\lambda}(\sin{\theta_{k}}+\sin{\widehat{\theta}_{k}})$ and
$\zeta_{\widehat{\theta}_{k}}=\frac{2\pi
d}{\lambda}2\sin{\widehat{\theta}_{k}}$.

###  E-B Analyzing the stationary points

In the gradient (E-A), the terms $\sum_{m=1}^{M}|h_{m,k}|^{2}$ and $\frac{2\pi
d}{\lambda}$ are always positive. Therefore, we focus on

$\displaystyle
g(\widehat{\theta}_{k})\triangleq\cos{\widehat{\theta}_{k}}\cdot\sum_{n=0}^{N-1}n\left[\sin{(\eta_{\widehat{\theta}_{k}}n)}-\sin{(\zeta_{\widehat{\theta}_{k}}n)}\right].$
(51)

The stationary points are the solutions to $g(\widehat{\theta}_{k})=0$.

Since $g(\widehat{\theta}_{k})$ is a summation of $N$ terms, we aim to find
its limit when $N\rightarrow\infty$. Specifically, define

$t(x)\triangleq
x\left[\sin{\left(\eta_{\widehat{\theta}_{k}}(N-1)x\right)}-\sin{\left(\zeta_{\widehat{\theta}_{k}}(N-1)x\right)}\right].$
(52)

Note that the relationship between $g(\widehat{\theta}_{k})$ and $t(x)$ is

$g(\widehat{\theta}_{k})=\cos{\widehat{\theta}_{k}}\cdot(N-1)\cdot\sum_{n=0}^{N-1}t\left(\frac{n}{N-1}\right).$
(53)

Meanwhile, we have

$\int_{0}^{1}t(x)dx=\lim_{N\to\infty}\frac{1}{N-1}\sum_{n=0}^{N-1}t\left(\frac{n}{N-1}\right).$
(54)

Note that the integral can be calculated analytically as

$\int_{0}^{1}t(x)dx=\frac{\cos{\left(\zeta_{\widehat{\theta}_{k}}(N-1)\right)}}{\zeta_{\widehat{\theta}_{k}}}-\frac{\cos{\left(\eta_{\widehat{\theta}_{k}}(N-1)\right)}}{\eta_{\widehat{\theta}_{k}}}.$
(55)

Finally, we have

$\displaystyle\\!\\!\lim_{N\to\infty}g(\widehat{\theta}_{k})=0$
$\displaystyle\Leftrightarrow$
$\displaystyle\\!\\!\lim_{N\to\infty}\cos{\widehat{\theta}_{k}}\cdot\sum_{n=0}^{N-1}t\left(\frac{n}{N-1}\right)=0$
$\displaystyle\Leftrightarrow$
$\displaystyle\\!\\!\cos{\widehat{\theta}_{k}}\cdot\int_{0}^{1}t(x)dx=0$
$\displaystyle\Leftrightarrow$
$\displaystyle\\!\\!\cos{\widehat{\theta}_{k}}\cdot\left[\frac{\cos{\left(\zeta_{\widehat{\theta}_{k}}(N-1)\right)}}{\zeta_{\widehat{\theta}_{k}}}-\frac{\cos{\left(\eta_{\widehat{\theta}_{k}}(N-1)\right)}}{\eta_{\widehat{\theta}_{k}}}\right]=0,$

which concludes the proof.

## References

  * [1] T. L. Marzetta, “Noncooperative cellular wireless with unlimited numbers of base station antennas,” _IEEE Transactions on Wireless Communications_ , vol. 9, no. 11, pp. 3590–3600, 2010.
  * [2] L. Liu, C. Oestges, J. Poutanen, K. Haneda, P. Vainikainen, F. Quitin, F. Tufvesson, and P. D. Doncker, “The COST 2100 MIMO channel model,” _IEEE Wireless Communications_ , vol. 19, no. 6, pp. 92–99, 2012.
  * [3] J. Li, Q. Zhang, X. Xin, Y. Tao, Q. Tian, F. Tian, D. Chen, Y. Shen, G. Cao, Z. Gao, and J. Qian, “Deep learning-based massive MIMO CSI feedback,” in _2019 18th International Conference on Optical Communications and Networks (ICOCN)_ , 2019, pp. 1–3.
  * [4] P. Stoica and A. Nehorai, “MUSIC, maximum likelihood and cramer-rao bound: further results and comparisons,” in _International Conference on Acoustics, Speech, and Signal Processing,_ , 1989, pp. 2605–2608 vol.4.
  * [5] R. Roy and T. Kailath, “ESPRIT – estimation of signal parameters via rotational invariance techniques,” _IEEE Transactions on acoustics, speech, and signal processing_ , vol. 37, no. 7, pp. 984–995, 1989.
  * [6] S. Marcos, A. Marsal, and M. Benidir, “The propagator method for source bearing estimation,” _Signal Processing_ , vol. 42, no. 2, pp. 121––138, Mar. 1995.
  * [7] P. Pal and P. P. Vaidyanathan, “Nested arrays: A novel approach to array processing with enhanced degrees of freedom,” _IEEE Transactions on Signal Processing_ , vol. 58, no. 8, pp. 4167–4181, 2010.
  * [8] X. Dai, X. Zhang, and Y. Wang, “Extended DOA-matrix method for DOA estimation via two parallel linear arrays,” _IEEE Communications Letters_ , vol. 23, no. 11, pp. 1981–1984, 2019.
  * [9] B. Li, S. Wang, J. Zhang, X. Cao, and C. Zhao, “Fast randomized-MUSIC for mm-Wave massive MIMO radars,” _IEEE Transactions on Vehicular Technology_ , vol. 70, no. 2, pp. 1952–1956, 2021.
  * [10] H. Liu, X. Yuan, and Y. J. Zhang, “Super-resolution blind channel-and-signal estimation for massive MIMO with one-dimensional antenna array,” _IEEE Transactions on Signal Processing_ , vol. 67, no. 17, pp. 4433–4448, 2019.
  * [11] P. Chen, Z. Chen, Z. Cao, and X. Wang, “A new atomic norm for DOA estimation with gain-phase errors,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 4293–4306, 2020.
  * [12] Z. Yang and L. Xie, “Exact joint sparse frequency recovery via optimization methods,” _IEEE Transactions on Signal Processing_ , vol. 64, no. 19, pp. 5145–5157, 2016.
  * [13] A. Garcia-Rodriguez and C. Masouros, “Low-complexity compressive sensing detection for spatial modulation in large-scale multiple access channels,” _IEEE Transactions on Communications_ , vol. 63, no. 7, pp. 2565–2579, 2015\.
  * [14] Y. Wang, Y. Zhang, Z. Tian, G. Leus, and G. Zhang, “Super-resolution channel estimation for arbitrary arrays in hybrid millimeter-wave massive MIMO systems,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 13, no. 5, pp. 947–960, 2019.
  * [15] R. W. Heath, N. González-Prelcic, S. Rangan, W. Roh, and A. M. Sayeed, “An overview of signal processing techniques for millimeter wave MIMO systems,” _IEEE Journal of Selected Topics in Signal Processing_ , vol. 10, no. 3, pp. 436–453, 2016.
  * [16] Z. Gao, L. Dai, S. Han, C.-L. I, Z. Wang, and L. Hanzo, “Compressive sensing techniques for next-generation wireless communications,” _IEEE Wireless Communications_ , vol. 25, no. 3, pp. 144–153, 2018.
  * [17] 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 Transactions on Signal Processing_ , vol. 68, pp. 764–779, 2020.
  * [18] M.-A. Badiu, T. L. Hansen, and B. H. Fleury, “Variational Bayesian inference of line spectra,” _IEEE Transactions on Signal Processing_ , vol. 65, no. 9, pp. 2247–2261, 2017.
  * [19] T. L. Hansen, B. H. Fleury, and B. D. Rao, “Superfast line spectral estimation,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 10, pp. 2511–2526, 2018.
  * [20] T. J. O’Shea, T. Roy, and N. West, “Approximating the void: Learning stochastic channel models from observation with variational generative adversarial networks,” in _2019 International Conference on Computing, Networking and Communications (ICNC)_ , 2019, pp. 681–686.
  * [21] H. Ye, G. Y. Li, B.-H. F. Juang, and K. Sivanesan, “Channel agnostic end-to-end learning based communication systems with conditional gan,” in _2018 IEEE Globecom Workshops (GC Wkshps)_ , 2018, pp. 1–5.
  * [22] V. Raj and S. Kalyani, “Backpropagating through the air: Deep learning at physical layer without channel models,” _IEEE Communications Letters_ , vol. 22, no. 11, pp. 2278–2281, 2018.
  * [23] Z. Qin, H. Ye, G. Y. Li, and B.-H. F. Juang, “Deep learning in physical layer communications,” _IEEE Wireless Communications_ , vol. 26, no. 2, pp. 93–99, 2019.
  * [24] H. He, S. Jin, C.-K. Wen, F. Gao, G. Y. Li, and Z. Xu, “Model-driven deep learning for physical layer communications,” _IEEE Wireless Communications_ , vol. 26, no. 5, pp. 77–83, 2019.
  * [25] C.-J. Chun, J.-M. Kang, and I.-M. Kim, “Deep learning-based channel estimation for massive mimo systems,” _IEEE Wireless Communications Letters_ , vol. 8, no. 4, pp. 1228–1231, 2019.
  * [26] J.-M. Kang, C.-J. Chun, and I.-M. Kim, “Deep learning based channel estimation for MIMO systems with received SNR feedback,” _IEEE Access_ , vol. 8, pp. 121 162–121 181, 2020.
  * [27] Z. Guo, K. Lin, X. Chen, and C.-Y. Chit, “Transfer learning for angle of arrivals estimation in massive mimo system,” in _2022 IEEE/CIC International Conference on Communications in China (ICCC)_ , 2022, pp. 506–511.
  * [28] L. Wu, Z.-M. Liu, and Z.-T. Huang, “Deep convolution network for direction of arrival estimation with sparse prior,” _IEEE Signal Processing Letters_ , vol. 26, no. 11, pp. 1688–1692, 2019.
  * [29] X. Su, P. Hu, Z. Liu, T. Liu, B. Peng, and X. Li, “Mixed near-field and far-field source localization based on convolution neural networks via symmetric nested array,” _IEEE Transactions on Vehicular Technology_ , vol. 70, no. 8, pp. 7908–7920, 2021.
  * [30] P. Pan, Y. Zhang, Z. Deng, and G. Wu, “Complex-valued frequency estimation network and its applications to superresolution of radar range profiles,” _IEEE Transactions on Geoscience and Remote Sensing_ , 2021.
  * [31] Y. Guo, Z. Zhang, Y. Huang, and P. Zhang, “DOA estimation method based on cascaded neural network for two closely spaced sources,” _IEEE Signal Processing Letters_ , vol. 27, pp. 570–574, 2020.
  * [32] A. Barthelme and W. Utschick, “A machine learning approach to DoA estimation and model order selection for antenna arrays with subarray sampling,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 3075–3087, 2021.
  * [33] P. Pan, Y. Zhang, Z. Deng, and W. Qi, “Deep learning-based 2-D frequency estimation of multiple sinusoidals,” _IEEE Transactions on Neural Networks and Learning Systems_ , 2021.
  * [34] L. Wan, Y. Sun, L. Sun, Z. Ning, and J. J. Rodrigues, “Deep learning based autonomous vehicle super resolution DOA estimation for safety driving,” _IEEE Transactions on Intelligent Transportation Systems_ , 2020.
  * [35] R. Akter, V.-S. Doan, T. Huynh-The, and D.-S. Kim, “RFDOA-Net: An efficient ConvNet for RF-based DOA estimation in UAV surveillance systems,” _IEEE Transactions on Vehicular Technology_ , 2021.
  * [36] S. Chakraborty and D. Saha, “Domain knowledge aided neural network for wireless channel estimation,” in _2021 IEEE Global Communications Conference (GLOBECOM)_ , 2021, pp. 1–6.
  * [37] X. Gao, S. Jin, C.-K. Wen, and G. Y. Li, “ComNet: Combination of deep learning and expert knowledge in OFDM receivers,” _IEEE Communications Letters_ , vol. 22, no. 12, pp. 2627–2630, 2018.
  * [38] J. Liao, J. Zhao, F. Gao, and G. Y. Li, “A model-driven deep learning method for massive MIMO detection,” _IEEE Communications Letters_ , vol. 24, no. 8, pp. 1724–1728, 2020.
  * [39] A. Barthelme and W. Utschick, “DoA estimation using neural network-based covariance matrix reconstruction,” _IEEE Signal Processing Letters_ , vol. 28, pp. 783–787, 2021.
  * [40] A. M. Elbir, “DeepMUSIC: Multiple signal classification via deep learning,” _IEEE Sensors Letters_ , vol. 4, no. 4, pp. 1–4, 2020.
  * [41] G. Izacard, S. Mohan, and C. Fernandez-Granda, “Data-driven estimation of sinusoid frequencies,” in _Proceedings of the 33rd International Conference on Neural Information Processing Systems (NeurIPS)_ , 2019, pp. 5127–5137.
  * [42] G. K. Papageorgiou, M. Sellathurai, and Y. C. Eldar, “Deep networks for direction-of-arrival estimation in low SNR,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 3714–3729, 2021.
  * [43] Y. Yuan, S. Wu, M. Wu, and N. Yuan, “Unsupervised learning strategy for direction-of-arrival estimation network,” _IEEE Signal Processing Letters_ , vol. 28, pp. 1450–1454, 2021.
  * [44] Y. Yang, N. Zou, E. Lin, F. Suo, and Z. Chen, “A neural network method for nonconvex optimization and its application on parameter retrieval,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 3383–3398, 2021.
  * [45] D. P. Kingma and M. Welling, “Auto-encoding variational Bayes,” _arXiv preprint arXiv:1312.6114_ , 2013.
  * [46] A. Caciularu and D. Burshtein, “Blind channel equalization using variational autoencoders,” in _2018 IEEE International Conference on Communications Workshops (ICC Workshops)_ , 2018, pp. 1–6.
  * [47] A. Caciularu and D. Burshtein, “Unsupervised linear and nonlinear channel equalization and decoding using variational autoencoders,” _IEEE Transactions on Cognitive Communications and Networking_ , vol. 6, no. 3, pp. 1003–1018, 2020.
  * [48] I. Khemakhem, D. Kingma, R. Monti, and A. Hyvarinen, “Variational autoencoders and nonlinear ICA: A unifying framework,” in _International Conference on Artificial Intelligence and Statistics_. PMLR, 2020, pp. 2207–2217.
  * [49] P. Jiang, C.-K. Wen, S. Jin, and G. Y. Li, “Dual CNN-based channel estimation for MIMO-OFDM systems,” _IEEE Transactions on Communications_ , vol. 69, no. 9, pp. 5859–5872, 2021.
  * [50] W.-K. Ma, T.-H. Hsieh, and C.-Y. Chi, “DOA estimation of quasi-stationary signals with less sensors than sources and unknown spatial noise covariance: A khatri–rao subspace approach,” _IEEE Transactions on Signal Processing_ , vol. 58, no. 4, pp. 2168–2180, 2010.
  * [51] W.-K. Ma, T.-H. Hsieh, and C.-Y. Chi, “Doa estimation of quasi-stationary signals via khatri-rao subspace,” in _2009 IEEE International Conference on Acoustics, Speech and Signal Processing_ , 2009, pp. 2165–2168.
  * [52] Z. Yang and L. Xie, “Enhancing sparsity and resolution via reweighted atomic norm minimization,” _IEEE Transactions on Signal Processing_ , vol. 64, no. 4, pp. 995–1006, 2015.
  * [53] H.-H. Chao and L. Vandenberghe, “Semidefinite representations of gauge functions for structured low-rank matrix decomposition,” _SIAM Journal on Optimization_ , vol. 27, no. 3, pp. 1362–1389, 2017.
  * [54] V. Lauinger, F. Buchali, and L. Schmalen, “Blind equalization and channel estimation in coherent optical communications using variational autoencoders,” _IEEE Journal on Selected Areas in Communications_ , vol. 40, no. 9, pp. 2529–2539, 2022.
  * [55] S. Y. Lee, “Gibbs sampler and coordinate ascent variational inference: A set-theoretical review,” _Communications in Statistics-Theory and Methods_ , pp. 1–21, 2021.
  * [56] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2016, pp. 770–778.
  * [57] B. D. Van Veen and K. M. Buckley, “Beamforming: a versatile approach to spatial filtering,” _IEEE ASSP Magazine_ , vol. 5, no. 2, pp. 4–24, 1988\.
  * [58] Q. He, L. Xiao, X. Zhong, and S. Zhou, “Increasing the sum-throughput of cells with a sectorization method for massive MIMO,” _IEEE Communications Letters_ , vol. 18, no. 10, pp. 1827–1830, 2014.
  * [59] K. Kreutz-Delgado, “The complex gradient operator and the cr-calculus,” _arXiv preprint arXiv:0906.4835_ , 2009.