# Modulating Intelligent Surfaces for Multi-User MIMO Systems: Beamforming and
Modulation Design

Haseeb Ur Rehman, Faouzi Bellili, Member, IEEE, Amine Mezghani, Member, IEEE,
and Ekram Hossain, Fellow, IEEE The authors are with the Department of
Electrical and Computer Engineering at the University of Manitoba, Canada
(emails<EMAIL_ADDRESS>{Faouzi.Bellili, Amine.Mezghani,
Ekram.Hossain}@umanitoba.ca). This work was supported by the Natural Sciences
and Engineering Research Council of Canada (NSERC) and Futurewei Technologies.

###### Abstract

This paper introduces a novel approach of utilizing the reconfigurable
intelligent surface (RIS) for joint data modulation and signal beamforming in
a multi-user downlink cellular network by leveraging the idea of backscatter
communication. We present a general framework in which the RIS, referred to as
modulating intelligent surface (MIS) in this paper, is used to: $i)$ beamform
the signals for a set of users whose data modulation is already performed by
the base station (BS), and at the same time, $ii)$ embed the data of a
different set of users by passively modulating the deliberately sent carrier
signals from the BS to the RIS. To maximize each user’s spectral efficiency, a
joint non-convex optimization problem is formulated under the sum minimum
mean-square error (MMSE) criterion. Alternating optimization is used to divide
the original joint problem into two tasks of: $i)$ separately optimizing the
MIS phase-shifts for passive beamforming along with data embedding for the BS-
and MIS-served users, respectively, and $ii)$ jointly optimizing the active
precoder and the receive scaling factor for the BS- and MIS-served users,
respectively. While the solution to the latter joint problem is found in
closed-form using traditional optimization techniques, the optimal phase-
shifts at the MIS are obtained by deriving the appropriate optimization-
oriented vector approximate message passing (OOVAMP) algorithm. Moreover, the
original joint problem is solved under both ideal and practical constraints on
the MIS phase shifts, namely, the unimodular constraint and assuming each MIS
element to be terminated by a variable reactive load. The proposed MIS-
assisted scheme is compared against state-of-the-art RIS-assisted wireless
communication schemes and simulation results reveal that it brings substantial
improvements in terms of system throughput while supporting a much higher
number of users.

###### Index Terms:

Reconfigurable intelligent surface (RIS), RIS-assisted multi-user MIMO, active
and passive beamforming, modulation, data encoding, backscatter communication,
modulating intelligent surface (MIS)

##  I Introduction

###  I-A Background

Reconfigurable intelligent surfaces (RISs) have recently emerged as a
promising technology for beyond-5G/6G wireless communications as it can
improve both the spectral and energy efficiencies of wireless systems. RISs
also offer the benefits of low-cost and easy integration into the currently
deployed wireless systems [1, 2, 3, 4]. Specifically, a RIS consists of an
array of reconfigurable passive reflective elements which enable one to
control the wireless propagation environment. RIS is an energy-efficient
technology since it allows to passively beamform the incoming signals without
the need for a power amplifier as required in traditional multiple-input
multiple-output (MIMO) base stations (BSs) [5, 6, 4, 7]. It does so by
suitably optimizing the phase shifts applied by each reflective element to
constructively combine the incoming signals so as to achieve improved received
power at the end users. The active transmit precoding at the BS together with
the RIS phase shifts (for passive beamforming) can be jointly optimized in
order to maximize the spectral-efficiency [8, 9, 10, 11, 7] or the energy-
efficiency of the system [12, 13]. Moreover, a RIS can also be utilized to
simultaneously perform passive beamforming and physical information transfer
[14, 15] (e.g., synchronization data or the channel state information [CSI]
estimated at the RIS). In [16], the authors presented a massive backscatter
wireless communication (MBWC) scheme to encode information on Wi-Fi signals
reflected by the RIS. So far, the vast majority of the work has focused on
utilizing the RIS for passive beamforming and the idea of using it also to
embed information in a multi-user cellular network has not been explored. In
this paper, by leveraging the idea of backscatter communication [16, 17, 18],
we propose a general framework in which the RIS is used for data embedding on
the reflected or re-emitted signals by changing the imdepence of each RIS
reflective element. In the proposed framework, the RIS can be used to either:
$i)$ perform passive beamforming for users served by a BS, or $ii)$ embed
information through backscatter communication, or $iii)$ do both
simultaneously. In this setting, the RIS phase shifts vary with the inevitable
changes in the propagation medium to perform passive beamforming for one set
of users, and also with every transmit symbol vector in order to modulate the
carrier by data for another set of users. We name the smart surface with such
capability as “modulating intelligent surface (MIS)” 111Such type of antenna
surfaces is also known as space-time-coding digital metasurfaces in the
antenna literature [19]..

###  I-B Motivation

The benefits of MISs, that can beamform and modulate signals at the same time,
are three-fold:

  1. 1.

In traditional purely reflective RIS-based schemes, the users’ received
signals are subject to severe attenuations stemming from the product path-loss
of the BS-RIS and RIS-user links since both links are used for data
communication. Using a MIS, however, allows the BS to focus all of the
transmit power (of a reference signal) towards the strongest path in the BS-
MIS link and the modulated signals are passively generated at the MIS by
appropriately designing its phase shifts. In this way, the information-bearing
signals undergo the path-loss of the MIS-user link only before reaching the
intended users.

  2. ii.

The total number of users which can be served is not limited by the number of
channel paths that are available in the BS-MIS link, but rather by the number
of reflective elements in the MIS. Therefore, more users than the number of BS
antennas can be served by allocating power towards a MIS having a high number
of antenna elements and then serving the users through the MIS.

  3. iii.

The MIS can serve users by recycling the incoming signals transmitted by the
BS without any RF chains, thus making the entire approach very cost-effective.

A practical implementation of the proposed MIS-based scheme requires the BS in
part to: $i)$ convey the optimized phase shifts to the MIS through a high
speed communication link, and $ii)$ radiate a carrier RF signal to provide the
MIS with the required power for passively embedding the data of a predefined
set of users. Ideally, the MIS must be installed at a location that is in
line-of-sight (LOS) to the BS to minimize the power loss. There are two
possible ways to transmit the optimized phase shifts to the MIS: $i)$ through
an optical fibre link between the BS and the MIS, $ii)$ through a high-speed
reliable wireless link (e.g., Terahertz wireless links) between the BS and the
MIS, if they are in close proximity to each other.

###  I-C Contributions

We consider a single-cell downlink MIMO system assisted by a single MIS which
is equipped with a large number of passive phase shifter elements. The MIS
helps the BS in better serving one portion of the users through passive
beamforming and also embeds the information-bearing data for the remaining
users on a separate carrier signal that it receives from the BS. In this
regard, we build on our recent work in [8] and propose a general method to
jointly optimize:

  * •

the active BS precoder,

  * •

the receive scaling factors for the BS- and MIS-served users,

  * •

the MIS phase shifts used for passive beamforming towards the BS-served users,
and

  * •

the data embedding on the reflected signal (i.e., modulation of the reflected
signal) for the MIS-served users.

The major contributions embodied by this paper are summarized as follows.

  * •

To the best of our knowledge, this is the first work that studies the use of
RISs for passive beamforming and data embedding at the same time in a multi-
user setup. We solve the problem of maximizing the users’ spectral-
efficiencies by jointly optimizing the transmit precoding matrix at the BS,
the receive scaling factor for the BS- and MIS-served users, and the MIS phase
shifts. By building on our recent work in [8], we formulate the joint
optimization problem under the sum minimum mean-square error (MMSE) criterion
in order to minimize the mean-square error (MSE) of the received symbols for
all users at the same time.

  * •

To solve the underlying joint optimization problem, we first split it using
alternating optimization [20] into two simpler sub-optimization tasks, one for
finding the optimum MIS phase shifts and the other for jointly optimizing the
active BS precoder and the receive scaling factors for both the BS- and MIS-
served users. The solution to the latter sub-optimization task is found in
closed form.

  * •

We apply the OOVAMP algorithm to optimize the MIS phase shifts under two
different constraints, namely, $i)$ the unimodular constraint, and $ii)$
assuming each MIS element to be terminated by a tunable reactive load.

  * •

We present various numerical results to compare the proposed scheme against
the standard approach in which a RIS is used for passive beamforming only
while active MMSE precoding is used at the BS for all users. In this context,
we consider two baseline techniques that rely on $i)$ semi-definite relaxation
(SDR) [21, 9], and $ii)$ OOVAMP-based alternating optimization [8] to find the
adequate RIS phase shifts. The latter approach boils down to a special case of
the herein proposed scheme when the number of MIS-served users is equal to
zero. The simulation results show that using MISs for joint beamforming and
information embedding significantly outperforms the classical schemes in which
RISs are used for passively beamforming the signals they receive from the BS.
We also study the resilience of the proposed scheme under CSI mismatches
stemming from imperfect CSI acquisition in practice.

###  I-D Paper Organization and Notations

The rest of this paper is organized as follows. The system model along with
the problem formulation for jointly optimizing the active BS precoder and the
MIS phase-shifts are discussed in Section II. Section III discusses the OOVAMP
algorithm and its constituent modules. In Section IV, we solve the
optimization problem at hand under multiple constraints using the underlying
OOVAMP algorithm. Lastly, numerical results are presented in Section V before
drawing out the conclusion in Section VI.

Notations: Lowercase letters (e.g., $r$) denote scalar variables. Vectors are
denoted by small boldface letters (e.g., $\mathbf{z}$) and the $k$-th element
of $\mathbf{z}$ is denoted as $z_{k}$. Exponent on a vector (e.g.,
$\mathbf{z}^{n}$) denotes component-wise exponentiation on every element of
the vector. Capital boldface letters (e.g., $\mathbf{A}$) are used to denote
matrices, while $a_{ik}$ and $\bm{a}_{i}$ stand, respectively, for the
$(i,k)$-th entry and the $i$-th column of $\mathbf{A}$. The zero matrix of
size $M\times N$ is denoted as $\mathbf{0}_{M\times N}$.
$\operatorname{\mathbb{C}}^{M\times N}$ stands for the set of matrices of size
$M\times N$ with complex elements and $\mathbf{A}^{-k}$ means
$\left(\mathbf{A}^{-1}\right)^{k}$. $\textnormal{Rank}(\mathbf{A})$ and
$\textnormal{Tr}(\mathbf{A})$, return, respectively, the rank and the trace of
any matrix $\mathbf{A}$. We also use $\|.\|_{\textnormal{2}},\
\|.\|_{\textnormal{F}},\ (.)^{*},\ (.)^{\mathsf{T}},\ (.)^{\mathsf{H}}$ to
denote the $\mathcal{L}_{\textnormal{2}}$ norm, Frobenius norm, the conjugate,
the transpose, and the conjugate transpose operators, respectively. The
operator $<.>$ returns the empirical average of all the elements/entries of
any vector or matrix. Moreover, $\textnormal{vec}(.)$ and
$\textnormal{unvec}(.)$ denote vectorization of a matrix and unvectorization
of a vector back to its original matrix form, respectively.
$\textnormal{Diag}(.)$ operates on a vector and generates a diagonal matrix by
placing that vector in the diagonal whereas $\textnormal{diag}(.)$ operates on
a matrix and returns its main diagonal in a vector. The statistical
expectation is denoted as $\operatorname{\mathbb{E}}\\{.\\}$. A random vector
with complex normal distribution is represented by
$\bm{\mathsf{x}}\sim\operatorname{\mathcal{C}}\operatorname{\mathcal{N}}(\mathbf{x};\mathbf{u},\mathbf{R})$,
where $\mathbf{u}$ and $\mathbf{R}$ denote its mean and covariance matrix,
respectively. Similarly, a random matrix with complex normal distribution is
represented by
$\bm{\mathsf{X}}\sim\operatorname{\mathcal{C}}\operatorname{\mathcal{M}}\operatorname{\mathcal{N}}(\mathbf{X};\mathbf{M},\mathbf{U},\mathbf{V})$,
where $\mathbf{M}$, $\mathbf{U}$ and $\mathbf{V}$ denote its mean and
covariance matrices among its rows and columns, respectively. The imaginary
unit is represented by ${\rm j}=\sqrt{-1}$ and the $\angle(.)$ operator
returns the angle of any complex number. The proportional relationship between
any two entities (functions or variables) is denoted by $\varpropto$. Lastly,
the operators $\otimes$, $\odot$ and $\ast$ stand for the Kronecker, the
Hadamard, and the column-wise Khatri-Rao products, respectively.

##  II System Model, Assumptions, and Problem Formulation

Consider a BS that is equipped with $N$ antenna elements that is serving a
total of $M$ single-antenna users (in the downlink) with the help of a MIS
that has $K>M$ reflective elements. Also consider a scheme in which the data
for $B$ out of the total $M$ users are directly modulated/encoded by the BS
(in baseband). Those $B$ users ($B<N$) are referred to as the BS-served users.
For the remaining $R=M-B$ users, the BS will simply send a known/reference
signal which will then be modulated by appropriately phase-shifting it using
the MIS reflective elements. For this reason, we call those $R$ users as the
MIS-served users although, strictly speaking, both types of users are being
served by the BS222Indeed, although being applied at the MIS, the information-
bearing phase shifts are designed centrally at the BS as function of the
users’ data.. The goal is to optimally design the MIS phase shifts not only to
modulate the data for the MIS-served users but also to passively beamform the
signals for the BS-served users. As illustrated by Fig. 1, each $b$-th BS-
served user has a direct link to the BS which is expressed by a channel vector
$\mathbf{h}_{\textnormal{b-ub},b}\in\operatorname{\mathbb{C}}^{N}$. The BS-
user channel vector for each $r$-th MIS-served user333Although the MIS-served
users do not receive data transmitted by the BS, they will experience
interference from the direct BS-user link of the BS-served users. is denoted
by $\mathbf{h}_{\textnormal{b-us},r}\in\operatorname{\mathbb{C}}^{N}$. The
surface-user channels for the $b$-th BS-served user and the $r$-th MIS-served
user are denoted, respectively, by
$\mathbf{h}_{\textnormal{s-ub},b}\in\operatorname{\mathbb{C}}^{K}$ and
$\mathbf{h}_{\textnormal{s-us},r}\in\operatorname{\mathbb{C}}^{K}$. Let
$\mathbf{H}_{\textnormal{b-s}}\in\operatorname{\mathbb{C}}^{K\times N}$ denote
the channel matrix of the MIMO MIS-BS link with
$\textnormal{Rank}(\mathbf{H}_{\textnormal{b-s}})\geq B$. The signal received
at the MIS is phase-shifted by a diagonal matrix
$\textnormal{Diag}(\bm{\upsilon})\in\operatorname{\mathbb{C}}^{K\times K}$,
where $\bm{\upsilon}\in\operatorname{\mathbb{C}}^{K}$ is the phase-shift
vector that might be subject to various constraints in practice. In
particular, under the unimodular constraint, we have
$\left\lvert\upsilon_{k}\right\rvert=1$ for $k=1,\cdots,K$, or equivalently,
$\upsilon_{k}=e^{{\rm j}\theta_{k}}$ for some phase shift
$\theta_{k}\in[0,2\pi]$. Let
$\mathbf{H}_{\textnormal{b-ub}}=\left[\mathbf{h}_{\textnormal{b-ub},1},\mathbf{h}_{\textnormal{b-ub},2},\cdots,\mathbf{h}_{\textnormal{b-ub},B}\right]\in\operatorname{\mathbb{C}}^{N\times
B}$,
$\mathbf{H}_{\textnormal{b-us}}=\left[\mathbf{h}_{\textnormal{b-us},1},\mathbf{h}_{\textnormal{b-us},2},\cdots,\mathbf{h}_{\textnormal{b-us},R}\right]\in\operatorname{\mathbb{C}}^{N\times
R}$,
$\mathbf{H}_{\textnormal{s-ub}}=[\mathbf{h}_{\textnormal{s-ub},1},\mathbf{h}_{\textnormal{s-ub},2},\cdots,\mathbf{h}_{\textnormal{s-ub},B}]\in\operatorname{\mathbb{C}}^{K\times
B}$, and
$\mathbf{H}_{\textnormal{s-us}}=[\mathbf{h}_{\textnormal{s-us},1},\mathbf{h}_{\textnormal{s-us},2},\cdots,\mathbf{h}_{\textnormal{s-us},R}]\in\operatorname{\mathbb{C}}^{K\times
R}$. For mathematical convenience, we gather the channels between the BS
(resp. MIS) and all users in a single matrix $\mathbf{H}_{\textnormal{b-u}}$
(resp. $\mathbf{H}_{\textnormal{s-u}}$):

$\displaystyle\mathbf{H}_{\textnormal{b-u}}$
$\displaystyle=\left[\mathbf{H}_{\textnormal{b-ub}}~{}\mathbf{H}_{\textnormal{b-us}}\right]\in\operatorname{\mathbb{C}}^{N\times
M},$ (1) $\displaystyle\mathbf{H}_{\textnormal{s-u}}$
$\displaystyle=\left[\mathbf{H}_{\textnormal{s-ub}}~{}\mathbf{H}_{\textnormal{s-us}}\right]~{}\in\operatorname{\mathbb{C}}^{K\times
M}.$ (2)

BSMIS $K$ elementsUser $1$User
$B$$\mathbf{h}_{\textnormal{b-ub},1}$$\mathbf{h}_{\textnormal{b-ub},B}$$\mathbf{H}_{\textnormal{b-s}}$$\mathbf{h}_{\textnormal{s-ub},1}$$\mathbf{h}_{\textnormal{s-ub},B}$$\vdots$$\vdots$$\cdots$$\cdots$$\ddots$User
$1$User
$R$$\mathbf{h}_{\textnormal{s-us},1}$$\mathbf{h}_{\textnormal{s-us},R}$

Figure 1: MIS-assisted multi-user MIMO system in which the MIS is concurrently
used for beamforming and data embedding.

Let $\bar{\mathsf{y}}_{m}$ be the noiseless signal received by the $m$-th user
and define
$\bar{\bm{\mathsf{y}}}\triangleq[\bar{\mathsf{y}}_{1},\bar{\mathsf{y}}_{2},\cdots,\bar{\mathsf{y}}_{M}]^{\mathsf{T}}$.
In the sequel, without loss of generality, we assume that
$\\{\bar{\mathsf{y}}_{1},\bar{\mathsf{y}}_{2},\cdots,\bar{\mathsf{y}}_{B}\\}$
and
$\\{\bar{\mathsf{y}}_{B+1},\bar{\mathsf{y}}_{B+2},\cdots,\bar{\mathsf{y}}_{M}\\}$
pertain to the BS-served and MIS-served users, respectively. Then,
$\bar{\bm{\mathsf{y}}}$ can be decomposed as
$\bar{\bm{\mathsf{y}}}\,=\,\bar{\bm{\mathsf{y}}}_{\textnormal{b}}\,+\,\bar{\bm{\mathsf{y}}}_{\textnormal{s}}$
where $\bar{\bm{\mathsf{y}}}_{\textnormal{b}}$ (resp.
$\bar{\bm{\mathsf{y}}}_{\textnormal{s}}$) is the information-bearing signal
intended to the BS-served (resp. MIS-served) users which are given by:

$\displaystyle\bar{\bm{\mathsf{y}}}_{\textnormal{b}}$
$\displaystyle~{}=~{}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}\bm{\mathsf{s}}_{\textnormal{b}}~{}+~{}\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\bm{\mathsf{s}}_{\textnormal{b}},$
(3) $\displaystyle\bar{\bm{\mathsf{y}}}_{\textnormal{s}}$
$\displaystyle~{}=~{}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon})\mathbf{H}_{\textnormal{b-s}}\left(\sqrt{P_{\textnormal{s}}}\mathbf{v}_{\textnormal{b}}\right).$
(4)

In (3),
$\bm{\mathsf{s}}_{\textnormal{b}}\sim\operatorname{\mathcal{C}}\operatorname{\mathcal{N}}(\mathbf{s};\mathbf{0},\mathbf{I}_{B})$
is the unknown symbol vector being transmitted by the BS (to the B-served
users), and
$\bm{\upsilon}\,=\,[\upsilon_{1},\upsilon_{2},\cdots,\upsilon_{K}]^{\mathsf{T}}\in\operatorname{\mathbb{C}}^{K}$
is a vector that gathers all the phase shifts used by the MIS. Moreover,
$\mathbf{F}\in\operatorname{\mathbb{C}}^{N\times B}$ is the active precoding
matrix that is used for beamforming purposes at the BS, which satisfies
$\left\lVert\mathbf{F}\right\rVert_{\textnormal{F}}^{2}=P_{\textnormal{b}}$,
where $P_{\textnormal{b}}$ is the fraction of power being allocated to the BS-
served users. In (4), $\sqrt{P_{\textnormal{s}}}\mathbf{v}_{\textnormal{b}}$
is a separate reference (which is constant and known) being transmitted by the
BS towards the MIS with
$\left\lVert\mathbf{v}_{\textnormal{b}}\right\rVert_{2}^{2}=1$ and
$P_{\textnormal{s}}$ is the fraction of transmit power being allocated to the
MIS-served users. The total transmit power is denoted by
$P=P_{\textnormal{b}}+P_{\textnormal{s}}$. Now, let
$\bm{\mathsf{w}}_{\textnormal{b}}\sim\operatorname{\mathcal{C}}\operatorname{\mathcal{N}}(\bm{\mathsf{w}};~{}\mathbf{0},~{}\sigma_{\textnormal{w}}^{2}\mathbf{I}_{B})$
and
$\bm{\mathsf{w}}_{\textnormal{s}}\sim\operatorname{\mathcal{C}}\operatorname{\mathcal{N}}(\bm{\mathsf{w}};~{}\mathbf{0},~{}\sigma_{\textnormal{w}}^{2}\mathbf{I}_{R})$
denote the additive white Gaussian noise (AWGN) vectors pertaining to the BS-
served and MIS-served users, respectively. Therefore, the noisy received
signal at all the users, $\bm{\mathsf{y}}$, is given by:

$\bm{\mathsf{y}}=\alpha_{\textnormal{b}}\bar{\bm{\mathsf{y}}}_{\textnormal{b}}+\alpha_{\textnormal{s}}\bar{\bm{\mathsf{y}}}_{\textnormal{s}}+\big{[}\alpha_{\textnormal{b}}\bm{\mathsf{w}}_{\textnormal{b}}^{\mathsf{T}},~{}\alpha_{\textnormal{s}}\bm{\mathsf{w}}_{\textnormal{s}}^{\mathsf{T}}\big{]}^{\mathsf{T}},$
(5)

wherein $\alpha_{\textnormal{b}}$ and $\alpha_{\textnormal{s}}$ are some real-
valued receive scaling factors444This is a common practice in the precoding
optimization literature [22, 23].. They are only utilized to facilitate the
optimization of the other variables (i.e., the BS precoder and the MIS phase
shifts) since the receivers can blindly estimate these scalars based on the
received symbol sequence [22, 23].

######  Remark.

The system model in (5) is an approximation of the exact system model which is
expressed as follows:

$\begin{bmatrix}\alpha_{\textnormal{b}}\mathbf{I}_{B}&\mathbf{0}_{B\times
R}\\\ \mathbf{0}_{R\times
B}&\alpha_{\textnormal{s}}\mathbf{I}_{R}\end{bmatrix}\Big{(}\bar{\bm{\mathsf{y}}}_{\textnormal{b}}+\bar{\bm{\mathsf{y}}}_{\textnormal{s}}+\big{[}\bm{\mathsf{w}}_{\textnormal{b}}^{\mathsf{T}},~{}\bm{\mathsf{w}}_{\textnormal{s}}^{\mathsf{T}}\big{]}^{\mathsf{T}}\Big{)}.$
(6)

The approximation will allow better tractability for the receive scaling
factors, $\alpha_{\textnormal{b}}$ and $\alpha_{\textnormal{s}}$, and the
precoding matrix, $\mathbf{F}$, by decoupling $\alpha_{\textnormal{s}}$ from
the other two variables in the joint optimization problem solved in Section
IV-B. Therefore, we use the approximate model in (5) to optimize the variables
and then use those optimized variables together with the exact system model in
(6) to compute the performance metrics such as sum-rate in Section V.

Consequently, by using (3) and (4) in (5), it follows that:

$\bm{\mathsf{y}}=\alpha_{\textnormal{b}}\left(\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}\bm{\mathsf{s}}_{\textnormal{b}}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\bm{\mathsf{s}}_{\textnormal{b}}\right)+\alpha_{\textnormal{s}}\left(\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon})\mathbf{H}_{\textnormal{b-s}}\sqrt{P_{\textnormal{s}}}\mathbf{v}_{\textnormal{b}}\right)+\big{[}\alpha_{\textnormal{b}}\bm{\mathsf{w}}_{\textnormal{b}}^{\mathsf{T}},~{}\alpha_{\textnormal{s}}\bm{\mathsf{w}}_{\textnormal{s}}^{\mathsf{T}}\big{]}^{\mathsf{T}}.$
(7)

Contrarily to the reference signal $\mathbf{v}_{\textnormal{b}}$, the
optimized MIS phase-shift vector varies in space (with changes in the channel)
and in time (with every transmit symbol vector). We now extend the system
model in (7) for a transmit block of length, $L$, as follows:

$\displaystyle\begin{split}\bm{\mathsf{Y}}&=\alpha_{\textnormal{b}}\Big{[}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{1})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}\bm{\mathsf{s}}_{\textnormal{b},1}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\bm{\mathsf{s}}_{\textnormal{b},1},\cdots,\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{L})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}\bm{\mathsf{s}}_{\textnormal{b},L}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\bm{\mathsf{s}}_{\textnormal{b},L}\Big{]}\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\alpha_{\textnormal{s}}\sqrt{P_{\textnormal{s}}}\Big{[}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{1})\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}},\cdots,\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{L})\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}}\Big{]}\,+\,\Big{[}\alpha_{\textnormal{b}}\bm{\mathsf{W}}_{\textnormal{b}}^{\mathsf{T}},~{}\alpha_{\textnormal{s}}\bm{\mathsf{W}}_{\textnormal{s}}^{\mathsf{T}}\Big{]}^{\mathsf{T}}.\end{split}$
(8)

We further gather the information symbols for the BS-served (resp. MIS-served)
users in one matrix denoted as
$\bm{\mathsf{S}}_{\textnormal{b}}\triangleq[\bm{\mathsf{s}}_{\textnormal{b},1},\cdots,\bm{\mathsf{s}}_{\textnormal{b},L}]\in\operatorname{\mathbb{C}}^{B\times
L}$ (resp. $\mathbf{S}_{\textnormal{s}}\in\operatorname{\mathbb{C}}^{R\times
L}$) and we define
$\bm{\mathsf{S}}\triangleq\big{[}\bm{\mathsf{S}}_{\textnormal{b}}^{\mathsf{T}},~{}\mathbf{S}_{\textnormal{s}}^{\mathsf{T}}\big{]}^{\mathsf{T}}$.
We also stack the phase-shift vectors being used at all time indices
$l=1,\cdots,L$ in one single matrix
$\bm{\Upsilon}=[\bm{\upsilon}_{1},~{}\bm{\upsilon}_{2},\cdots,\bm{\upsilon}_{L}]\in\operatorname{\mathbb{C}}^{K\times
L}$. The goal is to maximize each user’s signal-to-interference-plus-noise
ratio (SINR) by minimizing the error in its received symbols under the sum
MMSE criterion. A lower bound on the spectral efficiency for user $m$ can be
expressed in terms of the MMSE of its received symbol error [24] as follows:

$C_{m}^{\textnormal{MMSE}}=\log_{2}\left(\frac{1}{\textnormal{MMSE}_{m}}\right).$
(9)

The MSE of the received symbol for user $m$ and time index $l$ is given by
$\operatorname{\mathbb{E}}_{\mathsf{y}_{ml},\mathsf{s}_{ml}}\left\\{\left\lvert\mathsf{y}_{ml}-\mathsf{s}_{ml}\right\rvert^{2}\right\\}$,
and for $M$ users and a transmit block of length $L$, the sum symbol MSE is
given by:

$\displaystyle\sum_{m=1}^{M}\sum_{l=1}^{L}\operatorname{\mathbb{E}}_{\mathsf{y}_{ml},\mathsf{s}_{ml}}\left\\{\left\lvert\mathsf{y}_{ml}-\mathsf{s}_{ml}\right\rvert^{2}\right\\}\
=\
\operatorname{\mathbb{E}}_{\bm{\mathsf{Y}},\bm{\mathsf{S}}}\left\\{\left\lVert\bm{\mathsf{Y}}-\bm{\mathsf{S}}\right\rVert_{\textnormal{F}}^{2}\right\\}.$
(10)

Thus, our optimization problem under the sum MMSE criterion can be formulated
as follows:

$\displaystyle\displaystyle\operatorname*{arg\,min}_{\alpha_{\textnormal{s}},\alpha_{\textnormal{b}},\mathbf{F},{\mathbf{\Upsilon}}}\quad$
$\displaystyle\operatorname{\mathbb{E}}_{\bm{\mathsf{Y}},\bm{\mathsf{S}}}\left\\{\left\lVert\bm{\mathsf{Y}}-\bm{\mathsf{S}}\right\rVert_{\textnormal{F}}^{2}\right\\},$
(11a) subject to
$\displaystyle\left\lVert\mathbf{F}\right\rVert_{\textnormal{F}}^{2}=P_{\textnormal{b}},$
(11b) $\displaystyle|\upsilon_{kl}|=1,\quad
k=1,2,\cdots,K,~{}~{}l=1,2,\cdots,L.$ (11c)

We take the expectation involved in (11a) with respect to (w.r.t.) the random
matrices $\bm{\mathsf{S}}$, $\bm{\mathsf{W}}_{\textnormal{b}}$, and
$\bm{\mathsf{W}}_{\textnormal{s}}$ to further simplify the objective function
(see Appendix A) thereby resulting in the following optimization problem:

$\displaystyle\begin{split}\displaystyle\operatorname*{arg\,min}_{\alpha_{\textnormal{s}},\alpha_{\textnormal{b}},\mathbf{F},{\mathbf{\Upsilon}}}\quad&\bigg{\|}\alpha_{\textnormal{b}}\Big{[}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{1})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F},\cdots,\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{L})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\Big{]}\\\
&-\Big{[}[\mathbf{I}_{B,1},\ \mathbf{0}_{B\times
R,1}]^{\mathsf{T}},\cdots,[\mathbf{I}_{B,L},\ \mathbf{0}_{B\times
R,L}]^{\mathsf{T}}\Big{]}\bigg{\|}_{\textnormal{F}}^{2}\\\
&+\bigg{\|}\alpha_{\textnormal{s}}\sqrt{P_{\textnormal{s}}}\Big{[}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{1})\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}},\cdots,\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{L})\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}}\Big{]}-\left[\mathbf{0}_{L\times
B}\
\mathbf{S}_{\textnormal{s}}^{\mathsf{T}}\right]^{\mathsf{T}}\bigg{\|}_{\textnormal{F}}^{2}\\\
&+LB\sigma_{\textnormal{w}}^{2}\alpha_{\textnormal{b}}^{2}+LR\sigma_{\textnormal{w}}^{2}\alpha_{\textnormal{s}}^{2}.\end{split}$
(12a) s.t.
$\displaystyle\|\mathbf{F}\|_{\textnormal{F}}^{2}=P_{\textnormal{b}},$ (12b)
$\displaystyle|\upsilon_{kl}|=1,\quad k=1,2,\cdots,K,~{}~{}l=1,2,\cdots,L.$
(12c)

We shall denote the objective function in (12a) by
$f(\alpha_{\textnormal{s}},\alpha_{\textnormal{b}},\mathbf{F},{\mathbf{\Upsilon}})$
throughout the rest of the paper. The optimization problem in (12) is non-
convex due to the unimodular constraint on the MIS phase-shifts in (12c). We
aim to solve the problem by utilizing the OOVAMP algorithm in the same way as
done in [8].

The OOVAMP [8] is a low-complexity algorithm which is designed to solve
optimization problems involving a linear objective function and non-linear
constraints [25, 8]. The modular structure of OOVAMP decouples the constraint
from the objective function. This favorable property allows OOVAMP to optimize
the same objective function under different constraints by modifying the
OOVAMP module that satisfies the constraint (simple scalar functions). OOVAMP
has better convergence properties than other iterative algorithms (e.g.,
alternating direction method of multipliers [ADMM]) [26] because of its
automatic calculation of the optimal stepsize in every iteration.

##  III Optimization Oriented VAMP for Matrices

OOVAMP is a recently developed extension [8] of the standard max-sum VAMP
algorithm[25] which solves constrained optimization problems involving linear
objective functions under both linear and non-linear constraints. Moreover,
asymptotic optimality can be claimed for the computed solution under certain
mild conditions using state evolution arguments. In this section, we present
an extended version of the OOVAMP algorithm to optimize matrices involving
linear mixing. Given the knowledge of two matrices,
$\mathbf{A}\in\operatorname{\mathbb{C}}^{M\times N}$ and
$\mathbf{Z}\in\operatorname{\mathbb{C}}^{M\times K}$, the OOVAMP algorithm
solves the following optimization problem:

$\displaystyle\displaystyle\operatorname*{arg\,min}_{\mathbf{X}~{}\in~{}\operatorname{\mathbb{C}}^{N\times
K}}\quad$
$\displaystyle\left\lVert\mathbf{A}\mathbf{X}-\mathbf{Z}\right\rVert_{\textnormal{F}}^{2}$
(13a) s.t. $\displaystyle f_{ik}(x_{ik})=0\quad i=1,\cdots,N,\quad
k=1,\cdots,K.$ (13b)

The algorithm consists of the following two modules.

####  III-1 Linear MAP Estimator

At iteration $t$, the linear MAP (LMAP) estimator receives extrinsic
information (message) from the separable (i.e., entry-wise) MAP projector of
$\mathbf{X}$ in the form of a mean matrix, $\mathbf{R}_{t-1}$, and a common
scalar precision, $\gamma_{t-1}$. Then, under the Gaussian prior,
$\operatorname{\mathcal{C}}\operatorname{\mathcal{M}}\operatorname{\mathcal{N}}\Big{(}\mathbf{X};\mathbf{R}_{t-1},\gamma_{t-1}^{-1}\mathbf{I}_{N},\mathbf{I}_{K}\Big{)}$,
it computes the LMAP estimate, $\bar{\mathbf{X}}_{t}$, along with the
associated posterior precision, $\bar{\gamma}_{t}$, as follows:

$\displaystyle\bar{\mathbf{X}}_{t}$
$\displaystyle~{}=~{}\left(\mathbf{A}^{\mathsf{H}}\mathbf{A}+\gamma_{t-1}\mathbf{I}_{N}\right)^{-1}\left(\mathbf{A}^{\mathsf{H}}\mathbf{Z}+\gamma_{t-1}\mathbf{R}_{t-1}\right),$
(14) $\displaystyle\bar{\gamma}_{t}$
$\displaystyle~{}=~{}N\textnormal{Tr}\left(\left[\mathbf{A}^{\mathsf{H}}\mathbf{A}+\gamma_{t-1}\mathbf{I}_{N}\right]^{-1}\right)^{-1}.$
(15)

The extrinsic information on $\mathbf{X}$ is updated as:

$\operatorname{\mathcal{C}}\operatorname{\mathcal{M}}\operatorname{\mathcal{N}}(\mathbf{X};\bar{\mathbf{X}}_{t},\bar{\gamma}_{t}^{-1}\mathbf{I}_{N},\mathbf{I}_{K})/\operatorname{\mathcal{C}}\operatorname{\mathcal{M}}\operatorname{\mathcal{N}}(\mathbf{X};\mathbf{R}_{t-1},\gamma_{t-1}^{-1}\mathbf{I}_{N},\mathbf{I}_{K}),$

and then sent back in the form of a mean matrix,
$\widetilde{\mathbf{R}}_{t}=\left(\bar{\mathbf{X}}_{t}\bar{\gamma}_{t}-\mathbf{R}_{t-1}\gamma_{t-1}\right)/\left(\bar{\gamma}_{t}-\gamma_{t-1}\right)$,
and a scalar precision,
$\widetilde{\gamma}_{t}=\bar{\gamma}_{t}-\gamma_{t-1}$, to the separable MAP
projector of $\mathbf{X}$.

####  III-2 Separable MAP Projector

Because the constraint on $\mathbf{X}$ is component-wise, the constraint on
its entries, $x_{ik}$, is modeled as a prior with some precision,
$\gamma_{\textnormal{p}}$, i.e.,
$p_{\mathsf{x}}(x_{ik})\varpropto\exp\left(-\gamma_{\textnormal{p}}\left\lvert
f_{ik}(x_{ik})\right\rvert^{2}\right)$ with
$\gamma_{\textnormal{p}}\to\infty$, which results in the following prior
distribution on $\mathbf{X}$:

$p_{\bm{\mathsf{X}}}(\mathbf{X})=\displaystyle\prod_{i=1}^{N}\prod_{k=1}^{Q}p_{\mathsf{x}}(x_{ik}).$
(16)

This module computes the MAP estimate, $\widehat{\mathbf{X}}_{t}$, of
$\mathbf{X}$ from the joint distribution
$p_{\bm{\mathsf{X}}}(\mathbf{X})\operatorname{\mathcal{C}}\operatorname{\mathcal{M}}\operatorname{\mathcal{N}}\Big{(}\mathbf{X};\widetilde{\mathbf{R}}_{t},\widetilde{\gamma}^{-1}_{t}\mathbf{I}_{N},\mathbf{I}_{K}\Big{)}$.
The MAP estimate can be computed through a component-wise projector function
as follows:

$\widehat{x}_{ik,t}~{}=~{}g_{ik}(\widetilde{r}_{ik,t},\widetilde{\gamma}_{t})~{}\triangleq~{}\displaystyle\operatorname*{arg\,min}_{x_{ik}}\left[\widetilde{\gamma}_{t}\left\lvert
x_{ik}-\widetilde{r}_{ik,t}\right\rvert^{2}-\ln
p_{\mathsf{x}}(x_{ik})\right],$ (17)

or equivalently:

$g_{ik}(\widetilde{r}_{ik,t},\widetilde{\gamma}_{t})=\operatorname*{arg\,min}_{x_{ik}}\left[\widetilde{\gamma}_{t}\left\lvert
x_{ik}-\widetilde{r}_{ik,t}\right\rvert^{2}+\gamma_{\textnormal{p}}\left\lvert
f_{ik}(x_{ik})\right\rvert^{2}\right].$ (18)

The parameter $\gamma_{\textnormal{p}}$ in (18) accounts for the weight given
to the prior on $x_{i}$ inside the scalar MAP optimization. Therefore, taking
$\gamma_{p}\to\infty$ enforces the constraint. The derivative of the scalar
MAP projector w.r.t. $\widetilde{r}_{ik,t}$ is given by [25]:

$g_{ik}^{\prime}(\widetilde{r}_{ik,t},\widetilde{\gamma}_{t})~{}\triangleq~{}\frac{\partial
g_{ik}(\widetilde{r}_{ik,t},\widetilde{\gamma}_{t})}{\partial\widetilde{r}_{ik,t}}~{}=~{}\frac{1}{2}\left(\frac{\partial
g_{ik}\left(\widetilde{r}_{ik,t},\widetilde{\gamma}_{t}\right)}{\partial\Re\left\\{\widetilde{r}_{ik,t}\right\\}}-{\rm
j}\frac{\partial
g_{ik}\left(\widetilde{r}_{ik,t},\widetilde{\gamma}_{t}\right)}{\partial\Im\left\\{\widetilde{r}_{ik,t}\right\\}}\right)~{}=~{}\widetilde{\gamma}_{t}\widehat{\gamma}_{t},$
(19)

where $\widehat{\gamma}_{t}$ is the posterior precision. The matrix-valued
projector function and its derivative are defined as follows:

$\mathbf{G}(\widetilde{\mathbf{R}}_{t},\widetilde{\gamma}_{t})~{}\triangleq~{}\begin{bmatrix}g_{11}(\widetilde{r}_{11,t},\widetilde{\gamma}_{t})&g_{12}(\widetilde{r}_{12,t},\widetilde{\gamma}_{t})&\cdots&g_{1K}(\widetilde{r}_{1K,t},\widetilde{\gamma}_{t})\\\
g_{21}(\widetilde{r}_{21,t},\widetilde{\gamma}_{t})&g_{22}(\widetilde{r}_{22,t},\widetilde{\gamma}_{t})&\cdots&g_{2K}(\widetilde{r}_{2K,t},\widetilde{\gamma}_{t})\\\
\vdots&\vdots&\ddots&\vdots\\\
g_{N1}(\widetilde{r}_{N1,t},\widetilde{\gamma}_{t})&g_{N2}(\widetilde{r}_{N2,t},\widetilde{\gamma}_{t})&\cdots&g_{NK}(\widetilde{r}_{NK,t},\widetilde{\gamma}_{t})\end{bmatrix},$
(20)

$\mathbf{G}^{\prime}(\widetilde{\mathbf{R}}_{t},\widetilde{\gamma}_{t})~{}\triangleq~{}\begin{bmatrix}g^{\prime}_{11}(\widetilde{r}_{11,t},\widetilde{\gamma}_{t})&g^{\prime}_{12}(\widetilde{r}_{12,t},\widetilde{\gamma}_{t})&\cdots&g^{\prime}_{1K}(\widetilde{r}_{1K,t},\widetilde{\gamma}_{t})\\\
g^{\prime}_{21}(\widetilde{r}_{21,t},\widetilde{\gamma}_{t})&g^{\prime}_{22}(\widetilde{r}_{22,t},\widetilde{\gamma}_{t})&\cdots&g^{\prime}_{2K}(\widetilde{r}_{2K,t},\widetilde{\gamma}_{t})\\\
\vdots&\vdots&\ddots&\vdots\\\
g^{\prime}_{N1}(\widetilde{r}_{N1,t},\widetilde{\gamma}_{t})&g^{\prime}_{N2}(\widetilde{r}_{N2,t},\widetilde{\gamma}_{t})&\cdots&g^{\prime}_{NK}(\widetilde{r}_{NK,t},\widetilde{\gamma}_{t})\end{bmatrix}.$
(21)

Similar to the LMAP module, the MAP projector module computes an extrinsic
mean matrix,
$\mathbf{R}_{t}=\left(\widehat{\mathbf{X}}_{t}\widehat{\gamma}_{t}-\widetilde{\mathbf{R}}_{t}\widetilde{\gamma}_{t}\right)/\left(\widehat{\gamma}_{t}-\widetilde{\gamma}_{t}\right)$,
and a scalar precision,
$\gamma_{t}=\widehat{\gamma}_{t}-\widetilde{\gamma}_{t}$, and sends them back
to the LMAP module for the next iteration. The process is repeated until
convergence.

It is worth mentioning that the extrinsic parameters, i.e., the extrinsic mean
matrix and the scalar precision, calculated by each module act as a Gaussian
prior on the subsequent estimate of the adjacent module, thus making OOVAMP
parameter-free. Another major advantage of OOVAMP is that it decouples the
constraint from the objective function and also allows the projector function
to be separable. While the LMAP module optimizes the objective function with
no constraints, the latter are enforced by the projector function. This
modular property makes OOVAMP an attractive algorithm for solving optimization
problems involving linear mixing and under various component-wise constraints.
The algorithmic steps of OOVAMP are shown in Algorithm 1.

Algorithm 1 Optimization-oriented max-sum matrix VAMP

Given $\mathbf{A}\in\operatorname{\mathbb{C}}^{M\times N}$,
$\mathbf{Z}\in\operatorname{\mathbb{C}}^{M\times Q}$, a precision tolerance
$(\epsilon)$ and a maximum number of iterations $(T_{\textnormal{MAX}})$

1: Initialize $\mathbf{R}_{0}$, $\gamma_{\textnormal{0}}\geq 0$ and
$t\leftarrow 1$

2: repeat

3: // LMAP.

4:
$\bar{\mathbf{X}}_{t}=\left(\mathbf{A}^{\mathsf{H}}\mathbf{A}+\gamma_{t-1}\mathbf{I}_{N}\right)^{-1}\left(\mathbf{A}^{\mathsf{H}}\mathbf{Z}+\gamma_{t-1}\mathbf{R}_{t-1}\right)$

5:
$\bar{\gamma}_{t}=N\textnormal{Tr}\left(\left[\mathbf{A}^{\mathsf{H}}\mathbf{A}+\gamma_{t-1}\mathbf{I}_{N}\right]^{-1}\right)^{-1}$

6: $\widetilde{\gamma}_{t}=\bar{\gamma}_{t}-\gamma_{t-1}$

7:
$\widetilde{\mathbf{R}}_{t}=\widetilde{\gamma}_{t}^{-1}\left(\bar{\mathbf{X}}_{t}\bar{\gamma}_{t}-\mathbf{R}_{t-1}\gamma_{t-1}\right)$

8: // Separable MAP Projector

9:
$\widehat{\mathbf{X}}_{t}=\mathbf{G}\Big{(}\widetilde{\mathbf{R}}_{t},\widetilde{\gamma}_{t}\Big{)}$

10:
$\widehat{\gamma}_{t}=\widetilde{\gamma}_{t}^{-1}\left\langle\mathbf{G}^{\prime}\Big{(}\widetilde{\mathbf{R}}_{t},\widetilde{\gamma}_{t}\Big{)}\right\rangle$

11: $\gamma_{t}=\widehat{\gamma}_{t}-\widetilde{\gamma}_{t}$

12:
$\mathbf{R}_{t}=\gamma_{t}^{-1}\Big{(}\widehat{\gamma}_{t}\widehat{\mathbf{X}}_{t}-\widetilde{\gamma}_{t}\widetilde{\mathbf{R}}_{t}\Big{)}$

13: $t\leftarrow t+1$

14: until
$\left\lVert\widehat{\mathbf{X}}_{t}-\widehat{\mathbf{X}}_{t-1}\right\rVert_{2}^{2}\leq\epsilon\left\lVert\widehat{\mathbf{X}}_{t-1}\right\rVert_{2}^{2}$
or $t>T_{\textnormal{MAX}}$

15: return $\widehat{\mathbf{X}}_{t}$

##  IV OOVAMP-Based Solution for the Optimization Problem

In this section, we apply the OOVAMP algorithm, introduced in Section III, to
simultaneously optimize the matrix of phase shifters, ${\mathbf{\Upsilon}}$,
the optimal precoding matrix $\mathbf{F}$, and the scaling factors,
$\alpha_{\textnormal{b}}$ and $\alpha_{\textnormal{s}}$ (i.e., to solve the
optimization problem in (12)). We follow the optimization procedure presented
in [8], and decouple the joint optimization problem into two sub-problems
through alternating optimization. In one side, we optimize
${\mathbf{\Upsilon}}$ by utilizing OOVAMP and, on the other side, we find the
optimal transmit precoding, $\mathbf{F}$, and scalars
$\alpha_{\textnormal{b}}$ and $\alpha_{\textnormal{s}}$. It is a simple
iterative approach that optimizes a subset of all variables at a time while
fixing the other set of variables and the process is repeated until
convergence. More specifically, we divide the optimization problem in (11)
into the following two sub-optimization problems:

  1. 1.

$\displaystyle\widehat{\bm{\Upsilon}}=\displaystyle\operatorname*{arg\,min}_{{\mathbf{\Upsilon}}}\quad$
$\displaystyle
f(\alpha_{\textnormal{s}},\alpha_{\textnormal{b}},\mathbf{F},{\mathbf{\Upsilon}})$
(22a) s.t.
$\displaystyle|\upsilon_{kl}|=1,~{}~{}k=1,2,\cdots,K,~{}~{}l=1,2,\cdots,L.$
(22b)

  2. 2.

$\displaystyle\operatorname*{arg\,min}_{\alpha_{\textnormal{s}},\alpha_{\textnormal{b}},\mathbf{F}}\quad$
$\displaystyle
f(\alpha_{\textnormal{s}},\alpha_{\textnormal{b}},\mathbf{F},\widehat{\mathbf{\Upsilon}})$
(23a) s.t.
$\displaystyle\left\lVert\mathbf{F}\right\rVert_{\textnormal{F}}^{2}=P_{\textnormal{b}}.\mspace{210.0mu}$
(23b)

###  IV-A Optimizing the MIS Phase Shifts

Here, we derive the OOVAMP modules (i.e., LMAP estimator and separable
projector function) to solve the sub-optimization problem in (22) which is
restated explicitly as follows:

$\displaystyle\begin{split}\displaystyle\operatorname*{arg\,min}_{{\mathbf{\Upsilon}}}\quad&\bigg{\|}\alpha_{\textnormal{b}}\Big{[}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{1})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F},\cdots,\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{L})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\Big{]}\\\
&-\Big{[}[\mathbf{I}_{B,1},\ \mathbf{0}_{B\times
R,1}]^{\mathsf{T}},\cdots,[\mathbf{I}_{B,L},\ \mathbf{0}_{B\times
R,L}]^{\mathsf{T}}\Big{]}\bigg{\|}_{\textnormal{F}}^{2}\\\
&+\bigg{\|}\alpha_{\textnormal{s}}\sqrt{P_{\textnormal{s}}}\Big{[}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{1})\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}},\cdots,\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{L})\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}}\Big{]}-\left[\mathbf{0}_{L\times
B}\
\mathbf{S}_{\textnormal{s}}^{\mathsf{T}}\right]^{\mathsf{T}}\bigg{\|}_{\textnormal{F}}^{2}\\\
&+LB\sigma_{\textnormal{w}}^{2}\alpha_{\textnormal{b}}^{2}+LR\sigma_{\textnormal{w}}^{2}\alpha_{\textnormal{s}}^{2}\end{split}$
(24a) s.t.
$\displaystyle|\upsilon_{kl}|=1,~{}~{}k=1,2,\cdots,K,~{}~{}l=1,2,\cdots,L.$
(24b)

Next, we re-express the objective function in (24a) in a form that is similar
to the general OOVAMP objective function in (13a). In fact, by introducing the
following matrices:

$\displaystyle\mathbf{A}$
$\displaystyle=\alpha_{\textnormal{b}}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}},$
$\displaystyle(M\times K),$ (25) $\displaystyle\mathbf{B}$
$\displaystyle=(\mathbf{H}_{\textnormal{b-s}}\mathbf{F})^{\mathsf{T}},$
$\displaystyle(B\times K),$ (26) $\displaystyle\mathbf{D}$
$\displaystyle=\left[\mathbf{b}_{1}\otimes\mathbf{a}_{1},\cdots,\mathbf{b}_{K}\otimes\mathbf{a}_{K}\right],$
$\displaystyle(MB\times K),$ (27) $\displaystyle\mathbf{M}$
$\displaystyle=\sqrt{P_{\textnormal{s}}}\alpha_{\textnormal{s}}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}}),$
$\displaystyle(M\times K),$ (28) $\displaystyle\mathbf{X}$
$\displaystyle=\big{[}\textnormal{vec}\left([\mathbf{I}_{B},\
\mathbf{0}_{B\times
R}]^{\mathsf{T}}-\alpha_{\textnormal{b}}\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\right),\cdots,\textnormal{vec}\left([\mathbf{I}_{B},\
\mathbf{0}_{B\times
R}]^{\mathsf{T}}-\alpha_{\textnormal{b}}\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\right)\big{]},$
$\displaystyle(MB\times L),$ (29) $\displaystyle\mathbf{Z}$
$\displaystyle=\left[\mathbf{0}_{L\times B},\
\mathbf{S}_{\textnormal{s}}^{\mathsf{T}}\right]^{\mathsf{T}},$
$\displaystyle(M\times L),$ (30)

we show in Appendix B that the optimization problem in (24) can be rewritten
as follows:

$\displaystyle\displaystyle\operatorname*{arg\,min}_{{\mathbf{\Upsilon}}}\quad$
$\displaystyle\left\lVert\mathbf{D}{\bm{\Upsilon}}-\mathbf{X}\right\rVert_{\textnormal{F}}^{2}~{}+~{}\left\lVert\mathbf{M}{\bm{\Upsilon}}-\mathbf{Z}\right\rVert_{\textnormal{F}}^{2}~{}+~{}LB\sigma_{\textnormal{w}}^{2}\alpha_{\textnormal{b}}^{2}~{}+~{}LR\sigma_{\textnormal{w}}^{2}\alpha_{\textnormal{s}}^{2}$
(31a) s.t.
$\displaystyle|\upsilon_{kl}|=1,~{}~{}k=1,2,\cdots,K,~{}~{}l=1,2,\cdots,L.$
(31b)

The steps to derive both OOVAMP modules are detailed in the sequel.

####  IV-A1 LMAP Estimator

The LMAP module performs the minimization of the objection function in (31a)
under the Gaussian prior,
$\operatorname{\mathcal{C}}\operatorname{\mathcal{M}}\operatorname{\mathcal{N}}\Big{(}{\bm{\Upsilon}};\mathbf{R}_{t-1},\gamma_{t-1}^{-1}\mathbf{I}_{K},\mathbf{I}_{L}\Big{)}$,
by solving the following unconstrained optimization problem:

$\displaystyle\operatorname*{arg\,min}_{{\bm{\Upsilon}}}\quad\frac{1}{2}\left\lVert\mathbf{D}{\bm{\Upsilon}}-\mathbf{X}\right\rVert_{\textnormal{F}}^{2}+\frac{1}{2}\left\lVert\mathbf{M}{\bm{\Upsilon}}-\mathbf{Z}\right\rVert_{\textnormal{F}}^{2}+\frac{\gamma_{t-1}}{2}\left\lVert{\bm{\Upsilon}}-\mathbf{R}_{t-1}\right\rVert_{\textnormal{F}}^{2}.$
(32)

The solution (i.e., the LMAP estimate and the associated posterior precision)
to the optimization problem in (32) is given as follows:

$\displaystyle\bar{\bm{\Upsilon}}_{t}$
$\displaystyle~{}=~{}\left(\mathbf{D}^{\mathsf{H}}\mathbf{D}+\mathbf{M}^{\mathsf{H}}\mathbf{M}+\gamma_{t-1}\mathbf{I}_{K}\right)^{-1}\left(\mathbf{D}^{\mathsf{H}}\mathbf{X}+\mathbf{M}^{\mathsf{H}}\mathbf{Z}+\gamma_{t-1}\mathbf{R}_{t-1}\right),$
(33) $\displaystyle\bar{\gamma}_{t}$
$\displaystyle~{}=~{}K\textnormal{Tr}\left(\left[\mathbf{D}^{\mathsf{H}}\mathbf{D}+\mathbf{M}^{\mathsf{H}}\mathbf{M}+\gamma_{t-1}\mathbf{I}_{K}\right]^{-1}\right)^{-1}.$
(34)

####  IV-A2 Scalar MAP Projector

In this section, we introduce two projector functions defined in [8] to
satisfy the two types of constraint on the MIS reflection coefficients, i.e.,
$i)$ the unimodular constraint, and $ii)$ a practical constraint on the MIS
phase shifts in which each antenna element is terminated by a variable
reactive load. We refer the reader to [8] for more details on the derivation
and differences in the two projector functions. The projector function and its
derivative for the unimodular constraint is given as follows:

$\displaystyle g_{1,kl}\left(\widetilde{r}_{kl}\right)$
$\displaystyle~{}=~{}\widetilde{r}_{kl}\left\lvert\widetilde{r}_{kl}\right\rvert^{-1},$
(35) $\displaystyle g^{\prime}_{1,kl}\left(\widetilde{r}_{kl}\right)$
$\displaystyle~{}=~{}\frac{1}{2}\left\lvert\widetilde{r}_{kl}\right\rvert^{-1},$
(36)

where the derivative is taken according to (19).

To optimize the MIS phase shifts under a practical constraint, we consider a
reflective element that is combined with a tunable reactive load555This can be
implemented for instance by an antenna array composed of omni-directional
dipole elements loaded with the reactive elements in the absence of a ground
plane to allow for bidirectional beamforming and not just hemispherical
coverage. instead of an ideal phase-shifter666The value $1$ is the normalized
resistive part of the element impedance whereas $\chi_{kl}$ is the normalized
reactive part of the antenna plus reactive termination. Accordingly
$\upsilon_{kl}$ represents the induced current flowing across the antenna. We
assume the antenna elements to be uncoupled which holds approximately for
half-wavelength element spacing., i.e., $\upsilon_{kl}=-(1+{\rm
j}\chi_{kl})^{-1}$, where $\chi_{kl}\in\operatorname{\mathbb{R}}$ is a scalar
reactance value that must be optimized for each reflection coefficient. Under
the unimodular constraint, the idealistic MIS has a full field of view (FOV)
and the reflection coefficients correspond to ideal phase-shifters and are of
the form $\upsilon_{kl}=e^{{\rm j}\theta_{kl}}$, where777From a practical
standpoint, this assumption is difficult to realize. With the reactive-loading
constraint, the assumption of a MIS with full field of view (FOV) becomes more
acceptable. $\theta_{kl}\in[0,2\pi]$, whereas under the practical constraint
we have a restriction on the possible values of the MIS phase-shifts, i.e.,
$\angle-(1+{\rm
j}\chi)^{-1}\in\big{[}\textnormal{-}\frac{\pi}{2},\frac{\pi}{2}\big{]}$.
Moreover, the magnitude of each phase-shift under this constraint is always
less than $1$ for any $\chi\neq 0$. Practically, this introduces the phase-
dependent amplitude attenuation in the incident wave. The projector function
under this new constraint is defined as:

$g_{2,kl}\left(\widetilde{r}_{kl},\widetilde{\gamma}\right)~{}\triangleq~{}\displaystyle\operatorname*{arg\,min}_{\upsilon_{kl}}\left[\widetilde{\gamma}\left\lvert\upsilon_{kl}-\widetilde{r}_{kl}\right\rvert^{2}+\gamma_{\textnormal{p}}\left\lvert\upsilon_{kl}+\frac{1}{1+{\rm
j}\chi_{kl}^{\textnormal{opt}}}\right\rvert^{2}\right],$ (37)

where

$\chi_{kl}^{\textnormal{opt}}~{}=~{}g_{3}\left(\widetilde{r}_{kl}\right)~{}\triangleq~{}\displaystyle\operatorname*{arg\,min}_{\chi_{kl}}\left\lvert\widetilde{r}_{kl}+\frac{1}{1+{\rm
j}\chi_{kl}}\right\rvert^{2}.$ (38)

The solution to the optimization problem in (38), the projector function, and
its derivative are given by:

$\displaystyle g_{3}\left(\widetilde{r}_{kl}\right)$
$\displaystyle~{}=~{}\frac{1}{2\Im\left\\{\widetilde{r}_{kl}\right\\}}\left(1+2\Re\left\\{\widetilde{r}_{kl}\right\\}+\sqrt{\left(1+2\Re\left\\{\widetilde{r}_{kl}\right\\}\right)^{2}+4\Im\left\\{\widetilde{r}_{kl}\right\\}^{2}}\right),$
(39) $\displaystyle g_{2,kl}\left(\widetilde{r}_{kl}\right)$
$\displaystyle~{}=~{}-\left(1+{\rm
j}g_{3}\left(\widetilde{r}_{kl}\right)\right)^{-1},$ (40) $\displaystyle
g^{\prime}_{2,kl}\left(\widetilde{r}_{kl}\right)$
$\displaystyle~{}=~{}\left\lvert{\rm
j}g^{\prime}_{3}\left(\widetilde{r}_{kl}\right)\left(1+{\rm
j}g_{3}\left(\widetilde{r}_{kl}\right)\right)^{-2}\right\rvert.$ (41)

The matrix valued projector functions,
$\mathbf{G}_{1}(\widetilde{\mathbf{R}}_{t},\widetilde{\gamma}_{t})$ and
$\mathbf{G}_{2}(\widetilde{\mathbf{R}}_{t},\widetilde{\gamma}_{t})$, and their
derivatives,
$\mathbf{G}^{\prime}_{1}(\widetilde{\mathbf{R}}_{t},\widetilde{\gamma}_{t})$
and
$\mathbf{G}^{\prime}_{2}(\widetilde{\mathbf{R}}_{t},\widetilde{\gamma}_{t})$,
are obtained according to (20) and (21). Lastly, the constant transmitted
vector by the BS, $\mathbf{v}_{\textnormal{b}}$, is set to the right singular
vector of the matrix $\mathbf{H}_{\textnormal{b-s}}$ that corresponds to the
largest eigenvalue.

SeparableMAPProjectorLMAPEstimatorextext$\widehat{\bm{\Upsilon}}$$\mathbf{X}$$\mathbf{Z}$$\widetilde{\mathbf{R}},\widetilde{\gamma}$$\bar{\gamma}$$\bar{{\bm{\Upsilon}}}$$\mathbf{R},\gamma$$\widehat{{\bm{\Upsilon}}}$$\widehat{\gamma}$Precoder$\&$ScalarOptimizer$\widehat{\mathbf{F}}$$\widehat{\alpha}_{\textnormal{b}}$,
$\widehat{\alpha}_{\textnormal{s}}$

Figure 2: Block diagram of the proposed algorithm. The calculation of
extrinsic information is performed by the “ext” blocks.

###  IV-B Optimal Precoding and Scaling Factors

The receive scaling factor, $\alpha_{\textnormal{s}}$, is decoupled from the
other optimization variables, i.e., $\mathbf{F}$ and
$\alpha_{\textnormal{b}}$, in the objective function (23a). Therefore, it can
be optimized independently of the other two variables as follows:

$\displaystyle\begin{split}\displaystyle\operatorname*{arg\,min}_{\alpha_{\textnormal{s}}}\quad&\bigg{\|}\alpha_{\textnormal{s}}\sqrt{P_{\textnormal{s}}}\Big{[}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{1})\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}},\cdots,\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{L})\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}}\Big{]}-\left[\mathbf{0}_{L\times
B}\
\mathbf{S}_{\textnormal{s}}^{\mathsf{T}}\right]^{\mathsf{T}}\bigg{\|}_{\textnormal{F}}^{2}\\\
&+LR\sigma_{\textnormal{w}}^{2}\alpha_{\textnormal{s}}^{2}.\end{split}$ (42)

By defining the matrix:

$\mathbf{C}\triangleq\sqrt{P_{\textnormal{s}}}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}}){\bm{\Upsilon}},$
(43)

we rewrite (42) as:

$\displaystyle\operatorname*{arg\,min}_{\alpha_{\textnormal{s}}}\quad\left\lVert\alpha_{\textnormal{s}}\mathbf{C}-\left[\mathbf{0}_{L\times
B}\
\mathbf{S}_{\textnormal{s}}^{\mathsf{T}}\right]^{\mathsf{T}}\right\rVert_{\textnormal{F}}^{2}+LR\sigma_{\textnormal{w}}^{2}\alpha_{\textnormal{s}}^{2}.$
(44)

From (44), we establish the closed-form solution to the optimization problem
in (42) as follows:

$\alpha_{\textnormal{s}}^{\textnormal{opt}}=g_{4}(\mathbf{C})~{}\triangleq~{}\frac{\textnormal{Tr}\left(\mathbf{C}^{\mathsf{H}}\left[\mathbf{0}_{L\times
B}\
\mathbf{S}_{\textnormal{s}}^{\mathsf{T}}\right]^{\mathsf{T}}+\left[\mathbf{0}_{L\times
B}\
\mathbf{S}_{\textnormal{s}}^{\mathsf{H}}\right]\mathbf{C}\right)}{2\left(\left\lVert\mathbf{C}\right\rVert_{\textnormal{F}}^{2}+LR\sigma_{\textnormal{w}}^{2}\right)}.$
(45)

Now, the optimal precoding matrix, $\mathbf{F}$, and receive scaling factor,
$\alpha_{\textnormal{b}}$, are obtained as a solution to the following
optimization problem:

$\displaystyle\begin{split}\displaystyle\operatorname*{arg\,min}_{\alpha_{\textnormal{b}},\mathbf{F}}\quad&\bigg{\|}\alpha_{\textnormal{b}}\Big{[}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{1})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F},\cdots,\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{L})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\Big{]}\\\
&-\Big{[}[\mathbf{I}_{B,1},\ \mathbf{0}_{B\times
R,1}]^{\mathsf{T}},\cdots,[\mathbf{I}_{B,L},\ \mathbf{0}_{B\times
R,L}]^{\mathsf{T}}\Big{]}\bigg{\|}_{\textnormal{F}}^{2}+LB\sigma_{\textnormal{w}}^{2}\alpha_{\textnormal{b}}^{2}\end{split}$
(46a) s.t.
$\displaystyle\left\lVert\mathbf{F}\right\rVert_{\textnormal{F}}^{2}=P_{\textnormal{b}}.$
(46b)

By defining the matrices:

$\displaystyle\mathbf{K}$
$\displaystyle~{}\triangleq~{}\displaystyle\sum_{l=1}^{L}\left(\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}({\bm{\upsilon}}_{l})\mathbf{H}_{\textnormal{b-s}}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\right)^{\mathsf{H}}\left(\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}({\bm{\upsilon}}_{l})\mathbf{H}_{\textnormal{b-s}}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\right)$
(47) $\displaystyle\mathbf{E}$ $\displaystyle~{}\triangleq~{}[\mathbf{I}_{B}\
\mathbf{0}_{B\times
R}]\displaystyle\sum_{l=1}^{L}\left(\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}({\bm{\upsilon}}_{l})\mathbf{H}_{\textnormal{b-s}}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\right),$
(48)

the optimization problem in (46) becomes a constrained MMSE transmit precoding
optimization for MIMO systems. The problem can be solved jointly by Lagrange
optimization. We construct the Lagrangian function for the optimization
problem in (46) as follows:

$\operatorname{\mathcal{L}}(\mathbf{F},\alpha_{\textnormal{b}},\lambda)=\textnormal{Tr}\left(\alpha_{\textnormal{b}}^{2}\mathbf{K}\mathbf{F}\mathbf{F}^{\mathsf{H}}-\alpha_{\textnormal{b}}\mathbf{E}\mathbf{F}-\alpha_{\textnormal{b}}\mathbf{F}^{\mathsf{H}}\mathbf{E}^{\mathsf{H}}\right)+LB+LB\alpha_{\textnormal{b}}^{2}\sigma_{\textnormal{w}}^{2}+\lambda\left(\textnormal{Tr}\left(\mathbf{F}\mathbf{F}^{\mathsf{H}}\right)-P_{\textnormal{b}}\right),$
(49)

with $\lambda\in\operatorname{\mathbb{R}}$ being the Lagrange multiplier. The
closed-form solutions for optimal $\alpha_{\textnormal{b}}$ and $\mathbf{F}$
are given below and we refer the reader to [22] for more details:

$\displaystyle\alpha_{\textnormal{b}}^{\textnormal{opt}}$
$\displaystyle~{}=~{}g_{5}\left(\mathbf{K},~{}\mathbf{E}\right)~{}\triangleq~{}\sqrt{\frac{1}{P_{\textnormal{b}}}}\sqrt{\textnormal{Tr}\left(\left[\mathbf{K}+\frac{LB\sigma_{\textnormal{w}}^{2}\mathbf{I}_{N}}{P_{\textnormal{b}}}\right]^{-2}\mathbf{E}^{\mathsf{H}}\mathbf{E}\right)},$
(50) $\displaystyle\mathbf{F}^{\textnormal{opt}}$
$\displaystyle~{}=~{}g_{6}\left(\mathbf{K},~{}\mathbf{E}\right)~{}\triangleq~{}\frac{\sqrt{P_{\textnormal{b}}}\left[\mathbf{K}+\frac{LB\sigma_{\textnormal{w}}^{2}\mathbf{I}_{N}}{P_{\textnormal{b}}}\right]^{-1}\mathbf{E}^{\mathsf{H}}}{\sqrt{\textnormal{Tr}\left(\left[\mathbf{K}+\frac{LB\sigma_{\textnormal{w}}^{2}\mathbf{I}_{N}}{P_{\textnormal{b}}}\right]^{-2}\mathbf{E}^{\mathsf{H}}\mathbf{E}\right)}}.$
(51)

Algorithm 2 OOVAMP-based joint optimization algorithm

Given $\mathbf{H}_{\textnormal{s-u}}$, $\mathbf{H}_{\textnormal{b-u}}$,
$\mathbf{H}_{\textnormal{b-s}}$, $\mathbf{S}_{\textnormal{s}}$, a precision
tolerance $(\epsilon)$, and a maximum number of iterations
$(T_{\textnormal{MAX}})$

1: Initialize $\widehat{{\bm{\Upsilon}}}_{\textnormal{0}},\ \mathbf{R}_{0}$,
$\gamma_{0}\geq 0$ and $t\leftarrow 1$, and obtain
$\mathbf{v}_{\textnormal{b}}$ from $\mathbf{H}_{\textnormal{b-s}}$

2: Compute $\widehat{\mathbf{C}}_{0},~{}\widehat{\mathbf{K}}_{0}$ and
$\widehat{\mathbf{E}}_{0}$ by substituting
$\widehat{{\bm{\Upsilon}}}_{\textnormal{0}}$ into (43), (47) and (48).

3:
$\widehat{\alpha}_{\textnormal{s},0}=g_{4}\left(\widehat{\mathbf{C}}_{0}\right)$

4:
$\widehat{\alpha}_{\textnormal{b},0}=g_{5}\left(\widehat{\mathbf{K}}_{0},~{}\widehat{\mathbf{E}}_{0}\right)$

5:
$\widehat{\mathbf{F}}_{\textnormal{0}}=g_{6}\left(\widehat{\mathbf{K}}_{0},~{}\widehat{\mathbf{E}}_{0}\right)$

6: repeat

7: // LMAP Estimator

8: Compute matrices
$\mathbf{A},~{}\mathbf{B},~{}\mathbf{D},~{}\mathbf{M},~{}\mathbf{X}$ and
$\mathbf{Z}$ by substituting
$\widehat{\alpha}_{\textnormal{s},t-1},~{}\widehat{\alpha}_{\textnormal{b},t-1}$
and $\widehat{\mathbf{F}}_{t-1}$ into (25) to (30).

9:
$\bar{\bm{\Upsilon}}_{t}=\left(\mathbf{D}^{\mathsf{H}}\mathbf{D}+\mathbf{M}^{\mathsf{H}}\mathbf{M}+\gamma_{t-1}\mathbf{I}_{K}\right)^{-1}\left(\mathbf{D}^{\mathsf{H}}\mathbf{X}+\mathbf{M}^{\mathsf{H}}\mathbf{Z}+\gamma_{t-1}\mathbf{R}_{t-1}\right)$

10:
$\bar{\gamma}_{t}=K\textnormal{Tr}\left(\left[\mathbf{D}^{\mathsf{H}}\mathbf{D}+\mathbf{M}^{\mathsf{H}}\mathbf{M}+\gamma_{t-1}\mathbf{I}_{K}\right]^{-1}\right)^{-1}$

11: $\widetilde{\gamma}_{t}=\bar{\gamma}_{t}-\gamma_{t-1}$

12:
$\widetilde{\mathbf{R}}_{t}=\widetilde{\gamma}_{t}^{-1}\left(\bar{\bm{\Upsilon}}_{t}\bar{\gamma}_{t}-\mathbf{R}_{t-1}\gamma_{t-1}\right)$

13: // Separable MAP Projector

14:
$\widehat{{\bm{\Upsilon}}}_{t}=\mathbf{G}_{1}\left(\widetilde{\mathbf{R}}_{t}\right)$

15:
$\widehat{\gamma}_{t}=\widetilde{\gamma}_{t}^{-1}\left\langle\mathbf{G}^{\prime}_{1}\left(\widetilde{\mathbf{R}}_{t}\right)\right\rangle$.

16: $\gamma_{t}=\widehat{\gamma}_{t}-\widetilde{\gamma}_{t}$

17:
$\mathbf{R}_{t}=\gamma_{t}^{-1}\left(\widehat{\gamma}_{t}\widehat{{\bm{\Upsilon}}}_{t}-\widetilde{\gamma}_{t}\widetilde{\mathbf{R}}_{t}\right)$

18: //Find $\alpha_{\textnormal{s}},~{}\alpha_{\textnormal{b}}$ and
$\mathbf{F}$ through their closed-form solutions.

19: Compute $\widehat{\mathbf{C}}_{t},~{}\widehat{\mathbf{K}}_{t}$ and
$\widehat{\mathbf{E}}_{t}$ by substituting
$\widehat{{\bm{\Upsilon}}}_{\textnormal{t}}$ into (43), (47) and (48).

20:
$\widehat{\alpha}_{\textnormal{s},t}=g_{4}\left(\widehat{\mathbf{C}}_{t}\right)$

21:
$\widehat{\alpha}_{\textnormal{b},t}=g_{5}\left(\widehat{\mathbf{K}}_{t},~{}\widehat{\mathbf{E}}_{t}\right)$

22:
$\widehat{\mathbf{F}}_{\textnormal{t}}=g_{6}\left(\widehat{\mathbf{K}}_{t},~{}\widehat{\mathbf{E}}_{t}\right)$

23: $t\leftarrow t+1$

24: until $|E_{t}-E_{t-1}|<\epsilon E_{t-1}$ or $t>T_{\textnormal{MAX}}$

25: return $\widehat{{\bm{\upsilon}}}_{t},\ \widehat{\mathbf{F}}_{t},\
\widehat{\alpha}_{t}$.

Since we have now found the solution to both optimization problems in (22) and
(23), we can combine their solutions together into one algorithm. We define
the MSE at iteration $t$ as follows:

$E_{t}~{}\triangleq~{}f\big{(}\widehat{\alpha}_{\textnormal{s},t},~{}\widehat{\alpha}_{\textnormal{b},t},~{}\widehat{\mathbf{F}}_{t},~{}\widehat{\bm{\Upsilon}}_{t}\big{)}.$
(52)

The algorithm stops when $|E_{t}-E_{t-1}|<\epsilon E_{t-1}$, where
$\epsilon\in\operatorname{\mathbb{R}}_{+}$ is some precision tolerance. The
overall block diagram and the algorithmic steps are shown, respectively, in
Fig. 2 and Algorithm 2. The convergence and complexity of the OOVAMP-based
approach are discussed in [8]. Because of the monotone convergence theorem in
real analysis [27], Algorithm 2 is guaranteed to converge since the MSE is
minimized in every step and the objective function,
$f(\alpha_{\textnormal{s}},\alpha_{\textnormal{b}},\mathbf{F},{\mathbf{\Upsilon}})$,
is lower bounded by zero. The algorithm can be efficiently implemented by
exploiting matrix structures, and by using the singular value decomposition
(SVD) form of OOVAMP so that the computational complexity is of
$\operatorname{\mathcal{O}}(MNL(K+N))$.

##  V Numerical Results and Performance Analysis

###  V-A Simulation Model and Parameters

We present Monte-Carlo simulation results to assess the performance of the
proposed MIS-based communication scheme. We assume that the MIS is located at
a fixed distance of $500$ m from the BS and the users are spread uniformly at
a radial distance of $10$ m to $50$ m from the MIS. A path-based propagation
channel model, also known as parametric channel model [24], is used. Such a
model is more appropriate for systems with large antenna arrays. One key
parameter of such a channel model is the number of multi-path components of
the BS-MIS channel which governs the effect of channel correlation.

TABLE I: Simulation parameters with their notations and values Parameter |  Notation, Value |  Parameter |  Notation, Value
---|---|---|---
Number of channel paths in the MIS-BS link |  $Q_{\textnormal{MIS}}=10$ |  MIS-BS distance |  $d_{\textnormal{MIS}}=500$ m
Number of channel paths in the BS-user link |  $Q_{\textnormal{b-u}}=2$ |  User-BS distance |  $d=500$ m
Number of channel paths in the MIS-user link |  $Q_{\textnormal{s-u}}=2$ |  User-MIS distance |  $d^{\prime}\in[10,50]$ m
Path-loss exponent for the MIS-BS, MIS-user link |  $\eta=2.5$ |  Noise variance |  $\sigma_{\textnormal{w}}^{2}=-100$ dBm
Path-loss exponent for the BS-user link |  $\eta=3.7$ |  Channel path gains distribution |  $c_{q}\sim\operatorname{\mathcal{C}}\operatorname{\mathcal{N}}(0,1)$
Reference distance |  $d_{\textnormal{0}}=1$ m |  Path-loss at the reference distance |  $C_{\textnormal{0}}=-30$ dB

Assuming a uniform linear array with $N$ antennas at the BS and a square
uniform planar array with $K$ antenna elements at the MIS, the channel between
the MIS and the BS is generated according to:

$\mathbf{H}_{\textnormal{b-s}}~{}=~{}\sqrt{L(d_{\textnormal{MIS}})}\displaystyle\sum_{q=1}^{Q_{\textnormal{MIS}}}\mathsf{c}_{q}\mathbf{a}_{\textnormal{MIS}}(\varphi_{q},\psi_{q})\mathbf{a}_{\textnormal{BS}}(\phi_{q})^{\mathsf{T}}.$
(53)

The channel vectors for the link between each single antenna $m$-th user and
the MIS, and each $m$-th user and the BS are modeled, respectively, as
follows:

$\displaystyle\mathbf{h}_{\textnormal{s-u},m}$
$\displaystyle~{}=~{}\sqrt{L(d^{\prime}_{m})}\displaystyle\sum_{q=1}^{Q_{\textnormal{s-u}}}\mathsf{c}_{m,q}\mathbf{a}_{\textnormal{MIS}}(\varphi_{m,q},\psi_{m,q}),\quad
m=1,\cdots,M,$ (54) $\displaystyle\mathbf{h}_{\textnormal{b-u},m}$
$\displaystyle~{}=~{}\sqrt{L(d_{m})}\displaystyle\sum_{q=1}^{Q_{\textnormal{b-u}}}\mathsf{c}_{m,q}\mathbf{a}_{\textnormal{BS}}(\phi_{m,q}),\quad
m=1,\cdots,M.$ (55)

The term $L(d)=C_{\textnormal{0}}(d/d_{\textnormal{0}})^{\eta}$ in (53), (55),
(54) is the distance-dependent path-loss factor, where $C_{\textnormal{0}}$
denotes the path-loss at a reference distance $d_{\textnormal{0}}=1$ m, and
$\eta$ is the path-loss exponent. To account for channel correlation effects,
we set the number of multi-path components lower than the number of BS
antennas and the MIS antenna elements for the BS-MIS channel, thereby making
the channel matrix rank-deficient. Therefore, in the simulations, we set the
number of BS-served users lower than the rank of the BS-MIS channel matrix
$\mathbf{H}_{\textnormal{b-s}}$. For brevity, the details on the array
steering vectors, $\mathbf{a}_{\textnormal{MIS}}$ and
$\mathbf{a}_{\textnormal{BS}}$, are omitted, and interested readers are
referred to [8] for complete information on the channel model. Table I
describes the various simulation parameters along with their values used in
the simulations. Moreover, we split the total transmit power, $P$, between the
BS- and MIS-served users according to the share of each type of users, i.e.,
$P_{\textnormal{b}}=\frac{B}{M}P$ and $P_{\textnormal{s}}=\frac{R}{M}P$. The
results are averaged over $1000$ independent Monte-Carlo trials.

We use the sum-rate, $\widehat{C}$, for performance evaluation which is
defined as follows:

$\widehat{C}~{}=~{}\displaystyle\sum_{m=1}^{M}\log_{2}\left(\frac{1}{\textnormal{MMSE}_{m}}\right),$
(56)

where $\textnormal{MMSE}_{m}$ refers to the MSE of each $m$-th user’s received
symbol. Since the proposed approach of using the RIS as a modulating surface
(i.e., MIS) is novel, we benchmark it against the following two traditional
schemes wherein the RIS is merely used for passive beamforming purposes in
combination with MMSE precoding at the BS:

  * i.

Scheme 1: a multi-user MIMO system assisted by one RIS where the optimization
of the phase matrix is solved through alternating optimization and OOVAMP.

  * ii.

Scheme 2: a multi-user MIMO system assisted by one RIS where the optimization
of the phase matrix is performed using the semidefinite relaxation (SDR)
technique.

###  V-B Performance Results With Perfect CSI

We consider a typical urban or suburban environment where the BS is located
faraway from the users and has no LOS to them. However, the MIS is installed
at a location where an LOS component is present in the BS-MIS link but not in
the user-MIS link. We also set the number of BS antennas to $N=32$ and the
number of MIS reflective elements to $K=256$. Fig. 3(a), depicts the
achievable sum-rate versus the transmit power, $P$, for the different
considered transmission schemes. The total number of users is equal to $M=8$.
Here we consider two configurations for the proposed approach: $i)$ when all
the users are served by the MIS, and $ii)$ a hybrid case wherein the MIS and
the BS both serve $4$ users each. The configuration in which all the users are
served by the MIS significantly outperforms the scheme in which the RIS is
only used for beamforming (i.e., Scheme 1). This is because, the MIS-served
users do not suffer from the path-loss between the BS and the MIS. Moreover,
at low transmit power (e.g., $P=20$ dBm), the sum-rate for the users solely
served by the MIS is, respectively, two and four times the sum-rate of the BS-
served users for the beamforming-only RIS-assisted OOVAMP-based and SDR-based
approaches. For the hybrid configuration, the resulting sum-rate edges Scheme
1 but, it is lower than the case when all users are served by the MIS. This
confirms that it is more beneficial to serve the users by the MIS.

(a) $M=8,~{}N=32$, and
$K=256$.

(b) $M=8,\ N=32$, and
$P=25$ dBm.

Figure 3: (a) Sum-rate versus transmit power, and (b) Sum-rate versus the
number of MIS-served users.

Fig. 3(b) shows a plot of the sum-rate against the share of MIS-served users
among the total number of users. It is observed that the combined sum-rate
first decreases and then monotonically increases with the number of MIS-served
users. This is because of the presence of cross-user interference among the
MIS-served and the BS-served users since the MIS is performing both tasks,
i.e., beamforming to assist the BS and also data embedding to serve another
set of users. This implies that there is more loss than gain when the ratio of
the MIS-served users to the total number of users becomes small.

Fig. 4(a) illustrates the sum-rate versus the total number of users being
served. The users are solely served by the BS for one plot and by the MIS for
the other. Here we show the benefits of the approach of using the MIS as a
modulating surface to directly serve users. Although the sum-rate of the BS-
served users is higher when the number of users becomes small, the sum-rate of
the MIS-served users keeps increasing with the number of users while the sum-
rate of the BS-served users only increases up to the number of available
channel paths which is set to $10$. This is because, the number of users that
can be served by the MIS is independent of the number of BS antennas and the
correlation in the MIS-BS channel. The upper limit for the number of MIS-
served users is equal to the number of MIS antenna elements $K$.

(a) $N=32,~{}K=256$ and
$P=30$ dBm.

(b) $M=8,~{}N=32$, and
$K=256$.

Figure 4: (a) Sum-rate versus the number of users, and (b) Sum-rate versus
transmit power under imperfect CSI.

###  V-C Performance Results With Imperfect CSI

In this section, we asses the performance of the proposed scheme in the
presence of channel estimation errors. Specifically, we consider a scenario
where pilot training followed by MMSE estimation algorithms are used to
estimate the cascaded BS-MIS-users and the direct BS-users channels [28, 29].
We model the estimated channel matrix and vectors using the statistical CSI
error model proposed in [30, 31, 32] as follows:

$\widehat{\mathbf{H}}_{\textnormal{b-s}}=\kappa\mathbf{H}_{\textnormal{b-s}}+\sqrt{(1-\kappa^{2})L(d_{\textnormal{MIS}})}{\bm{\Delta}}_{\textnormal{b-s}},$
(57)
$\widehat{\mathbf{h}}_{\textnormal{b-u},m}=\kappa\mathbf{h}_{\textnormal{b-u},m}+\sqrt{(1-\kappa^{2})L(d_{m})}{\bm{\delta}}_{\textnormal{b-u},m},\quad
m=1,\cdots,M,$ (58)
$\widehat{\mathbf{h}}_{\textnormal{s-u},m}=\kappa\mathbf{h}_{\textnormal{s-u},m}+\sqrt{(1-\kappa^{2})L(d^{\prime}_{m})}{\bm{\delta}}_{\textnormal{s-u},m},\quad
m=1,\cdots,M,$ (59)

where $\kappa\in[0,1]$ denotes the channel estimation accuracy and
${\bm{\Delta}}_{\textnormal{b-s}}$, ${\bm{\delta}}_{\textnormal{b-u},m}$ and
${\bm{\delta}}_{\textnormal{s-u},m}$ follow the circularly symmetric complex
Gaussian (CSCG) distribution, i.e.,
$\textnormal{vec}({\bm{\Delta}}_{\textnormal{b-s}})\sim\operatorname{\mathcal{C}}\operatorname{\mathcal{N}}(\mathbf{0},\mathbf{1}_{N\times
N}\otimes\mathbf{I}_{K})$,
${\bm{\delta}}_{\textnormal{b-u},m}\sim\operatorname{\mathcal{C}}\operatorname{\mathcal{N}}(\mathbf{0},\mathbf{I}_{N})$
and
${\bm{\delta}}_{\textnormal{s-u},m}\sim\operatorname{\mathcal{C}}\operatorname{\mathcal{N}}(\mathbf{0},\mathbf{I}_{K})$.
We first optimize the variables
$\alpha_{\textnormal{s}},~{}\alpha_{\textnormal{b}},~{}\mathbf{F}$ and
${\bm{\Upsilon}}$ under imperfect CSI and then use the exact CSI matrices to
calculate the sum rate.

Fig. 4(b) plots the sum-rate versus transmit power for different values of the
channel estimation accuracy parameter $\kappa$. We also include plots for the
other schemes (i.e., Schemes 1 and 2) which use the RIS as a purely reflective
surface) under perfect CSI for reference. The results show the resilience of
the proposed MIS-based approach against small channel estimation errors. At
low SNR, it is observed that the proposed design with a low channel estimation
accuracy of $\kappa=0.90$ performs better than the baseline schemes even under
perfect CSIs. Moreover, the performance loss is negligible with a high channel
estimation accuracy value of $\kappa=0.99$.

##  VI Conclusion

We have presented a novel approach of employing the MIS for performing passive
beamforming and data embedding for the BS-served and the MIS-served users,
respectively, in a multi-user MIS-assisted downlink MIMO system. The
associated joint convex optimization problem has been formulated under the sum
MMSE criterion in order to maximize the users’ spectral efficiency.
Alternating minimization has been used to split the original optimization
problem into two tasks, i.e., separately optimizing the MIS phase shifts and
jointly optimizing the BS precoding and the receive scaling factors for the
BS- and MIS-served users. The optimal solution to the joint optimization
problem for precoding is found in closed form. We have optimized the MIS phase
shifts using OOVAMP by deriving the problem-specific modules of the OOVAMP
algorithm. Moreover, the original joint problem has been solved under both the
ideal and a practical constraint on the MIS phase shifts. Simulation results
illustrate highly superior system throughput performance of the proposed MIS-
based scheme over two baseline schemes in which the MIS is used for
beamforming purposes only. Moreover, the proposed approach can support more
number of users simultaneously than existing beamforming approaches wherein
the users are served by the BS only. Finally, the results under imperfect CSI
confirm that the performance remains nearly unchanged in the presence of small
channel estimation errors.

Three possible extensions to the proposed framework for MIS-assisted networks
are outlined as follows:

  * •

Optimal power allocation: In the present work, the power split between the BS-
and MIS-served users was not optimized. In fact, we simply split the power
between the BS- and MIS-served users according to the share of each type of
users. To utilize the transmit power efficiently, the joint optimization
problem can be extended to optimally allocate power between the BS- and MIS-
served users.

  * •

Optimal user assignment: Instead of splitting the users arbitrarily between
the BS and the MIS, the decision must be taken on the basis of CSI. The users
for which the direct BS-user link is strong enough to provide higher data rate
than if they were served by the MIS, can be served and beamformed by the BS
and the MIS, respectively. Similarly, the users that have a weak direct BS-
user link can be served through modulation by the MIS. User assignment can be
optimized to maximize system throughput.

  * •

Jointly optimized reference carrier signal: Presently, we assign the strongest
path in the BS-MIS link to transmit the reference carrier signal pertaining to
the MIS-served users. The remaining paths are used to serve the BS-served
users which results in lower sum-rate for BS-served users. To mitigate this
problem, the reference transmitted vector, $\mathbf{v}_{\textnormal{b}}$, can
also be jointly optimized with the other variables to enhance the overall
throughput of the system.

## Appendix A

We have the following function:

$f(\alpha_{\textnormal{s}},\alpha_{\textnormal{b}},\mathbf{F},{\mathbf{\Upsilon}})=\operatorname{\mathbb{E}}_{\bm{\mathsf{Y}},\bm{\mathsf{S}}}\left\\{\left\lVert\bm{\mathsf{Y}}-\bm{\mathsf{S}}\right\rVert_{\textnormal{F}}^{2}\right\\},$
(60)

where

$\displaystyle\begin{split}\bm{\mathsf{Y}}&=\alpha_{\textnormal{b}}\Big{[}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{1})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}\bm{\mathsf{s}}_{\textnormal{b},1}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\bm{\mathsf{s}}_{\textnormal{b},1},\cdots,\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{L})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}\bm{\mathsf{s}}_{\textnormal{b},L}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\bm{\mathsf{s}}_{\textnormal{b},L}\Big{]}\\\
&+\alpha_{\textnormal{s}}\sqrt{P_{\textnormal{s}}}\Big{[}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{1})\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}},\cdots,\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{L})\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}}\Big{]}\\\
&+\Big{[}\alpha_{\textnormal{b}}\bm{\mathsf{W}}_{\textnormal{b}}^{\mathsf{T}},~{}\alpha_{\textnormal{s}}\bm{\mathsf{W}}_{\textnormal{s}}^{\mathsf{T}}\Big{]}^{\mathsf{T}},\end{split}$
(61) $\displaystyle\bm{\mathsf{S}}$
$\displaystyle=\big{[}\bm{\mathsf{S}}_{\textnormal{b}}^{\mathsf{T}},~{}\mathbf{S}_{\textnormal{s}}^{\mathsf{T}}\big{]}^{\mathsf{T}}.$
(62)

We explicitly write the function in (60) and then expand it as follows:

$\displaystyle
f(\alpha_{\textnormal{s}},\alpha_{\textnormal{b}},\mathbf{F},{\mathbf{\Upsilon}})$
$\displaystyle=\operatorname{\mathbb{E}}_{\bm{\mathsf{Y}},\bm{\mathsf{S}}}\left\\{\textnormal{Tr}\left(\bm{\mathsf{Y}}^{\mathsf{H}}\bm{\mathsf{Y}}-\bm{\mathsf{Y}}^{\mathsf{H}}\bm{\mathsf{S}}-\bm{\mathsf{S}}^{\mathsf{H}}\bm{\mathsf{Y}}+\bm{\mathsf{S}}^{\mathsf{H}}\bm{\mathsf{S}}\right)\right\\}$
(63)
$\displaystyle\begin{split}&=\operatorname{\mathbb{E}}_{\bm{\mathsf{Y}},\bm{\mathsf{S}}_{\textnormal{b}}}\Big{\\{}\textnormal{Tr}\Big{(}\bm{\mathsf{Y}}^{\mathsf{H}}\bm{\mathsf{Y}}-\bm{\mathsf{Y}}^{\mathsf{H}}\big{[}\bm{\mathsf{S}}_{\textnormal{b}}^{\mathsf{T}},~{}\mathbf{0}_{L\times
R}\big{]}^{\mathsf{T}}-\bm{\mathsf{Y}}^{\mathsf{H}}\big{[}\mathbf{0}_{L\times
B},~{}\mathbf{S}_{\textnormal{s}}^{\mathsf{T}}\big{]}^{\mathsf{T}}-\big{[}\bm{\mathsf{S}}_{\textnormal{b}}^{\mathsf{H}},~{}\mathbf{0}_{L\times
R}\big{]}\bm{\mathsf{Y}}\\\ &-\big{[}\mathbf{0}_{L\times
B},~{}\mathbf{S}_{\textnormal{s}}^{\mathsf{H}}\big{]}\bm{\mathsf{Y}}+\big{[}\bm{\mathsf{S}}_{\textnormal{b}}^{\mathsf{H}},~{}\mathbf{0}_{L\times
R}\big{]}\big{[}\bm{\mathsf{S}}_{\textnormal{b}}^{\mathsf{T}},~{}\mathbf{0}_{L\times
R}\big{]}^{\mathsf{T}}+\big{[}\bm{\mathsf{S}}_{\textnormal{b}}^{\mathsf{H}},~{}\mathbf{0}_{L\times
R}\big{]}\big{[}\mathbf{0}_{L\times
B},~{}\mathbf{S}_{\textnormal{s}}^{\mathsf{T}}\big{]}^{\mathsf{T}}\\\
&+\big{[}\mathbf{0}_{L\times
B},~{}\mathbf{S}_{\textnormal{s}}^{\mathsf{H}}\big{]}\big{[}\bm{\mathsf{S}}_{\textnormal{b}}^{\mathsf{T}},~{}\mathbf{0}_{L\times
R}\big{]}^{\mathsf{T}}+\big{[}\mathbf{0}_{L\times
B},~{}\mathbf{S}_{\textnormal{s}}^{\mathsf{H}}\big{]}\big{[}\mathbf{0}_{L\times
B},~{}\mathbf{S}_{\textnormal{s}}^{\mathsf{T}}\big{]}^{\mathsf{T}}\Big{)}\Big{\\}}\end{split}$
$\displaystyle\begin{split}f(\alpha_{\textnormal{s}},\alpha_{\textnormal{b}},\mathbf{F},{\mathbf{\Upsilon}})&=\operatorname{\mathbb{E}}_{\bm{\mathsf{Y}},\bm{\mathsf{S}}_{\textnormal{b}}}\Big{\\{}\textnormal{Tr}\Big{(}\bm{\mathsf{Y}}^{\mathsf{H}}\bm{\mathsf{Y}}-\bm{\mathsf{Y}}^{\mathsf{H}}\big{[}\bm{\mathsf{S}}_{\textnormal{b}}^{\mathsf{T}},~{}\mathbf{0}_{L\times
R}\big{]}^{\mathsf{T}}-\big{[}\bm{\mathsf{S}}_{\textnormal{b}}^{\mathsf{H}},~{}\mathbf{0}_{L\times
R}\big{]}\bm{\mathsf{Y}}\\\
&-\bm{\mathsf{Y}}^{\mathsf{H}}\big{[}\mathbf{0}_{L\times
B},~{}\mathbf{S}_{\textnormal{s}}^{\mathsf{T}}\big{]}^{\mathsf{T}}-\big{[}\mathbf{0}_{L\times
B},~{}\mathbf{S}_{\textnormal{s}}^{\mathsf{H}}\big{]}\bm{\mathsf{Y}}+\big{[}\bm{\mathsf{S}}_{\textnormal{b}}^{\mathsf{H}},~{}\mathbf{0}_{L\times
R}\big{]}\big{[}\bm{\mathsf{S}}_{\textnormal{b}}^{\mathsf{T}},~{}\mathbf{0}_{L\times
R}\big{]}^{\mathsf{T}}\\\ &+\big{[}\mathbf{0}_{L\times
B},~{}\mathbf{S}_{\textnormal{s}}^{\mathsf{H}}\big{]}\big{[}\mathbf{0}_{L\times
B},~{}\mathbf{S}_{\textnormal{s}}^{\mathsf{T}}\big{]}^{\mathsf{T}}\Big{)}\Big{\\}}.\end{split}$
(64)

By defining the matrices:

$\displaystyle\bm{\mathsf{Y}}_{\textnormal{b}}$
$\displaystyle=\alpha_{\textnormal{b}}\Big{[}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{1})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}\bm{\mathsf{s}}_{\textnormal{b},1}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\bm{\mathsf{s}}_{\textnormal{b},1},\cdots,\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{L})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}\bm{\mathsf{s}}_{\textnormal{b},L}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\bm{\mathsf{s}}_{\textnormal{b},L}\Big{]},$
(65) $\displaystyle\mathbf{Y}_{\textnormal{e}}$
$\displaystyle=\alpha_{\textnormal{b}}\Big{[}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{1})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F},\cdots,\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{L})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\Big{]},$
(66) $\displaystyle\mathbf{Y}_{\textnormal{s}}$
$\displaystyle=\alpha_{\textnormal{s}}\sqrt{P_{\textnormal{s}}}\Big{[}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{1})\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}},\cdots,\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{L})\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}}\Big{]},$
(67) $\displaystyle\bm{\mathsf{W}}$
$\displaystyle=\Big{[}\alpha_{\textnormal{b}}\bm{\mathsf{W}}_{\textnormal{b}}^{\mathsf{T}},~{}\alpha_{\textnormal{s}}\bm{\mathsf{W}}_{\textnormal{s}}^{\mathsf{T}}\Big{]}^{\mathsf{T}},$
(68)

we expand the terms in (64) and take expectation w.r.t. the random matrices
$\bm{\mathsf{S}}_{\textnormal{b}}$ and $\bm{\mathsf{W}}$ thereby leading to:

$\displaystyle\begin{split}\operatorname{\mathbb{E}}_{\bm{\mathsf{Y}},\bm{\mathsf{S}}_{\textnormal{b}}}\big{\\{}\textnormal{Tr}\big{(}\bm{\mathsf{Y}}^{\mathsf{H}}\bm{\mathsf{Y}}\big{)}\big{\\}}&=\operatorname{\mathbb{E}}_{\bm{\mathsf{W}},\bm{\mathsf{S}}_{\textnormal{b}}}\big{\\{}\textnormal{Tr}\big{(}\bm{\mathsf{Y}}_{\textnormal{b}}^{\mathsf{H}}\bm{\mathsf{Y}}_{\textnormal{b}}+\bm{\mathsf{Y}}_{\textnormal{b}}^{\mathsf{H}}\mathbf{Y}_{\textnormal{s}}+\bm{\mathsf{Y}}_{\textnormal{b}}^{\mathsf{H}}\bm{\mathsf{W}}+\mathbf{Y}_{\textnormal{s}}^{\mathsf{H}}\bm{\mathsf{Y}}_{\textnormal{b}}+\mathbf{Y}_{\textnormal{s}}^{\mathsf{H}}\mathbf{Y}_{\textnormal{s}}\\\
&+\mathbf{Y}_{\textnormal{s}}^{\mathsf{H}}\bm{\mathsf{W}}+\bm{\mathsf{W}}^{\mathsf{H}}\bm{\mathsf{Y}}_{\textnormal{b}}+\bm{\mathsf{W}}^{\mathsf{H}}\mathbf{Y}_{\textnormal{s}}+\bm{\mathsf{W}}^{\mathsf{H}}\bm{\mathsf{W}}\big{)}\big{\\}}\end{split}$
(69)
$\displaystyle=\textnormal{Tr}\big{(}\mathbf{Y}_{\textnormal{e}}^{\mathsf{H}}\mathbf{Y}_{\textnormal{e}}+\mathbf{Y}_{\textnormal{s}}^{\mathsf{H}}\mathbf{Y}_{\textnormal{s}}\big{)}+LB\sigma_{\textnormal{w}}^{2}\alpha_{\textnormal{b}}^{2}+LR\sigma_{\textnormal{w}}^{2}\alpha_{\textnormal{s}}^{2},$
$\displaystyle\begin{split}\operatorname{\mathbb{E}}_{\bm{\mathsf{Y}},\bm{\mathsf{S}}_{\textnormal{b}}}\Big{\\{}\textnormal{Tr}\Big{(}\big{[}\bm{\mathsf{S}}_{\textnormal{b}}^{\mathsf{H}},~{}\mathbf{0}_{L\times
R}\big{]}\bm{\mathsf{Y}}\Big{)}\Big{\\}}&=\operatorname{\mathbb{E}}_{\bm{\mathsf{W}},\bm{\mathsf{S}}_{\textnormal{b}}}\Big{\\{}\textnormal{Tr}\Big{(}\big{[}\bm{\mathsf{S}}_{\textnormal{b}}^{\mathsf{H}},~{}\mathbf{0}_{L\times
R}\big{]}\bm{\mathsf{Y}}_{\textnormal{b}}+\big{[}\bm{\mathsf{S}}_{\textnormal{b}}^{\mathsf{H}},~{}\mathbf{0}_{L\times
R}\big{]}\mathbf{Y}_{\textnormal{s}}\\\
&+\big{[}\bm{\mathsf{S}}_{\textnormal{b}}^{\mathsf{H}},~{}\mathbf{0}_{L\times
R}\big{]}\bm{\mathsf{W}}\Big{)}\Big{\\}}\end{split}$ (70)
$\displaystyle=\textnormal{Tr}\bigg{(}\Big{[}[\mathbf{I}_{B,1},\
\mathbf{0}_{B\times R,1}]^{\mathsf{T}},\cdots,[\mathbf{I}_{B,L},\
\mathbf{0}_{B\times
R,L}]^{\mathsf{T}}\Big{]}^{\mathsf{H}}\mathbf{Y}_{\textnormal{e}}\bigg{)},$
$\displaystyle\begin{split}\operatorname{\mathbb{E}}_{\bm{\mathsf{Y}},\bm{\mathsf{S}}_{\textnormal{b}}}\Big{\\{}\textnormal{Tr}\Big{(}\big{[}\mathbf{0}_{L\times
B},~{}\mathbf{S}_{\textnormal{s}}^{\mathsf{H}}\big{]}\bm{\mathsf{Y}}\Big{)}\Big{\\}}&=\operatorname{\mathbb{E}}_{\bm{\mathsf{W}},\bm{\mathsf{S}}_{\textnormal{b}}}\Big{\\{}\textnormal{Tr}\Big{(}\big{[}\mathbf{0}_{L\times
B},~{}\mathbf{S}_{\textnormal{s}}^{\mathsf{H}}\big{]}\bm{\mathsf{Y}}_{\textnormal{b}}+\big{[}\mathbf{0}_{L\times
B},~{}\mathbf{S}_{\textnormal{s}}^{\mathsf{H}}\big{]}\mathbf{Y}_{\textnormal{s}}\\\
&+\big{[}\mathbf{0}_{L\times
B},~{}\mathbf{S}_{\textnormal{s}}^{\mathsf{H}}\big{]}\bm{\mathsf{W}}\Big{)}\Big{\\}}\end{split}$
(71) $\displaystyle=\textnormal{Tr}\Big{(}\big{[}\mathbf{0}_{L\times
B},~{}\mathbf{S}_{\textnormal{s}}^{\mathsf{H}}\big{]}\mathbf{Y}_{\textnormal{s}}\Big{)},$

and

$\displaystyle\begin{split}\operatorname{\mathbb{E}}_{\bm{\mathsf{S}}_{\textnormal{b}}}\bigg{\\{}\textnormal{Tr}\Big{(}\big{[}\bm{\mathsf{S}}_{\textnormal{b}}^{\mathsf{H}},~{}\mathbf{0}_{L\times
R}\big{]}\big{[}\bm{\mathsf{S}}_{\textnormal{b}}^{\mathsf{T}},~{}\mathbf{0}_{L\times
R}\big{]}^{\mathsf{T}}\Big{)}\bigg{\\}}&=\textnormal{Tr}\bigg{(}\Big{[}[\mathbf{I}_{B,1},\
\mathbf{0}_{B\times R,1}]^{\mathsf{T}},\cdots,[\mathbf{I}_{B,L},\
\mathbf{0}_{B\times R,L}]^{\mathsf{T}}\Big{]}^{\mathsf{H}}\\\
&\times\Big{[}[\mathbf{I}_{B,1},\ \mathbf{0}_{B\times
R,1}]^{\mathsf{T}},\cdots,[\mathbf{I}_{B,L},\ \mathbf{0}_{B\times
R,L}]^{\mathsf{T}}\Big{]}\bigg{)}.\end{split}$ (72)

The remaining two terms in (64) can be computed by following (70) and (71).
Finally, by substituting (69)$-$(72) into (64) and expressing the function in
the form of norms we obtain:

$\displaystyle\begin{split}f(\alpha_{\textnormal{s}},\alpha_{\textnormal{b}},\mathbf{F},{\mathbf{\Upsilon}})&=\bigg{\|}\alpha_{\textnormal{b}}\Big{[}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{1})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F},\cdots,\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{L})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\Big{]}\\\
&-\Big{[}[\mathbf{I}_{B,1},\ \mathbf{0}_{B\times
R,1}]^{\mathsf{T}},\cdots,[\mathbf{I}_{B,L},\ \mathbf{0}_{B\times
R,L}]^{\mathsf{T}}\Big{]}\bigg{\|}_{\textnormal{F}}^{2}\\\
&+\bigg{\|}\alpha_{\textnormal{s}}\sqrt{P_{\textnormal{s}}}\Big{[}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{1})\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}},\cdots,\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{L})\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}}\Big{]}-\left[\mathbf{0}_{L\times
B}\
\mathbf{S}_{\textnormal{s}}^{\mathsf{T}}\right]^{\mathsf{T}}\bigg{\|}_{\textnormal{F}}^{2}\\\
&+LB\sigma_{\textnormal{w}}^{2}\alpha_{\textnormal{b}}^{2}+LR\sigma_{\textnormal{w}}^{2}\alpha_{\textnormal{s}}^{2}.\end{split}$
(73)

## Appendix B

We have the following two norms:

$\displaystyle\begin{split}f_{1}&~{}=~{}\bigg{\|}\alpha_{\textnormal{b}}\Big{[}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{1})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F},\cdots,\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{L})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\Big{]}\\\
&~{}-~{}\Big{[}[\mathbf{I}_{B},\ \mathbf{0}_{B\times
R,1}]^{\mathsf{T}},\cdots,[\mathbf{I}_{B},\ \mathbf{0}_{B\times
R,L}]^{\mathsf{T}}\Big{]}\bigg{\|}_{\textnormal{F}}^{2},\end{split}$ (74)
$\displaystyle f_{2}$
$\displaystyle~{}=~{}\bigg{\|}\alpha_{\textnormal{s}}\sqrt{P_{\textnormal{s}}}\Big{[}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{1})\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}},\cdots,\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\bm{\upsilon}_{L})\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}}\Big{]}-\left[\mathbf{0}_{L\times
B}\
\mathbf{S}_{\textnormal{s}}^{\mathsf{T}}\right]^{\mathsf{T}}\bigg{\|}_{\textnormal{F}}^{2}.$
(75)

By defining the matrices,
$\mathbf{A}=\alpha_{\textnormal{b}}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\in\operatorname{\mathbb{C}}^{M\times
K}$ and
$\mathbf{B}=(\mathbf{H}_{\textnormal{b-s}}\mathbf{F})^{\mathsf{T}}\in\operatorname{\mathbb{C}}^{B\times
K}$, we rewrite (74) as follows:

$\displaystyle\begin{split}f_{1}&~{}=~{}\bigg{\|}\Big{[}\mathbf{A}\textnormal{Diag}(\bm{\upsilon}_{1})\mathbf{B}^{\mathsf{T}},\cdots,\mathbf{A}\textnormal{Diag}(\bm{\upsilon}_{L})\mathbf{B}^{\mathsf{T}}\Big{]}\\\
&~{}-~{}\Big{[}[\mathbf{I}_{B},\ \mathbf{0}_{B\times
R,1}]^{\mathsf{T}}-\alpha_{\textnormal{b}}\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F},\cdots,[\mathbf{I}_{B},\
\mathbf{0}_{B\times
R,L}]^{\mathsf{T}}-\alpha_{\textnormal{b}}\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\Big{]}\bigg{\|}_{\textnormal{F}}^{2}.\end{split}$
(76)

We then define a column-wise Khatri-Rao matrix,
$\mathbf{D}\in\operatorname{\mathbb{C}}^{MB\times K}$ and another matrix,
$\mathbf{X}\in\operatorname{\mathbb{C}}^{MB\times L}$, as follows:

$\displaystyle\mathbf{D}$
$\displaystyle=\left[\mathbf{b}_{1}\otimes\mathbf{a}_{1},\cdots,\mathbf{b}_{K}\otimes\mathbf{a}_{K}\right],$
(77) $\displaystyle\mathbf{X}$
$\displaystyle=\big{[}\textnormal{vec}\left([\mathbf{I}_{B},\
\mathbf{0}_{B\times
R,1}]^{\mathsf{T}}-\alpha_{\textnormal{b}}\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\right),\cdots,\textnormal{vec}\left([\mathbf{I}_{B},\
\mathbf{0}_{B\times
R,L}]^{\mathsf{T}}-\alpha_{\textnormal{b}}\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\right)\big{]}.$
(78)

Through vectorization, we have the following relation for the norm of a
matrix:

$\left\lVert\mathbf{D}\bm{\upsilon}-\textnormal{vec}(\mathbf{C})\right\rVert_{\textnormal{2}}^{2}=\left\lVert\mathbf{A}\textnormal{Diag}(\bm{\upsilon})\mathbf{B}^{\mathsf{T}}-\mathbf{C}\right\rVert_{\textnormal{F}}^{2}.$
(79)

By using the relation in (79) and substituting the matrices $\mathbf{D}$ and
$\mathbf{X}$ into (76), we get:

$f_{1}=\left\lVert\mathbf{D}[\bm{\upsilon}_{1},~{}\bm{\upsilon}_{2},\cdots,\bm{\upsilon}_{L}]-\mathbf{X}\right\rVert_{\textnormal{F}}^{2},$
(80)

or equivalently:

$f_{1}=\left\lVert\mathbf{D}\bm{\Upsilon}-\mathbf{X}\right\rVert_{\textnormal{F}}^{2}.$
(81)

Similarly, by defining the matrices:

$\displaystyle\mathbf{M}$
$\displaystyle=\sqrt{P_{\textnormal{s}}}\alpha_{\textnormal{s}}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\mathbf{H}_{\textnormal{b-s}}\mathbf{v}_{\textnormal{b}}),$
(82) $\displaystyle\mathbf{Z}$ $\displaystyle=\left[\mathbf{0}_{L\times B},\
\mathbf{S}_{\textnormal{s}}^{\mathsf{T}}\right]^{\mathsf{T}},$ (83)

and then substituting them into (75) one obtains:

$f_{2}=\left\lVert\mathbf{M}\bm{\Upsilon}-\mathbf{Z}\right\rVert_{\textnormal{F}}^{2}.$
(84)

## References

  * [1] T. S. Rappaport, Y. Xing, O. Kanhere, S. Ju, A. Madanayake, S. Mandal, A. Alkhateeb, and G. C. Trichopoulos, “Wireless Communications and applications above 100 GHz: Opportunities and challenges for 6G and beyond,” _IEEE Access_ , vol. 7, pp. 78 729–78 757, 2019.
  * [2] Q. Wu and R. Zhang, “Towards smart and reconfigurable environment: Intelligent reflecting surface aided Wireless Network,” _IEEE Communications Magazine_ , vol. 58, no. 1, pp. 106–112, 2020.
  * [3] S. Hu, F. Rusek, and O. Edfors, “Beyond Massive MIMO: The Potential of data transmission with large intelligent surfaces,” _IEEE Transactions on Signal Processing_ , vol. 66, no. 10, pp. 2746–2758, 2018.
  * [4] Y. Liu, X. Liu, X. Mu, T. Hou, J. Xu, M. Di Renzo, and N. Al-Dhahir, “Reconfigurable intelligent surfaces: Principles and opportunities,” _IEEE Communications Surveys Tutorials_ , pp. 1–1, 2021.
  * [5] H. Q. Ngo, E. G. Larsson, and T. L. Marzetta, “Energy and spectral efficiency of very large multiuser MIMO systems,” _IEEE Transactions on Communications_ , vol. 61, no. 4, pp. 1436–1449, 2013.
  * [6] E. G. Larsson, O. Edfors, F. Tufvesson, and T. L. Marzetta, “Massive MIMO for next generation wireless systems,” _IEEE Communications Magazine_ , vol. 52, no. 2, pp. 186–195, 2014.
  * [7] B. Guo, C. Sun, and M. Tao, “Two-way passive beamforming design for RIS-aided FDD communication systems,” in _2021 IEEE Wireless Communications and Networking Conference (WCNC)_ , 2021, pp. 1–6.
  * [8] H. U. Rehman, F. Bellili, A. Mezghani, and E. Hossain, “Joint active and passive beamforming design for IRS-assisted multi-User MIMO systems: A VAMP-based approach,” _IEEE Transactions on Communications_ , pp. 1–1, 2021\.
  * [9] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 11, pp. 5394–5409, 2019.
  * [10] E. Basar, M. Di Renzo, J. De Rosny, M. Debbah, M.-S. Alouini, and R. Zhang, “Wireless communications through reconfigurable intelligent surfaces,” _IEEE Access_ , vol. 7, pp. 116 753–116 773, 2019.
  * [11] X. Gan, C. Zhong, C. Huang, and Z. Zhang, “RIS-assisted multi-user MISO communications exploiting statistical CSI,” _IEEE Transactions on Communications_ , pp. 1–1, 2021.
  * [12] C. Huang, A. Zappone, G. C. Alexandropoulos, M. Debbah, and C. Yuen, “Reconfigurable intelligent surfaces for energy efficiency in wireless communication,” _IEEE Transactions on Wireless Communications_ , vol. 18, no. 8, pp. 4157–4170, 2019.
  * [13] L. You, J. Xiong, D. W. K. Ng, C. Yuen, W. Wang, and X. Gao, “Energy efficiency and spectral efficiency tradeoff in RIS-aided multiuser MIMO uplink transmission,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 1407–1421, 2021.
  * [14] W. Yan, X. Yuan, Z.-Q. He, and X. Kuai, “Passive beamforming and information transfer design for reconfigurable intelligent surfaces aided multiuser mimo systems,” _IEEE Journal on Selected Areas in Communications_ , vol. 38, no. 8, pp. 1793–1808, 2020.
  * [15] S. Lin, B. Zheng, G. C. Alexandropoulos, M. Wen, M. di Renzo, and F. Chen, “Joint passive beamforming and information transfer for RIS-empowered wireless communications,” in _GLOBECOM 2020 - 2020 IEEE Global Communications Conference_ , 2020, pp. 1–6.
  * [16] H. Zhao, Y. Shuang, M. Wei, T. J. Cui, P. D. Hougne, and L. Li, “Metasurface-assisted massive backscatter wireless communication with commodity Wi-Fi signals,” _Nature Communications_ , vol. 11, no. 3926, pp. 1–10, 2020.
  * [17] C. Xu, L. Yang, and P. Zhang, “Practical backscatter communication systems for battery-free Internet of Things: A tutorial and survey of recent research,” _IEEE Signal Processing Magazine_ , vol. 35, no. 5, pp. 16–27, 2018.
  * [18] C. Yao, Y. Liu, X. Wei, G. Wang, and F. Gao, “Backscatter technologies and the future of internet of things: Challenges and opportunities,” _Intelligent and Converged Networks_ , vol. 1, no. 2, pp. 170–180, 2020.
  * [19] L. Zhang, X. Chen, and S. L. et al., “Space-time-coding digital metasurfaces,” _Nat Commun_ , vol. 9, no. 4334, 2018.
  * [20] D. P. Bertsekas, “Nonlinear programming,” _Journal of the Operational Research Society_ , vol. 48, no. 3, pp. 334–334, 1997.
  * [21] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network: Joint active and passive beamforming design,” in _2018 IEEE Global Communications Conference (GLOBECOM)_ , 2018, pp. 1–6.
  * [22] M. Joham, W. Utschick, and J. A. Nossek, “Linear transmit processing in MIMO communications systems,” _IEEE Transactions on Signal Processing_ , vol. 53, no. 8, pp. 2700–2712, 2005.
  * [23] H. Jedda, A. Mezghani, A. L. Swindlehurst, and J. A. Nossek, “Precoding under instantaneous per-antenna peak power constraint,” in _2017 25th European Signal Processing Conference (EUSIPCO)_ , 2017, pp. 863–867.
  * [24] R. W. Heath Jr. and A. Lozano, _Foundations of MIMO Communication_. Cambridge University Press, 2018.
  * [25] S. Rangan, P. Schniter, and A. K. Fletcher, “Vector approximate message passing,” _IEEE Transactions on Information Theory_ , vol. 65, no. 10, pp. 6664–6684, 2019.
  * [26] A. Manoel, F. Krzakala, G. Varoquaux, B. Thirion, and L. Zdeborová, “Approximate message-passing for convex optimization with non-separable penalties,” 2018.
  * [27] J. Yeh, _Real Analysis_ , 3rd ed. WORLD SCIENTIFIC, 2014. [Online]. Available: https://www.worldscientific.com/doi/abs/10.1142/9037
  * [28] Q.-U.-A. Nadeem, H. Alwazani, A. Kammoun, A. Chaaban, M. Debbah, and M.-S. Alouini, “Intelligent reflecting surface-assisted multi-user MISO communication: Channel estimation and beamforming design,” _IEEE Open Journal of the Communications Society_ , vol. 1, pp. 661–680, 2020.
  * [29] Q.-U.-A. Nadeem, A. Kammoun, A. Chaaban, M. Debbah, and M.-S. Alouini, “Intelligent reflecting surface assisted wireless communication: Modeling and channel estimation,” 2019.
  * [30] W. Gifford, M. Win, and M. Chiani, “Diversity with practical channel estimation,” _IEEE Transactions on Wireless Communications_ , vol. 4, no. 4, pp. 1935–1947, 2005.
  * [31] R. Annavajjala, P. C. Cosman, and L. B. Milstein, “Performance analysis of linear modulation schemes With generalized diversity combining on Rayleigh fading channels with noisy channel estimates,” _IEEE Transactions on Information Theory_ , vol. 53, no. 12, pp. 4701–4727, 2007.
  * [32] J. Zhang, M. Kountouris, J. G. Andrews, and R. W. Heath, “Multi-mode transmission for the MIMO broadcast channel with imperfect channel state information,” _IEEE Transactions on Communications_ , vol. 59, no. 3, pp. 803–814, 2011.