# Joint Active and Passive Beamforming Design for IRS-Assisted Multi-User MIMO
Systems: A VAMP-Based Approach

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).

###### Abstract

This paper tackles the problem of joint active and passive beamforming
optimization for an intelligent reflective surface (IRS)-assisted multi-user
downlink multiple-input multiple-output (MIMO) communication system under both
ideal and practical IRS phase shifts. We aim to maximize the spectral
efficiency of the users by minimizing the sum mean square error (MSE) of the
users’ received symbols. For this, a joint non-convex optimization problem is
formulated under the sum minimum mean square error (MMSE) criterion.
Alternating minimization is used to break the original joint optimization
problem into the separate optimization of the active precoding matrix for the
base station (BS) and the matrix of phase shifts for the IRS. While the MMSE
active precoder is obtained in closed-form, the IRS phase shifts are optimized
iteratively using a modified version (developed in this paper) of the vector
approximate message passing (VAMP) algorithm. Moreover, the underlying joint
optimization problem is solved under two different models for the IRS phase
shifts, namely by assuming $i)$ a unimodular (i.e., ideal) constraint on the
reflection coefficients and $ii)$ a more practical reflection elements
termination by a variable reactive load (which inherently introduces the
phase-dependent amplitude attenuation in the IRS phase shifts). Simulation
results are presented to illustrate the performance of the proposed method
under both perfect and imperfect channel state information (CSI) and to show
the effect of the practical constraint on the system throughput. The results
validate the superiority of the proposed method over the state-of-the-art
techniques both in terms of throughput and computational complexity.

###### Index Terms:

Intelligent reflective surface (IRS), IRS-assisted multi-user MIMO, joint
active and passive beamforming, vector approximate message passing (VAMP),
optimization

##  I Introduction

###  I-A Background

The need for higher data rates in wireless communication is soaring. This
calls for innovative and economically viable communication technologies that
can keep up with the increasing network capacity requirement. Massive
multiple-input multiple-output (MIMO) technology can fulfill the network
capacity requirement for beyond fifth-generation (B5G) wireless networks [1,
2, 3]. The basic idea of massive MIMO is to equip the base stations (BSs) with
tens (if not hundreds) of antenna elements so as to simultaneously serve
multiple mobile devices using the same time/frequency resources. Despite the
many advantages of massive MIMO, its practical large-scale deployment is
hindered by the associated high hardware cost and energy consumption [4, 5].
Moreover, although millimeter wave (mmWave) communication benefits from
massive MIMO due to a symbiotic convergence of technologies, its practical use
is still limited by the less penetrative propagation characteristic of mmWave
signals in presence of blockages between the BS and the mobile device. [6].

One promising technology that has been introduced recently is intelligent
reflective surfaces (IRSs), also called reconfigurable intelligent surfaces
(RISs) [7, 8]. IRS is composed of a planar metasurface consisting of a large
number of passive reflective elements. Moreover, IRS does not require a power
amplifier for transmission which makes it an energy-efficient technology. This
allows the IRS to passively alter the signal propagation by reconfiguring the
phases of its reflective elements through a controller attached to the surface
[9]. Therefore, IRSs can be utilized to perform passive beamforming. Passive
beamforming refers to changing the IRS phases without actively powering the
IRS antenna elements as opposed to active beamforming at the BS so as to
improve the received power while reducing the interference for unintended
users, thereby enhancing the overall throughput of the network [10].
Practically, IRS deployment requires a large number of cost-effective phase
shifters (PSs) on a surface that can be easily integrated into a traditional
wireless network [11]. Due to the aforementioned reasons, IRS-assisted
communication has gained substantial research interest in the wireless
research community over the recent few years [10, 11, 12, 13, 14, 15, 16, 17,
18].

In [14], a multi-user multiple-input single-output (MISO) wireless system
assisted by a single IRS in the downlink configuration is studied. The authors
present a deep reinforcement learning (DRL)-based solution to jointly optimize
the IRS phase shifts and the BS precoding under different quality of service
constraints. In [18], authors tackle the problem of estimating the cascaded
BS-IRS-user channels for an IRS-assisted multi-user MISO system. The author
propose a pilot-based solution and improve its efficiency by exploiting the
fact that all users share the same BS-IRS channel. The same problem is solved
by utilizing the deep residual learning framework in [19]. In [20], an IRS-
assisted multi-cluster MISO system serving multiple users is considered
wherein the authors seek to minimize the transmit power under a minimum
signal-to-interference-plus-noise ratio (SINR) constraint by jointly
optimizing the IRS phase shifts and the transmit precoder. They tackle the
underlying problem through alternating direction method of multipliers (ADMM).
An IRS-aided MISO and MIMO system with discrete phase shifts for IRS elements
is also discussed in [12]. The authors formulate the problem of minimizing the
transmit power under minimum SINR constraint and jointly optimize the transmit
precoding and the IRS phase shifts in a mixed-integer non-linear programming
framework. In [17], a relatively more practical model for IRS reflection
coefficients is considered, and then a penalty-based algorithm is used to
optimize the phase matrix.

The vast majority of the existing work considers a MISO wireless system
assisted by a single or multiple IRSs serving a single user [13, 16, 17]. So
far, limited research has been conducted on IRS-aided multi-user MIMO systems.
Moreover, IRS reflection coefficients are often modeled as ideal phase
shifters and a realistic approach towards modeling reflection coefficients has
rarely been investigated. In fact, most of the existing methods are limited to
a single-phase shifter model, unimodular phase shifts being the most common
one, and hence they are not resilient to the various hardware impairments of
the IRS reflection elements [13, 12, 16, 10, 17].

###  I-B Contributions

In this paper, we consider a multi-user IRS-assisted single-cell downlink MIMO
system with a single IRS. The IRS is equipped with a large number of passive
phase shifters that aid the BS to serve a small number of users. We propose a
robust solution for the problem of jointly optimizing the active and passive
beamforming tasks under different models for the IRS reflection coefficients.
The main contributions embodied by this paper are as follows:

  * •

We solve the problem of maximizing the spectral efficiency of the users by
jointly optimizing the transmit precoding matrix at the BS and the reflection
coefficients at the IRS. To that end, we first formulate the joint
optimization problem under the sum MMSE criterion in order to minimize the MSE
of the received symbols for all users at the same time.

  * •

To solve the underlying joint optimization problem, we first split it using
alternate optimization [21] into two easier sub-optimization tasks of the
active precoder at the BS and the reflection coefficients at the IRS. The
precoding sub-optimization problem is similar to the MMSE transmit precoder
optimization for a traditional MIMO system, which can be solved in closed-form
through Lagrange optimization.

  * •

We modify and extend the existing VAMP algorithm [22] and propose a robust
technique to find locally optimal reflective coefficients for the IRS under
multiple constraints. Precisely, we find the sub-optimal phase matrix under
two different models for the reflection coefficients: $i)$ Under the
unimodular constraint on the IRS reflection coefficients and $ii)$ under a
practical constraint, where each IRS element is terminated by a tunable simple
reactive load.

  * •

We discuss the convergence and provide the order of complexity of the proposed
solution. We present various numerical results to compare the proposed
solution with the semi-definite relaxation (SDR) plus MMSE-based IRS
beamforming and precoding optimization approach [23, 10], an ADMM-based
solution, and a standalone massive MIMO system using MMSE precoder. The
results show that, the proposed solution: $i)$ outperforms both the SDR-based
and the ADMM-based solutions in terms of throughput while using the same
resources and being less computationally demanding, and $ii)$ achieves higher
throughput than a traditional massive MIMO system while using a significantly
smaller number of transmit antennas in typical propagation scenarios. We
illustrate the effect of practical phase shifts on the system throughout. We
also show the robustness of the proposed solution by assessing its performance
under imperfect CSI.

###  I-C Paper Organization and Notations

The rest of the paper is organized as follows: the system model along with the
problem formulation for jointly optimizing the active precoder and the
reflection coefficients are discussed in Section II. Section III briefly
introduces the VAMP algorithm and then extends it to solve optimization
problems. In Section IV, we solve the optimization problem at hand using the
proposed extended version of VAMP. In Section V, we further solve the
underlying optimization problem under the “simple reactive loading” constraint
on the IRS reflection elements. Exhaustive numerical results are shown in
Section VI. Finally, Section VII provides an analysis on the convergence and
computational complexity of the proposed solution.

Notations: Lowercase letters (e.g., $r$) denote scalar variables. The
uppercase italic letters (e.g., $N$) represent scalar constants. 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 exponent 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}$.
$\operatorname{\mathbb{C}}^{M\times N}$ denotes 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
$\mathbf{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. The imaginary unit is represented by ${\rm j}=\sqrt{-1}$ and the
$\angle(.)$ operator returns the angle of its complex argument. The
proportional relationship between any two entities (functions or variables) is
represented by $\varpropto$ operator. Lastly, the operators $\otimes$, $\odot$
and $\ast$ denote 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 serving $M\ (M<N)$
single-antenna users in the downlink. The BS is assisted by an IRS which has
$K\ (K>M)$ reflective elements. For each $m$-th user, we have a direct link to
the BS expressed by a channel vector
$\mathbf{h}_{\textnormal{b-u},m}\in\operatorname{\mathbb{C}}^{N}$. The channel
of the surface-user $m$ link is denoted by
$\mathbf{h}_{\textnormal{s-u},m}\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 IRS-BS link with
$\textnormal{Rank}(\mathbf{H}_{\textnormal{b-s}})\geq M$. The signal received
at the IRS is phase-shifted by a diagonal matrix
$\mathbf{\Upsilon}=\textnormal{Diag}(\bm{\upsilon})\in\operatorname{\mathbb{C}}^{K\times
K}$, where $\bm{\upsilon}\in\operatorname{\mathbb{C}}^{K}$ is the phase-shift
vector having unimodular elements, i.e.,
$\left\lvert\upsilon_{k}\right\rvert=1$ for $k=1,\cdots,K$. In other words,
for each reflection element, we have $\upsilon_{k}=e^{{\rm j}\theta_{k}}$ for
some phase shift $\theta_{k}\in[0,2\pi]$. The received signal for user $m$ can
be expressed as follows:

$\mathsf{y}_{m}=\alpha\left(\mathbf{h}_{\textnormal{s-u},m}^{\mathsf{H}}{\mathbf{\Upsilon}}\mathbf{H}_{\textnormal{b-s}}\displaystyle\sum_{m^{\prime}=1}^{M}\mathbf{f}_{m^{\prime}}\mathsf{s}_{m^{\prime}}+\mathbf{h}_{\textnormal{b-u},m}^{\mathsf{H}}\displaystyle\sum_{m^{\prime}=1}^{M}\mathbf{f}_{m^{\prime}}\mathsf{s}_{m^{\prime}}+\
\mathsf{w}\right),~{}~{}m=1,\cdots,M,$ (1)

where
$\mathsf{s}_{m}\sim\operatorname{\mathcal{C}}\operatorname{\mathcal{N}}(s;0,1)$
is the unknown transmit symbol,
$\mathsf{w}\sim\operatorname{\mathcal{C}}\operatorname{\mathcal{N}}(w;0,\sigma_{\textnormal{w}}^{2})$
denotes additive white Gaussian noise (AWGN), and
$\alpha\in\operatorname{\mathbb{R}}$ refers to the receiver scaling which is a
common practice in precoding optimization literature [24, 25]. Here,
$\mathbf{f}_{m}\in\operatorname{\mathbb{C}}^{N}$ for $m=1,\cdots,M$ are the
precoding vectors that are used for power allocation and beamforming purposes.
Let $\mathbf{F}=[\mathbf{f}_{1},\mathbf{f}_{2},\cdots,\mathbf{f}_{M}]$ be the
precoding matrix and let $P$ denote the total transmit power. By denoting
$\bm{\mathsf{s}}=[\mathsf{s}_{1},\mathsf{s}_{2},\cdots,\mathsf{s}_{M}]^{\mathsf{T}}$,
it follows that
$\operatorname{\mathbb{E}}\big{\\{}\|\mathbf{F}\bm{\mathsf{s}}\|^{2}\big{\\}}=P$.
Let
$\mathbf{H}_{\textnormal{b-u}}=\left[\mathbf{h}_{\textnormal{b-u},1},\mathbf{h}_{\textnormal{b-u},2},\cdots,\mathbf{h}_{\textnormal{b-u},M}\right]$
and
$\mathbf{H}_{\textnormal{s-u}}=[\mathbf{h}_{\textnormal{s-u},1},\mathbf{h}_{\textnormal{s-u},2},\cdots,\mathbf{h}_{\textnormal{s-u},M}]$.
Then, by stacking all the users’ signals in one vector
$\bm{\mathsf{y}}=[\mathsf{y}_{1},\mathsf{y}_{2},\cdots,\mathsf{y}_{M}]^{\mathsf{T}}$,
we can express the input-output relationship of the multi-user MIMO system as:

$\bm{\mathsf{y}}=\alpha\left(\underbrace{\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}{\mathbf{\Upsilon}}\mathbf{H}_{\textnormal{b-s}}\mathbf{F}\bm{\mathsf{s}}}_{\textnormal{Users-
IRS-
BS}}+\underbrace{\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F}\bm{\mathsf{s}}}_{\textnormal{Users-
BS}}+\ \bm{\mathsf{w}}\right).$ (2)

BSIRS $K$ reflectorsUser $1$User
$M$$\mathbf{h}_{\textnormal{b-u},1}$$\mathbf{h}_{\textnormal{b-u},M}$$\mathbf{h}_{\textnormal{s-u},1}$$\mathbf{h}_{\textnormal{s-u},M}$$\mathbf{H}_{\textnormal{b-s}}$$\vdots$$\vdots$$\cdots$$\cdots$$\ddots$

Figure 1: IRS-assisted multi-user MIMO system.

The overall effective channel matrix for all users is thus given by:

$\mathbf{H}^{\mathsf{H}}=\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}{\mathbf{\Upsilon}}\mathbf{H}_{\textnormal{b-s}}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}.$
(3)

We aim to minimize the received symbol error of each user under the MMSE
criterion, which consequently maximizes the user SINR. A lower bound on the
spectral efficiency for user $m$ can be expressed in terms of the MMSE of its
received symbol [26] as follows:

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

The MSE of the received symbol for user $m$ is given by
$\operatorname{\mathbb{E}}_{\mathsf{y}_{m},\mathsf{s}_{m}}\left\\{\left\lvert\mathsf{y}_{m}-\mathsf{s}_{m}\right\rvert^{2}\right\\}$,
and for $M$ users, the sum symbol MSE can be written as:

$\displaystyle\sum_{m=1}^{M}\operatorname{\mathbb{E}}_{\mathsf{y}_{m},\mathsf{s}_{m}}\left\\{\left\lvert\mathsf{y}_{m}-\mathsf{s}_{m}\right\rvert^{2}\right\\}\
=\
\operatorname{\mathbb{E}}_{\bm{\mathsf{y}},\bm{\mathsf{s}}}\left\\{\left\lVert\bm{\mathsf{y}}-\bm{\mathsf{s}}\right\rVert_{2}^{2}\right\\}.$
(5)

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

$\displaystyle\displaystyle\operatorname*{arg\,min}_{\alpha,\mathbf{F},{\mathbf{\Upsilon}}}\quad$
$\displaystyle\operatorname{\mathbb{E}}_{\bm{\mathsf{y}},\bm{\mathsf{s}}}\left\\{\|\bm{\mathsf{y}}-\bm{\mathsf{s}}\|_{2}^{2}\right\\},$
(6a) subject to
$\displaystyle\operatorname{\mathbb{E}}_{\bm{\mathsf{s}}}\left\\{\|\mathbf{F}\bm{\mathsf{s}}\|_{2}^{2}\right\\}=P,$
(6b) $\displaystyle\upsilon_{ik}=0,\quad i\neq k,$ (6c)
$\displaystyle|\upsilon_{ii}|=1,\quad i=1,2,\cdots,K.$ (6d)

######  Remark.

The objective function in (6a) leads to some fairness among the users by
ensuring that the MSE is minimized for each user. The lower bound on sum-
spectral-efficiency of $M$ users can be expressed in terms of the MMSE of the
users’ received symbols [26] as follows:

$\widehat{C}\ =\
\displaystyle\sum_{m=1}^{M}\log_{2}\left(\frac{1}{\textnormal{MMSE}_{m}}\right)\
=\
\log_{2}\left(\displaystyle\prod_{m=1}^{M}\frac{1}{\textnormal{MMSE}_{m}}\right).$
(7)

In other words, maximizing the sum-spectral-efficiency is equivalent to
minimizing the product MSE of all users. This can be achieved by minimizing
the MSE of the user with the strongest channel, thereby leading to a very
unfair solution. On the other hand, aiming for complete fairness results in a
very inefficient allocation of resources when it comes to the overall system
throughput. In this respect, the sum MMSE criterion is a good balance between
the two extremes. Since our aim is to maximize the spectral efficiency of each
user rather than the sum-spectral-efficiency, the MMSE criterion is a good fit
for our problem formulation.

The expectation involved in (6a) and (6b) is taken with respect to (w.r.t.)
the random vectors $\bm{\mathsf{s}}$ and $\bm{\mathsf{w}}$. Explicitly writing
the objective function in (6a) leads to:

$\operatorname{\mathbb{E}}_{\bm{\mathsf{w}},\bm{\mathsf{s}}}\big{\\{}\textnormal{Tr}\left(\alpha^{2}\bm{\mathsf{s}}^{\mathsf{H}}\mathbf{F}^{\mathsf{H}}\mathbf{H}\mathbf{H}^{\mathsf{H}}\mathbf{F}\bm{\mathsf{s}}-\alpha\bm{\mathsf{s}}^{\mathsf{H}}\mathbf{F}^{\mathsf{H}}\mathbf{H}\bm{\mathsf{s}}-\alpha\bm{\mathsf{s}}^{\mathsf{H}}\mathbf{H}^{\mathsf{H}}\mathbf{F}\bm{\mathsf{s}}+\bm{\mathsf{s}}\bm{\mathsf{s}}^{\mathsf{H}}\right)\\\
+\textnormal{Tr}\left(\alpha^{2}\bm{\mathsf{s}}^{\mathsf{H}}\mathbf{F}^{\mathsf{H}}\mathbf{H}\bm{\mathsf{w}}+\alpha^{2}\bm{\mathsf{w}}^{\mathsf{H}}\mathbf{H}^{\mathsf{H}}\mathbf{F}\bm{\mathsf{s}}-\alpha\bm{\mathsf{w}}^{\mathsf{H}}\bm{\mathsf{s}}-\alpha\bm{\mathsf{s}}^{\mathsf{H}}\bm{\mathsf{w}}+\alpha^{2}\bm{\mathsf{w}}^{\mathsf{H}}\bm{\mathsf{w}}\right)\big{\\}},$
(8)

thereby resulting in the following optimization problem.

$\displaystyle\displaystyle\operatorname*{arg\,min}_{\alpha,\mathbf{F},{\mathbf{\Upsilon}}}\quad$
$\displaystyle\left\lVert\alpha\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}{\bm{\Upsilon}}\mathbf{H}_{\textnormal{b-s}}\mathbf{F}-(\mathbf{I}_{M}-\alpha\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F})\right\rVert_{\textnormal{F}}^{2}+M\alpha^{2}\sigma_{\textnormal{w}}^{2},$
(9a) s.t. $\displaystyle\|\mathbf{F}\|_{\textnormal{F}}^{2}=P,$ (9b)
$\displaystyle\upsilon_{ik}=0,\quad i\neq k,$ (9c)
$\displaystyle|\upsilon_{ii}|=1,\quad i=1,2,\cdots,K.$ (9d)

The optimization problem in (9) is non-convex optimization problem due to the
unimodular constraint111Later, we will solve the same problem under another
constraint on the reflection coefficients. on the IRS phase shifts in (9d).
VAMP is a low-complexity algorithm which is designed to solve optimization
problems with a linear objective function and non-linear constraints [22].
VAMP has a modular structure that makes it possible to decouple the
constraints from the objective function. Therefore, the same objective
function can be minimized under different constraints by modifying the VAMP
module that satisfies the constraint (simple scalar functions). VAMP
automatically updates the stepsize at a per-iteration basis that leads to a
faster convergence compared to other iterative algorithms (e.g., ADMM) [27].
This favorable property makes VAMP tuning-free. The performance of VAMP can be
theoretically predicted to establish optimality through the statistical state
evolution framework [28, 22]. This is, however left to future work.

##  III Modified VAMP Algorithm for Constrained Optimization

Recently, message passing algorithms [22, 29, 30] have gained attention in
estimation theory because of their high performance and fast convergence.
Vector approximate message passing (VAMP) [22], in particular, is a low-
complexity algorithm that solves quadratic loss optimization of recovering a
vector from noisy linear measurements. In this section, we briefly discuss the
standard max-sum VAMP algorithm and we further modify it to solve the
constrained optimization problem at hand.

###  III-A Background on Max-Sum VAMP

Approximate message passing (AMP)-based computational techniques have gained a
lot of attention since their introduction within the compressed sensing
framework [29]. To be precise, AMP solves the standard linear regression
problem of recovering a vector $\mathbf{x}\in\operatorname{\mathbb{C}}^{N}$
from noisy linear observations:

$\bm{\mathsf{z}}=\mathbf{A}\mathbf{x}+\bm{\mathsf{w}},$ (10)

where $\mathbf{A}\in\operatorname{\mathbb{C}}^{M\times N}$ (with $M\ll N$) is
called sensing matrix and
$\bm{\mathsf{w}}\sim\operatorname{\mathcal{C}}\operatorname{\mathcal{N}}(\bm{\mathsf{w}};0,\gamma_{\textnormal{w}}^{-1}\mathbf{I}_{M})$,
with $\gamma_{\textnormal{w}}>0$, so that
$p_{\bm{\mathsf{z}}|\bm{\mathsf{x}}}(\mathbf{z}|\mathbf{x})=\operatorname{\mathcal{C}}\operatorname{\mathcal{N}}(\mathbf{z};\mathbf{A}\mathbf{x},\gamma_{\textnormal{w}}^{-1}\mathbf{I}_{M})$.
Interestingly, the performance of AMP under independent and identically
distributed (i.i.d.) Gaussian sensing matrices, $\mathbf{A}$, can be
rigorously tracked through scalar state evolution (SE) equations [31]. One
major drawback of AMP, however, is that it often diverges if the sensing
matrix, $\mathbf{A}$, is ill-conditioned or has a non-zero mean. To circumvent
this problem, vector AMP (VAMP) algorithm was proposed and rigorously analyzed
through SE equations in [22]. Although there is no theoretical guarantee that
VAMP will always converge, strong empirical evidence suggests that VAMP is
more resilient to badly conditioned sensing matrices given that they are
right-orthogonally invariant [22]. Consider the joint probability distribution
function (pdf) of $\bm{\mathsf{x}}$ and $\bm{\mathsf{z}}$,
$p_{\bm{\mathsf{x}},\bm{\mathsf{z}}}(\mathbf{x},\mathbf{z})$

$p_{\bm{\mathsf{x}},\bm{\mathsf{z}}}(\mathbf{x},\mathbf{z})=p_{\bm{\mathsf{x}}}(\mathbf{x})\operatorname{\mathcal{C}}\operatorname{\mathcal{N}}(\mathbf{z};\mathbf{A}\mathbf{x},\gamma_{\textnormal{w}}^{-1}\mathbf{I}_{M}).$
(11)

Here $p_{\bm{\mathsf{x}}}(\mathbf{x})$ is some prior distribution on the
vector $\mathbf{x}$ whose elements are assumed to be i.i.d. with a common
prior distribution, $p_{\mathsf{x}}(x)$, i.e.:

$p_{\bm{\mathsf{x}}}(\mathbf{x})=\displaystyle\prod_{i=1}^{N}p_{\mathsf{x}}(x_{i}).$
(12)

Max-sum VAMP can solve the following optimization problem:

$\widehat{\mathbf{x}}=\displaystyle\operatorname*{arg\,min}_{\mathbf{x}}\left\lVert\mathbf{z}-\mathbf{A}\mathbf{x}\right\rVert^{2},$
(13)

by finding the maximum a posteriori (MAP) estimate of $\mathbf{x}$ as follows:

$\widehat{\mathbf{x}}=\operatorname*{arg\,max}_{\mathbf{x}}p_{\bm{\mathsf{x}}|\bm{\mathsf{z}}}(\mathbf{x}|\mathbf{z}).$
(14)

The algorithm consists of the following two modules.

####  III-A1 Linear MAP/MMSE Estimator

At iteration $t$, the linear MAP estimator receives extrinsic information
(message) from the separable (i.e., element-wise) MAP denoiser of $\mathbf{x}$
in the form of a mean vector, $\mathbf{r}_{t-1}$, and a common scalar
precision, $\gamma_{t-1}$. Then, under the Gaussian prior,
$\operatorname{\mathcal{C}}\operatorname{\mathcal{N}}(\mathbf{x};\mathbf{r}_{t-1},\gamma_{t-1}^{-1}\mathbf{I}_{N})$,
it computes the linear MAP estimate, $\bar{\mathbf{x}}_{t}$, along with the
associated posterior precision, $\bar{\gamma}_{t}$, from the linear
observations, $\mathbf{z}=\mathbf{A}\mathbf{x}+\bm{\mathsf{w}}$ on
$\mathbf{x}$. Because we are dealing with Gaussian densities, the linear MAP
estimate is equal to the linear MMSE (LMMSE) and given as follows:

$\displaystyle\bar{\mathbf{x}}_{t}$
$\displaystyle~{}=~{}\left(\gamma_{\textnormal{w}}\mathbf{A}^{\mathsf{H}}\mathbf{A}+\gamma_{t-1}\mathbf{I}_{N}\right)^{-1}\left(\gamma_{\textnormal{w}}\mathbf{A}^{\mathsf{H}}\mathbf{z}+\gamma_{t-1}\mathbf{r}_{t-1}\right),$
(15) $\displaystyle\bar{\gamma}_{t}$
$\displaystyle~{}=~{}N\textnormal{Tr}\left(\left[\gamma_{\textnormal{w}}\mathbf{A}^{\mathsf{H}}\mathbf{A}+\gamma_{t-1}\mathbf{I}_{N}\right]^{-1}\right)^{-1}.$
(16)

The extrinsic information on $\mathbf{x}$ is updated as
$\operatorname{\mathcal{C}}\operatorname{\mathcal{N}}(\mathbf{x};\bar{\mathbf{x}}_{t},\bar{\gamma}_{t}^{-1}\mathbf{I}_{N})/\operatorname{\mathcal{C}}\operatorname{\mathcal{N}}(\mathbf{x};\mathbf{r}_{t-1},\gamma_{t-1}^{-1}\mathbf{I}_{N}),$
and then sent back in the form of a mean vector,
$\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
denoiser of $\mathbf{x}$. The SVD (singular value decomposition) form of VAMP
directly computes extrinsic mean vector $\widetilde{\mathbf{r}}_{t}$ and
scalar precision $\widetilde{\gamma}_{t}$, and can be readily obtained by
substituting
$\mathbf{A}=\mathbf{U}\textnormal{Diag}(\bm{\omega})\mathbf{V}^{\mathsf{H}}$
in (15) and (16).

####  III-A2 Separable MAP Denoiser of $\mathbf{x}$

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{N}}(\mathbf{x};\widetilde{\mathbf{r}}_{t},\widetilde{\gamma}^{-1}_{t}\mathbf{I}_{N})$.
Because $\mathbf{x}$ is i.i.d., the MAP estimate can be computed through a
component-wise denoising function as follows:

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

or equivalently,

$g_{1,i}(\widetilde{r}_{i,t},\widetilde{\gamma}_{t})=\displaystyle\operatorname*{arg\,min}_{x_{i}}\left[\widetilde{\gamma}\left\lvert
x_{i}-\widetilde{r}_{i,t}\right\rvert^{2}-\ln p_{\mathsf{x}}(x_{i})\right].$
(18)

The derivative of the scalar MAP denoiser w.r.t. $\widetilde{r}_{i,t}$ is
given by [22]:

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

where $\widehat{\gamma}_{t}$ is the posterior precision. The vector valued
denoiser function and its derivaive are defined as follows:

$\displaystyle\mathbf{g}_{1}(\widetilde{\mathbf{r}}_{t},\widetilde{\gamma}_{t})$
$\displaystyle~{}\triangleq~{}\big{[}g_{1,1}(\widetilde{r}_{1,t},\widetilde{\gamma}_{t}),\
g_{1,2}(\widetilde{r}_{2,t},\widetilde{\gamma}_{t}),\cdots,g_{1,N}(\widetilde{r}_{N,t},\widetilde{\gamma}_{t})\big{]}^{\mathsf{T}},$
(20)
$\displaystyle\mathbf{g}^{\prime}_{1}(\widetilde{\mathbf{r}}_{t},\widetilde{\gamma}_{t})$
$\displaystyle~{}\triangleq~{}\big{[}g^{\prime}_{1,1}(\widetilde{r}_{1,t},\widetilde{\gamma}_{t}),\
g^{\prime}_{1,2}(\widetilde{r}_{2,t},\widetilde{\gamma}_{t}),\cdots,g^{\prime}_{1,N}(\widetilde{r}_{N,t},\widetilde{\gamma}_{t})\big{]}^{\mathsf{T}}.$
(21)

Similar to the LMMSE module, the MAP denoiser module computes an extrinsic
mean vector,
$\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 LMMSE module for the next iteration. The process is repeated until
convergence.  It is worth mentioning that the extrinsic parameters, i.e., the
extrinsic mean vector and the scalar precision, calculated by each module act
as a Gaussian prior on the succeeding estimate of the adjacent module, thus
making VAMP parameter-free. Another key advantage of VAMP is that it decouples
the prior information, $p_{\bm{\mathsf{x}}}(\mathbf{x})$, and the
observations, $p_{\bm{\mathsf{x}}|\bm{\mathsf{z}}}(\mathbf{z}|\mathbf{x})$,
into two separate modules. Moreover, it also enables the denoising function to
be separable even if the elements of $\mathbf{x}$ are correlated in which case
the LMMSE module can easily incorporate such correlation information. The
steps of the standard max-sum VAMP algorithm are shown in Algorithm 1.

Algorithm 1 Max-sum VAMP SVD

Given $\mathbf{A}\in\operatorname{\mathbb{C}}^{M\times N}$,
$\mathbf{z}\in\operatorname{\mathbb{C}}^{M}$, 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: Compute economy-size SVD
$\mathbf{A}=\mathbf{U}\textnormal{Diag}(\bm{\omega})\mathbf{V}^{\mathsf{H}}$

3:
$R_{\textnormal{A}}=\textnormal{Rank}(\mathbf{A})=\textnormal{length}(\bm{\omega})$

4: Compute
$\widetilde{\mathbf{z}}=\textnormal{Diag}(\bm{\omega})^{-1}\mathbf{U}^{\mathsf{H}}\mathbf{z}$

5: repeat

6: // LMMSE SVD Form.

7:
$\mathbf{d}_{t}=\gamma_{\textnormal{w}}\textnormal{Diag}(\gamma_{\textnormal{w}}\bm{\omega}^{2}+\gamma_{t-1}\mathbf{1}_{R_{\textnormal{A}}})^{-1}\bm{\omega}^{2}$

8:
$\widetilde{\mathbf{r}}_{t}=\mathbf{r}_{t-1}+\frac{N}{R_{\textnormal{A}}}\mathbf{V}\textnormal{Diag}\left(\mathbf{d}_{t}/\langle\mathbf{d}_{t}\rangle\right)\left(\widetilde{\mathbf{z}}-\mathbf{V}^{\mathsf{H}}\mathbf{r}_{t-1}\right)$

9:
$\widetilde{\gamma}_{t}=\gamma_{t-1}\left\langle\mathbf{d}_{t}\right\rangle/\left(\frac{N}{R_{\textnormal{A}}}-\left\langle\mathbf{d}_{t}\right\rangle\right)$

10: // MAP Denoiser

11:
$\widehat{\mathbf{x}}_{t}=\mathbf{g}_{1}(\widetilde{\mathbf{r}}_{t},\widetilde{\gamma}_{t})$

12:
$\widehat{\gamma}_{t}=\left\langle\mathbf{g}_{1}^{\prime}(\widetilde{\mathbf{r}}_{t},\widetilde{\gamma}_{t})\right\rangle/\widetilde{\gamma}_{t}$

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

14:
$\mathbf{r}_{t}=(\widehat{\gamma}_{t}\widehat{\mathbf{x}}_{t}-\widetilde{\gamma}_{t}\widetilde{\mathbf{r}}_{t})/\gamma_{t}$

15: $t\leftarrow t+1$

16: 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}}$

17: return $\widehat{\mathbf{x}}_{t}$

###  III-B VAMP for Optimization

In this section, we explain how max-sum VAMP can be applied to constrained
optimization problems. Given the knowledge of three matrices
$\mathbf{A}\in\operatorname{\mathbb{C}}^{M\times N}$,
$\mathbf{B}\in\operatorname{\mathbb{C}}^{Q\times N}$ and
$\mathbf{Z}\in\operatorname{\mathbb{C}}^{M\times Q}$, the goal is to solve an
optimization problem of the form:

$\displaystyle\displaystyle\operatorname*{arg\,min}_{\mathbf{x}~{}\in~{}\operatorname{\mathbb{C}}^{N}}\quad$
$\displaystyle\left\lVert\mathbf{A}\textnormal{Diag}(\mathbf{x})\mathbf{B}^{\mathsf{T}}-\mathbf{Z}\right\rVert_{\textnormal{F}}^{2}$
(22a) s.t. $\displaystyle f_{i}(x_{i})=0\quad i=1,\cdots,N.$ (22b)

In the context of optimization, the observation matrix, $\mathbf{Z}$, is
considered as the desired output matrix and it is also assumed to be known.
Unlike the estimation problem in (13), we do not have a prior distribution on
$\mathbf{x}$. Yet, the optimization problem in (22) can be solved by modifying
the modules of standard max-sum VAMP.

####  III-B1 Extended LMMSE

Through vectorization, the objective function in (22a) can be written in the
same form as the quadratic objective function in (13) in the following way:

$\textnormal{vec}(\mathbf{Z})=(\mathbf{B}\otimes\mathbf{A})\textnormal{vec}(\textnormal{Diag}(\mathbf{x})).$
(23)

We then define a matrix, $\mathbf{D}\in\operatorname{\mathbb{C}}^{MQ\times
N}$, as follows:

$\mathbf{D}\triangleq\mathbf{B}\ast\mathbf{A}=\left[\mathbf{b}_{1}\otimes\mathbf{a}_{1},\cdots,\mathbf{b}_{K}\otimes\mathbf{a}_{K}\right].$
(24)

Then, the objective function in (22a) is equivalently expressed in a standard
form that is amenable to VAMP as follows:

$\displaystyle\displaystyle\operatorname*{arg\,min}_{\mathbf{x}~{}\in~{}\operatorname{\mathbb{C}}^{N}}\quad$
$\displaystyle\left\lVert\mathbf{D}\mathbf{x}-\textnormal{vec}(\mathbf{Z})\right\rVert_{\textnormal{2}}^{2}$
(25a) s.t. $\displaystyle f_{i}(x_{i})=0\quad i=1,\cdots,N.$ (25b)

The column-wise Khatri-Rao structure can be exploited to avoid taking SVD of
the large matrix $\mathbf{D}$ as explained in the sequel. Let
$\mathbf{A}=\mathbf{U}_{\textnormal{A}}\textnormal{Diag}(\bm{\omega}_{\textnormal{A}})\mathbf{V}_{\textnormal{A}}^{\mathsf{H}}$,
$\mathbf{B}=\mathbf{U}_{\textnormal{B}}\textnormal{Diag}(\bm{\omega}_{\textnormal{B}})\mathbf{V}_{\textnormal{B}}^{\mathsf{H}}$,
$\mathbf{D}=\mathbf{U}\textnormal{Diag}(\bm{\omega})\mathbf{V}^{\mathsf{H}}$
and
$\mathbf{V}_{\textnormal{BA}}=\left(\mathbf{V}_{\textnormal{B}}^{\mathsf{H}}\ast\mathbf{V}_{\textnormal{A}}^{\mathsf{H}}\right)^{\mathsf{H}}$.
By defining the normalization vector:

$\mathbf{v}_{\textnormal{n}}=\big{[}\left\lVert\mathbf{v}_{\textnormal{BA},1}\right\rVert_{2},\
\left\lVert\mathbf{v}_{\textnormal{BA},2}\right\rVert_{2},\cdots,\left\lVert\mathbf{v}_{\textnormal{BA},MQ}\right\rVert_{2}\big{]}^{\mathsf{T}},$
(26)

it can be shown that the SVD of the matrix $\mathbf{D}$ is given by:

$\mathbf{D}=\underbrace{\left(\mathbf{U}_{\textnormal{B}}\otimes\mathbf{U}_{\textnormal{A}}\right)}_{\mathbf{U}}\underbrace{\textnormal{Diag}\big{(}(\bm{\omega}_{\textnormal{B}}\otimes\bm{\omega}_{\textnormal{A}})\odot\mathbf{v}_{\textnormal{n}}\big{)}}_{\textnormal{Diag}(\bm{\omega})}\underbrace{\left(\mathbf{V}_{\textnormal{B}}^{\mathsf{H}}\ast\mathbf{V}_{\textnormal{A}}^{\mathsf{H}}\right)\odot\left(\mathbf{v}_{\textnormal{n}}^{-1}\mathbf{1}_{N}^{\mathsf{T}}\right)}_{\mathbf{V}^{\mathsf{H}}}.$
(27)

These steps can be easily incorporated in the Algorithm 1 accordingly. Similar
to the standard max-sum VAMP, at iteration $t$, the LMMSE module receives an
extrinsic mean vector, $\mathbf{r}_{t-1}$, and a scalar precision,
$\gamma_{t-1}$, from the separable MAP estimator. The SVD form of VAMP allows
for exploiting the Kronecker structure inside the algorithm to avoid any large
matrix multiplication. The product of a Kronecker-structured matrix and a
vector can be computed in an efficient way through reverse vectorization or
unvectorization by computing the product of three smaller matrices, and then
vectorizing the result. Therefore, line $4$ of Algorithm 1 can be modified as
follows:

$\widetilde{\mathbf{z}}~{}=~{}\textnormal{Diag}(\bm{\omega})^{-1}\mathbf{U}^{\mathsf{H}}\textnormal{vec}(\mathbf{Z})~{}=~{}\textnormal{Diag}(\bm{\omega})^{-1}\textnormal{vec}\left(\mathbf{U}_{\textnormal{A}}^{\mathsf{H}}\mathbf{Z}\mathbf{U}_{\textnormal{B}}^{*}\right).$
(28)

SeparableMAPProjectorLMMSEEstimatorextext$\widehat{\mathbf{x}}$$\mathbf{Z}$$\widetilde{\mathbf{r}},\widetilde{\gamma}$$\bar{\mathbf{x}}$$\bar{\gamma}$$\mathbf{r},\gamma$$\widehat{\mathbf{x}}$$\widehat{\gamma}$

Figure 2: Block diagram of batch VAMP for optimization. The calculation of
extrinsic information is performed by the “ext” blocks.

The steps for computing the extrinsic mean vector,
$\widetilde{\mathbf{r}}_{t}$, and the scalar precision,
$\widetilde{\gamma}_{t}$, remain unchanged and they are computed directly
without the need for computing the LMMSE estimate, $\bar{\mathbf{x}}_{t}$, and
the posterior precision, $\bar{\gamma}_{t}$. Hence, the only Kronecker product
required for the LMMSE is of the two vectors $\bm{\omega}_{\textnormal{B}}$
and $\bm{\omega}_{\textnormal{A}}$.

Algorithm 2 Max-sum VAMP SVD for optimization

Given $\mathbf{A}\in\operatorname{\mathbb{C}}^{M\times N}$,
$\mathbf{B}\in\operatorname{\mathbb{C}}^{Q\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: Select initial $\mathbf{r}_{0}$, $\gamma_{0}\geq 0$ and $t\leftarrow 1$

2: Compute economy-size SVD
$\mathbf{A}=\mathbf{U}_{\textnormal{A}}\textnormal{Diag}(\bm{\omega}_{\textnormal{A}})\mathbf{V}_{\textnormal{A}}^{\mathsf{H}}$

3: Compute economy-size SVD
$\mathbf{B}=\mathbf{U}_{\textnormal{B}}\textnormal{Diag}(\bm{\omega}_{\textnormal{B}})\mathbf{V}_{\textnormal{B}}^{\mathsf{H}}$

4: Compute
$\mathbf{V}_{\textnormal{BA}}=\left(\mathbf{V}_{\textnormal{B}}^{\mathsf{H}}\ast\mathbf{V}_{\textnormal{A}}^{\mathsf{H}}\right)^{\mathsf{H}}$

5: Compute normalization vector
$\mathbf{v}_{\textnormal{n}}=\big{[}\left\lVert\mathbf{v}_{\textnormal{BA},1}\right\rVert_{2},\
\left\lVert\mathbf{v}_{\textnormal{BA},2}\right\rVert_{2},\cdots,\left\lVert\mathbf{v}_{\textnormal{BA},MQ}\right\rVert_{2}\big{]}^{\mathsf{T}}$

6: Compute
$\mathbf{V}^{\mathsf{H}}=\mathbf{V}_{\textnormal{BA}}^{\mathsf{H}}\odot\left(\mathbf{v}_{\textnormal{n}}^{-1}\mathbf{1}_{N}^{\mathsf{T}}\right)$

7: Compute
$\bm{\omega}=(\bm{\omega}_{\textnormal{B}}\otimes\bm{\omega}_{\textnormal{A}})\odot\mathbf{v}_{\textnormal{n}}$

8: Compute
$\widetilde{\mathbf{z}}=\textnormal{Diag}(\bm{\omega})^{-1}\textnormal{vec}\left(\mathbf{U}_{\textnormal{A}}^{\mathsf{H}}\mathbf{Z}\mathbf{U}_{\textnormal{B}}^{*}\right)$

9:
$R_{\textnormal{BA}}=\textnormal{Rank}(\mathbf{B}\ast\mathbf{A})=\textnormal{length}(\bm{\omega})$

10: repeat

11: // LMMSE SVD Form.

12:
$\mathbf{d}_{t}=\gamma_{\textnormal{w}}\textnormal{Diag}(\gamma_{\textnormal{w}}\bm{\omega}^{2}+\gamma_{t-1}\mathbf{1}_{R_{\textnormal{BA}}})^{-1}\bm{\omega}^{2}$

13:
$\widetilde{\mathbf{r}}_{t}=\mathbf{r}_{t-1}+\frac{N}{R_{\textnormal{BA}}}\mathbf{V}\textnormal{Diag}\left(\mathbf{d}_{t}/\langle\mathbf{d}_{t}\rangle\right)\left(\widetilde{\mathbf{z}}-\mathbf{V}^{\mathsf{H}}\mathbf{r}_{t-1}\right)$

14:
$\widetilde{\gamma}_{t}=\gamma_{t-1}\left\langle\mathbf{d}_{t}\right\rangle/\left(\frac{N}{R_{\textnormal{BA}}}-\left\langle\mathbf{d}_{t}\right\rangle\right)$

15: // Separable MAP Projector.

16:
$\widehat{\mathbf{x}}_{t}=\mathbf{g}(\widetilde{\mathbf{r}}_{t},\widetilde{\gamma})$

17:
$\widehat{\gamma}_{t}=\widetilde{\gamma}_{t}^{-1}\left\langle\mathbf{g}^{\prime}(\widetilde{\mathbf{r}}_{t},\widetilde{\gamma})\right\rangle$

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

19:
$\mathbf{r}_{t}=\gamma_{t}^{-1}(\widehat{\gamma}_{t}\widehat{\mathbf{x}}_{t}-\widetilde{\gamma}_{t}\widetilde{\mathbf{r}}_{t})$

20: $t\leftarrow t+1$

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

22: return $\widehat{\mathbf{x}}_{t}$

####  III-B2 Scalar MAP Projector

Because the constraint on $\mathbf{x}$ is component-wise, we model the
constraint on its entries, $x_{i}$, as a prior with some precision,
$\gamma_{\textnormal{p}}$, i.e.,
$p_{\mathsf{x}}(x_{i})\varpropto\exp\left(-\gamma_{\textnormal{p}}\left\lvert
f_{i}(x_{i})\right\rvert^{2}\right)$ with $\gamma_{\textnormal{p}}\to\infty$.
We then define the scalar denoising function (now called projector function in
the context of optimization) as follows:

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

or equivalently:

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

The parameter $\gamma_{\textnormal{p}}$ in (30) 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. Taking the derivative of the
scalar projector function w.r.t. $\widetilde{r}_{i,t}$ as defined in equation
(19) yields:

$g_{i}^{\prime}(\widetilde{r}_{i,t},\widetilde{\gamma}_{t})=\widetilde{\gamma}_{t}\widehat{\gamma}_{t},$
(31)

where $\widehat{\gamma}_{t}$ is the posterior precision. The vector valued
projector function,
$\mathbf{g}(\widetilde{\mathbf{r}}_{t},\widetilde{\gamma})$, and its
derivative,
$\mathbf{g}^{\prime}(\widetilde{\mathbf{r}}_{t},\widetilde{\gamma})$, are
defined in the same way as (20) and (21) respectively. Similar to the denoiser
module, extrinsic information from the projector module is calculated in the
form of the mean vector,
$\mathbf{r}_{t}=\left(\widehat{\mathbf{x}}_{t}\widehat{\gamma}_{t}-\widetilde{\mathbf{r}}_{t}\widetilde{\gamma}_{t}\right)/(\widehat{\gamma}_{t}-\widetilde{\gamma}_{t})$,
and scalar precision,
$\gamma_{t}=\widehat{\gamma}_{t}-\widetilde{\gamma}_{t}$, which are then fed
to the LMMSE module. In an analogous way to sum-product VAMP, the max-sum VAMP
(for optimization) decouples the constraint from the objective function and
also enables the projector function to be separable. While the LMMSE module
optimizes the objective function with no constraints, the latter are enforced
by the projector function. This modular property makes VAMP a robust algorithm
for solving optimization problems in the presence of linear mixing and under
various component-wise constraints. The block diagram and the algorithmic
steps for the optimization-oriented VAMP are presented in Fig. 2 and Algorithm
2 respectively.

##  IV VAMP-Based Solution for the Joint Beamforming Problem

In this section, we apply the optimization-oriented VAMP algorithm, described
in Section III-B, to simultaneously optimize the matrix of phase shifters,
${\mathbf{\Upsilon}}$, as well as the optimal precoding matrix $\mathbf{F}$.
We decouple the joint optimization problem into two sub-problems through
alternate optimization. In one side we optimize ${\mathbf{\Upsilon}}$ by
utilizing the modified max-sum VAMP and, on the other side, we find the
optimal transmit precoding $\mathbf{F}$.

###  IV-A Alternating Optimization

We use alternating minimization which is the two-block version of the block
coordinate descent (BCD) algorithm. It is a simple iterative approach that
optimizes one variable at a time222Note here that a variable can be a scalar,
a vector, or a whole matrix. (while fixing the others) and the process is
repeated for every variable. Although it is hard to analytically establish the
optimality of the alternating minimization technique for non-convex
optimization problems, it is known that it performs really well for various
non-convex optimization problems especially for large system sizes [32, 33,
34, 35]. More specifically, we divide the optimization problem in (6) into the
following two sub-optimization problems:

  1. 1.

$\displaystyle\displaystyle\operatorname*{arg\,min}_{{\mathbf{\Upsilon}}}\quad$
$\displaystyle\operatorname{\mathbb{E}}_{\bm{\mathsf{y}},\bm{\mathsf{s}}}\big{\\{}\left\lVert\bm{\mathsf{y}}-\bm{\mathsf{s}}\right\rVert_{\textnormal{2}}^{2}\big{\\}}$
(32a) s.t. $\displaystyle\upsilon_{ik}=0,\quad i\neq k,$ (32b)
$\displaystyle|\upsilon_{ii}|=1,\quad i=1,2,\cdots,K.$ (32c)

  2. 2.

$\displaystyle\operatorname*{arg\,min}_{\alpha,\mathbf{F}}\quad$
$\displaystyle\operatorname{\mathbb{E}}_{\bm{\mathsf{y}},\bm{\mathsf{s}}}\big{\\{}\left\lVert\bm{\mathsf{y}}-\bm{\mathsf{s}}\right\rVert_{\textnormal{2}}^{2}\big{\\}}$
(33a) s.t.
$\displaystyle\operatorname{\mathbb{E}}_{\bm{\mathsf{s}}}\|\mathbf{F}\bm{\mathsf{s}}\|_{\textnormal{2}}^{2}=P.$
(33b)

Let us define the error at iteration $t$ as follows:

$E_{t}~{}\triangleq~{}\left\lVert\widehat{\alpha}_{t}\left(\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\widehat{{\mathbf{\Upsilon}}}_{t}\mathbf{H}_{\textnormal{b-s}}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\right)\widehat{\mathbf{F}}_{t}-\mathbf{I}_{M}\right\rVert_{\textnormal{F}}^{2}+M\widehat{\alpha}_{t}^{2}\sigma_{\textnormal{w}}^{2}.$
(34)

The algorithm stops iterating when $|E_{t}-E_{t-1}|<\epsilon E_{t-1}$, where
$\epsilon\in\operatorname{\mathbb{R}}_{+}$ is some precision tolerance. The
algorithmic steps for alternating minimization (after evaluating the
expectation) are shown in Algorithm 3.

Algorithm 3 Alternating minimization

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

1: Initialize $\widehat{{\mathbf{\Upsilon}}}_{\textnormal{0}}$ and
$t\leftarrow 1$.

2: repeat

3:
$\displaystyle[\widehat{\alpha}_{t},\widehat{\mathbf{F}}_{t}]=\displaystyle\operatorname*{arg\,min}_{\alpha,\mathbf{F}}$
$\displaystyle\quad\left\lVert\alpha\left(\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\widehat{{\mathbf{\Upsilon}}}_{t-1}\mathbf{H}_{\textnormal{b-s}}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\right)\mathbf{F}-\mathbf{I}_{M}\right\rVert_{\textnormal{F}}^{2}+M\alpha^{2}\sigma_{\textnormal{w}}^{2}$
s.t. $\displaystyle\quad\|\mathbf{F}\|_{\textnormal{F}}^{2}=P$

4:
$\displaystyle\widehat{{\mathbf{\Upsilon}}}_{t}=\operatorname*{arg\,min}_{{\mathbf{\Upsilon}}}\quad$
$\displaystyle\left\lVert\widehat{\alpha}_{t}\left(\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}{\mathbf{\Upsilon}}\mathbf{H}_{\textnormal{b-s}}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\right)\widehat{\mathbf{F}}_{t}-\mathbf{I}_{M}\right\rVert_{\textnormal{F}}^{2}$
s.t. $\displaystyle\upsilon_{ik}=0\quad i\neq k,$
$\displaystyle|\upsilon_{ii}|=1\quad i=1,2,\cdots,K$

5: $t\leftarrow t+1$

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

7: return $\widehat{{\mathbf{\Upsilon}}}_{t},\ \widehat{\mathbf{F}}_{t},\
\widehat{\alpha}_{t}$.

###  IV-B Optimization of the Phase Shifters Matrix

Here, we specialize optimization-oriented VAMP algorithm introduced in Section
III-B in order to optimize the phase matrix, ${\mathbf{\Upsilon}}$. Let us
restate the associated optimization after explicitly finding the expectation
in (32) as follows:

$\displaystyle\displaystyle\operatorname*{arg\,min}_{{\bm{\upsilon}}}\quad$
$\displaystyle\left\lVert\alpha\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}({\bm{\upsilon}})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}-(\mathbf{I}_{M}-\alpha\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F})\right\rVert_{\textnormal{F}}^{2}$
(35a) s.t. $\displaystyle|\upsilon_{i}|=1\quad i=1,2,\cdots,K.$ (35b)

The solution is obtained by setting
$\mathbf{A}=\alpha\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}$,
$\mathbf{B}=(\mathbf{H}_{\textnormal{b-s}}\mathbf{F})^{\mathsf{T}}$ and
$\mathbf{Z}=\mathbf{I}_{M}-\alpha\mathbf{H}_{\textnormal{b-u}}\mathbf{F}$ in
Algorithm 2 and then choosing a suitable projector function to satisfy the
constraints on the reflection coefficients. The unconstrained minimization of
the objective function in (35a) is performed by the LMMSE module. We define
the projector function that enforces the constraint on the reflection
coefficients as:

$g_{2,i}\left(\widetilde{r}_{i},\widetilde{\gamma}\right)~{}\triangleq~{}\displaystyle\operatorname*{arg\,min}_{\upsilon_{i}}\left[\widetilde{\gamma}\left\lvert\upsilon_{i}-\widetilde{r}_{i}\right\rvert^{2}+\gamma_{\textnormal{p}}\big{|}\left\lvert\upsilon_{i}\right\rvert-1\big{|}^{2}\right].$
(36)

Solving the optimization problem in (36) results in the following closed-form
expression for the underlying projector function:

$g_{2,i}\left(\widetilde{r}_{i},\widetilde{\gamma}\right)~{}=~{}\frac{\widetilde{\gamma}}{\widetilde{\gamma}+\gamma_{\textnormal{p}}}\widetilde{r}_{i}~{}+~{}\frac{\gamma_{\textnormal{p}}}{\widetilde{\gamma}+\gamma_{\textnormal{p}}}\widetilde{r}_{i}\left\lvert\widetilde{r}_{i}\right\rvert^{-1}.$
(37)

As $\gamma_{\textnormal{p}}\to\infty$, we have,
$\frac{\widetilde{\gamma}}{\widetilde{\gamma}+\gamma_{\textnormal{p}}}\to 0$
and
$\frac{\gamma_{\textnormal{p}}}{\widetilde{\gamma}+\gamma_{\textnormal{p}}}\to
1$. Therefore, the projector function simplifies to:

$g_{2,i}\left(\widetilde{r}_{i}\right)~{}=~{}\widetilde{r}_{i}\left\lvert\widetilde{r}_{i}\right\rvert^{-1}.$
(38)

The derivative of the projector function (38) w.r.t. $\widetilde{r}_{i}$ is
obtained according to equation (19) as follows:

$g^{\prime}_{2,i}\left(\widetilde{r}_{i}\right)=\frac{1}{2}\left\lvert\widetilde{r}_{i}\right\rvert^{-1}.$
(39)

Finally, the projector function, $\mathbf{g}_{2}(\widetilde{\mathbf{r}}_{t})$,
and its derivative $\mathbf{g}_{2}^{\prime}(\widetilde{\mathbf{r}}_{t})$ are
obtained by following (20) and (21) respectively.

SeparableMAPProjectorLMMSEEstimatorextext$\textnormal{Diag}(\widehat{{\bm{\upsilon}}})$$\mathbf{Z}$$\widetilde{\mathbf{r}},\widetilde{\gamma}$$\bar{\gamma}$$\bar{{\bm{\upsilon}}}$$\mathbf{r},\gamma$$\widehat{{\bm{\upsilon}}}$$\widehat{\gamma}$Precoder$\widehat{\mathbf{F}}$$\widehat{\alpha}$

Figure 3: Block diagram of the proposed algorithm.

###  IV-C Optimal Precoding

The sub-optimization problem in (33) is a constrained MMSE transmit precoding
optimization for traditional MIMO systems. It can be solved by jointly
optimizing $\mathbf{F}$ and $\alpha$ using Lagrange optimization. After
finding the expectation, we construct the Lagrangian function associated to
the problem in (33) as follows:

$\operatorname{\mathcal{L}}(\mathbf{F},\alpha,\lambda)=\left\lVert\alpha\mathbf{H}^{\mathsf{H}}\mathbf{F}-\mathbf{I}_{M}\right\rVert_{\textnormal{F}}^{2}+M\alpha^{2}\sigma_{\textnormal{w}}^{2}+\lambda(\textnormal{Tr}(\mathbf{F}\mathbf{F}^{\mathsf{H}})-P),$
(40)

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

$\alpha^{\textnormal{opt}}~{}=~{}g_{3}\left(\mathbf{H}\right)~{}\triangleq~{}\sqrt{\frac{1}{P}}\sqrt{\textnormal{Tr}\left(\left[\mathbf{H}\mathbf{H}^{\mathsf{H}}+\frac{M\sigma_{\textnormal{w}}^{2}\mathbf{I}_{N}}{P}\right]^{-2}\mathbf{H}\mathbf{H}^{\mathsf{H}}\right)}.$
(41)
$\mathbf{F}^{\textnormal{opt}}~{}=~{}g_{4}\left(\mathbf{H}\right)\triangleq\frac{\sqrt{P}\left[\mathbf{H}\mathbf{H}^{\mathsf{H}}+\frac{M\sigma_{\textnormal{w}}^{2}\mathbf{I}_{N}}{P}\right]^{-1}\mathbf{H}}{\sqrt{\textnormal{Tr}\left(\left[\mathbf{H}\mathbf{H}^{\mathsf{H}}+\frac{M\sigma_{\textnormal{w}}^{2}\mathbf{I}_{N}}{P}\right]^{-2}\mathbf{H}\mathbf{H}^{\mathsf{H}}\right)}}~{}=~{}\alpha^{\textnormal{opt}^{-1}}\left[\mathbf{H}\mathbf{H}^{\mathsf{H}}+\frac{M\sigma_{\textnormal{w}}^{2}\mathbf{I}_{N}}{P}\right]^{-1}\mathbf{H}.$
(42)

Note that, the scalar, $\alpha$, merely represents a scaling factor at the
receiver that is used to scale the incident signal as so to obtain the
transmitted constellation symbols and this a common practice in MMSE precoding
optimization [24, 25]. Choosing a common $\alpha$ for all users results in
better tractability and makes it possible to derive a closed-form solution for
the optimal $\alpha$. By inspecting the closed-form solution of the precoding
matrix, we observe that it is scaled by $\alpha^{\textnormal{opt}^{-1}}$. This
allows the transmitter to optimally scale all the transmit symbols based on
the available transmit power whereas the receiver upscales the received signal
plus noise to get back the original transmitted symbols while keeping the SNR
unaffected. It is also worth mentioning that the optimal scaling factor,
$\alpha^{\textnormal{opt}}$, is only utilized in optimizing the precoding
matrix since the receivers can blindly estimate this scalar based on the
received symbol sequence [24, 25].  Now that we have solved both sub-
optimization problems in (32) and (33), separately, we substitute their
solutions into Algorithm 3. The overall block diagram and algorithmic steps
are respectively shown in Fig. 3 and Algorithm 4.

######  Remark.

It is possible to include per-user data requirement by formulating the problem
under the weighted MMSE (WMMSE) criterion where we scale the MSEs of the users
with weights according to each user’s data requirement and then minimize the
sum MSE. To that end, we define a positive semi-definite real diagonal matrix,
$\mathbf{Q}$, containing user weights, $\\{q_{m}\\}_{m=1}^{M}$, in its
diagonal, i.e., $\mathbf{Q}=\textnormal{Diag}(q_{1},\cdots,q_{M})$. Then, the
optimization problem under the WMMSE criterion is given by:

$\displaystyle\displaystyle\operatorname*{arg\,min}_{\alpha,\mathbf{F},{\mathbf{\Upsilon}}}\quad$
$\displaystyle\operatorname{\mathbb{E}}_{\bm{\mathsf{y}},\bm{\mathsf{s}}}\left\\{\left\lVert\mathbf{Q}^{1/2}(\bm{\mathsf{y}}-\bm{\mathsf{s}})\right\rVert_{2}^{2}\right\\},$
(43a) subject to
$\displaystyle\operatorname{\mathbb{E}}_{\bm{\mathsf{s}}}\left\\{\|\mathbf{F}\bm{\mathsf{s}}\|_{2}^{2}\right\\}=P,$
(43b) $\displaystyle\upsilon_{ik}=0,\quad i\neq k,$ (43c)
$\displaystyle|\upsilon_{ii}|=1,\quad i=1,2,\cdots,K.$ (43d)

In this formulation, WMMSE precoding optimization is performed instead of the
ordinary MMSE precoding optimization, wherein the matrices $\mathbf{A}$,
$\mathbf{B}$ and $\mathbf{Z}$ are adjusted accordingly inside the VAMP part of
the Algorithm 4.

Algorithm 4 VAMP-based joint optimization algorithm

Given $\mathbf{H}_{\textnormal{s-u}}$, $\mathbf{H}_{\textnormal{b-u}}$,
$\mathbf{H}_{\textnormal{b-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$

2:
$\widehat{\mathbf{H}}_{0}=\left(\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\widehat{{\bm{\upsilon}}}_{\textnormal{0}})\mathbf{H}_{\textnormal{b-s}}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\right)^{\mathsf{H}}$

3:
$\widehat{\alpha}_{\textnormal{0}}=g_{3}\left(\widehat{\mathbf{H}}_{0}\right)$

4:
$\widehat{\mathbf{F}}_{\textnormal{0}}=g_{4}\left(\widehat{\mathbf{H}}_{0}\right)$

5: repeat

6: // LMMSE SVD Form.

7: Set
$\mathbf{A}=\widehat{\alpha}_{t-1}\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}$,
$\mathbf{B}=\left(\mathbf{H}_{\textnormal{b-s}}\widehat{\mathbf{F}}_{t-1}\right)^{\mathsf{T}}$
and
$\mathbf{Z}=\mathbf{I}_{M}-\widehat{\alpha}_{t-1}\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\widehat{\mathbf{F}}_{t-1}$.

8: Compute economy-size SVD
$\mathbf{A}=\mathbf{U}_{\textnormal{A}}\textnormal{Diag}(\bm{\omega}_{\textnormal{A}})\mathbf{V}_{\textnormal{A}}^{\mathsf{H}}$

9: Compute economy-size SVD
$\mathbf{B}=\mathbf{U}_{\textnormal{B}}\textnormal{Diag}(\bm{\omega}_{\textnormal{B}})\mathbf{V}_{\textnormal{B}}^{\mathsf{H}}$

10: Compute
$\mathbf{V}_{\textnormal{BA}}=\left(\mathbf{V}_{\textnormal{B}}^{\mathsf{H}}\ast\mathbf{V}_{\textnormal{A}}^{\mathsf{H}}\right)^{\mathsf{H}}$

11: Compute normalization vector
$\mathbf{v}_{\textnormal{n}}=\big{[}\left\lVert\mathbf{v}_{\textnormal{BA},1}\right\rVert_{2},\
\left\lVert\mathbf{v}_{\textnormal{BA},2}\right\rVert_{2},\cdots,\left\lVert\mathbf{v}_{\textnormal{BA},M^{2}}\right\rVert_{2}\big{]}^{\mathsf{T}}$

12: Compute
$\mathbf{V}^{\mathsf{H}}=\mathbf{V}_{\textnormal{BA}}^{\mathsf{H}}\odot\left(\mathbf{v}_{\textnormal{n}}^{-1}\mathbf{1}_{K}^{\mathsf{T}}\right)$

13: Compute
$\bm{\omega}=(\bm{\omega}_{\textnormal{B}}\otimes\bm{\omega}_{\textnormal{A}})\odot\mathbf{v}_{\textnormal{n}}$

14: Compute
$\widetilde{\mathbf{z}}={\bm{\omega}}^{-1}\odot\textnormal{vec}\left(\mathbf{U}_{\textnormal{A}}^{\mathsf{H}}\mathbf{Z}\mathbf{U}_{\textnormal{B}}^{*}\right)$

15:
$R_{\textnormal{BA}}=\textnormal{Rank}(\mathbf{B}\ast\mathbf{A})=\textnormal{length}(\bm{\omega})$

16:
$\mathbf{d}_{t}=\gamma_{\textnormal{w}}(\gamma_{\textnormal{w}}\bm{\omega}^{2}+\gamma_{t-1}{\mathbf{1}_{R_{\textnormal{BA}}}})^{-1}\odot\bm{\omega}^{2}$

17:
$\widetilde{\mathbf{r}}_{t}=\mathbf{r}_{t-1}+\frac{K}{R_{\textnormal{BA}}}\mathbf{V}\left(\frac{\mathbf{d}_{t}}{\langle\mathbf{d}_{t}\rangle}\odot\left(\widetilde{\mathbf{z}}-\mathbf{V}^{\mathsf{H}}\mathbf{r}_{t-1}\right)\right)$

18:
$\widetilde{\gamma}_{t}=\gamma_{t-1}\left\langle\mathbf{d}_{t}\right\rangle/\left(\frac{K}{R_{\textnormal{BA}}}-\left\langle\mathbf{d}_{t}\right\rangle\right)$

19: // Separable MAP Projector

20:
$\widehat{{\bm{\upsilon}}}_{t}=\mathbf{g}_{2}\left(\widetilde{\mathbf{r}}_{t}\right)$

21:
$\widehat{\gamma}_{t}=\widetilde{\gamma}_{t}^{-1}\left\langle\mathbf{g}^{\prime}_{2}\left(\widetilde{\mathbf{r}}_{t}\right)\right\rangle$.

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

23:
$\mathbf{r}_{t}=\gamma_{t}^{-1}\left(\widehat{\gamma}_{t}\widehat{{\bm{\upsilon}}}_{t}-\widetilde{\gamma}_{t}\widetilde{\mathbf{r}}_{t}\right)$

24: //Find $\alpha$ and $\mathbf{F}$ through their closed-form solutions.

25:
$\widehat{\mathbf{H}}_{t}=\left(\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}(\widehat{{\bm{\upsilon}}}_{t})\mathbf{H}_{\textnormal{b-s}}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\right)^{\mathsf{H}}$

26: $\widehat{\alpha}_{t}=g_{3}\left(\widehat{\mathbf{H}}_{t}\right)$

27: $\widehat{\mathbf{F}}_{t}=g_{4}\left(\widehat{\mathbf{H}}_{t}\right)$

28: $t\leftarrow t+1$

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

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

##  V Joint Beamforming Under Reactive Loading at the IRS

We consider a reflective element that is combined with a tunable reactive load
333This 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 shifter, i.e., 444The value
$1$ is the normalized resistive part of the element impedance whereas
$\chi_{i}$ is the normalized reactive part of the antenna plus reactive
termination. Accordingly $\upsilon_{i}$ represents the induced current flowing
across the antenna. We assume the antenna elements to be uncoupled which holds
approximately for half-wavelength element spacing. $\upsilon_{i}=-(1+{\rm
j}\chi_{i})^{-1}$, where $\chi_{i}\in\operatorname{\mathbb{R}}$ is a scalar
reactance value that has to be optimized for each reflection coefficient.
Under the unimodular constraint, the idealistic IRS has a full field of view
(FOV) and the reflection coefficients correspond to ideal phase-shifters and
are of the form $\upsilon_{i}=e^{{\rm j}\theta_{i}}$, where
$\theta_{i}\in[0,2\pi]$555In practice, this assumption is difficult from a
practical standpoint. With the reactive-loading constraint, the assumption of
an IRS with full FOV becomes more acceptable., whereas under the practical
constraint we have a restriction on the possible values of the IRS 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. We rewrite the objective
function under the new constraint on phases as follows:

$\displaystyle\displaystyle\operatorname*{arg\,min}_{{\bm{\upsilon}}}\quad$
$\displaystyle\left\lVert\alpha\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}\textnormal{Diag}({\bm{\upsilon}})\mathbf{H}_{\textnormal{b-s}}\mathbf{F}-(\mathbf{I}_{M}-\alpha\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F})\right\rVert_{\textnormal{F}}^{2}$
(44a) s.t. $\displaystyle\upsilon_{i}=\frac{-1}{1+{\rm j}\chi_{i}},\quad
i=1,2,\cdots,K.$ (44b)

To find the sub-optimal phase matrix under the new constraint, we change the
projector function accordingly as follows:

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

where

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

The optimization problem in (45) is a bi-level one [36]. The solution to (46)
is substituted in (45) which is then solved as ordinary MAP optimization. We
show in Appendix A that the solution to (46) is given by:

$g_{6}\left(\widetilde{r}_{i}\right)~{}=~{}\frac{1}{2\Im\left\\{\widetilde{r}_{i}\right\\}}\left(1+2\Re\left\\{\widetilde{r}_{i}\right\\}+\sqrt{\left(1+2\Re\left\\{\widetilde{r}_{i}\right\\}\right)^{2}+4\Im\left\\{\widetilde{r}_{i}\right\\}^{2}}\right).$
(47)

Substituting (47) back into (45) and solving the minimization leads to the
following result:

$g_{5,i}(\widetilde{r}_{i},\widetilde{\gamma})=\frac{\widetilde{\gamma}}{\widetilde{\gamma}+\gamma_{\textnormal{p}}}\widetilde{r}_{i}-\frac{\gamma_{\textnormal{p}}}{\widetilde{\gamma}+\gamma_{\textnormal{p}}}\left(1+{\rm
j}g_{6}\left(\widetilde{r}_{i}\right)\right)^{-1},$ (48)

where $\gamma_{\textnormal{p}}\to\infty.$ Thus, the projector function can be
expressed as:

$g_{5,i}\left(\widetilde{r}_{i}\right)=-\left(1+{\rm
j}g_{6}\left(\widetilde{r}_{i}\right)\right)^{-1},$ (49)

whose derivative is obtained as defined in equation (19) as follows:

$g^{\prime}_{5,i}\left(\widetilde{r}_{i}\right)={\rm
j}g^{\prime}_{6}\left(\widetilde{r}_{i}\right)\left(1+{\rm
j}g_{6}\left(\widetilde{r}_{i}\right)\right)^{-2},$ (50)

where

$g^{\prime}_{6}\left(\widetilde{r}_{i}\right)=\frac{1}{2}\left(\frac{\partial
g_{6}\left(\widetilde{r}_{i}\right)}{\partial\Re\left\\{\widetilde{r}_{i}\right\\}}-{\rm
j}\frac{\partial
g_{6}\left(\widetilde{r}_{i}\right)}{\partial\Im\left\\{\widetilde{r}_{i}\right\\}}\right).$
(51)

The partial derivatives involved in (51) are given by:

$\frac{\partial
g_{6}\left(\widetilde{r}_{i}\right)}{\partial\Re\left\\{\widetilde{r}_{i}\right\\}}=\Im\left\\{\widetilde{r}_{i}\right\\}^{-1}+\left(1+2\Re\left\\{\widetilde{r}_{i}\right\\}\right)\left(\Im\left\\{\widetilde{r}_{i}\right\\}\sqrt{\left(1+2\Re\left\\{\widetilde{r}_{i}\right\\}\right)^{2}+4\Im\left\\{\widetilde{r}_{i}\right\\}^{2}}\right)^{-1},$
(52)

and

$\frac{\partial
g_{6}\left(\widetilde{r}_{i}\right)}{\partial\Im\left\\{\widetilde{r}_{i}\right\\}}=-\left(1+2\Re\left\\{\widetilde{r}_{i}\right\\}\right)\left(2\Im\left\\{\widetilde{r}_{i}\right\\}^{2}\right)^{-1}\\\
-\left(1+2\Re\left\\{\widetilde{r}_{i}\right\\}\right)^{2}\left(2\Im\left\\{\widetilde{r}_{i}\right\\}^{2}\sqrt{\left(1+2\Re\left\\{\widetilde{r}_{i}\right\\}\right)^{2}+4\Im\left\\{\widetilde{r}_{i}\right\\}^{2}}\right)^{-1}.$
(53)

Since the derivative is required to be a real scalar, we take the absolute
value of the complex derivative and, therefore, we modify the derivative of
the projector function (50) as follows:

$g^{\prime}_{5,i}\left(\widetilde{r}_{i}\right)=\left\lvert{\rm
j}g^{\prime}_{6}\left(\widetilde{r}_{i}\right)\left(1+{\rm
j}g_{6}\left(\widetilde{r}_{i}\right)\right)^{-2}\right\rvert.$ (54)

Lastly, we obtain the vector valued projector function,
$\mathbf{g}_{5}(\widetilde{\mathbf{r}}_{t})$, and its derivative
$\mathbf{g}_{5}^{\prime}(\widetilde{\mathbf{r}}_{t})$ according to (20) and
(21) respectively and replace $\mathbf{g}_{2}(\widetilde{\mathbf{r}}_{t})$ and
$\mathbf{g}_{2}^{\prime}(\widetilde{\mathbf{r}}_{t})$ in lines $19$ and $20$
of Algorithm 4.

##  VI Numerical Results: Performance and Complexity Analysis

###  VI-A Simulation Model and Parameters

We present exhaustive Monte-Carlo simulation results to assess the performance
of the proposed algorithm. We assume that the IRS 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 IRS. A path-based propagation channel
model, also known as parametric channel model [26], is used. Such a model are
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-IRS
channel which governs the effect of channel correlation.

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

The channel between the IRS and the BS is generated according to:

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

Here, $Q_{\textnormal{IRS}}$ and $L(d_{\textnormal{IRS}})$ denote the number
of channel paths and the distance-dependent path-loss factor, respectively.
The vectors $\mathbf{a}_{\textnormal{BS}}(\phi)$ and
$\mathbf{a}_{\textnormal{IRS}}(\varphi,\psi)$ are the array response vectors
for the BS and the IRS, respectively. The coefficients $c_{q}$ in (55) denote
the path gains which are modeled by a complex normal distribution, i.e.,
$\mathsf{c}_{q}\sim\operatorname{\mathcal{C}}\operatorname{\mathcal{N}}(c_{q};0,1)$.
Assuming that a uniform linear array (ULA) with $N$ antennas is used at the
BS, we have $\mathbf{a}_{\textnormal{BS}}(\phi)=[1,e^{2\pi{\rm
j}\frac{d_{\textnormal{b}}}{\lambda}\cos\phi},\cdots,e^{2\pi{\rm
j}\frac{d_{\textnormal{b}}}{\lambda}(N-1)\cos\phi}]^{\mathsf{T}}$ wherein
$\lambda$, $\phi$, and $d_{\textnormal{b}}$ represent the wavelength, the
angle of departure (AOD), and the inter-antenna spacing at the BS,
respectively. The IRS is equipped with a (square) uniform planar array (UPA)
with $K$ antenna elements which are assumed to have a cosine embedded element
pattern. By defining the z-axis as the normal vector to the array, the array
response vector for the IRS is expressed as follows [37]:

$\mathbf{a}_{\textnormal{IRS}}(\varphi,\psi)~{}=~{}\sqrt{|\cos\varphi|}\begin{bmatrix}1\\\
e^{2\pi{\rm j}\frac{d_{\textnormal{s}}}{\lambda}\sin\varphi\sin\psi}\\\
\vdots\\\ e^{2\pi{\rm
j}\frac{d_{\textnormal{s}}}{\lambda}(\sqrt{K}-1)\sin\varphi\sin\psi}\end{bmatrix}\otimes\begin{bmatrix}1\\\
e^{2\pi{\rm j}\frac{d_{\textnormal{s}}}{\lambda}\sin\varphi\cos\psi}\\\
\vdots\\\ e^{2\pi{\rm
j}\frac{d_{\textnormal{s}}}{\lambda}(\sqrt{K}-1)\sin\varphi\cos\psi}\end{bmatrix}.$
(56)

Here $d_{\textnormal{s}}$ represents the inter-antenna spacing at the IRS
whereas $\varphi$ and $\psi$ are the angles of elevation and azimuth,
respectively. In simulations we set
$d_{\textnormal{b}}=d_{\textnormal{s}}=\lambda/2$. The angles $\psi_{q}$ and
$\phi_{q}$ are uniformly distributed in the interval $[0,2\pi)$ and the
$\varphi_{q}$’s are uniformly distributed in $[0,\pi)$. The channel of the
direct link between the BS and each $m$-th single-antenna user, with
$Q_{\textnormal{b-u}}$ paths, is modeled as follows:

$\mathbf{h}_{\textnormal{b-u},m}~{}=~{}\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.$ (57)

Similar to the IRS-BS channel,
$\mathsf{c}_{m,q}\sim\operatorname{\mathcal{C}}\operatorname{\mathcal{N}}(c_{m,q};0,1)$
and each angle $\phi_{m,q}$ is uniformly distributed in $[0,2\pi)$. The
channel vectors in (57) are assumed to be independent across all users.
Finally, the channel vector for the link between each $m$-th user, and the IRS
with $Q_{\textnormal{s-u}}$ channel paths, is modeled as follows:

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

The term $L(d)=C_{\textnormal{0}}(d/d_{\textnormal{0}})^{\eta}$ in (55), (57),
(58) 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. Moreover, to account for the line-of-sight
(LOS) component, the gain of one channel path is set to at least of $5$ dB
higher than the other path gains. To account for channel correlation effects,
we have set the number of multi-path components lower than the number of BS
antennas and the IRS antenna elements for the BS-IRS channel thereby making
the channel matrix rank-deficient. Therefore, in simulations we have set the
number of users lower than the rank of BS-IRS channel matrix
$\mathbf{H}_{\textnormal{b-s}}$. In the simulations, we fix
$d_{\textnormal{IRS}}=500$ m for the IRS-BS channel whereas the user-BS
distance, $d$, and the user-IRS distance, $d^{\prime}$, vary for each user
according to its location from the BS and the IRS, respectively. In all
simulations, we also set $C_{\textnormal{0}}=-30\ \textnormal{dB}$, $\eta=3.7$
(NLOS BS-user channel), $\eta=2.5$ (NLOS IRS-BS and IRS-user channels),
$Q_{\textnormal{b-u}}=2,\ Q_{\textnormal{s-u}}=2,\ \epsilon=10^{-3}$ and
$\sigma_{\textnormal{w}}^{2}=-100$ dBm. The results are averaged over $1000$
independent Monte Carlo simulations.

The following two scenarios are studied. First, we consider the case where
only the BS-IRS channel contains a LOS component. Then we consider the
scenario where both the BS-IRS and the IRS-user channels have a LOS component
but all the direct BS-user channels do not have a LOS component. The proposed
VAMP-based algorithm is compared against the following four different
configurations:

  * i.

A MIMO system assisted by one IRS where the SDR technique is used to optimize
the IRS reflection coefficients in combination with MMSE precoding.

  * ii.

A MIMO system assisted by one IRS where the joint optimization of the phase
matrix the and the precoding is solved through alternate optimization and
penalty-based ADMM.

  * iii.

A massive MIMO system with a large number of BS antennas with MMSE precoding.

  * iv.

An IRS-assisted MIMO system with unoptimized IRS phases and MMSE transmit
precoding.

###  VI-B Benchmarking Metrics

We use two metrics for performance evaluation, namely the sum-rate,
$\widehat{C}$, and the normalized root mean square error (NRMSE) which are
defined as follows:

$\widehat{C}~{}=~{}\displaystyle\sum_{m=1}^{M}\log_{2}\left(1+\frac{\left\lvert\mathbf{h}_{m}^{\mathsf{H}}\mathbf{f}_{m}\right\rvert^{2}}{\sigma_{\textnormal{w}}^{2}+\displaystyle\sum_{i\neq
m}\left\lvert\mathbf{h}_{m}^{\mathsf{H}}\mathbf{f}_{i}\right\rvert^{2}}\right),$
(59)

where
$\mathbf{h}_{m}^{H}=\mathbf{h}_{\textnormal{s-u},m}^{\mathsf{H}}{\mathbf{\Upsilon}}\mathbf{H}_{\textnormal{b-s}}+\mathbf{h}_{\textnormal{b-u},m}^{\mathsf{H}}$.

$\textnormal{NRMSE}(\alpha,{\mathbf{\Upsilon}},\mathbf{F})~{}\triangleq~{}\frac{1}{\sqrt{M}}\sqrt{\left\lVert\alpha\left(\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}{\mathbf{\Upsilon}}\mathbf{H}_{\textnormal{b-s}}+\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\right)\mathbf{F}-\mathbf{I}_{M}\right\rVert_{\textnormal{F}}^{2}+M\alpha^{2}\sigma_{\textnormal{w}}^{2}}.$
(60)

###  VI-C Performance Results With Perfect CSI

####  VI-C1 BS-IRS channel with LOS component

This situation is encountered in a typical urban or suburban environments
where the BS is located far away from the users and has no direct LOS
component. However, the IRS is installed at a location where a LOS component
is present in the BS-IRS link but not in the user-IRS link. Here we set the
number of users to $M=8$ and the number of BS antennas to $N=32$ for every
configuration except for massive MIMO for which we use $N=96$. Fig. 4(a),
depicts the achievable sum-rate versus the transmit power, $P$, for the
different considered transmission schemes. The proposed algorithm in this
scenario outperforms the massive MIMO system even with a significantly smaller
number of transmit antennas. VAMP automatically updates the stepsize at a per-
iteration basis that leads to a faster convergence compared to other iterative
algorithms Since the proposed algorithm is based on VAMP, it outperforms the
ADMM-based solution as it automatically updates the stepsize at a per-
iteration basis (by means of calculating extrinsic information at each step)
that leads to a faster convergence compared to ADMM, where the penalty
parameter must be manually chosen. As per the IRS-assisted configuration,
where one uses the SDR technique to optimize the IRS reflection coefficients,
a significant gap is observed between the achieved sum-rates as compared to
the proposed algorithm.

(a) $M=8$ and $K=256$.

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

(c) $M=8,\ N=32$, and $K=256$.

Figure 4: LOS IRS-BS channel: (a) Sum-rate versus transmit power, (b) Sum-rate
versus the number of IRS reflective elements, and (c) NRMSE versus the number
of iterations (BS-user link excluded).

Fig. 4(b) shows a plot of sum-rate against the number of IRS reflective
elements. It is observed that even with a small number of active transmit
antennas and merely ten paths between the IRS and the BS, the sum-rate for the
proposed solution keeps increasing with the number of reflective elements. In
contrast, the sum-rate saturates after a small gain when the IRS reflection
coefficients are not optimized. Compared to the ADMM-based solution and the
SDR method, the proposed algorithm shows higher throughput at every point.

The convergence of the proposed algorithm is investigated in Fig. 4(c) which
depicts the NRMSE as a function of the number of iterations. Here we exclude
the direct BS-user link to highlight the throughput of the BS-IRS-user link
after optimizing the IRS phase shifts. Observe that the major portion of the
gain is achieved in the first few iterations. The small number of iterations
required for convergence in combination with the low per-iteration complexity
makes the proposed algorithm very attractive from the practical implementation
point of view. The superiority of the proposed VAMP-based algorithm over the
ADMM-based approach stems from the feedback mechanism of VAMP. In fact, such
feedback controls the weight given to the update of $\bm{\Upsilon}$ at each
iteration compared to that of the preceding iteration. This is achieved with
the help of scalar precision parameters that act as weighting coefficients for
the $\bm{\Upsilon}$ updates that are computed in the current and the preceding
iteration. In addition to the plots shifting downward, the increase in
transmit power widens the gap between ADMM and the proposed VAMP-based
algorithm. This demonstrates that the latter utilizes the available transmit
power in a more efficient way than ADMM.

(a) $M=8$ and $K=256$.

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

(c) $M=8,\ N=32$, and $K=256$.

Figure 5: LOS IRS-user and BS-IRS channels: (a) Sum-rate versus transmit
power, (b) Sum-rate versus the number of IRS reflection elements, and (c)
NRMSE versus the number of iterations (BS-user link excluded).

####  VI-C2 BS-IRS and IRS-user channels with LOS components

Fig. 5(a) illustrates the sum-rate versus the transmit power for this
configuration. As expected, the results show that by adding a LOS component,
the use of an IRS together with the proposed joint beamforming optimization
solution yields considerably higher sum-rates compared to a massive MIMO
system with no IRS. Moreover, although the ADMM-based solution now matches the
performance of massive MIMO, the advantage of the proposed VAMP-based solution
over the ADMM- and SDR-based solutions is higher when compared to the NLOS
configuration.

The results in Fig. 5(b), i.e., sum-rate vs the number of IRS reflective
elements, also exhibit the same trends as in the NLOS scenario yet with a
broader gap between the curves, thereby, corroborating the superiority of the
proposed solution. Intuitively, the presence of a LOS component helps the
VAMP-based joint beamforming scheme to focus most of the transmit/reflected
energy in that direction. This is clearly depicted in Fig. 5(c), where the
NRMSE achieved by the proposed algorithm is approaching the NRMSE achieved by
the ADMM-based solution but at almost $5$ dB lower transmit power.

####  VI-C3 Practical IRS phase shifts

In this subsection, we assess the effect of replacing the unimodular
constraint on the reflection coefficients by reactively loaded omni-
directional elements. We use the same channel configuration as in Section
VI-C2. But, we rely on optimizing just the reactive part of the reflection
coefficients. Therefore, as portrayed by Fig. 6(b), the new constraint
decreases the throughput when compared with the ideal phase shifters setup.
However, the resulting sum-rate is still much higher than the one obtained by
using unoptimized IRS reflection coefficients. In fact, when the number of IRS
elements is higher than a certain value, the proposed approach with practical
phase shifts achieves higher throughput than both the SDR- and ADMM-based
solutions with ideal phase shifts. Similarly, due to having less room for
optimizing the reflection coefficients, Fig. 6(a) shows that the NRMSE
saturates sooner and at a higher value as compared to the case of a unimodular
constraint (i.e., ideal phase shifts). Nonetheless, even with the more
practical reactive load constraint, the resulting VAMP-based NRMSE is close to
the NRMSE achieved by ADMM with ideal phase shifters.

(a) $M=8,\ N=32$, $K=256$
and $P=30$ dBm.

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

(c) $M=8,\ N=32$,
and $K=256$.

Figure 6: LOS IRS-user and BS-IRS channels: (a) NRMSE versus the number of
iterations (BS-user link excluded), (b) Sum-rate versus IRS elements with
practical phase shifts, and (c) Sum-rate versus transmit power under imperfect
CSI.

###  VI-D Performance Results With Imperfect CSI

In this section, we measure the performance of the proposed solution 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-IRS-users and the direct BS-user channels [38, 39].
We model the estimated channel matrix and vectors using the statistical CSI
error model proposed in [40, 41, 42] as follows:

$\widehat{\mathbf{H}}_{\textnormal{b-s}}=\kappa\mathbf{H}_{\textnormal{b-s}}+\sqrt{(1-\kappa^{2})L(d_{\textnormal{IRS}})}{\bm{\Delta}}_{\textnormal{b-s}},$
(61)
$\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,$ (62)
$\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,$ (63)

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 matrices $\mathbf{F}$ and ${\bm{\Upsilon}}$ under
imperfect CSI and then use the exact CSI matrices to calculate the sum-rate.
Fig. 6(c) plots the sum-rate versus transmit power for different values of the
channel estimation accuracy parameter $\kappa$. We also include plots for the
other beamforming schemes under perfect CSI for reference. The results show
the resilience of the proposed VAMP-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.85$ performs better than the SDR
based approach and nearly as good as the ADMM-based approach under perfect
CSIs. Moreover, the performance loss with a high channel estimation accuracy
value of $\kappa=0.99$ is negligible.

##  VII Convergence, Optimality, and Complexity Analysis

According to the monotone convergence theorem in real analysis [43], a
monotonically decreasing sequence with a lower bound is convergent. In our
case, since the objective function,

$\left\lVert\alpha\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}{\bm{\Upsilon}}\mathbf{H}_{\textnormal{b-s}}\mathbf{F}-(\mathbf{I}_{M}-\alpha\mathbf{H}_{\textnormal{b-u}}^{\mathsf{H}}\mathbf{F})\right\rVert_{\textnormal{F}}^{2}+M\alpha^{2}\sigma_{\textnormal{w}}^{2}$
(64)

has a lower bound of zero, the proposed algorithm will always converge to a
solution if the MSE monotonically decreases in both steps of the algorithm,
i.e., the step of optimizing the phase-shifts (VAMP part) and the step of
optimizing the active precoding. For the latter, we have a closed-form optimal
solution. Therefore, it is necessary that the MSE decreases monotonically
inside the VAMP step in every iteration to guarantee the convergence of the
entire algorithm. In practice, most of the approximate message passing-based
algorithms (including VAMP) add damping steps inside the algorithm to avoid
any oscillations in the resultant MSE and thus, ensuring convergence [22]. The
lines 18 and 20 inside the VAMP part of the Algorithm 4 are, respectively,
replaced by the damped versions:

$\widetilde{\gamma}_{t}=\varrho\gamma_{t-1}\left\langle\mathbf{d}_{t}\right\rangle/\left(\frac{K}{R_{\textnormal{\tiny
BA}}}-\left\langle\mathbf{d}_{t}\right\rangle\right)+(1-\varrho)\widetilde{\gamma}_{t-1}.$
(65)
$\widehat{\bm{\upsilon}}_{t}=\varrho\mathbf{g}_{1}(\widetilde{\mathbf{r}}_{t})+(1-\varrho)\widehat{{\bm{\upsilon}}}_{t-1},$
(66)

for all iterations $t>1$ where $\varrho\in(0,1]$ is a suitably chosen damping
factor.

TABLE II: Comparison between the CPU execution time of the proposed VAMP-based algorithm, the ADMM-based algorithm and the SDR-based algorithm for different design configurations. The algorithms terminate when $\left\lvert\textnormal{NRMSE}_{t}-\textnormal{NRMSE}_{t-1}\right\rvert<10^{-3}\textnormal{NRMSE}_{t-1}$ or $t>100$. Design Parameters |  VAMP-based algorithm $\operatorname{\mathcal{O}}(MN(K+N))$ (msec) |  ADMM-based algorithm $\operatorname{\mathcal{O}}(MN(K+N))$ (msec) |  SDR-based algorithm $\operatorname{\mathcal{O}}(MN+K^{6})$ (msec)
---|---|---|---
$M=2,\ N=16,\hskip 7.0ptK=64$ |  $14$ |  $26$ |  $2100$
$M=4,\ N=32,\hskip 7.0ptK=256$ |  $104$ |  $340$ |  $12500$

The optimality of the proposed VAMP-based approach can be investigated through
statistical state evolution analysis of the proposed algorithm which we have
left for a future work. Please note that for non convex optimization problems
like optimizing the phase-shifts matrix under uni-modular constraint,
asymptotic (for large matrix sizes) optimality can be claimed for i.i.d.
matrices, if the proximal functions (projector function and its derivative
inside Algorithm 4) are shown to be Lipschitz continuous, and the state
evolution analysis reveals that the VAMP-based algorithm has only one fixed
point [22, 28]. For implementation purpose, we choose the maximum possible
value for precision tolerance, $\epsilon$, for which the mean square error
(MSE) approximately saturates before the algorithm is stopped. For the
proposed solution we have found out that $\epsilon=10^{-3}$ does the trick as
the MSE achieved by choosing any lower values for $\epsilon$ is approximately
equal to the MSE achieved by choosing $\epsilon=10^{-3}$. The maximum number
of iterations, $T_{\textnormal{max}}$, is a hardware-dependent parameter and
is manually chosen to have a limit on the number of iterations.

Note that, by utilizing the Kronecker structure, we avoid any large matrix
multiplication or even taking SVD of Kronecker or Khatri-Rao product of
matrices. Let $\mathbf{A}=\alpha\mathbf{H}_{\textnormal{s-u}}^{\mathsf{H}}$
and $\mathbf{B}=(\mathbf{H}_{\textnormal{b-u}}\mathbf{F})^{\mathsf{T}}$. For
our system model, the matrices $\mathbf{A}$ and $\mathbf{B}$ are of the same
size $M\times K$. Assuming that the matrices $\mathbf{A}$ and $\mathbf{B}$ are
of full rank, the complexity of the truncated SVDs of the matrices is of
$\operatorname{\mathcal{O}}(M^{2}K)$. The computational complexity of the
column-wise Khatri-Rao product in line $10$ and the following operations in
lines $11$ and $12$ of Algorithm 4 has a complexity of
$\operatorname{\mathcal{O}}(M^{2}K)$. The Kronecker product of two vectors in
line $13$ and the component-wise operations of vectors in lines $16$ and $17$
are of order $\operatorname{\mathcal{O}}(M^{2})$. The projector function and
its derivative has a complexity in the order of
$\operatorname{\mathcal{O}}(K)$. The functions $g_{3}\left(\mathbf{H}\right)$
and $g_{4}\left(\mathbf{H}\right)$ can be implemented efficiently by using the
matrix inversion lemma, thereby entailing a complexity of
$\operatorname{\mathcal{O}}(M^{3}+MN^{2})$. The complexity of all other matrix
multiplications elsewhere including the LMMSE part is of order
$\operatorname{\mathcal{O}}(MNK+M^{2}K)$. Therefore, the overall per-iteration
complexity of the algorithm is of order
$\operatorname{\mathcal{O}}(M^{3}+M^{2}K+MNK+MN^{2})$. Since $M<N$ and $M<K$
in our case, the overall per-iteration complexity simplifies to
$\operatorname{\mathcal{O}}(MN(K+N))$ or $\operatorname{\mathcal{O}}(MNK)$ for
$N\leq K$.

Table II provides a comparison of CPU (central processing unit) run time
between the VAMP-based approach, the ADMM-based approach and the SDR-based
approach for different design configurations. For comparison purpose, we
measure the time until the NRMSE saturates with a tolerance, $\epsilon$ ,
therefore we run the algorithms until
$\left\lvert\textnormal{NRMSE}_{t}-\textnormal{NRMSE}_{t-1}\right\rvert<\epsilon\textnormal{NRMSE}_{t-1}$
or $t>T_{\textnormal{MAX}}$, while setting $\epsilon=10^{-3}$ and
$T_{\textnormal{MAX}}=100$. We set the channel simulation parameters as in
Section VI-C2 with $P=30$ dBm. The algorithms are simulated on MATLAB R2020a
on a laptop PC having a Core i7-4720HQ processor and 16 GB of RAM with Windows
10 operating system. As expected, the simulation results confirm that the
proposed approach is significantly faster in terms of convergence time,
especially when there is a high number of IRS antenna elements.

##  VIII Conclusion

We investigated the problem of joint active and passive beamforming design for
an IRS-assisted downlink multi-user MIMO system under both ideal and practical
models for the IRS phase shifts. The associated joint non-convex optimization
has been formulated under sum-MMSE criterion. Using alternating minimization,
the joint optimization has been decomposed into two sub-optimization tasks,
i.e., optimizing the IRS phase shifts and the BS precoding matrix separately.
Regarding the phase shifts, we have presented a novel approach that relies on
the approximate message passing framework to solve the associated sub-
optimization problem. For this, we have first extended the traditional VAMP
algorithm, and then used the extended version to find a local optimum, for the
phase-shifts matrix under both ideal and practical constraints. The optimal
precoder at the BS, however, was found in closed-form using Lagrange
optimization. Simulation results illustrate the superiority of the proposed
approach over existing beamforming schemes (e.g., the SDR-and ADMM-based
approaches) both in terms of throughput and convergence speed. The results
also illustrate that the reduction in the throughput of the system under more
restrictive phase shifts is not significant. Moreover, it has been shown that
the performance of the proposed approach is largely unaffected by small
channel estimation errors. Its optimality can also be investigated through
state evolution analysis which is left for a follow-up work. Since the
proposed solution provides flexibility in terms of choosing the constraint on
the IRS reflection coefficients, it opens up the possibility of solving the
joint optimization problem using more physically consistent models for the IRS
reflection elements.

## Appendix A

We solve the following optimization problem:

$\displaystyle\operatorname*{arg\,min}_{\chi}\ f(\chi),$ (67)

where

$f(\chi)~{}\triangleq~{}\left\lvert\widetilde{r}+\frac{1}{1+{\rm
j}\chi}\right\rvert^{2},$ (68)

in which $\chi\in\operatorname{\mathbb{R}}$ and
$\widetilde{r}\in\operatorname{\mathbb{C}}$. Expanding the objective function,
we re-express it as follows:

$\displaystyle\operatorname*{arg\,min}_{\chi}\
\widetilde{r}^{*}\widetilde{r}+\frac{\widetilde{r}^{*}}{1+{\rm
j}\chi}+\frac{\widetilde{r}}{1-{\rm j}\chi}+\frac{1}{(1-{\rm j}\chi)(1+{\rm
j}\chi)}.$ (69)

Let $a\triangleq\Re\left\\{\widetilde{r}\right\\}$ and
$b\triangleq\Im\left\\{\widetilde{r}\right\\}$. We substitute $a$ and $b$ into
(69) and simplify the objective function as follows:

$\displaystyle\operatorname*{arg\,min}_{\chi}\
a^{2}+b^{2}+\frac{1+2a}{1+\chi^{2}}-\frac{2b\chi}{1+\chi^{2}}.$ (70)

By defining $c\triangleq(1+2a)$, we take the derivative w.r.t. $\chi$ and set
it to zero to obtain:

$f^{\prime}(\chi)~{}=~{}-\frac{2b(1-\chi^{2})}{(1+\chi^{2})^{2}}-\frac{2c\chi}{(1+\chi^{2})^{2}}~{}=~{}0.$
(71)

Simplifying (71) leads to:

$b\chi^{2}-c\chi-b~{}=~{}0.$ (72)

The roots of the quadratic equation in (72) are real and distinct and are
given by:

$\chi_{1}~{}=~{}\frac{c+\sqrt{c^{2}+4b^{2}}}{2b},$ (73)

and

$\chi_{2}=\frac{c-\sqrt{c^{2}+4b^{2}}}{2b}.$ (74)

where $b\neq 0$. By taking the second derivative of the objective function in
(70) w.r.t. $\chi$ and resorting to some straightforward algebraic
manipulations, we also obtain:

$f^{\prime\prime}(\chi)~{}=~{}\frac{2}{(1+\chi^{2})^{3}}(6b\chi-2b\chi^{3}+3c\chi^{2}-c).$
(75)

Substituting $\chi=\chi_{1}$ in (75) and simplifying the result yields:

$f^{\prime\prime}(\chi_{1})~{}=~{}\frac{1}{(1+\chi_{1}^{2})^{3}}\left(\frac{1}{b^{2}}\left(c^{3}+c^{2}\sqrt{c^{2}+4b^{2}}\right)+4\left(c+\sqrt{c^{2}+4b^{2}}\right)\right).$
(76)

Since $b\neq 0$, we have $c^{2}\sqrt{c^{2}+4b^{2}}>\left\lvert
c^{3}\right\rvert$ and $\sqrt{c^{2}+4b^{2}}>\left\lvert c\right\rvert$ which
implies that $f^{\prime\prime}(\chi_{1})>0$. Similarly, we have:

$f^{\prime\prime}(\chi_{2})~{}=~{}\frac{1}{(1+\chi_{2}^{2})^{3}}\left(\frac{1}{b^{2}}\left(c^{3}-c^{2}\sqrt{c^{2}+4b^{2}}\right)+4\left(c-\sqrt{c^{2}+4b^{2}}\right)\right)<0,\quad
b\neq 0.$ (77)

Thus, we choose:

$\chi^{\textnormal{opt}}~{}=~{}\chi_{1}~{}=~{}\frac{1+2a+\sqrt{(1+2a)^{2}+4b^{2}}}{2b},$
(78)

Interestingly, the solution $\chi_{1}$ results in the same sign for both
$\Im\left\\{-(1+{\rm j}\chi_{1})^{-1}\right\\}$ and
$\Im\left\\{\widetilde{r}\right\\}$.

## References

  * [1] F. Rusek, D. Persson, B. K. Lau, E. G. Larsson, T. L. Marzetta, O. Edfors, and F. Tufvesson, “Scaling Up MIMO: Opportunities and Challenges with Very Large Arrays,” _IEEE Signal Processing Magazine_ , vol. 30, no. 1, pp. 40–60, 2013.
  * [2] 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.
  * [3] T. L. Marzetta, “Noncooperative cellular wireless with unlimited numbers of base station antennas,” _IEEE Transactions on Wireless Communications_ , vol. 9, no. 11, pp. 3590–3600, 2010.
  * [4] S. Buzzi, C.-L. I, T. E. Klein, H. V. Poor, C. Yang, and A. Zappone, “A Survey of Energy-Efficient Techniques for 5G Networks and Challenges Ahead,” _IEEE Journal on Selected Areas in Communications_ , vol. 34, no. 4, pp. 697–709, 2016.
  * [5] S. Zhang, Q. Wu, S. Xu, and G. Y. Li, “Fundamental Green Tradeoffs: Progresses, Challenges, and Impacts on 5G Networks,” _IEEE Communications Surveys Tutorials_ , vol. 19, no. 1, pp. 33–56, 2017.
  * [6] X. Tan, Z. Sun, D. Koutsonikolas, and J. M. Jornet, “Enabling indoor mobile millimeter-wave networks based on smart reflect-arrays,” in _IEEE INFOCOM 2018 - IEEE Conference on Computer Communications_ , 2018, pp. 270–278.
  * [7] C. Liaskos, S. Nie, A. Tsioliaridou, A. Pitsillides, S. Ioannidis, and I. Akyildiz, “A new wireless communication paradigm through software-controlled metasurfaces,” _IEEE Communications Magazine_ , vol. 56, no. 9, pp. 162–169, 2018.
  * [8] 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.
  * [9] C. Pan, H. Ren, K. Wang, W. Xu, M. Elkashlan, A. Nallanathan, and L. Hanzo, “Multicell MIMO Communications Relying on Intelligent Reflecting Surfaces,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 8, pp. 5218–5233, 2020.
  * [10] 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.
  * [11] 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.
  * [12] Q. Wu and R. Zhang, “Beamforming Optimization for Wireless Network Aided by Intelligent Reflecting Surface With Discrete Phase Shifts,” _IEEE Transactions on Communications_ , vol. 68, no. 3, pp. 1838–1851, 2020.
  * [13] G. Zhou, C. Pan, H. Ren, K. Wang, and A. Nallanathan, “A Framework of Robust Transmission Design for IRS-Aided MISO Communications With Imperfect Cascaded Channels,” _IEEE Transactions on Signal Processing_ , vol. 68, pp. 5092–5106, 2020.
  * [14] H. Yang, Z. Xiong, J. Zhao, D. Niyato, L. Xiao, and Q. Wu, “Deep Reinforcement Learning-Based Intelligent Reflecting Surface for Secure Wireless Communications,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 1, pp. 375–388, 2021.
  * [15] X. Yu, D. Xu, and R. Schober, “Optimal Beamforming for MISO Communications via Intelligent Reflecting Surfaces,” in _2020 IEEE 21st International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)_ , 2020, pp. 1–5.
  * [16] P. Wang, J. Fang, X. Yuan, Z. Chen, and H. Li, “Intelligent Reflecting Surface-Assisted Millimeter Wave Communications: Joint Active and Passive Precoding Design,” _IEEE Transactions on Vehicular Technology_ , vol. 69, no. 12, pp. 14 960–14 973, 2020.
  * [17] S. Abeywickrama, R. Zhang, Q. Wu, and C. Yuen, “Intelligent Reflecting Surface: Practical Phase Shift Model and Beamforming Optimization,” _IEEE Transactions on Communications_ , vol. 68, no. 9, pp. 5849–5863, 2020\.
  * [18] B. Zheng, C. You, and R. Zhang, “Intelligent Reflecting Surface Assisted Multi-User OFDMA: Channel Estimation and Training Design,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 12, pp. 8315–8329, 2020\.
  * [19] C. Liu, X. Liu, D. W. K. Ng, and J. Yuan, “Deep residual learning for channel estimation in intelligent reflecting surface-assisted multi-user communications,” 2021.
  * [20] Y. Li, M. Jiang, Q. Zhang, and J. Qin, “Joint beamforming design in multi-cluster MISO NOMA reconfigurable intelligent surface-aided downlink communication networks,” _IEEE Transactions on Communications_ , vol. 69, no. 1, pp. 664–674, 2021.
  * [21] D. P. Bertsekas, “Nonlinear programming,” _Journal of the Operational Research Society_ , vol. 48, no. 3, pp. 334–334, 1997.
  * [22] 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.
  * [23] 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.
  * [24] 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.
  * [25] 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.
  * [26] R. W. Heath Jr. and A. Lozano, _Foundations of MIMO Communication_. Cambridge University Press, 2018.
  * [27] “approximate message-passing for convex optimization with non-separable penalties.”
  * [28] J. Barbier, F. Krzakala, N. Macris, L. Miolane, and L. Zdeborová, “Optimal errors and phase transitions in high-dimensional generalized linear models,” _Proceedings of the National Academy of Sciences_ , vol. 116, no. 12, pp. 5451–5460, 2019.
  * [29] D. L. Donoho, A. Maleki, and A. Montanari, “Message passing algorithms for compressed sensing: I. motivation and construction,” in _IEEE Information Theory Workshop on Information Theory (ITW 2010, Cairo)_ , 2010, pp. 1–5.
  * [30] C. M. Bishop, _Pattern Recognition and Machine Learning_. springer, 2006.
  * [31] M. Bayati and A. Montanari, “The dynamics of message passing on dense graphs, with applications to compressed sensing,” _IEEE Transactions on Information Theory_ , vol. 57, no. 2, pp. 764–785, 2011.
  * [32] Q. Li, Z. Zhu, and G. Tang, “Alternating Minimizations Converge to Second-Order Optimal Solutions,” in _Proceedings of the 36th International Conference on Machine Learning_ , ser. Proceedings of Machine Learning Research, K. Chaudhuri and R. Salakhutdinov, Eds., vol. 97. PMLR, 09–15 Jun 2019, pp. 3935–3943. [Online]. Available: http://proceedings.mlr.press/v97/li19n.html
  * [33] P. Jain, P. Netrapalli, and S. Sanghavi, “Low-Rank Matrix Completion Using Alternating Minimization,” in _Proceedings of the Forty-Fifth Annual ACM Symposium on Theory of Computing_. New York, NY, USA: Association for Computing Machinery, 2013, p. 665–674. [Online]. Available: https://doi.org/10.1145/2488608.2488693
  * [34] X. Yi, C. Caramanis, and S. Sanghavi, “Alternating Minimization for Mixed Linear Regression,” in _Proceedings of the 31st International Conference on Machine Learning_ , ser. Proceedings of Machine Learning Research, E. P. Xing and T. Jebara, Eds., vol. 32, no. 2. Bejing, China: PMLR, 22–24 Jun 2014, pp. 613–621. [Online]. Available: http://proceedings.mlr.press/v32/yia14.html
  * [35] P. Netrapalli, P. Jain, and S. Sanghavi, “Phase Retrieval Using Alternating Minimization,” _IEEE Transactions on Signal Processing_ , vol. 63, no. 18, pp. 4814–4826, 2015.
  * [36] A. Sinha, P. Malo, and K. Deb, “A review on bilevel optimization: From classical to evolutionary approaches and applications,” _IEEE Transactions on Evolutionary Computation_ , vol. 22, no. 2, pp. 276–295, 2018\.
  * [37] P. Hannan, “The element-gain paradox for a phased-array antenna,” _IEEE Transactions on Antennas and Propagation_ , vol. 12, no. 4, pp. 423–433, 1964\.
  * [38] 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.
  * [39] 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.
  * [40] 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.
  * [41] 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.
  * [42] 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.
  * [43] J. Yeh, _Real Analysis_ , 3rd ed. WORLD SCIENTIFIC, 2014. [Online]. Available: https://www.worldscientific.com/doi/abs/10.1142/9037