#  Intelligent Reflecting Surface Aided MIMO Broadcasting for Simultaneous
Wireless Information and Power Transfer

Cunhua Pan, Hong Ren, Kezhi Wang, Maged Elkashlan, Arumugam Nallanathan, ,
Jiangzhou Wang, _Fellow_ , _IEEE_ and Lajos Hanzo C. Pan, H. Ren, M. Elkashlan
and A. Nallanathan are with the School of Electronic Engineering and Computer
Science at Queen Mary University of London, London E1 4NS, U.K.
(e-mail:{c.pan, h.ren, maged.elkashlan, a.nallanathan}@qmul.ac.uk). K. Wang is
with Department of Computer and Information Sciences, Northumbria University,
UK. (e-mail: kezhi.wang@northumbria.ac.uk). J. Wang is with the School of
Engineering and Digital Arts, University of Kent, Canterbury, Kent, CT2 7NZ,
U.K. (e-mail: J.Z.Wang@kent.ac.uk). L. Hanzo is with the School of Electronics
and Computer Science, University of Southampton, Southampton, SO17 1BJ, U.K.
(e-mail: lh@ecs.soton.ac.uk).

###### Abstract

An intelligent reflecting surface (IRS) is invoked for enhancing the energy
harvesting performance of a simultaneous wireless information and power
transfer (SWIPT) aided system. Specifically, an IRS-assisted SWIPT system is
considered, where a multi-antenna aided base station (BS) communicates with
several multi-antenna assisted information receivers (IRs), while guaranteeing
the energy harvesting requirement of the energy receivers (ERs). To maximize
the weighted sum rate (WSR) of IRs, the transmit precoding (TPC) matrices of
the BS and passive phase shift matrix of the IRS should be jointly optimized.
To tackle this challenging optimization problem, we first adopt the classic
block coordinate descent (BCD) algorithm for decoupling the original
optimization problem into several subproblems and alternatively optimize the
TPC matrices and the phase shift matrix. For each subproblem, we provide a
low-complexity iterative algorithm, which is guaranteed to converge to the
Karush-Kuhn-Tucker (KKT) point of each subproblem. The BCD algorithm is
rigorously proved to converge to the KKT point of the original problem. We
also conceive a feasibility checking method to study its feasibility. Our
extensive simulation results confirm that employing IRSs in SWIPT beneficially
enhances the system performance and the proposed BCD algorithm converges
rapidly, which is appealing for practical applications.

###### Index Terms:

Intelligent Reflecting Surface (IRS), Large Intelligent Surface (LIS), SWIPT,
Energy Harvesting, MIMO.

##  I Introduction

Recently, intelligent reflecting surface (IRS)-assisted wireless communication
has received considerable research attention, since it is capable of
supporting cost-effective and energy-efficient high data rate communication
for next-generation communication systems [1, 2, 3]. In simple tangible terms,
an IRS is composed of a vast number of low-cost and passive reflective
components, each of which is capable of imposing a phase change on the signals
incident upon them. Thanks to the recent advances in meta-materials [4], it
has become feasible to reconfigure the phase shifts in real time. As a result,
the phase shifts of all reflective components can be collaboratively adjusted
for ensuring that the signals reflected from the IRS can be added
constructively or destructively at the receiver in order to beneficially steer
the signal component arriving from the base station (BS) for enhancing the
desired signal power or alternatively for suppressing the undesired signals,
such as interference. In contrast to conventional physical layer techniques
that are designed for accommodating the hostile time-varying wireless channels
[5, 6], IRSs constitute a new paradigm capable of ‘reprogramming’ the wireless
propagation environment into a more favorable transmission medium. Since the
reflective components are passive, they impose a much lower power consumption
than conventional relay-aided communication systems relying on active
transmission devices. Additionally, no thermal noise is imposed by the IRS,
since it directly reflects the incident signals without decoding or amplifying
them, which is in contrast to conventional relays. Furthermore, the reflective
phase arrays can be fabricated in small size and low weight, which enables
them to be easily coated in the buildings’ facade, ceilings, walls, etc.
Furthermore, as IRS is a complementary device, it can be readily integrated
into current wireless networks without modifying the physical layer
standardization, making it transparent to the users. To fully exploit the
benefits of IRS, the active beamforming at the BS and the passive beamforming
at the IRS should be jointly designed. However, the optimization variables are
coupled and the joint design leads to a complex optimization problem that is
difficult to solve.

Some innovative efforts have been devoted to the transceiver design when
integrating IRS into various wireless communication systems, including the
single-user scenarios of [7, 8, 9, 10, 11], the downlink multiple-user
scenarios of [12, 13, 14, 15], the physical layer security design of [16, 17,
18, 19, 20, 21], the mobile edge computing (MEC) networks of [22], multigroup
multicast networks of [23] and the multicell multiuser multiple-input
multiple-output (MIMO) case in [24]. Concretely, Wu _et al._ proposed joint
active and passive beamforming for a single-user scenario in [7], where
semidefinite relaxation (SDR) was proposed for optimizing the phase shift
matrix. However, its complexity is high since the number of optimization
variables increases quadratically with the number of phase shifts.
Additionally, the Gaussian random approximation employed leads to certain
performance loss. To resolve this issue, Yu _et al._ [8] proposed a pair of
efficient algorithms termed as fixed point iteration and manifold optimization
techniques, which can guarantee locally optimal solutions. As a further
advance, the authors of [9] considered realistic frequency-selective channels.
The phase shift design was studied in [10] when only statistical channel state
information (CSI) is available. A sophisticated phase shift model was derived
in [11], by taking into account a realistic amplitude-phase relationship. For
the multiuser case, the authors in [12] considered the total transmit power
minimization problem, while guaranteeing the users’ signal-to-interference-
plus-noise ratio (SINR) constraints. The associated energy efficiency
maximization problem was studied in [13] and zero-forcing beamforming was
adopted by the BS for simplifying the optimization problem. By contrast, a
weighted sum rate (WSR) maximization problem was considered in [14] and the
fairness issues were studied in [15]. The authors of [16, 17, 18] studied the
security issues of a single-user case, while the authors of [19, 20, 21]
considered multiple-user scenarios. In [22], the IRS was shown to be
beneficial in reducing the latency of MEC networks. In addition, the IRS can
help enhance the WSR performance for the multigroup multicast network in [23].
Most recently, we considered an IRS-assisted multicell MIMO communications
scenario [24], where we demonstrated that deploying an IRS at the cell edge is
also capable of mitigating the adjacent-cell interference. Channel state
information (CSI) is challenging to obtain in IRS-assisted communication
system due to its passive feature. There are some initial efforts to handle
this issue such as channel estimation and/or robust transmission design [25,
26, 27, 28]. Specifically, Huang _et al._ [25] proposed a deep learning method
for efficient online configuration of the phase shifts, where the phase values
can be immediately obtained by inputting the user location into the trained
deep neural network. A two-stage channel estimation method based on a sparse
matrix factorization and a matrix completion was proposed in [26]. A pair of
algorithms based on compressed sensing and deep learning were conceived by
Taha _et al._ [27] for tackling the challenging channel estimation issues of
IRS-assisted systems. Most recently, we first studied the robust beamforming
design for IRS-assisted communication systems in [28], where the imperfect
channel from an IRS to users was considered and the channel estimation error
was assumed to be within a bounded elliptical region.

Figure 1: An IRS-assisted SWIPT system.

On the other hand, information transmission enabled simultaneous wireless
information and power transfer (SWIPT) is an appealing technique for future
energy-hungry Internet-of-Things (IoTs) networks. Specifically, a base station
(BS) with constant power supply will transmit wireless signals to a set of
devices. Some devices intend to decode the information from the received
signal, which are termed as information receivers (IRs), while the others will
harvest the signal energy, which are called energy receiver (ER). In [29],
Zhang _et al._ studied the trade-off between the information rate attained and
the amount of harvested energy for a single-user MIMO system. In practice, a
typical ER such as a humidity sensor requires much higher energy for its
operation than that required by IRs. Due to severe channel attenuation, the
power received by the ERs is weak, which limits the maximum link-distance of
ERs. To mitigate this issue, we propose to deploy an IRS in the vicinity of
ERs to provide additional transmission links to support the ERs for enhancing
their harvested power as shown in Fig. 1, since there is a paucity of IRS-
assisted SWIPT contributions in the literature [30]. Explicitly in [30], the
weighted sum power maximization problem was studied by Wu and Zhang, who
proved that no dedicated energy-carrying signals were required for an IRS-
aided SWIPT system. The SDR method was adopted for solving the optimization
problem, which exhibits a high computational complexity as well as imposing a
performance degradation due to the associated rank-one extraction. However,
this method is not applicable when each user is equipped with multiple
antennas. Hence, in this paper we formulate a weighted sum rate (WSR)
maximization problem for the IRS-assisted SWIPT MIMO system of Fig. 1, in
which an IRS is installed in the vicinity of ERs for compensating the
associated power loss, while maximizing the WSR of distant IRs with the aid of
passive beamforming.

Against this background, the main contributions of this paper are summarized
as follows:

  1. 1.

We formulate the WSR maximization problem by jointly optimizing the transmit
precoding (TPC) matrices of the BS and those of the passive beamforming at the
IRS for our IRS-assisted SWIPT MIMO system subject to a non-convex unit-
modulus constraint imposed on the phase shifts, while simultaneously
satisfying the energy harvesting requirement of the ERs. To the best of our
knowledge, this is the first treatise considering the WSR maximization problem
of IRS-assisted SWIPT MIMO systems, which is much more challenging than the
weighted sum power minimization problem of [30] since the latter can be
readily transformed into a convex optimization problem. In contrast to the
multicell system of [24], an additional energy harvesting constraint is also
imposed in our current study, which further complicates the analysis.
Specifically, this constraint is non-convex and the optimization problem may
become infeasible. The WSR maximization problem is challenging to solve, since
the optimization variables are highly coupled and the data rate expressions of
the IRs are complex. To deal with this issue, we first reformulate the
original problem into an equivalent form by exploiting the equivalence between
the data rate and the weighted minimum mean-square error (WMMSE). Then, an
alternating optimization algorithm based on the popular block coordinate
descent (BCD) algorithm is proposed for alternatively updating the active TPC
matrices of the BS and the phase shift matrix of the IRS, which is rigorously
proved to converge to the Karush-Kuhn-Tucker (KKT) point of the original
optimization problem.

  2. 2.

For a given phase shift matrix, we then proceed by developing an iterative
algorithm based on the successive convex approximation (SCA) method and on the
Lagrangian dual decomposition method to derive a nearly closed-form solution
for the TPC matrices. A low-complexity bisection search method is proposed for
finding the optimal dual variables. The solutions generated by our iterative
algorithm are guaranteed to converge to the KKT point of the TPC optimization
problem.

  3. 3.

For the given TPC matrices, we formulate the phase shift optimization problem
as a non-convex quadratically constrained quadratic program (QCQP) subject to
an additional energy harvesting constraint by invoking some further matrix
manipulations. Then, a novel iterative algorithm based on the majorization-
minimization (MM) algorithm [31] and on the price-based method [32] is
developed for solving the QCQP. We strictly prove that the final solution
generated by the iterative algorithm is guaranteed to converge to the KKT
point of the phase shift optimization problem.

  4. 4.

The associated feasibility issue is also studied by formulating an alternative
optimization problem and an iterative algorithm is proposed for solving this
problem.

  5. 5.

Extensive simulation results are provided for verifying the performance
advantages of employing IRS in SWIPT in order to enhance the energy harvesting
performance. It is shown that the operating range of the ERs can be
dramatically expanded by placing IRSs in the ERs’ vicinity. Furthermore, the
BCD algorithm converges rapidly, and it is eminently suitable for practical
applications. Our simulation results also show that as expected, the path loss
exponent substantially affects the system’s performance and thus the location
of the IRS should be carefully chosen.

The remainder of this paper is organized as follows. In Section II, we
introduce the IRS-assisted SWIPT system model and our problem formulation. The
detailed algorithms used for solving the optimization problem are presented in
Section III. The feasibility issues of the original problem are discussed in
Section IV, followed by our extensive simulations and discussions in Section
V. Finally, our conclusions are provided in Section VI.

_Notations_ : For matrix $\bf{A}$, $\bf{A}^{*}$ and $\bf{A}^{\star}$ represent
the conjugate operator and converged solution, respectively.
${\rm{Re}}\\{a\\}$ represents the real part of a complex value $a$.
${{\mathbb{C}}^{M}}$ denotes the set of $M\times 1$ complex vectors.
${{\mathbb{E}}}\\{\cdot\\}$ denotes the expectation operation. For two
matrices $\bf A$ and $\bf B$, ${\bf{A}}\odot{\bf{B}}$ represents Hadamard
product of $\bf A$ and $\bf B$. ${\left\|{\bf{A}}\right\|_{F}}$,
${\rm{tr}}\left({\bf{A}}\right)$ and $\left|{\bf{A}}\right|$ denote the
Frobenius norm, trace operation and determinant of ${\bf{A}}$, respectively.
$\nabla{f_{\bf{x}}}\left({\bf{x}}\right)$ denotes the gradient of the function
$f$ with respect to (w.r.t.) the vector ${\bf{x}}$. ${\cal C}{\cal
N}({\bf{0}},{\bf{I}})$ represents a random vector following the distribution
of zero mean and unit variance matrix. $\arg\\{\cdot\\}$ means the extraction
of phase information. ${\rm{diag}}(\cdot)$ denotes the diagonalization
operation. ${\left(\cdot\right)^{*}}$, ${\left(\cdot\right)^{\rm{T}}}$ and
${\left(\cdot\right)^{\rm{H}}}$ denote the conjugate, transpose and Hermitian
operators, respectively. ${\rm{arg}}(\cdot)$ means the phase extraction
operation.

##  II System Model and Problem Formulation

###  II-A System Model

Consider the IRS-aided multiuser MIMO downlink of a SWIPT system operating
over the same frequency band both for data and energy transmission, as shown
in Fig. 1. Let us assume that there are $K_{I}$ IRs and $K_{E}$ ERs,
respectively. It is also assumed that the BS is equipped with $N_{B}\geq 1$
antennas, while each IR and ER is equipped with $N_{I}\geq 1$ and $N_{E}\geq
1$ antennas, respectively. Let us denote the sets of IRs and ERs as ${\cal
K}_{I}$ and ${\cal K}_{E}$, respectively. In general, low-power sensors
require a certain amount of power (e.g., 0.1 mW) for their real-time
operation. Due to the associated severe channel attenuation, the sensors
should be deployed sufficiently close to the BS, which limits their practical
implementation. To resolve this issue, we propose to employ an IRS, which has
$M$ reflective elements in the ERs’ vicinity for extending the operational
range of sensors, as shown in Fig. 1. Firstly, the IRS increases the energy
harvested by the ERs, and additionally it also assists in enhancing the signal
strength for distant IRs through careful phase shift optimization.

The number of data streams destined for each IR is assumed to be $d$,
satisfying $1\leq d\leq\min\\{N_{B},N_{I}\\}$. The signal transmitted from the
BS is given by

${{\bf{x}}}=\sum\limits_{k=1}^{K_{I}}{{\bf{F}}_{k}{\bf{s}}_{k}},$ (1)

where ${\bf{s}}_{k}\in\mathbb{C}^{d\times 1}$ is the $(d\times 1)$-element
data symbol vector designated for the $k$th IR satisfying
$\mathbb{E}\left[{{{\bf{s}}_{k}}{{{\bf{s}}_{k}^{H}}}}\right]={\bf{I}}_{d}$ and
$\mathbb{E}\left[{{{\bf{s}}_{i}}{{{\bf{s}}_{j}^{H}}}}\right]={\bf{0}},{\rm{for}}\
i\neq j$, while ${\bf{F}}_{k}\in\mathbb{C}^{N_{B}\times d}$ is the linear TPC
matrix used by the BS for the $k$th IR. Assuming non-dispersive narrow-band
transmission, the baseband equivalent channels spanning from the BS to the
IRS, from the BS to the $k$th IR, from the BS to the $l$th ER, from the IRS to
the $k$th IR, and finally from the IRS to the $l$th ER are modelled by the
matrices ${\bf{Z}}\in\mathbb{C}^{M\times N_{B}}$,
${\bf{H}}_{b,k}\in\mathbb{C}^{N_{I}\times N_{B}}$,
${\bf{G}}_{b,l}\in\mathbb{C}^{N_{E}\times N_{B}}$,
${\bf{H}}_{r,k}\in\mathbb{C}^{N_{I}\times M}$, and
${\bf{G}}_{r,l}\in\mathbb{C}^{N_{E}\times M}$, respectively. Let us denote the
diagonal reflection-coefficient matrix at the IRS by
${\bm{\Phi}}={\rm{diag}}\left\\{{{e^{j{\theta_{1}}}},\cdots,{e^{j{\theta_{m}}}},\cdots,{e^{j{\theta_{M}}}}}\right\\}$
111$j$ is the imaginary unit., where ${\theta_{m}}\in[0,2\pi]$ is the phase
shift of the $m$-th reflective element. Due to absorption and diffraction, the
signal power that has been reflected multiple times is ignored. As a result,
the signal received at the $k$th IR is given by

${{\bf{y}}_{I,k}}=\left({{{\bf{H}}_{b,k}}+{{\bf{H}}_{r,k}}{\bf{\Phi}}{\bf{Z}}}\right){\bf{x}}+{{\bf{n}}_{I,k}},\vspace{-0.3cm}$
(2)

where ${\bf{n}}_{I,k}$ is the $k$th IR’s noise vector satisfying ${\cal
C}{\cal N}\left({{\bf{0}},\sigma_{I}^{2}{{\bf{I}}_{N_{I}}}}\right)$.
Similarly, the signal received at the $l$th ER is given by

${{\bf{y}}_{E,l}}=\left({{{\bf{G}}_{b,l}}+{{\bf{G}}_{r,l}}{\bf{\Phi}}{\bf{Z}}}\right){\bf{x}}+{{\bf{n}}_{E,l}},\vspace{-0.2cm}$
(3)

where ${\bf{n}}_{E,l}$ is the $l$th ER’s noise vector obeying the distribution
of ${\cal C}{\cal N}\left({{\bf{0}},\sigma_{E}^{2}{{\bf{I}}_{N_{E}}}}\right)$.

We assume that all the CSIs are perfectly known at the BS, and the BS is
responsible for calculating the phase shifts of the IRS, which are then fed
back by them to the IRS controller through dedicated feedback channels. Given
this idealized and simplified assumption, the results obtained represent a
performance upper bound of how much performance gain can be achieved by an
IRS. Let us define the equivalent channel spanning from the BS to the $k$th IR
by
${\bf{\bar{H}}}_{k}\buildrel\Delta\over{=}{{{\bf{H}}_{b,k}}+{{\bf{H}}_{r,k}}{\bf{\Phi}}{\bf{Z}}}$.
Upon substituting $\bf{x}$ into (2), ${{\bf{y}}_{I,k}}$ can be rewritten as

${{\bf{y}}_{I,k}}={{{\bf{\bar{H}}}}_{k}}{{\bf{F}}_{k}}{{\bf{s}}_{k}}+\sum\limits_{i=1,i\neq
k}^{{K_{I}}}{{{{\bf{\bar{H}}}}_{k}}{{\bf{F}}_{i}}{{\bf{s}}_{i}}}+{{\bf{n}}_{I,k}}.\vspace{-0.2cm}$
(4)

Then, the achievable data rate (nat/s/Hz) of the $k$th IR is given by [33]

${R_{k}}\left({{\bf{F}},{\bm{\Phi}}}\right)={\log}\left|{{\bf{I}}+{{{\bf{\bar{H}}}}_{k}}{{\bf{F}}_{k}}{\bf{F}}_{k}^{\rm{H}}{\bf{\bar{H}}}_{k}^{\rm{H}}{\bf{J}}_{k}^{-1}}\right|,\vspace{-0.3cm}$
(5)

where ${\bf{F}}$ denotes the collection of TPC matrices, while
${{\bf{J}}_{k}}$ is the interference-plus-noise covariance matrix given by
${{\bf{J}}_{k}}=\sum\nolimits_{m=1,m\neq
k}^{K_{I}}{{{{\bf{\bar{H}}}}_{k}}{{\bf{F}}_{m}}{\bf{F}}_{m}^{\rm{H}}{\bf{\bar{H}}}_{k}^{\rm{H}}}+{\sigma_{I}^{2}}{\bf{I}}$.

On the other hand, due to the broadcast nature of wireless channels, the ERs
can extract energy from the electromagnetic wave. In general, the harvested
power is nonlinear over the received radio frequency (RF) power due to the
nonlinear RF-to-DC conversion, which depends on the input RF power level. This
nonlinear EH model has been characterized in [34], which is a complex function
of the RF power. Based on this nonlinear EH model, various transmission
designs have been proposed in [35] and [36]. However, there is still lack of a
general model that can accurately characterize this nonlinear relationship by
capturing all practical factors. Hence, for simplicity, we adopt the simple
linear EH model as widely used in the existing literature [29, 37, 38]. By
ignoring the noise power at the ERs, the total harvested power is proportional
to the total received power. Let us define the equivalent channel spanning
from the BS to the $l$th ER by
${\bf{\bar{G}}}_{l}\buildrel\Delta\over{=}{{{\bf{G}}_{b,l}}+{{\bf{G}}_{r,l}}{\bf{\Phi}}{\bf{Z}}}$.
Then, the total power harvested by the $l$th ER is

${Q_{i}}=\eta{\rm{tr}}\left({\sum\limits_{k=1}^{K_{I}}{{{{\bf{\bar{G}}}}_{l}}{{\bf{F}}_{k}}{\bf{F}}_{k}^{\rm{H}}{\bf{\bar{G}}}_{l}^{\rm{H}}}}\right),\vspace{-0.2cm}$
(6)

where $0<\eta\leq 1$ is the energy harvesting efficiency. In this paper, we
consider the constraint that the weighted sum of the power harvested by all
ERs should be higher than a predefined value, which is

$Q=\sum\limits_{l=1}^{{K_{E}}}{{\alpha_{l}}{Q_{l}}}={\rm{tr}}\left({\sum\limits_{k=1}^{K_{I}}{{\bf{F}}_{k}^{\rm{H}}{\bf{G}}{{\bf{F}}_{k}}}}\right)\geq\bar{Q},$
(7)

where
${\bf{G}}=\sum\nolimits_{l=1}^{{K_{E}}}{{\alpha_{l}}\eta{\bf{\bar{G}}}_{l}^{\rm{H}}{{{\bf{\bar{G}}}}_{l}}}$,
${\alpha_{l}}$ is the energy weighting factor of the $l$th ER, with a higher
value of ${\alpha_{l}}$ representing a higher priority for the $l$th ER than
for others. Finally, $\bar{Q}$ is the minimum harvested power threshold.

###  II-B Problem Formulation

Upon introducing the notations of $\phi_{m}={e^{j{\theta_{m}}}},\forall m$, we
have ${\bm{\Phi}}={\rm{diag}}\left\\{{\phi_{1},\cdots,\phi_{M}}\right\\}$.
Again, we aim for jointly optimizing the TPC matrices ${\bf{F}}$ and phase
shift matrix ${\bm{\Phi}}$ with the goal of maximizing the WSR of all IRs
subject to the total power budget, to the unit modulus of the phase shifters
and to the harvested power requirement. Then, this problem can be formulated
as follows:

$\displaystyle\mathop{\max}\limits_{{{\bf{F}},{\bm{\Phi}}}}\quad$
$\displaystyle{\sum\limits_{k=1}^{K_{I}}{{\omega_{k}}{R_{k}}\left({{\bf{F}},{\bm{\Phi}}}\right)}}$
(8a) s.t.
$\displaystyle\sum\limits_{k=1}^{K_{I}}{\left\|{{{\bf{F}}_{k}}}\right\|_{F}^{2}}\leq{P_{T}},$
(8b)
$\displaystyle{\rm{tr}}\left({\sum\limits_{k=1}^{K_{I}}{{\bf{F}}_{k}^{\rm{H}}{\bf{G}}{{\bf{F}}_{k}}}}\right)\geq\bar{Q},$
(8c) $\displaystyle\left|{{\phi_{m}}}\right|=1,m=1,\cdots,M,\vspace{-0.2cm}$
(8d)

where $\omega_{k}$ is the weighting factor controlling the scheduling priority
for each IR and ${P_{T}}$ is the power limit at the BS, while (8d) is the
unit-norm constraint imposed on the phase shifters.

As the IRS is passive and both the ERs and IRs are energy constrained, we
assume that this optimization problem is solved at the BS which posses the
knowledge of the CSI of all related links and other related parameters such as
$\bar{Q}$. After computing the phase shift values for the IRS, they are sent
to the IRS controller through dedicated control channels. Problem (8) is
difficult to solve, since the TPC matrices and the phase shifts are coupled.
If we remove the energy harvesting (EH) constraint, the problem reduces to the
WSR maximization problem recently studied in [24]. However, the additional EH
constraint makes the optimization more challenging to solve and the algorithms
developed in [24] cannot be directly applied for two reasons. Firstly, the EH
constraint is non-convex. Secondly, this problem may be infeasible due to the
conflicting constraints (8b) and (8c). In the following, we first conceive a
low-complexity algorithm to solve this problem by assuming that it is
feasible. Then, we study the feasibility of this problem.

##  III Low-Complexity Algorithm Development

In this section, we first transform Problem (8) into a more tractable one,
which allows the decoupling of the TPC matrices and of the phase shifts. Then,
the classic block coordinate descent (BCD) algorithm [33] is proposed for
solving the transformed problem.

###  III-A Reformulation of the Original Problem

To deal with the complex objective function, we reformulate Problem (8) by
employing the well-known WMMSE method [39]. The appealing feature of this
method is that it can transform the original complex problem into an
equivalent form, which facilitates the application of the BCD method.

Specifically, the linear decoding matrix ${\bf{U}}$ is applied to estimate the
signal vector ${{{\bf{\hat{s}}}}_{k}}$ for each IR

${{{\bf{\hat{s}}}}_{k}}={\bf{U}}_{k}^{\rm{H}}{{\bf{y}}_{I,k}},\forall k,$ (9)

where ${\bf{U}}_{k}\in\mathbb{C}^{N_{I}\times d}$ is the decoding matrix of
the $k$th IR. Then, the MSE matrix of the $k$th IR is given by

$\displaystyle{{\bf{E}}_{k}}$ $\displaystyle=$
$\displaystyle{\mathbb{E}_{{\bf{s,n}}}}\left[{\left({{{{\bf{\hat{s}}}}_{k}}-{{\bf{s}}_{k}}}\right){{\left({{{{\bf{\hat{s}}}}_{k}}-{{\bf{s}}_{k}}}\right)}^{H}}}\right]$
(10) $\displaystyle=$
$\displaystyle\\!\\!\left({{\bf{U}}_{k}^{\rm{H}}{{{\bf{\bar{H}}}}_{k}}{{\bf{F}}_{k}}-{\bf{I}}}\right)\\!\\!{\left({{\bf{U}}_{k}^{\rm{H}}{{{\bf{\bar{H}}}}_{k}}{{\bf{F}}_{k}}\\!-\\!{\bf{I}}}\right)^{\rm{H}}}+\\!\\!\\!\sum\limits_{m=1,m\neq
k}^{K_{I}}\\!\\!\\!{{\bf{U}}_{k}^{\rm{H}}{{{\bf{\bar{H}}}}_{k}}{{\bf{F}}_{m}}{\bf{F}}_{m}^{\rm{H}}{\bf{\bar{H}}}_{k}^{\rm{H}}{{\bf{U}}_{k}}}+{\sigma^{2}}{\bf{U}}_{k}^{\rm{H}}{{\bf{U}}_{k}},\forall
k\in{\cal K}_{I},$ (11)

where ${\bf{s}}$ and ${\bf{n}}$ denote the collections of data symbols and
noise vectors of all IRs, respectively.

By introducing a set of auxiliary matrices
${\bf{W}}=\\{{{\bf{W}}_{k}}\succeq{\bf{0}},\forall k\in{\cal K}_{I}\\}$ and
defining ${\bf{U}}=\\{{{\bf{U}}_{k}},\forall k\in{\cal K}_{I}\\}$, Problem (8)
can be reformulated as follows [33, 39]:

$\displaystyle\mathop{\max}\limits_{{{\bf{W}},{\bf{U}},{\bf{F}},{\bm{\Phi}}}}\quad$
$\displaystyle\sum\limits_{k=1}^{K_{I}}{{\omega_{k}}{h_{k}}\left({{\bf{W}},{\bf{U}},{\bf{F}},{\bm{\Phi}}}\right)}$
(12a) s.t. $\displaystyle(\ref{dehwifr}),(\ref{hhudiuwe}),(\ref{jiofjj}),$
(12b)

where ${h_{k}}\left({{\bf{W}},{\bf{U}},{\bf{F}},{\bm{\Phi}}}\right)$ is given
by

${h_{k}}\left({{\bf{W}},{\bf{U}},{\bf{F}},{\bm{\Phi}}}\right)={\log\left|{{{\bf{W}}_{k}}}\right|-{\rm{Tr}}\left({{{\bf{W}}_{k}}{{\bf{E}}_{k}}}\right)+d}.$
(13)

Although Problem (12) has more optimization variables than Problem (8), the
objective function (OF) in Problem (12) is much easier to handle, which allows
the BCD algorithm to solve this problem by iteratively obtaining one set of
variables while keeping the others fixed. Note that the decoding matrices
${\bf{U}}$ and the auxiliary matrices $\bf{W}$ only appear in the function
${h_{k}}\left({{\bf{W}},{\bf{U}},{\bf{F}},\bm{\Phi}}\right)$. Hence, the
optimal solution of ${\bf{U}}$ and $\bf{W}$ can be obtained while keeping the
other matrices fixed. Specifically, given ${\bm{\Phi}}$, ${\bf{W}}$, and
${\bf{F}}$, setting the first-order derivative of
${h_{k}}\left({{\bf{W}},{\bf{U}},{\bf{F}},\bm{\Phi}}\right)$ with respect to
${\bf{U}}_{k}$ and ${\bf{W}}_{k}$ to zero, we can obtain the optimal solution
of ${\bf{U}}_{k}$ and ${\bf{W}}_{k}$ respectively as follows

${\bf{U}}_{k}^{\star}={\left({{{\bf{J}}_{k}}+{{{\bf{\bar{H}}}}_{k}}{{\bf{F}}_{k}}{\bf{F}}_{k}^{\rm{H}}{\bf{\bar{H}}}_{k}^{\rm{H}}}\right)^{-1}}{{{\bf{\bar{H}}}}_{k}}{{\bf{F}}_{k}},{\bf{W}}_{k}^{\star}={\bf{E}}_{k}^{\star-1},$
(14)

where ${\bf{E}}_{k}^{\star}$ is obtained by inserting ${\bf{U}}_{k}^{\star}$
into the $k$th IR’s MSE matrix in (11), yielding

${\bf{E}}_{k}^{\star}={{\bf{I}}_{d}}-{\bf{F}}_{k}^{\rm{H}}{\bf{\bar{H}}}_{k}^{\rm{H}}{\left({\sum\limits_{m=1}^{{K_{I}}}{{{{\bf{\bar{H}}}}_{k}}{{\bf{F}}_{m}}{\bf{F}}_{m}^{\rm{H}}{\bf{\bar{H}}}_{k}^{\rm{H}}}+\sigma_{I}^{2}{\bf{I}}}\right)^{-1}}{{{\bf{\bar{H}}}}_{k}}{{\bf{F}}_{k}}.$
(15)

In the following, we focus our attention on the optimization of TPC matrices
${\bf{F}}$ and phase shifts ${\bm{\Phi}}$, when ${\bf{U}}$ and ${\bf{W}}$ are
given.

###  III-B Optimizing the Precoding Matrices ${\bf{F}}$

In this subsection, we aim to optimize the TPC matrices ${\bf{F}}$ with fixed
${\bf{W}},{\bf{U}}$ and $\bm{\Phi}$. By inserting ${\bf{E}}_{k}$ in (11) into
the OF of (12) and discarding the constant terms, the TPC matrices of our
optimization problem can be transformed as follows

$\displaystyle\mathop{\min}\limits_{{{\bf{F}}}}\quad\sum\limits_{k=1}^{K_{I}}{{\rm{tr}}\left({{\bf{F}}_{k}^{\rm{H}}{\bf{A}}{{\bf{F}}_{k}}}\right)}-\sum\limits_{k=1}^{K_{I}}{{\omega_{k}}{\rm{Tr}}\left({{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{{\bf{\bar{H}}}}_{k}}{{\bf{F}}_{k}}}\right)}-\sum\limits_{k=1}^{K_{I}}{{\omega_{k}}{\rm{tr}}\left({{{\bf{W}}_{k}}{\bf{F}}_{k}^{\rm{H}}{\bf{\bar{H}}}_{k}^{\rm{H}}{{\bf{U}}_{k}}}\right)}$
(16a) $\displaystyle\textrm{s.t.}\quad(\ref{dehwifr}),(\ref{hhudiuwe}),$ (16b)

where
${{\bf{A}}=\sum\nolimits_{m=1}^{K_{I}}{{\omega_{m}}}{\bf{\bar{H}}}_{m}^{\rm{H}}{{\bf{U}}_{m}}{{\bf{W}}_{m}}{\bf{U}}_{m}^{\rm{H}}{{{\bf{\bar{H}}}}_{m}}}$.

However, due to the non-convexity of the EH constraint, Problem (16) is still
non-convex. To resolve this issue, we observe that it can be viewed as a
difference of convex (d.c.) program, which can be efficiently solved by the
successive convex approximation (SCA) method [40]. In particular, we can
approximate it by its first-order Taylor expansion. By applying [41, Appendix
B] and Jensen’ inequality, we have

${\rm{tr}}\left({\sum\limits_{k=1}^{{K_{I}}}{{\bf{F}}_{k}^{\rm{H}}{\bf{G}}{{\bf{F}}_{k}}}}\right)\geq-{\rm{tr}}\left({\sum\limits_{k=1}^{{K_{I}}}{{\bf{F}}_{k}^{(n){\rm{H}}}{\bf{GF}}_{k}^{(n)}}}\right)+2{\mathop{\rm
Re}\nolimits}\left[{{\rm{tr}}\left({\sum\limits_{k=1}^{{K_{I}}}{{\bf{F}}_{k}^{(n){\rm{H}}}{\bf{G}}{{\bf{F}}_{k}}}}\right)}\right],$
(17)

where $\left\\{{{\bf{F}}_{k}^{(n)},\forall k}\right\\}$ is the solution
obtained from the previous iteration. Then, upon replacing the constraint (8c)
by the following constraint:

$2{\mathop{\rm
Re}\nolimits}\left[{{\rm{tr}}\left({\sum\limits_{k=1}^{{K_{I}}}{{\bf{F}}_{k}^{(n){\rm{H}}}{\bf{G}}{{\bf{F}}_{k}}}}\right)}\right]\geq\tilde{Q},$
(18)

where
$\tilde{Q}=\bar{Q}+{\rm{tr}}\left({\sum\nolimits_{k=1}^{{K_{I}}}{{\bf{F}}_{k}^{(n){\rm{H}}}{\bf{GF}}_{k}^{(n)}}}\right)$,
we may consider the following optimization problem:

$\displaystyle\mathop{\min}\limits_{{{\bf{F}}}}\quad\sum\limits_{k=1}^{K_{I}}{{\rm{tr}}\left({{\bf{F}}_{k}^{\rm{H}}{\bf{A}}{{\bf{F}}_{k}}}\right)}-\sum\limits_{k=1}^{K_{I}}{{\omega_{k}}{\rm{tr}}\left({{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{{\bf{\bar{H}}}}_{k}}{{\bf{F}}_{k}}}\right)}-\sum\limits_{k=1}^{K_{I}}{{\omega_{k}}{\rm{tr}}\left({{{\bf{W}}_{k}}{\bf{F}}_{k}^{\rm{H}}{\bf{\bar{H}}}_{k}^{\rm{H}}{{\bf{U}}_{k}}}\right)}$
(19a) $\displaystyle\textrm{s.t.}\quad(\ref{dehwifr}),(\ref{fsfrg}).$ (19b)

Since the OF is convex w.r.t. $\bf{F}$, and the constraints (8b) and (18) are
convex, Problem (19) constitutes a convex optimization problem, which can be
solved by standard convex solver packages, such as CVX [42]. However, the
resultant computational complexity is high. In the following, we provide a
low-complexity algorithm for obtaining a nearly optimal closed-form solution
by resorting to the Lagrangian dual decomposition method [43]. Since Problem
(19) is a convex problem and satisfies the slater’s condition222 According to
line 1 in Algorithm 2, the initial precoding matrix is initialized by the
solution obtained from Section IV. Assume the original problem is feasible.
Due to the randomness of channel matrices of $\bf{G}$ and $\bf{H}$, the
precoding matrix obtained in Section IV must be strictly larger than the
minimum EH requirement, i.e.,
${\rm{tr}}\left({\sum\nolimits_{k=1}^{K_{I}}{{\bf{F}}_{k}^{(0)\rm{H}}{\bf{G}}{{\bf{F}}_{k}^{(0)}}}}\right)>\bar{Q}$.
Then, based on [29], there must exist a strictly feasible solution, and thus
the slater’s condition holds. , the dual gap is zero and the optimal solution
can be obtained by solving its dual problem instead of its original one. We
first introduce the Lagrange multiplier $\lambda$ associated with the power
constraint, and derive the partial Lagrangian function of Problem (19) as
follows

$\begin{array}[]{l}{\cal
L}\left({{{\bf{F}}},{\lambda}}\right)=\sum\limits_{k=1}^{K_{I}}{{\rm{tr}}\left({{\bf{F}}_{k}^{\rm{H}}{{\bf{A}}}{{\bf{F}}_{k}}}\right)}-\sum\limits_{k=1}^{K_{I}}{{\omega_{k}}{\rm{tr}}\left({{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{{\bf{\bar{H}}}}_{k}}{{\bf{F}}_{k}}}\right)}-\sum\limits_{k=1}^{K_{I}}{{\omega_{k}}{\rm{tr}}\left({{{\bf{W}}_{k}}{\bf{F}}_{k}^{\rm{H}}{\bf{\bar{H}}}_{k}^{\rm{H}}{{\bf{U}}_{k}}}\right)}\\\
\qquad\qquad\qquad\quad+{\lambda}\sum\limits_{k=1}^{K_{I}}{\rm{tr}}\left({{\bf{F}}_{k}^{\rm{H}}{{\bf{F}}_{k}}}\right)-{\lambda}{P_{T}}.\end{array}$
(20)

The dual function can be obtained by solving the following problem

$g\left(\lambda\right)\buildrel\Delta\over{=}\mathop{\min}\limits_{\bf{F}}{\cal
L}({\bf{F}},\lambda)\ \ \ \ {\rm{s}}{\rm{.t}}{\rm{.}}\ (\ref{fsfrg}).$ (21)

Then, the dual problem is given by

$\displaystyle\mathop{\max}\limits_{{\lambda}}\quad g\left(\lambda\right)$
(22a) $\displaystyle\textrm{s.t.}\quad\lambda\geq 0.$ (22b)

Before solving the dual problem (22), we have to derive the expression of the
dual function $g\left(\lambda\right)$ by solving Problem (21) for a given
$\lambda$. By introducing the dual variable $\mu\geq 0$ associated with the
constraint (18), the Lagrangian function for Problem (21) is given by

$\begin{array}[]{l}{\cal
L}\left({{{\bf{F}}},{\mu}}\right)=\sum\limits_{k=1}^{K_{I}}{{\rm{tr}}\left({{\bf{F}}_{k}^{\rm{H}}\left({{{\bf{A}}}+{\lambda}{\bf{I}}}\right){{\bf{F}}_{k}}}\right)}-\sum\limits_{k=1}^{K_{I}}{{\omega_{k}}{\rm{tr}}\left({{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{{\bf{\bar{H}}}}_{k}}{{\bf{F}}_{k}}}\right)}-\sum\limits_{k=1}^{K_{I}}{{\omega_{k}}{\rm{tr}}\left({{{\bf{W}}_{k}}{\bf{F}}_{k}^{\rm{H}}{\bf{\bar{H}}}_{k}^{\rm{H}}{{\bf{U}}_{k}}}\right)}\\\
\qquad\qquad\qquad\quad+\mu\tilde{Q}-2\mu{\mathop{\rm
Re}\nolimits}\left[{{\rm{tr}}\left({\sum\limits_{k=1}^{{K_{I}}}{{\bf{F}}_{k}^{(n){\rm{H}}}{\bf{G}}{{\bf{F}}_{k}}}}\right)}\right]-{\lambda}{P_{T}}.\end{array}$
(23)

By setting the first-order derivative of ${\cal
L}\left({{{\bf{F}}},{\mu}}\right)$ w.r.t. ${{{\bf{F}}_{k}^{*}}}$ to the zero
matrix, we obtain the optimal solution of ${\bf{F}}_{k}$ as follows:

${{\bf{F}}_{k}^{\star}}(\mu)={\left({{\bf{A}}+\lambda{\bf{I}}}\right)^{{\dagger}}}\left({{\omega_{k}}{\bf{\bar{H}}}_{k}^{\rm{H}}{{\bf{U}}_{k}}{{\bf{W}}_{k}}+\mu{\bf{GF}}_{k}^{(n)}}\right),$
(24)

where ${(\cdot)^{\dagger}}$ denotes the matrix pseudoinverse. The value of
$\mu$ should be chosen for ensuring that the complementary slackness condition
for constraint (18) is satisfied:

$\mu\left(2{\mathop{\rm
Re}\nolimits}\left[{{\rm{tr}}\left({\sum\limits_{k=1}^{{K_{I}}}{{\bf{F}}_{k}^{(n){\rm{H}}}{\bf{G}}{{\bf{F}}_{k}^{\star}(\mu)}}}\right)}\right]-\tilde{Q}\right)=0.$
(25)

Hence, if the following condition holds

$2{\mathop{\rm
Re}\nolimits}\left[{{\rm{tr}}\left({\sum\limits_{k=1}^{{K_{I}}}{{\bf{F}}_{k}^{(n){\rm{H}}}{\bf{G}}{{\bf{F}}_{k}^{\star}(0)}}}\right)}\right]\geq\tilde{Q},$
(26)

the optimal solution of Problem (21) is given by
${{\bf{F}}_{k}^{\star}}(0),\forall k\in{\cal K}_{I}$. Otherwise, the optimal
$\mu$ is

$\mu=\frac{{\tilde{Q}-2{\mathop{\rm
Re}\nolimits}\left[{{\rm{tr}}\left({\sum\limits_{k=1}^{{K_{I}}}{{\omega_{k}}{\bf{F}}_{k}^{(n){\rm{H}}}{\bf{G}}{{\left({{\bf{A}}+\lambda{\bf{I}}}\right)}^{-1}}{\bf{\bar{H}}}_{k}^{\rm{H}}{{\bf{U}}_{k}}{{\bf{W}}_{k}}}}\right)}\right]}}{{2{\rm{tr}}\left({\sum\limits_{k=1}^{{K_{I}}}{{\bf{F}}_{k}^{(n){\rm{H}}}{\bf{G}}{{\left({{\bf{A}}+\lambda{\bf{I}}}\right)}^{-1}}{\bf{GF}}_{k}^{(n)}}}\right)}}.$
(27)

With the aid of the dual function, we may now commence the solution of the
dual problem (22) to find the optimal $\lambda$. Given $\lambda$, we denote
the optimal solution of Problem (21) by ${{\bf{F}}_{k}}(\lambda)$. The value
of $\lambda$ should be chosen for ensuring that the complementary slackness
condition for the power constraint is satisfied:

$\lambda\left({{\rm{tr}}\left(\sum\limits_{k=1}^{K_{I}}{{\bf{F}}_{k}^{\rm{H}}\left(\lambda\right){{\bf{F}}_{k}}\left(\lambda\right)}\right)-{P_{T}}}\right)=0.$
(28)

If the following condition holds:

${{\rm{tr}}\left(\sum\limits_{k=1}^{K_{I}}{{\bf{F}}_{k}^{\rm{H}}\left(0\right){{\bf{F}}_{k}}\left(0\right)}\right)\leq{P_{T}}},$
(29)

then the optimal solution is given by ${\bf{F}}_{k}(0)$. Otherwise, we have to
find $\lambda$ for ensuring that the following equation holds:

$P(\lambda)\buildrel\Delta\over{=}{\rm{tr}}\left(\sum\limits_{k=1}^{K_{I}}{{\bf{F}}_{k}^{\rm{H}}\left(\lambda\right){{\bf{F}}_{k}}\left(\lambda\right)}\right)={P_{T}}.$
(30)

Unfortunately, due to the complex expression of $\mu$ in (27), we are unable
to prove its monotonic nature by using the explicit expression of $P(\lambda)$
as in [24]. In the following lemma, we prove that $P(\lambda)$ is a
monotonically decreasing function of $\lambda$, which enables the bisection
search method to find $\lambda$.

Lemma 1: The total power $P(\lambda)$ is a monotonically decreasing function
of $\lambda$.

Proof: Please refer to Appendix A. $\Box$

Based on Lemma 1, the bisection search method can be used for finding the
solution of equation (30). In Algorithm 1, we provide the detailed steps of
solving Problem (19) for the case of $\lambda>0$. In each iteration of
Algorithm 1, we have to calculate ${{\bf{F}}_{k}^{\star}}(\mu)$ in (24), which
involves the calculation of ${\left({{\bf{A}}+\lambda{\bf{I}}}\right)}^{-1}$
at a complexity order of ${\cal O}(N_{B}^{3})$. If the total number of
iterations is $T$, then the total complexity of calculating
${\left({{\bf{A}}+\lambda{\bf{I}}}\right)}^{-1}$ is ${\cal O}(TN_{B}^{3})$,
which may be excessive. Here, we provide a method for reducing the
computational complexity. Specifically, as ${{\bf{A}}}$ is a non-negative
definite matrix, it can be decomposed as
${{\bf{A}}}={\bf{Q}}{\bm{\Lambda}}{{\bf{Q}}^{\rm{H}}}$ by using the singular
value decomposition (SVD), where
${\bf{Q}}{{\bf{Q}}^{\rm{H}}}={{\bf{Q}}^{\rm{H}}}{\bf{Q}}={{\bf{I}}_{{N_{T}}}}$
and ${\bm{\Lambda}}$ is a diagonal matrix with non-negative diagonal elements.
Then, we have
${\left({{\bf{A}}+\lambda{\bf{I}}}\right)}^{-1}={\bf{Q}}\left({\lambda{\bf{I}}+{\bm{\Lambda}}}\right)^{-1}{{\bf{Q}}^{\rm{H}}}$.
Hence, in each iteration, we only have to calculate the product of two
matrices, which has much lower complexity than calculating the inverse of the
matrix having the same dimension.

Algorithm 1 Bisection Search Method to Solve Problem (19)

1: Initialize the accuracy $\varepsilon$, the bounds $\lambda_{l}$ and
$\lambda_{u}$;

2: Calculate
$\lambda={{\left({{\lambda_{l}}+{\lambda_{u}}}\right)}\mathord{\left/{\vphantom{{\left({{\lambda_{l}}+{\lambda_{u}}}\right)}2}}\right.\kern-1.2pt}2}$;

3: If condition (26) is satisfied, $\mu$ is equal to zero. Otherwise, update
$\mu$ in (27);

4: Calculate $\\{{\bf{F}}_{k}(\lambda),\forall k\\}$ according to (24);

5: If $P(\lambda)\geq P_{T}$, set ${\lambda_{l}}={\lambda}$. Otherwise, set
${\lambda_{u}}={\lambda}$;

6: If $\left|{{\lambda_{l}}-{\lambda_{u}}}\right|\leq\varepsilon$, terminate.
Otherwise, go to step 2.

Based on the above discussions, in Algorithm 2 we provide the detailed steps
of the SCA algorithm conceived for solving Problem (16).

Algorithm 2 SCA Algorithm to Solve Problem (16)

1: Initialize the accuracy $\varepsilon$, the precoding matrices
${\bf{F}}^{(0)}$ from Section 2, the iteration index $n=0$, the maximum number
of iterations $n_{\rm{max}}$, calculate the OF value of Problem (16) as
$z({\bf{F}}^{(0)})$;

2: Calculate
$\tilde{Q}^{(n)}=\bar{Q}+{\rm{tr}}\left({\sum\nolimits_{k=1}^{{K_{I}}}{{\bf{F}}_{k}^{(n){\rm{H}}}{\bf{GF}}_{k}^{(n)}}}\right)$;

3: With $\tilde{Q}^{(n)}$, calculate $\\{{\bf{F}}_{k}^{(n+1)},\forall k\\}$ by
solving Problem (19) using Algorithm 1;

4: If $n\geq n_{\rm{max}}$ or
${{\left|{z({{\bf{F}}^{(n+1)}})-z({{\bf{F}}^{(n)}})}\right|}\mathord{\left/{\vphantom{{\left|{z({{\bf{F}}^{(n+1)}})-z({{\bf{F}}^{(n)}})}\right|}{\left|{z({{\bf{F}}^{(n+1)}})}\right|}}}\right.\kern-1.2pt}{\left|{z({{\bf{F}}^{(n+1)}})}\right|}}<\varepsilon$,
terminate. Otherwise, set $n\leftarrow n+1$ and go to step 2.

In the following, we show that Algorithm 2 converges to the KKT point of
Problem (16).

Theorem 1: The sequences of $\\{{\bf{F}}^{(n)},n=1,2,\cdots\\}$ generated by
Algorithm 2 converge to the KKT optimum point of Problem (16).

Proof: The proof is similar to that of [44] and hence it is omitted for
simplicity. $\Box$

Next, we briefly analyze the complexity of Algorithm 2. We assume that
$N_{B}\geq N_{I}\geq d$. In each iteration of Algorithm 2, the main complexity
contribution is the calculation of $\\{{\bf{F}}_{k}^{(n+1)},\forall k\\}$ by
using the bisection search method in Algorithm 1. In each iteration of
Algorithm 1, the main complexity lies in calculating $\bf{F}$ in (24), which
is on the order of ${\cal O}(K_{I}N_{B}^{3})$. The number of iterations
required for Algorithm 1 to converge is given by
${\log_{2}}\left({\frac{{{\lambda_{u}}-{\lambda_{l}}}}{\varepsilon}}\right)$.
Hence, the total complexity of Algorithm 1 is ${\cal
O}({\log_{2}}\left({\frac{{{\lambda_{u}}-{\lambda_{l}}}}{\varepsilon}}\right)K_{I}N_{B}^{3})$.
Then, the total complexity of Algorithm 2 is given by ${\cal
O}(n_{\rm{max}}{\log_{2}}\left({\frac{{{\lambda_{u}}-{\lambda_{l}}}}{\varepsilon}}\right)K_{I}N_{B}^{3})$.

###  III-C Optimizing the Phase Shift Matrix

In this subsection, we focus our attention on optimizing the phase shift
matrix ${\bm{\Phi}}$, while fixing the other parameters. Upon substituting
${\bf{E}}_{k}$ in (11) into (13) and removing the terms that are independent
of ${\bm{\Phi}}$, the phase shift optimization problem is formulated as:

$\displaystyle{\mathop{\min}\limits_{{\bm{\Phi}}}\quad\\!\\!\\!\\!\sum\limits_{k=1}^{{K_{I}}}{{\rm{tr}}\left(\\!{{\omega_{k}}{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{{\bf{\bar{H}}}}_{k}}{{\bf{\tilde{F}}}}{\bf{\bar{H}}}_{k}^{\rm{H}}{{\bf{U}}_{k}}}\\!\right)}\\!\\!-\\!{\sum\limits_{k=1}^{K_{I}}{{\rm{tr}}\left({\omega_{k}}{{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{{\bf{\bar{H}}}}_{k}}{{\bf{F}}_{k}}}\right)}}\\!-\\!{\sum\limits_{k=1}^{K_{I}}{{\rm{tr}}\left({\omega_{k}}{{{\bf{W}}_{k}}{\bf{F}}_{k}^{\rm{H}}{\bf{\bar{H}}}_{k}^{\rm{H}}{{\bf{U}}_{k}}}\right)}}}$
(31a) $\displaystyle\textrm{s.t.}\quad(\ref{hhudiuwe}),(\ref{jiofjj}),$ (31b)

where
${\bf{\tilde{F}}}=\sum\nolimits_{m=1}^{K_{I}}{{\bf{F}}_{m}}{\bf{F}}_{m}^{\rm{H}}$.

By substituting
${\bf{\bar{H}}}_{k}={{{\bf{H}}_{b,k}}+{{\bf{H}}_{r,k}}{\bf{\Phi}}{\bf{Z}}}$
into (31a), we have

$\begin{array}[]{l}\\!\\!{\omega_{k}}{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{{\bf{\bar{H}}}}_{k}}{{\bf{\tilde{F}}}}{\bf{\bar{H}}}_{k}^{\rm{H}}{{\bf{U}}_{k}}\\!=\\!{\omega_{k}}{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{\bf{H}}_{r,k}}{\bm{\Phi}}{{\bf{Z}}}{{\bf{\tilde{F}}}}{\bf{Z}}^{\rm{H}}{{\bm{\Phi}}^{\rm{H}}}{\bf{H}}_{r,k}^{\rm{H}}{{\bf{U}}_{k}}\\!+\\!{\omega_{k}}{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{\bf{H}}_{b,k}}{{\bf{\tilde{F}}}}{\bf{Z}}^{\rm{H}}{{\bm{\Phi}}^{\rm{H}}}{\bf{H}}_{r,k}^{\rm{H}}{{\bf{U}}_{k}}\\\
\qquad\qquad\quad\qquad\qquad\quad+{\omega_{k}}{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{\bf{H}}_{r,k}}{\bm{\Phi}}{{\bf{Z}}}{{\bf{\tilde{F}}}}{\bf{H}}_{b,k}^{\rm{H}}{{\bf{U}}_{k}}+{\omega_{k}}{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{\bf{H}}_{b,k}}{{\bf{\tilde{F}}}}{\bf{H}}_{b,k}^{\rm{H}}{{\bf{U}}_{k}},\end{array}$
(32)

and

${\omega_{k}}{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{{\bf{\bar{H}}}}_{k}}{{\bf{F}}_{k}}={\omega_{k}}{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{\bf{H}}_{r,k}}{\bm{\Phi}}{{\bf{Z}}}{{\bf{F}}_{k}}+{\omega_{k}}{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{\bf{H}}_{b,k}}{{\bf{F}}_{k}}.$
(33)

Let us define
${{\bf{B}}_{k}}\buildrel\Delta\over{=}{\omega_{k}}{\bf{H}}_{r,k}^{\rm{H}}{{\bf{U}}_{k}}{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{\bf{H}}_{r,k}}$,
${{\bf{C}}}\buildrel\Delta\over{=}{{\bf{Z}}}{{\bf{\tilde{F}}}}{\bf{Z}}^{\rm{H}}$
and
${{\bf{D}}_{k}}\buildrel\Delta\over{=}{\omega_{k}}{{\bf{Z}}}{\bf{\tilde{F}}}^{\rm{H}}{\bf{H}}_{b,k}^{\rm{H}}{{\bf{U}}_{k}}{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{\bf{H}}_{r,k}}$.
By using (32), we arrive at:

${\rm{tr}}\left({\omega_{k}}{{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{{\bf{\bar{H}}}}_{k}}{{\bf{\tilde{F}}}}{\bf{\bar{H}}}_{k}^{\rm{H}}{{\bf{U}}_{k}}}\right)={\rm{tr}}\left({{{\bm{\Phi}}^{\rm{H}}}{{\bf{B}}_{k}}{\bm{\Phi}}{{\bf{C}}}}\right)+{\rm{tr}}\left({{{\bm{\Phi}}^{\rm{H}}}{{\bf{D}}_{k}^{\rm{H}}}}\right)+{\rm{tr}}\left({{\bm{\Phi}}{\bf{D}}_{k}}\right)+{\rm{const}_{1}},$
(34)

where ${\rm{const}_{1}}$ is a constant term that is independent of
${\bm{\Phi}}$.

Similarly, by defining
${{\bf{T}}_{k}}\buildrel\Delta\over{=}{\omega_{k}}{{\bf{Z}}}{{\bf{F}}_{k}}{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{\bf{H}}_{r,k}}$,
from (33) we have

${\rm{tr}}\left({{\omega_{k}}{{\bf{W}}_{k}}{\bf{U}}_{k}^{\rm{H}}{{{\bf{\bar{H}}}}_{k}}{{\bf{F}}_{k}}}\right)={\rm{tr}}\left({\bm{\Phi}}{{\bf{T}}_{k}}\right)+{\rm{const}}_{2},$
(35)

where ${\rm{const}}_{2}$ is a constant term that is independent of
${\bm{\Phi}}$.

By defining
${{\bf{G}}_{b}}\buildrel\Delta\over{=}\sum\nolimits_{l=1}^{{K_{E}}}{{\alpha_{l}}\eta}{\bf{G}}_{b,l}^{\rm{H}}{{\bf{G}}_{b,l}}$,
${{\bf{G}}_{r}}\buildrel\Delta\over{=}\sum\nolimits_{l=1}^{{K_{E}}}{{\alpha_{l}}\eta}{\bf{G}}_{r,l}^{\rm{H}}{{\bf{G}}_{r,l}}$,
and
${{\bf{G}}_{br}}\buildrel\Delta\over{=}{\bf{Z\tilde{F}}}\sum\nolimits_{l=1}^{{K_{E}}}{{\alpha_{l}}\eta}{\bf{G}}_{b,l}^{\rm{H}}{{\bf{G}}_{r,l}}$,
the EH constraint in (8c) can be recast as follows:

${\rm{tr}}\left({{{\bm{\Phi}}^{\rm{H}}}{{\bf{G}}_{r}}{\bm{\Phi}}{\bf{C}}}\right)+{\rm{tr}}\left({{{\bm{\Phi}}^{\rm{H}}}{\bf{G}}_{b,r}^{\rm{H}}}\right)+{\rm{tr}}\left({{\bm{\Phi}}{{\bf{G}}_{br}}}\right)+{\rm{tr}}\left({{{\bf{G}}_{b}}{\bf{\tilde{F}}}}\right)\geq\bar{Q}.$
(36)

By inserting (34) and (35) into the OF of Problem (31) and removing the
constant terms, we have

$\displaystyle{\mathop{\min}\limits_{{\bm{\Phi}}}\quad{\rm{tr}}\left({{{\bm{\Phi}}^{\rm{H}}}{\bf{B}}{\bm{\Phi}}{\bf{C}}}\right)+{\rm{tr}}\left({{{\bm{\Phi}}^{\rm{H}}}{{\bf{V}}^{\rm{H}}}}\right)+{\rm{tr}}\left({{\bm{\Phi}}{\bf{V}}}\right)}$
(37a) $\displaystyle\textrm{s.t.}\quad(\ref{jiofjj}),(\ref{acvfeatgb}),$ (37b)

where ${\bf{B}}$ and ${\bf{V}}$ are given by
${\bf{B}}={\sum\nolimits_{k=1}^{K_{I}}{{{\bf{B}}_{k}}}}$ and
${\bf{V}}=\sum\nolimits_{k=1}^{K_{I}}{{{\bf{D}}_{k}}}-\sum\nolimits_{k=1}^{K_{I}}{{{\bf{T}}_{k}}}$,
respectively.

Upon denoting the collection of diagonal elements of ${\bm{\Phi}}$ by
${\bm{\phi}}={\left[{{\phi_{1}},\cdots,{\phi_{M}}}\right]^{\rm{T}}}$ and
adopting the matrix identity of [45, Eq. (1.10.6)], it follows that

${\rm{tr}}\left({{{\bm{\Phi}}^{\rm{H}}}{\bf{B}}{\bm{\Phi}}{\bf{C}}}\right)={{\bm{\phi}}^{\rm{H}}}\left({{\bf{B}}\odot{{\bf{C}}^{\rm{T}}}}\right){\bm{\phi}},{\rm{tr}}\left({{{\bm{\Phi}}^{\rm{H}}}{{\bf{G}}_{r}}{\bm{\Phi}}{\bf{C}}}\right)={{\bm{\phi}}^{\rm{H}}}\left({{\bf{G}}_{r}\odot{{\bf{C}}^{\rm{T}}}}\right){\bm{\phi}}.$
(38)

Upon denoting the collections of diagonal elements of ${\bf{V}}$ and
${\bf{G}}_{br}$ by
${\bf{v}}={\left[{{{\left[{\bf{V}}\right]}_{1,1}},\cdots,{{\left[{\bf{V}}\right]}_{M,M}}}\right]^{\rm{T}}}$
and
${\bf{g}}={\left[{{{\left[{\bf{G}}_{br}\right]}_{1,1}},\cdots,{{\left[{\bf{G}}_{br}\right]}_{M,M}}}\right]^{\rm{T}}}$,
we arrive at

${\rm{tr}}\left({{\bm{\Phi}}{\bf{V}}}\right)={{\bf{v}}^{\rm{T}}}{{\bm{\phi}}},{\rm{tr}}\left({{{\bm{\Phi}}^{\rm{H}}}{{\bf{V}}^{\rm{H}}}}\right)={\bm{\phi}}^{\rm{H}}{\bf{v}^{*}},{\rm{tr}}\left({{\bm{\Phi}}{{\bf{G}}_{br}}}\right)={{\bf{g}}^{\rm{T}}}{{\bm{\phi}}},{\rm{tr}}\left({{{\bm{\Phi}}^{\rm{H}}}{{{\bf{G}}_{br}^{\rm{H}}}}}\right)={\bm{\phi}}^{\rm{H}}{\bf{g}^{*}}.$
(39)

Moreover, the constraint (36) can be rewritten as

${{\bm{\phi}}^{\rm{H}}}\bm{\Upsilon}{\bm{\phi}}+2{\mathop{\rm
Re}\nolimits}\left\\{{{{\bm{\phi}}^{\rm{H}}}{\bf{g}^{*}}}\right\\}\geq{\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\frown$}}\over{Q}}},$ (40)

where we have ${\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\frown$}}\over{Q}}}=\bar{Q}-{\rm{Tr}}\left({{{\bf{G}}_{b}}{\bf{\tilde{F}}}}\right)$
and $\bm{\Upsilon}={{\bf{G}}_{r}\odot{{\bf{C}}^{\rm{T}}}}$. It can be verified
that ${\bf{G}}_{r}$ and ${{\bf{C}}^{\rm{T}}}$ are non-negative semidefinite
matrices. Then, according to [45], the Hadamard product
${{\bf{G}}_{r}\odot{{\bf{C}}^{\rm{T}}}}$ (or equivalently $\bm{\Upsilon}$) is
also a semidefinite matrix.

Thus, Problem (37) can be transformed as

$\displaystyle{\mathop{\min}\limits_{\bm{\phi}}\quad{{\bm{\phi}}^{\rm{H}}}{\bm{\Xi}}{\bm{\phi}}+2{\mathop{\rm
Re}\nolimits}\left\\{{{{\bm{\phi}}^{\rm{H}}}{\bf{v}^{*}}}\right\\}}$ (41a)
$\displaystyle\textrm{s.t.}\quad(\ref{jiofjj}),(\ref{acdfwefe}),$ (41b)

where we have $\bm{\Xi}={{\bf{B}}\odot{{\bf{C}}^{\rm{T}}}}$. Again, ${\bf{B}}$
can be verified to be a non-negative semidefinite matrix, and thus
${\bm{\Xi}}$ is a non-negative semidefinite matrix.

Due to the non-convex constraint (40), Problem (41) is difficult to solve. To
deal with this constraint, we again employ the SCA method [40]. Specifically,
since ${{\bm{\phi}}^{\rm{H}}}\bm{\Upsilon}{\bm{\phi}}$ is convex w.r.t.
${\bm{\phi}}$, its lower bound can be obtained as follows:

${{\bm{\phi}}^{\rm{H}}}\bm{\Upsilon}{\bm{\phi}}\geq-{{\bm{\phi}}^{{(n)}\rm{H}}}\bm{\Upsilon}{\bm{\phi}^{(n)}}+2{\mathop{\rm
Re}\nolimits}\left[{{\bm{\phi}}^{\rm{H}}}\bm{\Upsilon}{\bm{\phi}^{(n)}}\right],$
(42)

where ${{\bm{\phi}}^{{(n)}}}$ is obtained in the previous iteration. Then,
constraint (40) is replaced by the following constraint

$2{\mathop{\rm
Re}\nolimits}\left[{{{\bm{\phi}}^{\rm{H}}}\left({{\bf{g}^{*}}+{\bm{\Upsilon}}{{\bm{\phi}}^{(n)}}}\right)}\right]\geq{\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\frown$}}\over{Q}}+{{\bm{\phi}}^{(n){\rm{H}}}}\bm{\Upsilon}{{\bm{\phi}}^{(n)}}}\buildrel\Delta\over{=}\hat{Q},$
(43)

which is a linear constraint. Then, Problem (41) then becomes

$\displaystyle{\mathop{\min}\limits_{\bm{\phi}}\quad{{\bm{\phi}}^{\rm{H}}}{\bm{\Xi}}{\bm{\phi}}+2{\mathop{\rm
Re}\nolimits}\left\\{{{{\bm{\phi}}^{\rm{H}}}{\bf{v}^{*}}}\right\\}}$ (44a)
$\displaystyle\textrm{s.t.}\quad(\ref{jiofjj}),(\ref{sacvfgatr}).$ (44b)

In the following, we conceive the Majorization-Minimization (MM) algorithm
[31] for solving Problem (44). The key idea is to solve a challenging problem
by introducing a series of more tractable subproblems. Upon denoting the
objective function of Problem (44) by $f({\bm{\phi}})$, in the $(n+1)$th
iteration we have to find the upper bound of the OF, denoted as
$g({\bm{\phi}}|{\bm{\phi}}^{(n)})$, which should satisfy the following three
conditions:

$1)g({\bm{\phi}}^{(n)}|{\bm{\phi}}^{(n)})\\!=\\!f({\bm{\phi}}^{(n)});2){\left.{{\nabla_{\bm{\phi}^{*}}}g({\bm{\phi}}|{{\bm{\phi}}^{(n)}})}\right|_{{\bm{\phi}}={{\bm{\phi}}^{(n)}}}}={\left.{{\nabla_{\bm{\phi}^{*}}}f({\bm{\phi}})}\right|_{{\bm{\phi}}={{\bm{\phi}}^{(n)}}}};3)g({\bm{\phi}}|{\bm{\phi}}^{(n)})\\!\geq\\!f({\bm{\phi}}).$
(45)

Then, we solve the approximate subproblem defined by a more tractable new OF
$g({\bm{\phi}}|{\bm{\phi}}^{(n)})$. To find
$g({\bm{\phi}}|{\bm{\phi}}^{(n)})$, we introduce the following lemma [46].

Lemma 2: For any given ${\bm{\phi}}^{(n)}$, the following inequality holds for
any feasible ${\bm{\phi}}$:

${{\bm{\phi}}^{\rm{H}}}{\bm{\Xi}}{\bm{\phi}}\leq{{\bm{\phi}}^{\rm{H}}}{\bf{X}}{\bm{\phi}}-2{\mathop{\rm
Re}\nolimits}\left\\{{{{\bm{\phi}}^{\rm{H}}}\left({{\bf{X}}-{\bm{\Xi}}}\right){{\bm{\phi}}^{(n)}}}\right\\}+{\left({{{\bm{\phi}}^{(n)}}}\right)^{\rm{H}}}\left({{\bf{X}}-{\bm{\Xi}}}\right){{\bm{\phi}}^{(n)}}\buildrel\Delta\over{=}y({\bm{\phi}}|{{\bm{\phi}}^{(n)}}),$
(46)

where ${\bf{X}}={\lambda_{\rm{\max}}}{{\bf{I}}_{M}}$ and
${\lambda_{\rm{\max}}}$ is the maximum eigenvalue of ${\bm{\Xi}}$. $\Box$

Then, the function $g({\bm{\phi}}|{\bm{\phi}}^{(n)})$ can be constructed as
follows:

$g({\bm{\phi}}|{\bm{\phi}}^{(n)})=y({\bm{\phi}}|{{\bm{\phi}}^{(n)}})+2{\rm{Re}}\left\\{{{{\bm{\phi}}^{\rm{H}}}{{\bf{v}}^{*}}}\right\\},$
(47)

where $y({\bm{\phi}}|{{\bm{\phi}}^{(n)}})$ is defined in (46). The new OF
$g({\bm{\phi}}|{\bm{\phi}}^{(n)})$ is more tractable than the original OF
$f({\bm{\phi}})$. The subproblem to be solved is given by

$\displaystyle{\mathop{\min}\limits_{\bm{\phi}}\quad
g({\bm{\phi}}|{\bm{\phi}}^{(n)})}$ (48a)
$\displaystyle\textrm{s.t.}\quad(\ref{jiofjj}),(\ref{sacvfgatr}).$ (48b)

Since ${\bm{\phi}}^{\rm{H}}{\bm{\phi}}=M$, we have
${{\bm{\phi}}^{\rm{H}}}{\bf{X}}{\bm{\phi}}=M\lambda_{\rm{\max}}$, which is a
constant. By removing the other constants, Problem (48) can be rewritten as
follows:

$\displaystyle{\mathop{\max}\limits_{\bm{\phi}}\quad 2{\mathop{\rm
Re}\nolimits}\left\\{{{{\bm{\phi}}^{\rm{H}}{\bf{q}}^{(n)}}}\right\\}}$ (49a)
$\displaystyle\textrm{s.t.}\quad(\ref{jiofjj}),(\ref{sacvfgatr}),$ (49b)

where
${\bf{q}}^{(n)}=\left({{\lambda_{\rm{\max}}}{{\bf{I}}_{M}}-{\bm{\Xi}}}\right){{\bm{\phi}}^{(n)}}-{\bf{v}}^{*}.$
Due to the additional constraint (43), the optimal solution of Problem (49)
cannot be obtained as in [24]. Furthermore, due to the non-convex unit-modulus
constraint (8d), Problem (49) is a non-convex optimization problem. As a
result, the Lagrangian dual decomposition method developed for the convex
problem (19) is not applicable here, since the dual gap is not zero.

In the following, we propose a price mechanism for solving Problem (49) that
can obtain the globally optimal solution. Specifically, we consider the
following problem by introducing a non-negative price $p$ on the left hand
side of constraint (43):

$\displaystyle{\mathop{\max}\limits_{\bm{\phi}}\quad 2{\mathop{\rm
Re}\nolimits}\left\\{{{{\bm{\phi}}^{\rm{H}}{\bf{q}}^{(n)}}}\right\\}}+2p{\mathop{\rm
Re}\nolimits}\left[{{{\bm{\phi}}^{\rm{H}}}\left({{\bf{g}^{*}}+{\bm{\Upsilon}}{{\bm{\phi}}^{(n)}}}\right)}\right]$
(50a) $\displaystyle\textrm{s.t.}\quad(\ref{jiofjj}).$ (50b)

For a given $p$, the globally optimal solution is given by

$\bm{\phi}{(p)}={e^{j\arg\left({\bf{q}}^{(n)}+p\left({{\bf{g}^{*}}+{\bm{\Upsilon}}{{\bm{\phi}}^{(n)}}}\right)\right)}}.$
(51)

Our objective is to find a $p$ value for ensuring that the complementary
slackness condition for constraint (43) is satisfied:

$p\left(J(p)-\hat{Q}\right)=0,$ (52)

where $J(p)=2{\mathop{\rm
Re}\nolimits}\left[{{{\bm{\phi}}{(p)}^{\rm{H}}}\left({{\bf{g}^{*}}+{\bm{\Upsilon}}{{\bm{\phi}}^{(n)}}}\right)}\right]$.
To solve this equation, we consider two cases: 1) $p=0$; 2) $p>0$.

_Case I:_ In this case,
$\bm{\phi}{(0)}={e^{j\arg\left({\bf{q}}^{(n)}\right)}}$ has to satisfy
constraint (43). Otherwise, $p>0$.

_Case II:_ Since $p>0$, equation (52) holds only when $J(p)=\hat{Q}$. To solve
this equation, we first provide the following lemma.

Lemma 3: Function $J(p)$ is a monotonically increasing function of $p$.

Proof: The proof is similar to Lemma 1 and thus omitted. $\Box$

Based on Lemma 3, the bisection search method can be adopted for finding the
solution of $J(p)=\hat{Q}$. Based on the above discussions, we provide the
algorithm to solve Problem (49) in Algorithm 3. Although Problem (49) is a
non-convex problem, in the following theorem we prove that Algorithm 3 is
capable of finding the globally optimal solution.

Theorem 2: Algorithm 3 is capable of finding the globally optimal solution of
Problem (49) and thus also of Problem (48).

Proof: Please refer to Appendix B. $\Box$

Algorithm 3 Bisection Search Method to Solve Problem (49)

1: Calculate $J(0)$. If $J(0)\leq\hat{Q}$, terminate. Otherwise, go to step 2.

2: Initialize the accuracy $\varepsilon$, bounds $p_{l}$ and $p_{u}$;

3: Calculate
$p={{\left({{p_{l}}+{p_{u}}}\right)}\mathord{\left/{\vphantom{{\left({{p_{l}}+{p_{u}}}\right)}2}}\right.\kern-1.2pt}2}$;

4: Update $\bm{\phi}{(p)}$ in (51) and calculate $J(p)$;

5: If $J(p)\geq\hat{Q}$, set ${p_{u}}=p$; Otherwise, set ${p_{l}}=p$;

6: If $\left|{{p_{l}}-{p_{u}}}\right|\leq\varepsilon$, terminate; Otherwise,
go to step 3.

Based on the above, we now provide the details of solving Problem (31) in
Algorithm 4.

Algorithm 4 MM Combined with SCA Algorithm to Solve Problem (31)

1: Initialize the accuracy $\varepsilon$, the phase shifts
${\bm{\phi}}^{(0)}$, the iteration index to $n=0$, the maximum number of
iterations to $n_{\rm{max}}$, calculate the OF value of Problem (44) as
$f({\bm{\phi}}^{(0)})$;

2: Calculate $\hat{Q}^{(n)}={\mathord{\buildrel{\lower
3.0pt\hbox{$\scriptscriptstyle\frown$}}\over{Q}}+{{\bm{\phi}}^{(n){\rm{H}}}}\bm{\Upsilon}{{\bm{\phi}}^{(n)}}}$;

3: Calculate
${\bf{q}}^{(n)}=\left({{\lambda_{\rm{\max}}}{{\bf{I}}_{M}}-{\bm{\Xi}}}\right){{\bm{\phi}}^{(n)}}-{\bf{v}}^{*}$;

4: Update ${{\bm{\phi}}^{(n+1)}}$ by solving Problem (49) using Algorithm 3;

5: If $n\geq n_{\rm{max}}$ or
${{\left|{f({\bm{\phi}^{(n+1)}})-f({\bm{\phi}^{(n)}})}\right|}\mathord{\left/{\vphantom{{\left|{f({\bm{\phi}^{(n+1)}})-f({\bm{\phi}^{(n)}})}\right|}{f({\bm{\phi}^{(n+1)}})}}}\right.\kern-1.2pt}{f({\bm{\phi}^{(n+1)}})}}\leq\varepsilon$
holds, terminate; Otherwise, set $n\leftarrow n+1$ and go to step 2.

In the following theorem, we prove that the sequence of
$\\{{\bm{\phi}^{(n)}},n=1,2,\cdots\\}$ generated by Algorithm 4 converges to
the KKT-optimal point of Problem (31).

Theorem 3: The sequences of the OF value produced by Algorithm 4 are
guaranteed to converge, and the final solution satisfies the KKT point of
Problem (31).

Proof: Please refer to Appendix C. $\Box$

Let us now analyze the complexity of Algorithm 4. The complexity is dominated
by calculating ${{\bm{\phi}}^{(n+1)}}$ in step 4 using Algorithm 3. The
complexity mainly depends on calculating the maximum eigenvalue of
${\bm{\Xi}}$. Its complexity is on the order of ${\cal O}(M^{3})$. The number
of iterations required for Algorithm 3 is
${\log_{2}}\left({\frac{{{p_{u}}-{p_{l}}}}{\varepsilon}}\right)$. Then, the
total complexity of step 3 is ${\cal
O}({\log_{2}}\left({\frac{{{p_{u}}-{p_{l}}}}{\varepsilon}}\right)M^{3})$.
Hence, the total complexity of Algorithm 4 is given by ${\cal
O}(n_{\rm{max}}{\log_{2}}\left({\frac{{{p_{u}}-{p_{l}}}}{\varepsilon}}\right)M^{3})$.

###  III-D Overall Algorithm to Solve Problem (8)

Based on the above analysis, we provide the detailed steps of the BCD
algorithm to solve Problem (8) in Algorithm 5, where
$R({\bf{F}}^{(n)},{\bm{\phi}}^{(n)})$ denotes the OF value of Problem (8) in
the $n$th iteration.

Algorithm 5 Block Coordinate Descent Algorithm

1: Initialize iterative number $n=1$, maximum number of iterations
$n_{\rm{max}}$, feasible ${\bf{F}}^{(1)}$, ${\bm{\phi}}^{(1)}$, error
tolerance $\varepsilon$, calculate $R({\bf{F}}^{(1)},{\bm{\phi}}^{(1)})$,
calculate the optimal decoding matrices ${\bf{U}}^{(1)}$ and auxiliary
matrices ${\bf{W}}^{(1)}$ based on (14);

2: Given ${\bf{U}}^{(n)}$, ${\bf{W}}^{(n)}$ and ${\bm{\phi}}^{(n)}$, calculate
the optimal precoding matrices ${\bf{F}}^{(n+1)}$ by solving Problem (16)
using Algorithm 2;

3: Given ${\bf{U}}^{(n)}$, ${\bf{W}}^{(n)}$ and ${\bf{F}}^{(n+1)}$, calculate
the optimal ${\bm{\phi}}^{(n+1)}$ by solving Problem (31) using Algorithm 4;

4: Given ${\bf{F}}^{(n+1)}$ and ${\bm{\phi}}^{(n+1)}$, calculate the optimal
decoding matrices ${\bf{U}}^{(n+1)}$ in (14);

5: Given ${\bf{F}}^{(n+1)}$, ${\bf{U}}^{(n+1)}$ and ${\bm{\phi}}^{(n+1)}$,
calculate the optimal auxiliary matrices ${\bf{W}}^{(n+1)}$ in (14);

6: If $n\geq n_{\rm{max}}$ or
${{\left|{R({\bf{F}}^{(n+1)},{\bm{\phi}}^{(n+1)})-R({\bf{F}}^{(n)},{\bm{\phi}}^{(n)})}\right|}\mathord{\left/{\vphantom{{\left|{{\rm{Obj}}({{\bf{F}}^{(n+1)}},{{\bm{\theta}}^{(n+1)}})-R({\bf{F}}^{(n)},{\bm{\phi}}^{(n)})}\right|}{R({\bf{F}}^{(n)},{\bm{\phi}}^{(n)})}}}\right.\kern-1.2pt}{R({\bf{F}}^{(n+1)},{\bm{\phi}}^{(n+1)})}}<\varepsilon$,
terminate. Otherwise, set $n\leftarrow n+1$ and go to step 2.

The following theorem shows the convergence and solution properties of
Algorithm 5.

Theorem 4: The OF value sequence
$\\{R({{\bf{F}}^{(n)}},{\bm{\phi}}^{(n)}{\rm{)}},n=1,2,\cdots\\}$ generated by
Algorithm 5 is guaranteed to converge, and the final solution satisfies the
KKT conditions of Problem (8).

Proof: Please refer to Appendix D. $\Box$

The complexity of Algorithm 5 mainly depends on that of Step 2 and Step 3, the
complexity of which has been analyzed in the above subsections. In specific,
the total complexity of step 2 and step 3 are respectively given by ${\cal
O}(n_{1}^{\rm{max}}{\log_{2}}\left({\frac{{{\lambda_{u}}-{\lambda_{l}}}}{\varepsilon}}\right)K_{I}N_{B}^{3})$
and ${\cal
O}(n_{2}^{\rm{max}}{\log_{2}}\left({\frac{{{p_{u}}-{p_{l}}}}{\varepsilon}}\right)M^{3})$,
where $n_{1}^{\rm{max}}$ and $n_{2}^{\rm{max}}$ denote the number of
iterations for Algorithm 2 and Algorithm 4 to converge. Denote the total
number of iterations of Algorithm 5 as $N_{\rm{max}}$. Then, the overall
complexity of Algorithm 5 is given by
$O\left({{N_{\max}}\left({n_{1}^{{\rm{max}}}{{\log}_{2}}\left({\frac{{{\lambda_{u}}-{\lambda_{l}}}}{\varepsilon}}\right){K_{I}}N_{B}^{3}+n_{2}^{{\rm{max}}}{{\log}_{2}}\left({\frac{{{p_{u}}-{p_{l}}}}{\varepsilon}}\right){M^{3}}}\right)}\right)$.
Additionally, the simulation results show that Algorithm 5 converges rapidly,
which demonstrates the low complexity of this algorithm.

##  IV Feasibility Check for Problem (8)

Due to the conflicting EH and limited transmit power constraints, Problem (8)
may be infeasible. Hence, we have to first check whether Problem (8) is
feasible or not. To this end, we construct the following optimization problem:

$\displaystyle\mathop{\max}\limits_{{{\bf{F}},{\bm{\Phi}}}}\
{\rm{tr}}\left({\sum\limits_{k=1}^{K_{I}}{{\bf{F}}_{k}^{\rm{H}}{\bf{G}}{{\bf{F}}_{k}}}}\right)$
(53a) $\displaystyle\textrm{s.t.}\quad(\ref{dehwifr}),(\ref{jiofjj}).$ (53b)

If the optimal OF value is larger than $\bar{Q}$, Problem (8) is feasible.
Otherwise, it is infeasible. As the TPC matrices and phase shift matrix are
coupled, the globally optimal solution is difficult to obtain. In the
following, we can obtain a suboptimal solution by alternatively optimizing the
TPC matrices and phase shifts.

For a given phase shift matrix, the TPC matrix optimization problem is given
by

$\displaystyle\mathop{\max}\limits_{{{\bf{F}}}}\
{\rm{tr}}\left({\sum\limits_{k=1}^{K_{I}}{{\bf{F}}_{k}^{\rm{H}}{\bf{G}}{{\bf{F}}_{k}}}}\right)$
(54a) $\displaystyle\textrm{s.t.}\quad(\ref{dehwifr}).$ (54b)

Upon denoting the maximum eigenvalue and the corresponding eigenvector of
${\bf{G}}$ by $\chi$ and $\bf{b}$ respectively, the optimal solution can be
readily obtained as
${{\bf{F}}_{k}}=\left[{\sqrt{{p_{k}}}{\bf{b}},{{\bf{0}}_{{N_{B}}\times(d-1)}}}\right],\forall
k=1,\cdots,K_{I}$, where $\sum\nolimits_{k=1}^{{K_{I}}}{{p_{k}}}={P_{T}}$ and
$p_{k}\geq 0,\forall k=1,\cdots,K_{I}$. Without loss of generality, we can set
$p_{i}=P_{T}/K_{I},\forall i\in{\cal K}_{I}$. The OF value is given by $\chi
P_{T}$. In this case, the optimal TPC matrix represents the optimal energy
beamforming, which is the same as that for the single-antenna IR case of [38].

For a given TPC matrix $\bf{F}$, the phase shift optimization problem is
formulated as:

$\displaystyle\mathop{\max}\limits_{{{\bm{\phi}}}}\
{{\bm{\phi}}^{\rm{H}}}\bm{\Upsilon}{\bm{\phi}}+2{\mathop{\rm
Re}\nolimits}\left\\{{{{\bm{\phi}}^{\rm{H}}}{\bf{g}^{*}}}\right\\}$ (55a)
$\displaystyle\textrm{s.t.}\quad(\ref{jiofjj}),$ (55b)

where $\bm{\Upsilon}$ and $\bf{g}$ are defined in the above section. The OF is
convex w.r.t. ${\bm{\phi}}$, and maximizing a convex function is a d.c
program. Hence, it can be solved by using the SCA method by approximating
${{\bm{\phi}}^{\rm{H}}}\bm{\Upsilon}{\bm{\phi}}$ as its first-order Taylor
expansion, details of which are omitted.

Finally, alternatively solve Problem (54) and (55) until the OF is larger than
$\bar{Q}$.

##  V Simulation Results

In this section, we provide simulation results for demonstrating the benefits
of applying IRS to SWIPT systems, as seen in Fig. 2, where there are four ERs
and two IRs. The ERs and IRs are uniformly and randomly scattered in a circle
centered at $(x_{\rm{ER}},0)$ and $(x_{\rm{IR}},0)$ with radius 1 m and 4 m,
respectively. The IRS is located at $(x_{\rm{IRS}},2)$. In the simulations, we
assume that the IRS is just above the ERs and thus we set
$x_{\rm{ER}}=x_{\rm{IRS}}$. The large-scale path loss is modeled in dB as

$PL=P{L_{0}}{\left({\frac{D}{{{D_{0}}}}}\right)^{-\alpha}},$ (56)

where ${\rm{P}}{{\rm{L}}_{0}}$ is the path loss at the reference distance
$D_{0}$, $D$ is the link length in meters, and $\alpha$ is the path loss
exponent. Here, we set $D_{0}=1$ and ${\rm{P}}{{\rm{L}}_{0}}=-30\rm{dB}$. The
path loss exponents of the BS-IRS, IRS-ER, IRS-IR, BS-IR and BS-ER links are
respectively set as $\alpha_{\rm{BSIRS}}=2.2$, $\alpha_{\rm{IRSER}}=2.2$,
$\alpha_{\rm{IRSIR}}=2.4$, $\alpha_{\rm{BSIR}}=3.6$ and
$\alpha_{\rm{BSER}}=3.6$. Unless otherwise stated, the other parameters are
set as follows: Channel bandwidth of 1 MHz, noise power density of $-160$
dBm/Hz, $N_{B}=4$, $N_{I}=N_{E}=2$, $d=2$, $\bar{Q}=2\times{10^{-4}}$ W,
$\eta=0.5$, $M=50$, $P_{T}=10\ W$, weight factors $\omega_{k}=1,\forall
k\in{\cal K}_{I},{\alpha_{l}}=1,\forall l\in{\cal K}_{E}$, $x_{\rm{ER}}=5$ m,
and $x_{\rm{IR}}=400$ m. The following results are obtained by averaging over
100 random locations and channel generations. Due to the severe blockage and
long distance, the channels from the BS and the IRS to the IRs are assumed to
be Rayleigh fading. However, as the BS, the ERs and the IRS are close to each
other, the small-scale channels are assumed to be Rician fading. In
particular, the small-scale channels from the IRS to the ERs are denoted as:

${{\bf{\tilde{G}}}_{r,l}}=\sqrt{\frac{{{\beta_{{\rm{irser}}}}}}{{{\beta_{{\rm{irser}}}}+1}}}{\bf{\tilde{G}}}_{r,l}^{{\rm{LoS}}}+\sqrt{\frac{1}{{{\beta_{{\rm{irser}}}}+1}}}{\bf{\tilde{G}}}_{r,l}^{{\rm{NLoS}}},l=1,\cdots,{K_{E}},$
(57)

where $\beta_{{\rm{irser}}}$ is the Rician factor,
${\bf{\tilde{G}}}_{r,l}^{{\rm{LoS}}}$ is the deterministic line of sight
(LoS), and ${\bf{\tilde{G}}}_{r,l}^{{\rm{NLoS}}}$ is the non-LoS (NLoS)
component that is Rayleigh fading. The LoS component
${\bf{\tilde{G}}}_{r,l}^{{\rm{LoS}}}$ can be modeled as
${\bf{\tilde{G}}}_{r,l}^{{\rm{LoS}}}={{\bf{a}}_{{N_{E}}}}\left({\vartheta_{{\rm{irser}},l}^{AoA}}\right){\bf{a}}_{M}^{H}\left({\vartheta_{{\rm{irser}},l}^{AoD}}\right)$,
where ${{\bf{a}}_{{N_{E}}}}\left({\vartheta_{{\rm{irser}},l}^{AoA}}\right)$ is
defined as

${{\bf{a}}_{{N_{E}}}}\left({\vartheta_{{\rm{irser}},l}^{AoA}}\right)={\left[{1,{e^{j\frac{{2\pi
d}}{\lambda}\sin\vartheta_{{\rm{irser}},l}^{AoA}}},\cdots,{e^{j\frac{{2\pi
d}}{\lambda}({N_{E}}-1)\sin\vartheta_{{\rm{irser}},l}^{AoA}}}}\right]^{T}}$
(58)

and

${{\bf{a}}_{M}}\left({\vartheta_{{\rm{irser}},l}^{AoD}}\right)={\left[{1,{e^{j\frac{{2\pi
d}}{\lambda}\sin\vartheta_{{\rm{irser}},l}^{AoD}}},\cdots,{e^{j\frac{{2\pi
d}}{\lambda}(M-1)\sin\vartheta_{{\rm{irser}},l}^{AoD}}}}\right]^{T}}.$ (59)

In (58) and (59), $d$ is the antenna separation distance, $\lambda$ is the
wavelength, $\vartheta_{{\rm{irser}},l}^{AoD}$ is the angle of departure and
$\vartheta_{{\rm{irser}},l}^{AoA}$ is the angle of arrival. It is assumed that
$\vartheta_{{\rm{irser}},l}^{AoD}$ and $\vartheta_{{\rm{irser}},l}^{AoA}$ are
randomly distributed within $[0,2\pi]$. For simplicity, we set
$d/\lambda=1/2$. The small-scale channels from the BS to the ERs and the IRS
are similarly defined. For simplicity, the Rician factors for all Rician
fading channels are assumed to be the same as $\beta=3$.

Figure 2: The simulated IRS-aided SWIPT MIMO communication scenario.

Figure 3: Maximum harvested power achieved by various schemes.

Figure 4: Convergence behaviour of the BCD algorithm.

We first compare the maximum power harvested by various schemes in Fig. 4.
Specifically, we solve the EH maximization problem (53) by using the
feasibility check method in Section IV. Additionally, we also present the
results without using IRS. Fig. 4 shows the maximum EH power versus the ER
circle center location $x_{\rm{EH}}$. As expected, the EH power gleaned by all
schemes decreases, when the ERs are far away from the BS. As expected, more
power can be harvested with the aid of IRS than that without IRS, especially
when the number of phase shifters $M$ is large. This is mainly due to the fact
that an additional strong link is reflected by the IRS, which can be harvested
by the ERs. This figure also shows that the IRS is effective in expanding the
operational range of ERs. For example, when the harvested power limit is
$\bar{Q}=2\times 10^{-4}$ W, the maximum operational range of the system
without IRS is only 5.5 m, while the system having $M=40$ phase shifters can
operate for distances up to 9 m.

In Fig. 4, we study the convergence behaviour of the BCD algorithm for
different numbers of phase shifters $M$. It is observed from Fig. 4 that the
WSR achieved for various $M$ values increases monotonically with the number of
iterations, which verifies Theorem 4. Additionally, the BCD algorithm
converges rapidly and in general a few iterations are sufficient for the BCD
algorithm to achieve a large portion of the converged WSR. This reflects the
low complexity of the BCD algorithm, which is appealing for practical
applications.

Figure 5:  WSR versus the number of phase shifters.

Figure 6: WSR versus the harvested power requirement $\bar{Q}$.

In the following, we compare our proposed BCD algorithms to a pair of
benchmark schemes: 1)‘No-IRS’: In this scheme, there is no IRS to assist the
transmission as in conventional systems; 2) ‘Fixed Phase’: In this method, the
phase shifts are fixed at the solutions obtained by solving the harvested
power maximization problem (53), while they are not optimized, when using the
BCD algorithm by removing Step 3 of the BCD algorithm. When any of the methods
fails to obtain a feasible solution, its achievable WSR is set to zero.

In Fig. 6, we first study the impact of the number of phase shifters $M$ on
the performance of various algorithms. As expected, the WSR achieved by all
the algorithms - except for the No-IRS method - increases with $M$, since a
higher degree of freedom can be exploited for optimizing the system
performance. By carefully optimizing the phase shifts for maximizing the WSR,
the BCD algorithm significantly outperforms the fixed-phase scheme.
Additionally, the performance gain increases with $M$, which emphasizes the
importance of optimizing the phase shifts. By employing the IRS in our SWIPT
system, the WSR obtained by the BCD algorithm becomes drastically higher than
that of No-IRS. For example, when $M=60$, the WSR performance gain is up to 10
bit/s/Hz. These results demonstrate that introducing the IRS into our SWIPT
system is a promising technique of enhancing the system performance.

In Fig. 6, the impact of harvested power requirement $\bar{Q}$ is
investigated. It is seen from this figure that the WSR achieved by all the
algorithms decreases upon increasing $\bar{Q}$, because the probability of
infeasibility increases, which in turn reduces the average WSR value. We also
find that the WSR obtained by the No-IRS scheme decreases more rapidly than
that of the other two IRS-aided transmission schemes. The WSR of the No-IRS is
approaching zero when $\bar{Q}=4\times 10^{-4}$ W, while those relying on IRSs
achieve a WSR gain in excess of 20 bit/s/Hz. It is observed again that the BCD
algorithm performs better than the fixed-phase scheme, but the gap narrows
with the increase of $\bar{Q}$. This can be explained as follows. With the
increase of $\bar{Q}$, both the TPC matrices and the phase shifts should be
designed for maximizing the power harvested at the ERs, and thus the final
solutions of the fixed-phase and BCD method will become the same.

Figure 7:  WSR versus the IRS-related path loss exponent $\alpha_{\rm{IRS}}$.

Figure 8: WSR versus the location of ER circle center $x_{\rm{ER}}$.

The above results are obtained for $\alpha_{\rm{BSIRS}}=2.2$,
$\alpha_{\rm{IRSER}}=2.2$, $\alpha_{\rm{IRSIR}}=2.4$ based on the assumption
that the IRS relies on an obstacle-free scenario. In practice, this ideal
scenario is seldom encountered. Hence, it is imperative to investigate the
impact of
$\alpha_{\rm{IRS}}\buildrel\Delta\over{=}\alpha_{\rm{BSIRS}}=\alpha_{\rm{IRSER}}=\alpha_{\rm{IRSIR}}$
on the system performance, which is shown in Fig. 8. Observe from this figure
that the WSR achieved by the algorithms using IRS decreases drastically with
$\alpha_{\rm{IRS}}$. When $\alpha_{\rm{IRS}}=3$, the WSR-performance gain of
our algorithm over the No-IRS scenario is only 7 bit/s/Hz, because upon
increasing $\alpha_{\rm{IRS}}$, the signal power reflected from the IRS
becomes weaker. Hence, the benefits of the IRS can be eroded. This provides an
important engineering design insight: the location of IRS should be carefully
considered for finding an obstacle-free scenario associated with a low
$\alpha_{\rm{IRS}}$.

In Fig. 8, we study the impact of ER locations on the system performance. As
expected, the WSR achieved by all the schemes decreases with $x_{\rm{IRS}}$,
since the ERs become more distant from the BS and the signals gleaned from
both the BS and IRS become weaker. The WSR achieved by the No-IRS approaches
zero when $x_{\rm{IRS}}=8$ m, hence this method cannot reach the energy
transmission target of the ERs. The proposed algorithm is again observed to
significantly outperform the other two algorithms, especially when the ERs are
close to the BS.

Figure 9:  WSR versus the location of IR circle center $x_{\rm{IR}}$.

Finally, the impact of IR locations is investigated in Fig. 9. It is observed
that the WSR achieved by all the algorithms decreases with $x_{\rm{IR}}$ since
the IRs become farther away from the BS when increasing $x_{\rm{IR}}$. The
proposed algorithm is shown to achieve nearly the WSR gain of 10 bit/s/Hz over
the No-IRS when $x_{\rm{IR}}=100\ m$, and the WSR gain slightly increases with
$x_{\rm{IR}}$. This means that the IRS is more advantageous when the IRs are
far away from the BS, and the IRS can provide one additional favorable link.

##  VI Conclusions

In this paper, we have invoked an IRS in a SWIPT MIMO system for enhancing the
performance of both the ERs and IRs. By carefully adjusting the phase shifts
at the IRS, the signal reflected by the IRS can be added constructively at
both the ERs and IRs. We considered the WSR maximization problem of IRs, while
guaranteeing the energy harvesting requirements of the ERs and the associated
non-convex unit-modulus constraints. We conceived a BCD algorithm for
alternatively optimizing the TPC matrices at the BS and the phase shift matrix
at the IRSs. For each subproblem, a low-complexity iterative algorithm was
proposed, which guarantees to be at worst locally optimal. Our simulation
results demonstrated that the IRS enhances the performance of the SWIPT system
and that the proposed algorithm converges rapidly, hence it is eminently
suitable for practical implementations.

This paper assumes perfect CSI at the BS, which is challenging to obtain. For
the future work, we will consider the robust transmission design for the IRS-
aided SWIPT system, where the CSI is assumed to be imperfectly known. In
addition, how to design the discrete phase shifts will be left for future
work.

##  Appendix A Proof of Lemma 1

We consider a pair of variables $\lambda$ and $\lambda^{\prime}$, where
$\lambda>\lambda^{\prime}$. Let ${\bf{F}(\lambda)}$ and
${\bf{F}(\lambda^{\prime})}$ be the optimal solutions of Problem (21) with
$\lambda$ and $\lambda^{\prime}$, respectively. Since ${\bf{F}(\lambda)}$ is
the optimal solution of Problem (21) with $\lambda$, we have

${\cal L}[{\bf{F}(\lambda)},\lambda]\leq{\cal
L}[{\bf{F}(\lambda^{\prime})},\lambda].$ (A.1)

Similarly, we have

${\cal L}[{\bf{F}(\lambda^{\prime})},\lambda^{\prime}]\leq{\cal
L}[{\bf{F}(\lambda)},\lambda^{\prime}].$ (A.2)

By adding these two inequalities and simplifying them, we have
$\left({\lambda-\lambda^{\prime}}\right)P(\lambda)\leq\left({\lambda-\lambda^{\prime}}\right)P(\lambda^{\prime})$.
Since $\lambda>\lambda^{\prime}$, we have $P(\lambda)\leq
P(\lambda^{\prime})$, which completes the proof.

##  Appendix B Proof of Theorem 2

Denote the globally optimal solution of Problem (49) by $\bm{\phi}^{\star}$.
According to [43], for a non-convex optimization problem, all its locally
optimal solutions (including the globally optimal solution) should satisfy the
Karush-Kuhn-Tucker (KKT) optimality conditions, one of which is the
complementary slackness condition for constraint (43):

$\lambda^{\star}\left(2{\mathop{\rm
Re}\nolimits}\left[{{{\bm{\phi}}^{\star\rm{H}}}\left({{\bf{g}^{*}}+{\bm{\Upsilon}}{{\bm{\phi}}^{(n)}}}\right)}\right]-\hat{Q}\right)=0,\vspace{-0.2cm}$
(B.1)

where $\lambda^{\star}$ is the corresponding optimal dual variable. We
consider two cases: 1) $\lambda^{\star}=0$; 2) $\lambda^{\star}>0$.

The first case means that constraint (43) is not tight in the optimum. Then,
the optimal solution can be obtained as
$\bm{\phi}^{\star}={e^{j\arg\left({\bf{q}}^{(n)}\right)}}$, which is equal to
$\bm{\phi}{(0)}$. Hence, Algorithm 3 achieves the optimal solution of Problem
(49).

For the second case, the following equality should hold:

$2{\mathop{\rm
Re}\nolimits}\left[{{{\bm{\phi}}^{\star\rm{H}}}\left({{\bf{g}^{*}}+{\bm{\Upsilon}}{{\bm{\phi}}^{(n)}}}\right)}\right]=\hat{Q}.$
(B.2)

We prove the second case by using the method of contradiction. Denote the
optimal $p$ obtained by Algorithm 3 as $p^{\star}$, and the corresponding
$\bm{\phi}$ as $\bm{\phi}(p^{\star})$. Then, we have

$2{\mathop{\rm
Re}\nolimits}\left[{{\bm{\phi}(p^{\star})^{\rm{H}}}\left({{\bf{g}^{*}}+{\bm{\Upsilon}}{{\bm{\phi}}^{(n)}}}\right)}\right]=\hat{Q}.$
(B.3)

Let us assume that $\bm{\phi}(p^{\star})$ is not the globally optimal solution
of Problem (49). Then, we have

$2{\mathop{\rm
Re}\nolimits}\left\\{{{\bm{\phi}(p^{\star})^{\rm{H}}{\bf{q}}^{(n)}}}\right\\}<2{\mathop{\rm
Re}\nolimits}\left\\{{{\bm{\phi}^{\star\rm{H}}{\bf{q}}^{(n)}}}\right\\}.$
(B.4)

Since $\bm{\phi}(p^{\star})$ is the globally optimal solution of Problem (50)
when $p=p^{\star}$, we have

$\\!\\!\\!2{\mathop{\rm
Re}\nolimits}\left\\{{{{\bm{\phi}}(p^{\star})^{\rm{H}}{\bf{q}}^{(n)}}}\right\\}\\!+\\!2p^{\star}{\mathop{\rm
Re}\nolimits}\left[{{{\bm{\phi}}(p^{\star})^{\rm{H}}}\left({{\bf{g}^{*}}\\!+\\!{\bm{\Upsilon}}{{\bm{\phi}}^{(n)}}}\right)}\right]\\!\geq\\!2{\mathop{\rm
Re}\nolimits}\left\\{{{{\bm{\phi}}^{\star\rm{H}}{\bf{q}}^{(n)}}}\right\\}\\!+\\!2p^{\star}{\mathop{\rm
Re}\nolimits}\left[{{{\bm{\phi}}^{\star\rm{H}}}\left({{\bf{g}^{*}}\\!+\\!{\bm{\Upsilon}}{{\bm{\phi}}^{(n)}}}\right)}\right].\vspace{-0.2cm}$
(B.5)

By substituting (B.2) and (B.3) into (B.5), we have

$2{\mathop{\rm
Re}\nolimits}\left\\{{{{\bm{\phi}}(p^{\star})^{\rm{H}}{\bf{q}}^{(n)}}}\right\\}\geq
2{\mathop{\rm
Re}\nolimits}\left\\{{{{\bm{\phi}}^{\star\rm{H}}{\bf{q}}^{(n)}}}\right\\},\vspace{-0.2cm}$
(B.6)

which contradicts (B.4). Hence, the solution obtained by Algorithm 3 is the
globally optimal solution of Problem (49). Since Problem (48) is equivalent to
Problem (49), the proof is complete.

##  Appendix C Proof of Theorem 3

Let us define the following functions:

$\displaystyle T({\bm{\phi}})$ $\displaystyle\buildrel\Delta\over{=}$
$\displaystyle{{\bm{\phi}}^{\rm{H}}}\bm{\Upsilon}{\bm{\phi}}+2{\mathop{\rm
Re}\nolimits}\left\\{{{{\bm{\phi}}^{\rm{H}}}{\bf{g}^{*}}}\right\\}+{\rm{tr}}\left({{{\bf{G}}_{b}}{\bf{\tilde{F}}}}\right),$
(C.1) $\displaystyle\bar{T}({\bm{\phi}}|{\bm{\phi}}^{(n)})$
$\displaystyle\buildrel\Delta\over{=}$
$\displaystyle-{{\bm{\phi}^{(n)}}^{\rm{H}}}\bm{\Upsilon}{\bm{\phi}^{(n)}}+2{\mathop{\rm
Re}\nolimits}\left[{{{\bm{\phi}}^{\rm{H}}}\left({{\bf{g}^{*}}+{\bm{\Upsilon}}{{\bm{\phi}}^{(n)}}}\right)}\right]+{\rm{tr}}\left({{{\bf{G}}_{b}}{\bf{\tilde{F}}}}\right).$
(C.2)

It can be verified that
$T({\bm{\phi}}^{(n)})=\bar{T}({\bm{\phi}^{(n)}}|{\bm{\phi}}^{(n)})$.

We first show that the solution sequence
$\\{{\bm{\phi}}^{(n)},n=1,2,\cdots\\}$ is feasible for Problem (31). The unit-
modulus constraint is guaranteed in (51). We only have to check the EH
constraint in (8c). Note that ${\bm{\phi}}^{(n+1)}$ is a feasible solution of
Problem (49), and thus satisfies constraint (42). Hence, we have
$\bar{T}({\bm{\phi}}^{(n+1)}|{\bm{\phi}}^{(n)})\geq\bar{Q}$. By using
inequality (42), we have
$T({\bm{\phi}}^{(n+1)})\geq\bar{T}({\bm{\phi}}^{(n+1)}|{\bm{\phi}}^{(n)})$.
Then, $T({\bm{\phi}}^{(n+1)})\geq\bar{Q}$ holds, which means that the sequence
of ${\bm{\phi}}^{(n+1)}$ satisfies the EH constraint in (8c).

Now, we show that the OF value sequence
$\\{f({\bm{\phi}}^{(n)}),n=1,2,\cdots\\}$ is monotonically decreasing. Based
on Theorem 2, the globally optimal solution $\bm{\Phi}$ to Problem (48) can be
obtained. Then, we have $g({\bm{\phi}}^{(n+1)}|{\bm{\phi}}^{(n)})\leq
g({\bm{\phi}^{(n)}}|{\bm{\phi}}^{(n)})$. According to the first condition in
(45), we have
$g({\bm{\phi}}^{(n)}|{\bm{\phi}}^{(n)})\\!=\\!f({\bm{\phi}}^{(n)})$. Hence, we
have $g({\bm{\phi}}^{(n+1)}|{\bm{\phi}}^{(n)})\leq f({\bm{\phi}}^{(n)})$. By
using the third condition of (45), we have
$g({\bm{\phi}}^{(n+1)}|{\bm{\phi}}^{(n)})\geq f({\bm{\phi}}^{(n+1)})$. As a
result, we have $f({\bm{\phi}}^{(n)})\geq f({\bm{\phi}}^{(n+1)})$.
Additionally, the OF must have a lower bound due to the unit-modulus
constraint. Hence, the OF value sequence
$\\{f({\bm{\phi}}^{(n)}),n=1,2,\cdots\\}$ is guaranteed to converge.

Now, we prove that the converged solution satisfies the KKT conditions of
Problem (31). Let us denote the converged solution by
$\\{{\bm{\phi}}^{\star}\\}$. Since ${\bm{\phi}}^{\star}$ is the globally
optimal solution of Problem (48), it must satisfy the KKT conditions of
Problem (48). Specifically, the Lagrange function of Problem (48) is given by

${\cal
L}({\bm{\phi}},\nu,\bm{\tau})=g({\bm{\phi}}|{\bm{\phi}}^{\star})+\nu\left(\hat{Q}-2{\mathop{\rm
Re}\nolimits}\left[{{{\bm{\phi}}^{\rm{H}}}\left({{\bf{g}^{*}}+{\bm{\Upsilon}}{{\bm{\phi}}^{\star}}}\right)}\right]\right)+\sum\limits_{m=1}^{M}{{\tau_{m}}\left({\left|{{\phi_{m}}}\right|-1}\right)},\vspace{-0.2cm}$
(C.3)

where $\nu$ and ${\bm{\tau}}=\\{\tau_{1},\cdots,\tau_{M}\\}$ are the
corresponding dual variables. Then, there must exist a $\nu^{\star}$ and
${\bm{\tau}}^{\star}=\\{\tau_{1}^{\star},\cdots,\tau_{M}^{\star}\\}$ for
ensuring that the following conditions are satisfied:

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!{{{\nabla_{{\bm{\phi}^{*}}}}{\cal
L}({\bm{\phi}},\nu,\bm{\tau})|_{{\bm{\phi}}={\bm{\phi}}^{\star}}}}\\!\\!=\\!\\!{{{\nabla_{\bm{\phi}^{*}}}g({\bm{\phi}}|{{\bm{\phi}}^{\star}})|_{\bm{\phi}=\bm{\phi}^{\star}}}}\\!\\!-\\!\\!\nu^{\star}\left({{\bf{g}^{*}}\\!+\\!{\bm{\Upsilon}}{{\bm{\phi}}^{\star}}}\right)\\!\\!+\\!\\!{\sum\limits_{m=1}^{M}\\!\\!{\tau_{m}^{\star}\\!\left({{\nabla_{\bm{\phi}^{*}}}\left|{{\phi_{m}}}\right|}\right)|_{\bm{\phi}=\bm{\phi}^{\star}}}}\\!\\!\\!$
$\displaystyle=$ $\displaystyle\\!\\!\\!{\bf{0}},$ (C.4)
$\displaystyle\\!\\!\\!\\!\nu^{\star}\left(\\!\hat{Q}\\!-\\!2{\mathop{\rm
Re}\nolimits}\left[\\!{{{\bm{\phi}^{\star}}^{\rm{H}}}\left({{\bf{g}^{*}}\\!+\\!{\bm{\Upsilon}}{{\bm{\phi}}^{\star}}}\right)}\\!\right]\\!\right)$
$\displaystyle\\!\\!\\!\\!\\!=\\!\\!\\!\\!\\!$ $\displaystyle 0,$ (C.5)
$\displaystyle{\tau_{m}^{\star}}\left({\left|{{\phi_{m}^{\star}}}\right|-1}\right)$
$\displaystyle\\!\\!\\!\\!\\!=\\!\\!\\!\\!\\!$ $\displaystyle 0,\forall
m.\vspace{-0.2cm}$ (C.6)

According to the second condition of (45), we have

${\nabla_{\bm{\phi}^{*}}}g({\bm{\phi}}|{{\bm{\phi}}^{\star}})|_{{\bm{\phi}}={\bm{\phi}}^{\star}}=\nabla_{\bm{\phi}^{*}}f({\bm{\phi}})|_{{\bm{\phi}}={\bm{\phi}}^{\star}}.\vspace{-0.2cm}$
(C.7)

Upon denoting the OF of Problem (31) as $\varphi({\bm{\phi}})$, which is the
same as $f({\bm{\phi}})$ except that $\varphi({\bm{\phi}})$ has more
constants, we have
$\nabla_{\bm{\phi}^{*}}f({\bm{\phi}})|_{{\bm{\phi}}={\bm{\phi}}^{\star}}=\nabla_{\bm{\phi}^{*}}\varphi({\bm{\phi}})|_{{\bm{\phi}}={\bm{\phi}}^{\star}}$.
Combining with (C.7), we have
${\nabla_{\bm{\phi}^{*}}}g({\bm{\phi}}|{{\bm{\phi}}^{\star}})|_{{\bm{\phi}}={\bm{\phi}}^{\star}}=\nabla_{\bm{\phi}^{*}}\varphi({\bm{\phi}})|_{{\bm{\phi}}={\bm{\phi}}^{\star}}$.
By substituting it into (C.4), we arrive at

$\nabla_{\bm{\phi}^{*}}\varphi({\bm{\phi}})|_{{\bm{\phi}}={\bm{\phi}}^{\star}}-\nu^{\star}\left({{\bf{g}^{*}}+{\bm{\Upsilon}}{{\bm{\phi}}^{\star}}}\right)+{\sum\limits_{m=1}^{M}{\tau_{m}^{\star}\left({{\nabla_{\bm{\phi}^{*}}}\left|{{\phi_{m}}}\right|}\right)|_{\bm{\phi}=\bm{\phi}^{\star}}}}={\bf{0}}.\vspace{-0.2cm}$
(C.8)

It can be checked that the set of equations (C.5), (C.6) and (C.8) constitutes
exactly the KKT conditions of Problem (31). Hence, the proof is complete.

##  Appendix D Proof of Theorem 4

Let us define the OF of Problem (12) as

$h\left({{\bf{W}},{\bf{U}},{\bf{F}},{\bm{\Phi}}}\right)\buildrel\Delta\over{=}\sum\limits_{k=1}^{K_{I}}{{\omega_{k}}{h_{k}}\left({{\bf{W}},{\bf{U}},{\bf{F}},{\bm{\Phi}}}\right)}.\vspace{-0.2cm}$
(D.1)

It can be readily verified that the sequence of solutions
$\\{{{\bf{F}}^{(n)}},{\bm{\phi}}^{(n)}\\}$ generated by Algorithm 5 is always
feasible for Problem (8). The monotonic property of Algorithm 5 can be
similarly proved by using the method of [33].

In the following, we prove that the converged solution satisfies the KKT
conditions of Problem (8). Let us denote the converged solution as
$\\{{\bf{W}}^{\star},{\bf{U}}^{\star},{\bf{F}}^{\star},{\bm{\Phi}}^{\star}\\}$.

According to Theorem 1, ${\bf{F}}^{\star}$ is the KKT-optimum point of Problem
(16). Upon denoting the OF of Problem (16) as
$z({\bf{F}},{\bm{\Phi}^{\star}})$, the Lagrange function of Problem (16) is
given by

${\cal
L}({\bf{F}},\lambda,\mu)=z({\bf{F}},{\bm{\Phi}^{\star}})+\lambda\left(\sum\limits_{k=1}^{K_{I}}{\left\|{{{\bf{F}}_{k}}}\right\|_{F}^{2}}-{P_{T}}\right)+\mu\left(\bar{Q}-{{\rm{tr}}\left({\sum\limits_{k=1}^{{K_{I}}}{{\bf{F}}_{k}^{{\rm{H}}}{\bf{G}}{{\bf{F}}_{k}}}}\right)}\right),$
(D.2)

where $\lambda$ and $\mu$ are the corresponding dual variables. Then, there
must exist a $\lambda^{\star}$ and $\mu^{\star}$ for ensuring that the
following conditions are satisfied 333For simplicity, the prime constraints
are omitted. :

$\displaystyle{{\left.{{\nabla_{{\bf{F}}_{k}^{*}}}\cal
L({\bf{F}},\lambda,\mu)}\right|_{{{\bf{F}}_{k}}={\bf{F}}_{k}^{\star}}}}\\!=\\!{{\left.{{\nabla_{{\bf{F}}_{k}^{*}}}z({\bf{F}},{\bm{\Phi}^{\star}})}\right|_{{{\bf{F}}_{k}}={\bf{F}}_{k}^{\star}}}}+\lambda^{\star}{\bf{F}}_{k}^{\star}-\mu^{\star}{\bf{GF}}_{k}^{\star}$
$\displaystyle\\!=\\!$ $\displaystyle{\bf{0}},\forall k\in{\cal K}_{I},$ (D.3)
$\displaystyle\lambda^{\star}\left(\sum\limits_{k=1}^{K_{I}}{\left\|{{{\bf{F}}_{k}^{\star}}}\right\|_{F}^{2}}-{P_{T}}\right)$
$\displaystyle\\!\\!=\\!\\!$ $\displaystyle 0,$ (D.4)
$\displaystyle\mu^{\star}\left(\bar{Q}-{{\rm{tr}}\left({\sum\limits_{k=1}^{{K_{I}}}{{\bf{F}}_{k}^{\star{\rm{H}}}{\bf{G}}{{\bf{F}}_{k}^{\star}}}}\right)}\right)$
$\displaystyle\\!\\!=\\!\\!$ $\displaystyle 0.$ (D.5)

Furthermore, it can be readily checked that

${\left.{{\nabla_{{{\bf{F}}_{k}^{*}}}}h\left({{{\bf{W}}^{\star}},{{\bf{U}}^{\star}},{\bf{F}},{\bm{\Phi}^{\star}}}\right)}\right|_{{{\bf{F}}_{k}}={\bf{F}}_{k}^{\star}}}={\left.{{\nabla_{{{\bf{F}}_{k}^{*}}}}z({\bf{F}},{\bm{\Phi}^{\star}})}\right|_{{{\bf{F}}_{k}}={\bf{F}}_{k}^{\star}}},\forall
k\in{\cal K}_{I}.\vspace{-0.2cm}$ (D.6)

To expound a little further, we have the following chain of inequalities:

$\displaystyle{\left.{{\nabla_{{{\bf{F}}_{k}^{*}}}}h_{k}\left({{{\bf{W}}^{\star}},{{\bf{U}}^{\star}},{\bf{F}},{\bm{\Phi}^{\star}}}\right)}\right|_{{{\bf{F}}_{k}}={\bf{F}}_{k}^{\star}}}$
(D.7) $\displaystyle=$
$\displaystyle-{\rm{tr}}\left({{\bf{W}}_{k}^{\star}\left({{{\left.{{\nabla_{{\bf{F}}_{k}^{*}}}{{\bf{E}}_{k}}\left({{{\bf{U}}^{\star}},{\bf{F}},{\bm{\Phi}^{\star}}}\right)}\right|}_{{{\bf{F}}_{k}}={\bf{F}}_{k}^{\star}}}}\right)}\right)$
(D.8) $\displaystyle=$
$\displaystyle-{\rm{tr}}\left({{{\left({{{\bf{E}}_{k}}\left({{{\bf{U}}^{\star}},{{\bf{F}}^{\star}},{\bm{\Phi}^{\star}}}\right)}\right)}^{-1}}\left({{{\left.{{\nabla_{{\bf{F}}_{k}^{*}}}{{\bf{E}}_{k}}\left({{{\bf{U}}^{\star}},{\bf{F}},{\Phi^{\star}}}\right)}\right|}_{{{\bf{F}}_{k}}={\bf{F}}_{k}^{\star}}}}\right)}\right)$
(D.9) $\displaystyle=$
$\displaystyle{\left.{\left({{\nabla_{{\bf{F}}_{k}^{*}}}\log\left|{{{\left({{{\bf{E}}_{k}}\left({{{\bf{U}}^{\star}},{\bf{F}},{\bm{\Phi}^{\star}}}\right)}\right)}^{-1}}}\right|}\right)}\right|_{{{\bf{F}}_{k}}={\bf{F}}_{k}^{\star}}}$
(D.10) $\displaystyle=$
$\displaystyle{\left.{{\nabla_{{\bf{F}}_{k}^{*}}}{R_{k}}({\bf{F}},{\bm{\Phi}^{\star}})}\right|_{{{\bf{F}}_{k}}={\bf{F}}_{k}^{\star}}},$
(D.11)

where (D.8) follows from the chain rule, and the final equality follows from
applying the Woodbury matrix identity to (15). Combining (D.11) with (D.6), we
have

${\left.{{\nabla_{{{\bf{F}}_{k}^{*}}}}z({\bf{F}},{\bm{\Phi}^{\star}})}\right|_{{{\bf{F}}_{k}}={\bf{F}}_{k}^{\star}}}={\left.{{\nabla_{{\bf{F}}_{k}^{*}}}{R_{k}}({\bf{F}},{\bm{\Phi}^{\star}})}\right|_{{{\bf{F}}_{k}}={\bf{F}}_{k}^{\star}}}.$
(D.12)

By substituting (D.12) into (D.3), we arrive at

${\left.{{\nabla_{{\bf{F}}_{k}^{*}}}{R_{k}}({\bf{F}},{\bm{\Phi}^{\star}})}\right|_{{{\bf{F}}_{k}}={\bf{F}}_{k}^{\star}}}+\lambda^{\star}{\bf{F}}_{k}^{\star}-\mu^{\star}{\bf{GF}}_{k}^{\star}={\bf{0}},\forall
k\in{\cal K}_{I}.$ (D.13)

According to Theorem 3, ${\bm{\phi}}^{\star}$ satisfies the KKT conditions of
Problem (31), and thus the set of equations (C.5), (C.6) and (C.8) hold.

Furthermore, it can be readily verified that

${\nabla_{{\bm{\phi}^{*}}}}h\left({\bf{W}}^{\star},{\bf{U}}^{\star},{\bf{F}}^{\star},{\bm{\Phi}}\right)|_{{\bm{\phi}}={\bm{\phi}}^{\star}}=\nabla_{\bm{\phi}^{*}}\varphi({\bm{\phi}})|_{{\bm{\phi}}={\bm{\phi}}^{\star}}.$
(D.14)

By using similar derivations as in (D.7)-(D.11), we can prove that

${\nabla_{{\bm{\phi}^{*}}}}h\left({\bf{W}}^{\star},{\bf{U}}^{\star},{\bf{F}}^{\star},{\bm{\Phi}}\right)|_{{\bm{\phi}}={\bm{\phi}}^{\star}}={\left.{{\nabla_{\bm{\phi}^{*}}}{R_{k}}(\bm{\phi},{\bf{F}}^{\star})}\right|_{{\bm{\phi}}=\bm{\phi}^{\star}}}.$
(D.15)

Hence, we have

$\nabla_{\bm{\phi}^{*}}\varphi({\bm{\phi}})|_{{\bm{\phi}}={\bm{\phi}}^{\star}}={\left.{{\nabla_{\bm{\phi}^{*}}}{R_{k}}(\bm{\phi},{\bf{F}}^{\star})}\right|_{{\bm{\phi}}=\bm{\phi}^{\star}}}.$
(D.16)

By substituting (D.16) into (C.8), we arrive at:

${\left.{{\nabla_{\bm{\phi}^{*}}}{R_{k}}(\bm{\phi},{\bf{F}}^{\star})}\right|_{{\bm{\phi}}=\bm{\phi}^{\star}}}-\nu^{\star}\left({{\bf{g}^{*}}\\!+\\!{\bm{\Upsilon}}{{\bm{\phi}}^{\star}}}\right)+{\sum\limits_{m=1}^{M}{\tau_{m}^{\star}\\!\left({{\nabla_{\bm{\phi}^{*}}}\left|{{\phi_{m}}}\right|}\right)|_{\bm{\phi}=\bm{\phi}^{\star}}}}={\bf{0}}.$
(D.17)

Then, the set of equations (D.13), (D.4), (D.5), (D.17), (C.5), and (C.6)
constitute exactly the KKT conditions of Problem (8).

## References

  * [1] M. Di Renzo, M. Debbah, D.-T. Phan-Huy, A. Zappone, M.-S. Alouini, C. Yuen, V. Sciancalepore, G. C. Alexandropoulos, J. Hoydis, H. Gacanin _et al._ , “Smart radio environments empowered by reconfigurable AI meta-surfaces: an idea whose time has come,” _EURASIP J. Wireless Commun. Networking_ , vol. 2019, no. 1, p. 129, 2019.
  * [2] Q. Wu and R. Zhang, “Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network.” [Online]. Available: https://arxiv.org/abs/1905.00152
  * [3] J. Zhang, E. Björnson, M. Matthaiou, D. W. K. Ng, H. Yang, and D. Love, “Multiple antenna technologies for beyond 5G.” [Online]. Available: https://arxiv.org/abs/1910.00092
  * [4] T. J. Cui, M. Q. Qi, X. Wan, J. Zhao, and Q. Cheng, “Coding metamaterials, digital metamaterials and programmable metamaterials,” _Light: Science & Applications_, vol. 3, no. 10, p. e218, 2014.
  * [5] J. Zhang, S. Chen, Y. Lin, J. Zheng, B. Ai, and L. Hanzo, “Cell-free massive MIMO: A new next-generation paradigm,” _IEEE Access_ , vol. 7, pp. 99 878–99 888, 2019.
  * [6] J. Zhang, L. Dai, Z. He, B. Ai, and O. A. Dobre, “Mixed-ADC/DAC multipair massive MIMO relaying systems: Performance analysis and power optimization,” _IEEE Trans. Commun._ , vol. 67, no. 1, pp. 140–153, 2018\.
  * [7] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network: Joint active and passive beamforming design,” in _2018 IEEE Global Communications Conference (GLOBECOM)_ , Dec. 2018, pp. 1–6.
  * [8] X. Yu, D. Xu, and R. Schober, “MISO wireless communication systems via intelligent reflecting surfaces.” [Online]. Available: https://arxiv.org/abs/1904.12199
  * [9] Y. Yang, B. Zheng, S. Zhang, and R. Zhang, “Intelligent reflecting surface meets OFDM: Protocol design and rate maximization.” [Online]. Available: https://arxiv.org/abs/1906.09956
  * [10] Y. Han, W. Tang, S. Jin, C. Wen, and X. Ma, “Large intelligent surface-assisted wireless communication exploiting statistical CSI,” _IEEE Trans. Veh. Technol._ , 2019.
  * [11] S. Abeywickrama, R. Zhang, and C. Yuen, “Intelligent reflecting surface: Practical phase shift model and beamforming optimization.” [Online]. Available: https://arxiv.org/abs/1907.06002
  * [12] 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.
  * [13] C. Huang, A. Zappone, G. C. Alexandropoulos, M. Debbah, and C. Yuen, “Reconfigurable intelligent surfaces for energy efficiency in wireless communication,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 8, pp. 4157–4170, Aug. 2019.
  * [14] H. Guo, Y. Liang, J. Chen, and E. Larsson, “Weighted sum-rate optimization for intelligent reflecting surface enhanced wireless networks.” [Online]. Available: https://arxiv.org/abs/1905.07920
  * [15] Q. Nadeem, A. Kammoun, A. Chaaban, M. Debbah, and M. Alouini, “Large intelligent surface assisted MIMO communications.” [Online]. Available: https://arxiv.org/abs/1903.08127
  * [16] X. Yu, D. Xu, and R. Schober, “Enabling secure wireless communications via intelligent reflecting surfaces.” [Online]. Available: https://arxiv.org/abs/1904.09573
  * [17] M. Cui, G. Zhang, and R. Zhang, “Secure wireless communication via intelligent reflecting surface,” _IEEE Wireless Commun. Lett._ , 2019.
  * [18] H. Shen, W. Xu, W. Xu, S. Gong, Z. He, and C. Zhao, “Secrecy rate maximization for intelligent reflecting surface assisted multi-antenna communications,” _IEEE Commun. Lett._ , pp. 1–1, 2019.
  * [19] J. Chen, Y. Liang, Y. Pei, and H. Guo, “Intelligent reflecting surface: A programmable wireless environment for physical layer security.” [Online]. Available: https://arxiv.org/abs/1905.03689
  * [20] D. Xu, X. Yu, Y. Sun, D. W. K. Ng, and R. Schober, “Resource allocation for secure IRS-assisted multiuser MISO systems.” [Online]. Available: https://arxiv.org/abs/1907.03085
  * [21] X. Guan, Q. Wu, and R. Zhang, “Intelligent reflecting surface assisted secrecy communication via joint beamforming and jamming.” [Online]. Available: https://arxiv.org/abs/1907.12839
  * [22] T. Bai, C. Pan, Y. Deng, M. Elkashlan, A. Nallanathan, and L. Hanzo, “Latency minimization for intelligent reflecting surface aided mobile edge computing.” [Online]. Available: https://arxiv.org/abs/1910.07990
  * [23] G. Zhou, C. Pan, H. Ren, K. Wang, and A. Nallanathan, “Intelligent reflecting surface aided multigroup multicast MISO communication systems.” [Online]. Available: https://arxiv.org/abs/1909.04606
  * [24] C. Pan, H. Ren, K. Wang, W. Xu, M. Elkashlan, A. Nallanathan, and L. Hanzo, “Intelligent reflecting surface for multicell MIMO communications.” [Online]. Available: https://arxiv.org/abs/1907.10864
  * [25] C. Huang, G. C. Alexandropoulos, C. Yuan, and M. Debbah, “Indoor signal focusing with deep learning designed reconfigurable intelligent surfaces.” [Online]. Available: https://arxiv.org/abs/1905.07726
  * [26] Z. He and X. Yuan, “Cascaded channel estimation for large intelligent metasurface assisted massive MIMO.” [Online]. Available: https://arxiv.org/abs/1905.07948
  * [27] A. Taha, M. Alrabeiah, and A. Alkhateeb, “Enabling large intelligent surfaces with compressive sensing and deep learning.” [Online]. Available: https://arxiv.org/abs/1904.10136
  * [28] G. Zhou, C. Pan, H. Ren, K. Wang, M. D. Renzo, and A. Nallanathan, “Robust beamforming design for intelligent reflecting surface aided MISO communication systems.” [Online]. Available: https://arxiv.org/abs/1911.06237
  * [29] R. Zhang and C. K. Ho, “MIMO broadcasting for simultaneous wireless information and power transfer,” _IEEE Trans. Wireless Commun._ , vol. 12, no. 5, pp. 1989–2001, May 2013.
  * [30] Q. Wu and R. Zhang, “Weighted sum power maximization for intelligent reflecting surface aided SWIPT.” [Online]. Available: https://arxiv.org/abs/1907.05558
  * [31] Y. Sun, P. Babu, and D. P. Palomar, “Majorization-minimization algorithms in signal processing, communications, and machine learning,” _IEEE Trans. Signal Process._ , vol. 65, no. 3, pp. 794–816, Feb. 2017.
  * [32] C. Pan, W. Xu, J. Wang, H. Ren, W. Zhang, N. Huang, and M. Chen, “Pricing-based distributed energy-efficient beamforming for MISO interference channels,” _IEEE J. Sel. Areas Commun._ , vol. 34, no. 4, pp. 710–722, Apr. 2016.
  * [33] C. Pan, H. Zhu, N. J. Gomes, and J. Wang, “Joint precoding and RRH selection for user-centric green MIMO C-RAN,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 5, pp. 2891–2906, May 2017.
  * [34] E. Boshkovska, D. W. K. Ng, N. Zlatanov, and R. Schober, “Practical non-linear energy harvesting model and resource allocation for SWIPT systems,” _IEEE Commun. Lett._ , vol. 19, no. 12, pp. 2082–2085, 2015.
  * [35] D. Mishra, G. C. Alexandropoulos, and S. De, “Energy sustainable IoT with individual QoS constraints through MISO SWIPT multicasting,” _IEEE Internet Things J._ , vol. 5, no. 4, pp. 2856–2867, Aug. 2018.
  * [36] K. Xiong, B. Wang, and K. R. Liu, “Rate-energy region of SWIPT for MIMO broadcasting under nonlinear energy harvesting model,” _IEEE Trans. Wireless Commun._ , vol. 16, no. 8, pp. 5147–5161, Aug. 2017.
  * [37] F. Wang, J. Xu, X. Wang, and S. Cui, “Joint offloading and computing optimization in wireless powered mobile-edge computing systems,” _IEEE Trans. Wireless Commun._ , vol. 17, no. 3, pp. 1784–1797, Mar. 2018.
  * [38] J. Xu, L. Liu, and R. Zhang, “Multiuser MISO beamforming for simultaneous wireless information and power transfer,” _IEEE Trans. Signal Process._ , vol. 62, no. 18, pp. 4798–4810, Sep. 2014.
  * [39] Q. Shi, M. Razaviyayn, Z.-Q. Luo, and C. He, “An iteratively weighted MMSE approach to distributed sum-utility maximization for a MIMO interfering broadcast channel,” _IEEE Trans. Signal Process._ , vol. 59, no. 9, pp. 4331–4340, Sep. 2011.
  * [40] C. Pan, H. Ren, M. Elkashlan, A. Nallanathan, and L. Hanzo, “The non-coherent ultra-dense C-RAN is capable of outperforming its coherent counterpart at a limited fronthaul capacity,” _IEEE J. Sel. Areas Commun._ , vol. 36, no. 11, pp. 2549–2560, Nov. 2018.
  * [41] C. Pan, H. Ren, M. Elkashlan, A. Nallanathan, and L. Hanzo, “Robust beamforming design for ultra-dense user-centric C-RAN in the face of realistic pilot contamination and limited feedback.” [Online]. Available: https://arxiv.org/abs/1804.03990
  * [42] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1,” 2014.
  * [43] S. Boyd and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [44] C. Pan, W. Xu, W. Zhang, J. Wang, H. Ren, and M. Chen, “Weighted sum energy efficiency maximization in ad hoc networks,” _IEEE Wireless Commun. Lett._ , vol. 4, no. 3, pp. 233–236, Jun. 2015.
  * [45] X.-D. Zhang, _Matrix analysis and applications_. Cambridge University Press, 2017.
  * [46] J. Song, P. Babu, and D. P. Palomar, “Sequence design to minimize the weighted integrated and peak sidelobe levels,” _IEEE Trans. Signal Process._ , vol. 64, no. 8, pp. 2051–2064, Apr. 2016.