# Robust Design for Intelligent Reflecting Surface Assisted MIMO-OFDMA
Terahertz Communications

Wanming Hao, Gangcan Sun, Ming Zeng, Zhengyu Zhu, Zheng Chu, Octavia A. Dobre,
Pei Xiao, W. Hao, G. Sun, and Z. Zhu are with the School of Information
Engineering, Zhengzhou University, Zhengzhou 450001, China. (E-mail: {iewmhao,
iegcsun, iezyzhu}@zzu.edu.cn)M. Zeng and O. A. Dobre are with the Faculty of
Engineering and Applied Science, Memorial University, St. Johns, NL A1B 3X9,
Canada. (E-mail: {mzeng, odobre}@mun.ca)Z. Chu and P. Xiao are with the 5G
Innovation Center, Institute of Communication Systems, University of Surrey,
Guildford GU2 7XH, U.K. (Email: {zheng.chu<EMAIL_ADDRESS>
###### Abstract

Recently, terahertz (THz) communication has drawn considerable attention as
one of the promising technologies for the future wireless communications
owning to its ultra-wide bandwidth. Nonetheless, one major obstacle that
prevents the actual deployment of THz lies in its inherent huge attenuation.
Intelligent reflecting surface (IRS) and multiple-input multiple-output (MIMO)
represent two effective solutions for compensating the large pathloss in THz
systems. In this paper, we consider an IRS-aided multi-user THz MIMO system
with orthogonal frequency division multiple access, where the sparse radio
frequency chain antenna structure is adopted for reducing the power
consumption. The objective is to maximize the weighted sum rate via jointly
optimizing the hybrid analog/digital beamforming at the base station and
reflection matrix at the IRS. Since the analog beamforming and reflection
matrix need to cater all users and subcarriers, it is difficult to directly
solve the formulated problem, and thus, an alternatively iterative
optimization algorithm is proposed. Specifically, the analog beamforming is
designed by solving a MIMO capacity maximization problem, while the digital
beamforming and reflection matrix optimization are both tackled using
semidefinite relaxation technique. Considering that obtaining perfect channel
state information (CSI) is a challenging task in IRS-based systems, we further
explore the case with the imperfect CSI for the channels from the IRS to
users. Under this setup, we propose a robust beamforming and reflection matrix
design scheme for the originally formulated non-convex optimization problem.
Finally, simulation results are presented to demonstrate the effectiveness of
the proposed algorithms.

###### Index Terms:

Hybrid beamforming, Intelligent Reflecting Surfaces, THz, Multiple-input
multiple-output.

##  I Introduction

With the rapid proliferation of various novel applications, such as virtual
reality, augmented reality, and telemedicine, the data rate demands in future
wireless communication are expected to grow explosively [1]. As such, the
sub-6 Gigahertz (GHz) and millimeter-wave (mmWave) may not be able to support
these bandwidth-hungry applications. That being said, terahertz (THz)
communication (0.1-10 THz) has been regarded as a promising technology to deal
with the above problem due to its ultra-wide bandwidth [2, 3]. However, there
are two major shortcomings for THz communications, namely severe signal
attenuation and poor diffraction[4].

Multiple-input multiple-output (MIMO) has been recognized as an effective
technology to enhance the THz signal strength owing to the high beamforming
gain. Indeed, it has been shown that the signal strength grows linearly with
the number of antennas at the base station (BS) [5]. Meanwhile, the small
wavelength in THz makes it easy to pack more antennas together, and form a
massive MIMO array. This way, the problem of severe signal antenuation of THz
can be substantially relieved. Nonetheless, the property of poor diffraction
still makes THz vulnerable to blocking obstacles that break the line-of-sight
(LoS) links. To address this problem, intelligent reflect surface (IRS) can be
deployed to create additional links [6, 7], and thus, enhance the THz systems.
Being equipped with a large number of reconfigurable passive elements [8, 9],
the IRS can reflect the incident signals to any direction via adjusting the
phase shifts. As a result, when there is no direct link between the
transmitter and receiver, communication can still be realized via building a
reflective link with the help of the IRS as shown in Fig. 1. Therefore,
incorporating MIMO and IRS into the THz communication can effectively enhance
the signal reception and reduce the probability of signal blockage.

In this paper, we study a multi-user IRS-aided THz MIMO system, where the BS
employs sparse RF chain structure for lowering the circuit power consumption
[10]. Meanwhile, considering that the wideband THz signals may suffer from
frequency selective fading, orthogonal frequency division multiple (OFDM) is
also adopted. Based on this system model, we design the hybrid analog/digital
beamforming at the BS and the reflection matrix at the IRS for maximizing the
weighted sum rate under perfect and imperfect channel state information (CSI).

###  I-A Related Works

The MIMO THz communication has become a research hotspot in recent years.
Considering the large signal attenuation, Lin et al. study the indoor short
range MIMO THz communications [11, 12]. The authors propose a hybrid
analog/digital beamforming to maximize the energy efficiency of the system.
Busari et al. consider three hybrid beamforming array structures, namely fully
connection, subconnection and overlapped subarray [13]. Then, a single-path
THz channel model is used to investigate the performance of the system under
different array structures. Additionally, due to the ultra-wide bandwidth,
frequency selective hybrid beamforming design in THz system is necessary. For
example, Tan and Dai first analyze the array gain loss in the wideband THz
system and then propose a time delay network to obtain the near-optimal array
gain[14]. However, the complexity of the considered system is prohibitively
high. Yuan et al. build a 3-D wideband THz channel model and propose a two-
stage hybrid analog/digital beamforming for maximizing the capacity of the
system [15]. After that, the imperfect CSI is also considered and a robust
beamforming design scheme is developed.

In parallel, IRS has attracted great attention in the past two years owning to
its ability to enable cost-effective and energy-efficient communications. Wu
and Zhang provide a basic IRS communication system model in [8], based upon
which the joint active beamforming at the BS and passive beamforming at the
IRS is designed to minimize the system power consumption. In addition, Ning et
al. propose to apply THz to IRS [16], and consider the beam training and
hybrid analog/digital beamforming. They propose two effective hierarchical
codebooks and beamforming design schemes to obtain the near-optimal
performance. To study the performance of IRS in frequency-selective fading
channels, Zhang et al. consider a MIMO-OFDM system [17], where only one common
set of IRS reflective matrix is designed for all subcarriers. Based on this, a
new alternative optimization algorithm is proposed. Yang et al. investigate
the channel estimation and beamforming design problem in the IRS-based OFDM
system [18], and propose a practical transmission protocol as well as channel
estimation scheme. On this basis, a strategy of jointly optimizing power
allocation and the reflection matrix is developed for maximizing the
achievable rate.

Although the THz and IRS techniques have been investigated in the literature,
e.g., in [11, 12, 13, 14, 15, 8, 19, 16, 17, 18], most of them do not consider
the hybrid beamforming at the BS for IRS communication [8, 19, 16, 17, 18]. In
fact, in a THz-based IRS communication system, the BS should employ a sparse
RF antenna structure for reducing the power consumption and the multiple
subcarriers transmission technology should be adopted for overcoming the
frequency selection channel fading. In this case, how to design the hybrid
analog/digital beamforming at the BS and reflection matrix at the IRS catering
to all subchannels will be challenging. In addition, how to obtain the perfect
CSI remains a non-trivial task for IRS-based reflection links. For the direct
link from the BS to users, the CSI can be readily estimated by conventional
channel estimation methods. For the indirect link from the BS to the IRS, the
CSI is also relatively easy to obtain since the locations of IRS and BS are
fixed. However, the accurate CSIs of reflection links from the IRS to users
are usually difficult to obtain due to the mobility of users. However, [8, 19,
16, 17, 18] all assume perfect CSI. Although Zhou et al. investigate the
robust beamforming design in an IRS system [20], the conventional multiple
antenna structure and single carrier scenario are considered.

###  I-B Main Contributions

To the best of our knowledge, this is the first work to consider the hybrid
analog/digital beamforming in the IRS-aided THz MIMO-OFDMA under imperfect
CSI, and the main contributions of this paper include:

  * •

We construct an IRS-aided THz MIMO-OFDMA communication system, where the BS
employs sparse RF chain structure for reducing the circuit power consumption.
First, we investigate the joint optimization of the hybrid beamforming at the
BS and reflection matrix at the IRS for maximizing the weighted sum rate under
perfect CSI.

  * •

To solve the formulated non-trivial problem, we first initialize the
reflection matrix. Since all subcarriers share one analog beamforming matrix,
we ignore the multi-user interference and obtain the analog beamforming by
solving the corresponding MIMO capacity optimization problem. We subsequently
reformulate a multi-user weighted sum rate maximization problem to optimize
the digital beamforming. With the help of successive convex approximation
(SCA) and semidefinite relaxation (SDR) techniques, we propose an iterative
algorithm to solve the digital beamforming that mitigates the multi-user
interference. Next, we formulate the reflection matrix optimization problem
under given hybrid analog/digital beamforming, and an iterative algorithm is
proposed to solve it. The above procedure is repeated until convergence.

  * •

Next, we assume that the perfect CSIs of reflection links cannot be obtained,
and there exists bounded estimation error. We apply the same method to solve
the analog beamforming. For the digital beamforming and reflection matrix, we
develop a robust optimization scheme for the weighted sum rate optimization
problem relying on the $\mathcal{S}$-Procedure and the convex approximation
techniques. Finally, our simulation results demonstrate the effectiveness of
the proposed algorithms.

We organize the rest of this paper as follows. The system model and weighted
sum rate optimization problem are introduced in Section II. An alternatively
iterative optimization algorithm is designed in Section III. The imperfect
CSIs from the IRS to users are considered and the corresponding optimization
algorithm is developed in Section IV. Simulation results are presented in
Section VI. Finally, conclusions are drawn in Section VII.

Notations: We use the following notations throughout this paper: $(\cdot)^{T}$
and $(\cdot)^{H}$ denote the transpose and Hermitian transpose, respectively,
$\|\cdot\|$ is the Frobenius norm, ${\mathbb{C}}^{x\times y}$ means the space
of $x\times y$ complex matrix, Re($\cdot$) and Tr($\cdot$) denote real number
operation and trace operation, respectively, and Diag($a_{1},\ldots,a_{n}$) is
a diagonal matrix. $\angle(\cdot)$ represents the phase of a complex number.

Fig. 1: The IRS system model. Fig. 2: The sparse RF chain structure at the BS.

##  II System Model and Problem Formulation

In this section, we first describe the IRS-aided THz MIMO-OFDMA system model
and antenna structure. Next, we present the THz channel transmission model and
corresponding parameters. Finally, we formulate the weighted sum rate
maximization problem.

###  II-A System Model

We consider an IRS-aided THz multi-user MIMO system with OFDMA as shown in
Fig. 1, where the BS is equipped with $N_{\rm{TX}}$ antennas and
$N_{\rm{RF}}\;(N_{\rm{RF}}\leq N_{\rm{TX}})$ RF chains. The diagram of the
sparse RF chain at the BS is illustrated in Fig. 2. We assume that there are
no direct links between BS and users due to the occlusion of walls or other
obstacles, and the users can only receive the reflected signals from IRS. Let
$N_{\rm{IRS}}$, $M$ and $K$ denote the number of IRS elements, users and
subcarriers, respectively.

The received signal on the $k$th subcarrier at the $m$th user can be expressed
as

$\displaystyle
y_{m}[k]={\bf{G}}_{m}[k]{\bf{F}}{\bf{v}}_{m}[k]x_{m}[k]+\sum_{j\neq
m}^{M}{\bf{G}}_{m}[k]{\bf{F}}{\bf{v}}_{j}[k]x_{j}[k]+n_{m}[k],$ (1)

where
${\bf{G}}_{m}[k]=G_{t}G_{r}\eta_{k}{\widehat{\bf{g}}}_{m}[k]{\bf\Phi}{\widehat{\bf{H}}}[k]$,
with $G_{t}$ and $G_{r}$ as the transmit and receive antenna gains,
respectively, and $\eta_{k}$ as the pathloss compensation factor [16].
${\widehat{\bf{g}}}_{m}[k]\in{\mathbb{C}}^{1\times N_{\rm{IRS}}}$ denotes the
channel vector from IRS to the $m$th user on the $k$th subcarrier,
${\bf\Phi}\in{\mathbb{C}}^{N_{\rm{IRS}}\times N_{\rm{IRS}}}$ is the reflection
coefficient matrix with
${\bf\Phi}={\rm{diag}}\\{\phi_{1},...,\phi_{N_{\rm{IRS}}}\\}$,
${\widehat{\bf{H}}}[k]\in{\mathbb{C}}^{N_{\rm{IRS}}\times N_{\rm{TX}}}$
represents the channel matrix from BS to IRS on the $k$th subcarrier,
${\bf{F}}\in{\mathbb{C}}^{N_{\rm{TX}}\times N_{\rm{RF}}}$ is the analog
beamforming matrix with ${\bf{F}}=[{\bf{f}}_{1},...,{\bf{f}}_{N_{\rm{RF}}}]$,
${\bf{v}}_{m}[k]\in{\mathbb{C}}^{N_{\rm{RF}}\times 1}$ and $x_{m}[k]$ denote
the digital beamforming and transmit signal for the $m$th user on the $k$th
subcarrier, respectively, $n_{m}[k]$ is the independent and identically
distributed (i.i.d.) additive white Gaussian noise (AWGN) with zero-mean and
variance $N_{0}$. In (1), the first term is the designed signal, while the
second term is the multi-user interference that must be mitigated by designing
proper digital beamforming and reflection matrix.

Next, we present the THz channel model. Let $f_{c}$ and $B$, respectively,
represent the central frequency and bandwidth. Then, the frequency band of the
$k$th subcarrier can be expressed as
$f_{k}=f_{c}+\frac{B}{K}(k-1-\frac{K-1}{2}),k=1,2,...,K$. Although there are a
few scattering components in THz communication, their power are much lower
(more than 20 dB) than that of LoS component [21], and thus, we only consider
the LoS component and ignore the other scattering components. Accordingly, the
channel matrix ${\widehat{\bf{H}}}[k]$ can be expressed as

$\displaystyle{\widehat{\bf{H}}}[k]=q(f_{k},d){{\bf{H}}}[k],$ (2)

where $q(f_{k},d)$ is the complex path gain satisfying

$\displaystyle q(f_{k},d)=\frac{c}{4\pi fd}e^{-\frac{1}{2}\tau(f_{k})d},$ (3)

where $c$ stands for the speed of light, $\tau(f_{k})$ represents the medium
absorption factor and $d$ is the distance from the BS to IRS [22].
${{\bf{H}}}[k]$ can be expressed as

$\displaystyle{{\bf{H}}}[k]={\bf{a}}_{r}(\theta_{k}){\bf{a}}_{t}^{H}(\varphi_{k}),$
(4)

with ${\bf{a}}_{t}(\theta_{k})$ and ${\bf{a}}_{r}(\varphi_{k})$, respectively,
as the antenna array response vector of the transmitter and receiver, namely

$\displaystyle{\bf{a}}_{t}(\theta_{k})=\frac{1}{\sqrt{N_{\rm{TX}}}}\left[1,e^{j\pi\theta_{k}},e^{j2\pi\theta_{k}},\cdots,e^{j(N_{\rm{TX}}-1)\pi\theta_{k}}\right]^{T},$
(5a)
$\displaystyle{\bf{a}}_{r}(\varphi_{k})=\frac{1}{\sqrt{N_{\rm{IRS}}}}\left[1,e^{j\pi\varphi_{k}},e^{j2\pi\varphi_{k}},\cdots,e^{j(N_{\rm{IRS}}-1)\pi\varphi_{k}}\right]^{T}.$
(5b)

Here, $\theta_{k}=2d_{0}f_{k}\sin(\phi_{t})/c$ and
$\varphi_{k}=2d_{0}f_{k}\sin(\phi_{r})/c$, $d_{0}$ denotes the antenna
distance, and $\phi_{t}/\phi_{r}\in[-\pi/2,\pi/2]$ are, respectively, angle of
departure (AoD) and angle of arrival (AoA). Similarly, ${\bf{g}}_{m}[k]$ can
be expressed as

$\displaystyle{\widehat{\bf{g}}}_{m}[k]=q(f_{k},d_{m}){{\bf{g}}}_{m}[k],$ (6)

where
${{\bf{g}}}_{m}[k]=\frac{1}{\sqrt{N_{\rm{IRS}}}}\left[1,e^{j\pi\varphi_{k,m}},e^{j2\pi\varphi_{k,m}},\cdots,e^{j(N_{\rm{IRS}}-1)\pi\varphi_{k,m}}\right]$,
and $q(f_{k},d_{m})$ is defined as

$\displaystyle q(f_{k},d_{m})=\frac{c}{4\pi
f_{k}d_{m}}e^{-\frac{1}{2}\tau(f_{k})d_{m}},$ (7)

with $d_{m}$ as the distance from the IRS to the $m$th user. The BS-IRS-user
$m$ link channel can such be expressed as

$\displaystyle{\bf{G}}_{m}[k]=u_{m}[k]{{\bf{g}}}_{m}[k]{\bf\Phi}{{\bf{H}}}[k],$
(8)

where $u_{m}[k]=G_{t}G_{r}\eta_{k}q(f,d_{m})q(f_{k},d_{m})$. According to
[22], the cascaded path loss of the BS-IRS-user link should satisfy

$\displaystyle\eta_{k}q(f_{k},d)q(f_{k},d_{m})=\frac{\chi
c}{8\sqrt{\pi^{3}}f_{k}dd_{m}}e^{-\frac{1}{2}\tau(f_{k})(d+d_{m})},$ (9)

where $\chi$ is the IRS element gain.

Finally, we rewrite (1) as

$\displaystyle
y_{m}[k]=u_{m}[k]{{\bf{g}}}_{m}[k]{\bf\Phi}{{\bf{H}}}[k]{\bf{F}}{\bf{v}}_{m}[k]x_{m}[k]+\sum\nolimits_{j\neq
m}^{M}u_{m}[k]{{\bf{g}}}_{m}[k]{\bf\Phi}{{\bf{H}}}[k]{\bf{F}}{\bf{v}}_{j}[k]x_{j}[k]+n_{m}[k].$
(10)

###  II-B Problem Formulation

By employing (10), the achievable rate of the $m$th user on the $k$th
subcarrier can be expressed as

$\displaystyle
R_{m}[k]=\frac{B}{K}\log\left(1\\!+\\!\frac{\left|u_{m}[k]{{\bf{g}}}_{m}[k]{\bf\Phi}{{\bf{H}}}[k]{\bf{F}}{\bf{v}}_{m}[k]\right|^{2}}{\sum\nolimits_{j\neq
m}^{M}\left|u_{m}[k]{{\bf{g}}}_{m}[k]{\bf\Phi}{{\bf{H}}}[k]{\bf{F}}{\bf{v}}_{j}[k]\right|^{2}\\!+\\!BN_{0}/K}\right),$
(11)

and thus, the achievable sum rate for the $m$th user can be written as

$\displaystyle R_{m}=\sum\nolimits_{k=1}^{K}R_{m}[k].$ (12)

Next, we formulate the weighted sum rate maximization problem as follows

$\displaystyle\underset{\left\\{{\bf\Phi},{\bf{F}},{\bf{v}}_{m}[k]\right\\}}{\rm{max}}\;\sum\nolimits_{m=1}^{M}\alpha_{m}R_{m}$
(13a) $\displaystyle{\rm{s.t.}}$
$\displaystyle\;|\phi_{i}|=1,i\in\\{1,\cdots,N_{\rm{IRS}}\\},$ (13b)
$\displaystyle\;\sum\nolimits_{M=1}^{M}\sum\nolimits_{k=1}^{K}||{\bf{F}}{\bf{v}}_{m}[k]||^{2}\leq
P_{\rm{max}},$ (13c)
$\displaystyle\;{\bf{F}}(i,j)=1/\sqrt{N_{\rm{TX}}},i\in\\{1,\cdots,N_{\rm{TX}}\\},j\in\\{1,\cdots,N_{\rm{RF}}\\},$
(13d)

where $\alpha_{m}$ denotes the weight of the $m$th user, (13b) is the uni-
modular constraint for each reflection coefficient $\phi_{i}$, (13c) is the
sum transmit power constraint, and (13d) is the amplitude constraint of analog
beamforming. The objective of (13) is to jointly optimize the reflection
matrix ${\bf\Phi}$ and hybrid analog/digital beamforming ${\bf{F}}$ and
${\bf{v}}_{m}[k]$ for maximizing the weighed sum rate of the system.
Obviously, (13) is a non-convex optimization problem due to the non-convex
objective function (13a) and constrains (13b), (13d). Finding the optimal
solution is a challenging task, and we propose an effective alternatively
iterative optimization algorithm to deal with it.

##  III Solution of The Weighted Sum Rate Optimization Problem

In this section, we propose an alternatively iterative optimization algorithm.
First, we consider the hybrid analog/digital beamforming of ${\bf{F}}$,
${\bf{v}}_{m}[k]$ under given reflection matrix ${\bf\Phi}$. We formulate a
MIMO-OFDM sum rate maximization problem and use it to obtain the analog
beamforming matrix ${\bf{F}}$. Next, we transform the original problem into a
SDP and solve the digital beamforming ${\bf{v}}_{m}[k]$ by the SDR technique.
Lastly, we solve the reflection matrix ${\bf\Phi}$ according to the obtained
hybrid beamforming ${\bf{F}}$ and ${\bf{v}}_{m}[k]$. The above procedure is
repeated until convergence.

###  III-A Optimization of ${\bf{F}}$ and ${\bf{v}}_{m}[k]$ under Fixed
${\bf\Phi}$

Under given ${\bf\Phi}$, the original problem can be transformed as

$\displaystyle\underset{\left\\{{\bf{F}},{\bf{v}}_{m}[k]\right\\}}{\rm{max}}\;\sum_{m=1}^{M}\sum_{k=1}^{K}a_{m}\log\left(1\\!+\\!\frac{\left|{\hat{\bf{h}}}_{m}[k]{\bf{F}}{\bf{v}}_{m}[k]\right|^{2}}{\sum\nolimits_{j\neq
m}^{M}\left|{\hat{\bf{h}}}_{m}[k]{\bf{F}}{\bf{v}}_{j}[k]\right|^{2}\\!+\\!\delta^{2}}\right)$
(14a) $\displaystyle{\rm{s.t.}}$
$\displaystyle\;\rm(\ref{OptA2}),\rm(\ref{OptA3}),$ (14b)

where ${\hat{\bf{h}}}_{m}[k]=u_{m}[k]{{\bf{g}}}_{m}[k]{\bf\Phi}{{\bf{H}}}[k]$,
$a_{m}=B\alpha_{m}/K$, and $\delta^{2}=BN_{0}/K$.

Problem (14) is still difficult to solve since each common element of the
analog beamforming ${\bf{F}}$ needs to cater all users and subcarriers.
Furthermore, the data streams of different users may own different priority
weights, which leads to more complicated analog beamforming design. For
simplicity, we assume that all users have the same weights and it can be
regarded as a MIMO-OFDM system by neglecting the inter-user interference. To
solve the analog beamforming ${\bf{F}}$, we first define
${\hat{\bf{H}}}[k]=[{\hat{\bf{h}}}_{1}[k]^{T},...,{\hat{\bf{h}}}_{M}[k]^{T}]^{T}$,
and reformulate a MIMO-OFDM sum rate maximization problem as follows:

$\displaystyle\underset{\left\\{{\bf{F}},{\bf{V}}[k]\right\\}}{\rm{max}}\;\frac{1}{K}\sum_{k=1}^{K}\log\left({\bf{I}}\\!+\\!\frac{{\hat{\bf{H}}}[k]{\bf{F}}{\bf{V}}[k]{\bf{V}}[k]^{H}{\bf{F}}^{H}{{\hat{\bf{H}}}[k]^{H}}}{\\!\delta^{2}}\right)$
(15a) $\displaystyle{\rm{s.t.}}$
$\displaystyle\;||{\bf{F}}{\bf{V}}[k]||^{2}\leq
P_{\rm{max}}/K,\rm(\ref{OptA3}).$ (15b)

Here, we consider the transmit power constraint for each subcarrier to obtain
the low bound. Note that (15) is only used for designing ${\bf{F}}$. For a
given ${\bf{F}}$, the optimal digital beamforming matrix can be calculated as
[23]

$\displaystyle{\bf{V}}[k]=({\bf{F}}^{H}{\bf{F}})^{-{1}/{2}}{\bf{U}}_{e}[k]{\bf{\Gamma}}_{e}[k],$
(16)

where ${\bf{U}}_{e}[k]$ is the set of right singular vector according to the
$N_{\rm{RF}}$ largest singular value of
${\hat{\bf{H}}}[k]{\bf{F}}({\bf{F}}^{H}{\bf{F}})^{-{1}/{2}}$, and
${\bf{\Gamma}}_{e}[k]$ is a diagonal matrix of the power allocated to the data
streams on each subcarrier. Here, we assume that there are $N_{\rm{RF}}$ data
streams. In addition, it is obvious that ${\bf{f}}_{i}^{H}{\bf{f}}_{i}=1$,
while ${\bf{f}}_{i}^{H}{\bf{f}}_{j}\ll 1\;(i\neq j)$ with high probability for
a large $N_{\rm{TX}}$. Therefore, the analog beamforming satisfies
${\bf{F}}^{H}{\bf{F}}\approx{\bf{I}}$ that can always be approximated as
proportional to the identity matrix, namely
${\bf{F}}^{H}{\bf{F}}\propto{\bf{I}}$. Moreover, in the high or moderate
signal-to-noise ratio (SNR) regime, the equal power allocation scheme for all
streams on each subcarrier can be adopted without significant performance
degradation, namely ${\bf{\Gamma}}_{e}[k]\propto{\bf{I}}$ [24]. As a result,
the digital beamforming matrix can be approximated as
${\bf{V}}[k]\approx\lambda{\bf{U}}_{e}[k]$, where
$\lambda=\sqrt{P_{\rm{max}}/(KN_{\rm{TX}}N_{\rm{RF}})}$. Based on the above
analysis, we have

$\displaystyle\frac{1}{K}\sum_{k=1}^{K}\log\left({\bf{I}}\\!+\\!\frac{{\hat{\bf{H}}}[k]{\bf{F}}{\bf{V}}[k]{\bf{V}}[k]^{H}{\bf{F}}^{H}{{\hat{\bf{H}}}[k]^{H}}}{\\!\delta^{2}}\right)$
(17) $\displaystyle=$
$\displaystyle\frac{1}{K}\sum_{k=1}^{K}\log\left({\bf{I}}\\!+\\!\frac{\lambda}{\\!\delta^{2}}{\hat{\bf{H}}}[k]{\bf{F}}{\bf{U}}_{e}[k]{\bf{U}}_{e}[k]^{H}{\bf{F}}^{H}{{\hat{\bf{H}}}[k]^{H}}\right)$
$\displaystyle=$
$\displaystyle\frac{1}{K}\sum_{k=1}^{K}\log\left({\bf{I}}\\!+\\!\frac{\lambda}{\\!\delta^{2}}{\bf{F}}^{H}{\hat{\bf{H}}}[k]^{H}{{{\hat{\bf{H}}}[k]}\bf{F}}\right).$

Finally, we can obtain the upper bound of (17) using the Jensen’s inequality
as

$\displaystyle\frac{1}{K}\sum_{k=1}^{K}\log\left({\bf{I}}\\!+\\!\frac{\lambda}{\\!\delta^{2}}{\bf{F}}^{H}{\hat{\bf{H}}}[k]^{H}{{{\hat{\bf{H}}}[k]}\bf{F}}\right)\leq\log\left({\bf{I}}\\!+\\!\frac{\lambda}{\\!\delta^{2}}{\bf{F}}^{H}{{{{\bf{\Sigma}}}}\bf{F}}\right),$
(18)

where
${{{\bf{\Sigma}}}}=\frac{1}{K}\sum_{k=1}^{K}\left({\hat{\bf{H}}}[k]^{H}{\hat{\bf{H}}}[k]\right)$
and the analog beamforming matrix can be obtained by solving the following
problem

$\displaystyle{\bf{F}}^{\star}=\underset{{\bf{F}}(i,j)=1/\sqrt{N_{\rm{TX}}}}{\rm{arg\;max}}\;\;\;\;\;\log\left({\bf{I}}\\!+\\!\frac{\lambda}{\\!\delta^{2}}{\bf{F}}^{H}{{{{\bf{\Sigma}}}}\bf{F}}\right).$
(19)

Since ${{{\bf{\Sigma}}}}$ is a hermitian matrix, its singular value
decomposition (SVD) can be written as
${{{\bf{\Sigma}}}}={\bf{S}}{\bf{\Lambda}}{\bf{S}}^{H}.$ Therefore, the
solution of (19) can be given by
${\bf{F}}^{\star}(i,j)=\frac{1}{\sqrt{N_{\rm{TX}}}}e^{\angle\left({\bf{S}}_{1:N_{\rm{BF}}}(i,j)\right)}$,
where ${\bf{S}}_{1:N_{\rm{BF}}}$ denotes the first $N_{\rm{BF}}$ columns of
${\bf{S}}$.

Note that although (16) provides a digital beamforming solution, it is only
suitable for the single user case. Therefore, after obtaining the analog
beamforming matrix ${\bf{F}}$, problem (14) can be transformed as follows for
solving digital beamforming

$\displaystyle\underset{\left\\{{\bf{v}}_{m}[k]\right\\}}{\rm{max}}\;\sum_{m=1}^{M}\sum_{k=1}^{K}a_{m}\log\left(1\\!+\\!\frac{\left|{\bar{\bf{h}}}_{m}[k]{\bf{v}}_{m}[k]\right|^{2}}{\sum\nolimits_{j\neq
m}^{M}\left|{\bar{\bf{h}}}_{m}[k]{\bf{v}}_{j}[k]\right|^{2}\\!+\\!\delta^{2}}\right)$
(20a) $\displaystyle\;\;\;\;\;\;{\rm{s.t.}}$
$\displaystyle\;\sum\nolimits_{M=1}^{M}\sum\nolimits_{k=1}^{K}||{\bf{F}}{\bf{v}}_{m}[k]||^{2}\leq
P_{\rm{max}},$ (20b)

where ${\bar{\bf{h}}}_{m}[k]={\hat{\bf{h}}}_{m}[k]{\bf{F}}$. We define
${\bar{\bf{H}}}_{m}[k]={\bar{\bf{h}}}_{m}[k]^{H}{\bar{\bf{h}}}_{m}[k]$ and
${\bf{V}}_{m}[k]={\bf{v}}_{m}[k]{\bf{v}}_{m}[k]^{H}$, and by using an
auxiliary variable $t_{m,k}$, we reformulate (20) as

$\displaystyle\underset{\left\\{{\bf{V}}_{m}[k]\right\\}}{\rm{max}}\;\sum_{m=1}^{M}\sum_{k=1}^{K}a_{m}\log\left(1\\!+\\!t_{m,k}\right)$
(21a) $\displaystyle\;\;\;\;\;\;\;\;\;{\rm{s.t.}}$
$\displaystyle\;t_{m,k}\leq\frac{{\rm{Tr}}({\bar{\bf{H}}}_{m}[k]{\bf{V}}_{m}[k])}{\sum\nolimits_{j\neq
m}^{M}{\rm{Tr}}({\bar{\bf{H}}}_{m}[k]{\bf{V}}_{j}[k])\\!+\\!\delta^{2}},$
(21b)
$\displaystyle\;\sum\nolimits_{M=1}^{M}\sum\nolimits_{k=1}^{K}{\rm{Tr}}({\bf{F}}^{H}{\bf{F}}{\bf{V}}_{m}[k])\leq
P_{\rm{max}},$ (21c)
$\displaystyle\;{\rm{rank}}({\bf{V}}_{m}[k])=1,{\bf{V}}_{m}[k]\succeq 0.$
(21d)

It is obvious that (21) is a non-convex optimization problem due to (21b) and
(21d). To cope with (21b), we introduce an auxiliary variable $b_{m,k}$ and
transform it as

$\displaystyle
t_{m,k}b_{m,k}\leq{\rm{Tr}}({\bar{\bf{H}}}_{m}[k]{\bf{V}}_{m}[k]),$ (22a)
$\displaystyle b_{m,k}\geq\sum\nolimits_{j\neq
m}^{M}{\rm{Tr}}({\bar{\bf{H}}}_{m}[k]{\bf{V}}_{j}[k])\\!+\\!\delta^{2}.$ (22b)

Now, we only need to deal with $t_{m,k}b_{m,k}$. According to [25], the upper
bound of $t_{m,k}b_{m,k}$ can be obtained as

$\displaystyle\frac{t_{m,k}^{(n)}}{2b_{m,k}^{(n)}}b_{m,k}^{2}+\frac{b_{m,k}^{(n)}}{2t_{m,k}^{(n)}}t_{m,k}^{2}\geq
t_{m,k}b_{m,k},$ (23)

where $t_{m,k}^{(n)}$ and $b_{m,k}^{(n)}$ are the values of $t_{m,k}$ and
$b_{m,k}$ at the $n$th iteration, respectively. Consequently, we transform
(22a) into the following convex constraints

$\displaystyle\frac{t_{m,k}^{(n)}}{2b_{m,k}^{(n)}}b_{m,k}^{2}+\frac{b_{m,k}^{(n)}}{2t_{m,k}^{(n)}}t_{m,k}^{2}\leq{\rm{Tr}}({\bar{\bf{H}}}_{m}[k]{\bf{V}}_{m}[k]).$
(24)

Finally, (21) can be recast as the following SDP problem

$\displaystyle\underset{\left\\{t_{m,k},b_{m,k},{\bf{V}}_{m}[k]\right\\}}{\rm{max}}\;\sum_{m=1}^{M}\sum_{k=1}^{K}a_{m}\log\left(1\\!+\\!t_{m,k}\right)$
(25a)
$\displaystyle\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{s.t.}}$
$\displaystyle\;{\rm{(\ref{OptE2}),(\ref{eq23B}),(\ref{eq25})}},$ (25b)
$\displaystyle\;{\rm{rank}}({\bf{V}}_{m}[k])=1,{\bf{V}}_{m}[k]\succeq 0.$
(25c)

Since the rank-one constraint is non-convex, we need to drop it and formulate
a SDR problem that can be solved by existing convex solvers such as the CVX
toolbox. Summarily, to obtain the digital beamforming matrix
${\bf{V}}_{m}[k]$, we need to iteratively solve (25). Specifically, we first
initialize the auxiliary variables $b_{m,k}^{\rm{o}}$, $t_{m,k}^{\rm{o}}$ and
solve (25) for obtaining the optimal $b_{m,k}^{\star}$, $t_{m,k}^{\star}$ and
${\bf{V}}_{m}[k]^{\star}$. Next, $b_{m,k}^{\rm{o}}$ and $t_{m,k}^{\rm{o}}$ are
updated with the obtained $b_{m,k}^{\star}$ and $t_{m,k}^{\star}$, and then,
we resolve (25). The above procedure is repeated until the results converge or
the iteration number reaches its maximum value. In addition, since the SDR
problem of (25) is a convex optimization problem, the solutions are optimal
for each iteration. Therefore, iteratively solving (25) and updating variables
increase or at least maintain the value of the objective function [26, 27].
Given the limited transmit power, the designed iterative algorithm guarantees
the value of the objective function to be a monotonically non-decreasing
sequence with an upper bound, and it converges to a stationary solution that
is at least a local optimal.

For solving (25), we remove the rank-one constraint
${\rm{rank}}({\bf{V}}_{m}[k])=1$. To explore the characteristic of the
solutions, we first give the following theorem for the obtained digital
beamforming ${\bf{V}}^{\star}_{m}[k]$.

###### Theorem 1

For a large number of BS antenna $N_{\rm{TX}}$, the obtained digital
beamforming ${\bf{V}}^{\star}_{m}[k]$ satisfies
${\rm{rank}}({\bf{V}}^{\star}_{m}[k])=1$.

###### Proof:

Refer to Appendix A. ∎

In fact, even for a medium number of BS antennas, such as $N_{\rm{TX}}=16$, we
find that the optimal ${\bf{V}}^{\star}_{m}[k]$ always satisfies
${\rm{rank}}({\bf{V}}^{\star}_{m}[k])=1$. If
${\rm{rank}}({\bf{V}}^{\star}_{m}[k])=1$, the optimal
${\bf{v}}^{\star}_{m}[k]$ can be recovered by the eigenvalue decomposition,
namely
${\bf{V}}^{\star}_{m}[k]={\bf{v}}^{\star}_{m}[k]^{H}{\bf{v}}^{\star}_{m}[k]$.
If the obtained ${\bf{V}}^{\star}_{m}[k]$ is not a rank-one matrix, the
Gaussian randomization technique is typically used to obtain a rank-one
solution [28]. Specifically, we first generate a random vector
$\hat{{\bf{v}}}_{m}[k]$ satisfying
$\hat{{\bf{v}}}_{m}[k]\sim{\mathcal{CN}}({\bf{0}},{\bf{V}}^{\star}_{m}[k])$.
Then, we apply $\hat{{\bf{v}}}_{m}[k]$ to problem (25) and check its
feasibility. Here, we need to independently generate enough feasible
beamforming vectors $\hat{{\bf{v}}}_{m}[k]$ and select the optimal one
$\hat{{\bf{v}}}_{m}[k]^{\star}$ from all the random vectors. Based on this,
the final digital beamforming vector can be approximated as
${{\bf{v}}}_{m}[k]^{\star}=\hat{{\bf{v}}}_{m}[k]^{\star}$.

Now, we analyze the computational complexity for solving problem (25). Let
$\epsilon$ be the iteration accuracy; then, the number of iterations is on the
order of $\sqrt{KMN_{\rm{RF}}+5KM+1}\ln(1/\epsilon)$ [29]. For problem (25),
there are $KM$ linear matrix inequality (LMI) constraints of size
$N_{\rm{RF}}$, $3KM+1$ LMI constraints of size $1$, and $KM$ second-order cone
constraints. Therefore, the complexity of solving problem (25) is on the order
of
$\sqrt{\Delta_{1}}\zeta\ln(1/\epsilon)(\Delta_{2}+\zeta\Delta_{3}+\zeta^{2})$,
where $\Delta_{1}=KMN_{\rm{RF}}+5KM+1$, $\Delta_{2}=KMN_{\rm{RF}}^{3}+12KM+1$,
$\Delta_{3}=KMN_{\rm{RF}}^{2}+3KM+1$, and $\zeta=KMN_{\rm{RF}}^{2}$.

###  III-B Optimization of $\bf{\Phi}$ under Fixed ${\bf{F}}$ and
${\bf{V}}_{m}[k]$

After obtaining the hybrid analog/digital beamforming ${\bf{F}}$ and
${\bf{v}}_{m}[k]$, we consider the reflection matrix of IRS here and transform
(13) into the following optimization problem

$\displaystyle\underset{\left\\{{\bf{\Phi}}\right\\}}{\rm{max}}\;\sum_{m=1}^{M}\sum_{k=1}^{K}a_{m}\log\left(1\\!+\\!\frac{\left|{{\bf{g}}}_{m}[k]{\bf\Phi}{\bf{z}}_{m}[k]\right|^{2}}{\sum\nolimits_{j\neq
m}^{M}\left|{{\bf{g}}}_{m}[k]{\bf\Phi}{\bf{z}}_{j}[k]\right|^{2}\\!+\\!\delta^{2}}\right)$
(26a)
$\displaystyle{\rm{s.t.}}\;\;|\phi_{i}|=1,i\in\\{1,\cdots,N_{\rm{IRS}}\\},$
(26b)

where ${\bf{z}}_{m}[k]=u_{m}[k]{{\bf{H}}}[k]{\bf{F}}{\bf{v}}_{m}[k]$. Let
${{\bf{g}}}_{m}[k]{\bf\Phi}{\bf{z}}_{m}[k]={\bf{c}}_{m,k}{\bm{\phi}}$, where
${\bm{\phi}}=[\phi_{1},...,\phi_{N_{\rm{IRS}}}]^{T}$ and
${\bf{c}}_{m}[k]={{\bf{g}}}_{m}[k]{\rm{diag}}({\bf{z}}_{m}[k])$. Thus, we can
reformulate (26) as

$\displaystyle\underset{\left\\{{\bm{\phi}}\right\\}}{\rm{max}}\;\sum_{m=1}^{M}\sum_{k=1}^{K}a_{m}\log\left(1\\!+\\!\frac{\left|{\bf{c}}_{m}[k]{\bm{\phi}}\right|^{2}}{\sum\nolimits_{j\neq
m}^{M}\left|{\bf{c}}_{j}[k]{\bm{\phi}}\right|^{2}\\!+\\!\delta^{2}}\right)$
(27a)
$\displaystyle{\rm{s.t.}}\;\;|\phi_{i}|=1,i\in\\{1,\cdots,N_{\rm{IRS}}\\}.$
(27b)

Similarly, we define ${\bf{C}}_{m}[k]={\bf{c}}_{m}[k]^{H}{\bf{c}}_{m}[k]$ and
${\bf{\Omega}}={\bm{\phi}}{\bm{\phi}}^{H}$, and transform (27) into the
following optimization problem

$\displaystyle\underset{\left\\{{\bm{\Omega}}\right\\}}{\rm{max}}\;\sum_{m=1}^{M}\sum_{k=1}^{K}a_{m}\log\left(1\\!+\\!\frac{{\rm{Tr}}({\bf{C}}_{m}[k]{\bf{\Omega}})}{\sum\nolimits_{j\neq
m}^{M}{\rm{Tr}}({\bf{C}}_{j}[k]{\bf{\Omega}})\\!+\\!\delta^{2}}\right)$ (28a)
$\displaystyle{\rm{s.t.}}\;\;{\bf{\Omega}}(i,i)=1,i\in\\{1,\cdots,N_{\rm{IRS}}\\},$
(28b)
$\displaystyle\;\;\;\;\;\;\;{\rm{rank}}({\bf{\Omega}})=1,{\bf{\Omega}}\succeq
0.$ (28c)

One can observe that (28) has a similar form to (21), and we can adopt the
same iterative method proposed in Section. III-A to solve it. We omit the
details due to the space limitation and directly write the SDR problem as
follows:

$\displaystyle\underset{\left\\{t^{\prime}_{m,k},b^{\prime}_{m,k},{\bf{\Omega}}\right\\}}{\rm{max}}\;\sum_{m=1}^{M}\sum_{k=1}^{K}a_{m}\log\left(1\\!+\\!t^{\prime}_{m,k}\right)$
(29a) $\displaystyle{\rm{s.t.}}$
$\displaystyle\;b^{\prime}_{m,k}\geq\sum\nolimits_{j\neq
m}^{M}{\rm{Tr}}({{\bf{C}}}_{j}[k]{\bf{\Omega}})\\!+\\!\delta^{2},$ (29b)
$\displaystyle\;\frac{{t^{\prime}}_{m,k}^{(n)}}{2{b^{\prime}}_{m,k}^{(n)}}{b^{\prime}}_{m,k}^{2}+\frac{{b^{\prime}}_{m,k}^{(n)}}{2{t^{\prime}}_{m,k}^{(n)}}{t^{\prime}}_{m,k}^{2}\leq{\rm{Tr}}({\bf{C}}_{m}[k]{\bf{\Omega}}),$
(29c) $\displaystyle\;{\bf{\Omega}}(i,i)=1,{\bf{\Omega}}\succeq 0.$ (29d)

Here, we still need to alternatively solve problem (29) to obtain the
reflection matrix ${\bf{\Phi}}$. If the obtained ${\bm{\phi}}^{\star}$
satisfies ${\rm{rank}}({\bf{\Omega}}^{\star})=1$, the optimal
${\bm{\phi}}^{\star}$ can be recovered by the eigenvector of
${\bf{\Omega}}^{\star}={\bm{\phi}}^{\star}{\bm{\phi}}^{\star H}$, and the
reflection coefficient matrix ${\bf{\Phi}}^{\star}$ can be expressed as
${\bf\Phi}^{\star}={\rm{diag}}\\{\phi_{1}^{\star},...,\phi_{N_{\rm{IRS}}}^{\star}\\}$.
Otherwise, the Gaussian randomization technique can be used to obtain a rank-
one solution [28].

Now, we analyze the computational complexity for solving problem (29). Let
$\epsilon$ be the iteration accuracy; then, the number of iterations is on the
order of $\sqrt{2N_{\rm{IRS}}+5KM}\ln(1/\epsilon)$. For problem (29), there
are 1 LMI constraint of size $N_{\rm{IRS}}$, $3KM+N_{\rm{IRS}}$ LMI constraint
of size $1$, and $KM$ second-order cone constraints. Therefore, the complexity
of solving problem (29) is on the order of
$\sqrt{\Delta_{1}}\zeta\ln(1/\epsilon)(\Delta_{2}+\zeta\Delta_{3}+\zeta^{2})$,
where $\Delta_{1}=2N_{\rm{IRS}}+5KM$,
$\Delta_{2}=N_{\rm{IRS}}^{3}+N_{\rm{IRS}}+12KM$,
$\Delta_{3}=N_{\rm{IRS}}^{2}+N_{\rm{IRS}}+3KM$, and $\zeta=N_{\rm{IRS}}^{2}$.

Finally, we summarize the proposed alteratively iterative optimization scheme
in Algorithm 1.

1 Initialize the reflection matrix ${\bf{\Phi}}^{(0)}$, iteration number $r=1$
and maximum iteration number $r_{\rm{max}}$.

2 repeat

3 Obtain the analog beamforming ${\bf{F}}^{(r)}$ according to (19).

4 Initialize variables $t_{m,k}^{(0)}$, $b_{m,k}^{(0)}$, iteration number
$r^{\prime}=1$ and maximum iteration number $r^{\prime}_{\rm{max}}$.

5 repeat

6 Obtain $t_{m,k}^{\star}$, $b_{m,k}^{\star}$ and ${\bf{v}}_{m}[k]^{\star}$ by
solving (25).

7 Update variables $t_{m,k}^{(r^{\prime})}\leftarrow t_{m,k}^{\star}$,
$b_{m,k}^{(r^{\prime})}\leftarrow b_{m,k}^{\star}$.

8 Update $r^{\prime}\leftarrow r^{\prime}+1$.

9 until  _$r^{\prime}=r^{\prime}_{\rm{max}}$ or Convergence_;

10 Initialize variables ${t^{\prime}}_{m,k}^{(0)}$,
${b^{\prime}}_{m,k}^{(0)}$, iteration number $r{{}^{\prime\prime}}=1$ and
maximum iteration number $r{{}^{\prime\prime}}_{\rm{max}}$.

11 repeat

12 Obtain ${t^{\prime}}_{m,k}^{\star}$, ${b^{\prime}}_{m,k}^{\star}$ and
${\bf{\Omega}}^{\star}$ by solving (29).

13 Update variables
${t^{\prime}}_{m,k}^{(r{{}^{\prime\prime}})}\leftarrow{t^{\prime}}_{m,k}^{\star}$,
${b^{\prime}}_{m,k}^{(r{{}^{\prime\prime}})}\leftarrow{b^{\prime}}_{m,k}^{\star}$.

14 Update $r{{}^{\prime\prime}}\leftarrow r{{}^{\prime\prime}}+1$.

15 until  _$r{{}^{\prime\prime}}=r{{}^{\prime\prime}}_{\rm{max}}$ or
Convergence_;

16 Obtain ${\bf{\Phi}}^{\star}$ according to ${\bf{\Omega}}^{\star}$.

17 Update ${\bf{\Phi}}^{(r)}\leftarrow{\bf{\Phi}}^{\star}$.

18 Update $r\leftarrow r+1$.

19

20until  _$r=r_{\rm{max}}$ or Convergence_;

Obtain the analog beamforming ${\bf{\Phi}}^{(r)}$, digital beamforming
${\bf{v}}_{m}[k]^{(r)}$ and reflection matrix ${\bf{\Phi}}^{(r)}$.

Algorithm 1 The Proposed Alternatively Iterative Optimization Algorithm.

##  IV Extension to Imperfect CSIs from IRS to Users

Due to the mobility of users, it is difficult to obtain the perfect CSIs from
the IRS to users. Therefore, in this section, we assume that there exists
channel estimation error for the reflection links between IRS and users,
namely ${{\bf{g}}}_{m}[k]={\tilde{\bf{g}}}_{m}[k]+\triangle{{\bf{g}}}_{m}[k]$,
where ${\tilde{\bf{g}}}_{m}[k]$ denotes the estimated CSI and
$\triangle{{\bf{g}}}_{m}[k]$ is the estimation error. Here, we assume that the
estimation error is upper bounded with
${\tilde{\bf{g}}}_{m}[k]{\tilde{\bf{g}}}_{m}[k]^{H}\leq\varepsilon$. Next, we
redesign the hybrid analog/digital beamforming and reflection matrix with
imperfect CSI as follows

$\displaystyle\underset{\left\\{{\bf\Phi},{\bf{F}},{\bf{v}}_{m}[k]\right\\}}{\rm{max}}\;\sum\nolimits_{m=1}^{M}\alpha_{m}R_{m}$
(30a) $\displaystyle{\rm{s.t.}}$
$\displaystyle\;|\phi_{i}|=1,i\in\\{1,\cdots,N_{\rm{IRS}}\\},$ (30b)
$\displaystyle\;\sum\nolimits_{M=1}^{M}\sum\nolimits_{k=1}^{K}||{\bf{F}}{\bf{v}}_{m}[k]||^{2}\leq
P_{\rm{max}},$ (30c)
$\displaystyle\;{{\bf{g}}}_{m}[k]={\tilde{\bf{g}}}_{m}[k]+\triangle{{\bf{g}}}_{m}[k],$
(30d)
$\displaystyle\;{\tilde{\bf{g}}}_{m}[k]{\tilde{\bf{g}}}_{m}[k]^{H}\leq\varepsilon,$
(30e)
$\displaystyle\;{\bf{F}}(i,j)=1/\sqrt{N_{\rm{TX}}},i\in\\{1,\cdots,N_{\rm{TX}}\\},j\in\\{1,\cdots,N_{\rm{RF}}\\}.$
(30f)

Similar to the proposed Algorithm 1, we design an alternatively iterative
optimization scheme to solve it.

###  IV-A Optimization of ${\bf{F}}$ and ${\bf{v}}_{m}[k]$ under Fixed
${\bf\Phi}$

In fact, it is extremely difficult to design analog beamforming ${\bf{F}}$
under channel estimation error, and thus, we ignore the estimation error and
define ${{\bf{g}}}_{m}[k]\triangleq{\tilde{\bf{g}}}_{m}[k]$ for simplicity,
namely
${\hat{\bf{h}}}_{m}[k]=u_{m}[k]{\tilde{\bf{g}}}_{m}[k]{\bf\Phi}{{\bf{H}}}[k]$.
In this way, we can adopt the same scheme used in Section III-A to obtain the
analog beamforming ${\bf{F}}$. Next, we directly solve the digital beamforming
${\bf{v}}_{m}[k]$. After obtaining ${\bf{F}}$, we define
$u_{m}[k]{{\bf{g}}}_{m}[k]{\bf\Phi}{{\bf{H}}}[k]{\bf{F}}\triangleq{{\bf{g}}}_{m}[k]{\bf{\Xi}}_{m}[k]$,
where ${\bf{\Xi}}_{m}[k]=u_{m}[k]{\bf\Phi}{{\bf{H}}}[k]{\bf{F}}$. Thus, the
received signal of the $m$th user on the $k$th subcarrier can be expressed as

$\displaystyle
y_{m}[k]=({\tilde{\bf{g}}}_{m}[k]+\triangle{{\bf{g}}}_{m}[k]){\bf{\Xi}}_{m}[k]{\bf{v}}_{m}[k]x_{m}[k]+\sum\nolimits_{j\neq
m}^{M}({\tilde{\bf{g}}}_{m}[k]+\triangle{{\bf{g}}}_{m}[k]){\bf{\Xi}}_{m}[k]{\bf{v}}_{j}[k]x_{j}[k]+n_{m}[k],$
(31)

and the achievable rate is given by

$\displaystyle
R_{m}[k]=\frac{B}{K}\log\left(1\\!+\\!\frac{\left|({\tilde{\bf{g}}}_{m}[k]\\!+\\!\triangle{{\bf{g}}}_{m}[k]){\bf{\Xi}}_{m}[k]{\bf{v}}_{m}[k]\right|^{2}}{\sum\nolimits_{j\neq
m}^{M}\Upsilon_{m}[k]\\!+\\!BN_{0}/K}\right),$ (32)

where
$\Upsilon_{m}[k]=\left|({\tilde{\bf{g}}}_{m}[k]\\!+\\!\triangle{{\bf{g}}}_{m}[k]){\bf{\Xi}}_{m}[k]{\bf{v}}_{j}[k]\right|^{2}$.
By introducing an auxiliary variable $\tau_{m,k}$, the original problem (30)
can be recast as

$\displaystyle\underset{\left\\{{\bf{v}}_{m}[k]\right\\}}{\rm{max}}\;\sum_{m=1}^{M}\sum_{k=1}^{K}a_{m}\log\left(1\\!+\tau_{m,k}\right)$
(33a) $\displaystyle{\rm{s.t.}}$
$\displaystyle\;\tau_{m,k}\leq\frac{\left|({\tilde{\bf{g}}}_{m}[k]\\!+\\!\triangle{{\bf{g}}}_{m}[k]){\bf{\Xi}}_{m}[k]{\bf{v}}_{m}[k]\right|^{2}}{\sum\nolimits_{j\neq
m}^{M}\Upsilon_{m}[k]\\!+\\!BN_{0}/K},$ (33b)
$\displaystyle\;\;\triangle{{\bf{g}}}_{m}[k]\triangle{{\bf{g}}}_{m}[k]^{H}\leq\varepsilon,$
(33c)
$\displaystyle\;\sum\nolimits_{M=1}^{M}\sum\nolimits_{k=1}^{K}||{\bf{F}}{\bf{v}}_{m}[k]||^{2}\leq
P_{\rm{max}}.$ (33d)

By comparing the weighted sum rate maximization problem under perfect CSI (21)
and imperfect CSI (33), one can observe that the original scheme used to solve
(21) can not be directly used to solve (33) due to the uncertain term
$\triangle{{\bf{g}}}_{m}[k]$. To cope with (33), we first give the following
Lemma 1 referred to as $\mathcal{S}$-Procedure [30].

###### Lemma 1

Define the function

$\displaystyle
f_{i}({\bf{x}})={\bf{x}}{\bf{Q}}_{i}{\bf{x}}^{H}+2{\rm{Re}}\\{{\bf{p}}_{i}{\bf{x}}^{H}\\}+e_{i},i\in\\{1,2\\},$
(34)

where ${\bf{x}}\in{\mathbb{C}}^{1\times O}$,
${\bf{Q}}_{i}\in{\mathbb{C}}^{O\times O}$,
${\bf{p}}_{i}\in{\mathbb{C}}^{1\times O}$, and $e_{i}\in\mathbb{R}$ with $O$
representing any integer, and thus, the expression $f_{1}({\bf{x}})\leq
0\Rightarrow f_{2}({\bf{x}})\leq 0$ holds if and only if there exists a
$\beta$ satisfying

$\displaystyle\beta\left[\begin{array}[]{ccc}{\bf{Q}}_{1}&{\bf{p}}_{1}^{H}\\\
{\bf{p}}_{1}&e_{1}\end{array}\right]-\left[\begin{array}[]{ccc}{\bf{Q}}_{2}&{\bf{p}}_{2}^{H}\\\
{\bf{p}}_{2}&e_{2}\end{array}\right]\succeq{\bf{0}}.$ (39)

Next, we have

$\displaystyle\begin{aligned}
&\left|({\tilde{\bf{g}}}_{m}[k]\\!+\\!\triangle{{\bf{g}}}_{m}[k]){\bf{\Xi}}_{m}[k]{\bf{v}}_{m}[k]\right|^{2}\\\
=&\left(({\tilde{\bf{g}}}_{m}[k]\\!+\\!\triangle{{\bf{g}}}_{m}[k]){\bf{\Xi}}_{m}[k]{\bf{v}}_{m}[k]\right)\left(({\tilde{\bf{g}}}_{m}[k]\\!+\\!\triangle{{\bf{g}}}_{m}[k]){\bf{\Xi}}_{m}[k]{\bf{v}}_{m}[k]\right)^{H}\\\
=&\triangle{{\bf{g}}}_{m}[k]{\bf{\Xi}}_{m}[k]{\bf{V}}_{m}[k]{\bf{\Xi}}_{m}[k]^{H}\triangle{{\bf{g}}}_{m}[k]^{H}+2{\rm{Re}}({\tilde{\bf{g}}}_{m}[k]{\bf{\Xi}}_{m}[k]{\bf{V}}_{m}[k]{\bf{\Xi}}_{m}[k]^{H}\triangle{{\bf{g}}}_{m}[k]^{H})\\\
&+{\tilde{\bf{g}}}_{m}[k]{\bf{\Xi}}_{m}[k]{\bf{V}}_{m}[k]{\bf{\Xi}}_{m}[k]^{H}{\tilde{\bf{g}}}_{m}[k]^{H}.\end{aligned}$
(40)

Similarly, $\sum_{j\neq m}\Upsilon_{m}[k]$ can be expressed as

$\displaystyle\begin{aligned} &\sum\nolimits_{j\neq
m}\Upsilon_{m}[k]=\sum\nolimits_{j\neq
m}\left|({\tilde{\bf{g}}}_{m}[k]\\!+\\!\triangle{{\bf{g}}}_{m}[k]){\bf{\Xi}}_{m}[k]{\bf{v}}_{j}[k]\right|^{2}\\\
=&\triangle{{\bf{g}}}_{m}[k]{\bf{\Xi}}_{m}[k]\left(\sum\nolimits_{j\neq
m}{\bf{V}}_{j}[k]\right){\bf{\Xi}}_{m}[k]^{H}\triangle{{\bf{g}}}_{m}[k]^{H}+2{\rm{Re}}({\tilde{\bf{g}}}_{m}[k]{\bf{\Xi}}_{m}[k]\left(\sum\nolimits_{j\neq
m}{\bf{V}}_{j}[k]\right){\bf{\Xi}}_{m}[k]^{H}\triangle{{\bf{g}}}_{m}[k]^{H})\\\
&+{\tilde{\bf{g}}}_{m}[k]{\bf{\Xi}}_{m}[k]\left(\sum\nolimits_{j\neq
m}{\bf{V}}_{j}[k]\right){\bf{\Xi}}_{m}[k]^{H}{\tilde{\bf{g}}}_{m}[k]^{H}.\end{aligned}$
(41)

Next, we transform (33b) into the following three constraints

$\displaystyle\mu_{m,k}\leq\left|({\tilde{\bf{g}}}_{m}[k]\\!+\\!\triangle{{\bf{g}}}_{m}[k]){\bf{\Xi}}_{m}[k]{\bf{v}}_{m}[k]\right|^{2},$
(42a) $\displaystyle\tau_{m,k}\theta_{m,k}\leq\mu_{m,k},$ (42b)
$\displaystyle\theta_{m,k}\geq\sum\nolimits_{j\neq
m}^{M}\Upsilon_{m}[k]\\!+\\!BN_{0}/K.$ (42c)

Combining (40), (33c) and Lemma 1, (42a) can be transformed into the following
convex LMI constraint

$\displaystyle\left[\begin{array}[]{ccc}\beta_{m,k}{\bf{I}}\\!+\\!\hat{{\bf{V}}}_{m}[k]&({\tilde{\bf{g}}}_{m}[k]\hat{{\bf{V}}}_{m}[k])^{H}\\\
{\tilde{\bf{g}}}_{m}[k]\hat{{\bf{V}}}_{m}[k]&\varepsilon\\!+\\!{\tilde{\bf{g}}}_{m}[k]\hat{{\bf{V}}}_{m}[k]{\tilde{\bf{g}}}_{m}[k]^{H}-\mu_{m,k}\end{array}\right]\\!\succeq\\!{\bf{0}},$
(45)

where
$\hat{{\bf{V}}}_{m}[k]={\bf{\Xi}}_{m}[k]{\bf{V}}_{m}[k]{\bf{\Xi}}_{m}[k]^{H}$.
The non-convex constraint (42b) can be expressed as the following convex
constraint

$\displaystyle\frac{\tau_{m,k}^{(n)}}{2\theta_{m,k}^{(n)}}\theta_{m,k}^{2}+\frac{\theta_{m,k}^{(n)}}{2\tau_{m,k}^{(n)}}\tau_{m,k}^{2}\leq\mu_{m,k},$
(46)

where $\tau_{m,k}^{(n)}$ and $\theta_{m,k}^{(n)}$ are the values of
$\tau_{m,k}$ and $\theta_{m,k}$ at the $n$th iteration, respectively. In
addition, combining (41), (33c) and Lemma 1, (42c) can be transformed into the
following convex LMI constraint

$\displaystyle\left[\begin{array}[]{ccc}\beta_{m,k}{\bf{I}}\\!-\\!\check{{\bf{V}}}_{m}[k]&-({\tilde{\bf{g}}}_{m}[k]\check{{\bf{V}}}_{m}[k])^{H}\\\
-{\tilde{\bf{g}}}_{m}[k]\check{{\bf{V}}}_{m}[k]&\varepsilon\\!-\\!{\tilde{\bf{g}}}_{m}[k]\check{{\bf{V}}}_{m}[k]{\tilde{\bf{g}}}_{m}[k]^{H}+\theta_{m,k}-BN_{0}/K\end{array}\right]\\!\succeq\\!{\bf{0}},$
(49)

where $\check{{\bf{V}}}_{m}[k]={\bf{\Xi}}_{m}[k]\left(\sum\nolimits_{j\neq
m}{\bf{V}}_{j}[k]\right){\bf{\Xi}}_{m}[k]^{H}$.

Finally, we transform (28) into the following SDP problem

$\displaystyle\underset{\left\\{\tau_{m,k},\mu_{m,k},\theta_{m,k},\left\\{{\bf{V}}_{m}[k]\right\\}\right\\}}{\rm{max}}\;\sum_{m=1}^{M}\sum_{k=1}^{K}a_{m}\log\left(1\\!+\tau_{m,k}\right)$
(50a) $\displaystyle{\rm{s.t.}}$
$\displaystyle\;{\rm{(\ref{eq40}),(\ref{eq41}),(\ref{eq42})}},$ (50b)
$\displaystyle\;\sum\nolimits_{M=1}^{M}\sum\nolimits_{k=1}^{K}{\rm{Tr}}({\bf{F}}^{H}{\bf{F}}{\bf{V}}_{m}[k])\leq
P_{\rm{max}},$ (50c)
$\displaystyle\;{\rm{rank}}({\bf{V}}_{m}[k])=1,{\bf{V}}_{m}[k]\succeq 0.$
(50d)

It is obvious that (50) can be solved with the convex optimization toolbox
such as CVX by removing the rank-one constraint. Likewise, the optimal digital
beamforming matrix ${\bf{V}}_{m}[k]$ is obtained by iteratively solving
problem (50).

To handle problem (50), we remove the rank-one constraint
${\rm{rank}}({\bf{V}}_{m}[k])=1$. If problem (50) does not yield a rank-one
solution, namely ${\rm{rank}}({\bf{V}}_{m}[k])>1$, the Gaussian randomization
technique is used to obtain a rank-one solution [28].

Now, we analyze the computational complexity for solving problem (50). Let
$\epsilon$ be the iteration accuracy; then, the number of iterations is on the
order of $\sqrt{3KMN_{\rm{RF}}+7KM}\ln(1/\epsilon)$[29]. For problem (50),
there are $2KM$ LMI constraints of size $N_{\rm{IRS}}+1$, $KM$ LMI constraints
of size $N_{\rm{IRS}}$, $3KM+1$ LMI constraints of size 1, and $KM$ second-
order cone constraints. Therefore, the complexity of solving problem (50) is
on the order of
$\sqrt{\Delta_{1}}\zeta\ln(1/\epsilon)(\Delta_{2}+\zeta\Delta_{3}+\zeta^{2})$,
where $\Delta_{1}=3KMN_{\rm{RF}}+7KM$,
$\Delta_{2}=KM(2(N_{\rm{RF}}+1)^{3}+N_{\rm{RF}}^{3}+11)+1$,
$\Delta_{3}=2KM(N_{\rm{RF}}+1)^{2}+KMN_{\rm{RF}}^{2}+3KM+1$, and
$\zeta=KMN_{\rm{RF}}^{2}$.

###  IV-B Optimization of $\bf{\Phi}$ under Fixed ${\bf{F}}$ and
${\bf{V}}_{m}[k]$

After obtaining the analog beamforming ${\bf{F}}$ and digital beamforming
${\bf{v}}_{m}[k]$, we reconsider the reflection matrix $\bf{\Phi}$ with
imperfect CSIs. According to the definition in Section III-B, we have

$\displaystyle
u_{m}[k]{{\bf{g}}}_{m}[k]{\bf\Phi}{{\bf{H}}}[k]{\bf{F}}{\bf{v}}_{m}[k]\triangleq{{\bf{g}}}_{m}[k]{\bf\Phi}{{\bf{z}}}_{m}[k].$
(51)

Next, we define ${{\bf{Z}}}_{m}[k]\triangleq{\rm{diag}}({{\bf{z}}}_{m}[k])$,
and thus, we have
${{\bf{g}}}_{m}[k]{\bf\Phi}{{\bf{z}}}_{m}[k]={{\bf{g}}}_{m}[k]{{\bf{Z}}}_{m}[k]{\bm{\phi}}$.
After that, the achievable rate can be expressed as

$\displaystyle
R_{m}[k]=\frac{B}{K}\log\left(1\\!+\\!\frac{\left|({\tilde{\bf{g}}}_{m}[k]\\!+\\!\triangle{{\bf{g}}}_{m}[k]){{\bf{Z}}}_{m}[k]{\bm{\phi}}\right|^{2}}{\sum\nolimits_{j\neq
m}^{M}\Upsilon^{\prime}_{m}[k]\\!+\\!BN_{0}/K}\right),$ (52)

where
$\Upsilon^{\prime}_{m}[k]=\left|({\tilde{\bf{g}}}_{m}[k]\\!+\\!\triangle{{\bf{g}}}_{m}[k]){{\bf{Z}}}_{j}[k]{\bm{\phi}}\right|^{2}$.
It can be found that (52) has a similar expression with (32), and thus, the
same scheme can be used to solve the reflection matrix $\bf{\Phi}$. Therefore,
we omit the details and directly formulate the following SDR problem as

$\displaystyle\underset{\left\\{\tau^{\prime}_{m,k},\mu^{\prime}_{m,k},\theta^{\prime}_{m,k},\bf{\Omega}\right\\}}{\rm{max}}\;\sum_{m=1}^{M}\sum_{k=1}^{K}a_{m}\log\left(1\\!+\tau^{\prime}_{m,k}\right)$
(53a) $\displaystyle{\rm{s.t.}}$
$\displaystyle\;\left[\begin{array}[]{ccc}\beta^{\prime}_{m,k}{\bf{I}}\\!+\\!\hat{{\bf{\Omega}}}_{m}[k]&({\tilde{\bf{g}}}_{m}[k]\hat{{\bf{\Omega}}}_{m}[k])^{H}\\\
{\tilde{\bf{g}}}_{m}[k]\hat{{\bf{\Omega}}}_{m}[k]&\varepsilon\\!+\\!{\tilde{\bf{g}}}_{m}[k]\hat{{\bf{\Omega}}}_{m}[k]{\tilde{\bf{g}}}_{m}[k]^{H}-\mu^{\prime}_{m,k}\end{array}\right]\\!\succeq\\!{\bf{0}},$
(53d)
$\displaystyle\frac{{\tau^{\prime}}_{m,k}^{(n)}}{2{\theta^{\prime}}_{m,k}^{(n)}}{\theta^{\prime}}_{m,k}^{2}+\frac{{\theta^{\prime}}_{m,k}^{(n)}}{2{\tau^{\prime}}_{m,k}^{(n)}}{\tau^{\prime}}_{m,k}^{2}\leq\mu^{\prime}_{m,k},$
(53e)
$\displaystyle\left[\begin{array}[]{ccc}\beta^{\prime}_{m,k}{\bf{I}}\\!-\\!\check{{\bf{\Omega}}}_{m}[k]&-({\tilde{\bf{g}}}_{m}[k]\check{{\bf{\Omega}}}_{m}[k])^{H}\\\
-{\tilde{\bf{g}}}_{m}[k]\check{{\bf{\Omega}}}_{m}[k]&\varepsilon\\!-\\!{\tilde{\bf{g}}}_{m}[k]\check{{\bf{\Omega}}}_{m}[k]{\tilde{\bf{g}}}_{m}[k]^{H}+\theta^{\prime}_{m,k}-BN_{0}/K\end{array}\right]\\!\succeq\\!{\bf{0}},$
(53h) $\displaystyle\;{\bf{\Omega}}(i,i)=1,{\bf{\Omega}}\succeq 0,$ (53i)

where $\tau_{m,k}$, $\mu^{\prime}_{m,k}$ and $\theta^{\prime}_{m,k}$ are the
introduced auxiliary variables,
$\hat{{\bf{\Omega}}}_{m}[k]={\bf{Z}}_{m}[k]{\bf{\Omega}}{\bf{Z}}_{m}[k]^{H}$,
and $\check{{\bf{\Omega}}}_{m}[k]=\sum\nolimits_{j\neq
m}{\bf{\Xi}}_{j}[k]{\bf{\Omega}}{\bf{\Xi}}_{j}[k]^{H}$. We can use the CVX
toolbox to solve (53). Likewise, if problem (29) does not yield a rank-one
solution, namely ${\rm{rank}}({\bf{\Omega}}^{\star})\neq 1$. The Gaussian
randomization technique can be used to obtain a rank-one solution [28].

Now, we analyze the computational complexity for solving problem (53). Let
$\epsilon$ be the iteration accuracy; then, the number of iterations is on the
order of $\sqrt{KMN_{\rm{IRS}}+7KM}\ln(1/\epsilon)$. For problem (53), there
are $2KM$ LMI constraints of size $N_{\rm{IRS}}+1$, 1 LMI constraint of size
$N_{\rm{IRS}}$, $4KM$ LMI constraints of size $1$, and $KM$ second-order cone
constraints. Therefore, the complexity of solving problem (53) is on the order
of
$\sqrt{\Delta_{1}}\zeta\ln[1/\epsilon)(\Delta_{2}+12KM+1+\zeta\Delta_{2}+\zeta^{2}]$,
where $\Delta_{1}=KMN_{\rm{IRS}}+7KM$,
$\Delta_{2}=2KM(N_{\rm{IRS}}+1)^{3}+N_{\rm{IRS}}^{3}+13KM$,
$\Delta_{3}=2KM(N_{\rm{RF}}+1)^{2}+N_{\rm{RF}}^{2}+4KM$, and
$\zeta=N_{IRS}^{2}$.

##  V Numerical Results

In this section, simulation results are presented to evaluate the performance
of the proposed schemes in IRS-aided THz MIMO-OFDMA systems. Due to the severe
pathloss in THz, we consider a short distance communication scenario as shown
in Fig. 3, where users are located within a circle with $1.5$ m radius. The
AoD/AOA follows the uniform distribution within $[-\pi/2,\pi/2]$, and the
antenna spacing is assumed to be half wavelength. The default simulation
parameters are listed in Table I, and they are used in simulation unless
otherwise specified.

Fig. 3: The location distribution in the IRS-aided system. TABLE I: Default Parameters Used in Simulations. Parameters | Value
---|---
Number of antennas | $N_{\rm{TX}}=64$
Number of RF chains | $N_{\rm{RF}}=4$
Number of reflection elements | $N_{\rm{IRS}}=4$
Central frequency |  $f_{c}=340$ [GHz]
Bandwidth |  $B=20$ [GHz]
Number of subcarriers | $K=16$
Number of users | $M=2$
Transmit antenna gain | $G_{t}=4\\!+\\!20\log_{10}(\sqrt{N_{\rm{TX}}})$
Receive antenna gain | $G_{r}=1$
Absorption coefficient | $0.0033/m$
Speed of light | $c=3\times 10^{8}$
Noise variance |  $N_{0}=-174$ [dBm/Hz]
Fig. 4: The rate versus iteration for solving the digital beamforming. Fig. 5:
The rate versus iteration for solving the reflection matrix. Fig. 6: The rate
versus iteration for the proposed Algorithm 1. Fig. 7: The rate versus the
allowable maximum transmit power. Fig. 8: The rate versus the number of
antennas. Fig. 9: The rate versus the allowable maximum transmit power under
different numbers of IRS reflection elements. Fig. 10: The rate versus the
allowable maximum transmit power with different weights. Fig. 11: Rate
comparison under different weights.

Figs. 4 and 5 show the convergence performance of the proposed inner iterative
algorithm for solving the digital beamforming and reflection matrix,
respectively, i.e., Line $5\sim 9$ and Line $11\sim 15$ in Algorithm 1. Here,
we set $\varepsilon=0$, $P_{\rm{max}}=0$ dBm and
$\alpha_{m}=1\;(m\in\\{1,2\\})$. The legend “$n$th iteration” in Figs. 4 and 5
stands for the outer iteration number. One can observe that the inner
iterative algorithm tends to converge after 5 iterations for each outer
iteration. In addition, it can be found that the gap is small between the
$2$nd and $3$rd iterations, but large between the $1$st and $2$nd iteration.
This means that outer iterative loop (i.e., Algorithm 1) also converges
rapidly, and we will elaborate this details later in Fig. 6.

The convergence performance of Algorithm 1 under different estimation errors
is plotted in Fig. 6, where we set the maximum transmit power $P_{\rm{max}}=4$
dBm and $\alpha_{m}=1\;(m\in\\{1,2\\})$. It is clear that the rate tends to
stabilize after 3 iterations, which demonstrates the fast convergence of the
proposed algorithm. In addition, it is easy to understand that the rate is low
for a large estimation error as shown in Fig. 6, where $\varepsilon=0$ means
perfect CSIs between IRS and users.

Fig. 7 shows the rate versus $P_{\rm{max}}$ under different estimation errors,
where we set $\alpha_{m}=1\;(m\in\\{1,2\\})$. Meanwhile, we also plot the rate
under fully the digital structure, namely each antenna is connected to each RF
chain. It is clear that the rate under the fully digital structure is higher
than that under the sparse RF chain structure for the same condition, while
the circuit power consumption is very high for the former. This is also one of
the reasons for which the sparse RF chain structure is usually adopted when
the ultra high frequency carrier is applied. In addition, one can observe that
the rate increases with $P_{\rm{max}}$.

We plot the rate under different BS antennas in Fig. 8, where we set
$P_{\rm{max}}=4$ dBm and $\alpha_{m}=1\;(m\in\\{1,2\\})$. It is obvious that
the rate increases with the number of antennas, but with a decreasing slope.
In addition, we also plot the rate versus $P_{\rm{max}}$ for different numbers
of IRS reflection elements in Fig. 9, with $\varepsilon=0$. It is observed
again that a large number of IRS reflection elements leads to a higher rate.
This is because a higher beamforming gain can be achieved when there are
either more antennas at the BS or more reflection elements at the IRS.

To compare the system performance under different users’ weights, we set the
user weight $\alpha=[1\;\;1]$ and $\alpha=[1\;\;3]$ as shown in Figs. 10 and
11. $\alpha=[1\;\;1]$ means that user 1 and user 2 have the same weight,
whereas $\alpha=[1\;\;3]$ means that user 2 has a high priority than user 1.
One can observe that the rate with $\alpha=[1\;\;3]$ is always lower than that
with $\alpha=[1\;\;1]$. Furthermore, we find that the rate of user 1 is higher
than that of user 2 with $\alpha=[1\;\;1]$, while their rates are almost the
same for $\alpha=[1\;\;3]$. This means when $\alpha=[1\;\;1]$, all users have
the same priority, and the system always optimizes the power or beamforming to
maximize the sum rate. However, when user 2 has higher priority, the system
allocates more power to user 2. In this case, although it can improve the rate
of user 2, both the rate of user 1 and the sum rate are reduced. In practice,
we can set different weights according to different quality of service
requirements of the users.

##  VI Conclusion

In this paper, we have considered an IRS-aided THz MIMO-OFDMA system, where
the BS is equipped with a sparse RF chain structure. First, we have proposed a
joint hybrid analog/digital beamforming and reflection matrix design to
maximize the weighted sum rate under perfect CSIs. Next, considering the
imperfect CSIs from the IRS to users, we have redesigned a robust joint
optimization algorithm. From simulation results, we have found that the
channel estimation error has a large impact on the system sum rate. In
addition, one user with a high weight can improve this user’s rate, while the
system sum rate decreases. Consequently, channel estimation schemes and users
weight selection are important considerations for the design of practical
systems.

##  Appendix A Proof of Theorem 1

First, we give the Lagrangian function of (25) without the rank-one constraint
as

$\displaystyle\begin{aligned}
&F(t_{m,k},b_{m,k},\xi,\psi_{m,k},\nu_{m,k},{\bf{V}}_{m}[k])=\sum_{m=1}^{M}\sum_{k=1}^{K}a_{m}\log\left(1\\!+\\!t_{m,k}\right)\\\
+&\xi\left(P_{\rm{max}}-\sum\nolimits_{M=1}^{M}\sum\nolimits_{k=1}^{K}{\rm{Tr}}({\bf{F}}^{H}{\bf{F}}{\bf{V}}_{m}[k])\right)\\\
+&\sum_{m=1}^{M}\sum_{k=1}^{K}\psi_{m,k}\left(b_{m,k}-\left(\sum\nolimits_{j\neq
m}^{M}{\rm{Tr}}({\bar{\bf{H}}}_{m}[k]{\bf{V}}_{j}[k])\\!+\\!\delta^{2}\right)\right)\\\
+&\sum_{m=1}^{M}\sum_{k=1}^{K}\nu_{m,k}\left({\rm{Tr}}({\bar{\bf{H}}}_{m}[k]{\bf{V}}_{m}[k])-\left(\frac{t_{m,k}^{(n)}}{2b_{m,k}^{(n)}}b_{m,k}^{2}+\frac{b_{m,k}^{(n)}}{2t_{m,k}^{(n)}}t_{m,k}^{2}\right)\right)\\\
+&\sum_{m=1}^{M}\sum_{k=1}^{K}{\rm{Tr}}(\Theta_{m,k}{\bf{V}}_{m}[k]).\end{aligned}$
(54)

where $\xi,\psi_{m,k},\nu_{m,k}$, $\Theta_{m,k}$, respectively, represent the
Lagrange multipliers corresponding the constraints (21c), (22b), (24) and
(25c). Since the relaxed SDP problem (25) is convex, and the gap between the
primal problem and its dual problem is zero, namely it satisfies the Slater’s
condition [31]. Therefore, the Karaush-Kuhn-Tucker (KKT) conditions are
necessary and sufficient for the optimal solutions of problem (25) with rank-
one constraint. Next, we give the KKT conditions related to the optimal
digital beamforming ${\bf{V}}_{m}[k]^{\star}$ as:

$\displaystyle\xi^{\star}{\bf{F}}^{H}{\bf{F}}+\sum_{j\neq
m}^{M}\sum_{k=1}^{K}\psi_{m,k}^{\star}{\bar{\bf{H}}}_{j}[k]-\nu_{m,k}^{\star}{\bar{\bf{H}}}_{m}[k]$
$\displaystyle={\bf{\Theta}}_{m,k}^{\star},$ (55a)
$\displaystyle{\bf{\Theta}}_{m,k}^{\star}{\bf{V}}_{m}[k]^{\star}$
$\displaystyle={\bf{0}},$ (55b) $\displaystyle{\bf{\Theta}}_{m,k}^{\star}$
$\displaystyle\succeq{\bf{0}},$ (55c)

where $\xi^{\star}$, $\psi_{m,k}^{\star}$, $\nu_{m,k}^{\star}$ and
${\bf{\Theta}}_{m,k}^{\star}$ are the optimal Lagrange multipliers. The analog
beamforming can be expressed as
${\bf{F}}=[{\bf{f}}_{1},\cdots,{\bf{f}}_{\rm{RF}}]$, and we have
${\bf{f}}_{i}^{H}{\bf{f}}_{i}=1$ and ${\bf{f}}_{i}^{H}{\bf{f}}_{j}\ll
1\;(i\neq j)$ for a large $N_{\rm{TX}}$ as our analysis for the analog
beamforming matrix in Section III-A. In this way, for a large $N_{\rm{TX}}$,
we can obtain ${\bf{F}}^{H}{\bf{F}}\approx{\bf{I}}$, and
${\bf{F}}^{H}{\bf{F}}$ is a full rank matrix, namely
${\rm{rank}}({\bf{F}}^{H}{\bf{F}})=N_{\rm{RF}}$. Because $\xi^{\star}>0$,
$\psi_{m,k}^{\star}>0$, we define

$\displaystyle{\bf{Y}}=\xi^{\star}{\bf{F}}^{H}{\bf{F}}+\sum_{j\neq
m}^{M}\sum_{k=1}^{K}\psi_{m,k}^{\star}{\bar{\bf{H}}}_{j}[k],$ (56)

and thus, ${\bf{Y}}$ is a positive-definite matrix which has full rank with
${\rm{rank}}({\bf{Y}})=N_{\rm{RF}}$. Based on this, we have

Therefore, we can claim that the rank of ${\bf{\Theta}}_{m,k}^{\star}$ is
either $N_{\rm{RF}}$ or $N_{\rm{RF}}-1$. If
${\rm{rank}}({\bf{\Theta}}_{m,k}^{\star})=N_{\rm{RF}}$, according to (55b),
the optimal ${\bf{V}}_{m}[k]^{\star}={\bf{0}}$, which means that the BS does
not transmit any signal. Thus, we have
${\rm{rank}}({\bf{\Theta}}_{m,k}^{\star})=N_{\rm{RF}}-1$, and the null space
of ${\bf{\Theta}}_{m,k}^{\star}$ is one dimensional. Meanwhile, (55b) means
that ${\bf{V}}_{m}[k]^{\star}$ must lie in the null-space of
${\bf{\Theta}}_{m,k}^{\star}$, and we have
${\rm{rank}}({\bf{V}}_{m}[k]^{\star})=1$ and the proof is completed.

## References

  * [1] Cisco, “Cisco annual internet report (20182023),” _White Paper_ , 2020.
  * [2] H. Sarieddeen, M. Alouini, and T. Y. Al-Naffouri, “Terahertz-band ultra-massive spatial modulation MIMO,” _IEEE J. Sel. Areas Commun._ , vol. 37, no. 9, pp. 2040–2052, Sep. 2019.
  * [3] R. Zhang, W. Hao, G. Sun, and S. Yang, “Hybrid precoding design for wideband THz massive MIMO-OFDM systems with beam squint,” _IEEE Sys. J._ , pp. 1–4, to be published, 2020.
  * [4] S. Priebe and T. Kurner, “Stochastic modeling of THz indoor radio channels,” _IEEE Trans. Wireless Commun._ , vol. 12, no. 9, pp. 4445–4455, Sep. 2013.
  * [5] L. Lu, G. Y. Li, A. L. Swindlehurst, A. Ashikhmin, and R. Zhang, “An overview of massive mimo: Benefits and challenges,” _IEEE J. Sel. Topics Signal Proc._ , vol. 8, no. 5, pp. 742–758, May 2014.
  * [6] M. Zeng, X. Li, G. Li, W. Hao, and O. A. Dobre, “Sum rate maximization for IRS-assisted uplink NOMA,” _arXiv e-prints_ , p. arXiv:2004.10791, 2020.
  * [7] M. Zeng, E. B. Mohamed, O. A. Dobre, P. Fortier, and Q. V. Pham, “Energy-efficient resource allocation for IRS-assisted multi-antenna uplink systems,” _arXiv e-prints_ , p. arXiv:2007.10002, 2020.
  * [8] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 11, pp. 5394–5409, Nov. 2019.
  * [9] Z. Chu, W. Hao, P. Xiao, and J. Shi, “Intelligent reflecting surface aided multi-antenna secure transmission,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 1, pp. 108–112, Jan. 2020.
  * [10] X. Gao, L. Dai, S. Han, C. I, and R. W. Heath, “Energy-efficient hybrid analog and digital precoding for mmwave MIMO systems with large antenna arrays,” _IEEE J. Sel. Areas Commun._ , vol. 34, no. 4, pp. 998–1009, Apr. 2016.
  * [11] C. Lin and G. Y. Li, “Adaptive beamforming with resource allocation for distance-aware multi-user indoor terahertz communications,” _IEEE Trans. Wireless Commun._ , vol. 63, no. 8, pp. 2985–2995, Aug. 2015.
  * [12] ——, “Energy-efficient design of indoor mmwave and sub-THz systems with antenna arrays,” _IEEE Trans. Wireless Commun._ , vol. 15, no. 7, pp. 4660–4672, Jul. 2016.
  * [13] S. A. Busari, K. M. S. Huq, S. Mumtaz, J. Rodriguez, Y. Fang, D. C. Sicker, S. Al-Rubaye, and A. Tsourdos, “Generalized hybrid beamforming for vehicular connectivity using THz massive MIMO,” _IEEE Trans. Veh. Tech._ , vol. 68, no. 9, pp. 8372–8383, Sep. 2019.
  * [14] J. Tan and L. Dai, “Delay-phase precoding for THz massive MIMO with beam split,” in _Proc. IEEE GLOBECOM_ , 2019, pp. 1–6.
  * [15] H. Yuan, N. Yang, K. Yang, C. Han, and J. An, “Hybrid beamforming for terahertz MIMO-OFDM systems over frequency selective fading,” _arXiv e-prints_ , p. arXiv:1910.05967, 2019.
  * [16] B. Ning, Z. Chen, W. Chen, Y. Du, and J. Fang, “Terahertz multi-user massive mimo with intelligent reflecting surface: Beam training and hybrid beamforming,” _arXiv e-prints_ , p. arXiv:1912.11662, 2019.
  * [17] S. Zhang and R. Zhang, “Capacity characterization for intelligent reflecting surface aided MIMO communication,” _IEEE J. Sel. Areas Commun._ , to be published, 2020.
  * [18] Y. Yang, B. Zheng, S. Zhang, and R. Zhang, “Intelligent reflecting surface meets OFDM: Protocol design and rate maximization,” _IEEE Trans. Commun._ , vol. 68, no. 7, pp. 4522–4535, Jul. 2020.
  * [19] Z. Zhang and L. Dai, “A joint precoding framework for wideband reconfigurable intelligent surface-aided cell-free network,” _arXiv e-prints_ , p. arXiv:2002.03744, 2020.
  * [20] G. Zhou, C. Pan, H. Ren, K. Wang, M. Di Renzo, and A. Nallanathan, “Robust beamforming design for intelligent reflecting surface aided MISO communication systems,” _IEEE Wireless Commun. Lett._ , to be published, 2020\.
  * [21] S. Priebe, M. Kannicht, M. Jacob, and T. Kürner, “Ultra broadband indoor channel measurements and calibrated ray tracing propagation modeling at THz frequencies,” _J. Commun. and Net._ , vol. 15, no. 6, pp. 547–558, Jun. 2013.
  * [22] W. Tang, M. Z. Chen, X. Chen, J. Y. Dai, Y. Han, M. Di Renzo, Y. Zeng, S. Jin, Q. Cheng, and T. J. Cui, “Wireless communications with reconfigurable intelligent surface: Path loss modeling and experimental measurement,” _arXiv e-prints_ , p. arXiv:1911.05326, 2019.
  * [23] S. Park, A. Alkhateeb, and R. W. Heath, “Dynamic subarrays for hybrid precoding in wideband mmwave MIMO systems,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 5, pp. 2907–2920, May 2017.
  * [24] W. Shen, X. Bu, X. Gao, C. Xing, and L. Hanzo, “Beamspace precoding and beam selection for wideband millimeter-wave MIMO relying on lens antenna arrays,” _IEEE Trans. Signal Process._ , vol. 67, no. 24, pp. 6301–6313, Dec. 2019.
  * [25] P. Song, G. Scutari, F. Facchinei, and L. Lampariello, “D3m: Distributed multi-cell multigroup multicasting,” in _Proc. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2016, pp. 3741–3745.
  * [26] Q. Zhang, Q. Li, and J. Qin, “Robust beamforming for nonorthogonal multiple-access systems in MISO channels,” _IEEE Trans. Veh. Technol._ , vol. 65, no. 12, pp. 10 231–10 236, Dec. 2016.
  * [27] Z. Yang, C. Pan, W. Xu, Y. Pan, M. Chen, and M. Elkashlan, “Power control for multi-cell networks with non-orthogonal multiple access,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 2, pp. 927–942, Feb. 2018.
  * [28] Z. Luo, W. Ma, A. M. So, Y. Ye, and S. Zhang, “Semidefinite relaxation of quadratic optimization problems,” _IEEE Signal Process. Mag._ , vol. 27, no. 3, pp. 20–34, Mar. 2010.
  * [29] K. Wang, A. M. So, T. Chang, W. Ma, and C. Chi, “Outage constrained robust transmit optimization for multiuser MISO downlinks: Tractable approximations by conic optimization,” _IEEE Trans. Signal Process._ , vol. 62, no. 21, pp. 5690–5705, Nov. 2014.
  * [30] F. Zhou, Z. Li, J. Cheng, Q. Li, and J. Si, “Robust AN-aided beamforming and power splitting design for secure MISO cognitive radio with SWIPT,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 4, pp. 2450–2464, Apr. 2017.
  * [31] S. Boyd and L. Vandenberghe, “Convex optimization,” _Cambridge University Press_ , 2004.