# Parallel and Distributed Hybrid Beamforming for Multicell Millimeter Wave
Full Duplex

Chandan Kumar Sheemar, Symeon Chatzinotas, Dirk Slock and Jorge Querol
Chandan Kumar Sheemar, Symeon Chatzinotas and Jorge Querol are with the SnT
department at the University of Luxembourg
(email:{chandankumar.sheemar,symeon.chatzinotas,jorge.querol}@uni,.lu). Dirk
Slock is with the communication systems department at EURECOM, Sophia
Antipolis, France (email:slock@eurecom.fr).

###### Abstract

This article presents two novel hybrid beamforming (HYBF) designs for a
multicell massive multiple-input-multiple-output (mMIMO) millimeter wave
(mmWave) full duplex (FD) system under limited dynamic range (LDR). Firstly, a
novel centralized HYBF (C-HYBF) scheme based on the minorization-maximization
(MM) method is presented. However, C-HYBF presents many drawbacks such as high
computational complexity, massive communication overhead to transfer complete
channel state information (CSI) to the central node (CN) every channel
coherence time (CCT), and requirement of expensive computational resources for
joint optimization. To overcome these drawbacks, we present a very low
complexity, per-link parallel and distributed HYBF (P$\&$D-HYBF) scheme based
on cooperation. Due to the per-link independent decomposition, it enables each
FD base station (BS) to solve its local sub-problems independently and in
parallel on multiple processors, which leads to significant reduction in the
communication overhead. It requires that each FD BS cooperates by exchanging
information about its beamformers, which allows each FD BS to adapt its
beamformers correctly, leading to a negligible performance loss compared to
C-HYBF. Simulation results show that both designs achieve similar performance
and outperform the fully digital half duplex (HD) system with only a few
radio-frequency (RF) chains.

##  I Introduction

Cellular communication networks are continuously evolving to keep up with the
escalating demands of the newly emerging wireless data services and maintain
seamless connectivity. To accommodate future data growth, the use of ultra-
high-frequency bands such as the millimeter wave (mmWave) band has recently
attracted considerable interest in the research community [1, 2]. In contrast
to the radio frequency (RF) and microwave bands, mmWave offers a $200$ greater
spectrum. However, mmWave communication systems suffer from severe propagation
losses, tackled with a massive number of antennas, which due to shorter
wavelength, can be packed in a tiny area.

In parallel to the development of mmWave half duplex (HD) systems has been the
progression of full duplex (FD) technology which enables the devices to
simultaneously transmit and receive in the same frequency band, which
theoretically doubles the spectral efficiency. Beyond that, FD can be
beneficial to improve security, enable advanced joint communication and
sensing, reduce end-to-end delays, and enable more flexible utilization of the
vast mmWave spectrum [3]. However, FD systems suffer from self-interference
(SI), which can be $90-100$ dB higher than the received signal power. Advanced
SI cancellation (SIC) techniques are crucial to mitigate the SI power and make
FD feasible. Although the promising two-fold gain of the ideal FD systems in
the spectrum efficiency is very much desirable, in practice, the maximum
achievable SIC levels depend on the non-ideal hardware, including power
amplifiers (PAs), analog-to-digital-converters (ADCs), digital-to-analog-
converters, mixers and low noise PAs. The overall effect of the non-ideal
hardware in the transmit and receive RF chains can be captured with the well-
established limited dynamic range (LDR) noise model [4, 5], which allows to
investigate the effective achievable gains in a practical FD systems.

###  I-A Prior Work and Motivation

Hybrid beamforming (HYBF) is a promising signal processing tool to deal with
SI in mmWave and enable FD operation with fewer RF chains than the massive
number of antennas. In [6], the authors presented a novel HYBF design for a
point-to-point massive multiple-input multiple-output (mMIMO) mmWave FD system
to optimize the sum-spectral efficiency while keeping the signal level at the
input of the ADCs under control. In [7], a new algorithm for the design of
constant-amplitude analog precoders and combiners of an FD mmWave single-
stream bidirectional link is presented. In [8], HYBF for an FD mmWave MIMO
interference channel is presented. In [9], a two-timescale HYBF approach for
the mmWave FD relay is presented. A novel robust HYBF under imperfect CSI is
presented in [10]. In [11], HYBF for integrated access and backhaul is
presented. In [12], the authors presented a novel HYBF design for amplify and
forwards mmWave FD relay. In [13], a novel HYBF design for a single-antenna
multi-user mmWave FD system with $1$-bit phase resolution is proposed. In
[14], the achievable gain of a mmWave FD system with one uplink (UL) and one
downlink (DL) user only under the limited receive dynamic range is
investigated. In [15], we presented a practical HYBF design for a single-cell
mMIMO mmWave FD system. In [16], the authors proposed HYBF for a mmWave FD
joint communication and sensing system.

Although the aforementioned papers investigated HYBF for mmWave FD systems in
many different scenario, we remark that the literature does not contribute to
the multicell mmWave FD systems with multi-antenna UL and DL users.
Furthermore, no distributed beamforming design for the FD systems with UL and
DL users, either for sub-$6\leavevmode\nobreak\ $GHz or mmWave, has ever been
presented. To the best of our knowledge, the only exiting solutions for
distributed beamforming in FD are available for sub-$6$ GHz, limited only to
FD relays [17, 18, 19] and point-to-point FD systems [20]. It is to be noted
that FD systems require significant additional channel state information (CSI)
compared to the HD systems, as shown in Figures 2-2. Centralized solutions are
undesirable in general, but especially for FD due to large amount of CSI, as
they require full CSI to be transferred to the central node (CN) every channel
coherence time (CCT), e.g., based on a multi-hop communication, if the CN is
located far from the network. The CN has to optimize all the variables
jointly, which leads to high hardware cost for the computational resources and
large execution times. After the optimization, all the variables need to be
communicated back to the FD BSs via feedback, and finally, each FD BS must
communicate the optimized variables to its users. For any centralized design,
this procedure must be executed at the millisecond scale, which is
prohibitive. Distributed designs [21, 22, 23] are crucial to solve the
aforementioned challenges by allowing each FD BS to optimize its own variables
in a fully distributed fashion based on the local information, which can lead
to significant reduction in the communication overhead [24, 25].

Figure 1: multicell MIMO HD system.

Figure 2: multicell MIMO FD system.

###  I-B Main Contributions

To fill the aforementioned research gaps, in this paper, we present novel HYBF
designs for a multicell mMIMO mmWave FD system with multi-antenna UL and DL
users to maximize the weighted sum rate (WSR). The FD BSs and users are
assumed to have fully connected hybrid and digital processing capabilities,
respectively. The analog stage is assumed to be quantized and non-ideal
hardware is assumed for all the devices, which is modelled with the LDR noise
model [5, 15].

The main contributions of our work are as follows:

  * •

Although centralized solutions are undesirable due to their inherent
challenges at the implementation level, they exhibit superior performance.
Hence, they can serve as a benchmark to compare the performance loss of
distributed schemes [26]. To this end, we present the first-ever novel
centralized HYBF (C-HYBF) scheme for the considered scenario using the
minorization-maximization (MM) [27], relying on alternating optimization.

  * •

MM method has been adopted for over a decade as a promising tool for
beamforming [27]. However, its full potential has been undiscovered so far. By
identifying the hidden structure obtained after applying MM, we show how to
perform per-link independent optimization [24, 23], thus paving the path
towards per-link parallel and distributed (P$\&$D)-HYBF. Our proposed novel
P$\&$D-HYBF design allows decomposing the global optimization problem into
very low complexity per-link independent subproblems, which can be solved
independently on multiple processors by each FD BS in parallel. This
eliminates all the implementation challenges of C-HYBF. Our design is based on
cooperation, which requires the FD BSs to cooperate by exchanging information
about their beamformers with the neighbouring cells.

  * •

P$\&$D-HYBF is non-trivial as the analog beamformers and analog combiners are
shared between the DL and UL users in the same cell, respectively. Moreover,
the beamformers of the DL users are subject to a coupled total sum-power
constraint imposed at the FD BSs, which is also affected by the analog
beamformers. Our presented P$\&$D-HYBF framework also shows how to handle
these coupling constraints to make per-link independent optimization feasible.

  * •

Computational analysis shows that P$\&$D-HYBF exhibits very low complexity,
increasing only linearly with the network size and density, in contrast to the
C-HYBF scheme, which has quadratic dependence on both and requires extremely
high computational complexity. This makes P$\&$D-HYBF highly scalable, and
thus it has the potential to enable the deployment of low-cost distributed
computational processors over a large multicell FD network.

  * •

Simulation results show that P$\&$D-HYBF exhibits negligible performance loss
due to distributed implementation. Both designs achieve similar performance
gains and outperform the conventional fully digital HD system with only a few
RF chains and with very low phase resolution at the analog stage of the FD
BSs. It is also shown that P$\&$D-HYBF, due to low complexity, requires
considerably less execution time compared to C-HYBF.

_Paper Organization:_ The rest of the paper is organized as follows. We first
present the system model and problem formulation in Section II. Then, the MM
optimization method and C-HYBF design are presented in Sections III and IV,
respectively. Section V presents the P$\&$D-HYBF design based on cooperation.
Finally, Sections VI and VII present the simulation results and conclusions,
respectively.

_Mathematical Notations:_ Boldface lower and upper case characters denote
vectors and matrices, respectively.
$\mathbb{E}\\{\cdot\\},\mbox{Tr}\\{\cdot\\},(\cdot)^{H}$, $\otimes$,
$\mathbf{I}$, and $\mathbf{D}_{d}$ denote expectation, trace, conjugate
transpose, Kronecker product, identity matrix and the $d$ generalized dominant
eigenvectors (GDEs) selection matrix, respectively. Vector of zeros of size
$M$ is denoted as $\mathbf{0}_{M\times 1}$, $\mbox{vec}(\mathbf{X})$ stacks
the column of $\mathbf{X}$ into $\mathbf{x}$, $\mbox{unvec}(\mathbf{x})$
reshapes $\mathbf{x}$ into $\mathbf{X}$, and $\angle\mathbf{X}$ returns the
phasors of matrix $\mathbf{X}$. $\mbox{Cov}(\cdot)$ and diag$(\cdot)$ denote
the covariance and diagonal matrices, respectively, and $svd(\mathbf{X})$
returns the singular value decomposition (SVD) of $\mathbf{X}$. The element of
matrix $\mathbf{X}$ at the m-th row and n-th column is denoted as
$\mathbf{X}(m,n)$.

##  II System Model

Let $\mathcal{B}=\\{1,....,B\\}$ denote the set containing the indices of $B$
FD BSs serving in $B$ cells. Let $\mathcal{D}_{b}=\\{1,...,D_{b}\\}$ and
$\mathcal{U}_{b}=\\{1,...,U_{b}\\}$ denote the sets containing the indices of
$D_{b}$ DL and $U_{b}$ UL multi-antenna HD users communicating with BS
$b\in\mathcal{B}$. The DL user $j_{b}\in\mathcal{D}_{b}$ and UL user
$k_{b}\in\mathcal{U}_{b}$ are assumed to have $N_{j_{b}}$ receive and
$M_{k_{b}}$ transmit antennas, respectively. The FD BS $b\in\mathcal{B}$ is
assumed to have $M_{b}^{RF}$ and $N_{b}^{RF}$ transmit and receive RF chains,
respectively, and $M_{b}$ and $N_{b}$ transmit and receive antennas,
respectively. We denote with
$\mathbf{V}_{j_{b}}\in\mathbb{C}^{M_{b}^{RF}\times d_{j_{b}}}$ and
$\mathbf{U}_{k_{b}}\in\mathbb{C}^{M_{k_{b}}\times d_{k_{b}}}$ the digital
beamformers for the white unitary variance data streams
$\mathbf{s}_{j_{b}}\in\mathbb{C}^{d_{j_{b}}\times 1}$ and
$\mathbf{s}_{k_{b}}\in\mathbb{C}^{d_{k_{b}}\times 1}$ transmitted for DL user
$j_{b}\in\mathcal{D}_{b}$ and from UL user $k_{b}\in\mathcal{U}_{b}$,
respectively. Let $\mathbf{W}_{b}\in\mathbb{C}^{M_{b}\times M_{b}^{RF}}$ and
$\mathbf{F}_{b}\in\mathbb{C}^{N_{b}^{RF}\times N_{b}}$ denote the fully
connected analog beamformer and analog combiner for FD BS $b\in\mathcal{B}$,
respectively. Let
$\mathcal{P}_{b}=\\{1,\;e^{i2\pi/n_{b}},\;e^{i4\pi/n_{b}},...,e^{i2\pi
n_{b}-1/n_{b}}\\}$ denote the set of $n_{b}$ possible discrete values that the
unit-modulus elements of $\mathbf{W}_{b}$ and $\mathbf{F}_{b}$ can assume. Let
$\mathbb{Q}_{b}(\cdot)$ denote the quantizer function to quantize the infinite
resolution unit-modulus elements of $\mathbf{W}_{b}$ and $\mathbf{F}_{b}$ such
that $\mathbb{Q}_{b}(\angle\mathbf{W}_{b}(m,n))$ and
$\mathbb{Q}_{b}(\angle\mathbf{F}_{b}(m,n)\in\mathcal{P}_{b},\forall m,n$.

The users and the FD BSs are assumed to suffer from the LDR noise due to non-
ideal hardware. It is denoted as $\mathbf{c}_{k_{b}}$ and $\mathbf{e}_{j_{b}}$
for the UL user $k_{b}\in\mathcal{U}_{b}$ and DL user
$j_{b}\in\mathcal{D}_{b}$, respectively, modelled as [5]

$\mathbf{c}_{k_{b}}\sim\mathcal{CN}(\mathbf{0}_{M_{k_{b}}\times
1},k_{k_{b}}\;\mbox{diag}(\mathbf{U}_{k_{b}}\mathbf{U}_{k_{b}}^{H})),\quad\mathbf{e}_{j_{b}}\sim\mathcal{CN}(\mathbf{0}_{N_{j_{b}}\times
1},\beta_{j_{b}}\;\mbox{diag}(\mathbf{\Phi}_{j_{b}})),$ (1)

where $k_{k_{b}}\ll 1$, $\beta_{j_{b}}\ll
1,\mathbf{\Phi}_{j_{b}}=\mbox{Cov}(\mathbf{r}_{j_{b}})$ and
$\mathbf{r}_{j_{b}}$ denotes the undistorted received signal for DL user
$j_{b}\in\mathcal{D}_{b}$. Let $\mathbf{c}_{b}$ and $\mathbf{e}_{b}$ denote
the transmit and receive LDR noise for a hybrid FD BS $b\in\mathcal{B}$,
respectively, which can be modelled as [15]

$\mathbf{c}_{b}\sim\mathcal{CN}(\mathbf{0}_{M_{b}\times
1},k_{b}\;\mbox{diag}(\sum_{n_{b}\in\mathcal{D}_{b}}\mathbf{W}_{b}\mathbf{V}_{n_{b}}\mathbf{V}_{n_{b}}^{H}\mathbf{W}_{b}^{H})),\quad\mathbf{e}_{b}\sim\mathcal{CN}({\mathbf{0}_{N_{b}^{RF}\times
1},\beta_{b}\;\mbox{diag}(\mathbf{\Phi}_{b})}),$ (2)

with $k_{b}\ll 1$, $\beta_{b}\ll
1,\mathbf{\Phi}_{b}=\mbox{Cov}(\mathbf{r}_{b})$ and $\mathbf{r}_{b}$ denotes
the undistorted received signal by FD BS $b\in\mathcal{B}$ after the analog
combiner $\mathbf{F}_{b}$. Let $\mathbf{n}_{b}$ and $\mathbf{n}_{j_{b}}$
denote the thermal noise for FD BS $b$ and DL user $j_{b}$, respectively,
modelled as

$\mathbf{n}_{b}\sim\mathcal{CN}(\mathbf{0}_{N_{b}\times
1},\sigma_{b}^{2}\mathbf{I}),\quad\quad\mathbf{n}_{j_{b}}\sim\mathcal{CN}(\mathbf{0}_{N_{j_{b}}\times
1},\sigma_{j_{b}}^{2}\mathbf{I}),$ (3)

where $\sigma_{b}^{2}$ and $\sigma_{j_{b}}^{2}$ denote the noise variances.

###  II-A Channel Modelling

We assume perfect CSI and let
$\mathbf{H}_{j_{b}}\in\mathbb{C}^{N_{j_{b}}\times M_{b}}$ and
$\mathbf{H}_{k_{b}}\in\mathbb{C}^{N_{b}\times M_{k_{b}}}$ denote the direct
channels responses between the DL user $j_{b}\in\mathcal{D}_{b}$ and UL user
$k_{b}\in\mathcal{U}_{b}$, respectively, and their serving FD BS
$b\in\mathcal{B}$. Let $\mathbf{H}_{j_{b},k_{b}}\in\mathbb{C}^{N_{j_{b}}\times
M_{k_{b}}}$ and $\mathbf{H}_{j_{b},k_{c}}\in\mathbb{C}^{N_{j_{b}}\times
M_{k_{c}}}$ denote the in-cell UL cross-interference (CI), i.e., generated
from opposite transmission directions, channel response between the DL user
$j_{b}\in\mathcal{D}_{b}$ and UL user $k_{b}\in\mathcal{U}_{b}$ and the out-
cell UL CI channel response between the DL user $j_{b}\in\mathcal{D}_{b}$ and
UL user $k_{c}\in\mathcal{U}_{c}$, respectively, with $b\neq c$. Let
$\mathbf{H}_{j_{b},c}\in\mathbb{C}^{N_{j_{b}}\times M_{c}}$ and
$\mathbf{H}_{b,k_{c}}\in\mathbb{C}^{N_{b}\times M_{k_{c}}}$ denote the
interference channels responses from FD BS $c\in\mathcal{B}$ to DL user
$j_{b}\in\mathcal{D}_{b}$ and from UL user $k_{c}\in\mathcal{U}_{c}$ to FD BS
$b$, respectively, with $c\neq b$. Let
$\mathbf{H}_{b,c}\in\mathbb{C}^{N_{b}\times M_{c}}$ and
$\mathbf{H}_{b,b}\in\mathbb{C}^{N_{b}\times M_{b}}$ denote the DL CI channel
response from FD BS $c\in\mathcal{B}$ to FD BS $b\in\mathcal{B}$, with $c\neq
b$, and the SI channel response for FD BS $b\in\mathcal{B}$, respectively.

In mmWave, channel response $\mathbf{H}_{k_{b}}$ can be modelled as

$\displaystyle\mathbf{H}_{k_{b}}=\sqrt{\frac{{1}}{{N_{k_{b}}}}}\sum_{n=1}^{N_{k_{b}}^{p}}$
$\displaystyle\alpha_{k_{b}}^{n}\mathbf{a}_{r}^{b}\;(\phi_{k_{b}}^{n})\;\mathbf{a}_{t}^{k_{b}}(\theta_{k_{b}}^{n})^{T}.$
(4)

The scalars $N_{k_{b}}^{p}$ and $\alpha_{k_{b}}^{n}$ denote the number of
paths and a complex Gaussian random variable with amplitudes and phases
distributed according to the Rayleigh and uniform distribution, respectively.
The vectors $\mathbf{a}_{r}^{b}(\phi_{k_{b}}^{n})$ and
$\mathbf{a}_{t}^{k_{b}}(\theta_{k_{b}}^{n})^{T}$ denote the receive and
transmit antenna array response for FD BS $b\in\mathcal{B}$ and UL user
$k_{b}\in\mathcal{U}_{b}$, respectively, with the angle of arrival (AoA)
$\phi_{k_{b}}^{n}$ and angle of departure (AoD) $\theta_{k_{b}}^{n}$,
respectively. The channel responses
$\mathbf{H}_{j_{b}},\mathbf{H}_{j_{b},k_{b}},\mathbf{H}_{j_{b},k_{c}},\mathbf{H}_{j_{b},c}$
and $\mathbf{H}_{b,k_{c}}$ can be modelled similarly as
$\eqref{channel_Model}$ and the SI channel
$\mathbf{H}_{b,b}\in\mathbb{C}^{N_{b}\times M_{b}}$ can be modelled as [15]

Table I: Main Notations

$M_{b}^{RF}$ |  Transmit RF chains for FD BS $b$
---|---
$N_{b}^{RF}$ |  Receive RF chains for FD BS $b$
$M_{b}$ |  Transmit antennas at the FD BS $b$
$N_{b}$ |  Receive antennas at the FD BS $b$
$M_{k_{b}}$ |  Transmit antennas for the UL user $k_{b}$
$N_{j_{b}}$ |  Receive antennas for the DL user $j_{b}$
$\mathbf{U}_{k_{b}}$ |  Digital beamformer for UL user $k_{b}$
$\mathbf{V}_{j_{b}}$ |  Digital beamformer for DL user $j_{b}$
$\mathbf{W}_{b}$ |  Analog beamformer for BS $b$
$\mathbf{F}_{b}$ |  Analog combiner at BS b
$\mathbf{H}_{j_{b}}$ |  Direct channel for DL user $j_{b}$
$\mathbf{H}_{k_{b}}$ |  Direct channel for UL user $k_{b}$
$\mathbf{H}_{j_{b},k_{b}}$ |  CI channel from $k_{b}$ to $j_{b}$
$\mathbf{H}_{j_{b},k_{c}}$ |  CI channel from $k_{c}$ to $j_{b}$
$\mathbf{H}_{j_{b},c}$ |  Interference channel from BS $c$ to $j_{b}$
$\mathbf{H}_{b,k_{c}}$ |  Interference channel from user $k_{c}$ to the BS $b$
$\mathbf{H}_{b,c}$ |  BS-to-BS interference channel from BS $c$ to the BS $b$
$\mathbf{H}_{b,b}$ |  SI channel for BS $b$

$\mathbf{H}_{b}=\sqrt{\frac{\kappa_{b}}{\kappa_{b}+1}}\mathbf{H}_{b}^{L}+\sqrt{\frac{1}{\kappa_{b}+1}}\mathbf{H}_{b}^{R},$
(5)

where $\mathbf{H}_{b}^{R}$ denotes the reflected components which can be
modelled as (4) and $\mathbf{H}_{b}^{L}$ is the line of sight (LoS) channel,
modelled as

$\mathbf{H}_{b}^{L}(m,n)=\frac{\rho_{b}}{r_{m,n}}e^{-j2\pi\frac{r_{m,n}}{\lambda}}.$
(6)

The scalars $\kappa_{b}$, $\rho_{b}$, $r_{m,n}$ and $\lambda$ denote the
Rician factor, the power normalization constant to assure
$\mathbb{E}(||\mathbf{H}_{b}^{L}(m,n)||_{F}^{2})=M_{b}N_{b}$, the distance
between $m$-th receive and $n$-th transmit antenna and the wavelength,
respectively. The aforementioned main notations are summarized in Table I.

###  II-B Problem Formulation

Let $\mathbf{y}_{j_{b}}$ and $\mathbf{y}_{k_{b}}$ denote the signals received
by the DL user $j_{b}\in\mathcal{D}_{b}$ and by the FD BS $b\in\mathcal{B}$
from UL user $k_{b}\in\mathcal{U}_{b}$ after the analog combiner
$\mathbf{F}_{b}$, respectively, given as

$\displaystyle\mathbf{y}_{j_{b}}=\mathbf{H}_{j_{b}}($
$\displaystyle\sum_{n_{b}\in\mathcal{D}_{b}}\mathbf{W}_{b}\mathbf{V}_{n_{b}}\mathbf{s}_{n_{b}}+\mathbf{c}_{b})+\mathbf{e}_{j_{b}}+\mathbf{n}_{j_{b}}+\sum_{k_{b}\in\mathcal{U}_{b}}\mathbf{H}_{j_{b},k_{b}}(\mathbf{U}_{k_{b}}\mathbf{s}_{k_{b}}+\mathbf{c}_{k_{b}})$
(7) $\displaystyle+\sum_{\begin{subarray}{c}c\in\mathcal{B},\\\ c\neq
b\end{subarray}}[\mathbf{H}_{j_{b},c}(\sum_{n_{c}\in\mathcal{D}_{c}}\mathbf{W}_{c}\mathbf{V}_{n_{c}}\mathbf{s}_{n_{c}}+\mathbf{c}_{c})+\sum_{k_{c}\in\mathcal{U}_{c}}\mathbf{H}_{j_{b},k_{c}}(\mathbf{U}_{k_{c}}\mathbf{s}_{k_{c}}+\mathbf{c}_{k_{c}})],$
$\displaystyle\mathbf{y}_{k_{b}}={\mathbf{F}_{b}}^{H}$
$\displaystyle(\sum_{k_{b}\in\mathcal{U}_{b}}\mathbf{H}_{k_{b}}(\mathbf{U}_{k_{b}}\mathbf{s}_{k_{b}}+\mathbf{c}_{k_{b}})+\mathbf{n}_{b}+\mathbf{H}_{b,b}(\sum_{j_{b}\in\mathcal{D}_{b}}\mathbf{W}_{b}\mathbf{V}_{j_{b}}\mathbf{s}_{j_{b}}+\mathbf{c}_{b})$
(8) $\displaystyle+\sum_{\begin{subarray}{c}c\in\mathcal{B},\\\ c\neq
b\end{subarray}}[\mathbf{H}_{b,c}(\sum_{j_{c}\in\mathcal{D}_{c}}\mathbf{W}_{c}\mathbf{V}_{j_{c}}\mathbf{s}_{j_{c}}+\mathbf{c}_{c})+\sum_{k_{c}\in\mathcal{U}_{c}}\mathbf{H}_{b,k_{c}}(\mathbf{U}_{k_{c}}\mathbf{s}_{k_{c}}+\mathbf{c}_{k_{c}})]+\mathbf{e}_{b}.$

Let $\overline{k_{b}}$, $\overline{j_{b}}$ and $\overline{b}$ denote the
indices in the sets $\mathcal{U}_{b}$, $\mathcal{D}_{b}$ and $\mathcal{B}$
without the elements $k_{b}$, $j_{b}$ and $b$, respectively. Let
$\mathbf{T}_{k_{b}}\triangleq\mathbf{U}_{k_{b}}\mathbf{U}_{k_{b}}^{H}$ and
$\mathbf{Q}_{j_{b}}\triangleq\mathbf{W}_{b}\mathbf{V}_{j_{b}}\mathbf{V}_{j_{b}}^{H}\mathbf{W}_{b}^{H}$
denote the transmit covariance matrices of UL user $k_{b}\in\mathcal{U}_{b}$
and of FD BS $b\in\mathcal{B}$ intended for its DL user
$j_{b}\in\mathcal{D}_{b}$, respectively. Let ($\mathbf{R}_{k_{b}}$)
$\mathbf{R}_{\overline{k_{b}}}$ and ($\mathbf{R}_{j_{b}}$)
$\mathbf{R}_{\overline{j_{b}}}$ denote the (signal plus) interference plus
noise covariance matrices received by the FD BS $b\in\mathcal{B}$ from UL user
$k_{b}\in\mathcal{U}_{b}$ and by the DL user $j_{b}\in\mathcal{D}_{b}$,
respectively. The matrices $\mathbf{R}_{k_{b}}$ and $\mathbf{R}_{j_{b}}$ can
be written as (9) and $\mathbf{R}_{\overline{k_{b}}}$ and
$\mathbf{R}_{\overline{j_{b}}}$ can be obtained as
$\mathbf{R}_{\overline{k_{b}}}=\mathbf{R}_{k_{b}}-\mathbf{H}_{j_{b}}\mathbf{Q}_{j_{b}}\mathbf{H}_{j_{b}}^{H}$
and
$\mathbf{R}_{\overline{j_{b}}}=\mathbf{R}_{j_{b}}-{\mathbf{F}_{b}}^{H}\mathbf{H}_{k_{b}}\mathbf{T}_{k_{b}}\mathbf{H}_{k_{b}}^{H}{\mathbf{F}_{b}}$,
respectively.

$\displaystyle\mathbf{R}_{j_{b}}=$
$\displaystyle\mathbf{H}_{j_{b}}\mathbf{Q}_{j_{b}}\mathbf{H}_{j_{b}}^{H}+\mathbf{H}_{j_{b}}(\sum_{\begin{subarray}{c}n_{b}\in\mathcal{D}_{b}\\\
n_{b}\neq
j_{b}\end{subarray}}\mathbf{Q}_{n_{b}})\mathbf{H}_{j_{b}}^{H}+\mathbf{H}_{j_{b}}k_{b}\mbox{diag}(\sum_{n_{b}\in\mathcal{D}_{b}}\mathbf{Q}_{n_{b}})\mathbf{H}_{j_{b}}^{H}+\sum_{k_{b}\in\mathcal{U}_{b}}\mathbf{H}_{j_{b},k_{b}}(\mathbf{T}_{k_{b}}+k_{k_{b}}\mbox{diag}(\mathbf{T}_{k_{b}}))$
(9a)
$\displaystyle\mathbf{H}_{j_{b},k_{b}}^{H}+\sum_{\begin{subarray}{c}c\in\mathcal{B}\\\
c\neq
b\end{subarray}}\mathbf{H}_{j_{b},c}(\sum_{n_{c}\in\mathcal{D}_{c}}\mathbf{Q}_{n_{c}}+k_{c}\mbox{diag}(\mathbf{Q}_{n_{c}}))\mathbf{H}_{j_{b},c}^{H}+\sum_{\begin{subarray}{c}c\in\mathcal{B}\\\
c\neq
b\end{subarray}}\sum_{k_{c}\in\mathcal{U}_{c}}\mathbf{H}_{j_{b},k_{c}}(\mathbf{T}_{k_{c}}+k_{k_{c}}\mbox{diag}(\mathbf{T}_{k_{c}}))\mathbf{H}_{j_{b},k_{c}}^{H}$
$\displaystyle+\beta_{j_{b}}\mbox{diag}(\Phi_{j_{b}})+\sigma_{j_{b}}^{2}\mathbf{I}_{N_{j_{b}}},$
$\displaystyle\mathbf{R}_{k_{b}}=$
$\displaystyle{\mathbf{F}_{b}}^{H}(\mathbf{H}_{k_{b}}\mathbf{T}_{k_{b}}\mathbf{H}_{k_{b}}^{H}+\sum_{\begin{subarray}{c}m_{b}\in\mathcal{U}_{b}\\\
m_{b}\neq
k_{b}\end{subarray}}\mathbf{H}_{m_{b}}\mathbf{T}_{m_{b}}\mathbf{H}_{m_{b}}^{H}+\sum_{m_{b}\in\mathcal{U}_{b}}k_{m_{b}}\mathbf{H}_{m_{b}}\mbox{diag}(\mathbf{T}_{m_{b}})\mathbf{H}_{m_{b}}^{H}+\mathbf{H}_{b,b}\sum_{j_{b}\in\mathcal{D}_{b}}(\mathbf{Q}_{j_{b}}$
(9b)
$\displaystyle+k_{b}\mbox{diag}(\mathbf{Q}_{j_{b}}))\mathbf{H}_{b,b}^{H}+\sum_{\begin{subarray}{c}c\in\mathcal{B}\\\
c\neq
b\end{subarray}}\mathbf{H}_{b,c}\sum_{j_{c}\in\mathcal{D}_{c}}(\mathbf{Q}_{j_{c}}+k_{c}\mbox{diag}(\mathbf{Q}_{j_{c}}))\mathbf{H}_{b,c}^{H}+\sum_{\begin{subarray}{c}c\in\mathcal{B}\\\
c\neq
b\end{subarray}}\sum_{k_{c}\in\mathcal{U}_{c}}\mathbf{H}_{b,k_{c}}(\mathbf{T}_{k_{c}}$
$\displaystyle+k_{k_{c}}\mbox{diag}(\mathbf{T}_{k_{c}}))\mathbf{H}_{b,k_{c}}^{H}+\sigma_{b}^{2}\mathbf{I}_{n_{b}}){\mathbf{F}_{b}}+\beta_{b}\mbox{diag}(\mathbf{\Phi}_{b}),$

The WSR maximization problem for HYBF in the considered scenario with $J_{b}$
DL and $U_{b}$ UL multi-antenna users $\forall b\in\mathcal{B}$, under the
joint sum-power, unit-modulus and discrete phase-shifters constraints can be
stated as

$\displaystyle\underset{\begin{subarray}{c}\mathbf{U},\mathbf{V},\\\
\mathbf{W}_{b}\mathbf{F}^{RF}\end{subarray}}{\max}\sum_{b\in\mathcal{B}}\Big{[}\sum_{k_{b}\in\mathcal{U}_{b}}w_{k_{b}}\mbox{lndet}(\mathbf{R}_{\overline{k_{b}}}^{-1}\mathbf{R}_{k_{b}})+\sum_{j_{b}\in\mathcal{D}_{b}}w_{j_{b}}\mbox{lndet}(\mathbf{R}_{\overline{j_{b}}}^{-1}\mathbf{R}_{j_{b}})\Big{]}$
(10a) $\text{s.t.}\quad\mbox{Tr}(\mathbf{U}_{k_{b}}\mathbf{U}_{k_{b}}^{H})\leq
p_{k_{b}},\quad\forall k_{b}\in\mathcal{U}_{b},$ (10b)
$\mbox{Tr}(\sum_{j\in\mathcal{D}_{b}}\mathbf{W}_{b}\mathbf{V}_{j}\mathbf{V}_{j}^{H}\mathbf{W}_{b}^{H})\leq
p_{b},\quad\forall b\in\mathcal{B},$ (10c)
$\mathbf{W}_{b}(m,n)\;\&\;\mathbf{F}_{b}(m,n)\in\mathcal{P}_{b},\quad\;\forall
m,n.$ (10d)

The scalars $w_{k_{b}}$ and $w_{j_{b}}$ denote rate weights for UL user
$k_{b}\in\mathcal{U}_{b}$ and DL user $j_{b}\in\mathcal{D}_{b}$, respectively,
and the scalars $p_{k_{b}}$ and $p_{b}$ denote the sum-power constraint for UL
user $k_{b}\in\mathcal{U}_{b}$ and FD BS $b\in\mathcal{B}$, respectively. The
collections of digital beamformers in UL and DL are denoted as $\mathbf{U}$
and $\mathbf{V}$, respectively, and the collections of analog beamformers and
combiners are denoted as $\mathbf{W}_{b}$ and $\mathbf{F}_{b}$, respectively.

##  III Minorization-Maximization

Problem (10a) is non-concave in the transmit covariance matrices
$\mathbf{T}_{k_{b}}$ and $\mathbf{Q}_{j_{b}}$ due to the interference terms,
and finding its global optimum is very challenging. To find its sub-optimal
solution based on alternating optimization, we leverage the minorization-
maximization (MM) method [27], which allows reformulating (10a) with its
minorizer using the difference-of-convex (DC) functions [27].

Let $\mbox{WR}_{k_{b}}^{UL}$ and $\mbox{WR}_{j_{b}}^{DL}$ denote the weighted
rate (WR) of users $k_{b}\in\mathcal{U}_{b}$ and $j_{b}\in\mathcal{D}_{b}$,
respectively, and let $\mbox{WSR}_{\overline{k_{b}}}^{UL}$ and
$\mbox{WSR}_{\overline{j_{b}}}^{DL}$ denote the WSR of users outside the cell
$b$ in UL and DL, respectively. The dependence of the global WSR in (10a) on
the aforementioned terms can be highlighted as

$\displaystyle\mbox{WSR}=\mbox{WR}_{k_{b}}^{UL}+\mbox{WSR}_{\overline{k_{b}}}^{UL}+\mbox{WR}_{j_{b}}^{DL}+\mbox{WSR}_{\overline{j_{b}}}^{DL}+\mbox{WSR}_{\overline{b}}^{UL}+\mbox{WSR}_{\overline{b}}^{DL}$
(11)

in which the WSR in UL and DL for FD BS $b\in\mathcal{B}$ is given as
$\mbox{WSR}_{b}^{UL}=\mbox{WR}_{k_{b}}^{UL}+\mbox{WSR}_{\overline{k_{b}}}^{UL}$
and
$\mbox{WSR}_{b}^{DL}=\mbox{WR}_{j_{b}}^{DL}+\mbox{WSR}_{\overline{j_{b}}}^{DL}$,
respectively. Considering the dependence of the transmit covariance matrices
on the global WSR, only $\mbox{WR}_{k_{b}}^{UL}$ is concave in
$\mathbf{T}_{k_{b}}$ and $\mbox{WSR}_{\overline{k_{b}}}^{UL}$,
$\mbox{WSR}_{b}^{DL}$, $\mbox{WSR}_{\overline{b}}^{UL}$ and
$\mbox{WSR}_{\overline{b}}^{DL}$ are non concave in $\mathbf{T}_{k_{b}}$ due
to interference. Similarly, only $\mbox{WR}_{j_{b}}^{DL}$ is concave in
$\mathbf{Q}_{j_{b}}$ and $\mbox{WSR}_{\overline{j_{b}}}^{DL}$,
$\mbox{WSR}_{b}^{UL}$,$\mbox{WSR}_{\overline{b}}^{UL}$,$\mbox{WSR}_{\overline{b}}^{DL}$
are non concave in $\mathbf{Q}_{j_{b}}$. As a linear function is
simultaneously convex and concave, DC programming introduces the first order
Taylor series expansion of $\mbox{WSR}_{\overline{k_{b}}}^{UL}$,
$\mbox{WSR}_{b}^{DL}$, $\mbox{WSR}_{\overline{b}}^{UL}$ and
$\mbox{WSR}_{\overline{b}}^{DL}$ in $\mathbf{T}_{k_{b}}$, around
$\hat{\mathbf{T}}_{k_{b}}$ (i.e. around all $\mathbf{T}_{k_{b}}$), and for
$\mbox{WSR}_{\overline{j_{b}}}^{DL}$,
$\mbox{WSR}_{b}^{UL}$,$\mbox{WSR}_{\overline{b}}^{UL}$,$\mbox{WSR}_{\overline{b}}^{DL}$
around $\hat{\mathbf{Q}}_{j_{b}}$ (i.e. around all $\mathbf{Q}_{j_{b}}$). Let
$\hat{\mathbf{T}}$ and $\hat{\mathbf{Q}}$ denote the sets containing all such
$\hat{\mathbf{T}}_{k_{b}}$ and $\hat{\mathbf{Q}}_{j_{b}}$, respectively. The
tangent expressions with respect to the non-concave terms for
$\mathbf{T}_{k_{b}}$ can be written by computing the following gradients

$\hat{\mathbf{G}}_{\overline{k_{b}},b}^{UL}=-\frac{\partial\mbox{WSR}_{\overline{k_{b}}}^{UL}}{\partial\mathbf{T}_{k_{b}}}|_{\hat{\mathbf{T}},\hat{\mathbf{Q}}},\;\hat{\mathbf{G}}_{{k_{b}},b}^{DL}=-\frac{\partial\mbox{WSR}_{b}^{DL}}{\partial\mathbf{T}_{k_{b}}}|_{\hat{\mathbf{T}},\hat{\mathbf{Q}}},\hat{\mathbf{G}}_{k_{b},\overline{b}}^{UL}=-\frac{\partial\mbox{WSR}_{\overline{b}}^{UL}}{\partial\mathbf{T}_{k_{b}}}|_{\hat{\mathbf{T}},\hat{\mathbf{Q}}},\;\hat{\mathbf{G}}_{{k_{b}},\overline{b}}^{DL}=-\frac{\partial\mbox{WSR}_{\overline{b}}^{DL}}{\partial\mathbf{T}_{k_{b}}}|_{\hat{\mathbf{T}},\hat{\mathbf{Q}}},$
(12a)

which allow to write the following minorizers
$\underline{\mbox{WSR}}_{\overline{k_{b}}}^{UL},$
$\underline{\mbox{WSR}}_{b}^{DL},$
$\underline{\mbox{WSR}}_{\overline{b}}^{UL}$ and
$\underline{\mbox{WSR}}_{\overline{b}}^{DL}$ with respect to
$\mathbf{T}_{k_{b}}$. Similarly, for the transmit covariance matrix
$\mathbf{Q}_{j_{b}}$, we have the gradients

$\hat{\mathbf{G}}_{j_{b},b}^{UL}=-\frac{\partial\mbox{WSR}_{b}^{UL}}{\partial\mathbf{Q}_{j_{b}}}|_{\hat{\mathbf{T}},\hat{\mathbf{Q}}},\;\hat{\mathbf{G}}_{\overline{j_{b}},b}^{DL}=-\frac{\partial\mbox{WSR}_{\overline{j_{b}}}^{DL}}{\partial\mathbf{Q}_{j_{b}}}|_{\hat{\mathbf{T}},\hat{\mathbf{Q}}},\hat{\mathbf{G}}_{j_{b},\overline{b}}^{UL}=-\frac{\partial\mbox{WSR}_{\overline{b}}^{UL}}{\partial\mathbf{Q}_{j_{b}}}|_{\hat{\mathbf{T}},\hat{\mathbf{Q}}},\;\hat{\mathbf{G}}_{{j_{b}},\overline{b}}^{DL}=-\frac{\partial\mbox{WSR}_{\overline{b}}^{DL}}{\partial\mathbf{Q}_{j_{b}}}|_{\hat{\mathbf{T}},\hat{\mathbf{Q}}},$
(13a)

which allow to write the minorizers $\underline{\mbox{WSR}}_{b}^{UL},$
$\underline{\mbox{WSR}}_{\overline{j_{b}}}^{DL},$
$\underline{\mbox{WSR}}_{\overline{b}}^{UL}$ and
$\underline{\mbox{WSR}}_{\overline{b}}^{DL}$ with respect to
$\mathbf{Q}_{j_{b}}$. The gradients (12) and (13) can be computed by applying
the matrix differentiation properties, and the result Lemma 3 [15], and they
are reported in Table II. The tangent expressions constitute a touching lower
bound for the original WSR cost function. Hence, the DC functions constitute a
MM approach, regardless of the restatement of the transmit covariance matrices
$\mathbf{T}_{k_{b}}$ and $\mathbf{Q}_{j_{b}}$ as a function of the
beamformers.

Table II: Gradients expressions to construct the minorized WSR cost function.

$\hat{\mathbf{G}}_{\overline{k_{b}},b}^{UL}$ |  $\sum_{\begin{subarray}{c}m_{b}\in\mathcal{U}_{b}\\\ m_{b}\neq k_{b}\end{subarray}}w_{m_{b}}[\mathbf{H}_{m_{b}}^{H}\mathbf{F}_{b}(\mathbf{R}_{\overline{m_{b}}}^{-1}-\mathbf{R}_{m_{b}}^{-1}+\beta_{b}\mbox{diag}(\mathbf{R}_{\overline{m_{b}}}^{-1}-\mathbf{R}_{m_{b}}^{-1})){\mathbf{F}_{b}}^{H}\mathbf{H}_{m_{b}}+k_{m_{b}}\mbox{diag}(\mathbf{H}_{m_{b}}^{H}\mathbf{F}_{b}(\mathbf{R}_{\overline{m_{b}}}^{-1}-\mathbf{R}_{m_{b}}^{-1}){\mathbf{F}_{b}}^{H}\mathbf{H}_{m_{b}})].$
---|---
$\hat{\mathbf{G}}_{{k_{b}},b}^{DL}$ |  $\sum_{\begin{subarray}{c}j_{b}\in\mathcal{D}_{b}\end{subarray}}w_{j_{b}}[\mathbf{H}_{j_{b},k_{b}}^{H}(\mathbf{R}_{\overline{j_{b}}}^{-1}-\mathbf{R}_{j_{b}}^{-1}+\beta_{j_{b}}\mbox{diag}(\mathbf{R}_{\overline{j_{b}}}^{-1}-\mathbf{R}_{j_{b}}^{-1}))\mathbf{H}_{j_{b},k_{c}}+k_{k_{b}}\mbox{diag}(\mathbf{H}_{j_{b},k_{b}}^{H}(\mathbf{R}_{\overline{j_{b}}}^{-1}-\mathbf{R}_{j_{b}}^{-1})\mathbf{H}_{j_{b},k_{b}})]$
$\hat{\mathbf{G}}_{{k_{b}},\overline{b}}^{UL}$ |  $\sum_{\begin{subarray}{c}c\in\mathcal{B}\\\ c\neq b\end{subarray}}\sum_{\begin{subarray}{c}k_{c}\in\mathcal{U}_{c}\end{subarray}}w_{k_{c}}[\mathbf{H}_{c,k_{b}}^{H}{\mathbf{F}_{c}}(\mathbf{R}_{\overline{k_{c}}}^{-1}-\mathbf{R}_{k_{c}}^{-1}+\beta_{c}\mbox{diag}(\mathbf{R}_{\overline{k_{c}}}^{-1}-\mathbf{R}_{k_{c}}^{-1})){\mathbf{F}_{c}}^{H}\mathbf{H}_{c,k_{b}}+k_{k_{b}}\mbox{diag}(\mathbf{H}_{c,k_{b}}^{H}\mathbf{F}_{c}(\mathbf{R}_{\overline{k_{c}}}^{-1}-\mathbf{R}_{k_{c}}^{-1}){\mathbf{F}_{c}}^{H}\mathbf{H}_{c,k_{b}})]$
$\hat{\mathbf{G}}_{{k_{b}},\overline{b}}^{DL}$ |  $\sum_{\begin{subarray}{c}c\in\mathcal{B}\\\ c\neq b\end{subarray}}\sum_{\begin{subarray}{c}j_{c}\in\mathcal{D}_{c}\end{subarray}}w_{j_{c}}[\mathbf{H}_{j_{c},k_{b}}^{H}(\mathbf{R}_{\overline{j_{c}}}^{-1}-\mathbf{R}_{j_{c}}^{-1}+\beta_{j_{c}}\mbox{diag}(\mathbf{R}_{\overline{j_{c}}}^{-1}-\mathbf{R}_{j_{c}}^{-1}))\mathbf{H}_{c,k_{b}}+k_{k_{b}}\mbox{diag}(\mathbf{H}_{c,k_{b}}^{H}(\mathbf{R}_{\overline{j_{c}}}^{-1}-\mathbf{R}_{j_{c}}^{-1})\mathbf{H}_{c,k_{b}})]$
$\hat{\mathbf{G}}_{{j_{b}},b}^{UL}$ |  $\sum_{\begin{subarray}{c}k_{b}\in\mathcal{U}_{b}\end{subarray}}w_{k_{b}}[\mathbf{H}_{b,b}^{H}\mathbf{F}_{b}(\mathbf{R}_{\overline{k_{b}}}^{-1}-\mathbf{R}_{k_{b}}^{-1}+\beta_{b}\mbox{diag}(\mathbf{R}_{\overline{k_{b}}}^{-1}-\mathbf{R}_{k_{b}}^{-1})){\mathbf{F}_{b}}^{H}\mathbf{H}_{b,b}+k_{k_{b}}\mbox{diag}(\mathbf{H}_{b,b}^{H}{\mathbf{F}_{b}}(\mathbf{R}_{\overline{k_{b}}}^{-1}-\mathbf{R}_{k_{b}}^{-1}){\mathbf{F}_{b}}^{H}\mathbf{H}_{b,b})]$
$\hat{\mathbf{G}}_{\overline{j_{b}},b}^{DL}$ |  $\sum_{\begin{subarray}{c}l_{b}\in\mathcal{D}_{b}\\\ l_{b}\neq j_{b}\end{subarray}}w_{l_{b}}[\mathbf{H}_{l_{b}}^{H}(\mathbf{R}_{\overline{l_{b}}}^{-1}-\mathbf{R}_{l_{b}}^{-1}+\beta_{l_{b}}\mbox{diag}(\mathbf{R}_{\overline{l_{b}}}^{-1}-\mathbf{R}_{l_{b}}^{-1}))\mathbf{H}_{l_{b}}^{H}+k_{c}\mbox{diag}(\mathbf{H}_{l_{b}}^{H}(\mathbf{R}_{\overline{l_{b}}}^{-1}-\mathbf{R}_{l_{b}}^{-1})\mathbf{H}_{l_{b}})]$
$\hat{\mathbf{G}}_{j_{b},\overline{b}}^{UL}$ |  $\sum_{\begin{subarray}{c}c\in\mathcal{B}\\\ c\neq b\end{subarray}}\sum_{\begin{subarray}{c}k_{c}\in\mathcal{U}_{c}\end{subarray}}w_{k_{c}}[\mathbf{H}_{c,b}^{H}{\mathbf{F}_{c}}(\mathbf{R}_{\overline{k_{c}}}^{-1}-\mathbf{R}_{k_{c}}^{-1}+\beta_{c}\mbox{diag}(\mathbf{R}_{\overline{k_{c}}}^{-1}-\mathbf{R}_{k_{c}}^{-1})){\mathbf{F}_{c}}^{H}\mathbf{H}_{c,b}+k_{b}\mbox{diag}(\mathbf{H}_{c,b}^{H}{\mathbf{F}_{c}}(\mathbf{R}_{\overline{k_{c}}}^{-1}-\mathbf{R}_{k_{c}}^{-1}){\mathbf{F}_{c}}^{H}\mathbf{H}_{c,b})]$
$\hat{\mathbf{G}}_{j_{b},\overline{b}}^{DL}$ |  $\sum_{\begin{subarray}{c}c\in\mathcal{B}\\\ c\neq b\end{subarray}}\sum_{\begin{subarray}{c}j_{c}\in\mathcal{D}_{c}\end{subarray}}w_{j_{c}}[\mathbf{H}_{j_{c},b}^{H}(\mathbf{R}_{\overline{j_{c}}}^{-1}-\mathbf{R}_{j_{c}}^{-1}+\beta_{j_{c}}\mbox{diag}(\mathbf{R}_{\overline{j_{c}}}^{-1}-\mathbf{R}_{j_{c}}^{-1}))\mathbf{H}_{j_{c},b}+k_{b}\mbox{diag}(\mathbf{H}_{j_{c},b}^{H}(\mathbf{R}_{\overline{j_{c}}}^{-1}-\mathbf{R}_{j_{c}}^{-1})\mathbf{H}_{j_{c},b}]$

Let $\lambda_{k_{b}}$ and $\psi_{b}$ denote the Lagrange multipliers
associated with the sum-power constraint for UL user $k_{b}\in\mathcal{U}_{b}$
and FD BS $b\in\mathcal{B}$, respectively. Hereafter, for notational
convenience, we define the following matrices

$\mathbf{\Sigma}_{k_{b}}^{1}=\mathbf{H}_{k_{b}}^{H}\mathbf{F}_{b}\mathbf{R}_{\overline{k_{b}}}^{-1}{\mathbf{F}_{b}}^{H}\mathbf{H}_{k_{b}},\quad\mathbf{\Sigma}_{j_{b}}^{1}=\mathbf{H}_{j_{b}}^{H}\mathbf{R}_{\overline{j_{b}}}^{-1}\mathbf{H}_{j_{b}},$
(14a)
$\mathbf{\Sigma}_{k_{b}}^{2}=\hat{\mathbf{G}}_{\overline{k_{b}},b}^{UL}+\hat{\mathbf{G}}_{{k_{b}},b}^{DL}+\hat{\mathbf{G}}_{k_{b},\overline{b}}^{UL}+\hat{\mathbf{G}}_{k_{b},\overline{b}}^{DL}+\lambda_{k_{b}}\mathbf{I},\;\;\mathbf{\Sigma}_{j_{b}}^{2}=\hat{\mathbf{G}}_{j_{b},b}^{UL}+\hat{\mathbf{G}}_{\overline{j_{b}},b}^{DL}+\hat{\mathbf{G}}_{j_{b},\overline{b}}^{UL}+\hat{\mathbf{G}}_{j_{b},\overline{b}}^{DL}+\psi_{b}\mathbf{I}.$
(14b)

By considering the minorized WSR constructed with the gradients (12)-(13),
ignoring the constant terms and the unit-modulus and quantization constraints
(10d), and augmenting the minorized WSR only with the sum-power constraints
leads to the following Lagrangian

$\displaystyle\mathcal{L}=$
$\displaystyle\sum_{b\in\mathcal{B}}\big{[}\sum_{k_{b}\in\mathcal{U}_{b}}(w_{k_{b}}\mbox{lndet}(\mathbf{I}+\mathbf{U}_{k_{b}}^{H}\mathbf{\Sigma}_{k_{b}}^{1}\mathbf{U}_{k_{b}})-\mbox{Tr}(\mathbf{U}_{k_{b}}^{H}\mathbf{\Sigma}_{k_{b}}^{2}\mathbf{U}_{k_{b}})+\lambda_{k_{b}}p_{k_{b}})$
(15)
$\displaystyle+\sum_{j_{b}\in\mathcal{D}_{b}}(w_{j_{b}}\mbox{lndet}(\mathbf{I}+\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{1}{\mathbf{W}_{b}}\mathbf{V}_{j_{b}})-\mbox{Tr}(\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{2}{\mathbf{W}_{b}}\mathbf{V}_{j_{b}})+\psi_{b}p_{b})\big{]}.$

We note that the constraints on the analog part, ommittied in (31), will be
incorporated later.

##  IV Centralized Hybrid Beamforming

This section presents a novel C-HYBF design based on alternating optimization
to solve $\eqref{Largrangian}$ to a local optimum. Hereafter, different sub-
sections are dedicated for the optimization of different variables and at each
step complete information of the other variables is summarized in the
gradients, which are updated at each iteration.

###  IV-A Digital Beamforming

To optimize the digital beamformers $\mathbf{U}_{k_{b}}$ and
$\mathbf{V}_{j_{b}}$ we take the derivatives of (31) with respect to their
conjugates, which yield the following Karush–Kuhn–Tucker (KKT) conditions

$\displaystyle\mathbf{\Sigma}_{k_{b}}^{1}\mathbf{U}_{k_{b}}(\mathbf{I}+\mathbf{U}_{k_{b}}^{H}\mathbf{\Sigma}_{k_{b}}^{1}\mathbf{U}_{k_{b}})^{-1}-\mathbf{\Sigma}_{k_{b}}^{2}\mathbf{U}_{k_{b}}=0,$
(16a) $\displaystyle{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{1}$
$\displaystyle{\mathbf{W}_{b}}\mathbf{V}_{j_{b}}(\mathbf{I}+\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{1}{\mathbf{W}_{b}}\mathbf{V}_{j_{b}})^{-1}-{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{2}{\mathbf{W}_{b}}\mathbf{V}_{j_{b}}=0.$
(16b)

######  Theorem 1.

The WSR maximizing digital beamformers $\mathbf{U}_{k_{b}}$ and
$\mathbf{V}_{j_{b}}$, when the remaining variables are fixed, can be computed
as the GDE solution of the pair of the following matrices

$\mathbf{U}_{k_{b}}=\mathbf{D}_{d_{k_{b}}}(\mathbf{\Sigma}_{k_{b}}^{1},\mathbf{\Sigma}_{k_{b}}^{2}),$
(17a)
$\mathbf{V}_{j_{b}}=\mathbf{D}_{d_{j_{b}}}({\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{1}{\mathbf{W}_{b}},{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{2}{\mathbf{W}_{b}}),$
(17b)

where the matrix $\mathbf{D}_{d_{k_{b}}}$ ($\mathbf{D}_{d_{k_{b}}}$) selects
$d_{k_{b}}(d_{k_{b}})$ GDEs.

###### Proof.

The proof is provided in Appendix A. ∎

Note that GDEs provide optimal beamforming directions but not the powers.
Therefore, we consider scaling them to unit-norm columns by preserving only
optimized beamforming directions and design the optimal power allocation
scheme separately. We remark that such operation does not affect the
optimality of the optimized digital beamformers.

###  IV-B Analog Beamforming

In this section, we present a novel design for the analog beamformer
$\mathbf{W}_{b}$ for FD BS $b\in\mathcal{B}$ in a multicell scenario, assuming
the remaining variables to be fixed. We note that herein we focus solely on
the fully connected structure, however, the proposed method can be easily
extended to the partially connected case. By considering only the dependence
of the WSR on the unconstrained analog beamformer $\mathbf{W}_{b}$, we have to
solve the following problem

$\displaystyle\underset{\begin{subarray}{c}\mathbf{W}_{b}\end{subarray}}{\max}\sum_{j_{b}\in\mathcal{D}_{b}}$
$\displaystyle[w_{j_{b}}\mbox{lndet}(\mathbf{I}+\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{1}\mathbf{W}_{b}\mathbf{V}_{j_{b}})-\mbox{Tr}(\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{2}\mathbf{W}_{b}\mathbf{V}_{j_{b}})].$
(18)

We take its derivative with respect to the conjugate of $\mathbf{W}_{b}$ which
leads to the following KKT condition

$\displaystyle\sum_{j_{b}\in\mathcal{D}}$
$\displaystyle(\mathbf{\Sigma}_{j_{b}}^{1}\mathbf{W}_{b}\mathbf{V}_{j_{b}}\mathbf{V}_{j_{B}}^{H}(\mathbf{I}+\mathbf{V}_{j_{b}}\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{1}\mathbf{W}_{b})^{-1}-\mathbf{\Sigma}_{j_{b}}^{2}\mathbf{W}_{b}\mathbf{V}_{j_{b}}\mathbf{V}_{j_{b}}^{H})=0.$
(19)

######  Theorem 2.

The vectorized unconstrained analog beamformer $\mathbf{W}_{b}$ which is
common to all the DL users in set $\mathcal{D}_{b}$, can be optimized as one
GDE solution of the pair of the sum of following matrices

$\displaystyle\mbox{vec}(\mathbf{W}_{b})=\mathbf{D}_{1}(\sum_{j_{b}\in\mathcal{D}_{b}}(\mathbf{V}_{j_{b}}$
$\displaystyle\mathbf{V}_{j_{b}}^{H}(\mathbf{I}+\mathbf{V}_{j_{b}}\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{1}\mathbf{W}_{b})^{-1})^{T}\otimes\mathbf{\Sigma}_{j_{b}}^{1},\;\sum_{j_{b}\in\mathcal{D}_{b}}(\mathbf{V}_{j_{b}}\mathbf{V}_{j_{b}}^{H})^{T}\otimes(\mathbf{\Sigma}_{j_{b}}^{2})).$
(20)

###### Proof.

The proof is provided in Appendix B ∎

The result stated in Theorem 2 provides the optimized vectorized unconstrained
analog beamformer. Operation $\mbox{unvec}(\mbox{vec}(\mathbf{W}_{b}))$ is
required to reshape it into correct dimensions. To projects its elements on
the set $\mathcal{P}_{b}$, we normalize its amplitude to unit-norm with
$\angle\cdot$ and pass it through the quantizer such that
$\mathbf{W}_{b}=\mathbb{Q}_{b}(\angle\mathbf{W}_{b}(m,n))\in\mathcal{P}_{b},\forall
m,n$. We remark that such operation result into optimality loss, which depends
on the resolution of the phase-shifts, and will be evaluated in in Section VI.

###  IV-C Analog Combining

Optimization of the analog combiner $\mathbf{F}_{b}$ is straightforward
compared to the analog beamformer. Note that the analog combiners do not
appear in the trace operators of (31) but only in
$\mathbf{\Sigma}_{k_{b}}^{1},\forall b,$ as they do not generate any
interference towards other links. Therefore, to optimize $\mathbf{F}_{b}$ we
can directly consider the original problem (10a) which is purely concave with
respect to $\mathbf{F}_{b}$.

The objective of the analog combiner $\mathbf{F}_{b}$ is to combine the
received covariance matrices at the antenna level such that the WSR is
maximized. Let $(\mathbf{R}_{k_{b}}^{a})$ $\mathbf{R}_{\overline{k_{b}}}^{a}$
denote the (signal plus) interference and noise covariance matrices received
at the antennas of the FD BS $b\in\mathcal{B}$ to be combined with
$\mathbf{F}_{b}$. Given $\mathbf{R}_{k_{b}}^{a}$ and
$\mathbf{R}_{\overline{k_{b}}}^{a}$, the matrices $\mathbf{R}_{k_{b}}$ and
$\mathbf{R}_{\overline{k_{b}}}$ can be recovered as
$\mathbf{R}_{k_{b}}={\mathbf{F}_{b}}^{H}\mathbf{R}_{k_{b}}^{a}\mathbf{F}_{b}$
and
$\mathbf{R}_{\overline{k_{b}}}={\mathbf{F}_{b}}^{H}\mathbf{R}_{\overline{k_{b}}}^{a}\mathbf{F}_{b}$.
Problem (10a), with respect to the unconstrained analog combiner
$\mathbf{F}_{b}$, by using the properties of the logarithm function can be
restated as

$\displaystyle\underset{\begin{subarray}{c}\mathbf{F}_{b}\end{subarray}}{\max}\sum_{k_{b}\in\mathcal{U}_{b}}w_{k_{b}}[$
$\displaystyle\mbox{lndet}({\mathbf{F}_{b}}^{H}\mathbf{R}_{k_{b}}^{a}\mathbf{F}_{b})-\mbox{lndet}({\mathbf{F}_{b}}^{H}\mathbf{R}_{\overline{k_{b}}}^{a}\mathbf{F}_{b})].$
(21)

Solving (21), which is purely concave, leads to the following optimal analog
combiner

$\mathbf{F}_{b}=\mathbf{D}_{N_{b}^{RF}}(\sum_{k_{b}\in\mathcal{U}_{b}}w_{k_{b}}\mathbf{R}_{k_{b}}^{a},\sum_{k_{b}\in\mathcal{U}_{b}}w_{k_{b}}\mathbf{R}_{\overline{k_{b}}}^{a}),$
(22)

where the matrix $\mathbf{D}_{N_{b}^{RF}}$ selects GDEs equal to the number of
receive RF chains $N_{b}^{RF}$ at the FD BS $b\in\mathcal{B}$. To meet the
constraints for (22), we normalize its amplitudes with the operator
$\angle\cdot$ and pass it through the quantizer such that
$\mathbf{F}_{b}=\mathbb{Q}_{b}(\angle\mathbf{F}_{b}(m,n))\in\mathcal{P}_{b}$.
Recall that this result into optimality loss, which will be evaluated in
Section VI.

###  IV-D Optimal Power Allocation

Let $\mathbf{P}_{k_{b}}$ and $\mathbf{P}_{j_{b}}$ denote the stream power
matrices for the UL user $k_{b}\in\mathcal{U}_{b}$ and DL user
$j_{b}\in\mathcal{D}_{b}$, respectively, to be included in the beamformers
$\mathbf{U}_{k_{b}}$ and $\mathbf{V}_{k_{b}}$, respectively. Given the
normalized digital beamformers with the unit-norm columns, the power
allocation problems for $\mathbf{P}_{k_{b}}$ and $\mathbf{P}_{j_{b}}$ can be
formally stated as

$\displaystyle\underset{\mathbf{P}_{k_{b}}}{\max}$ $\displaystyle\quad
w_{k_{b}}\mbox{lndet}(\mathbf{I}+\mathbf{U}_{k_{b}}^{H}\mathbf{\Sigma}_{k_{b}}^{1}\mathbf{U}_{k_{b}}\mathbf{P}_{k_{b}})-\mbox{Tr}(\mathbf{U}_{k_{b}}^{H}\mathbf{\Sigma}_{k_{b}}^{2}\mathbf{U}_{k_{b}}\mathbf{P}_{k_{b}}),$
(23a) $\displaystyle\underset{\mathbf{P}_{j_{b}}}{\max}$ $\displaystyle
w_{j_{b}}\mbox{lndet}(\mathbf{I}+\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{1}{\mathbf{W}_{b}}\mathbf{V}_{j_{b}}\mathbf{P}_{j_{b}})-\mbox{Tr}(\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{2}{\mathbf{W}_{b}}\mathbf{V}_{j_{b}}\mathbf{P}_{j_{b}}).$
(23b)

Recall that right multiplying the unit-norm beamformers with powers is
equivalent to not normalized digital beamformers which also includes the
optimal power allocation. Solving the problems stated above leads to the
following optimal power allocation scheme

$\mathbf{P}_{k_{b}}=(w_{k_{b}}(\mathbf{U}_{k_{b}}^{H}\mathbf{\Sigma}_{k_{b}}^{2}\mathbf{U}_{k_{b}})^{-1}-(\mathbf{U}_{k_{b}}^{H}\mathbf{\Sigma}_{k_{b}}^{1}\mathbf{U}_{k_{b}})^{-1})^{+},$
(24a)
$\displaystyle\mathbf{P}_{j_{b}}=(w_{j_{b}}(\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{2}{\mathbf{W}_{b}}\mathbf{V}_{j_{b}})^{-1}-(\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{1}{\mathbf{W}_{b}}\mathbf{V}_{j_{b}})^{-1})^{+},$
(24b)

where $(\mathbf{X})^{+}=max\\{\mathbf{0},\mathbf{X}\\}$. Note that given the
other variables fixed, the power allocation is optimal as each problem is made
of a concave part and a linear part which is simultaneously convex and
concave. Given the optimal stream powers, we can search for the Lagrange
multipliers satisfying the total sum-power constraint. Let $\mathbf{P}^{DL}$
and $\mathbf{P}^{UL}$ denote the collection of powers in DL and UL,
respectively, and let $\mathbf{\Lambda}$ and $\mathbf{\Psi}$ denote the
collection of multipliers for $\lambda_{k_{b}}$ and $\psi_{b}$, respectively.
Given $\eqref{optimal_pow}$, consider the dependence of the Lagrangian only on
the multipliers and powers as
$\mathcal{L}(\mathbf{\Lambda},\mathbf{\Psi},\mathbf{P}^{DL},\mathbf{P}^{UL})$,
obtained by including the power matrices $\mathbf{P}_{k_{b}}$ and
$\mathbf{P}_{j_{b}}$ in (31).

The multipliers in $\mathbf{\Lambda}$ and $\mathbf{\Psi}$ should be such that
the Lagrangian is finite and the values of multipliers are strictly positive,
i.e.,

$\displaystyle\underset{\mathbf{\Psi},\mathbf{\Lambda}}{\min}$
$\displaystyle\underset{\mathbf{P}^{DL},\mathbf{P}^{UL}}{\max}\quad\mathcal{L}(\mathbf{\Lambda},\mathbf{\Psi},\mathbf{P}^{DL},\mathbf{P}^{UL}),$
(25)
$\displaystyle\quad\quad\mbox{s.t.}\quad\quad\mathbf{\Psi},\mathbf{\Lambda}\succeq
0.$

The dual function

$\max_{\mathbf{P}^{DL},\mathbf{P}^{UL}}\quad\mathcal{L}(\mathbf{\Lambda},\mathbf{\Psi},\mathbf{P}^{DL},\mathbf{P}^{UL})$
(26)

is the pointwise supremum of a family of functions of
$\mathbf{\Psi},\mathbf{\Lambda}$, it is convex [28] and the globally optimal
values for $\mathbf{\Psi}$ and $\mathbf{\Lambda}$ can be found by using any of
the numerous convex-optimization techniques. In this work, we adopt the
Bisection method. Let $\underline{\lambda_{k_{b}}},\underline{\psi_{b}}$ and
$\overline{\psi_{b}},\overline{\lambda_{k_{b}}}$ denote the upper and lower
bounds for searching the multipliers $\psi_{b}$ and $\lambda_{k_{b}}$,
respectively, and let $[0,\lambda_{k_{b}}^{max}]$ and $[0,\psi_{b}^{max}]$
denote their search range. Note also that as the GDE solution is computed
given fixed multipliers, doing water-filling for the powers while searching
for the multipliers leads to non diagonal power matrices. Hence, consider a
SVD of the powers as
$[\mathbf{L}_{\mathbf{P}_{i}}^{svd},\mathbf{D}_{\mathbf{P}_{i}}^{svd},\mathbf{R}_{\mathbf{P}_{i}}^{svd}]=svd(\mathbf{P}_{i})$,
with $\mathbf{P}_{i}=\mathbf{P}_{k_{b}}$ or
$\mathbf{P}_{i}=\mathbf{P}_{j_{b}}$, and the matrices
$\mathbf{L}_{\mathbf{P}_{i}}^{svd},\mathbf{D}_{\mathbf{P}_{i}}^{svd}$ and
$\mathbf{R}_{\mathbf{P}_{i}}^{svd}$ denote the left unitary, diagonal and
right unitary matrices obtained from the SVD of $\mathbf{P}_{i}$. The diagonal
structure of the power matrices while searching for the multipliers can be re-
established as $\mathbf{P}_{i}=\mathbf{D}_{i}$, with $i\in\mathcal{U}_{b}$ or
$i\in\mathcal{D}_{b}$.

By using the closed form expressions derived above, the complete alternating
optimization based C-HYBF procedure to solve (10a) is formally stated in
Algorithm 1.

Algorithm 1 Centralized Hybrid Beamforming

Given: The CSI and rate weights.
Initialize:
$\mathbf{W}_{b},\mathbf{F}_{b},\mathbf{V}_{j_{b}},\mathbf{U}_{k_{b}},\quad\forall
j_{b}\;\&\;\forall k_{b}$.
Set:
$\underline{\lambda_{k_{b}}}=0,\overline{\lambda_{k_{b}}}=\lambda_{k_{b}}^{max},\underline{\psi_{b}}=0,\overline{\psi_{b}}=\psi_{b}^{max}$,
$\forall k_{b}\;\&\;\forall b$
Repeat until convergence

for $b=1:B$

Compute $\mathbf{W}_{b}$ with (20), do $\mbox{unvec}(\mathbf{W}_{b})$ and get
$\angle\mathbf{W}_{b}$

for: $j_{b}=1:D_{b}$

Compute
$\hat{\mathbf{G}}_{j_{b},b}^{UL},\hat{\mathbf{G}}_{{j_{b}},b}^{DL},\hat{\mathbf{G}}_{j_{b},\overline{b}}^{UL},\hat{\mathbf{G}}_{j_{b},\overline{b}}^{DL}$
from Table II Optimize $\mathbf{V}_{j_{b}}$ with (17a) and normalize its
columns

Next $j_{b}$

Repeat until convergence

set $\psi_{b}=(\underline{\psi_{b}}+\overline{\psi_{b}})/2$

for $j_{b}=1:D_{b}$

Compute $\mathbf{P}_{j_{b}}$ with (24), do SVD, set
$\mathbf{P}_{j_{b}}=\mathbf{D}_{P_{j_{b}}}^{svd}$ Set
$\mathbf{Q}_{j_{b}}=\mathbf{W}_{b}\mathbf{V}_{j_{b}}\mathbf{P}_{j_{b}}\mathbf{V}_{j_{b}}^{H}\mathbf{W}_{b}^{H}$

Next $j_{b}$

if constraint for $\psi_{b}$ is violated

set $\underline{\psi_{b}}=\psi_{b}$

else $\overline{\psi_{b}}=\psi_{b}$

for: $k_{b}=1:K_{b}$

Compute
$\hat{\mathbf{G}}_{k_{b},b}^{UL},\hat{\mathbf{G}}_{{k_{b}},b}^{DL},\hat{\mathbf{G}}_{k_{b},\overline{b}}^{UL},\hat{\mathbf{G}}_{k_{b},\overline{b}}^{DL}$
from Table II.

Optimize $\mathbf{U}_{k_{b}}$ with (17b) and normalize its columns

Repeat until convergence

set
$\lambda_{k_{b}}=(\underline{\lambda_{k_{b}}}+\overline{\lambda_{k_{b}}})/2$

Compute $\mathbf{P}_{k_{b}}$ with (24), do SVD, set
$\mathbf{P}_{k_{b}}=\mathbf{D}_{P_{k_{b}}}^{svd}$ Set
$\mathbf{T}_{{}_{b}}=\mathbf{U}_{k_{b}}\mathbf{P}_{k_{b}}\mathbf{U}_{k_{b}}^{H}$

if constraint for $\lambda_{k_{b}}$ is violated

set $\underline{\lambda_{k_{b}}}=\lambda_{k_{b}}$

else $\overline{\lambda_{k_{b}}}=\lambda_{k_{b}}$

Next $k_{b}$

Next $b$

Quantize $\mathbf{W}_{b}$ and $\mathbf{F}_{b},$ with
$\mathbb{Q}_{b}(\cdot),\forall b$

###  IV-E Convergence of C-HYBF

The convergence of Algorithm 1 can be proved by using the minorization theory
[27], alternating or cyclic optimization [27], Lagrange dual function [28],
saddle-point interpretation [28] and KKT conditions [28]. For the WSR cost
function (10a), we construct its minorizer, which is a touching lower bound
for (10a), hence we can write

WSR
$\displaystyle\geq\underline{\mbox{WSR}}=\underline{\mbox{WR}}_{k_{b},b}^{UL}+\underline{\mbox{WSR}}_{\overline{k_{b}},b}^{UL}+\underline{\mbox{WR}}_{j_{b},b}^{DL}+\underline{\mbox{WSR}}_{\overline{j_{b}},b}^{DL}+\underline{\mbox{WSR}}_{\overline{b}}^{DL}+\underline{\mbox{WSR}}_{\overline{b}}^{UL}.$
(27)

The minorized WSR, which is concave in $\mathbf{T}_{k_{b}}$ and
$\mathbf{Q}_{j_{b}}$, has the same gradient of the original WSR maximization
problem (10a), hence the KKT conditions are not affected. Reparameterizing
$\mathbf{T}_{k_{b}}$ or $\mathbf{Q}_{j_{b}}$ in terms of
$\mathbf{W}_{b},\mathbf{V}_{j_{b}},\forall j_{b}\in\mathcal{D}_{b}$, or
$\mathbf{U}_{k_{b}},\forall k_{b}\in\mathcal{U}_{b}$, respectively,
augementing the minorized WSR cost function with the Lagrange multipliers and
power constraints leads to (31). By incorporating further the power matrices
we get to
$\mathcal{L}(\mathbf{\Lambda},\mathbf{\Psi},\mathbf{P}^{DL},\mathbf{P}^{UL})$.
Every alternating update of the $\mathcal{L}$ for the variables
$\mathbf{W}_{b},\mathbf{F}_{b},\forall
b\in\mathcal{B},\mathbf{V}_{j_{b}},\forall{j_{b}}\in\mathcal{D}_{b},\mathbf{U}_{k_{b}},\forall
k_{b}\in\mathcal{U}_{b},\mathbf{P}_{k_{b}},\mathbf{P}_{j_{b}},\lambda_{k_{b}}$
and $\psi_{b}$, leads to a monotonic increase of the WSR, which assures
convergence. For the KKT conditions, at the convergence point, the gradients
of $\mathcal{L}$ for $\mathbf{V}_{j_{b}},\mathbf{W}_{b},\mathbf{U}_{k_{b}}$ or
$\mathbf{P}_{k_{b}},\mathbf{P}_{j_{b}}$ correspond to the gradients of the
Lagrangian of the original problem (10a), and hence the sub-optimal solution
for the minorized WSR matches the sub-optimal solution of the original
problem. For the fixed analog and digital beamformers, $\mathcal{L}$ is
concave in powers, hence we have strong duality for the saddle point, i.e.,

$\max_{\mathbf{P}^{DL},\mathbf{P}^{UL}}\min_{\mathbf{\Lambda},\mathbf{\Psi}}\quad\mathcal{L}(\mathbf{\Lambda},\mathbf{\Psi},\mathbf{P}^{UL},\mathbf{P}^{DL}).$
(28)

Let $\mathbf{X}^{*}$ and $x^{*}$ denote the optimal solution for matrix
$\mathbf{X}$ or scalar $x$ at the convergence, respectively. As each iteration
leads to a monotonic increase in the WSR and the power are updated by
satisfying the sum-power constraint, at the convergence point, the solution of
the optimization problem

$\min_{\mathbf{\Lambda},\mathbf{\Psi}}\quad\mathcal{L}(\mathbf{V}_{j_{b}}^{*},{\mathbf{W}_{b}}^{*},{\mathbf{F}_{b}}^{*},\mathbf{U}_{b}^{*},{\mathbf{P}^{DL}}^{*},{\mathbf{P}^{UL}}^{*},\mathbf{\Lambda},\mathbf{\Psi})$
(29)

satisfies the KKT conditions for the powers in $\mathbf{P}^{DL}$ and
$\mathbf{P}^{UL}$ and the complementary slackness conditions

$\psi_{b}^{*}\,(p_{b}-\sum_{j_{b}\in\mathcal{D}_{b}}\mbox{Tr}({\mathbf{W}_{b}}^{*}\mathbf{V}_{j_{b}}^{*}\mathbf{P}_{j_{b}}^{*}\mathbf{V}_{j_{b}}^{*\,H}{\mathbf{W}_{b}}^{*\,H}))=0,\quad\lambda_{k_{b}}^{*}\,(p_{k_{b}}-\mbox{Tr}(\mathbf{U}_{k_{b}}^{*}\mathbf{P}_{k_{b}}^{*}\mathbf{U}_{k_{b}}^{*\,H}))=0,$
(30a)

with the individual factors in the products being non-negative.

##  V Parallel and Distributed Implementation

As well-discussed in Section I, implementation of Algorithm $1$ in a real-time
FD network is infeasible in practice as it requires centralized
implementation. To overcome its drawbacks, we present the concept of
P$\&$D-HYBF for mmWave based on cooperation among FD BSs. Such design enables
per-link independent optimization of the beamformers for each FD BS, thus
avoiding full CSI exchange with the CN every CCT and communicating back the
optimized variables to the FD BSs at the millisecond scale, leading to
significant reduction in the communication overhead.

To develop the P$\&$D-HYBF scheme, we assume the following:

  1. 1.

there exist a feedback link among the FD BSs and they cooperate by exchanging
information about the digital beamformers, analog beamformers and analog
combiners via the feedback link;

  2. 2.

local CSI is accessible by each FD BS;

  3. 3.

each FD BS has multiple computational processors dedicated for UL and DL;

  4. 4.

the computations take place at the FD BSs in each cell in a synchronized
manner.

Recall that the MM optimization technique allowed to write the Lagrangian of
the original WSR problem (10a) as

$\displaystyle\mathcal{L}=$
$\displaystyle\sum_{b\in\mathcal{B}}\big{[}\sum_{k_{b}\in\mathcal{U}_{b}}(w_{k_{b}}\mbox{lndet}(\mathbf{I}+\mathbf{U}_{k_{b}}^{H}\mathbf{\Sigma}_{k_{b}}^{1}\mathbf{U}_{k_{b}})-\mbox{Tr}(\mathbf{U}_{k_{b}}^{H}\mathbf{\Sigma}_{k_{b}}^{2}\mathbf{U}_{k_{b}})+\lambda_{k_{b}}p_{k_{b}})$
(31)
$\displaystyle+\sum_{j_{b}\in\mathcal{D}_{b}}(w_{j_{b}}\mbox{lndet}(\mathbf{I}+\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{1}{\mathbf{W}_{b}}\mathbf{V}_{j_{b}})-\mbox{Tr}(\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{2}{\mathbf{W}_{b}}\mathbf{V}_{j_{b}})+\psi_{b}p_{b})\big{]},$

from which it is evident that when optimizing one variable, assuming the
remaining variables to be fixed, only the gradients appearing in
$\mathbf{\Sigma}_{k_{b}}^{2}$ or $\mathbf{\Sigma}_{j_{b}}^{2}$ associated with
that particular link are required. Therefore, the gradients summarize complete
information about all the remaining interfering links in the network. From a
practical point-of-view, the gradients for each link take into account the
interference generated towards all the other links, and hence limit greedy
behaviour while optimizing its beamforming directions. However,
$\eqref{Largrangian}$ is coupled among different links as the covariance
matrices of other users directly appear in the gradients, which vary at the
update of each beamformer.

To decouple (31) into local per-link independent optimization sub-problems, we
assume that each FD BS has some memory to save information. Hereafter,
overline will emphasize that the variables are only local and saved in the
memory of each FD BS. We introduce the following local variables

$\overline{\mathbf{L}_{k_{b}}^{In}}=\hat{\mathbf{G}}_{\overline{k_{b}},b}^{UL}+\hat{\mathbf{G}}_{k_{b},b}^{DL},\;\overline{\mathbf{L}_{k_{b}}^{Out}}=\hat{\mathbf{G}}_{k_{b},\overline{b}}^{UL}+\hat{\mathbf{G}}_{k_{b},\overline{b}}^{DL},\quad\forall
k_{b}\in\mathcal{U}_{b},$ (32a)
$\overline{\mathbf{L}_{j_{b}}^{In}}=\hat{\mathbf{G}}_{j_{b},b}^{UL}+\hat{\mathbf{G}}_{\overline{j_{b}},b}^{DL},\;\;\overline{\mathbf{L}_{j_{b}}^{Out}}=\hat{\mathbf{G}}_{j_{b},\overline{b}}^{UL}+\hat{\mathbf{G}}_{j_{b},\overline{b}}^{DL},\quad\forall
j_{b}\in\mathcal{D}_{b},$ (32b)

to be saved in the memory of each FD BS $b\in\mathcal{B}$. The local variables
$\overline{\mathbf{L}_{k_{b}}^{In}}$ and $\overline{\mathbf{L}_{k_{b}}^{out}}$
save information about the overall interference generated inside and outside
the cell by the beamformer of UL user $k_{b}\in\mathcal{U}_{b}$, respectively.
Similarly, the local variables $\overline{\mathbf{L}_{j_{b}}^{In}}$ and
$\overline{\mathbf{L}_{j_{b}}^{Out}}$ save information about the interference
generated in the same cell and in the neighbouring cells by the FD BS
$b\in\mathcal{B}$, respectively, while serving its DL user
$j_{b}\in\mathcal{D}_{b}$. Note that each FD BS can update the in-cell local
variables $\overline{\mathbf{L}_{k_{b}}^{In}}$ and
$\overline{\mathbf{L}_{j_{b}}^{In}}$ by itself. A feedback from the
neighboring BSs is required only to update the out-cell variables
$\overline{\mathbf{L}_{k_{b}}^{Out}}$ and
$\overline{\mathbf{L}_{j_{b}}^{Out}}$. To save information about the
interference-plus-noise covariance matrices at the RF chains and antenna level
(for the analog combiner), we define the following local variables

$\overline{\mathbf{R}_{\overline{j_{b}}}^{-1}}=\mathbf{R}_{\overline{j_{b}}}^{-1},\quad\forall
j_{b},$ (33a)
$\overline{\mathbf{R}_{\overline{k_{b}}}^{-1}}=\mathbf{R}_{\overline{k_{b}}}^{-1},\quad\overline{\mathbf{R}_{k_{b}}^{a}}=\mathbf{R}_{k_{b}}^{a},\quad\overline{\mathbf{R}_{\overline{k_{b}}}^{a}}=\mathbf{R}_{\overline{k_{b}}}^{a},\quad\forall
k_{b}.$ (33b)

For notational compactness, similar to (14), we also define the following
variables

$\mathbf{Z}_{k_{b}}^{1}=\mathbf{H}_{k_{b}}^{H}\mathbf{F}_{b}\overline{\mathbf{R}_{\overline{k_{b}}}^{-1}}{\mathbf{F}_{b}}^{H}\mathbf{H}_{k_{b}},\quad\mathbf{Z}_{k_{b}}^{2}=\overline{\mathbf{L}_{k_{b}}^{In}}+\overline{\mathbf{L}_{k_{b}}^{out}}+\lambda_{k_{b}}\mathbf{I},\quad\forall
k_{b},$ (34a)
$\mathbf{Z}_{j_{b}}^{1}=\mathbf{H}_{j_{b}}^{H}\overline{\mathbf{R}_{\overline{j_{b}}}^{-1}}\mathbf{H}_{j_{b}},\quad\;\;\mathbf{Z}_{j_{b}}^{2}=\overline{\mathbf{L}_{j_{b}}^{In}}+\overline{\mathbf{L}_{j_{b}}^{out}}+\psi_{b}\mathbf{I},\quad\forall
j_{b},$ (34b)

which are function only of the fixed local variables. By replacing the
gradients with the fixed local variables, the Lagrangian (31) can be rewritten
as a function of the saved local variables as

$\displaystyle\mathcal{L}=\sum_{b\in\mathcal{B}}$
$\displaystyle[\sum_{k_{b}\in\mathcal{U}_{b}}(w_{k_{b}}\mbox{lndet}(\mathbf{I}+\mathbf{U}_{k_{b}}^{H}\mathbf{Z}_{k_{b}}^{1}\mathbf{U}_{k_{b}})-\mbox{Tr}(\mathbf{U}_{k_{b}}^{H}\mathbf{Z}_{k_{b}}^{2}\mathbf{U}_{k_{b}})+\lambda_{k_{b}}p_{k_{b}})$
(35)
$\displaystyle+\sum_{j_{b}\in\mathcal{D}_{b}}(w_{j_{b}}\mbox{lndet}(\mathbf{I}+\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{Z}_{j_{b}}^{1}\mathbf{W}_{b}\mathbf{V}_{j_{b}})-\mbox{Tr}(\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{Z}_{j_{b}}^{1}\mathbf{W}_{b}\mathbf{V}_{j_{b}})+\psi_{b}p_{b})\big{]}.$

In contrast to (31), (35) becomes fully decoupled per-cell as it is a function
only of the local variables, which are fixed for each FD BS and will be
updated solely when a feedback is received from the neighbouring BSs. However,
note that in UL and DL, the optimization of the analog combiner
$\mathbf{F}_{b}$ in $\mathbf{Z}_{k_{b}}^{1}$ and analog beamformer
$\mathbf{W}_{b}$, $\forall b\in\mathcal{B}$, is still coupled as they are
common among the UL and DL users in the same cell, respectively. Moreover, the
analog beamformer $\mathbf{W}_{b}$ also affects the total transmit power of
each FD BS, posing a serious challenge for enabling per-link independent
optimization in DL from (35) and making per-link independent optimization
infeasible. Handling of the coupling constraints and P$\&$D optimization for
HYBF from (35) is discussed in the following.

###  V-A Per-Link Independent Sub-Problems in UL

Each UL user has its own sum-power constraint but the analog combiner
$\mathbf{F}_{b}$, appearing in $\textbf{Z}_{k_{b}}^{1}$, is common among all
the UL users in the same cell. To decouple the optimization, we assume that FD
BS $b\in\mathcal{B}$ updates $\mathbf{F}_{b}$ only after updating all the
digital beamformers $\mathbf{U}_{k_{b}},\forall k_{b}\in\mathcal{U}_{b}$.
Given this assumption and fixed local variables, UL WSR maximization problem
for each FD BS reduces into three layers of sub-problems. At the bottom layer,
FD BS $b\in\mathcal{B}$ has to solve independent sub-problems to update
$\mathbf{U}_{k_{b}}$, whose optimization is fully decoupled and therefore can
be executed in parallel $\forall k_{b}$. At the middle layer, FD BS
$b\in\mathcal{B}$ has to independently update the stream power matrix
$\mathbf{P}_{k_{b}}$ while searching the multiplier $\lambda_{k_{b}}$,
$\forall k_{b}$, for its independent power constraint. Finally, at the top
layer, once the two-layer UL sub-problems are solved, only one update of the
common analog combiner is required. Fig. 4 highlights the idea of the proposed
per-link decomposition for the UL WSR for FD BS $b\in\mathcal{B}$ into three
sub-layers, and the sub-problems at each layer must solved from the bottom to
the top.

Due to per-link independent decomposition of MM obtained with fixed local
variables, the Lagrangian for the UL user $k_{b}\in\mathcal{U}_{b}$ with
independent sum-power constraint $p_{k_{b}}$, can be written as

$\displaystyle\mathcal{L}_{k_{b}}=w_{k_{b}}\mbox{lndet}(\mathbf{I}+\mathbf{U}_{k_{b}}^{H}\mathbf{Z}_{k_{b}}^{1}\mathbf{U}_{k_{b}})-\mbox{Tr}(\mathbf{U}_{k_{b}}^{H}\mathbf{Z}_{k_{b}}^{2}\mathbf{U}_{k_{b}})+\lambda_{k_{b}}p_{k_{b}},$
(36)

in which for the bottom layer the analog combiner $\mathbf{F}_{b}$ in
$\mathbf{Z}_{k_{b}}^{1}$ and the powers are fixed. To optimize
$\mathbf{U}_{k_{b}}$, a derivative of (36) can be taken, which leads to a
similar KKT condition as (16a), with $\mathbf{\Sigma}_{k_{b}}^{i}$ replaced
with $\mathbf{Z}_{k_{b}}^{i},\forall i$. By following a similar proof of
Appendix A, it can be easily shown that the WSR maximizing
$\mathbf{U}_{k_{b}}$ for (36) can be computed as

$\displaystyle\mathbf{U}_{k_{b}}=\mathbf{D}_{d_{k_{b}}}(\mathbf{Z}_{k_{b}}^{1},\mathbf{Z}_{k_{b}}^{2}).$
(37)

Note that (37) can be computed in parallel by the multi-processor FD BS
$b\in\mathcal{B}$, $\forall k_{b}\in\mathcal{U}_{b}$ as the common part
$\mathbf{F}_{b}$ will be updated at the top layer and information about the
interference generated towards other links is captured in the local variables,
equals to the gradients from the previous feedback.

At the middle layer, the power optimization remains decoupled as each UL user
has its own sum-power constraint and the local variables are fixed. To find
the optimal $\mathbf{P}_{k_{b}}$ in parallel $\forall k_{b}$, we first
consider the normalization of the columns of (37) to unit-norm $\forall
k_{b}$, and the independent power allocation problem for P$\&$D-HYBF in UL can
be formally stated as

$\underset{\mathbf{P}_{k_{b}}}{\max}\quad[w_{k_{b}}\mbox{lndet}(\mathbf{I}+\mathbf{U}_{k_{b}}^{H}\mathbf{Z}_{k_{b}}^{1}\mathbf{U}_{k_{b}}\mathbf{P}_{k_{b}})-\mbox{Tr}(\mathbf{U}_{k_{b}}^{H}\mathbf{Z}_{k_{b}}^{2}\mathbf{U}_{k_{b}}\mathbf{P}_{k_{b}})].$
(38)

Solving (38) independently $\forall k_{b}$ yields the following parallel power
allocation scheme

$\displaystyle\mathbf{P}_{k_{b}}=(w_{k_{b}}($
$\displaystyle\mathbf{U}_{k}^{H}\mathbf{Z}_{k_{b}}^{1}\mathbf{U}_{k_{b}})^{-1}-(\mathbf{U}_{k_{b}}^{H}\mathbf{Z}_{k_{b}}^{2}\mathbf{U}_{k_{b}})^{-1})^{+},$
(39)

which can be computed while searching for the multiplier $\lambda_{k_{b}}$
associated with its independent sum-power constraint in parallel $\forall
k_{b}$. Note that the power allocation problem is made of a concave part and a
linear part with the fixed local variables and therefore simultaneously
concave and the solution (39) results into optimal power allocation scheme.
Multiplier $\lambda_{k_{b}}$, $\forall k_{b}$, should be such that (36) is
finite and the value of $\lambda_{k_{b}}$ should be strictly positive. It can
be searched by solving the following problem in parallel

$\displaystyle\underset{\lambda_{k_{b}}}{\min}$
$\displaystyle\underset{\mathbf{P}_{k_{b}}}{\max}\quad\mathcal{L}_{k_{b}}(\lambda_{k_{b}},\mathbf{P}_{k_{b}}),$
(40) $\displaystyle\quad\quad\mbox{s.t.}\quad\quad\lambda_{k_{b}}\succeq 0,$

while independently allocating the powers at the middle layer $\forall k_{b}$.
The dual function

$\underset{\mathbf{P}_{k_{b}}}{\max}\quad\mathcal{L}_{k_{b}}(\lambda_{k_{b}},\mathbf{P}_{k_{b}}),$
(41)

is convex [28] and can be solved with the Bisection method, as for the C-HYBF
scheme. If $\mathbf{P}_{k_{b}}$ becomes non-diagonal, its diagonal structure
can be reestablished as
$\mathbf{P}_{k_{b}}=\mathbf{D}_{\mathbf{P}_{k_{b}}}^{svd}$, where
$\mathbf{D}_{\mathbf{P}_{k_{b}}}^{svd}$ is a diagonal matrix obtained from SVD
of the non-diagonal $\mathbf{P}_{k_{b}}$.

At the top layer, one update of $\mathbf{F}_{b}$ is required $\forall
b\in\mathcal{B}$. Note that simultaneous variation in parallel of the
beamformers $\mathbf{U}_{k_{b}}$ and powers $\mathbf{P}_{k_{b}}$, $\forall
k_{b},$ at the bottom and middle layer vary the received covariance matrices.
This information should be updated in the local variables
$\overline{\mathbf{R}_{k_{b}}^{a}}$ and
$\overline{\mathbf{R}_{\overline{k_{b}}}^{a}}$ at the antenna level, which
$\mathbf{F}_{b}$ should combine (similar to (22)). As each FD BS
$b\in\mathcal{B}$ has complete information about the optimized variables at
the middle and bottom layers, it can use it to update first
$\overline{\mathbf{R}_{k_{b}}^{a}}$ and
$\overline{\mathbf{R}_{\overline{k_{b}}}^{a}},\forall
k_{b}\in\mathcal{U}_{b}$. As the WSR is fully concave with respect to the
analog combiner $\mathbf{F}_{b}$, to optimize it we consider the following
optimization problem

Figure 3: Decomposition of the UL WSR into three layers of sub-problems
$\forall b\in\mathcal{B}$.

Figure 4: Decomposition of the DL WSR into three layers of sub-problems
$\forall b\in\mathcal{B}$.

$\displaystyle\underset{\begin{subarray}{c}\mathbf{F}_{b}\end{subarray}}{\max}\sum_{k_{b}\in\mathcal{U}_{b}}w_{k_{b}}[$
$\displaystyle\mbox{lndet}({\mathbf{F}_{b}}^{H}\overline{\mathbf{R}_{k_{b}}^{a}}\mathbf{F}_{b})-\mbox{lndet}({\mathbf{F}_{b}}^{H}\overline{\mathbf{R}_{\overline{k_{b}}}^{a}}\mathbf{F}_{b})],$
(42)

where the local variables $\overline{\mathbf{R}_{k_{b}}^{a}}$ and
$\overline{\mathbf{R}_{\overline{k_{b}}}^{a}}$ have been recently updated by
using the information from the middle and bottom layers. Problem (42) is fully
concave and solving its leads to the following optimal unconstrained analog
combiner

$\mathbf{F}_{b}=\mathbf{D}_{N_{b}^{RF}}(\sum_{k_{b}\in\mathcal{U}_{b}}w_{k_{b}}\overline{\mathbf{R}_{k_{b}}^{a}},\sum_{k_{b}\in\mathcal{U}_{b}}w_{k_{b}}\overline{\mathbf{R}_{\overline{k_{b}}}^{a}}).$
(43)

To meet the unit-modulus and quantization constraints, we normalize its
amplitudes with $\angle\cdot$ to unit-norm and quantize it as
$\mathbf{F}_{b}=\mathbb{Q}_{b}(\angle\mathbf{F}_{b})\in\mathcal{P}_{b}$. Note
that such operation results into optimality loss, which strictly depends on
the phase resolution of the phase shifters and will be investigated in the
simulations.

###  V-B Per-Link Independent Sub-Problems in DL

Decomposition of the DL WSR is more challenging due to the coupled sum-power
constraint among the DL users in the set $\mathcal{D}_{b},\forall b$.
Moreover, the analog combiner $\mathbf{W}_{b}$ is also common between the DL
users in the same cell and thus affects the total transmit power. To introduce
per-link independent decomposition in DL, we assume that each FD BS
$b\in\mathcal{B}$ first updates the digital beamformers for the DL users with
unit-norm columns, while keeping the Lagrange multiplier $\psi_{b}$ and the
analog beamformer $\mathbf{W}_{b}$ fixed. Furthermore, the powers are included
afterwards, while searching the common multiplier $\psi_{b}$. Given this
assumption, the DL WSR problem, for each FD BS $b\in\mathcal{B}$, decomposes
into three layers of sub-problems. At the bottom layer, each FD BS has to
update the DL beamformers $\mathbf{V}_{j_{b}}$, $\forall j_{b}$, and normalize
its columns to unit-norm, in parallel. At the middle layer, one update of the
analog beamformer $\mathbf{W}_{b}$ is required. Finally, at the top layer, we
have to search for the Lagrange multiplier $\psi_{b}$ satisfying the coupled
sum-power constraint and update the power matrices $\mathbf{P}_{j_{b}}$ for
the DL users in parallel $\forall j_{b}$. Fig. 4 shows the decomposition of
the DL WSR into three layers of sub-problems, which must be solved from the
bottom to the top.

For FD BS $b\in\mathcal{B}$, the Lagrangian for the DL WSR can be written as

$\displaystyle\mathcal{L}_{b}^{DL}=$
$\displaystyle\sum_{j_{b}\in\mathcal{D}_{b}}(w_{j_{b}}\mbox{lndet}(\mathbf{I}+\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{Z}_{j_{b}}^{1}\mathbf{W}_{b}\mathbf{V}_{j_{b}})-\mbox{Tr}(\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{Z}_{j_{b}}^{2}\mathbf{W}_{b}\mathbf{V}_{j_{b}})+\psi_{b}p_{b}).$
(44)

In (44), for the bottom layer, as $\psi_{b}$, $\mathbf{W}_{b}$,
$\mathbf{Z}_{j_{b}}^{i}$ are fixed, optimization of the digital beamformers
remains decoupled. To optimize $\mathbf{V}_{j_{b}}$ a derivative can be taken,
which will lead to a similar KKT condition as (16b), but function of fixed
$\mathbf{Z}_{j_{b}}^{i}$ instead of $\mathbf{\Sigma}_{j_{b}}^{i}$. By
following a similar proof for Appendix A, it can be easily shown that the WSR
maximizing $\mathbf{V}_{j_{b}}$ can be computed as

$\displaystyle\mathbf{V}_{j_{b}}=\mathbf{D}_{d_{j_{b}}}({\mathbf{W}_{b}}^{H}\mathbf{Z}_{j_{b}}^{1}\mathbf{W}_{b},{\mathbf{W}_{b}}^{H}\mathbf{Z}_{j_{b}}^{2}{\mathbf{W}_{b}}).$
(45)

We consider the normalization of the columns of $\mathbf{V}_{j_{b}}$ to unit-
norm in parallel $\forall j_{b}$ such that optimal power allocation could be
included at the top layer. Once the parallel update of the digital beamformers
$\mathbf{V}_{j_{b}},\forall j_{b}$, has been made, at the middle layer, FD BS
$b\in\mathcal{B}$ has to optimize the analog combiner $\mathbf{W}_{b}$. By
considering the unconstrained analog combiner, each FD BS $b\in\mathcal{B}$
has to independently solve the following optimization problem at the middle
layer

$\displaystyle\underset{\begin{subarray}{c}\mathbf{W}_{b}\end{subarray}}{\max}\sum_{j_{b}\in\mathcal{D}_{b}}$
$\displaystyle(w_{j_{b}}\mbox{lndet}(\mathbf{I}+\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{Z}_{j_{b}}^{1}\mathbf{W}_{b}\mathbf{V}_{j_{b}})-\mbox{Tr}(\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{Z}_{j_{b}}^{2}\mathbf{W}_{b}\mathbf{V}_{j_{b}})).$
(46)

Note that each FD BS has complete information about the digital beamformers
optimized at the bottom layer, which must be first used to update
$\overline{\mathbf{R}_{\overline{j_{b}}}^{-1}}$ and
$\overline{\mathbf{L}_{j_{b}}^{In}}$ appearing in $\mathbf{Z}_{j_{b}}^{1}$ and
$\mathbf{Z}_{j_{b}}^{2}$ in (46), respectively, $\forall j_{b}$. To optimize
$\mathbf{W}_{b}$ a derivative of (46) can be taken, which will lead to a
similar KKT condition as (19), function of $\mathbf{Z}_{j_{b}}^{i}$ instead of
$\mathbf{\Sigma}_{j_{b}}^{i},\forall i$. By following a similar proof of
Appendix B, it can be easily shown that $\mathbf{W}_{b}$ can be optimized as

$\displaystyle\mbox{vec}(\mathbf{W}_{b})=\mathbf{D}_{1}($
$\displaystyle\sum_{j_{b}\in\mathcal{D}_{b}}(\mathbf{V}_{j_{b}}\mathbf{V}_{j_{b}}^{H}(\mathbf{I}+\mathbf{V}_{j_{b}}\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{Z}_{j_{b}}^{1}\mathbf{W}_{b})^{-1})^{T}\otimes\mathbf{Z}_{j_{b}}^{1},\;\sum_{j_{b}\in\mathcal{D}_{b}}(\mathbf{V}_{j_{b}}\mathbf{V}_{j_{b}}^{H})^{T}\otimes\mathbf{Z}_{j_{b}}^{2}).$
(47)

The analog combiner $\mathbf{W}_{b}$ optimized according to (47) is
unconstrained and vectorized. Therefore, we do
unvec$(\mbox{vec}(\mathbf{W}_{b}))$ to shape it into correct dimensions,
normalize the amplitude with $\angle\cdot$ and quantize it such that
$\mathbf{W}_{b}=\mathbb{Q}(\angle\mathbf{W}_{b})\in\mathcal{P}_{b}$.

For the top layer, the optimal stream power allocation can be included while
searching the multiplier $\psi_{b}$ to satisfy the sum-power constraint
$p_{b}$. Assuming the multiplier $\psi_{b}$ to be fixed, which is captured in
$\mathbf{Z}_{j_{b}}^{2}$, the power optimization problem $\forall
j_{b}\in\mathcal{D}_{b}$ can be stated as

$\displaystyle\underset{\mathbf{P}_{j_{b}}}{\max}$ $\displaystyle
w_{j_{b}}\mbox{lndet}(\mathbf{I}+\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{Z}_{j_{b}}^{1}\mathbf{W}_{b}\mathbf{V}_{j_{b}}\mathbf{P}_{j_{b}})-\mbox{Tr}(\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{Z}_{j_{b}}^{2}\mathbf{W}_{b}\mathbf{V}_{j_{b}}\mathbf{P}_{j_{b}}).$
(48)

In (48), the update of power matrix $\mathbf{P}_{j_{b}},\forall j_{b}$,
remains independent and the multiplier $\psi_{b}$ must be updated based on the
sum of the transmit covariance matrices
$\mbox{Tr}(\sum_{j_{b}}{\mathbf{W}_{b}}^{H}\mathbf{V}_{j_{b}}^{H}\mathbf{P}_{j_{b}}\mathbf{V}_{j_{b}}\mathbf{W}_{b})$,
once all the power matrices $\mathbf{P}_{j_{b}}$ are updated in parallel.
Solving $\eqref{power_DL_parallel_problem}$ in parallel $\forall j_{b}$ leads
to the following optimal power allocation scheme

$\displaystyle\mathbf{P}_{j_{b}}=($ $\displaystyle
w_{j_{b}}(\mathbf{V}_{j_{b}}^{H}\mathbf{W}_{b}^{H}\mathbf{Z}_{j_{b}}^{1}\mathbf{W}_{b}\mathbf{V}_{j_{b}})^{-1}-(\mathbf{V}_{j_{b}}^{H}\mathbf{W}_{b}^{H}\mathbf{Z}_{j_{b}}^{2}\mathbf{W}_{b}\mathbf{V}_{j_{b}})^{-1})^{+}.$
(49)

As a final step, the multiplier $\psi_{b}$ can be searched with the Bisection
method similar to (40), and while doing so, the optimal power allocation for
the DL user in the set $\mathcal{D}_{b}$ can be computed in parallel according
to (49).

Being a cooperative approach, after each synchronized iteration, each FD BS
must exchange information about the updated digital beamformers, analog
beamformer and combiner such that the neighbouring BSs could update their
local variables, function of the gradients which depends on the beamformers,
in a synchronized manner. Continuous sharing of the information makes each FD
BS aware of the generated interference at each iteration and hence enables
correct adaption of the beamformers during the optimization. It is shown in
Section VI that due to correct adaptation, P$\&$D-HYBF suffers negligible
performance loss compared to C-HYBF. Based on the per-link decomposition
proposed above, the complete procedure to execute the cooperative P$\&$D-HYBF
for WSR maximization in a multicell mmWave FD system is given in Algorithm
$2$.

Algorithm 2 Parallel and Distributed Hybrid Beamforming

Given: The rate weights, CSI and multiple processors in UL and DL $\forall b$.
Initialize: $\mathbf{W}_{b},\mathbf{V}_{j_{b}},\mathbf{U}_{k_{b}},\forall
j_{b}\in\mathcal{D}_{b},\forall k_{b}\in\mathcal{U}_{b}$ in each cell.
Repeat until convergence

$\forall b\in\mathcal{B}$, share $\mathbf{W}_{b},\mathbf{F}_{b}$ and
$\mathbf{U}_{k_{b}},\mathbf{V}_{j_{b}},\forall k_{b},\forall j_{b}$ with the
neighbouring FD BSs.

In parallel $\forall b$ ($\forall k_{b}\in\mathcal{U}_{b},\forall
j_{b}\in\mathcal{D}_{b}$)

Update $\overline{\mathbf{L}_{j_{b}}^{In}},\overline{\mathbf{L}_{k_{b}}^{In}}$
and $\overline{\mathbf{L}_{j_{b}}^{Out}}$ and
$\overline{\mathbf{L}_{k_{b}}^{Out}}$ in (34) , from the memory and based on
the feedback, respectively.

Update $\overline{\mathbf{R}_{\overline{k_{b}}}^{-1}}$ and
$\overline{\mathbf{R}_{\overline{j_{b}}}}$ in (33) from the memory.

Solve in parallel $\forall b\in\mathcal{B}$

Parallel DL for FD BS $b$

Set:  $\underline{\psi_{b}}=0,\overline{\psi_{b}}=\psi_{b}^{max}$.

Bottom layer: Compute $\mathbf{V}_{j_{b}},\forall j_{b}$ in parallel with (45)
and normalize it to unit-norm columns

Update $\overline{\mathbf{R}_{\overline{j_{b}}}^{-1}}$ and
$\overline{\mathbf{L}_{j_{b}}^{In}},\forall j_{b}$ from the memory

Middle layer: Compute $\mathbf{W}_{b}$ with (47), do unvec and get
$\angle\mathbf{W}_{b}$

Top layer: Repeat until convergence

set $\psi_{b}=(\underline{\psi_{b}}+\overline{\psi_{b}})/2$

In parallel $\forall j_{b}$

Compute $\mathbf{P}_{j_{b}}$ (49), do SVD, set
$\mathbf{P}_{j_{b}}=\mathbf{D}_{\mathbf{P}_{j_{b}}}$

Set
$\mathbf{Q}_{j_{b}}=\mathbf{W}_{b}\mathbf{V}_{j_{b}}\mathbf{P}_{j_{b}}\mathbf{V}_{j_{b}}^{H}\mathbf{W}_{b}^{H}$

if constraint for $\psi_{b}$ is violated

set $\underline{\psi_{b}}=\psi_{b}$ ,

else $\overline{\psi_{b}}=\psi_{b}$

Parallel UL for FD BS $b$

Set:
$\underline{\lambda_{k_{b}}}=0,\overline{\lambda_{k_{b}}}=\lambda_{k_{b}}^{max},\;\forall
k_{b}$.

Bottom layer: Compute $\mathbf{U}_{k_{b}},\forall k_{b},$ (37) in parallel and
normalize it

Middle layer: Repeat until convergence in parallel $\forall k_{b}$

set
$\lambda_{k_{b}}=(\underline{\lambda_{k_{b}}}+\overline{\lambda_{k_{b}}})/2$

Compute $\mathbf{P}_{k_{b}}$ in parallel with (39), do SVD and set
$\mathbf{P}_{k_{b}}=\mathbf{D}_{P_{k_{b}}}$

Set
$\mathbf{T}_{{}_{b}}=\mathbf{U}_{k_{b}}\mathbf{P}_{k_{b}}\mathbf{U}_{k_{b}}^{H}$

if constraint for $\lambda_{k_{b}}$ is violated

set $\underline{\lambda_{k_{b}}}=\lambda_{k_{b}}$

else $\overline{\lambda_{k_{b}}}=\lambda_{k_{b}}$

Update
$\overline{\mathbf{R}_{k_{b}}^{a}},\overline{\mathbf{R}_{\overline{k_{b}}}^{a}}$,
$\forall k_{b}$.

Top layer: Compute $\mathbf{F}_{b}$ with (43) and get $\angle\mathbf{F}_{b}$.

Repeat

Quantize $\mathbf{W}_{b}$ and $\mathbf{F}_{b},$ with
$\mathbb{Q}_{b}(\cdot),\forall b$.

_Remark 1:_ Recall that C-HYBF would require complete CSI to the CN every CCT
to optimize the beamformers. Then the optimization takes place jointly,
leading to significant computational complexity which is analyzed in Section
VI, and then the optimized variables need to be communicated back to each FD
BS, which then has to communicate them to its users at the millisecond scale.
The presented P$\&$D-HYBF requires that each FD BS is only aware of the local
CSI, i.e., the interfering links appearing in the gradients, which results in
significant reduction in the communication overhead. Moreover, as the
computations are fully distributed in parallel on different multi-processors
FD BSs, which reduces the computational complexity considerably and leading to
faster convergence.

_Remark 2:_ A feedback from the neighbouring FD BSs is required only to update
the out-cell local variables $\overline{\mathbf{L}_{j_{b}}^{Out}},\forall
j_{b}$ and $\overline{\mathbf{L}_{k_{b}}^{Out}},\forall k_{b}$. To further
reduce the overhead, if the interference and CI channels among different cells
vary slowly, the FD BSs can consider updating them only when these channels
have changed significantly compared to the last feedback stage.

###  V-C On the Convergence of P$\&$D-HYBF

The convergence proof for P$\&$D-HYBF follows similarly from the proof stated
for the C-HYBF scheme. Compared to C-HYBF scheme, the local variables have a
different type of information saved for each communication link which dictate
the gradients. Computing the beamformers given the information shared from the
neighbouring BSs as the GDEs increase the WSR at each iteration for each link.
However, the increase is different compared to C-HYBF as the local variables’
information differs, thus resulting for the P$\&$D-HYBF to converge to a
different local optimum.

###  V-D Computational Complexity Analysis

In this section, we present the per-iteration computational complexity for the
C-HYBF and P$\&$D-HYBF schemes. For such purpose, equal number of users in DL
and UL in each cell, i.e., $D_{b}=D$ and $U_{b}=U$, $\forall b\in\mathcal{B}$
with $|\mathcal{B}|=B$, is assumed. Moreover, the number of antennas in each
cell for the FD BSs, UL and DL users is also assumed to be the same.

Let us consider only the computational complexity of C-HYBF denoted with
$\mathcal{C}_{C}$, by ignoring its huge communication overhead to transfer
complete CSI every CCT and communicating back the optimized variables. Its one
iteration consists in updating $BD$ and $BU$ digital beamformers for the DL
and UL users, respectively, and $B$ analog beamformers and $B$ analog
combiners by the CN. Optimization of these many variables result into the
following computational complexity

$\displaystyle\mathcal{C}_{C}=$
$\displaystyle\mathcal{O}(B^{2}U^{2}{N_{b}^{RF}}^{3}+B^{2}UDN_{j_{b}}^{3}+B^{2}D^{2}N_{j_{b}}^{3}+B^{2}DU{N_{b}^{RF}}^{3}+B{M_{b}^{RF}}^{2}M_{b}^{2}+B{N_{b}^{RF}}N_{b}^{2}$
(50) $\displaystyle+BDd_{j_{b}}{M_{b}^{RF}}^{2}+BDd_{k_{b}}N_{k_{b}}^{2}).$

For P$\&$D-HYBF, different computational processors have different
computational burden and therefore we consider comparing its worst-case
complexity. Namely, let $Z_{DL}$ and $Z_{UL}$ denote the number of processors
dedicated for DL and UL by each FD BS, respectively. It is necessary to
distinguish two case: 1) $Z_{DL}=D$, $Z_{UL}=U$, and 2) $Z_{DL}<D$,
$Z_{UL}<U,\forall b$. The first case considers fully parallel implementation
with the number of processors equal to the number of users. For such a case,
the worst-case complexity in UL and DL is given for the processors which makes
one update of the digital beamformer in UL and the analog combiner, and one
update of the digital beamformer in DL and the analog beamformer,
respectively, in each cell. Let $\mathcal{C}_{P\&D}^{DL}$ and
$\mathcal{C}_{P\&D}^{UL}$ denote the worst-case complexity in DL and UL for
P$\&$D-HYBF, respectively. For the fully parallel implementation (case 1), the
complexity results to be

$\displaystyle\mathcal{C}_{P\&D}^{DL}=\mathcal{O}(BDN_{j_{b}}^{3}+BU{N_{b}^{RF}}^{3}+d_{j_{b}}{M_{b}^{RF}}^{2}+{M_{b}^{RF}}^{2}M_{b}^{2}),$
(51a)
$\mathcal{C}_{P\&D}^{UL}=\mathcal{O}(BU{N_{b}^{RF}}^{3}+BDN_{j_{b}}^{3}+d_{k_{b}}N_{k_{b}}^{2}+{N_{b}^{RF}}N_{b}^{2})$
(51b)

The second case considers that the number of processors dedicated is less than
the number of users. In such a case, each processor may have to update $K$ and
$N$ digital beamformers in DL and UL, respectively, before updating the analog
part. In such a case, the worst-case complexity in DL and UL is given for the
processors which update $K$ digital beamformers and the analog beamformer, and
$N$ digital beamformers and the analog combiner, respectively, in each cell.
In such a case, the worst case complexity in UL and DL for P$\&$D-HYBF result
to be

$\displaystyle\mathcal{C}_{P\&D}^{DL}=\mathcal{O}(KBDN_{j_{b}}^{3}+KBU{N_{b}^{RF}}^{3}+Kd_{j_{b}}{M_{b}^{RF}}^{2}+{M_{b}^{RF}}^{2}M_{b}^{2})$
(52a)
$\mathcal{C}_{P\&D}^{UL}=\mathcal{O}(NBU{N_{b}^{RF}}^{3}+NBDN_{j_{b}}^{3}+Nd_{k_{b}}N_{k_{b}}^{2}+{N_{b}^{RF}}N_{b}^{2})$
(52b)

We can see that the worst-case complexity of P$\&$D-HYBF is significantly less
compared to C-HYBF, scaling only linearly as a function of network size $B$
and density, i.e., $D$ and $U$. On the other hand, C-HYBF not only requires
significant computational complexity, but it also scales quadratically with
the network size and density. Note that C-HYBF requires significant additional
overhead to transfer full CSI to the CN, which is ignored in the computational
analysis but makes C-HYBF even more undesirable. Finally, we also remark that
higher computational complexity leads to higher power consumption at the CN
for processing [29]. P$\&$D-HYBF can enable the deployment of low-cost multi-
processor FD BSs, and processors will require minimal power consumption for
optimization.

##  VI Simulation Results

This section presents simulation results to evaluate the performance of the
proposed C-HYBF and P$\&$D-HYBF schemes. For comparison, we consider the
following benchmark schemes:

  * •

A centralized _Fully Digital FD_ scheme with the LDR noise.

  * •

A centralized _Fully Digital HD_ scheme with LDR noise, serving the UL and DL
users by separating the resources in times. It is neither affected by the SI
nor by the CI.

To compare the performance with a fully digital HD system, we define the
additional gain in terms of percentage for an FD system over an HD system as

$Gain=\frac{\mbox{WSR}_{FD}-\mbox{WSR}_{HD}}{\mbox{WSR}_{HD}}\times 100[\%],$
(53)

where $\mbox{WSR}_{FD}$ and $\mbox{WSR}_{HD}$ are the network WSR for the FD
and HD system, respectively. We assume the same SNR level for all the FD BSs,
defined as $\mbox{SNR}=p_{b}/\sigma_{b}^{2}$, with transmit power $p_{b}$ and
thermal noise variance $\sigma_{b}^{2}$. We assume that the UL users and FD
BSs transmit the same amount of power, i.e., $p_{k_{b}}=p_{b},\forall k_{b}$.
The thermal noise level for DL users is set as
$\sigma_{j_{b}}^{2}=\sigma_{b}^{2},\forall j_{b}$. The total transmit power is
normalized to $1$, and we choose the thermal noise variance to meet the
desired SNR. We simulate a multicell network of $B=2$ cells, with each FD BS
serving one DL and one UL user. P$\&$D-HYBF is evaluated on a computer
consisting of $4$ computational processors, equal to the number of users in
the network, i.e., fully parallel implementation, by using the parfor function
in Matlab. The FD BSs are assumed to have $M_{b}=100$ transmit and $N_{b}=60$
receive antennas, and the number of RF chains in transmission and reception is
chosen as $M_{b}^{RF}=N_{b}^{RF}=32,16,10$ or $8$. The phase-shifters are
assumed to be quantized with a uniform quantizer $\mathbb{Q}_{b}(\cdot)$ of
$10$ or $4$ bits. The DL and UL users are assumed to have
$N_{j_{b}}=N_{k_{b}}=5$ antennas and are served with $d_{j_{b}}=d_{k_{b}}=2$
data streams. The number of paths for each device is chosen to be
$N_{k_{b}}^{p}=3$, and the AOA $\theta_{k_{b}}^{n}$ and AOD $\phi_{k_{b}}^{n}$
are assumed to be uniformly distributed in the interval
$\mathcal{U}\sim[-30^{\circ},30^{\circ}],\forall j_{b},k_{b}$. We assume
uniform-linear arrays (ULAs) for the FD BSs and users. For the FD BSs, the
transmit and the receive array are assumed to be separated with distance
$D_{b}=20\leavevmode\nobreak\ $cm with a relative angle
$\Theta_{b}=90^{\circ}$ and $r_{m,n}$ in (5) is set given $D_{b}$ and
$\Theta_{b}$. The Rician factor is chosen to be $\kappa_{b}=1$ and the rate
weights are set to be $w_{k_{b}}=w_{j_{b}}=1$. Digital beamformers are
initialized as the dominant eigenvectors of the channel covariance matrices of
each user. The analog beamformers and combiners are initialized as the
dominant eigenvectors of the sum of the channel covariance matrices across all
the DL and UL users, respectively. The results reported herein are averaged
over $100$ channel realizations.

Figure 5: Typical convergence behaviour of the proposed schemes.

Figure 6: Execution time for the C-HYBF and the P$\&$D HYBF schemes with
$32$-RF chains.

Figure 7: Average WSR as a function of the LDR noise with
SNR$=20\leavevmode\nobreak\ $dB.

Figure 8: Average WSR as a function of the LDR noise with
SNR$=20\leavevmode\nobreak\ $dB.

Fig. 6 shows the typical convergence behaviour of the proposed schemes
compared to the fully digital FD system. It is evident that the proposed
schemes experience a different increase in the WSR at each iteration and
converge to different local optimum points, but the difference in the achieved
WSR is negligible. We also remark that as the P$\&$D-HYBF converges in a few
iterations, a minimal amount of communication overhead to share information
about the updated variables is required. Fig. 6 shows the typical execution
time to run the C-HYBF and the P$\&$D-HYBF with $32$ RF chains. We can see
that the former requires significant computational time as it can only update
different variables iteratively based on alternating optimization, one after
the other. Transferring full CSI to the CN and communicating back the
optimized variables to all cells will add significant additional time, which
is ignored in the simulations. For the latter, computation of the local
variables takes place in parallel at each FD BS which has to compute only the
variables associated with its users based on the feedback on multiple
processors. It is evident that the P$\&$D-HYBF requires $\sim 1/21$ and $\sim
1/2.3$ less time in UL and DL, respectively than the average execution time of
C-HYBF. The complexity of P$\&$D-HYBF in DL is dominated by the computation of
one large GDE to update the vectorized analog beamformer, which has complexity
$\mathcal{O}({M_{b}^{RF}}^{2}M_{b}^{2})$. In UL, complexity is dominated by
the computation of the analog combiner, which is only
$\mathcal{O}({N_{b}^{RF}}N_{b}^{2})$. The execution time to solve one sub-
problem for the bottom layers in UL and DL is also shown, which is negligible
compared to the average execution time of C-HYBF. From the complexity analysis
reported above, we can expect the execution time of C-HYBF to scale
quadratically as the number of users or cells increases. On the other hand,
P$\&$D-HYBF requires significantly less execution time and is expected to
increase only linearly as the network grows.

Figure 9: Average WSR as a function of the SNR with LDR noise
$\kappa_{k_{b}}=-80\leavevmode\nobreak\ $dB.

Figure 10: Average WSR as a function of the SNR with LDR noise
$\kappa_{k_{b}}=-80\leavevmode\nobreak\ $dB.

Fig. 8 shows the average WSR achieved with both schemes as a function of the
LDR noise with $32$ or $16$ RF chains and $10$ or $4$ bits phase-resolution.
It is visible that the P$\&$D-HYBF performs very close to the C-HYBF scheme
with the same number of RF chains and phase resolution. Fully digital FD
achieves $\sim 83\%$ of additional gain than the fully digital HD for any LDR
noise level. For a low LDR noise level $k_{b}<-80$ dB, C-HYBF and P$\&$D-HYBF
with $32$ RF chains achieve $\sim 74,55$ and $\sim 71,54\%$ additional gain
and with $16$ RF chains they achieve $\sim 67\%,48\%$ and $\sim 64\%,47\%$
additional gain with $10,4$ bits phase-resolution, respectively. We can also
see that when the LDR noise variance increases, the achievable WSR for both
the FD and HD systems decreases considerably. Fig. 8 shows the average WSR as
function of the LDR noise with only $12$ or $10$ RF chains and with $10$ or
$4$ bits phase-resolution. For LDR noise $k_{b}\leq 80\leavevmode\nobreak\
$dB, C-HYBF and P$\&$D-HYBF with $12$ RF chains achieve $\sim 60,43\%$ and
$\sim 57,43\%$ additional gain and with $10$ RF chains they achieve additional
gain of $\sim 58,38\%$ and $\sim 53,37\%$ with $10,4$ bit phase-resolution,
respectively. Fig. 8-8 exhibits the detrimental effect of the LDR noise on the
maximum achievable WSR with HYBF for both the schemes for mmWave FD and
motivates to deploy RF circuitry, which generates small LDR noise. It is to be
noted that P$\&$D-HYBF achieves similar performance as the C-HYBF scheme at
any LDR noise level.

Fig. 10 shows the average WSR as a function of the SNR with $32$ and
$16\leavevmode\nobreak\ $RF chains and with $10$ or $4$ bit phase-resolution
affected with LDR noise $k_{b}=-80\leavevmode\nobreak\ $dB, in comparison with
the benchmark schemes. We can see that a fully digital FD system achieves
$\sim 94\%$ and $\sim 82\%$ additional gain at low and high SNR, respectively.
With $32$ RF chains and $10$ bit phase-resolution, the C-HYBF scheme achieves
$\sim 79\%$ gain at all the SNR levels and the P$\&$D-HYBF achieves $\sim
77\%$ and $\sim 68\%$ gain at low and high SNR, respectively. As the phase-
resolution decreases to $4$-bits, we can see that the loss in WSR compared to
the $10$-bit phase-resolution case is much more evident at high SNR. Still,
with $16$ RF chains and $10$ or $4$ bit phase-resolution, both schemes
significantly outperform the fully digital HD scheme for any SNR level. Fig.
10 shows the average WSR as a function of the SNR with same LDR noise level as
in Fig. 10, i.e., $k_{b}=-80\leavevmode\nobreak\ $dB, but with $10$ or $12$ RF
chains and $10$ or $4$ bit phase-resolution. The achieved average WSR presents
a similar behaviour as in the case of a high number of RF chains, and it is
visible that the proposed schemes significantly outperform the fully digital
HD system also with very low number of RF chains and phase-resolution.
Moreover, P$\&$D-HYBF achieves similar performance as the C-HYBF scheme
regardless of the phase resolution and number of RF chains.

Figure 11: Average WSR as function of the SNR with LDR noise
$\kappa_{k_{b}}=-40\leavevmode\nobreak\ $dB.

Figure 12: Average WSR as a function of the SNR with LDR noise
$\kappa_{k_{b}}=-40\leavevmode\nobreak\ $dB.

Fig. 12 shows the achieved average WSR as a function of the SNR with LDR noise
$k_{b}=-40\leavevmode\nobreak\ $dB, which reflects highly non-ideal RF
circuitry. We can see that when the LDR noise dominates, reduction of the
thermal noise variance has negligible effect on the effective signal-to-LDR-
plus-thermal-noise ratio (SLNR). Therefore, dominance of the LDR noise
variance acts as a ceiling to the effective SLNR ratio which limits the
achievable WSR, tending to saturate at SNR$=10$ dB. We can also see that with
a large LDR noise level, C-HYBF and P$\&$D-HYBF still perform similarly with
the same phase-resolution and RF chains. At high SNR, both schemes achieve
higher WSR with $16$ RF chains and $10$ bit phase-resolution than the case of
$32$ RF chains and $4$ bit phase-resolution. Finally, Fig. 12 shows the
average WSR as a function of the SNR with only $10$ or $12$ RF chains and with
$10$ or $4$ bit phase-resolution, clearly showing that the proposed schemes
can still significantly outperform the full digital HD system. Fig. 12 also
shows that both schemes with $10$ RF chains and $10$ bit phase resolution are
more robust to the LDR noise than the case of $12$ RF chains and $4$ bit
phase-resolution, thus motivating to deploy high resolution phase-shifters at
the analog stage rather than deploying more number of RF chains, which not
only requires higher cost but also it will not yield higher gains than the
high resolution phase shifters case.

##  VII Conclusion

This article presented two HYBF schemes for WSR maximization in multicell
mmWave mMIMO FD systems. Firstly, a C-HYBF scheme based on alternating
optimization is presented. However, C-HYBF requires massive communication
overhead to exchange information between the FD network and the CN. Moreover,
very high computational power is required to optimize numerous variables
jointly. To overcome these drawbacks, a very low-complexity P$\&$D-HYBF design
is proposed, which enables each FD BS to solve its local per-link independent
sub-problems simultaneously on different computational processors, which
drastically reduces the communication overhead. Its complexity scales only
linearly as a function of the network size, making it highly scalable and
enabling the deployment of low-cost computational processors. Simulation
results show that the proposed HYBF designs achieve similar average WSR and
significantly outperform the centralized fully digital HD systems with only a
few RF chains.

##  Appendix A Proof of Theorem 1

We firstly proceed with the proof by considering the minorized WSR for the
digital beamformer $\mathbf{V}_{j_{b}}$, and the proof for the digital
beamformer $\mathbf{U}_{k_{b}}$ follows similarly. For the digital beamformer
$\mathbf{V}_{j_{b}}$, when the remaining variables are fixed, the following
optimization problem can be considered

$\displaystyle\underset{\begin{subarray}{c}\mathbf{V}_{j_{b}}\end{subarray}}{\max}$
$\displaystyle[w_{j_{b}}\mbox{lndet}(\mathbf{I}+\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{1}{\mathbf{W}_{b}}\mathbf{V}_{j_{b}})-\mbox{Tr}(\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{2}{\mathbf{W}_{b}}\mathbf{V}_{j_{b}})].$
(54)

To prove the result stated in Theorem 1, (56) must be simplified such that the
Hadamard’s inequality applies as in Proposition 1 [30] or Theorem 1 [31]. The
Cholesky decomposition of the matrix
$({\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{2}\mathbf{W}_{b})$ can be
written as $\mathbf{L}_{j_{b}}\mathbf{L}_{j_{b}}^{H}$ where
$\mathbf{L}_{j_{b}}$ is the lower triangular Cholesky factor. We can define
$\widetilde{\mathbf{V}_{j_{b}}}=\mathbf{L}_{j_{b}}^{H}\mathbf{V}_{j_{b}}$,
which allows to write (56) as

$\displaystyle\underset{\begin{subarray}{c}\mathbf{V}_{j_{b}}\end{subarray}}{\max}$
$\displaystyle[w_{j_{b}}\mbox{lndet}(\mathbf{I}+\widetilde{\mathbf{V}_{j_{b}}}^{H}\mathbf{L}_{j_{b}}^{-1}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{1}{\mathbf{W}_{b}}\mathbf{L}_{j_{b}}^{-H}\widetilde{\mathbf{V}_{j_{b}}})-\mbox{Tr}(\widetilde{\mathbf{V}_{j_{b}}}^{H}\widetilde{\mathbf{V}_{j_{b}}})].$
(55)

Let $\mathbf{A}_{j_{b}}\mathbf{\Lambda}_{j_{b}}\mathbf{A}_{j_{b}}^{H}$ be the
eigen decomposition of
$\mathbf{L}_{j_{b}}^{-1}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{1}{\mathbf{W}_{b}}\mathbf{L}_{j_{b}}^{-H}$,
where $\mathbf{A}_{j_{b}}$ and $\mathbf{\Lambda}_{j_{b}}$ are the unitary and
diagonal power matrices. Let
$\mathbf{D}_{j_{b}}=\mathbf{A}_{j_{b}}^{H}\widetilde{\mathbf{V}_{j_{b}}}^{H}\widetilde{\mathbf{V}_{j_{b}}}\mathbf{A}_{j_{b}}$,
and the problem (55) can be restated as

$\displaystyle\underset{\begin{subarray}{c}\mathbf{D}_{j_{b}}\end{subarray}}{\max}$
$\displaystyle[w_{j_{b}}\mbox{lndet}(\mathbf{I}+\mathbf{D}_{j_{b}}\mathbf{\Lambda}_{j_{b}})-\mbox{Tr}(\mathbf{D}_{j_{b}})].$
(56)

By using the Hadamard’s inequality [Page 233 [32]], it can be easily shown
that the optimal $\bm{D}_{j_{b}}$ must be diagonal. Let
$\mathbf{V}_{j_{b}}=\mathbf{L}_{j_{b}}^{-H}\mathbf{A}_{j_{b}}\mathbf{D}_{j_{b}}^{\frac{1}{2}}$
and we can write

${\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{1}{\mathbf{W}_{b}}=\mathbf{L}_{j_{b}}\mathbf{L}_{j_{b}}^{H}\mathbf{L}_{j_{b}}^{-H}\mathbf{A}_{j_{b}}\mathbf{D}_{j_{b}}^{\frac{1}{2}}\mathbf{\Lambda}_{j_{b}}={\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{2}{\mathbf{W}_{b}}\mathbf{V}_{j_{b}}\mathbf{\Lambda}_{j_{b}}$
(57)

which proves the result. A similar proof can be followed to prove also that
the optimal digital beamformer $\mathbf{U}_{k_{b}}$ when the remaining
variable are fixed must be a GDE.

##  Appendix B Proof of Theorem 2

To optimize the analog beamformer $\mathbf{W}_{b}$ given the remaining
variables to be fixed, we consider the following optimization problem

$\displaystyle\underset{\begin{subarray}{c}\mathbf{W}_{b}\end{subarray}}{\max}\sum_{j_{b}\in\mathcal{D}_{b}}$
$\displaystyle[w_{j_{b}}\mbox{lndet}(\mathbf{I}+\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{1}\mathbf{W}_{b}\mathbf{V}_{j_{b}})-\mbox{Tr}(\mathbf{V}_{j_{b}}^{H}{\mathbf{W}_{b}}^{H}\mathbf{\Sigma}_{j_{b}}^{2}\mathbf{W}_{b}\mathbf{V}_{j_{b}})].$
(58)

To prove the result for the analog beamformer, note that a similar reasoning
as for the digital beamformer in Appendix A does not apply directly. This is
due to the fact that the KKT condition for the analog beamformer (19) have the
form $\bm{A}_{1}\mathbf{W}_{b}\bm{A}_{2}=\bm{B}_{1}\mathbf{W}_{b}\bm{B}_{2}$
and thus not resolvable for $\mathbf{W}_{b}$. To solve it for the analog
beamformer $\mathbf{W}_{b}$, we can apply the identity
$\mbox{vec}(\bm{A}\bm{X}\bm{B})=\bm{B}^{T}\otimes\bm{A}\mbox{vec}(\bm{X})$
[33], which allows to rewrite (58) as a function of the vectorized analog
beamformers $\mathbf{W}_{b}$, leading to a resolvable KKT condition for
$\mbox{vec}(\mathbf{W}_{b})$. A similar proof from Appendix A, as done in
Theorem 3 [15], can be applied to the vectorized analog beamformer
$\mathbf{W}_{b}$, and show that the optimal vectorized analog beamformer is a
GDE of the sum of the matrices stated in (20).

## References

  * [1] C. K. Sheemar, “Hybrid beamforming techniques for massive mimo full duplex radio systems,” Ph.D. dissertation, EURECOM, 2022.
  * [2] X. Wang, L. Kong, F. Kong, F. Qiu, M. Xia, S. Arnon, and G. Chen, “Millimeter wave communication: A comprehensive survey,” _IEEE Communications Surveys & Tutorials_, vol. 20, no. 3, pp. 1616–1653, 2018.
  * [3] C. B. Barneto, S. D. Liyanaarachchi, M. Heino, T. Riihonen, and M. Valkama, “Full duplex radio/radar technology: The enabler for advanced joint communication and sensing,” _IEEE Wireless Commun._ , vol. 28, no. 1, pp. 82–88, Feb. 2021.
  * [4] A. C. Cirik, R. Wang, Y. Hua, and M. Latva-aho, “Weighted sum-rate maximization for full-duplex mimo interference channels,” _IEEE Transactions on communications_ , vol. 63, no. 3, pp. 801–815, 2015.
  * [5] B. P. Day, A. R. Margetts, D. W. Bliss, and P. Schniter, “Full-duplex bidirectional MIMO: Achievable rates under limited dynamic range,” _IEEE Trans. Signal Process._ , vol. 60, no. 7, pp. 3702–3713, Apr. 2012\.
  * [6] J. Palacios, J. Rodriguez-Fernandez, and N. González-Prelcic, “Hybrid precoding and combining for full-duplex millimeter wave communication,” in _IEEE Global Commun. Conf. (GLOBECOM)_ , Dec. 2019, pp. 1–6.
  * [7] R. López-Valcarce and N. González-Prelcic, “Analog beamforming for full-duplex millimeter wave communication,” in _IEEE 16th Int. Symp. Wirel. Commun. Syst. (ISWCS)_ , Aug. 2019, pp. 687–691.
  * [8] C. K. Sheemar and D. Slock, “Hybrid beamforming and combining for millimeter wave full duplex massive MIMO interference channel,” _arXiv preprint arXiv:2108.00465_ , 2021.
  * [9] Y. Cai, K. Xu, A. Liu, M. Zhao, B. Champagne, and L. Hanzo, “Two-timescale hybrid analog-digital beamforming for mmwave full-duplex MIMO multiple-relay aided systems,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 9, pp. 2086–2103, Jun. 2020.
  * [10] Z. Luo, L. Zhao, L. Tonghui, H. Liu, and R. Zhang, “Robust hybrid precoding/combining designs for full-duplex millimeter wave relay systems,” _IEEE Trans. Veh. Technol._ , vol. 70, no. 9, pp. 9577–9582, Jul. 2021.
  * [11] C. K. Sheemar and D. Slock, “Massive MIMO mmwave full duplex relay for IAB with limited dynamic range,” in _IEEE 11th IFIP Int. Conf. on New Techn., Mob. and Sec. (NTMS)_ , Apr. 2021, pp. 1–5.
  * [12] R. López-Valcarce and N. González-Prelcic, “Beamformer design for full-duplex amplify-and-forward millimeter wave relays,” in _IEEE 16th Int. Symp. Wirel. Commun. Syst. (ISWCS)_ , Oct. 2019, pp. 86–90.
  * [13] J. M. B. da Silva, A. Sabharwal, G. Fodor, and C. Fischione, “1-bit phase shifters for large-antenna full-duplex mmwave communications,” _IEEE Trans. Wirel. Commun._ , vol. 19, no. 10, pp. 6916–6931, Jul. 2020.
  * [14] I. P. Roberts, J. G. Andrews, and S. Vishwanath, “Hybrid beamforming for millimeter wave full-duplex under limited receive dynamic range,” _arXiv preprint arXiv:2012.11647_ , 2020.
  * [15] C. K. Sheemar, C. K. Thomas, and D. Slock, “Practical hybrid beamforming for millimeter wave massive mimo full duplex with limited dynamic range,” _IEEE Open Journal of the Communications Society_ , vol. 3, pp. 127–143, 2022\.
  * [16] M. A. Islam, G. C. Alexandropoulos, and B. Smida, “Integrated sensing and communication with millimeter wave full duplex hybrid beamforming,” _arXiv preprint arXiv:2201.05240_ , 2022.
  * [17] J.-H. Lee and O.-S. Shin, “Full-duplex relay based on distributed beamforming in multiuser mimo systems,” _IEEE transactions on vehicular technology_ , vol. 62, no. 4, pp. 1855–1860, 2012.
  * [18] X. Xu, X. Chen, M. Zhao, S. Zhou, C.-Y. Chi, and J. Wang, “Power-efficient distributed beamforming for full-duplex mimo relaying networks,” _IEEE Transactions on Vehicular Technology_ , vol. 66, no. 2, pp. 1087–1103, 2016.
  * [19] J.-H. Lee and O.-S. Shin, “Distributed beamforming approach to full-duplex relay in multiuser mimo transmission,” in _2012 IEEE Wireless Communications and Networking Conference Workshops (WCNCW)_. IEEE, 2012, pp. 278–282.
  * [20] T. M. Kim, H. J. Yang, and A. J. Paulraj, “Distributed sum-rate optimization for full-duplex MIMO system under limited dynamic range,” _IEEE Signal Process. Lett._ , vol. 20, no. 6, pp. 555–558, Feb. 2013.
  * [21] O. Tervo, H. Pennanen, D. Christopoulos, S. Chatzinotas, and B. Ottersten, “Distributed optimization for coordinated beamforming in multicell multigroup multicast systems: Power minimization and sinr balancing,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 1, pp. 171–185, 2017\.
  * [22] D. P. Palomar and M. Chiang, “A tutorial on decomposition methods for network utility maximization,” _IEEE Journal on Selected Areas in Communications_ , vol. 24, no. 8, pp. 1439–1451, 2006.
  * [23] G. Scutari, F. Facchinei, and L. Lampariello, “Parallel and distributed methods for constrained nonconvex optimization—part i: Theory,” _IEEE Trans. Signal_ , vol. 65, no. 8, pp. 1929–1944, Dec. 2016.
  * [24] G. Scutari, F. Facchinei, P. Song, D. P. Palomar, and J.-S. Pang, “Decomposition by partial linearization: Parallel optimization of multi-agent systems,” _IEEE Trans. Signal_ , vol. 62, no. 3, pp. 641–656, Nov. 2013.
  * [25] Y. Yang, G. Scutari, D. P. Palomar, and M. Pesavento, “A parallel decomposition method for nonconvex stochastic multi-agent optimization problems,” _IEEE Trans. Signal_ , vol. 64, no. 11, pp. 2949–2964, Feb. 2016\.
  * [26] C. K. Sheemar, L. Badia, and S. Tomasin, “Game-theoretic mode scheduling for dynamic TDD in 5G systems,” _IEEE Commun. Lett._ , vol. 25, no. 7, pp. 2425 – 2429, Apr. 2021.
  * [27] P. Stoica and Y. Selen, “Cyclic minimizers, majorization techniques, and the expectation-maximization algorithm: A refresher,” _IEEE Signal Process. Mag._ , vol. 21, no. 1, pp. 112–114, Feb. 2004.
  * [28] S. Boyd, S. P. Boyd, and L. Vandenberghe, _Convex Optimization_. Cambridge university press, 2004.
  * [29] D. Schneider, “Could supercomputing turn to signal processors (again)?” _IEEE Spectrum_ , vol. 49, no. 10, pp. 13–14, 2012.
  * [30] S.-J. Kim and G. B. Giannakis, “Optimal resource allocation for MIMO ad hoc cognitive radio networks,” _IEEE Trans. Inf. Theory_ , vol. 57, no. 5, pp. 3117–3131, May 2011.
  * [31] D. Hoang and R. A. Iltis, “Noncooperative eigencoding for MIMO ad hoc networks,” _IEEE Trans. Signal_ , vol. 56, no. 2, pp. 865–869, 2008.
  * [32] T. M. Cover, _Elements of Information Theory_. John Wiley & Sons, 1999.
  * [33] J. R. Magnus and H. Neudecker, _Matrix Differential Calculus with Applications in Statistics and Econometrics_. John Wiley & Sons, 2019.