# Efficient Channel Estimation for RIS-Aided MIMO Communications with Unitary
Approximate Message Passing

Yabo Guo, Peng Sun, Zhengdao Yuan, Chongwen Huang, Qinghua Guo, , Zhongyong
Wang, and Chau Yuen The work of Y. Guo, P. Sun, Z. Yuan and Z. Wang was
supported by National Natural Science Foundation of China (61901417) and
Science and Technology Research Project of Henan Province (212102210173,
212102210556, 202102210313, 212102210566). (Corresponding authors: Qinghua Guo
and Zhongyong Wang.)Y. Guo, P. Sun and Z. Wang are with the School of
Information Engineering, Zhengzhou University, Zhengzhou 450002, China
(e-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>zywangzzu@gmail.com).Z. Yuan
is with the Artificial Intelligence Technology Engineering Research Center,
Open University of Henan, and also with the School of Information Engineering,
Zhengzhou University, Zhengzhou 450002, China (e-mail:
yuan_zhengdao@163.com).C. Huang is with the College of Information Science and
Electronic Engineering, Zhejiang University, Hangzhou 310007, China, and
Zhejiang Provincial Key Lab of Information Processing, Communication and
Networking (IPCAN), Hangzhou 310007, China, and the International Joint
Innovation Center, Zhejiang University, Haining 314400, China (e-mail:
chongwenhuang@zju.edu.cn).Q. Guo is with the School of Electrical, Computer
and Telecommunications Engineering, University of Wollongong, Wollongong, NSW
2522, Australia (e-mail: qguo@uow.edu.au).C. Yuen is with the Engineering
Product Development (EPD) Pillar, Singapore University of Technology and
Design, Singapore 487372 (e-mail: yuenchau@sutd.edu.sg).

###### Abstract

Reconfigurable intelligent surface (RIS) is very promising for wireless
networks to achieve high energy efficiency, extended coverage, improved
capacity, massive connectivity, etc. To unleash the full potentials of RIS-
aided communications, acquiring accurate channel state information is crucial,
which however is very challenging. For RIS-aided multiple-input and multiple-
output (MIMO) communications, the existing channel estimation methods have
computational complexity growing rapidly with the number of RIS units $N$
(e.g., in the order of $N^{2}$ or $N^{3}$) and/or have special requirements on
the matrices involved (e.g., the matrices need to be sparse for algorithm
convergence to achieve satisfactory performance), which hinder their
applications. In this work, instead of using the conventional signal model in
the literature, we derive a new signal model obtained through proper
vectorization and reduction operations. Then, leveraging the unitary
approximate message passing (UAMP), we develop a more efficient channel
estimator that has complexity linear with $N$ and does not have special
requirements on the relevant matrices, thanks to the robustness of UAMP. These
facilitate the applications of the proposed algorithm to a general RIS-aided
MIMO system with a larger $N$. Moreover, extensive numerical results show that
the proposed estimator delivers much better performance and/or requires
significantly less number of training symbols, thereby leading to notable
reductions in both training overhead and latency.

###### Index Terms:

Reconfigurable intelligent surface (RIS), channel estimation, approximate
message passing (AMP).

##  I Introduction

Reconfigurable intelligent surface (RIS) is a reconfigurable planar array with
a massive number of passive reflecting elements, which is capable of altering
the wireless propagation environment to achieve desired channel responses. RIS
has been recognized as a promising technology in future wireless
communications. With the aid of RIS, wireless networks are able to achieve
high energy efficiency, improve the system capacity and radio coverage,
enhance massive connectivity, etc [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,
14, 15, 16].

The research on RIS has attracted tremendous attention and many works have
been conducted to explore the potentials of RIS-aided communications. Under
the assumption of perfect channel state information (CSI), energy-efficient
designs were studied in [1] with the assist of RIS for wireless
communications. The works in [2] and [8] studied RIS-aided secrecy
communications. The work in [15] investigated the use of RIS to enhance the
energy harvesting and information transmission capabilities of the internet of
things systems. Energy efficient transmission with distributed RISs was
studied in [6]. The works in [14] and [17] studied RIS-aided non-othogonal
multiple access. RIS-aided communications were extended to millimeter-wave and
Terahertz communications in [18, 16, 9]. The works in [7] and [19] maximize
the weighted sum rate for both single and multi-cells networks. To unleash the
potentials of RIS aided communications, efficient accurate CSI acquisition is
crucial [5, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30].

In this work, we focus on channel estimation of RIS-aided multiple-input
multiple-output (MIMO) communications, which is challenging especially when
the number of RIS unit $N$ is large. In [25], a two-stage algorithm that
includes a sparse matrix factorization stage and a matrix completion stage was
proposed for the estimation of cascaded (transmitter-RIS and RIS-receiver)
channel, where the RIS phase matrix needs to be sparse and the channel between
RIS and BS is required to be a low-rank matrix. The work in [26] proposed a
trilinear semi-blind cascaded channel estimation problem, in which the
receiver estimates the channel coefficients and the transmitted signals
jointly, and a message passing algorithm was developed. The algorithm has
relatively low complexity, but it also requires that the relevant matrices are
sparse for the convergence of the algorithm, thereby achieving good
performance. A message-passing algorithm was proposed in [27] to estimate the
cascaded channels by exploiting the information on the slow-varying channel
components and the hidden channel sparsity, and its complexity growing with
$N^{2}$. The sparseness of the matrices may be achieved by using a sparse RIS
phase matrix, which however is not optimal in terms of channel estimation and
transmission, or by restricting the applications to special scenarios with
sparse channel matrices. The work in [28] proposed a three-phase pilot-based
channel estimation framework for RIS-assisted uplink multiuser communications,
and the complexity of the method is cubic in the number of RIS units $N$. In
[29], leveraging the parallel factor (PARAFAC) decomposition to represent high
dimensional tensors that involve unknown channels in different unfolded forms,
channel estimation methods based on alternating least squares (ALS) and vector
approximate message passing (VAMP) were proposed to estimate the channels from
BS to RIS and from RIS to users alternatively. Although there are no
requirements on the sparseness of the involved matrices, its complexity grows
with $N^{3}$, which is a concern for a large or even a moderate $N$. In
summary, the existing methods have special requirements on the involved
matrices and/or have the scalability issue as their computational complexity
rapidly grows with $N$, which hinder their applications.

To overcome the problems of the existing methods, leveraging the unitary
approximate message passing (UAMP) algorithm [31, 32, 33], we design a new RIS
channel estimator in this work. As a variant of the AMP algorithm [34], UAMP
achieves remarkably improved robustness by using a unitary transformation,
which enables it to deal with a linear reverse problem with a general (or
tough) system transfer matrix while with low complexity [32], [33]. In this
work, instead of using the conventional signal model for RIS channel
estimation in the literature, we derive a new signal model through proper
vectorization and reduction operations, which reformulates the channel
estimation to a structured signal recovery problem. Then a factor graph
representation is developed, and a message passing algorithm is derived, where
UAMP plays a crucial role. Thanks to the low complexity and robustness of
UAMP, the proposed algorithm is very efficient, which has complexity linear
with the number of RIS units $N$, and does not have special requirements on
the relevant matrices. These enable the applications of the proposed algorithm
to a general RIS aided MIMO system with a larger $N$. The Cramér-Rao lower
bound (CRLB) for the channel estimation considered is derived to serve as a
performance benchmark. Extensive numerical results show that, with much lower
complexity, the proposed channel estimator significantly outperforms existing
ones in terms of channel estimation performance and/or the training overhead.

The remainder of the paper is organized as follows. In Section II, we
introduce the system model and problem formulation of RIS channel estimation.
In Section III, the problem is reformulated, a new signal model is obtained
and the problem is represented as a factor graph. Then the UAMP based message
passing algorithm is developed in Section III-B, and the CRLB for channel
estimation is derived in Section V. Numerical results are provided in Section
VI, followed by conclusions in Section VII.

Notations-Boldface lower-case and upper-case letters denote vectors and
matrices, respectively. Superscripts $\bm{A}^{H}$ and $\bm{A}^{T}$ represent
conjugate transpose and transpose, respectively, and $\bm{A}^{*}$ represents
the conjugate of $\bm{A}$. A Gaussian distribution of $x$ with mean $\hat{x}$
and variance $\nu_{x}$ is denoted by
$\mathcal{N}\left(x;\hat{x},\nu_{x}\right).$ Notations $\otimes$ and $\odot$
represent the Kronecker and Khatri-Rao products, respectively. The relation
$f(x)=cg(x)$ for some positive constant $c$ is written as $f(x)\propto g(x)$.
We use $\bm{a}\cdot\bm{b}$ and $\bm{a}\cdot/\bm{b}$ to represent the element-
wise product and division between vectors $\bm{a}$ and $\bm{b}$, respectively.
The notation $\bm{a}^{.-1}$ denotes the element-wise inverse operation to
vector $\bm{a}$. We use $|\bm{A}|^{2}$ to denote element-wise magnitude
squared operation for $\bm{A}$, $\|\bm{a}\|$ to denote the $l_{2}$ norm of
$\bm{a}$, and $\|\bm{A}\|_{F}$ to denote the Frobenius norm of $\bm{A}$. The
notation $<\bm{a}>$ denotes the average operation for $\bm{a}$, i.e., the sum
of the elements of $\bm{a}$ divided by its length. We use $\mathbf{1}$ and
$\mathbf{0}$ to denote an all-one vector and an all-zero vector with a proper
length, respectively. The notation $\textbf{Diag}(\bm{a})$ represents a
diagonal matrix with $\bm{a}$ as its diagonal and $\bm{I}_{n}$ donates a
$n\times n$ identity matrix.

##  II System Model and Problem Formulation

Figure 1: Illustration of RIS-aided MIMO uplink transmission.

We consider RIS-aided MIMO uplink transmission as shown in Fig. 1, where the
BS equipped with $M$ antennas receives signals from $K$ users. A RIS with $N$
passive reflecting elements is equipped between the BS and users, each user
equipped with a single antenna (the extension to multiple antennas is
straightforward). The RIS is attached to the facade of a building in the
vicinity of users. Due to the highly attenuation caused by unfavorable
propagation environments such as tall buildings, the direct propagation path
between the BS and users is neglected [29].

We denote the channel matrix between the BS and the RIS by
$\bm{G}\in\mathbb{C}^{M\times N}$, and use
$\bm{H}\triangleq[\bm{h}_{1}^{u},\ldots,\bm{h}_{K}^{u}]\in\mathbb{C}^{N\times
K}$ to represent the channel matrix between the RIS and the $K$ users, where
$\bm{h}_{k}^{u}$ denote the channel vector from the $k$-th user to the RIS.
With the $l$-th ($l=1,2,\ldots,L$) RIS phase configuration, the $K$
consecutive received signals $\bm{Y}_{l}\in\mathbb{C}^{M\times K}$ is given by

$\displaystyle\bm{Y}_{l}=\bm{G}\textbf{Diag}(\bm{\Phi}_{l,:})\bm{H}\bm{X}+\bm{W}_{l},$
(1)

where $\bm{\Phi}_{l,:}$ is the $l$-th row of the RIS phase matrix
$\bm{\Phi}\in\mathbb{C}^{L\times N}$ (part of the RIS phase configurations)
and $\bm{W}_{l}$ models the zero mean complex additive white Gaussian noise
(AWGN) with precision $\beta$ (i.e., variance $\beta^{-1}$), and
$\bm{X}\in\mathbb{C}^{K\times K}$ denotes the transmitted orthogonal training
matrix from the users, i.e., $\bm{X}\bm{X}^{H}=\bm{I}_{K}$. Right-multiplying
both sides of $\eqref{eq:recvYl}$ by $\bm{X}^{H}$ leads to

$\displaystyle\bm{\tilde{Y}}_{l}=\bm{G}\textbf{Diag}(\bm{\Phi}_{l,:})\bm{H}+\bm{\tilde{W}}_{l},$
(2)

where $\bm{\tilde{Y}}_{l}\triangleq\bm{Y}_{l}\bm{X}^{H}\in\mathbb{C}^{M\times
K}$, and
$\bm{\tilde{W}}_{l}\triangleq\bm{W}_{l}\bm{X}^{H}\in\mathbb{C}^{M\times K}$.
As $\bm{X}$ is a unitary matrix, entries of $\bm{\tilde{W}}_{l}$ are still
zero-mean white Gaussian with the same noise precision $\beta$.

During the training process, $L$ configurations of the RIS are used, leading
to $L$ matrices $\\{\bm{\tilde{Y}}_{l},l=1,...,L\\}$, based which we aim to
estimate the channel matrices. As the number of RIS unit $N$ can be large, low
complexity of channel estimators is crucial for practical implementation. The
cubic or quadratic complexity (in terms of $N$) of existing channel estimation
algorithms can be a concern in their applications. During the training
process, the RIS undergoes $L$ configurations. In order to reduce the training
overhead and communication latency, a small $L$ is desirable. We address the
above challenges by reformulating the RIS channel estimation to a structured
signal recovery problem with a new signal model and develop an efficient
message passing algorithm. In particular, we do not make any special
requirements on the matrices $\bm{G}$, $\bm{H}$ or $\bm{\Phi}$.

##  III Problem Reformulation and Factor Graph Representation

###  III-A New Model for RIS Channel Estimation

Instead of using model $\eqref{eq:recvYl_nopilot}$ directly, we reformulate it
to a new model through vectorization and reduction, which leads to a
structured signal recovery problem. Vectorizing $\eqref{eq:recvYl_nopilot}$
leads to

$\displaystyle\textbf{vec}(\bm{\tilde{Y}}_{l})=(\bm{H}^{T}\otimes\bm{G})\textbf{vec}\left(\textbf{Diag}(\bm{\Phi}_{l,:})\right)+\textbf{vec}(\bm{\tilde{W}}_{l}),$
(3)

where the vector $\textbf{vec}\left(\textbf{Diag}(\bm{\Phi}_{l,:})\right)$ can
be represented as

$\textbf{vec}\left(\textbf{Diag}(\bm{\Phi}_{l,:})\right)=[{\phi}_{l,1},\mathbf{0}_{N}^{T},{\phi}_{l,2},\mathbf{0}_{N}^{T},\ldots,\mathbf{0}_{N}^{T},{\phi}_{l,N}]^{T},$
(4)

and $\phi_{l,n}$ donates the $n$-th ($n=1,2,\ldots,N$) element of
$\bm{\Phi}_{l,:}$. Here we note that, the vector
$\textbf{vec}\left(\textbf{Diag}(\bm{\Phi}_{l,:})\right)$ is highly sparse,
and the non-zeros elements of
$\textbf{vec}\left(\textbf{Diag}(\bm{\Phi}_{l,:})\right)$ are separated by
all-zero vector $\mathbf{0}_{N}$. This can be exploited to significantly
reduce the dimension of $\eqref{eq:vecYl_nopilot}$. We drop the zero elements
in $\textbf{vec}\left(\textbf{Diag}(\bm{\Phi}_{l,:})\right)$ and the
corresponding columns in $\bm{H}^{T}\otimes\bm{G}$, then (3) can be reduced to

$\bm{\tilde{y}}_{l}=\bm{\tilde{S}}\bm{\Phi}^{T}_{l,:}+\bm{\tilde{w}}_{l},$ (5)

where
$\bm{\tilde{y}}_{l}\triangleq\textbf{vec}(\bm{\tilde{Y}}_{l})\in\mathbb{C}^{KM\times
1}$, $\bm{\tilde{S}}\triangleq\bm{H}^{T}\odot\bm{G}\in\mathbb{C}^{KM\times
N}$,
$\bm{\tilde{S}}\bm{\Phi}^{T}_{l,:}=(\bm{H}^{T}\otimes\bm{G})\textbf{vec}\left(\textbf{Diag}(\bm{\Phi}_{l,:})\right)$,
and
$\bm{\tilde{w}}_{l}\triangleq\textbf{vec}(\bm{\tilde{W}}_{l})\in\mathbb{C}^{KM\times
1}$. By stacking $\\{\bm{\tilde{y}}_{l},l=1,...,L\\}$ into a matrix, we have
$\bm{\tilde{Y}}=\bm{\tilde{S}}\bm{\Phi}^{T}+\bm{\tilde{W}}$, which is
rewritten as

$\bm{Y}=\bm{\Phi}\bm{S}+\bm{W}.$ (6)

where $\bm{Y}=\bm{\tilde{Y}}^{T}$, $\bm{W}=\bm{\tilde{W}}^{T}$and

$\bm{S}=\bm{\tilde{S}}^{T}=(\bm{H}^{T}\odot\bm{G})^{T}\in\mathbb{C}^{N\times
KM}.$ (7)

Now we can see from (6) that, the channel estimation is reformulated as the
recovery of the matrix $\bm{S}$, which admits the structure (7). We can also
treat $\bm{S}$ as an intermediate variable, as our aim is to estimate $\bm{H}$
and $\bm{G}$.

###  III-B Probabilistic and Factor Graph Representation

We consider recovering the signal using the message passing techniques, in
particular leveraging UAMP to achieve low complexity while with high
robustness. We first represent the problem in a probabilistic form. It is
noted that UAMP works with a unitary transform of the linear observation model
(6), which is crucial to achieving high robustness. So, to facilitate the use
of UAMP later, we first carry out a unitary transformation to (6) based on the
singular value decomposition (SVD) $\bm{\Phi}=\bm{U\Lambda V}$, leading to

$\bm{R}=\bm{\Psi}\bm{S}+\bm{\overline{W}},$ (8)

where $\bm{R}=\bm{U}^{H}\bm{Y}$, $\bm{\Psi}=\bm{U}^{H}\bm{\Phi}=\bm{\Lambda
V}$ and $\bm{\overline{W}}=\bm{U}^{H}\bm{W}$. Since $\bm{U}$ is unitary,
entries in $\bm{\overline{W}}$ are still AWGN with precision $\beta$.

Note that $\bm{\tilde{S}}=[\bm{\tilde{s}}_{1},\ldots,\bm{\tilde{s}}_{N}]$,
$\bm{H}^{T}=[\bm{h}_{1},\ldots,\bm{h}_{N}]$ and
$\bm{G}=[\bm{g}_{1},\ldots,\bm{g}_{N}]$. Then, according to
$\bm{\tilde{S}}=\bm{H}^{T}\odot\bm{G}$, we have

$\bm{\tilde{s}}_{n}=\bm{h}_{n}\otimes\bm{g}_{n},$ (9)

where $\bm{h}_{n}=[{h}_{1,n},\ldots,{h}_{K,n}]^{T}$,
$\bm{g}_{n}=[{g}_{1,n},\ldots,{g}_{M,n}]^{T}$ and
$\bm{\tilde{s}}_{n}=[{\tilde{s}}_{1,1,n},\ldots,{\tilde{s}}_{m,k,n},\ldots,{\tilde{s}}_{M,K,n}]^{T}$
with ${\tilde{s}}_{m,k,n}={h}_{k,n}{g}_{m,n}$. Let $J=KM$, and note that
$\bm{R}=[\bm{r}_{1},\ldots,\bm{r}_{J}]$,
$\bm{S}=[\bm{s}_{1},\ldots,\bm{s}_{J}]$, and
$\bm{\overline{W}}=[\bm{w}_{1},\ldots,\bm{w}_{J}]$. Define an auxiliary
variable $\bm{Z}\triangleq[\bm{z}_{1},\ldots,\bm{z}_{J}]$ with
$\bm{z}_{j}=\bm{\Psi}\bm{s}_{j}$. Then, the joint distribution of
$\bm{H},\bm{G},\bm{S},\bm{\tilde{S}},\bm{Z}$ and $\beta$ given $\bm{R}$ can be
factorized as

$\displaystyle p$
$\displaystyle(\bm{H},\bm{G},\bm{S},\bm{\tilde{S}},\bm{Z},\beta|\bm{R})$
$\displaystyle\propto$
$\displaystyle~{}p(\bm{R}|\bm{Z},\beta)p(\bm{\tilde{S}}|\bm{S})p(\bm{Z}|\bm{S})p(\bm{\tilde{S}}|\bm{H},\bm{G})p(\bm{H})p(\bm{G})p(\beta)$
$\displaystyle=$
$\displaystyle~{}p(\bm{\tilde{S}}|\bm{S})p(\beta){\prod}_{n}p(\bm{\tilde{s}}_{n}|\bm{h}_{n},\bm{g}_{n})p(\bm{h}_{n})p(\bm{g}_{n}){\prod}_{j}p(\bm{r}_{j}|\bm{z}_{j},\beta)$
$\displaystyle\qquad\qquad\times p(\bm{z}_{j}|\bm{s}_{j})$
$\displaystyle\triangleq$
$\displaystyle~{}f_{\bm{S}}(\bm{\tilde{S}},\bm{S})f_{\beta}(\beta){\prod}_{n}f_{\bm{\tilde{s}}_{n}}(\bm{\tilde{s}}_{n},\bm{h}_{n},\bm{g}_{n})f_{\bm{h}_{n}}(\bm{h}_{n})f_{\bm{g}_{n}}(\bm{g}_{n})$
$\displaystyle\qquad\qquad\times{\prod}_{j}f_{\bm{r}_{j}}(\bm{z}_{j},\beta)f_{\bm{z}_{j}}(\bm{z}_{j},\bm{s}_{j}).$
(10)

where the involved distributions are listed in Table I. To facilitate the
factor graph representation of the factorization in (10), local functions
(factors) are defined, and the correspondence between the distributions and
local functions are also shown in Table I. It is noted that $p(\bm{H})$ and
$p(\bm{G})$ represent the priors for the channel matrices $\bm{H}$ and
$\bm{G}$, respectively. When no priors are available for the channel matrices,
the priors can be set to be a non-informative one, e.g.,
$\rho_{h}=\rho_{g}=+\infty$ in Table I. The factor graph representation is
depicted in Fig. 2.

TABLE I: Factors and distributions in (10). Factor |  Distribution |  Function
---|---|---
$f_{\beta}$ |  $p(\beta)$ |  $\propto\beta^{-1}$
$f_{\bm{h}_{n}}$ |  $p(\bm{h}_{n})$ |  $\mathcal{CN}(\bm{h}_{n};\bm{0}_{K},\rho_{h}\bm{I}_{K})$
$f_{\bm{g}_{n}}$ |  $p(\bm{g}_{n})$ |  $\mathcal{CN}(\bm{g}_{n};\bm{0}_{M},\rho_{g}\bm{I}_{M})$
$f_{\bm{r}_{j}}$ |  $p(\bm{r}_{j}|\bm{z}_{j},\beta)$ |  $\mathcal{N}(\bm{r}_{j};\bm{z}_{j},\beta^{-1}\bm{I}_{L})$
$f_{\bm{z}_{j}}$ |  $p(\bm{z}_{j}|\bm{s}_{j})$ |  $\delta(\bm{z}_{j}-\bm{\Psi}\bm{s}_{j})$
$f_{\bm{\tilde{s}}_{n}}$ |  $p(\bm{\tilde{s}}_{n}|\bm{h}_{n},\bm{g}_{n})$ |  $\delta(\bm{\tilde{s}}_{n}-\bm{h}_{n}\otimes\bm{g}_{n})$
$f_{\bm{S}}$ |  $p(\bm{\tilde{S}}|\bm{S})$ |  $\delta(\bm{\tilde{S}}^{T}-\bm{S})$
$\bm{z}_{1}$$\bm{\vdots}$$\bm{z}_{J}$$f_{\bm{z}_{1}}$$f_{\bm{z}_{J}}$$\bm{s}_{1}$$\bm{s}_{J}$$f_{\bm{\tilde{s}}_{1}}$$f_{\bm{\tilde{s}}_{N}}$$\bm{\tilde{s}}_{1}$$\bm{\vdots}$$\bm{\tilde{s}}_{N}$$\bm{g}_{1}$$\bm{g}_{N}$$f_{\bm{r}_{1}}$$f_{\bm{r}_{J}}$$\beta$$f_{\beta}$$f_{\bm{S}}$$\bm{h}_{1}$$\bm{\vdots}$$\bm{h}_{N}$$f_{\bm{h}_{1}}$$f_{\bm{g}_{1}}$$f_{\bm{h}_{N}}$$f_{\bm{g}_{N}}$
Figure 2:  Factor graph representation of (10).

Our aim is to find the posteriori distributions $p(\bm{H}|\bm{R})$ and
$p(\bm{G}|\bm{R})$ and their estimates in terms of a posteriori means, i.e.,
$\hat{\bm{H}}=\mathbb{E}\left\\{\bm{H}|\bm{R}\right\\}$ and
$\hat{\bm{G}}=\mathbb{E}\left\\{\bm{G}|\bm{R}\right\\}$. We note that it is
difficult to find the exact a posteriori distributions, and approximate
inference has to be resorted. In next section, we will develop a low
complexity message passing algorithm to find their approximations. To
facilitate the message passing algorithm design, a scalar factor graph
representation is shown in Fig. 3.

Figure 3: Scalar factor graph representation of (10).

##  IV UAMP Based Message Passing Algorithm For RIS Channel Estimation

In this section, we develop an efficient message passing algorithm for
approximate inference, where UAMP is incorporated to deal with the most
computationally intensive part of message passing, which is crucial to
achieving low complexity while with high robustness. The message passing
algorithm carries out an iterative process, where each iteration involves a
forward message passing process and backward message passing process in the
graph shown in Fig. 2 or Fig. 3. We use $m_{A\rightarrow B}(x)$ to denote a
message passed from node $A$ to node $B$, which is a function of $x$. For
Gaussian messages, the arrows above its mean and variance indicate the
direction of the message passing. In addition, we use $b(x)$ to denote the
belief of a variable $x$. Note that, if a forward computation requires
backward messages, the relevant messages in the previous iteration is used by
default.

###  IV-A Forward Message Passing

In the forward direction, according to the rules of variational message
passing [35], we have

$\displaystyle m_{f_{\bm{r}_{j}}\rightarrow\bm{z}_{j}}(\bm{z}_{j})$
$\displaystyle\propto\exp\left\\{\int_{\beta}\mathfrak{b}(\beta)\log
f_{\bm{r}_{j}}\right\\}$
$\displaystyle\propto\mathcal{N}\left(\bm{z}_{j};\bm{r}_{j},\hat{\beta}^{-1}\right),$
(11)

where

$\displaystyle\mathfrak{b}(\beta)$ $\displaystyle\propto
m_{f_{\bm{r}_{j}}\rightarrow\beta}(\beta)f_{\beta}$
$\displaystyle\propto\beta^{LJ-1}\exp\left\\{\sum_{j=1}-\beta\left(\|\bm{r}_{j}-\hat{\bm{z}_{j}}\|^{2}+\mathbf{1}^{T}\bm{\nu}_{\bm{z}_{j}}\right)\right\\},$
(12)

and

$\hat{\beta}=\int_{\beta}\beta\mathfrak{b}(\beta)=\frac{LJ}{\sum_{j=1}(\|\bm{r}_{j}-\hat{\bm{z}_{j}}\|^{2}+\mathbf{1}_{L}^{T}\bm{\nu}_{\bm{z}_{j}})}.$
(13)

It is noted that in the above equation, $\hat{\bm{z}}_{j}$ and the message
$m_{f_{\bm{r}_{j}}\rightarrow\beta}(\beta)$ are required to compute
$\hat{\beta}$, which are obtained from the last iteration and their
computations are delayed to (71) and (74).

The Gaussian form of the message in (IV-A) suggests the following model

$\bm{r}_{j}=\bm{z}_{j}+\bm{w}_{j},$ (14)

where the noise $\bm{w}_{j}$ is Gaussian with mean zero and precision
$\hat{\beta}$. This allows seamless integration with the forward recursion of
UAMP. Here, assume that the mean and variance of $\bm{s}_{j}$ are available,
i.e., $\hat{\bm{s}}_{j}$ and $\nu_{\bm{s}_{j}}$. Specially, we assume
$\bm{s}_{j}$ have a common variance $\nu_{\bm{s}_{j}}$, and the corresponding
computation will be detailed later. Following UAMP, we define a vector
$\bm{\psi}\in\mathbb{C}^{L\times 1}$ as

$\bm{\psi}=|\bm{\Psi}|^{2}\bm{1}_{N}.$ (15)

Then we calculate two vectors $\bm{\nu_{p_{j}}}$ and $\bm{p}_{j}$ as

$\bm{\nu_{p_{j}}}=\bm{\psi}\nu_{\bm{{s}}_{j}},$ (16)
$\bm{p}_{j}=\bm{\Psi}\hat{{\bm{s}}}_{j}-\bm{\nu_{p_{j}}}\cdot\bm{\mu}_{j},$
(17)

where $\bm{\mu}_{j}$ is a vector that is computed in last iteration. According
to UAMP, we update the intermediate vectors $\bm{\nu}_{\bm{\mu}_{j}}$ and
$\bm{\mu}_{j}$ by

$\displaystyle\bm{\nu}_{\bm{\mu}_{j}}$
$\displaystyle=\mathbf{1}_{L}./(\bm{\nu}_{\bm{p}_{j}}+\hat{\beta}^{-1}\mathbf{1}_{L}),$
(18) $\displaystyle\bm{\mu}_{j}$
$\displaystyle=\bm{\nu}_{\bm{\mu}_{j}}\cdot(\bm{r}_{j}-\bm{p}_{j}),$ (19)

Then compute vectors $\bm{\nu}_{\bm{q}_{j}}$ and $\bm{q}_{j}$ with

$\displaystyle\bm{\nu}_{\bm{q}_{j}}$
$\displaystyle=\mathbf{1}_{N}./|\bm{\Psi}^{H}|^{2}\bm{\nu}_{\bm{\mu}_{j}},$
(20) $\displaystyle\bm{q}_{j}$
$\displaystyle=\hat{\bm{s}}_{j}+\bm{\nu}_{\bm{q}_{j}}.\bm{\Psi}^{H}\bm{\mu}_{j},$
(21)

The message $\bm{q}_{j}$ and $\bm{\nu}_{\bm{q}_{j}}$ are the mean and variance
of $\bm{s}_{j}$. According to the belief propagation derivation of (U)AMP, we
have

$m_{\bm{s}_{j}\rightarrow
f_{\bm{s}}}(\bm{s}_{j})=\mathcal{N}(\bm{s}_{j};\bm{q}_{j},\textbf{Diag}(\bm{\nu}_{\bm{q}_{j}})).$
(22)

Stack $\bm{q}_{j}$ and $\bm{\nu}_{\bm{q}_{j}}$ into matrices as

$\displaystyle\bm{Q}=\left[\bm{q}_{1},\ldots,\bm{q}_{J}\right],$ (23)
$\displaystyle\bm{\nu}_{\bm{Q}}=\left[\bm{\nu}_{{\bm{q}}_{1}},\ldots,\bm{\nu}_{{\bm{q}}_{J}}\right],$
(24)

Due to the deterministic relationship between $\bm{S}$ and $\bm{\tilde{S}}$,
i.e., $\bm{\tilde{S}}^{T}=\bm{S}$, as shown in Table I, we have

$\displaystyle\bm{\tilde{Q}}=\bm{Q}^{T}=\left[\bm{\tilde{q}}_{1},\ldots,\bm{\tilde{q}}_{N}\right],$
(25)
$\displaystyle\bm{\nu}_{\bm{\tilde{Q}}}=\bm{\nu}^{T}_{\bm{Q}}=[\bm{\nu}_{\bm{\tilde{q}}_{1}},\ldots,\bm{\nu}_{\bm{\tilde{q}}_{N}}],$
(26)

and the vectors $\bm{\tilde{q}}_{n}\in\mathbb{C}^{J\times 1}$ and
$\bm{\nu}_{\bm{\tilde{q}}_{n}}\in\mathbb{C}^{J\times 1}$ can be divided into
$K$ length-$M$ vectors, i.e.,

$\bm{\tilde{q}}_{n}=\left[\bm{\tilde{q}}^{T}_{1,n},\ldots,\bm{\tilde{q}}^{T}_{K,n}\right]^{T},$
(27)
$\bm{\nu}_{\bm{\tilde{q}}_{n}}=[\bm{\nu}^{T}_{\bm{\tilde{q}}_{1,n}},\ldots,\bm{\nu}^{T}_{\bm{\tilde{q}}_{K,n}}]^{T}.$
(28)

Let
${\nu}_{\bm{\tilde{q}}_{k,n}}=\langle{\bm{\nu}_{\bm{\tilde{q}}_{k,n}}}\rangle$,
and the message $m_{\bm{\tilde{s}}_{n}\rightarrow
f_{\bm{\tilde{s}}_{n}}}(\bm{\tilde{s}}_{n})$ can be expressed as

$m_{\bm{\tilde{s}}_{n}\rightarrow
f_{\bm{\tilde{s}}_{n}}}(\bm{\tilde{s}}_{n})=\mathcal{N}(\bm{\tilde{s}}_{n};\bm{\tilde{q}}_{n},\textbf{Diag}(\bm{{\nu}}^{\prime}_{\bm{\tilde{q}}_{n}})),$
(29)

where

$\bm{{\nu}}^{\prime}_{\bm{\tilde{q}}_{n}}=\left[{\nu}_{\bm{\tilde{q}}_{1,n}},\ldots,{\nu}_{\bm{\tilde{q}}_{K,n}}\right]^{T}\otimes\bm{1}_{M}.$
(30)

Next, we compute the message at $f_{\bm{\tilde{s}}_{n}}$, $\bm{h}_{n}$ and
$\bm{g}_{n}$. It is noted that $\bm{\tilde{s}}_{n}\in\mathbb{C}^{J\times 1}$
also can be divided into $K$ length-$M$ vectors and expressed as

$\bm{\tilde{s}}_{n}=[\bm{\tilde{s}}_{1,n}^{T},\ldots,\bm{\tilde{s}}_{K,n}^{T}]^{T},$
(31)

where
$\bm{\tilde{s}}_{k,n}\triangleq[{\tilde{s}}_{1,k,n},\ldots,{\tilde{s}}_{M,k,n},]^{T}$.
We further factorize the function
$f_{\bm{\tilde{s}}_{n}}(\bm{\tilde{s}}_{n},\bm{h}_{n},\bm{g}_{n})$ as

$f_{\bm{\tilde{s}}_{n}}(\bm{\tilde{s}}_{n},\bm{h}_{n},\bm{g}_{n})=\prod_{m,k}f_{{\tilde{s}}_{m,k,n}}({h}_{k,n},{g}_{m,n}),$
(32)

and the factor $f_{\tilde{s}_{m,k,n}}({h}_{k,n},{g}_{m,n})$ is shown in Fig.
4.

$\tilde{s}_{m^{\prime},k,n}$$\tilde{s}_{m,k,n}$$\tilde{s}_{m,k^{\prime},n}$$f_{\tilde{s}_{m^{\prime},k,n}}$$f_{\tilde{s}_{m,k,n}}$$f_{\tilde{s}_{m,k^{\prime},n}}$$\bm{\vdots}$$\bm{\vdots}$$h_{k,n}$$g_{m,n}$$f_{{h_{k,n}}}$$f_{\bm{g_{m,n}}}$
Figure 4:  Factor graph representation of $f_{{\tilde{s}}_{m,k,n}}$.

With the definition
$\bm{\tilde{q}}_{k,n}\triangleq[{\tilde{q}}_{1,k,n},\ldots,{\tilde{q}}_{M,k,n}]^{T}$,
(29) implies that

$\displaystyle m_{\tilde{s}_{m,k,n}\rightarrow
f_{\tilde{s}_{m,k,n}}}\left(\tilde{s}_{m,k,n}\right)=\mathcal{N}\left(\tilde{s}_{m,k,n};\tilde{q}_{m,k,n},{\nu}_{\bm{\tilde{q}}_{k,n}}\right),$
(33)

and we note that the factor
$f_{\tilde{s}_{m,k,n}}=\delta\left(\tilde{s}_{m,k,n}-h_{k,n}g_{m,n}\right)$ as
shown in Table I. To compute the message $m_{f_{\tilde{s}_{m,k,n}}\rightarrow
g_{m,n}}\left(g_{m,n}\right)$ with belief propagation at factor node
$f_{\tilde{s}_{m,k,n}}$, we need to integrate out $s_{m,k,n}$ and $h_{k,n}$.
However, due to the multiplication of $g_{m,n}$ and $h_{k,n}$, the message
will be intractable even the incoming message $m_{h_{k,n}\rightarrow
f_{\tilde{s}_{m,k,n}}}\left(h_{k,n}\right)$ is Gaussian. To solve this, we
first apply belief propagation and eliminate the variable $\tilde{s}_{m,k,n}$
to get an intermediate function node
$\tilde{f}_{\tilde{s}_{m,k,n}}(h_{k,n},g_{m,n})$, i.e.,

$\displaystyle\tilde{f}_{\tilde{s}_{m,k,n}}(h_{k,n},g_{m,n})$
$\displaystyle\\!=\\!\\!\int_{\tilde{s}_{m,k,n}}\\!\\!\\!\\!m_{\tilde{s}_{m,k,n}\rightarrow
f_{\tilde{s}_{m,k,n}}}\\!\\!\left(\tilde{s}_{m,k,n}\right)\cdot
f_{\tilde{s}_{m,k,n}}$
$\displaystyle=\mathcal{N}\left({\tilde{s}}_{m,k,n};\tilde{q}_{m,k,n},{\nu}_{\bm{\tilde{q}}_{k,n}}\right).$
(34)

This turns the function node $f_{\tilde{s}_{m,k,n}}$ with the hard constraint
$\delta\left(\tilde{s}_{m,k,n}-h_{k,n}g_{m,n}\right)$ to a soft function node,
enabling the use of variational inference to handle $h_{k,n}$ and $g_{m,n}$.
With the intermediate local function
$\tilde{f}_{\tilde{s}_{m,k,n}}(h_{k,n},g_{m,n})$, we can compute the outgoing
message $m_{f_{\tilde{s}_{m,k,n}}\rightarrow g_{m,n}}\left(g_{m,n}\right)$ as

$\displaystyle m_{f_{\tilde{s}_{m,k,n}}\rightarrow
g_{m,n}}\left(g_{m,n}\right)$
$\displaystyle\propto\exp\left\\{\int_{h_{k,n}}\mathfrak{b}\left(h_{k,n}\right)\log\tilde{f}_{\tilde{s}_{m,k,n}}\right\\}$
$\displaystyle\propto\mathcal{N}(g_{m,n};\overset{\scriptscriptstyle\rightarrow}{g}_{m,k,n},\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{m,k,n}}),$
(35)

where

$\displaystyle\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{m,k,n}}=\frac{{\nu}_{\bm{\tilde{q}}_{k,n}}}{\left|\hat{h}_{k,n}\right|^{2}+\nu_{h_{k,n}}},$
(36)
$\displaystyle\overset{\scriptscriptstyle\rightarrow}{g}_{m,k,n}=\frac{\tilde{q}_{m,k,n}\hat{h}_{k,n}^{*}}{\left|\hat{h}_{k,n}\right|^{2}+\nu_{h_{k,n}}},$
(37)

with $\hat{h}_{k,n}$ and $\nu_{h_{k,n}}$ being the approximate a posteriori
mean and variance of ${h}_{k,n}$, which are computed in (52) and (51). With
belief propagation and referring to Fig. 4, the message $m_{g_{m,n}\rightarrow
f_{{g_{m,n}}}}\left(g_{m,n}\right)$ can be represented as

$\displaystyle m_{g_{m,n}\rightarrow
f_{{g_{m,n}}}}\left(g_{m,n}\right)=\mathcal{N}\left(g_{m,n};\overset{\scriptscriptstyle\rightarrow}{g}_{m,n},\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{m,n}}\right),$
(38)

with

$\displaystyle\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{m,n}}=1/\sum_{k=1}^{K}\frac{1}{\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{m,k,n}}},$
(39)
$\displaystyle\overset{\scriptscriptstyle\rightarrow}{g}_{m,n}=\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{m,n}}\sum_{k=1}^{K}\frac{\overset{\scriptscriptstyle\rightarrow}{g}_{m,k,n}}{\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{m,k,n}}}.$
(40)

So, the marginal of $g_{m,n}$ can be expressed as

$\displaystyle\mathfrak{b}\left(g_{m,n}\right)$
$\displaystyle=m_{g_{m,n}\rightarrow
f_{{g}_{m,n}}}\left(g_{m,n}\right)f_{{g_{m,n}}}$
$\displaystyle\propto\mathcal{N}\left(g_{m,n};\hat{g}_{m,n},\nu_{g_{m,n}}\right),$
(41)

with

$\displaystyle\nu_{g_{m,n}}=\frac{\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{m,n}}\rho_{g}}{\rho_{g}+\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{m,n}}},$
(42)
$\displaystyle\hat{g}_{m,n}=\frac{\overset{\scriptscriptstyle\rightarrow}{g}_{m,n}\rho_{g}}{\rho_{g}+\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{m,n}}}.$
(43)

Similarly, we can compute the message $m_{f_{\tilde{s}_{m,k,n}}\rightarrow
h_{k,n}}$ as

$\displaystyle m_{f_{\tilde{s}_{m,k,n}}\rightarrow
h_{k,n}}\left(h_{k,n}\right)\propto\mathcal{N}\left(h_{k,n};\overset{\scriptscriptstyle\rightarrow}{h}_{m,k,n},\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{m,k,n}}\right),$
(44)

where

$\displaystyle\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{m,k,n}}=\frac{{\nu}_{\bm{\tilde{q}}_{k,n}}}{\left|\hat{g}_{m,n}\right|^{2}+\nu_{g_{m,n}}},$
(45)
$\displaystyle\overset{\scriptscriptstyle\rightarrow}{h}_{m,k,n}=\frac{\tilde{q}_{m,k,n}\hat{g}_{m,n}^{*}}{\left|\hat{g}_{m,n}\right|^{2}+\nu_{g_{m,n}}},$
(46)

with $\hat{g}_{m,n}$ and $\nu_{g_{m,n}}$ being the approximate a posteriori
mean and variance of ${g}_{m,n}$, which are computed in (43) and (42). With
belief propagation, the message $m_{h_{k,n}\rightarrow
f_{{h_{k,n}}}}\left(h_{k,n}\right)$ can be represented as

$\displaystyle m_{h_{k,n}\rightarrow
f_{{h_{k,n}}}}\left(h_{k,n}\right)=\mathcal{N}\left(h_{k,n};\overset{\scriptscriptstyle\rightarrow}{h}_{k,n},\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{k,n}}\right),$
(47)

with

$\displaystyle\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{k,n}}=1/\sum_{m=1}^{M}\frac{1}{\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{m,k,n}}},$
(48)
$\displaystyle\overset{\scriptscriptstyle\rightarrow}{h}_{k,n}=\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{k,n}}\sum_{m=1}^{M}\frac{\overset{\scriptscriptstyle\rightarrow}{h}_{m,k,n}}{\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{m,k,n}}}.$
(49)

The marginal of $h_{k,n}$ can be expressed as

$\displaystyle\mathfrak{b}\left(h_{k,n}\right)$
$\displaystyle=m_{h_{k,n}\rightarrow
f_{{h}_{k,n}}}\left(h_{k,n}\right)f_{{h_{k,n}}}$
$\displaystyle\propto\mathcal{N}\left(h_{k,n};\hat{h}_{k,n},\nu_{h_{k,n}}\right),$
(50)

with

$\displaystyle\nu_{h_{k,n}}=\frac{\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{k,n}}\rho_{h}}{\rho_{h}+\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{k,n}}},$
(51)
$\displaystyle\hat{h}_{k,n}=\frac{\overset{\scriptscriptstyle\rightarrow}{h}_{k,n}\rho_{h}}{\rho_{h}+\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{k,n}}}.$
(52)

This is the end of forward message passing.

###  IV-B Backward message passing

Next, we elaborate the backward message passing. The backward message from
$h_{k,n}$ to $f_{\tilde{s}_{m,k,n}}$ can be expressed as

$\displaystyle m_{h_{k,n}\rightarrow
f_{\tilde{s}_{m,k,n}}}\left(h_{k,n}\right)=\frac{\mathfrak{b}\left(h_{k,n}\right)}{m_{f_{\tilde{s}_{m,k,n}\rightarrow
h_{k,n}}}\left(h_{k,n}\right)},$ (53)

They are represented collectively as $m_{\bm{h}_{n}\rightarrow
f_{\bm{\tilde{s}}_{n}}}\left(\bm{h}_{n}\right)$, which is Gaussian with mean
$\overset{\scriptscriptstyle\leftarrow}{\bm{h}}_{n}$ and variance
$\textbf{Diag}(\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{h}_{n}})$.
With the factor graph shown in Fig. 4, the mean and variance can be computed
as

$\displaystyle\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{h}_{n}}=\left(\left(\mathbf{1}_{K}./\bm{\nu}_{\bm{h}_{n}}\right)\otimes\mathbf{1}_{M}-\mathbf{1}_{J}./\overset{\scriptscriptstyle\rightarrow}{\bm{\nu}}_{\bm{h}_{n}}\right)^{.-1},$
(54)
$\displaystyle\overset{\scriptscriptstyle\leftarrow}{\bm{h}}_{n}=\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{h}_{n}}\cdot\left((\hat{\bm{h}}_{n}./\bm{\nu}_{\bm{h}_{n}})\otimes\mathbf{1}_{M}-\overset{\scriptscriptstyle\rightarrow}{\bm{h}}_{n}./\overset{\scriptscriptstyle\rightarrow}{\bm{\nu}}_{\bm{h}_{n}}\right),$
(55)

where $\bm{\nu}_{\bm{h}_{n}}=[\nu_{h_{1,n}},\ldots,\nu_{h_{K,n}}]^{T}$,
$\hat{\bm{h}}_{n}=[\hat{h}_{1,n},\ldots,\hat{h}_{K,n}]^{T}$,
$\left[\overset{\scriptscriptstyle\rightarrow}{\bm{\nu}}_{\bm{h}_{n}}\right]_{(k-1)M+m}=\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{m,k,n}}$
and
$\left[{\overset{\scriptscriptstyle\rightarrow}{\bm{h}}_{n}}\right]_{(k-1)M+m}=\overset{\scriptscriptstyle\rightarrow}{h}_{m,k,n}$.

Similarly, the message $m_{\bm{g}_{n}\rightarrow
f_{\bm{s}_{n}}}\left(\bm{g}_{n}\right)$ is also Gaussian with mean
$\overset{\scriptscriptstyle\leftarrow}{\bm{g}}_{n}$ and variance
$\textbf{Diag}(\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{g}_{n}})$,
which can be computed as

$\displaystyle\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{g}_{n}}=\left(\mathbf{1}_{K}\otimes\left(\mathbf{1}_{M}./\bm{\nu}_{\bm{g}_{n}}\right)-\mathbf{1}_{J}./\overset{\scriptscriptstyle\rightarrow}{\bm{\nu}}_{\bm{g}_{n}}\right)^{.-1},$
(56)
$\displaystyle\overset{\scriptscriptstyle\leftarrow}{\bm{g}}_{n}=\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{g}_{n}}\cdot\left(\mathbf{1}_{K}\otimes(\hat{\bm{g}}_{n}./\bm{\nu}_{\bm{g}_{n}})-\overset{\scriptscriptstyle\rightarrow}{\bm{g}}_{n}./\overset{\scriptscriptstyle\rightarrow}{\bm{\nu}}_{\bm{g}_{n}}\right),$
(57)

where
$\bm{\nu}_{\bm{g}_{n}}=[{\nu}_{{g}_{1,n}},\ldots,{\nu}_{{g}_{M,n}}]^{T}$,
$\hat{\bm{g}}_{n}=[\hat{g}_{1,n},\ldots,\hat{g}_{M,n}]^{T}$,
$\left[\overset{\scriptscriptstyle\rightarrow}{\bm{\nu}}_{\bm{g}_{n}}\right]_{(m-1)K+k}=\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{m,k,n}}$
and
$\left[{\overset{\scriptscriptstyle\rightarrow}{\bm{g}}_{n}}\right]_{(m-1)K+k}=\overset{\scriptscriptstyle\rightarrow}{g}_{m,k,n}$.
Then, the backward message
$m_{f_{\bm{\tilde{s}}_{n}}\rightarrow\bm{\tilde{s}}_{n}}(\bm{\tilde{s}}_{n})=\mathcal{N}(\bm{\tilde{s}}_{n};\overset{\scriptscriptstyle\leftarrow}{\bm{\tilde{s}}}_{n},\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{\tilde{s}}_{n}})$
with

$\displaystyle\overset{\scriptscriptstyle\leftarrow}{\bm{\tilde{s}}}_{n}=\overset{\scriptscriptstyle\leftarrow}{\bm{h}}_{n}\cdot\overset{\scriptscriptstyle\leftarrow}{\bm{g}}_{n},$
(58)
$\displaystyle\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{\tilde{s}}_{n}}=|\overset{\scriptscriptstyle\leftarrow}{\bm{h}}_{n}|^{2}\cdot\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{g}_{n}}+\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{h}_{n}}\cdot|\overset{\scriptscriptstyle\leftarrow}{\bm{g}}_{n}|^{2}+\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{h}_{n}}\cdot\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{g}_{n}},$
(59)

where
$\overset{\scriptscriptstyle\leftarrow}{\bm{\tilde{s}}}_{n}=[\overset{\scriptscriptstyle\leftarrow}{\bm{\tilde{s}}}_{1,n}^{T},\ldots,\overset{\scriptscriptstyle\leftarrow}{\bm{\tilde{s}}}_{K,n}^{T}]^{T}$
and
$\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{\tilde{s}}_{n}}=[\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{\tilde{s}}_{1,n}}^{T},\ldots,\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{\tilde{s}}_{K,n}}^{T}]^{T}$.
Following (U)AMP, the backward message is combined with message
$m_{\bm{\tilde{s}}_{n}\rightarrow
f_{\bm{\tilde{s}}_{n}}}\left(\bm{\tilde{s}}_{n}\right)$, i.e.,

$\displaystyle\bm{\nu}_{\bm{\tilde{s}}_{k,n}}=\left(1/\nu_{\bm{q}_{k,n}}\mathbf{1}_{M}+\mathbf{1}_{M}./\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{\tilde{s}}_{k,n}}\right)^{.-1},$
(60)
$\displaystyle\hat{\bm{\tilde{s}}}_{k,n}=\bm{\nu}_{\bm{\tilde{s}}_{k,n}}\cdot\left(1/\nu_{\bm{\tilde{s}}_{k,n}}\hat{\bm{q}}_{k,n}+\overset{\scriptscriptstyle\leftarrow}{\bm{\tilde{s}}}_{k,n}./\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{\tilde{s}}_{k,n}}\right),$
(61)

Stack $\bm{\nu}_{\bm{\tilde{s}}_{k,n}}$ and $\hat{\bm{\tilde{s}}}_{k,n}$ into

$\displaystyle\bm{\nu}_{\bm{\tilde{S}}}=[\bm{\nu}_{\bm{\tilde{s}}_{1}},\ldots,\bm{\nu}_{\bm{\tilde{s}}_{N}}],$
(62)
$\displaystyle\hat{\bm{\tilde{S}}}=[\hat{\bm{\tilde{s}}}_{1},\ldots,\hat{\bm{\tilde{s}}}_{N}],$
(63)

where

$\displaystyle\bm{\nu}_{\bm{\tilde{s}}_{n}}=[\bm{\nu}_{\bm{\tilde{s}}_{1,n}}^{T},\ldots,\bm{\nu}_{\bm{\tilde{s}}_{K,n}}^{T}]^{T},$
(64)
$\displaystyle\hat{\bm{\tilde{s}}}_{n}=[\hat{\bm{\tilde{s}}}_{1,n}^{T},\ldots,\hat{\bm{\tilde{s}}}_{K,n}^{T}]^{T}.$
(65)

Further, we can obtain $\bm{\nu}_{\bm{S}}=\bm{\nu}_{\bm{\tilde{S}}}^{T}$ and
$\hat{\bm{S}}=\hat{\bm{\tilde{S}}}^{T}$ with

$\displaystyle\bm{\nu}_{\bm{S}}=[\bm{\nu}_{{\bm{s}}_{1}},\ldots,\bm{\nu}_{{\bm{s}}_{J}}],$
(66) $\displaystyle{\nu}_{\bm{s}_{j}}=\langle{\bm{\nu}_{\bm{s}_{j}}}\rangle,$
(67) $\displaystyle\hat{\bm{S}}=[\hat{\bm{s}}_{1},\ldots,\hat{{\bm{s}}}_{J}],$
(68) $\displaystyle\hat{{\bm{s}}}_{j}=[\hat{s}_{j,1},\ldots,\hat{s}_{j,N}],$
(69)

According to the belief propagation derivation of (U)AMP, it holds that

$m_{\bm{z}_{j}\rightarrow
f_{\bm{r}_{j}}}(\bm{z}_{j})=m_{f_{\bm{z}_{j}}\rightarrow\bm{z}_{j}}(\bm{z}_{j})=\mathcal{N}(\bm{z}_{j};\bm{p}_{j},\textbf{Diag}(\bm{\nu_{p_{j}}})),$
(70)

where $\bm{p}_{j}$ and $\bm{\nu_{p_{j}}}$ are computed in the forward
direction. It is noted that the factor node $f_{\bm{r}_{j}}$ connects the
variable node $\beta$. According to the rules of the variational message
passing, the message $m_{f_{\bm{r}_{j}}\rightarrow\beta}(\beta)$ can be
expressed as

$m_{f_{\bm{r}_{j}}\rightarrow\beta}(\beta)\propto\exp\left\\{\sum_{j=1}^{J}\int_{\bm{z}_{j}}\mathfrak{b}(\bm{z}_{j})\log
f_{\bm{r}_{j}}\right\\},$ (71)

where $\mathfrak{b}(\bm{z}_{j})$ is the approximate marginal of $\bm{z}_{j}$,
which can be expressed as

$\displaystyle\mathfrak{b}(\bm{z}_{j})$ $\displaystyle\propto
m_{f_{\bm{r}_{j}}\rightarrow\bm{z}_{j}}(\bm{z}_{j})m_{\bm{z}_{j}\rightarrow
f_{\bm{r}_{j}}}(\bm{z}_{j})$
$\displaystyle=\mathcal{N}\left(\bm{z}_{j};\hat{\bm{z}}_{j},\textbf{Diag}\left(\bm{\nu}_{\bm{z}_{j}}\right)\right),$
(72)

where

$\bm{\nu}_{\bm{z}_{j}}=\mathbf{1}_{L}./\left(\mathbf{1}_{L}\cdot/\bm{\nu}_{\bm{p}_{j}}+\hat{\beta}\mathbf{1}_{L}\right),$
(73)
$\hat{\bm{z}}_{j}=\bm{\nu}_{\bm{z}_{j}}\cdot\left(\bm{p}_{j}\cdot/\bm{\nu}_{\bm{p}_{j}}+\hat{\beta}\bm{r}_{j}\right),$
(74)

with $\hat{\beta}$ being the approximate a posteriori mean of the noise
precision that is obtained with (13). It is noted that in the above
derivation, the message $m_{f_{\bm{r}_{j}}\rightarrow\bm{z}_{j}}(\bm{z}_{j})$
is required, which is Gaussian, i.e.,
$m_{f_{\bm{r}_{j}}\rightarrow\bm{z}_{j}}(\bm{z}_{j})=\mathcal{N}(\bm{z}_{j};\bm{r}_{j},\hat{\beta}^{-1})$
, and its derivation is shown in (IV-A). Then, it is not hard to show that the
message

$m_{f_{\bm{r}_{j}}\rightarrow\beta}(\beta)\propto\beta^{LJ}\exp\left\\{\sum^{J}_{j=1}-\beta\left(\|\bm{r}_{j}-\hat{\bm{z}_{j}}\|^{2}+\mathbf{1}_{L}^{T}\bm{\nu}_{\bm{z}_{j}}\right)\right\\}.$
(75)

This is the end of backward message passing.

The message passing algorithm is summarized in Algorithm 1 and it can be
terminated when it reaches a maximum number of iteration or the difference
between the estimates of two consecutive iterations is less than a threshold.

Algorithm 1 UAMP-Based Channel Estimation for RIS-Aided MIMO System

Input: A feasible $\bm{\Phi}$, $\epsilon>0$ and the maximum number of
iteration $I_{max}$.
Initialize:
$\hat{h}_{k,n},\nu_{h_{k,n}}=1,\hat{{\bm{s}}}_{j}=\bm{0},{\nu}_{{\bm{s}}_{j}}=1,\bm{\mu}_{j}=\bm{0}$,
$\forall k,n,j$, and $\hat{\beta}=1$.
Repeat:

1: update noise precision $\hat{\beta}$ with (13);

2: $\forall j$: update $\bm{\nu_{p_{j}}}$ and $\bm{p}_{j}$ with (16) and (17);

3: $\forall j$: update $\bm{\nu}_{\bm{\mu}_{j}}$ and $\bm{\mu}_{j}$ with (18)
and (19);

4: $\forall j$: update $\bm{\nu}_{\bm{q}_{j}}$ and $\bm{q}_{j}$ with (20) and
(21);

5: $\forall n$: update $\bm{{\nu}}^{\prime}_{\bm{\tilde{q}}_{n}}$ and
$\bm{\tilde{q}}_{n}$ with (30) and (27);

6: $\forall m,k,n$: update
$\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{m,k,n}}$ and
$\overset{\scriptscriptstyle\rightarrow}{g}_{m,k,n}$ with (36) and (37);

7: $\forall m,n$: update
$\overset{\scriptscriptstyle\rightarrow}{\nu}_{g_{m,n}}$ and
$\overset{\scriptscriptstyle\rightarrow}{g}_{m,n}$ with (39) and (40);

8: $\forall m,n$: update $\nu_{g_{m,n}}$ and $\hat{g}_{m,n}$ with (42) and
(43);

9: $\forall m,k,n$: update
$\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{m,k,n}}$ and
$\overset{\scriptscriptstyle\rightarrow}{h}_{m,k,n}$ with (45) and (46);

10: $\forall k,n$: update
$\overset{\scriptscriptstyle\rightarrow}{\nu}_{h_{k,n}}$ and
$\overset{\scriptscriptstyle\rightarrow}{h}_{k,n}$ with (48) and (49);

11: $\forall k,n$: update $\nu_{h_{k,n}}$ and $\hat{h}_{k,n}$ with (51) and
(52);

12: $\forall n$: update
$\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{h}_{n}}$ and
$\overset{\scriptscriptstyle\leftarrow}{\bm{h}}_{n}$ with (54) and (55);

13: $\forall n$: update
$\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{g}_{n}}$ and
$\overset{\scriptscriptstyle\leftarrow}{\bm{g}}_{n}$ with (56) and (57);

14: $\forall n$: update
$\overset{\scriptscriptstyle\leftarrow}{\bm{\nu}}_{\bm{\tilde{s}}_{n}}$ and
$\overset{\scriptscriptstyle\leftarrow}{\bm{\tilde{s}}}_{n}$ with (59) and
(58);

15: $\forall k,n$: update $\bm{\nu}_{\bm{\tilde{s}}_{k,n}}$ and
$\hat{\bm{\tilde{s}}}_{k,n}$ with (60) and (61);

16: $\forall j$: update ${\nu}_{{\bm{s}}_{j}}$ and ${\hat{\bm{s}}_{j}}$ with
(67) and (69);

17: $\forall j$: update $\bm{\nu}_{\bm{z}_{j}}$ and $\hat{\bm{z}}_{j}$ with
(73) and (74);

18: Construct $\hat{\bm{H}}=[\hat{\bm{h}}_{1},\ldots,\hat{\bm{h}}_{N}]^{T}$
with $\hat{\bm{h}}_{n}=[\hat{{h}}_{1,n},\ldots,\hat{{h}}_{K,n}]^{T}$ and
$\hat{\bm{G}}=[\hat{\bm{g}}_{1},\ldots,\hat{\bm{g}}_{N}]$ with
$\hat{\bm{g}}_{n}=[\hat{{g}}_{1,n},\ldots,\hat{{g}}_{M,n}]^{T}$.

Until $||\hat{\bm{H}}-\bm{H}||_{F}^{2}\left\|\bm{H}\right\|_{F}^{-2}<\epsilon$
and $||\hat{\bm{G}}-\bm{G}||_{F}^{2}\left\|\bm{G}\right\|_{F}^{-2}<\epsilon$
or the number of iteration is more than $I_{max}$.

Output: $\hat{\beta}$, $\hat{\bm{H}}$ and $\hat{\bm{G}}$ that are the
estimations of ${\beta}$, ${\bm{H}}$ and ${\bm{G}}$, respectively.

###  IV-C Computational Complexity Analysis

We analyze the computational complexity of the proposed algorithm and compare
it with that of sate-of-the-art algorithms. The UAMP-based message passing
algorithm needs pre-processing, i.e., performing a single economic SVD for
$\bm{\Phi}$ and unitary transformation, and the complexity is
$\mathcal{O}(NL^{2})$. It noted that the SVD can be carried out offline and
there is no matrix inversion involved in Algorithm 1. Also note that the
formulated problem is a multiple measurement vector one. The complexity of the
proposed algorithm is dominated by the computation of $\bm{p}$ in step 1,
which requires $\mathcal{O}(NLKM+LKM)$, and the computations of $\nu_{\bm{q}}$
and $\bm{q}$ in step 6, which require $\mathcal{O}(NLKM+NKM)$ and
$\mathcal{O}(NLKM)$, respectively. It can be shown that the overall complexity
of the algorithm is $\mathcal{O}(NLKM)+\mathcal{O}(NKM)$ per iteration. As we
consider a general RIS-aided MIMO system without any special requirements on
the (channel or RIS phase) matrices, the most relevant algorithms for
comparison are the ALS-based algorithm and VAMP-based algorithm in [29]. The
complexity of the ALS-based algorithm is
$\mathcal{O}(NLKM+L(M+K)N^{2})+\mathcal{O}(N^{3})$ per iteration, and that of
the VAMP-based algorithm is $\mathcal{O}((M+K)N^{3}+(M+K)N^{2})$ per
iteration. It is noted that $N$ is the number of RIS units, which can be much
larger than $M$, $K$ and $L$. From the analysis, we can see that the
complexity of the proposed algorithm, which is linear with $N$, is
significantly smaller than that of the ALS or VAMP-based algorithm.

TABLE II: Complexity Comparison Algorithm |  Complexity
---|---
UAMP |  $\mathcal{O}(NLKM)+\mathcal{O}(NKM)$
ALS |  $\mathcal{O}(NLKM+L(M+K)N^{2})+\mathcal{O}(N^{3})$
VAMP |  $\mathcal{O}((M+K)N^{3}+(M+K)N^{2})$
$\displaystyle\bm{\mathcal{P}}=\begin{bmatrix}\mathbb{E}\left\\{\left(\frac{\partial
f}{\partial h_{1,1}}\right)\\!\left(\frac{\partial f}{\partial
h_{1,1}^{*}}\right)\right\\}&\ldots&\mathbb{E}\left\\{\left(\frac{\partial
f}{\partial h_{1,1}}\right)\left(\frac{\partial f}{\partial
h_{K,N}^{*}}\right)\right\\}&\mathbb{E}\left\\{\left(\frac{\partial
f}{\partial h_{1,1}}\right)\\!\left(\frac{\partial f}{\partial
g_{1,1}^{*}}\right)\right\\}&\ldots&\mathbb{E}\left\\{\left(\frac{\partial
f}{\partial h_{1,1}}\right)\\!\left(\frac{\partial f}{\partial
g_{M,N}^{*}}\right)\right\\}\\\ \vdots&&\vdots&\vdots&&\vdots\\\
\mathbb{E}\left\\{\left(\frac{\partial f}{\partial
h_{K,N}}\right)\\!\left(\frac{\partial f}{\partial
h_{1,1}^{*}}\right)\right\\}&\ldots&\mathbb{E}\left\\{\left(\frac{\partial
f}{\partial h_{K,N}}\right)\\!\left(\frac{\partial f}{\partial
h_{K,N}^{*}}\right)\right\\}&\mathbb{E}\left\\{\left(\frac{\partial
f}{\partial h_{K,N}}\right)\\!\left(\frac{\partial f}{\partial
g_{1,1}^{*}}\right)\right\\}&\ldots&\mathbb{E}\left\\{\left(\frac{\partial
f}{\partial h_{K,N}}\right)\left(\frac{\partial f}{\partial
g_{M,N}^{*}}\right)\right\\}\\\ \\\ \mathbb{E}\left\\{\left(\frac{\partial
f}{\partial g_{1,1}}\right)\\!\left(\frac{\partial f}{\partial
h_{1,1}^{*}}\right)\right\\}&\ldots&\mathbb{E}\left\\{\left(\frac{\partial
f}{\partial g_{1,1}}\right)\left(\frac{\partial f}{\partial
h_{K,N}^{*}}\right)\right\\}&\mathbb{E}\left\\{\left(\frac{\partial
f}{\partial g_{1,1}}\right)\\!\left(\frac{\partial f}{\partial
g_{1,1}^{*}}\right)\right\\}&\ldots&\mathbb{E}\left\\{\left(\frac{\partial
f}{\partial g_{1,1}}\right)\\!\left(\frac{\partial f}{\partial
g_{M,N}^{*}}\right)\right\\}\\\ \vdots&&\vdots&\vdots&&\vdots\\\
\mathbb{E}\left\\{\left(\frac{\partial f}{\partial
g_{M,N}}\right)\\!\left(\frac{\partial f}{\partial
h_{1,1}^{*}}\right)\right\\}&\ldots&\mathbb{E}\left\\{\left(\frac{\partial
f}{\partial g_{M,N}}\right)\\!\left(\frac{\partial f}{\partial
h_{K,N}^{*}}\right)\right\\}&\mathbb{E}\left\\{\left(\frac{\partial
f}{\partial g_{M,N}}\right)\\!\left(\frac{\partial f}{\partial
g_{1,1}^{*}}\right)\right\\}&\ldots&\mathbb{E}\left\\{\left(\frac{\partial
f}{\partial g_{M,N}}\right)\left(\frac{\partial f}{\partial
g_{M,N}^{*}}\right)\right\\}\\\ \end{bmatrix}.$ (86)

##  V Cramér-Rao Lower Bound

In this section, we derive the CRLB for the RIS channel estimation, which is
used to serve as another performance benchmark, besides the ALS and VAMP-based
algorithms.

We firstly rewrite the system model as

$\displaystyle\bm{\tilde{Y}}=(\bm{H}^{T}\odot\bm{G})\bm{\Phi}^{T}+\bm{\tilde{W}},$
(76)

and define a complex parameter $\bm{\theta}\in\mathbb{C}^{2N(M+K)\times 1}$,
which includes all of unknown complex parameters in $\bm{H}$ and $\bm{G}$ as

$\displaystyle\bm{\theta}\triangleq\left[\bm{h}_{1}^{T},\ldots,\bm{h}_{N}^{T},\bm{g}_{1}^{T},\ldots,\bm{g}_{N}^{T},\bm{h}_{1}^{H},\ldots,\bm{h}_{N}^{H},\bm{g}_{1}^{H},\ldots,\bm{g}_{N}^{H}\right]^{T}.$
(77)

The likelihood function of $\bm{\tilde{Y}}$ can be expressed as

$\displaystyle
p(\bm{\tilde{Y}};\bm{\theta})\\!=\\!\left(\pi\sigma^{2}\right)^{-KML}\\!\exp\\!\left\\{\\!-{\sigma^{-2}}\sum_{l=1}^{L}\left\|\bm{\tilde{y}}_{l}\\!-\\!(\bm{H}^{T}\\!\odot\\!\bm{G})\bm{\Phi}^{T}_{l,:}\right\|^{2}\right\\}\\!,$
(78)

and the logarithm of likelihood function can be expressed as

$\displaystyle\ln(p(\bm{\tilde{Y}};\bm{\theta}))\\!\\!=\\!\\!-KML\ln\left(\\!\pi\sigma^{2}\right)\\!\\!-\\!\\!{\sigma^{-2}}\sum_{l=1}^{L}\left\|\bm{\tilde{y}}_{l}\\!-\\!(\bm{H}^{T}\\!\\!\odot\\!\bm{G})\bm{\Phi}^{T}_{l,:}\right\|^{2}\\!\\!,$
(79)

Then, define $f_{\bm{\theta}}\triangleq\ln(p(\bm{\tilde{Y}};\bm{\theta}))$ and
the Fisher information matrix (FIM)
$\bm{\mathcal{J}}_{\bm{\theta}}\in\mathbb{C}^{(2N(M+K))\times(2N(M+K))}$ can
be obtained by

$\displaystyle\bm{\mathcal{J}}_{\bm{\theta}}=\mathbb{E}\left\\{\left(\frac{\partial
f_{\bm{\theta}}}{\partial\bm{\theta}}\right)\left(\frac{\partial
f_{\bm{\theta}}}{\partial\bm{\theta}}\right)^{H}\right\\}.$ (80)

The partial derivatives of $f_{\bm{\theta}}$ with respect to $\bm{\theta}$ can
be expressed as

$\displaystyle\frac{\partial f_{\bm{\theta}}}{\partial
h_{k,n}}\\!=\\!{\sigma^{-2}}\sum_{l=1}^{L}\\!\left\\{[\bm{\Phi}^{T}_{l,:}]_{n}\bm{g}^{T}_{n}\left(\\!\\!([\bm{Y}^{{}^{\prime}}_{l}]_{:,k}))^{*}\\!\\!-\\!((\bm{H}^{T})_{k,:}\\!\odot\\!\bm{G})^{*}(\bm{\Phi}^{T}_{l,:})^{*}\right)\\!\right\\},$
(81) $\displaystyle\frac{\partial f_{\bm{\theta}}}{\partial
g_{m,n}}\\!=\\!{\sigma^{-2}}\sum_{l=1}^{L}\left\\{[\bm{\Phi}^{T}_{l,:}]_{n}\bm{h}^{T}_{n}\left((\bm{Y}^{{}^{\prime}}_{l})_{m,:})^{H}\\!-\\!(\bm{H}^{T}\odot\bm{G}_{m,:})^{*}(\bm{\Phi}^{T}_{l,:})^{*}\right)\right\\},$
(82) $\displaystyle\frac{\partial f_{\bm{\theta}}}{\partial
h_{k,n}^{*}}=\left(\frac{\partial f_{\bm{\theta}}}{\partial
h_{k,n}}\right)^{*},\qquad\frac{\partial f_{\bm{\theta}}}{\partial
g_{m,n}^{*}}=\left(\frac{\partial f_{\bm{\theta}}}{\partial
g_{m,n}}\right)^{*}.$ (83)

where $[\bm{\Phi}^{T}_{l,:}]_{n}$ is the $n$-th element of
$\bm{\Phi}^{T}_{l,:}$,
$\bm{Y}^{{}^{\prime}}_{l}\triangleq[\bm{\tilde{y}}_{1,l},\ldots,\bm{\tilde{y}}_{K,l}]\in\mathbb{C}^{M\times
K}$ and $\bm{\tilde{y}}_{k,l}\in\mathbb{C}^{M\times 1}$ represents the sub-
vectors of $\bm{\tilde{y}}_{l}$ as

$\displaystyle\bm{\tilde{y}}_{l}=[\bm{\tilde{y}}_{1,l}^{T},\ldots,\bm{\tilde{y}}_{K,l}^{T}]^{T}.$
(84)

Hence, the FIM $\bm{\mathcal{J}}_{\bm{\theta}}$ can be expressed as

$\displaystyle\bm{\mathcal{J}}_{\bm{\theta}}=\left[\begin{array}[]{cc}\bm{\mathcal{P}}&\bm{0}\\\
\bm{0}&\bm{\mathcal{P}}^{*}\end{array}\right],$ (87)

where $\bm{\mathcal{P}}$ is shown in (86) at the top of next page and the size
of $\bm{\mathcal{P}}$ and $\bm{0}$ are $(N(M+K))\times(N(M+K))$. The inverse
of the FIM of $\bm{\theta}$ gives, under some regularity conditions, a lower
bound for the augmented covariance matrix of an unbiased estimator of
$\bm{\theta}$ as

$\displaystyle\bm{\mathcal{J}}_{\bm{\theta}}^{-1}=\left[\begin{array}[]{cc}\bm{\mathcal{P}}^{-1}&\bm{0}\\\
\bm{0}&(\bm{\mathcal{P}}^{-1})^{*}\end{array}\right],$ (89)

where $\bm{\mathcal{P}}^{-1}$ is given by

$\displaystyle\bm{\mathcal{P}}^{-1}=\left[\begin{array}[]{cc}\bm{\Omega}_{\bm{H}}&\bm{\mathcal{I}}\\\
\bm{\mathcal{I}}^{H}&\bm{\Omega}_{\bm{G}}\end{array}\right],$ (92)

and $\bm{\Omega}_{\bm{H}}\in\mathbb{C}^{KN\times KN}$ and
$\bm{\Omega}_{\bm{G}}\in\mathbb{C}^{MN\times MN}$ are the CRLB matrices for
the estimates of $\bm{H}$ and $\bm{G}$, respectively, and
$\bm{\mathcal{I}}\in\mathbb{C}^{KN\times MN}$ represents the remaining sub-
matrices.

Furthermore, $\bm{\mathcal{P}}$ can be divided into four sub-matrices as

$\displaystyle\bm{\mathcal{P}}=\left[\begin{array}[]{cc}\bm{\mathcal{P}}_{\bm{HH}}&\bm{\mathcal{P}}_{\bm{HG}}\\\
\bm{\mathcal{P}}_{\bm{HG}}^{H}&\bm{\mathcal{P}}_{\bm{GG}}\end{array}\right],$
(95)

where $\bm{\mathcal{P}}_{\bm{HH}}\in\mathbb{C}^{KN\times KN}$,
$\bm{\mathcal{P}}_{\bm{GG}}\in\mathbb{C}^{MN\times MN}$ and
$\bm{\mathcal{P}}_{\bm{HG}}\in\mathbb{C}^{KN\times MN}$, as shown in (86).
According to the formula for inverse of a partitioned Hermitian matrix in
[36], we can obtain $\bm{\Omega}_{\bm{H}}$ and $\bm{\Omega}_{\bm{G}}$ as
following

$\displaystyle\bm{\Omega}_{\bm{H}}=$
$\displaystyle(\bm{\mathcal{P}}_{\bm{HH}}-\bm{\mathcal{P}}_{\bm{HG}}\bm{\mathcal{P}}_{\bm{GG}}^{-1}\bm{\mathcal{P}}_{\bm{HG}}^{H})^{-1},$
(96) $\displaystyle\bm{\Omega}_{\bm{G}}=$
$\displaystyle(\bm{\mathcal{P}}_{\bm{GG}}-\bm{\mathcal{P}}_{\bm{HG}}^{H}\bm{\mathcal{P}}_{\bm{HH}}^{-1}\bm{\mathcal{P}}_{\bm{HG}})^{-1}.$
(97)

So the CRLB of $\bm{H}$ and $\bm{G}$ can be donated as

$\displaystyle\mathrm{CRLB}_{\bm{H}}=\frac{\mathrm{trace}(\bm{\Omega}_{\bm{H}})}{KN},$
(98)
$\displaystyle\mathrm{CRLB}_{\bm{G}}=\frac{\mathrm{trace}(\bm{\Omega}_{\bm{G}})}{MN}.$
(99)

##  VI Simulation Results

In this section, we provide extensive numerical experiments to demonstrate the
superior performance of the proposed UAMP-based channel estimation algorithm.
For comparison, we also include the ALS-based and VAMP-based channel
estimation algorithms in [29]. The threshold $\epsilon=10^{-3}$ and the
maximum number of iterations is set to 30. The entries of $\bm{H}$ and
$\bm{G}$ are independently drawn from a complex Gaussian distribution with
zero mean and unit variance. The scaling ambiguity of the estimation is
eliminated in the calculation of the normalized mean square error (NMSE). The
SNR (in dB) is defined as

$\displaystyle{\rm SNR}=$ $\displaystyle
10\log_{10}\left(\frac{N\mathbb{E}\left\\{\left\|\bm{\Phi}\right\|_{F}^{2}\right\\}}{L\beta^{-1}}\right).$
(100)

In order to estimate the RIS channels during the training phase, the RIS goes
through $L$ configurations, which leads to the RIS phase matrix $\bm{\Phi}$ in
(6). In the simulations, we consider two methods. One is that the RIS units
are turned on or off randomly, leading to a matrix $\bm{\Phi}$ with entries 1
or 0 [37], which is called a binary matrix hereafter. In the simulations we
assume that each entry in $\bm{\Phi}$ takes 1 or 0 with the same probability.
The other one is that the phases of the RIS units are set to some discrete
values, and in particular, the phase matrix $\bm{\Phi}$ is part of the DFT
matrix (called partial DFT matrix), as in [22].

We evaluate the performance of estimators in terms of the NMSE of estimated
channel matrices $\hat{\bm{H}}$ and $\hat{\bm{G}}$. The NMSE performance of
various estimators versus SNR with $L=N=K=M=64$ is shown in Fig. 5, where Fig.
5 (a) and (b) are for partial DFT matrix and binary phase matrix,
respectively. It is observed that the proposed UAMP-based algorithm
significantly outperforms the ALS-based and VAMP-based channel estimators for
both $\bm{H}$ and $\bm{G}$, especially in the case that $\bm{\Phi}$ is a
binary matrix. Note that the NMSEs of $\bm{H}$ and $\bm{G}$ are very similar.
To keep the figures clear, we only show the NMSE performance of $\bm{H}$ in
the subsequent simulation results.

Figure 5:  NMSE performance of the estimators versus SNR, where $L=N=K=M=64$.
(a) Partial DFT matrix; (b) Binary matrix. Figure 6: NMSE performance of the
estimators versus SNR for $N=K=M=64$ and different $L$ (a) Partial DFT matrix;
(b) Binary matrix. Figure 7: NMSE performance of the estimators versus $L$,
where $N=K=M=64$ and SNR=20dB. (a) Partial DFT matrix; (b) Binary matrix.

In Fig. 6, with $N=K=M=64$, we compare the NMSE performance versus SNR of the
estimators with different values of $L$. According to the results, the
performance of the UAMP-based method is significantly better than that of the
ALS and VAMP-based methods, especially when $L$ is relatively small. As
smaller $L$ (the number of RIS phase configurations needed for channel
estimation) is highly desirable to reduce the training overhead and latency,
next we vary the value of $L$ and examine the performance of the estimators.
The results are shown in Fig. 7, where the SNR is set to 20dB, and $\bm{\Phi}$
is partial DFT matrix in (a) and binary matrix in (b). It can be seen that,
with the increase of $L$, the performance of all estimators improves, as
expected. However, the performance improvement of the ALS and VAMP-based
estimators is very slow with $L$. We can also see from Fig. 7 (a) that, the
UAMP-based estimator with $L=16$ achieves the same performance of ALS and
VAMP-based estimators with $L=32$. According to Fig. 7 (b), the UAMP-based
estimator with $L=16$ even outperforms the ALS and VAMP-based estimators with
$L=32$. The results demonstrate that the use of the proposed algorithm can
lead to a huge reduction in training overheads.

Figure 8: NMSE performance of $\hat{\bm{H}}$ versus SNR with different $N$,
where $L=K=M=32$. (a) Partial DFT matrix; (b) Binary matrix. Figure 9:
Performance of the estimators and CRLB, where $L=K=M=N=16$. (a) Partial DFT
matrix; (b) Binary matrix. Figure 10: NMSE performance of $\hat{\bm{H}}$ and
$\hat{\bm{G}}$ versus the iteration number, where $L=20$ and $K=M=N=32$.
Figure 11: Noise variance estimation of the proposed algorithm.

We next examine the impact of the number of RIS units $N$ on the performance
of channel estimation, where we set $L=K=M=32$. As a smaller $N$ leads to a
less number of channel coefficients to be estimated, for a fixed $L$ the
performance of the estimators improves with the decrease of $N$. In Fig. 8, we
can find that, in the case of binary matrix, the performance of the UAMP-based
estimator with $N=128$ is even better than that of the ALS and VAMP-based
estimators with $N=32$, which again demonstrates the superior performance of
the proposed one. The performance of the UAMP-based algorithm, the ALS
algorithm and the CRLB are shown in Fig. 9, where $L=K=M=N=16$. We can see
that the performance of UAMP-based algorithm is almost the same as the CRLB,
which is significantly better than that of the ALS algorithm.

From the above results, we find that the performance of the proposed UAMP-
based algorithm consistently show good performance for both partial DFT matrix
and binary matrix. In contrast, the ALS and VAMP-based algorithms exhibit
significantly worse performance in the case of binary matrix, compared to the
partial DFT matrix.

We also investigate the performance of the estimation with iteration number
for different $\bm{\Phi}$ when the SNR is 20dB. The results are shown in Fig.
10. It can be seen that the proposed algorithm converges fast in different
cases, especially when $\bm{\Phi}$ is a partial DFT matrix. It is worth
mentioning that the estimation of the noise precision (reciprocal of the
variance) is incorporated in the UAMP-based channel estimator, so no separate
noise power estimator is needed. In Fig. 11, we compare the estimated noise
power and its true value, where $\bm{\Phi}$ is a partial DFT matrix. In Fig.
11 (a), the estimate of the noise variance with the iteration number is shown,
where the SNR is 20dB. We can see that the convergence is fast. The results in
Fig. 11 (b) show that the proposed algorithm provides accurate noise variance
estimates for a wide range of SNRs.

##  VII Conclusions

In this paper, we have addressed the issue of channel estimation in RIS-aided
MIMO communications. Through vectorization and reduction, we obtain a new
signal model for channel estimation, based on which a message passing based
algorithm is developed, leveraging UAMP. Compared to the state-of-the-art
algorithms, the proposed algorithm does not have any special requirements on
the matrices involved, and it shows significant advantages in computational
complexity, estimation performance and training overhead. Extensive numerical
results demonstrate the merits of the proposed algorithm.

## References

  * [1] 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.
  * [2] X. Guan, Q. Wu, and R. Zhang, “Intelligent reflecting surface assisted secrecy communication via joint beamforming and jamming,” _CoRR_ , vol. abs/1907.12839, 2019. [Online]. Available: http://arxiv.org/abs/1907.12839
  * [3] Q. Nadeem, A. Kammoun, A. Chaaban, M. Debbah, and M. Alouini, “Asymptotic analysis of large intelligent surface assisted MIMO communication,” _CoRR_ , vol. abs/1903.08127, 2019. [Online]. Available: http://arxiv.org/abs/1903.08127
  * [4] Y. Liang, R. Long, Q. Zhang, J. Chen, H. V. Cheng, and H. Guo, “Large intelligent surface/antennas (LISA): making reflective radios smart,” _CoRR_ , vol. abs/1906.06578, 2019. [Online]. Available: http://arxiv.org/abs/1906.06578
  * [5] Q. Wu and R. Zhang, “Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network,” _IEEE Communications Magazine_ , vol. 58, no. 1, pp. 106–112, 2020.
  * [6] Z. Yang, M. Chen, W. Saad, W. Xu, M. Shikh-Bahaei, H. V. Poor, and S. Cui, “Energy-efficient wireless communications with distributed reconfigurable intelligent surfaces,” _IEEE Transactions on Wireless Communications_ , pp. 1–1, 2021.
  * [7] H. Guo, Y.-C. Liang, J. Chen, and E. G. Larsson, “Weighted sum-rate maximization for intelligent reflecting surface enhanced wireless networks,” in _2019 IEEE Global Communications Conference (GLOBECOM)_ , 2019, pp. 1–6.
  * [8] Y. Liang, J. Yang, W. Xie, M. Hasna, and M. Renzo, “Secrecy performance analysis of RIS-aided wireless communication systems,” _IEEE Transactions on Vehicular Technology_ , vol. PP, no. 99, pp. 1–1, 2020.
  * [9] 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.
  * [10] E. Basar, “Reconfigurable intelligent surface-based index modulation: A new beyond mimo paradigm for 6G,” _IEEE Transactions on Communications_ , vol. 68, no. 5, pp. 3187–3196, 2020.
  * [11] C. Huang, S. Hu, G. C. Alexandropoulos, A. Zappone, C. Yuen, R. Zhang, M. D. Renzo, and M. Debbah, “Holographic MIMO surfaces for 6G wireless networks: Opportunities, challenges, and trends,” _IEEE Wireless Communications_ , vol. 27, no. 5, pp. 118–125, 2020.
  * [12] A. Zappone, M. Di Renzo, F. Shams, X. Qian, and M. Debbah, “Overhead-aware design of reconfigurable intelligent surfaces in smart radio environments,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 1, pp. 126–141, 2021.
  * [13] L. Yang, Y. Yang, M. O. Hasna, and M.-S. Alouini, “Coverage, probability of SNR gain, and DOR analysis of RIS-aided communication systems,” _IEEE Wireless Communications Letters_ , vol. 9, no. 8, pp. 1268–1272, 2020\.
  * [14] 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.
  * [15] Z. Zhu, Z. Li, Z. Chu, G. Sun, W. Hao, P. Liu, and I. Lee, “Resource allocation for intelligent reflecting surface assisted wireless powered iot systems with power splitting,” _IEEE Transactions on Wireless Communications_ , pp. 1–1, 2021.
  * [16] C. Huang, Z. Yang, G. C. Alexandropoulos, K. Xiong, L. Wei, C. Yuen, Z. Zhang, and M. Debbah, “Multi-hop ris-empowered terahertz communications: A drl-based hybrid beamforming design,” _IEEE Journal on Selected Areas in Communications_ , vol. 39, no. 6, pp. 1663–1677, 2021.
  * [17] T. Hou, Y. Liu, Z. Song, X. Sun, Y. Chen, and L. Hanzo, “Reconfigurable intelligent surface aided noma networks,” _IEEE Journal on Selected Areas in Communications_ , vol. 38, no. 11, pp. 2575–2588, 2020.
  * [18] I. F. Akyildiz, C. Han, and S. Nie, “Combating the distance problem in the millimeter wave and terahertz frequency bands,” _IEEE Communications Magazine_ , vol. 56, no. 6, pp. 102–108, 2018.
  * [19] 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.
  * [20] D. Mishra and H. Johansson, “Channel estimation and low-complexity beamforming design for passive intelligent surface assisted MISO wireless energy transfer,” in _ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2019, pp. 4659–4663.
  * [21] B. Zheng and R. Zhang, “Intelligent reflecting surface-enhanced ofdm: Channel estimation and reflection optimization,” _IEEE Wireless Communications Letters_ , vol. 9, no. 4, pp. 518–522, 2020.
  * [22] T. L. Jensen and E. De Carvalho, “An optimal channel estimation scheme for intelligent reflecting surfaces based on a minimum variance unbiased estimator,” in _ICASSP 2020 - 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2020, pp. 5000–5004.
  * [23] 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.
  * [24] A. Taha, M. Alrabeiah, and A. Alkhateeb, “Enabling large intelligent surfaces with compressive sensing and deep learning,” _IEEE Access_ , vol. 9, pp. 44 304–44 321, 2021.
  * [25] Z.-Q. He and X. Yuan, “Cascaded channel estimation for large intelligent metasurface assisted massive MIMO,” _IEEE Wireless Communications Letters_ , vol. 9, no. 2, pp. 210–214, 2020.
  * [26] Z.-Q. He, H. Liu, X. Yuan, Y.-J. A. Zhang, and Y.-C. Liang, “Semi-blind cascaded channel estimation for reconfigurable intelligent surface aided massive mimo,” 2021.
  * [27] H. Liu, X. Yuan, and Y.-J. A. Zhang, “Matrix-calibration-based cascaded channel estimation for reconfigurable intelligent surface assisted multiuser MIMO,” _IEEE Journal on Selected Areas in Communications_ , vol. 38, no. 11, pp. 2621–2636, 2020.
  * [28] Z. Wang, L. Liu, and S. Cui, “Channel estimation for intelligent reflecting surface assisted multiuser communications: Framework, algorithms, and analysis,” _IEEE Transactions on Wireless Communications_ , vol. 19, no. 10, pp. 6607–6620, 2020.
  * [29] L. Wei, C. Huang, G. C. Alexandropoulos, C. Yuen, Z. Zhang, and M. Debbah, “Channel estimation for RIS-empowered multi-user MISO wireless communications,” _IEEE Transactions on Communications_ , vol. 69, no. 6, pp. 4144–4157, 2021.
  * [30] L. Wei, C. Huang, Q. Guo, Z. Zhang, M. Debbah, and C. Yuen, “Bidirectional approximate message passing for ris-assisted multi-user miso communications,” in _2021 IEEE 94th Vehicular Technology Conference (VTC2021-Fall)_ , 2021, pp. 1–6.
  * [31] Q. Guo and J. Xi, “Approximate message passing with unitary transformation,” _CoRR_ , vol. abs/1504.04799, 2015. [Online]. Available: http://arxiv.org/abs/1504.04799
  * [32] Z. Yuan, Q. Guo, and M. Luo, “Approximate message passing with unitary transformation for robust bilinear recovery,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 617–630, 2021.
  * [33] M. Luo, Q. Guo, M. Jin, Y. C. Eldar, D. Huang, and X. Meng, “Unitary approximate message passing for sparse bayesian learning,” _IEEE Transactions on Signal Processing_ , vol. 69, pp. 6023–6039, 2021.
  * [34] D. L. Donoho, A. Maleki, and A. Montanari, “Message passing algorithms for compressed sensing: I. motivation and construction,” in _2010 IEEE Information Theory Workshop on Information Theory (ITW 2010, Cairo)_ , 2010, pp. 1–5.
  * [35] J. Dauwels, “On variational message passing on factor graphs,” in _2007 IEEE International Symposium on Information Theory_ , 2007, pp. 2546–2550.
  * [36] X. Liu and N. Sidiropoulos, “Cramer-rao lower bounds for low-rank decomposition of multidimensional arrays,” _IEEE Transactions on Signal Processing_ , vol. 49, no. 9, pp. 2074–2086, 2001.
  * [37] D. Mishra and H. Johansson, “Channel estimation and low-complexity beamforming design for passive intelligent surface assisted MISO wireless energy transfer,” in _ICASSP 2019 - 2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2019, pp. 4659–4663.