# RIS-Aided Angular-Based Hybrid Beamforming Design in mmWave Massive MIMO
Systems

Ibrahim Yildirim$\ast$,$\bullet$, Asil Koc$\circ$, Ertugrul Basar$\ast$, Tho
Le-Ngoc$\circ$ $\ast$CoreLab, Department of Electrical and Electronics
Engineering, Koç University, Sariyer 34450, Istanbul, Turkey
$\bullet$Faculty of Electrical and Electronics Engineering, Istanbul Technical
University, Sariyer 34469, Istanbul, Turkey.
$\circ$Department of Electrical and Computer Engineering, McGill University,
Montreal, QC, Canada
Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This paper proposes a reconfigurable intelligent surface (RIS)-aided and
angular-based hybrid beamforming (AB-HBF) technique for the millimeter wave
(mmWave) massive multiple-input multiple-output (MIMO) systems. The proposed
RIS-AB-HBF architecture consists of three stages: (i) RF beamformer, (ii)
baseband (BB) precoder/combiner, and (iii) RIS phase shift design. First, in
order to reduce the number of RF chains and the channel estimation overhead,
RF beamformers are designed based on the 3D geometry-based mmWave channel
model using slow time-varying angular parameters of the channel. Second, a BB
precoder/combiner is designed by exploiting the reduced-size effective channel
seen from the BB stages. Then, the phase shifts of the RIS are adjusted to
maximize the achievable rate of the system via the nature-inspired particle
swarm optimization (PSO) algorithm. Illustrative simulation results
demonstrate that the use of RISs in the AB-HBF systems has the potential to
provide more promising advantages in terms of reliability and flexibility in
system design.

###### Index Terms:

Reconfigurable intelligent surface, massive MIMO, millimeter-wave, hybrid
beamforming, low CSI overhead, particle swarm optimization.

##  I Introduction

In envisioned sixth-generation (6G) wireless communication systems, enormously
abundant spectrum at high frequencies is needed to be exploited in order to
support eccentric applications which impose extreme performance requirements
[1]. Although millimeter-wave (mmWave) transmission is one of the cornerstones
of fifth-generation (5G) and beyond communication systems, promising
technologies such as massive multiple-input multiple-output (MIMO) and
reconfigurable intelligent surface (RIS) have been introduced to combat the
high path loss and blockage susceptibility at these frequencies [2, 3, 4].
RISs, which are composed of a low-cost and passive reflecting elements, have
recently attracted great interest by both academia and industry to improve the
capacity and the transmission performance of wireless communication systems
[5, 6, 7]. ††This work was supported by Scientific and Technological Research
Council of Turkey (TUBITAK) under grant no 120E401.

Beamforming (i.e., precoding in transmission and combining in reception) is an
essential signal processing technique to ensure reliable communications
between the transmitter and receiver. Although the conventional MIMO systems
widely considers single-stage fully-digital beamforming, it causes two crucial
challenges for the massive MIMO systems: (i) high hardware cost due to a
single dedicated power-hungry radio frequency (RF) chain per antenna, (ii)
large channel estimation overhead size [8]. Hybrid beamforming (HBF) is
introduced as a promising solution, which interconnects analog RF-stage and
digital baseband (BB)-stage via employing limited number of RF chains [9].
Furthermore, in [10, 11, 12], angular-based HBF (AB-HBF) is developed for also
reducing the instantaneous channel estimation overhead size. Particularly, it
builds RF-stage via slow time-varying angular information and BB-stage via
reduced dimensional channel state information (CSI).

Recently, RIS-assisted HBF studies have been introduced, in which RISs and
mmWave massive MIMO systems are combined to achieve the targeted requirements
for next-generation communication systems [13, 14, 15]. In [13], the authors
improved the spectral efficiency performance by jointly optimizing the RIS
phase shifts and the hybrid precoder/combiner via a two-stage manifold
optimization-based algorithm. While the hardware cost and energy consumption
are reduced with the HBF structure proposed in [13], the direct path between
the transmitter and the receiver is not taken into account. In [14], the
authors considered an RIS-assisted HBF scheme for broadband MIMO system at
mmWave frequency by proposing a geometric mean decomposition-based beamforming
to avoid the complex bit/power allocation. Moreover, the authors jointly
designed the RIS phase shift matrix and hybrid beamformer in [15] for the RIS-
assisted hybrid mmWave MIMO systems in the presence of a direct link between
the transmitter and the receiver. Despite improvements in error performance
and achievable rate, the existing RIS-assisted HBF studies require full CSI in
hybrid beamformer design.

Figure 1: RIS-assisted hybrid mmWave massive MIMO system model.

In this work, we introduce an RIS-aided AB-HBF (RIS-AB-HBF) system requiring
low CSI overhead for mmWave massive MIMO systems. By considering the 3D
geometry-based mmWave channel model, a three-stage RIS-aided HBF system design
is proposed: (i) RF beamformers, (ii) BB precoder/combiner, and (iii) RIS
phase shift design. First, the transmit/receive RF beamformers are designed by
using the slow time-varying angle-of-departure/angle-of-arrival (AoD/AoA)
information of the channel, while the BB precoder and combiner are constructed
exploiting the reduced-size effective channel seen from the BB-stage. Then,
the RIS phase reflection matrix is designed via particle swarm optimization
(PSO) to maximize the achievable rate of massive MIMO systems. Our
comprehensive numerical results demonstrate that the RISs significantly
enhance the achievable rate performance of the AB-HBF mmWave systems by
alleviating the disruptive effects of the wireless propagation environment as
well as reducing the CSI overhead.

##  II System Model

In this section, we present the system and channel models of the RIS-assisted
HBF scheme for mmWave massive MIMO systems as well as provide their working
principles.

###  II-A System Model and Problem Formulation

We consider an RIS-assisted hybrid mmWave massive MIMO architecture shown in
Fig. 1, where there is a single RIS between the transmitter/receiver pair. As
shown in Fig.1, the transmission is carried out via an RIS and non-RIS direct
link between the transmitter (Tx) and the receiver (Rx). We assume that the
Tx, Rx and RIS are in the form of uniform planar array (UPA). Here, the real-
time controlled RIS is equipped with $M_{I}=M_{I_{x}}\times M_{I_{y}}$ dynamic
reflecting elements that can be adjusted based on channel phases while the Tx
and Rx are equipped with $M_{T}=M_{T_{x}}\times M_{T_{y}}$ and
$M_{R}=M_{R_{x}}\times M_{R_{y}}$ antennas, respectively. The number of
antennas along the $x$ and $y$-axis of the corresponding UPA are respectively
denoted as $M_{i_{x}}$ and $M_{i_{y}}$ for $i\in\left\\{I,T,R\right\\}$.
$N_{S}$ data streams are sent from the Tx via digital BB precoder (${\bf
B}_{t}$) and analog RF beamformer (${\bf F}_{t}$), which are connected with
$N_{T}$ RF chains. At the Rx, the received signal is combined via the analog
RF beamformer (${\bf F}_{r}$) and digital BB combiner (${\bf B}_{r}$), which
are connected with $N_{R}$ RF chains. Here we assume that the number of data
streams is constraint as $N_{S}\leq\min\left(N_{T},N_{R}\right)$ and the
number of RF chains are much smaller than both the number of transmit and
receive antennas to exploit the benefit of HBF systems (i.e. $N_{T}\ll M_{T}$
and $N_{R}\ll M_{R}$).

The precoded signal ${\bf s}\in\mathbb{C}^{M_{T}}$ at the Tx is expressed as:

$\displaystyle{\bf s}={{\bf{F}}_{t}}{{\bf{B}}_{t}}{\bf{x}},$ (1)

where ${\bf F}_{t}\in\mathbb{C}^{M_{T}\times N_{T}}$ represents the transmit
RF beamformer, ${{\bf B}_{t}\in\mathbb{C}^{N_{T}\times N_{S}}}$ represents the
BB precoder and ${\bf x}\in\mathbb{C}^{N_{S}}$ is the data symbol vector
satisfying $\mathbb{E}\left\\{{\bf x}{\bf x}^{H}\right\\}={\bf I}_{N_{s}}$.
The power constraint is also ensured for the transmitted power $P_{T}$ as
$\mathbb{E}\big{\\{}\big{\|}{\bf s}\big{\|}_{2}^{2}\big{\\}}\leq P_{T}$.

The direct channel between the Tx and Rx is represented by ${\bf
H}_{TR}\in\mathbb{C}^{M_{R}\times M_{T}}$, while the RIS-assisted link
consists of two components: ${\bf H}_{TI}\in\mathbb{C}^{M_{I}\times M_{T}}$ is
the channel from the Tx to RIS and ${\bf H}_{IR}\in\mathbb{C}^{M_{R}\times
M_{I}}$ is the channel from the RIS to Rx.
${\bf\Theta}=\text{diag}{(e^{j\Omega_{1}},\dots,e^{j\Omega_{M_{I}}})}\in\mathbb{C}^{M_{I}\times
M_{I}}$ is the matrix of RIS element response with $\Omega_{i}$ being the
phase shift introduced by $i$th reflecting element of the RIS for
$i=1,\dots,M_{I}$. At the receiver, the signal received from the RIS-assisted
and direct link is obtained as:

$\displaystyle{\bf{y}}=$ $\displaystyle{\bf{Hs}}+{\bf{v}}=({\bf
H}_{IR}{\bf\Theta}{\bf H}_{TI}+{\bf H}_{TR})\bf{s}+{\bf{v}}$ (2)
$\displaystyle=$ $\displaystyle{{({\bf H}_{IR}{\bf\Theta}{\bf H}_{TI}+{\bf
H}_{TR})}}{{\bf{F}}_{t}}{{\bf{B}}_{t}}{\bf{x}}+{\bf{v}},$

where ${\bf H}={\bf H}_{IR}{\bf\Theta}{\bf H}_{TI}+{\bf
H}_{TR}\in\mathbb{C}^{M_{R}\times M_{T}}$ stands for the end-to-end channel
matrix and ${{\bf{v}}\sim\mathcal{CN}\left({\bf 0},\sigma_{v}^{2}{\bf
I}_{M_{R}}\right)}$ represents the additive white Gaussian noise. By using (1)
and (2), the received signal after combining is written by:

$\displaystyle\tilde{\bf y}={\bf
B}_{r}\bm{\mathcal{H}}{\bf{B}}_{t}{\bf{x}}+{\bf B}_{r}{\bf F}_{r}{\bf{v}},$
(3)

where $\bm{\mathcal{H}}={\bf F}_{r}{{({\bf H}_{IR}{\bf\Theta}{\bf H}_{TI}+{\bf
H}_{TR})}}{{\bf{F}}_{t}}\in\mathbb{C}^{N_{r}\times N_{t}}$ is the effective
channel matrix seen from the BB-stages, ${\bf F}_{r}\in\mathbb{C}^{N_{R}\times
M_{T}}$ represents the receive RF beamformer, ${\bf
B}_{r}\in\mathbb{C}^{N_{S}\times N_{R}}$ denotes the BB combiner.

In this work, our main aim is to efficiently design the transmit RF beamformer
and BB precoder at the Tx (i.e., ${\bf F}_{t}$, ${\bf B}_{t}$), the receive RF
beamformer and BB combiner at the Rx (i.e., ${\bf F}_{r}$, ${\bf B}_{r}$), and
the phase shifts at the RIS (i.e. $\bf{\Theta}$) to maximize the achievable
rate of the system as follows:

$\displaystyle\max_{\left\\{{\bf F}_{t},{\bf F}_{r},{\bf B}_{t},{\bf
B}_{r},\bf{\Theta}\right\\}}$ $\displaystyle R=$
$\displaystyle{\log_{2}}\left|{{{{\bf{I}}_{{N_{s}}}}+{\bf R}_{{w}}^{-1}}{\bf
B}_{r}\bm{\mathcal{H}}{\bf B}_{t}{\bf B}_{t}^{H}\bm{\mathcal{H}}^{H}{\bf
B}_{r}^{H}}\right|,$ (4a) s.t. $\displaystyle~{}\bm{\mathcal{H}}={\bf
F}_{r}{\bf HF}_{r}\in\mathbb{C}^{N_{r}\times N_{t}},$ (4d)
$\displaystyle\left|{\bf{F}}_{t}\left(m,n\right)\right|=\frac{1}{\sqrt{M_{t}}},\left|{\bf{F}}_{r}\left(m,n\right)\right|=\frac{1}{\sqrt{M_{r}}},~{}\forall\left(m,n\right),$
$\displaystyle\mathbb{E}\big{\\{}{\big{\|}{\bf{x}}\big{\|}_{2}^{2}}\big{\\}}={\rm{tr}}\left({{\bf{B}}_{t}^{H}{\bf{F}}_{t}^{H}{{\bf{F}}_{t}}{{\bf{B}}_{t}}}\right)\leq{P_{T}},$

where ${\bf R}_{{w}}=\sigma_{v}^{2}{\bf B}_{r}{\bf F}_{r}{\bf F}_{r}^{H}{\bf
B}_{r}^{H}$ is the covariance matrix of the noise ${\bf{w}}={\bf B}_{r}{\bf
F}_{r}{\bf{v}}$. Since phase shifters are used in the implementation of RF
beamformer, the constraint of the constant module is ensured as in (4d).

###  II-B Channel Model

In mmWave transmission, contrary to the conventional rich-scattering
propagation environment, only a few spatial paths are experienced due to
scattering limited propagation conditions. Hence, the mmWave channel will be
sparse in the angular domain [16]. We assume that ${\bf H}_{l}$
($l\in\left\\{TR,TI,IR\right\\}$) is a sum of the contributions of $C_{l}$
scattering clusters, each of which contribute $Z_{l,c}$ ($c=1,\dots,C_{l}$)
propagation paths to the corresponding channel matrix. Each sub-channel is
assumed to have $Z_{l}=\sum\nolimits_{c=1}^{C}Z_{l,c}$ paths in total between
the terminals.

By following the UPA structure [10] and the 3D geometry-based mmWave channel
model [17], the channel matrix between the terminals is defined as follows:

$\displaystyle{\bf H}_{l}=$
$\displaystyle\sum\limits_{c=1}^{C_{l}}{\sum\limits_{p=1}^{Z_{l,c}}g_{c_{p}}^{l}}{{{\bm{\phi}}}_{l,r}}\left({{\gamma_{x,r}^{l,c_{p}}},{\gamma_{y,r}^{l,c_{p}}}}\right){\bm{\phi}}_{l,t}^{H}\big{(}{{\gamma_{x,t}^{l,c_{p}}},{\gamma_{y,t}^{l,c_{p}}}}\big{)}$
(5) $\displaystyle=$ $\displaystyle{{\bf{\Phi}}_{l}^{r}}{\bf
G}_{l}{{\bf{\Phi}}_{l}^{t}},\>l\in\left\\{TR,TI,IR\right\\},$

where $g_{c_{p}}^{l}\sim\mathcal{CN}\left(0,\frac{\beta^{l}}{P_{l}}\right)$
stands for the complex path gain of the $p^{th}$ path within the $c^{th}$
cluster of the corresponding channel with the path index of
$c_{p}=p+\sum\nolimits_{u=1}^{c-1}Z_{l,u}$, and
${\bm{\phi}_{l,r}}\big{(}{\gamma_{x,r}^{l,c_{p}},\gamma_{y,r}^{l,c_{p}}}\big{)}$
is the receive phase response vector while the transmit phase response vector
is denoted by
${\bm{\phi}_{l,t}}\big{(}{\gamma_{x,t}^{l,c_{p}},\gamma_{y,t}^{l,c_{p}}}\big{)}$.
Here, $\beta^{l}$ represents the path loss of the corresponding path, which is
inversely proportional to the terminal distances ($d_{l}$) as
$\beta^{l}\propto d_{l}^{-\eta}$ with $\eta$ being path loss exponent for
$l\in\left\\{TR,TI,IR\right\\}$. Moreover, ${\bf
G}_{l}=\textrm{diag}\left(g_{1}^{l},\cdots,g_{P}^{l}\right)\in\mathbb{C}^{Z_{l}\times
Z_{l}}$ represents the diagonal matrix including the complex path gains,
${{\bf{\Phi}}_{l}^{r}}$ and ${{\bf{\Phi}}_{l}^{t}}$ denote the receive and
transmit phase response matrices for the corresponding channel, respectively.

Then, the transmit and receive phase response vectors are defined as:

$\displaystyle{\bm{\phi}_{ij,t}}\left({{\gamma_{x}^{l},\gamma_{y}^{l}}}\right)$
$\displaystyle=\big{[}{1,{e^{j2\pi d{{\gamma_{x}^{l}}}}},\cdots,{e^{j2\pi
d\left({{M_{i_{x}}}-1}\right){{\gamma_{x}^{l}}}}}}\big{]}^{T}\vspace{-1ex}$
(6) $\displaystyle\otimes\big{[}{1,{e^{j2\pi
d{{\gamma_{y}^{l}}}}},\cdots,{e^{j2\pi
d\left({{M_{i_{y}}}-1}\right){{\gamma_{y}^{l}}}}}}\big{]}^{T},$
$\displaystyle{\bm{\phi}_{ij,r}}\left({{\gamma_{x}^{l},\gamma_{y}^{l}}}\right)$
$\displaystyle=\big{[}{1,{e^{j2\pi d{{\gamma_{x}^{l}}}}},\cdots,{e^{j2\pi
d\left({{M_{j_{x}}}-1}\right){{\gamma_{x}^{l}}}}}}\big{]}^{T}$ (7)
$\displaystyle\otimes\big{[}{1,{e^{j2\pi
d{{\gamma_{y}^{l}}}}},\cdots,{e^{j2\pi
d\left({{M_{j_{y}}}-1}\right){{\gamma_{y}^{l}}}}}}\big{]}^{T},$

where $(i,j)\in\left\\{(T,R),(T,I),(I,R)\right\\}$ and $d$ denotes the antenna
spacing normalized by wavelength. From (5), (6) and (7), the slow time-varying
transmit and receive phase response matrices for each sub-channel are
respectively written by:

$\displaystyle{{\bf{\Phi}}_{l}^{t}}=\left[{\begin{array}[]{*{20}{c}}{{\bm{\phi}_{l,t}^{H}}\big{(}{{\gamma_{x,r}^{l,1}}},{{\gamma_{y,r}^{l,1}}}\big{)}}\\\
{\vdots}\\\
{{\bm{\phi}_{l,t}^{H}}\big{(}{{\gamma_{x,r}^{l,P_{l}}}},{{\gamma_{y,r}^{l,P_{l}}}}\big{)}}\end{array}}\right],~{}{{\bf{\Phi}}_{l}^{r}}=\left[{\begin{array}[]{*{20}{c}}{{\bm{\phi}_{l,r}^{H}}\big{(}{{\gamma_{x,r}^{l,1}}},{{\gamma_{y,r}^{l,1}}}\big{)}}\\\
{\vdots}\\\
{{\bm{\phi}_{l,r}^{H}}\big{(}{{\gamma_{x,r}^{l,P_{l}}}},{{\gamma_{y,r}^{l,P_{l}}}}\big{)}}\end{array}}\right]^{H}.$
(8)

By considering the mean and spread parameters of AoA/AoD in Table I, at the
arrival side of each sub-channel $l$ for $l\in\left\\{TR,TI,IR\right\\}$, the
coefficients reflecting the elevation AoA and azimuth AoA of the $p^{th}$ path
within the $c^{th}$ cluster are respectively defined as
${\gamma_{x,r}^{l,c_{p}}}=\sin\big{(}{{\theta_{r,c_{p}}^{l}}}\big{)}\cos\big{(}{{\psi_{r,c_{p}}^{l}}}\big{)}$
and
${\gamma_{y,r}^{l,c_{p}}}=\sin\big{(}{{\theta_{r,c_{p}}^{l}}}\big{)}\sin\big{(}{{\psi_{r,c_{p}}^{l}}}\big{)}$,
where
${\theta_{r,c_{p}}^{l}}\in\big{[}{{\theta_{r,c}^{l}}-{\delta_{r,c}^{l,\theta}}},{{\theta_{r,c}^{l}}+{\delta_{r,c}^{l,\theta}}}\big{]}$
and
${\psi_{r,c_{p}}^{l}}\in\big{[}{{\psi_{r,c}^{l}}-{\delta_{r,c}^{l,\psi}}},{{\psi_{r,c}^{l}}+{\delta_{r,c}^{l,\psi}}}\big{]}$.
At the departure side of each sub-channel $l$ for
$l\in\left\\{TR,TI,IR\right\\}$, the elevation AoD and azimuth AoD of the
$p^{th}$ path within the $c^{th}$ cluster are respectively defined as
${\gamma_{x,t}^{l,c_{p}}}=\sin\big{(}{{\theta_{t,c_{p}}^{l}}}\big{)}\cos\big{(}{{\psi_{t,c_{p}}^{l}}}\big{)}$
and
${\gamma_{y,t}^{l,c_{p}}}=\sin\big{(}{{\theta_{t,c_{p}}^{l}}}\big{)}\sin\big{(}{{\psi_{t,c_{p}}^{l}}}\big{)}$
, where
${\theta_{t,c_{p}}^{l}}\in\big{[}{{\theta_{t,c}^{l}}-{\delta_{t,c}^{l,\theta}}},{{\theta_{t,c}^{l}}+{\delta_{t,c}^{l,\theta}}}\big{]}$
and
${\psi_{t,c_{p}}^{l}}\in\big{[}{{\psi_{t,c}^{l}}-{\delta_{t,c}^{l,\psi}}},{{\psi_{t,c}^{l}}+{\delta_{t,c}^{l,\psi}}}\big{]}$.

TABLE I: AoA and AoD parameters of the $l^{\text{th}}$ channel for $l\in\left\\{TR,TI,IR\right\\}$. $\theta_{r,c}^{l}$ | Mean elevation AoA | $\theta_{t,c}^{l}$ | Mean elevation AoD
---|---|---|---
$\delta_{r,c}^{l,\theta}$ | Elevation AoA spread | $\delta_{t,c}^{l,\theta}$ | Elevation AoD spread
$\psi_{r,c}^{l}$ | Mean azimuth AoA | $\psi_{t,c}^{l}$ | Mean azimuth AoD
$\delta_{r,c}^{l,\psi}$ | Azimuth AoA spread | $\delta_{t,c}^{l,\psi}$ | Azimuth AoD spread

##  III RIS-Aided Angular-Based Hybrid Beamforming

In this section, we respectively provide the design of RF beamformers, BB
precoder/combiner and phase shift matrix of the RIS for the proposed RIS-aided
AB-HBF system. Our primary goal is to maximize the achievable rate of the
system as well as reduce the CSI overhead size and the hardware complexity.
From this point of view, first, the design procedures of RF beamformers is
introduced based on slow-time varying angular parameters. Afterwards, BB
precoder/combiner is designed by the singular value decomposition (SVD) and
water filling algorithm using effective channel with reduced size. Finally,
the phase shifts of the RIS, which enhances the system capacity, is obtained
by applying the PSO method.

###  III-A RF Beamformers

In order to design RF-stage efficiently, the range of the azimuth and
elevation AoA/AoD angles of the existing paths in the clusters is first
required to be appropriately defined as in [11]. Therefore, AoA and AoD angle
supports are respectively defined based on (5) as:

$\displaystyle{\Upsilon_{\textrm{AoA}}}$ $\displaystyle=$
$\displaystyle\Big{\\{}{{{\left[{{\gamma_{x}},{\gamma_{y}}}\right]}}={{\left[\sin\left(\theta\right){\cos\left(\psi\right),\sin\left(\theta\right)\sin\left(\psi\right)}\right]}}}\hskip
0.86108pt\Big{|}\hskip
0.86108pt\theta\in{\bm{\theta}}_{r},\psi\in{\bm{\psi}}_{r}\Big{\\}},\vspace{-0.6ex}\vspace{0.5ex}$
(9a) $\displaystyle\Upsilon_{\textrm{AoD}}$ $\displaystyle=$
$\displaystyle\Big{\\{}{{{\left[{{\gamma_{x}},{\gamma_{y}}}\right]}}={{\left[\sin\left(\theta\right){\cos\left(\psi\right),\sin\left(\theta\right)\sin\left(\psi\right)}\right]}}}\hskip
0.86108pt\Big{|}\hskip
0.86108pt\theta\in{\bm{\theta}}_{t},\psi\in{\bm{\psi}}_{t}\Big{\\}},\vspace{-0.6ex}$
(9b)

where ${\bm{\theta}}_{r}={\bm{\theta}}_{r}^{TI}\cup{\bm{\theta}}_{r}^{TR}$,
${\bm{\psi}}_{r}={\bm{\psi}}_{r}^{TI}\cup{\bm{\psi}}_{r}^{TR}$,
${\bm{\theta}}_{t}={\bm{\theta}}_{t}^{IR}\cup{\bm{\theta}}_{t}^{TR}$, and
${\bm{\psi}}_{t}={\bm{\psi}}_{t}^{IR}\cup{\bm{\psi}}_{t}^{TR}$. Here, the
elevation and azimuth AoA supports for the corresponding sub-channels are
respectively defined as:
${\bm{\theta}}_{r}^{l}=\cup_{c=1}^{C_{l}}\big{[}{{\theta_{r,c}^{l}}-{\delta_{r,c}^{l,\theta}}},{{\theta_{r,c}^{l}}+{\delta_{r,c}^{l,\theta}}}\big{]}$
and
${\bm{\psi}}_{r}^{l}=\cup_{c=1}^{C_{l}}\big{[}{{\psi_{r,c}^{l}}-{\delta_{r,c}^{l,\psi}}},{{\psi_{r,c}^{l}}+{\delta_{r,c}^{l,\psi}}}\big{]}$
for $l\in\left\\{IR,TR\right\\}$. Then, the elevation and azimuth AoD supports
for the corresponding sub-channels are respectively obtained as:
${\bm{\theta}}_{t}^{l}=\cup_{c=1}^{C_{l}}\big{[}{{\theta_{t,c}^{l}}-{\delta_{t,c}^{l,\theta}}},{{\theta_{t,c}^{l}}+{\delta_{t,c}^{l,\theta}}}\big{]}$
and
${\bm{\psi}}_{t}^{l}=\cup_{c=1}^{C_{l}}\big{[}{{\psi_{t,c}^{l}}-{\delta_{t,c}^{l,\psi}}},{{\psi_{t,c}^{l}}+{\delta_{t,c}^{l,\psi}}}\big{]}$
for $l\in\left\\{TI,TR\right\\}$.

By using (2) and (5), the effective channel seen from the BB-stages can be
reorganized as:

$\displaystyle\bm{\mathcal{H}}=$ $\displaystyle{\bf F}_{r}{\bf H}{\bf F}_{t},$
(10) $\displaystyle=$ $\displaystyle{\bf F}_{r}({{\bf{\Phi}}_{IR}^{r}}{\bf
G}_{IR}{{\bf{\Phi}}_{IR}^{t}}{\bf\Theta}{{\bf{\Phi}}_{{TI}}^{r}}{\bf
G}_{TI}{{\bf{\Phi}}_{TI}^{t}}$ $\displaystyle+{{\bf{\Phi}}_{{TR}}^{r}}{\bf
G}_{TR}{{\bf{\Phi}}_{TR}^{t}}){\bf F}_{t}.$

In order to optimize the receive (transmit) beamforming gain and take
advantage of all degrees of freedom enabled by the channel, the subspace
spanned by ${{\bf{\Phi}}_{IR}^{r}}$ and ${{\bf{\Phi}}_{TR}^{r}}$
(${{\bf{\Phi}}_{TI}^{t}}$ and ${{\bf{\Phi}}_{TR}^{t}}$) should be encompassed
in the columns of ${\bf F}_{r}$ (${\bf F}_{t}$): $\textrm{Span}\left({\bf
F}_{r}\right)\subset\left(\textrm{Span}\left({{\bf{\Phi}}_{IR}^{r}}\right)\cup\textrm{Span}\left({{\bf{\Phi}}_{TR}^{r}}\right)\right)$
and $\textrm{Span}\left({\bf
F}_{t}\right)\subset\left(\textrm{Span}\left({{\bf{\Phi}}_{TI}^{t}}\right)\cup\textrm{Span}\left({{\bf{\Phi}}_{TR}^{t}}\right)\right)$
should be satisfied. Afterwards, the columns of the RF precoder matrix are
created using the transmit steering vector
${{\bf{a}}_{t}}\left({{\gamma_{x}},{\gamma_{y}}}\right)=\frac{1}{{\sqrt{{M_{t}}}}}{\big{[}{1,{e^{j2\pi
d\gamma_{x}}},\cdots,{e^{j2\pi
d\left({{M_{T_{x}}}-1}\right)\gamma_{x}}}}\big{]}^{T}}\otimes{\big{[}{1,{e^{j2\pi
d\gamma_{y}}},\cdots,{e^{j2\pi
d\left({{M_{T_{y}}}-1}\right)\gamma_{y}}}}\big{]}^{T}}\in\mathbb{C}^{M_{T}}$.
The quantized angle pairs are defined as follows to obtain $M_{T}$ orthogonal
transmit steering vectors as the minimum number of angle pairs spanning the
entire elevation and azimuth space [10, 11]:
${{\kappa_{x,t}^{m}}=\frac{2m-1}{{{M_{T_{x}}}}}}-1$ for
$m=1,\cdots,{M_{T_{x}}}$ and
${{\kappa_{y,t}^{n}}=\frac{2n-1}{{{M_{T_{y}}}}}}-1$ for
$n=1,\cdots,{M_{T_{y}}}$. Thus, while minimizing the use of the RF chain by
covering the entire AoD support, the quantized angle pairs within the AoD
support are written as:

$\left({\kappa_{x,t}^{m},\kappa_{y,t}^{n}}\right)~{}\Big{|}\Big{.}~{}{\gamma_{x}}\in{\bm{\kappa}}_{x,t}^{m},{\gamma_{y}}\in{\bm{\kappa}}_{y,t}^{n},\left({{\gamma_{x}},{\gamma_{y}}}\right)\in\Upsilon_{\textrm{AoD}},\vspace{-1ex}$
(11)

where
${\bm{\kappa}}_{x,t}^{m}=\big{[}{\kappa_{x,t}^{m}-\frac{1}{{{M_{T_{x}}}}},\kappa_{x,t}^{m}+\frac{1}{{{M_{T_{x}}}}}}\big{]}$
and
${\bm{\kappa}}_{y,t}^{n}=\big{[}{\kappa_{y,t}^{n}-\frac{1}{{{M_{T_{y}}}}},\kappa_{y,t}^{n}+\frac{1}{{{M_{T_{y}}}}}}\big{]}$
respectively represent the boundary of $\kappa_{x,t}^{m}$ and
$\kappa_{y,t}^{n}$. Therefore, the RF precoder with $N_{T}$ quantized angle-
pairs is constructed based on (9b) as:

${{\bf{F}}_{t}}=\big{[}{{\bf{a}}_{t}\big{(}{\kappa^{m_{1}}_{x,t},\kappa^{n_{1}}_{y,t}}\big{)},\cdots,{\bf{a}}_{t}\big{(}{\kappa^{m_{N_{T}}}_{x,t},\kappa^{n_{N_{T}}}_{y,t}}\big{)}}\big{]}\in\mathbb{C}^{M_{T}\times
N_{T}}.\vspace{-1ex}$ (12)

Following the similar procedure for the RF combiner ($\bf{F}_{r}$), the
receive steering vector is obtained as:
${{\bf{a}}_{r}}\left({{\gamma_{x}},{\gamma_{y}}}\right)=\frac{1}{{\sqrt{{M_{R}}}}}{\big{[}{1,{e^{j2\pi
d\gamma_{x}}},\cdots,{e^{j2\pi
d\left({{M_{R_{x}}}-1}\right)\gamma_{x}}}}\big{]}^{T}}\otimes{\big{[}{1,{e^{j2\pi
d\gamma_{y}}},\cdots,{e^{j2\pi
d\left({{M_{R_{y}}}-1}\right)\gamma_{y}}}}\big{]}^{T}}\in\mathbb{C}^{M_{R}}$.
The quantized angle pairs are given by
${{\kappa_{x,r}^{m}}=\frac{2m-1}{{{M_{R_{x}}}}}}-1$ for
$m=1,\cdots,{M_{R_{x}}}$ and
${{\kappa_{y,r}^{n}}=\frac{2n-1}{{{M_{R_{y}}}}}}-1$ for
$n=1,\cdots,{M_{R_{y}}}$.

In order to minimize the use of the RF chain in the receiver side by covering
the entire AoA support, the quantized angle pairs within the AoA support are
obtained as:

$\left({\kappa_{x,r}^{m},\kappa_{y,r}^{n}}\right)~{}\Big{|}\Big{.}~{}{\gamma_{x}}\in{\bm{\kappa}}_{x,r}^{m},{\gamma_{y}}\in{\bm{\kappa}}_{y,r}^{n},\left({{\gamma_{x}},{\gamma_{y}}}\right)\in\Upsilon_{\textrm{AoA}},\vspace{-1ex}$
(13)

where
${\bm{\kappa}}_{x,t}^{m}=\big{[}{\kappa_{x,t}^{m}-\frac{1}{{{M_{T_{x}}}}},\kappa_{x,t}^{m}+\frac{1}{{{M_{T_{x}}}}}}\big{]}$
and
${\bm{\kappa}}_{y,r}^{n}=\big{[}{\kappa_{y,r}^{n}-\frac{1}{{{M_{R_{y}}}}},\kappa_{y,r}^{n}+\frac{1}{{{M_{R_{y}}}}}}\big{]}$
respectively denote the boundary of $\kappa_{x,r}^{m}$ and $\kappa_{y,r}^{n}$.
Finally, the RF combiner with $N_{R}$ quantized angle-pairs is constructed
fulfilling AoA support in (9a) as:

${{\bf{F}}_{r}}=\big{[}{{\bf{a}}_{r}\big{(}{\kappa^{m_{1}}_{x,r},\kappa^{n_{1}}_{y,r}}\big{)},\cdots,{\bf{a}}_{r}\big{(}{\kappa^{m_{N_{R}}}_{x,r},\kappa^{n_{N_{R}}}_{y,r}}\big{)}}\big{]}\in\mathbb{C}^{M_{R}\times
N_{R}}.\vspace{-0.5ex}$ (14)

It should be noted that the unitary property for RF beamformers (i.e. ${\bf
F}_{t}^{H}{\bf F}_{t}={\bf I}_{N_{T}}$, ${\bf F}_{r}{\bf F}_{r}^{H}={\bf
I}_{N_{R}}$) and the constant modulus condition provided in (4d) are satisfied
in the RF beamformer design stage.

###  III-B Baseband Precoder/Combiner

After the design of the RF beamformers, this subsection focuses on the BB
precoder/combiner design using the low-dimensional effective channel matrix
given in (10). The singular value decomposition (SVD) of effective channel
matrix$\bm{\mathcal{H}}$ is obtained as $\bm{\mathcal{H}}={\bf U\Sigma
V}^{H}$, where ${\bf
U}\in\mathbb{C}^{N_{R}\times\textrm{rank}\left(\bm{\mathcal{H}}\right)}$
represents a tall unitary matrix,
${\bf\Sigma}=\textrm{diag}\big{(}\sigma_{1}^{2},\cdots,\sigma_{\textrm{rank}\left(\bm{\mathcal{H}}\right)}^{2}\big{)}\in\mathbb{R}^{\textrm{rank}\left(\bm{\mathcal{H}}\right)\times\textrm{rank}\left(\bm{\mathcal{H}}\right)}$
denotes a diagonal matrix of the singular values in the descending order and
${\bf V}\in\mathbb{C}^{N_{T}\times\textrm{rank}\left(\bm{\mathcal{H}}\right)}$
stands for a tall unitary matrix. Here, ${\bf V}$ can be partitioned as ${\bf
V}=\left[{\bf V}_{1},{\bf V}_{2}\right]$ for ${\bf
V}_{1}\in\mathbb{C}^{N_{T}\times N_{S}}$ under
$\textrm{rank}\left(\bm{\mathcal{H}}\right)\geq N_{S}$.

Similar to [18, 11], the optimal solution for the BB precoder can be
represented as:

${\bf B}_{t}={\bf V}_{1}{\bf\Gamma}^{\frac{1}{2}}\in\mathbb{C}^{N_{T}\times
N_{T}},$ (15)

where
${\bf\Gamma}=\textrm{diag}\left(\Gamma_{1},\cdots,\Gamma_{N_{S}}\right)\in\mathbb{R}^{N_{S}\times
N_{S}}$ represents a transmit power matrix with
$\sum\nolimits_{n=1}^{N_{S}}\Gamma_{n}={P_{T}}$ for meeting the power
constraint given in (4d). By using the water filling power control solution
[19], the powers are allocated as:

${\Gamma_{n}}={\left({\mu-\frac{\sigma_{v}^{2}}{{P_{T}\sigma_{n}^{2}}}}\right)^{+}},{\textrm{
s.t.}}\sum\limits_{n=1}^{N_{S}}{\left({\mu-\frac{\sigma_{v}^{2}}{{P_{T}\sigma_{n}^{2}}}}\right)^{+}}=P_{T}.\vspace{-1ex}$
(16)

Then, considering the given RF beamformer (${\bf F}_{t}$/${\bf F}_{r}$), phase
response matrix (${\bf\Theta}$) and BB precoder (${\bf B}_{t}$), we strive to
design the BB combiner (${\bf B}_{r}$) that minimize the mean-squared-error
(MSE) between the transmitted and processed received signals. Based on (3),
the MSE expression can therefore be stated as:

$\displaystyle{\bf{e}_{\rm{MSE}}}=$
$\displaystyle\mathbb{E}\big{\\{}{{{\big{\|}{\tilde{\bf{y}}-{\bf{x}}}\big{\|}}^{2}}}\big{\\}}={\rm{tr}}\big{(}\mathbb{E}{\big{\\{}{\tilde{\bf{y}}{{\tilde{\bf{y}}}^{H}}-{\bf{x}}{{\tilde{\bf{y}}}^{H}}-\tilde{\bf{y}}{{\bf{x}}^{H}}+{\bf{x}}{{\bf{x}}^{H}}}\big{\\}}}\big{)}$
(17) $\displaystyle=$ $\displaystyle\rm{tr}\Big{(}{{\bf
B}_{r}\bm{\mathcal{H}}{{\bf{B}}_{t}}{{\bf{B}}_{t}^{H}}\bm{\mathcal{H}}^{H}}{\bf
B}_{r}^{H}+{\bf B}_{r}{\bf R}_{{q_{r}}}{\bf B}_{r}^{H}$
$\displaystyle-{{{\bf{B}}_{t}^{H}}\bm{\mathcal{H}}^{H}{\bf B}_{r}^{H}}+{\bf
I}_{N_{s}}\Big{.}\Big{)}.$

For the purpose of finding the optimal solution for (17), the differentiation
rules for the complex-valid matrices are applied. Then, the derivative of
$\bf{e}_{\textrm{MSE}}$ with respect to BB combiner (${\bf B}_{r}$) is
calculated as:

$\displaystyle\frac{{\partial\bf{e}_{{\rm{MS}}{{\rm{E}}{}}}}}{{\partial{{\bf{B}}_{r}}}}$
$\displaystyle=2{\bf
B}_{r}{\bm{\mathcal{H}}{{\bf{B}}_{t}}{{\bf{B}}_{t}^{H}}\bm{\mathcal{H}}^{H}}-2{{{\bf{B}}_{t}^{H}}\bm{\mathcal{H}}^{H}}.$
(18)

Finally, the BB combiner meeting the minimum MSE criterion in (17) can be
obtained as follows:

${\bf{B}}_{r}={\bf{B}}_{t}^{H}\bm{\mathcal{H}}^{H}\left(\bm{\mathcal{H}}{{\bf{B}}_{t}}{{\bf{B}}_{t}^{H}}\bm{\mathcal{H}}^{H}\right)^{-1}.\vspace{-1ex}$
(19)

###  III-C RIS Phase Shift Design

In this subsection, phase shifts of the RIS elements are adjusted by
optimizing the achievable rate of the system via the PSO algorithm. In PSO,
each particle reconfigures its position at every step by moving through
multidimensional space with its own experience and the experience of its peers
towards an optimum solution by the whole swarm [20]. Under given RF
beamformers (${\bf F}_{t},{\bf F}_{r}$) and BB precoder/combiner (${\bf
B}_{t}$, ${\bf B}_{r}$), the optimization problem for maximizing the
achievable rate can be formulated as:

$\displaystyle\max_{\left\\{\bf{\Theta}\right\\}}{}$ $\displaystyle
R\left({\bf F}_{t},{\bf F}_{r},{\bf B}_{t},{\bf B}_{r},\bf{\Theta}\right)$
(20) s.t.
$\displaystyle{\bf\Theta}=\text{diag}{(e^{j\Omega_{1}},\dots,e^{j\Omega_{M_{I}}})}$
$\displaystyle\mathbb{E}\big{\\{}{\big{\|}{\bf{x}}\big{\|}_{2}^{2}}\big{\\}}={\rm{tr}}\left({{\bf{B}}_{t}^{H}{\bf{F}}_{t}^{H}{{\bf{F}}_{t}}{{\bf{B}}_{t}}}\right)\leq{P_{T}}$
$\displaystyle\Omega_{i}\in[0,2\pi],\forall i,\vspace{-1ex}$

where $R\left({\bf F}_{t},{\bf F}_{r},{\bf B}_{t},{\bf
B}_{r},\bf{\Theta}\right)$ is given in (4a). In our PSO solution, the phase
shifts of the RIS elements ($\Omega_{i}$, $i\in{1,\dots,M_{I}}$) are modeled
as the positions of the particles, and the objective function in (20) is
considered as the fitness function of PSO.

$N_{\textrm{par}}$ particles are randomly initialized with position ${\bf
p}\in\mathbb{C}^{M_{I}}$ and velocity ${\bf v}\in\mathbb{C}^{M_{I}}$. Each
particle has a memory for its local best position (${\bf
p}_{\textrm{best},i}^{(q)}$) and the global best position (${\bf
p}_{G}^{(q)}$), where the superiority of a position is evaluated by the given
fitness function. Here, ${\bf p}_{\textrm{best},i}^{(q)}$ is the best position
that the $i^{\textrm{th}}$ particle has achieved while ${\bf
p}_{G}^{\textrm{(t)}}$ is the best position achieved by any of the particles
in the entire swarm in the $q^{\textrm{th}}$ iteration. For each iteration,
the velocity and position of each particle are updated according to:

$\displaystyle{\bf v}_{i}^{(q+1)}=$ $\displaystyle\zeta{\bf
v}_{i}^{(q)}+u_{1}\rho_{1}\big{(}{\bf p}_{\textrm{best},i}^{(q)}-{\bf
p}_{i}^{(q)}\big{)}+u_{2}\rho_{2}\big{(}{\bf p}_{G}^{(q)}-{\bf
p}_{i}^{(q)}\big{)},$ (21) $\displaystyle{\bf p}_{i}^{(q+1)}=$
$\displaystyle{\bf v}_{i}^{(q+1)}+{\bf p}_{i}^{(q)},$

where $\zeta$ is the inertia weight of velocity, $\rho_{1}$ and $\rho_{2}$ are
the acceleration constants, and $u_{1}$ and $u_{2}$ are uniformly distributed
random numbers in $[0,1]$. Through the iterations, the particles move from
${\bf p}_{i}^{(q)}$ to ${\bf p}_{i}^{(q+1)}$. Finally, after $Q$ iterations,
the RIS phase shift matrix is determined at the last iteration as
${\bf\Theta}=\textrm{diag}\big{(}\exp\big{(}{j{\bf
p}_{G}^{(Q)}}\big{)}\big{)}$. The summary of the proposed RIS-aided AB-HBF
system design including three stages is precisely explained in Algorithm 1.

Algorithm 1 Proposed RIS-aided AB-HBF System Design

0: $\bm{\mathcal{H}}$, $\Upsilon_{\textrm{AoA}}$, $\Upsilon_{\textrm{AoD}}$,
$P_{T}$, $Q$, $N_{\textrm{par}}$, $\zeta$, $\rho_{1}$, $\rho_{2}$

0: ${\bf F}_{t}$, ${\bf F}_{r}$, ${\bf B}_{t}$, ${\bf B}_{r}$, ${\bf\Theta}$

Build ${\bf F}_{t}$ and ${\bf F}_{r}$ via (12) and (14), respectively.

for $q=1:Q$ do

Calculate each particle’s velocity ${\bf v}_{i}^{(q)}$ via (21).

Calculate each particle’s position ${\bf p}_{i}^{(q)}$ via (21).

Find each particle’s achievable rate via (20).

Update personal best ${\bf p}_{\textrm{best},i}^{(q)}$ and global best ${\bf
p}_{{G}}^{(q)}$.

end for

Set RIS phase-shift matrix as
${\bf\Theta}=\textrm{diag}\big{(}\exp\big{(}{j{\bf
p}_{G}^{(Q)}}\big{)}\big{)}$.

Build ${\bf B}_{t}$ and ${\bf B}_{r}$ via (15) and (19), respectively.

##  IV Numerical Results

In this section, we provide comprehensive numerical results to investigate the
performance of the proposed RIS-aided AB-HBF system. Here, Table II presents
the common parameters for all simulations considering the geometry-based
mmWave channel given in (5). Due to the simplicity, we assume that all
subchannels have $Z_{l,c}=5$ paths within a single cluster $C_{l}=1$ for
$l\in\left\\{TR,TI,IR\right\\}$ and $c=1$. Moreover, the azimuth and elevation
AoA/AoD spreads for all subchannels are assumed to be the same as
$\delta_{i,c}^{l,j}=$$$ for $l\in\left\\{TR,TI,IR\right\\}$,
$i=\in\left\\{t,r\right\\}$, $j\in\left\\{\theta,\psi\right\\}$ and $c=1$. The
spacing between two antennas/reflectors is also $d=0.5$ due to the half-
wavelength separation. The top view of the simulation setup, including the
terminal distances, is also illustrated in Fig. 2. Here, the Tx-Rx channel is
assumed to be nearly blocked under non-line-of-sight (NLOS) conditions, while
the RIS-assisted cascade channel is assumed to have LOS conditions.

Figure 2: Top view of simulation setup for the proposed RIS-aided hybrid MIMO system. TABLE II: Simulation parameters. # of transmit antennas | $M_{T}=8\times 8=64$
---|---
# of receive antennas | $M_{R}=4\times 4=16$
# of RIS elements | $M_{I}=16\times 16=256$
# of RF Chains |  $N_{T}=6$, $N_{R}=2$
# of data streams | $N_{S}=2$
# of clusters for ${\bf H}_{l}$ |  $C_{l}=1$ for $l\in\left\\{TR,TI,IR\right\\}$
# of paths for ${\bf H}_{l}$ |  $Z_{l,c}=5$ for $l\in\left\\{TR,TI,IR\right\\}$
Noise PSD |  $-174$ dBm/Hz
Channel Bandwidth |  $10$ kHz
PSO # of iterations | $T=200$
PSO # of particles | $N_{\textrm{par}}=100$
Path loss exponent |  $\eta=2.3$ for LOS, $\eta=4.5$ for NLOS [21]
Mean elevation AoA/AoD | $\theta_{i,c}^{TR}=$$$ , $\theta_{i,c}^{TI}=$$$, $\theta_{i,c}^{IR}=$$$, $i\in\left\\{t,r\right\\}$
Mean azimuth AoA/AoD |  $\phi_{i,c}^{TR}=$$$ , $\phi_{i,c}^{TI}=$$$ , $\phi_{i,c}^{IR}=$$$, $i\in\left\\{t,r\right\\}$

Fig. 3(a) compares the achievable rate performance for the AB-HBF system with
or without RIS for varying Tx-RIS separation under fixed Tx and Rx positions.
In this comparison, the RIS-assisted AB-HBF schemes with random and constant
phase adjustment are also considered as a benchmark system for $d_{TR}=200$
and $d_{V}=5$ m. In this scenario, it is observed that RIS significantly
improves system performance when the direct link between the Tx and Rx is
nearly blocked. For $d_{1}=100$ m, the PSO-based RIS-assisted AB-HBF system
provides $25$ dB superiority to the conventional AB-HBF system, while $4.5$ dB
superiority to the random phase-adjusted RIS-AB-HBF system. Fig. 3(a) also
provides practical perspective for the effective positioning of the RIS. For
$d_{1}=20$ m, the achievable rate will be measured as $10.47$ bps/Hz, while it
will approximately reduce to one-third, measuring $3.24$ bps/Hz for
$d_{1}=100$ m. In Fig. 3(b), the achievable rate performance of the AB-HBF
system is illustrated depending on the position of Rx changing along the
$x$-axis under the fixed Tx and RIS positions. Here, the fixed distances
between the Tx and RIS are given as $d_{1}=10$ m and $d_{V}=20$ m. However,
doubling the $d_{TR}$ ($25$ m to $50$ m) will result in a performance loss of
$22.85\%$ in the PSO-based RIS-AB-HBF system, $28.2\%$ in the random phase RIS
AB-HBF, $57.03\%$ in the constant phase RIS-AB-HBF, and $76.53\%$ in the
conventional AB-HBF. From the presented results in Fig. 3, we conclude that
the presence of an RIS in AB-HBF systems alleviates the deteriorated channel
conditions by providing an additional reliable transmission link alternative
to the blocked Tx-Rx link. Thus, a remarkable increase in achievable rate
performance is achieved by effectively positioning the RIS even in the absence
of a LOS Tx-Rx channel.

Figure 3: Comparison of achievable rate performance for AB-HBF system with or
without RIS for varying (a) $d_{1}$ and (b) $d_{TR}$.

In Fig. 4, we present achievable rate comparison results for PSO-based, random
and constant phase adjustment for the RIS-aided AB-HBF system at $100$ and
$400$ m $d_{TR}$ values with increasing $P_{T}$. Here, the distances for this
setup are assumed as $d_{V}=5$ and $d_{2}=4d_{1}$. Fig. 4 gives an useful
insights for the different operating modes of an RIS. Even if the phase shifts
of the RIS are set to a fixed value, a notable increase in the system
performance can be achieved by providing an additional LOS path for the
transmission. However, when the received signal quality is below a certain
threshold value, the transmission performance can be increased by activating
the PSO-assisted phase adjustment algorithm. From this point of view, the use
of RISs in AB-HBF systems may has the potential to provide more promising
advantages in terms of reliability and flexibility in system design.

Figure 4: Achievable rate comparison of different phase adjustment methods in
RIS-aided AB-HPC systems for varying $M_{I}$.

Finally, Fig. 5 demonstrates the achievable rate comparison of the AB-HPC
system in the presence of an RIS under varying $M_{I}$ for $d_{TR}=100$,
$d_{1}=10$, $d_{2}=80$ and $d_{V}=10$ m. As observed from Fig. 5, increasing
$M_{I}$ significantly enhance the effect of the RIS in transmission both in
the case of PSO-based and random phase adjustment. While the constant phase
adjustment contributes to the transmission by providing a reliable LOS path,
it will become saturated at increasing $M_{I}$ values. Although the random
phase adjustment also leads to an remarkable increase in the system
performance, the effect of increasing $M_{I}$ is more evident for the PSO-
based RIS-AB-HBF system. For instance, the PSO-based scenario provides a $3$
dB advantage over the random phase adjustment case for $M_{I}=10$, while
provides $4.5$ dB increase for $M_{I}=100$.

##  V Conclusions

In this paper, a RIS-aided AB-HBF system has been introduced for the hybrid
mmWave massive MIMO system. In this context, we have investigated a number of
RIS-aided AB-HBF systems in terms of achievable rate performance to find the
most efficient architecture where the RIS and HBF technologies work in harmony
with each other. It has been shown via extensive numerical results that a
remarkable increase in achievable rate performance is achieved by effectively
positioning the RIS even in the absence of a LOS channel between the Tx and
Rx. The use of RISs in the AB-HBF systems may have the potential to provide
more favorable advantages in terms of reliability and flexibility in system
design.

Figure 5: Achievable rate performance of the AB-HBF system with and without
RIS for varying $M_{I}$.

## References

  * [1] W. Jiang et al., “The road towards 6G: A comprehensive survey,” _IEEE Open J. Commun. Society_ , vol. 2, pp. 334–366, 2021.
  * [2] E. G. Larsson _et al._ , “Massive MIMO for next generation wireless systems,” _IEEE Commun. Mag._ , vol. 52, no. 2, pp. 186–195, Feb. 2014.
  * [3] E. Basar _et al._ , “Reconfigurable intelligent surfaces for future wireless networks: A channel modeling perspective,” _IEEE Wireless Commun._ , vol. 28, no. 3, pp. 108–114, 2021.
  * [4] E. Basar et al., “Wireless communications through reconfigurable intelligent surfaces,” _IEEE Access_ , pp. 116 753–116 773, Sep. 2019\.
  * [5] C. Huang et al., “Reconfigurable intelligent surfaces for energy efficiency in wireless communication,” _IEEE Trans. Wireless Commun._ , Aug. 2019.
  * [6] I. Yildirim _et al._ , “Modeling and analysis of reconfigurable intelligent surfaces for indoor and outdoor applications in future wireless networks,” _IEEE Trans. Commun._ , vol. 69, no. 2, pp. 1290–1301, 2021.
  * [7] I. Yildirim _et al._ , “Hybrid RIS-empowered reflection and decode-and-forward relaying for coverage extension,” _IEEE Commun. Lett._ , vol. 25, no. 5, pp. 1692–1696, 2021.
  * [8] N. Fatema _et al._ , “Massive MIMO linear precoding: A survey,” _IEEE Syst. J._ , vol. 12, no. 4, pp. 3920–3931, Dec. 2017.
  * [9] I. Ahmed _et al._ , “A survey on hybrid beamforming techniques in 5G: Architecture and system model perspectives,” _IEEE Commun. Surveys Tuts._ , vol. 20, no. 4, pp. 3060–3097, 4th Quart. 2018.
  * [10] A. Koc _et al._ , “3D angular-based hybrid precoding and user grouping for uniform rectangular arrays in massive MU-MIMO systems,” _IEEE Access_ , vol. 8, pp. 84 689–84 712, May 2020.
  * [11] A. Koc _et al._ , “Hybrid millimeter-wave massive MIMO systems with low CSI overhead and few-bit DACs/ADCs,” in _2020 IEEE 92nd Veh. Technol. Conf. (VTC2020-Fall)_ , Dec. 2020, pp. 1–5.
  * [12] A. Koc _et al._ , “Swarm intelligence based power allocation in hybrid massive MIMO systems,” in _2021 IEEE Wireless Commun. and Netw. Conf. (WCNC)_ , Mar. 2021, pp. 1–7.
  * [13] P. Wang et al., “Joint transceiver and large intelligent surface design for massive MIMO mmWave systems,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 2, pp. 1052–1064, 2021.
  * [14] K. Ying et al., “GMD-based hybrid beamforming for large reconfigurable intelligent surface assisted millimeter-wave massive MIMO,” _IEEE Access_ , vol. 8, pp. 19 530–19 539, 2020.
  * [15] S. H. Hong et al., “Hybrid beamforming for intelligent reflecting surface aided millimeter wave MIMO systems,” _IEEE Trans. Wireless Commun (Early Access)_ , 2022.
  * [16] A. Alkhateeb _et al._ , “MIMO precoding and combining solutions for millimeter-Wave systems,” _IEEE Commun. Mag._ , vol. 52, no. 12, pp. 122–131, Dec. 2014.
  * [17] X. Cheng _et al._ , “Communicating in the real world: 3D MIMO,” _IEEE Wireless Commun._ , vol. 21, no. 4, pp. 136–144, Aug. 2014.
  * [18] S. Park _et al._ , “Dynamic subarrays for hybrid precoding in wideband mmWave MIMO systems,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 5, pp. 2907–2920, Mar. 2017.
  * [19] F. Sohrabi _et al._ , “Hybrid analog and digital beamforming for mmWave OFDM large-scale antenna arrays,” _IEEE J. Sel. Areas Commun._ , vol. 35, no. 7, pp. 1432–1443, July 2017.
  * [20] W. Tong _et al._ , “Artificial intelligence for vehicle-to-everything: A survey,” _IEEE Access_ , vol. 7, pp. 10 823–10 843, 2019.
  * [21] A. I. Sulyman et al., “Directional radio propagation path loss models for millimeter-wave wireless networks in the 28-, 60-, and 73-GHz bands,” _IEEE Trans. Wireless Commun._ , vol. 15, no. 10, pp. 6939–6947, 2016.


# Uncertainty Quantification of Deep Learning for Spatiotemporal Data:
Challenges and Opportunities

Wenchong He<EMAIL_ADDRESS>University of FloridaGainesvilleFLUSA and  Zhe Jiang
<EMAIL_ADDRESS>University of FloridaGainesvilleFLUSA

###### Abstract.

With the advancement of GPS, remote sensing, and computational simulations,
large amounts of geospatial and spatiotemporal data are being collected at an
increasing speed. Such emerging spatiotemporal big data assets, together with
the recent progress of deep learning technologies, provide unique
opportunities to transform society. However, it is widely recognized that deep
learning sometimes makes unexpected and incorrect predictions with unwarranted
confidence, causing severe consequences in high-stake decision-making
applications (e.g., disaster management, medical diagnosis, autonomous
driving). Uncertainty quantification (UQ) aims to estimate a deep learning
model’s confidence. This paper provides a brief overview of UQ of deep
learning for spatiotemporal data, including its unique challenges and existing
methods. We particularly focus on the importance of uncertainty sources. We
identify several future research directions for spatiotemporal data.

††conference: 2nd KDD Workshop on Uncertainty Reasoning and Quantification in
Decision Making; August 07, 2023; Long Beach, CA

##  1\. Introduction

With the advancement of GPS, remote sensing, and computational simulations,
large amounts of geospatial and spatiotemporal data are being collected at an
increasing speed (Shekhar et al., 2015). Such emerging spatiotemporal big data
assets, together with the recent progress of deep learning technologies,
provide unique opportunities to transform society in broad applications. For
example, deep learning is widely used to process spatiotemporal data from
radar or lidar sensors and video cameras to monitor road conditions, detect
pedestrians, and navigate through traffic (Zhao et al., 2019). In disaster
management, deep learning systems have been developed to analyze satellite or
drone imagery to enhance situational awareness during deadly hurricane flood
disasters (Jiang et al., 2019). Although deep learning is known for higher
prediction accuracy compared with many traditional machine learning
techniques, it is widely recognized that they can sometimes make unexpected
and incorrect predictions with unwarranted confidence, particularly in complex
real-world environments (Reichstein et al., 2019). This can have serious
consequences in high-stakes applications like autonomous driving (Choi et al.,
2019), medical diagnosis (Begoli et al., 2019), and disaster response (Alam et
al., 2017). Therefore, uncertainty quantification is essential for a deep
learning model to be aware of its limitations and avoid overconfident
predictions.

Applications: One important application is disaster response. Deep learning
has been used to predict the track of hurricanes or estimate flooded areas,
the results of which directly impact decision-makers in planning evacuation
and rescue efforts. Thus, in hurricane tracking, scientists often provide not
only the most likely point of landfall but also provide a “cone of
uncertainty” across other likely points of impact and future trajectories of
the storm. Similarly, in autonomous driving, complex environmental factors
like extreme weather or the ambiguous appearance of nearby vehicles can fool a
deep learning model to ignore an obstacle and cause traffic crashes. In the
medical domain, deep learning has been widely used for medical image analysis,
clinical diagnosis, and treatment planning. Overconfident predictions can not
only cause unnecessary medical expenses but also endanger patient life (Loftus
et al., 2022).

Challenges: Uncertainty quantification of deep learning for spatiotemporal
data poses unique challenges due to their special data characteristics. First,
spatiotemporal data violate the common assumption that samples follow an
identical and independent distribution. Instead, implicit dependency
structures exist in continuous space and time (e.g., spatial and temporal
autocorrelation, and temporal dynamics) (Shekhar et al., 2015; Jiang, 2018; He
et al., 2022c). Thus, the uncertainty quantification process should be aware
of such a dependency structure. Second, spatiotemporal data have various
spatial, temporal, and spectral resolutions and diverse sources of noise and
errors (e.g., sensor noise, obstacles, and atmospheric effects in remote
sensing signals (Licata and Mehta, 2022), GPS errors). Analyzing such data
often requires the co-registration of different layers (e.g., points, lines,
polygons, geo-rasters) into the same spatial reference system. The process is
subject to registration uncertainty due to GPS errors or annotation mistakes
in map generation (Jiang et al., 2022; He et al., 2022a). Third,
spatiotemporal data are heterogeneous (non-stationary), i.e., the data
distribution often varies across different regions or time periods (Jiang et
al., 2019). Thus, a deep learning model trained in one region (or time) may
not generalize well to another region (or time). Spatiotemporal non-stationary
requires characterizing uncertainty due to out-of-distribution samples
(Shekhar et al., 2015). This issue is particularly important when spatial
observation samples are sparsely distributed, causing uncertainty when
inferring the observations at other locations in continuous space (Hengl et
al., 2017). Moreover, modeling such uncertainty needs to consider sample
density both in the non-spatial feature space and in the geographic space.
Fourth, the large volume of spatiotemporal data requires efficient computation
for prediction and uncertainty quantification, especially considering the need
of capturing long-range and complex dependencies. Finally, spatiotemporal
phenomena are often governed by physical laws (e.g., water flow, temperature
diffusion), and it is crucial to consider the inherent physics knowledge of
the system when quantifying prediction uncertainty for spatiotemporal systems.

This paper does not aim to provide a thorough review of uncertainty
quantification methods for deep learning, as this can be found in several
existing surveys (Gawlikowski et al., 2021; He and Jiang, 2023). Instead, the
paper provides a brief overview of uncertainty quantification methods for
spatiotemporal data. We first categorize the sources of uncertainty into two
types: data uncertainty and model uncertainty. We delve into these sources and
explore how they can be effectively represented in neural network models. We
briefly summarize existing methodologies for quantifying uncertainty in deep
neural networks, taking into account the perspective of uncertainty sources.
We analyze the advantages and disadvantages of these methodologies, shedding
light on their applicability and limitations in the context of spatiotemporal
data. Finally, we identify several future research directions related to
uncertainty quantification for spatiotemporal data.

##  2\. Types of uncertainty source

In this section categorize the sources of uncertainty in GeoAI into two types:
data uncertainty and model uncertainty. We discuss the potential sources of
each type and its representation.

###  2.1. Data uncertainty

####  2.1.1. Source of data uncertainty

(a) Classes distribution with clean separable boundary (b) Entropy of samples
in (a) (c) Classes distribution with ambiguous boundary (d) Entropy of samples
in (c)

Figure 1. Data uncertainty visualization examples (Different colors represent
samples in different classes)

Data uncertainty, also known as aleatoric uncertainty, arises from the
inherent randomness, noise, or overlapping feature distribution within the
data, which cannot be eliminated even with additional training data (Yarin,
2016). Randomness or noise in data can stem from various factors during data
acquisition. Uncertainty in the data acquisition process can be attributed to
instrument errors, inadequate data sampling frequency, and transmission errors
(Hariri et al., 2019). Complex environmental conditions, such as poor weather,
can hinder the proper functioning of data acquisition devices. Additionally,
uncertainty can be influenced by inappropriate sampling methods, data storage,
data representation techniques, and interpolation methods during data
processing (Yarin, 2016).

In the case of spatiotemporal data collected from space and airborne platforms
like CubeSats and UAVs, data uncertainty can arise not only from sensor errors
in data acquisition devices but also from the discretization of data in
digital formats despite the underlying phenomenon being continuous (Cheng et
al., 2014). The uncertainty in representing an object’s movement is also
impacted by the sampling frequency, or sample rate, at which location samples
are collected (Pfoser and Jensen, 1999). Furthermore, uncertainty within the
data can accumulate from multiple sources and propagate into the model.

####  2.1.2. Data uncertainty representation

In a classification problem, the data uncertainty arises from the complexity
of the data and the structure of the class boundary in the feature space. The
uncertainty of the class variable, given a specific input instance, is
quantified by the entropy of the true class distribution
$\mathcal{H}[p(y|\boldsymbol{x})]$ where $y$ is the class label and
$\boldsymbol{x}$ is the sample feature. The entropy reflects the randomness of
the class distribution due to feature overlap among samples from different
classes as shown in Fig. 1.

For regression problems, data uncertainty arises from the inherent noise or
variability in the data generation process. The observation noise, denoted as
$\epsilon(\boldsymbol{x})$, is added to the true function $f(\boldsymbol{x})$
to obtain the target variable $y$. There are two types of noise: homeostatic
noise, which assumes a constant observation noise across all inputs, and
heteroscedastic noise, where the observation noise varies as a function of the
input $\epsilon(\boldsymbol{x})\sim p(\epsilon|\boldsymbol{x})$. The
heteroscedastic noise model is useful when the noise level differs among
samples.

###  2.2. Model uncertainty

####  2.2.1. Source of model uncertainty

Model uncertainty, also known as aleatoric uncertainty, encompasses the
uncertainty in a model’s predictions resulting from imperfections in the model
training process. In the context of spatiotemporal data, model uncertainty can
be attributed to three primary sources: uncertainty in model architecture,
uncertainty in model parameters, and uncertainty due to dataset distribution
mismatch. Uncertainty in model architecture arises from the lack of
understanding regarding the most suitable model architecture for a given
geospatial dataset. For instance, in deep learning models, uncertainty may
exist regarding the optimal number of neural network layers and neurons in
each layer, as an overly complex model can lead to overfitting. Uncertainty in
model parameters stems from unknown optimal parameter values. This uncertainty
can arise due to factors such as an improper training strategy, limited
geospatial training instances, or convergence to a local optimum rather than
the global minimum of the loss function. Therefore, the model’s weight values
may not accurately represent the true optimal solution. The last type of model
uncertainty is caused by dataset distribution drift, where the distribution of
test samples differs from that of the training dataset. This issue, known as
out-of-distribution (OOD) data, is not uncommon in real-world spatiotemporal
deployments, as the test cases often involve complex and diverse scenarios.

####  2.2.2. Model uncertainty representation

Representing model uncertainty in spatiotemporal data poses challenges due to
its multiple sources. Different methods can be adopted to estimate and
represent uncertainty associated with each type. For uncertainty stemming from
model parameters, Bayesian neural networks (BNNs) are commonly used (Jospin et
al., 2022). BNNs assume a prior distribution over the model parameters and aim
to infer the posterior distribution to reflect parameter uncertainty.
Uncertainty arising from model architectures can be estimated using deep
ensembles. This approach involves constructing an ensemble of neural network
architectures, training each model separately, and generating predictions that
form a distribution on the target variable. The variance of these predictions
serves as an estimation of prediction uncertainty. Uncertainty resulting from
dataset distribution mismatch can be assessed by considering the proximity of
new test samples to the training samples. As the test sample deviates further
from the distribution of training data, the model uncertainty increases.

In summary, model uncertainty in spatiotemporal data arises from
misspecifications in model architectures, parameters, and dataset
distributions. Depending on the specific application, one type of uncertainty
may dominate over the other, necessitating tailored methods to address it. The
recognition and appropriate handling of these uncertainties are crucial for
robust analysis and interpretation of data and model outputs in various
geospatial fields.

##  3\. Methodolies

In the following, the main intuitions and approaches of the three types are
presented and their main advantages and disadvantages are discussed.

###  3.1. Model uncertainty

####  3.1.1. Bayesian Neural Networks:

The Bayesian neural network (BNN) incorporates a prior distribution
$p(\boldsymbol{\theta})$ on the neural network parameters and learns the
posterior distribution $p(\boldsymbol{\theta}|\mathbf{X},\mathbf{Y})$ based on
the training dataset in Eq. 1. However, analytically solving for this
posterior distribution is intractable, and approximation methods are necessary
for prediction in BNN.

(1) $\footnotesize
p(\boldsymbol{\theta}|\mathbf{X},\mathbf{Y})=\frac{p(\mathbf{Y}|\mathbf{X},\boldsymbol{\theta})p(\boldsymbol{\theta})}{p(\mathbf{Y}|\mathbf{X})}$

One approach is to select an approximation $q_{\phi}(\boldsymbol{\theta})$
from a parameterized class of distributions $\mathcal{Q}$ to approximate the
posterior. Popular optimization methods for selecting
$q_{\phi}(\boldsymbol{\theta})$ include variational inference (Blei et al.,
2017) and Laplace approximation (Friston et al., 2007), both of which impose
assumptions and restrictions on the form of the approximated posterior.
However, these restrictions can lead to inaccuracies in predictions and
uncertainty quantification. Another widely used approach is the MC dropout
method (Gal and Ghahramani, 2016), which is simple and easy to implement. It
demonstrates that optimizing a neural network with a dropout layer is
equivalent to approximating a BNN using variational inference on a parametric
Bernoulli distribution. Uncertainty estimation is obtained by computing the
variance across multiple stochastic forward predictions with different dropout
masks. However, MC dropout tends to be less calibrated than other baseline
uncertainty quantification methods in many benchmark datasets.

####  3.1.2. Ensemble models

Ensemble models are a powerful approach that involves combining multiple
neural network models during the prediction process. By aggregating the
predictions of individual models, an output distribution is formed. The
variability in predictions among the ensemble models can serve as an indicator
of model uncertainty, where a higher variance implies a greater degree of
uncertainty. To capture this uncertainty stemming from various factors,
several strategies for constructing ensembles can be employed, such as
bootstrapping or combining different neural network architectures
(Lakshminarayanan et al., 2017).

####  3.1.3. Sample-density aware neural networks

The mentioned approaches do not effectively handle model uncertainty arising
from low sample density, where samples lying far from the training set support
may lead to overly confident predictions. To address this, various approaches
have been developed to create sample density-aware neural networks capable of
capturing model uncertainty in such scenarios. These approaches include
Gaussian process modeling (e.g., kriging (Montero et al., 2015)) for spatial
data, deep Gaussian processes, and distance-aware neural networks. Distance-
aware neural networks, inspired by Gaussian process models, aim to
characterize uncertainty based on the density of sample features. They utilize
the feature extraction capabilities of deep neural networks (DNNs) to learn a
hidden representation $h(\boldsymbol{x})$ that reflects meaningful distances
within the data manifold (Van Amersfoort et al., 2020).

###  3.2. Data uncertainty

Generally speaking, data uncertainty is represented by the entropy of the
distribution $p(y|\boldsymbol{x},\boldsymbol{\theta})$, where
$\boldsymbol{\theta}$ is the neural network parameters. Existing approaches
for learning this distribution can be classified into deep discriminative
models and generative models.

####  3.2.1. Deep discriminative models

Deep discriminative models can be further categorized as parametric or non-
parametric based on the format of the distribution. To quantify data
uncertainty, a discriminative model directly outputs a predictive distribution
using a neural network. The distribution can be represented by a parametric
model, which assumes a parameterized family of probability distributions
(e.g., Gaussian or mixture Gaussian) whose parameters (e.g., mean and
variance) are predicted by the neural network (Kendall and Gal, 2017).
Alternatively, a non-parametric model does not make any assumptions about the
underlying distributions and outputs a prediction interval (Pearce et al.,
2018). Prediction intervals provide a lower and upper bound $[y_{l},y_{u}]$,
within which the ground truth $y$ is expected to fall with a prescribed
confidence level of $1-\alpha$ (i.e., $p(y\in[y_{l},y_{u}])>1-\alpha$).
However, a drawback is that standard optimization strategies may not be
applicable.

####  3.2.2. Deep generative models

Deep generative models (DGMs) are capable of learning the intractable data
distribution $p_{\text{data}}(\boldsymbol{x})$ in the high-dimensional feature
space $\mathcal{X}\in\mathbb{R}^{n}$ from a large number of independent and
identically distributed observed samples. To quantify DNN data uncertainty,
the basic idea is to employ DGMs to learn the predictive distribution
$p(y|\boldsymbol{x})$ using the conditional deep generative model (cDGM) (Sohn
et al., 2015). Uncertainty quantification models based on cDGMs aim to learn a
conditional density over the prediction $y$, given the input feature
$\boldsymbol{x}$. This involves learning a model
$g_{\boldsymbol{\theta}}(\boldsymbol{z},\cdot):\mathbb{X}\rightarrow\mathbb{Y}$,
where the generative model $g(\boldsymbol{z},\boldsymbol{x})$ with
$\boldsymbol{z}\sim p(\boldsymbol{z})$ approximates the true unknown
distribution $p_{\text{true}}(y|\boldsymbol{x})$. The variability of the
prediction distribution $p(y|\boldsymbol{x})$ is encoded in the latent
variable $\boldsymbol{z}$ and the generative model. During inference, for any
$\boldsymbol{x}\in\mathbb{X}$, we can generate $m$ samples with
$y_{i}=g_{\boldsymbol{\theta}}(\boldsymbol{z}_{i},\boldsymbol{x})$, where
$\boldsymbol{z}_{i}\sim p(\boldsymbol{z})$. From these samples
$\\{y_{i}\\}_{i=1}^{m}$, we can quantify the prediction uncertainty by
measuring their variability.

###  3.3. The combination of model and data uncertainty

Many frameworks have been proposed to jointly consider both data and model
uncertainty, aiming for more accurate uncertainty quantification. One
straightforward approach is to select one method from each category and
combine them within a single framework. For instance, Bayesian neural network
(BNN) models like Monte-Carlo dropout or ensemble models can be merged with
prediction distributions, where uncertainty can be obtained from the total
variance of the prediction. However, such simple combinations often introduce
significant computational complexity. Alternatively, other methods have
developed evidential deep learning (Malinin and Gales, 2018), which integrates
both data and model uncertainty within a single deterministic model using
evidence theory. This approach offers computational efficiency, but it
requires the design of new optimization strategies and may not be suitable for
certain network architectures.

##  4\. Future Directions

###  4.1. Spatial Sample Density and Nonstationarity

As discussed in the introduction, one major challenge of spatiotemporal data
is spatiotemporal non-stationarity. The data distribution can vary from one
region (time) to another region (time). Such phenomena can be characterized as
out-of-distribution (OOD) data or spatiotemporal outliers. Given a training
data distribution $p(x)$, the OOD data are those samples that are either
unlikely under the training data distribution or outside the support of
$p(x)$. Accurate detection of OOD samples is of paramount importance in
spatiotemporal model generalizability. Because the space is continuous and the
boundary of individual homogeneous sub-regions is implicit, the model needs to
learn such spatial patterns in order to quantify uncertainty due to
nonstationarity. Another relevant source of spaital uncertainty is due to
sparse training samples in the geographical space. This is often due to
limited sensor observations. Traditionally, Gaussian process has been widely
used to quantify such spatial uncertainty in continuous space. However, for
deep neural network models, new techniques are needed that consider sample
density both in the non-spatial feature space and in the geographic space.

###  4.2. Spatial Imaging and inverse problem

The goal of the imaging process is to reconstruct an unknown image from
measurements, which is an inverse problem commonly used in medical imaging
(e.g., magnetic resonance imaging and X-ray computed tomography) and
geophysical inversion (e.g., seismic inversion) (Fessler, 2010). However, this
process is challenging due to the limited and noisy information used to
determine the original image, leading to structured uncertainty and
correlations between nearby pixels in the reconstructed image (Kendall and
Gal, 2017). To overcome this issue, current research in uncertainty
quantification of inverse problems employs conditional deep generative models,
such as cVAE, cGAN, and conditional normalizing flow models (Dorta et al.,
2018). These methods utilize a low-dimensional latent space for image
generation but may overlook unique data characteristics, such as structural
constraints from domain physics in certain types of image data, such as remote
sensing images, MRI images, or geological subsurface images (Sun et al.,
2021). The use of physics-informed models may improve uncertainty
quantification in these cases. It’s promising to incorporate the physics
constraints for quantifying the uncertainty associated with the imaging
process.

###  4.3. UQ for physics-aware DNN models

Many applications in the field of spatiotemporal data minging involve physical
systems that can be described using principles such as partial differential
equations (PDEs) governing diffusion processes. While deep learning has proven
effective in extracting complex patterns from data, standard deep neural
networks lack explicit incorporation of the underlying physics (Kashinath et
al., 2021; He et al., 2022b). To address this challenge, physics-informed
neural networks (PINNs) have been developed to integrate physical principles
and domain knowledge into deep learning for consistent predictions
(Krishnapriyan et al., 2021). PINNs encode the governing equations as residual
losses to guide the optimization of neural networks (Karniadakis et al.,
2021), enabling the incorporation of physical constraints during model
training.

However, many physical systems exhibit non-deterministic or unknown underlying
principles. Chaotic and stochastic behavior is inherent in these systems, as
exemplified by weather forecasting and climate prediction, where small
perturbations can lead to significant differences in the predicted state
(known as the ”butterfly effect”) (Kashinath et al., 2021). Even with
deterministic initial conditions, the governing equations (e.g., PDEs) driving
the system may be stochastic, described by stochastic differential equations.
Uncertainty quantification is crucial for improving the reliability of
predictions, especially under distribution shifts. The uncertainty in modeling
physical systems can arise from multiple sources. First, the initial and
boundary conditions may be non-deterministic, and the system itself may
exhibit chaotic behavior (Wang et al., 2019). Second, the underlying physical
principles may not be fully known, or the parameters of the governing
equations may be stochastic, leading to violations of conservation laws in
imperfect systems (Yang and Perdikaris, 2019).

These cases highlight the presence of inherent data uncertainty in stochastic
physical systems. Probabilistic models offer a natural way to model
distributions and incorporate stochasticity and uncertainty into neural
networks. However, quantifying uncertainty in PINNs poses specific challenges.
It requires simultaneously considering the physical principles and their
uncertainties. Incorporating physical constraints can help mitigate data and
model uncertainty (He and Jiang, 2020). Various sources of uncertainty may
arise, including randomness in the physical system itself, measurement errors,
and limited knowledge of the governing equations. One potential approach to
uncertainty quantification in PINNs involves building probabilistic neural
networks that propagate uncertainty from multiple sources based on deep
generative models for structured outputs (Oberdiek et al., 2022).

##  5\. Conclusion

In this paper, we provides a brief overview of UQ of deep learning for
spatiotemporal data, including its unique challenges and existing methods. We
particularly focus on the importance of uncertainty sources. We also identify
several future research directions related to spatiotemporal data.

## References

  * (1)
  * Alam et al. (2017) Firoj Alam, Muhammad Imran, and Ferda Ofli. 2017.  Image4act: Online social media image processing for disaster response. In _Proceedings of the 2017 IEEE/ACM international conference on advances in social networks analysis and mining 2017_. 601–604.
  * Begoli et al. (2019) Edmon Begoli, Tanmoy Bhattacharya, and Dimitri Kusnezov. 2019\.  The need for uncertainty quantification in machine-assisted medical decision making.  _Nature Machine Intelligence_ 1, 1 (2019), 20–23.
  * Blei et al. (2017) David M Blei, Alp Kucukelbir, and Jon D McAuliffe. 2017\.  Variational inference: A review for statisticians.  _Journal of the American statistical Association_ 112, 518 (2017), 859–877.
  * Cheng et al. (2014) Reynold Cheng, Tobias Emrich, Hans-Peter Kriegel, Nikos Mamoulis, Matthias Renz, Goce Trajcevski, and Andreas Züfle. 2014\.  Managing uncertainty in spatial and spatio-temporal data. In _2014 IEEE 30th International Conference on Data Engineering_. IEEE, 1302–1305.
  * Choi et al. (2019) Jiwoong Choi, Dayoung Chun, Hyun Kim, and Hyuk-Jae Lee. 2019\.  Gaussian yolov3: An accurate and fast object detector using localization uncertainty for autonomous driving. In _Proceedings of the IEEE/CVF International Conference on Computer Vision_. 502–511.
  * Dorta et al. (2018) Garoe Dorta, Sara Vicente, Lourdes Agapito, Neill DF Campbell, and Ivor Simpson. 2018\.  Structured uncertainty prediction networks. In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_. 5477–5485.
  * Fessler (2010) Jeffrey A Fessler. 2010\.  Model-based image reconstruction for MRI.  _IEEE signal processing magazine_ 27, 4 (2010), 81–89.
  * Friston et al. (2007) Karl Friston, Jérémie Mattout, Nelson Trujillo-Barreto, John Ashburner, and Will Penny. 2007.  Variational free energy and the Laplace approximation.  _Neuroimage_ 34, 1 (2007), 220–234.
  * Gal and Ghahramani (2016) Yarin Gal and Zoubin Ghahramani. 2016.  Dropout as a bayesian approximation: Representing model uncertainty in deep learning. In _international conference on machine learning_. PMLR, 1050–1059.
  * Gawlikowski et al. (2021) Jakob Gawlikowski, Cedrique Rovile Njieutcheu Tassi, Mohsin Ali, Jongseok Lee, Matthias Humt, Jianxiang Feng, Anna Kruspe, Rudolph Triebel, Peter Jung, Ribana Roscher, et al. 2021\.  A survey of uncertainty in deep neural networks.  _arXiv preprint arXiv:2107.03342_ (2021).
  * Hariri et al. (2019) Reihaneh H Hariri, Erik M Fredericks, and Kate M Bowers. 2019\.  Uncertainty in big data analytics: survey, opportunities, and challenges.  _Journal of Big Data_ 6, 1 (2019), 1–16.
  * He and Jiang (2020) Wenchong He and Zhe Jiang. 2020.  Semi-supervised learning with the em algorithm: A comparative study between unstructured and structured prediction.  _IEEE Transactions on Knowledge and Data Engineering_ 34, 6 (2020), 2912–2920.
  * He and Jiang (2023) Wenchong He and Zhe Jiang. 2023.  A Survey on Uncertainty Quantification Methods for Deep Neural Networks: An Uncertainty Source Perspective.  _arXiv preprint arXiv:2302.13425_ (2023).
  * He et al. (2022a) Wenchong He, Zhe Jiang, Marcus Kriby, Yiqun Xie, Xiaowei Jia, Da Yan, and Yang Zhou. 2022a.  Quantifying and Reducing Registration Uncertainty of Spatial Vector Labels on Earth Imagery. In _Proceedings of the 28th ACM SIGKDD Conference on Knowledge Discovery and Data Mining_. 554–564.
  * He et al. (2022b) Wenchong He, Arpan Man Sainju, Zhe Jiang, Da Yan, and Yang Zhou. 2022b.  Earth Imagery Segmentation on Terrain Surface with Limited Training Labels: A Semi-supervised Approach based on Physics-Guided Graph Co-Training.  _ACM Transactions on Intelligent Systems and Technology (TIST)_ 13, 2 (2022), 1–22.
  * He et al. (2022c) Wenchong He, Minh N Vu, Zhe Jiang, and My T Thai. 2022c.  An explainer for temporal graph neural networks. In _GLOBECOM 2022-2022 IEEE Global Communications Conference_. IEEE, 6384–6389.
  * Hengl et al. (2017) Tomislav Hengl, Jorge Mendes de Jesus, Gerard BM Heuvelink, Maria Ruiperez Gonzalez, Milan Kilibarda, Aleksandar Blagotić, Wei Shangguan, Marvin N Wright, Xiaoyuan Geng, Bernhard Bauer-Marschallinger, et al. 2017\.  SoilGrids250m: Global gridded soil information based on machine learning.  _PLoS one_ 12, 2 (2017), e0169748.
  * Jiang (2018) Zhe Jiang. 2018\.  A survey on spatial prediction methods.  _IEEE Transactions on Knowledge and Data Engineering_ 31, 9 (2018), 1645–1664.
  * Jiang et al. (2022) Zhe Jiang, Wenchong He, Marcus Stephen Kirby, Arpan Man Sainju, Shaowen Wang, Lawrence V Stanislawski, Ethan J Shavers, and E Lynn Usery. 2022\.  Weakly Supervised Spatial Deep Learning for Earth Image Segmentation Based on Imperfect Polyline Labels.  _ACM Transactions on Intelligent Systems and Technology (TIST)_ 13, 2 (2022), 1–20.
  * Jiang et al. (2019) Zhe Jiang, Arpan Man Sainju, Yan Li, Shashi Shekhar, and Joseph Knight. 2019.  Spatial ensemble learning for heterogeneous geographic data with class ambiguity.  _ACM Transactions on Intelligent Systems and Technology (TIST)_ 10, 4 (2019), 1–25.
  * Jospin et al. (2022) Laurent Valentin Jospin, Hamid Laga, Farid Boussaid, Wray Buntine, and Mohammed Bennamoun. 2022\.  Hands-on Bayesian neural networks—A tutorial for deep learning users.  _IEEE Computational Intelligence Magazine_ 17, 2 (2022), 29–48.
  * Karniadakis et al. (2021) George Em Karniadakis, Ioannis G Kevrekidis, Lu Lu, Paris Perdikaris, Sifan Wang, and Liu Yang. 2021.  Physics-informed machine learning.  _Nature Reviews Physics_ 3, 6 (2021), 422–440.
  * Kashinath et al. (2021) K Kashinath, M Mustafa, A Albert, JL Wu, C Jiang, S Esmaeilzadeh, K Azizzadenesheli, R Wang, A Chattopadhyay, A Singh, et al. 2021\.  Physics-informed machine learning: case studies for weather and climate modelling.  _Philosophical Transactions of the Royal Society A_ 379, 2194 (2021), 20200093.
  * Kendall and Gal (2017) Alex Kendall and Yarin Gal. 2017.  What uncertainties do we need in bayesian deep learning for computer vision?  _Advances in neural information processing systems_ 30 (2017).
  * Krishnapriyan et al. (2021) Aditi Krishnapriyan, Amir Gholami, Shandian Zhe, Robert Kirby, and Michael W Mahoney. 2021.  Characterizing possible failure modes in physics-informed neural networks.  _Advances in Neural Information Processing Systems_ 34 (2021), 26548–26560.
  * Lakshminarayanan et al. (2017) Balaji Lakshminarayanan, Alexander Pritzel, and Charles Blundell. 2017.  Simple and scalable predictive uncertainty estimation using deep ensembles.  _Advances in neural information processing systems_ 30 (2017).
  * Licata and Mehta (2022) Richard J Licata and Piyush M Mehta. 2022.  Uncertainty quantification techniques for data-driven space weather modeling: thermospheric density application.  _Scientific Reports_ 12, 1 (2022), 7256.
  * Loftus et al. (2022) Tyler J Loftus, Benjamin Shickel, Matthew M Ruppert, Jeremy A Balch, Tezcan Ozrazgat-Baslanti, Patrick J Tighe, Philip A Efron, William R Hogan, Parisa Rashidi, Gilbert R Upchurch Jr, et al. 2022\.  Uncertainty-aware deep learning in healthcare: a scoping review.  _PLOS digital health_ 1, 8 (2022), e0000085.
  * Malinin and Gales (2018) Andrey Malinin and Mark Gales. 2018.  Predictive uncertainty estimation via prior networks.  _Advances in neural information processing systems_ 31 (2018).
  * Montero et al. (2015) José-María Montero, Gema Fernández-Avilés, and Jorge Mateu. 2015.  _Spatial and spatio-temporal geostatistical modeling and kriging_.  John Wiley & Sons.
  * Oberdiek et al. (2022) Philipp Oberdiek, Gernot A Fink, and Matthias Rottmann. 2022\.  UQGAN: A Unified Model for Uncertainty Quantification of Deep Classifiers trained via Conditional GANs.  _arXiv preprint arXiv:2201.13279_ (2022).
  * Pearce et al. (2018) Tim Pearce, Alexandra Brintrup, Mohamed Zaki, and Andy Neely. 2018.  High-quality prediction intervals for deep learning: A distribution-free, ensembled approach. In _International conference on machine learning_. PMLR, 4075–4084.
  * Pfoser and Jensen (1999) Dieter Pfoser and Christian S Jensen. 1999.  Capturing the uncertainty of moving-object representations. In _International Symposium on Spatial Databases_. Springer, 111–131.
  * Reichstein et al. (2019) Markus Reichstein, Gustau Camps-Valls, Bjorn Stevens, Martin Jung, Joachim Denzler, and Nuno Carvalhais. 2019\.  Deep learning and process understanding for data-driven Earth system science.  _Nature_ 566, 7743 (2019), 195–204.
  * Shekhar et al. (2015) Shashi Shekhar, Zhe Jiang, Reem Y Ali, Emre Eftelioglu, Xun Tang, Venkata MV Gunturi, and Xun Zhou. 2015.  Spatiotemporal data mining: A computational perspective.  _ISPRS International Journal of Geo-Information_ 4, 4 (2015), 2306–2338.
  * Sohn et al. (2015) Kihyuk Sohn, Honglak Lee, and Xinchen Yan. 2015.  Learning structured output representation using deep conditional generative models.  _Advances in neural information processing systems_ 28 (2015).
  * Sun et al. (2021) Jian Sun, Kristopher A Innanen, and Chao Huang. 2021\.  Physics-guided deep learning for seismic inversion with hybrid training and uncertainty analysis.  _Geophysics_ 86, 3 (2021), R303–R317.
  * Van Amersfoort et al. (2020) Joost Van Amersfoort, Lewis Smith, Yee Whye Teh, and Yarin Gal. 2020.  Uncertainty estimation using a single deep deterministic neural network. In _International conference on machine learning_. PMLR, 9690–9700.
  * Wang et al. (2019) Bin Wang, Jie Lu, Zheng Yan, Huaishao Luo, Tianrui Li, Yu Zheng, and Guangquan Zhang. 2019.  Deep uncertainty quantification: A machine learning approach for weather forecasting. In _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_. 2087–2095.
  * Yang and Perdikaris (2019) Yibo Yang and Paris Perdikaris. 2019.  Adversarial uncertainty quantification in physics-informed neural networks.  _J. Comput. Phys._ 394 (2019), 136–152.
  * Yarin (2016) Gal Yarin. 2016\.  Uncertainty in deep learning.  _University of Cambridge, Cambridge_ (2016).
  * Zhao et al. (2019) Ling Zhao, Yujiao Song, Chao Zhang, Yu Liu, Pu Wang, Tao Lin, Min Deng, and Haifeng Li. 2019.  T-gcn: A temporal graph convolutional network for traffic prediction.  _IEEE Transactions on Intelligent Transportation Systems_ 21, 9 (2019), 3848–3858.