# Channel Estimation for IRS-aided Multiuser Communications with Reduced Error
Propagation

Yi Wei, Ming-Min Zhao, Min-Jian Zhao, and Yunlong Cai The authors are with the
College of Information Science and Electronic Engineering, Zhejiang
University, Hangzhou 310027, China (email: {21731133, zmmblack, mjzhao,
ylcai}@zju.edu.cn).

###### Abstract

Intelligent reflecting surface (IRS) has emerged as a promising paradigm to
improve the capacity and reliability of a wireless communication system by
smartly reconfiguring the wireless propagation environment. To achieve the
promising gains of IRS, the acquisition of the channel state information (CSI)
is essential, which however is practically difficult since the IRS does not
employ any transmit/receive radio frequency (RF) chains in general and it has
limited signal processing capability. In this paper, we study the uplink
channel estimation problem for an IRS-aided multiuser single-input multi-
output (SIMO) system. The existing channel estimation approach for IRS-aided
multiuser systems mainly consists of three phases, where the direct channels
from the base station (BS) to all the users, the reflected channel from the BS
to a typical user via the IRS, and the other reflected channels are estimated
sequentially based on the estimation results of the previous phases. However,
this approach will lead to a serious error propagation issue, i.e., the
channel estimation errors in the first and second phases will deteriorate the
estimation performance in the second and third phases. To resolve this
difficulty, we propose a novel two-phase channel estimation (2PCE) strategy
which is able to alleviate the negative effects caused by error propagation
and enhance the channel estimation performance with the same amount of channel
training overhead as in the existing approach. Specifically, in the first
phase, the direct and reflected channels associated with a typical user are
estimated simultaneously by varying the reflection patterns at the IRS, such
that the estimation errors of the direct channel associated with this typical
user will not affect the estimation of the corresponding reflected channel. In
the second phase, we estimate the CSI associated with the other users and
demonstrate that by properly designing the pilot symbols of the users and the
reflection patterns at the IRS, the direct and reflected channels associated
with each user can also be estimated simultaneously, which helps to reduce the
error propagation. Moreover, the asymptotic mean squared error (MSE) of the
proposed 2PCE strategy is analyzed when the least-square (LS) channel
estimation method is employed, and we show that the 2PCE strategy can
outperform the existing approach. Finally, extensive simulation results are
presented to validate the effectiveness of our proposed channel estimation
strategy.

###### Index Terms:

Intelligent reflecting surface, channel estimation, multiuser communications,
single-input multiple-output (SIMO).

##  I Introduction

To meet the ever-increasing demand for high-speed and low-latency data
applications, various advanced technologies, e.g., massive multiple-input
multiple-output (MIMO) communication, millimeter wave (mmWave) communication
and ultra-dense network (UDN), have been proposed and extensively studied to
improve the spectral efficiency of wireless systems [1, 2]. These technologies
have led to fundamental changes in the design of the fifth generation (5G)
cellular network and brought significant performance improvement, however,
they also incur high energy consumption and hardware cost, due to the large
number of active nodes/antennas/radio-frequency (RF) chains used. Recently,
intelligent reflecting surface (IRS) has been proposed as a promising new
paradigm to alleviate the above issues and enhance the system performance [3,
4, 5, 6]. In general, IRS is a man-made planar surface consisting of a large
number of passive and low-cost reflecting elements, each of which can be
independently adjusted to change the amplitude and/or phase of the incident
signal by an external software controller. As a result, by smartly
coordinating these reflecting elements, IRS is able to create a favorable
wireless signal propagation environment to improve the wireless communication
coverage, throughput, and energy efficiency substantially. Besides, different
from the conventional half-duplex amplify-and-forward (AF) relay, IRS can
achieve high beamforming gains by intelligently adjusting its reflection
coefficients at different reflecting elements in a full-duplex manner, with
low hardware and energy cost.

As such, IRS has attracted significant attention recently and various works
have been conducted in the literature to show the effectiveness of IRS in
enhancing the system performance, see e.g., [7, 8, 9, 10]. However, to reap
the performance gains offered by IRS, the acquisition of channel state
information (CSI) is of significant importance since the IRS is generally not
equipped with any transmit/receive RF chains and thus not capable of
performing complex baseband signal processing tasks. This issue becomes more
stringent when the numbers of base station (BS) antennas, IRS reflecting
elements and users are large, since the number of channel coefficients that
are required to be estimated will become considerably large. Although there
are several studies in the literature proposed to design the IRS passive
beamforming based on the statistical CSI and employ the two-timescale
transmission protocol to reduce the channel training overhead [11], how to
estimate the full CSI efficiently is still an essential problem since various
works require full CSI for IRS reflection optimization and obtaining a larger
number of channel samples is also an indispensable part for statistical CSI
estimation.

Recently, various strategies have been proposed to tackle the channel
estimation problem efficiently for IRS-aided wireless communication systems
[12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22]. Specifically, for the single-
user system, the studies of [12],[13] proposed binary reflection controlled
channel estimation schemes, where the user-BS channel (referred to as the
_direct channel_) is estimated first when all the reflecting elements are
turned off, then only one reflecting element is switched on at each subsequent
time slot to estimate the cascaded user-IRS-BS channel (referred to as the
_reflected channel_) associated with this reflecting element. In this
strategy, as only a small portion of the IRS reflecting elements is switched
on at each time, the channel estimation accuracy may be degraded as the IRS’s
large aperture gain is not exploited. To address this issue, the works [14]
and [15] improved the binary reflection strategy by considering the full
reflection of the IRS at all time slots and designing a discrete Fourier
transform (DFT) based IRS training reflection matrix (also known as
_reflection pattern_). Further, the work [17] studied the reflection pattern
design problem under the practical constraint that the phase shift at each
reflecting element only takes a finite number of discrete values. When massive
MIMO system is considered, the authors in [18] exploited the low-rank
structure of the channel matrix and formulated the reflected channel
estimation problem as a combined sparse matrix factorization and matrix
completion problem. As for the IRS-aided mmWave MIMO systems, the works [20,
21] studied the channel estimation problem by exploiting the channel sparsity.
On the other hand, for the multiuser system, the work [19] proposed a
compressive sensing (CS) based channel estimation method with reduced channel
training overhead, based on the assumption that the reflected channels are
sparse. For more general channel models, the work [22] proposed a novel three-
phase channel estimation (3PCE) framework to estimate the direct and reflected
channels, where the strong correlation among the reflected channels associated
with different users is exploited to reduce the channel training overhead.
Specifically, since each IRS reflecting element reflects the signals from
different users to the BS via the same IRS-BS channel, the reflected channel
of an arbitrary user is a scaled version of that of the other users and only
the scaling factors (scalars), rather than the whole channels (vectors), are
needed to be estimated. Under this framework, only
$K+N+\max(K-1,\lceil\frac{(K-1)N}{M}\rceil)$ pilot symbols are required to
estimate the full CSI which consists of $KM+KMN$ channel coefficients, instead
of $K+KN$ pilot symbols if this correlation is not considered, where $K$, $N$
and $M$ denote the numbers of users, IRS reflecting elements, and BS antennas,
respectively. However, these exist a serious error propagation issue in this
3PCE framework, i.e., the channel estimation errors in the first and second
phases will deteriorate the estimation performance in the second and third
phases.

Motivated by the above, we consider an IRS-aided multiuser uplink SIMO
communication system in this paper where multiple single-antenna users
communicate with a multi-antenna BS with the assist of an IRS, and propose a
new channel estimation strategy to alleviate the error propagation issue in
the existing 3PCE strategy [22]. Our main idea is to reduce the channel
estimation phases by estimating the direct and reflected channels associated
with each user simultaneously, such that the error propagation between
different phases can be well-controlled. Specifically, our proposed strategy
which is referred to as the two-phase channel estimation (2PCE) strategy
mainly consists of two estimation phases, where the first phase aims to
estimate the direct and reflected channels associated with a typical user and
the second phase is dedicated to estimate the other channels based on the
knowledge of this typical user’s CSI. We first show that the proposed 2PCE
strategy is able to perfectly estimate all the channel coefficients with the
theoretically minimal pilot sequence length, under the ideal case without
receive noise at the BS. Then, for the practical case with receive noise at
the BS, the least-square (LS) channel estimation solutions are derived in both
phases. The asymptotic MSE of the proposed 2PCE strategy is analyzed when the
number of transmit antennas at the BS becomes large to evaluate the
performance of our proposed 2PCE strategy. The analysis shows that the
proposed strategy is able to outperform the existing 3PCE strategy with the
same channel training overhead. Finally, extensive numerical results are
presented to validate the effectiveness of the proposed 2PCE strategy under
various numbers of BS antennas, IRS reflecting elements and users, and we show
that the performance of the proposed 2PCE strategy is better than that of the
3PCE strategy, even when the number of BS antennas is not large. The impacts
of various system parameters, such as the channel training signal-to-noise
ratio (SNR), the path-loss exponents, the channel Rician factors and the
correlation coefficients on the channel estimation performance are
investigated and useful insights are drawn.

The rest of this paper is organized as follows. Section II presents the
considered IRS-aided uplink multiuser SIMO system model. Section III
introduces the details of the proposed 2PCE strategy under the ideal case
without receive noise at the BS. Theoretical analysis on the asymptotic MSE of
the proposed 2PCE strategy is presented in Section IV and comparison with the
existing 3PCE strategy is provided. Numerical results are presented in Section
V and finally, Section VI concludes this paper.

Notation: Scalars, vectors and matrices are respectively denoted by lower
(upper) case, boldface lower case and boldface upper case letters. The symbol
$[\mathbf{X}]_{n,m}$ denotes the $(n,m)$-th entry of the matrix $\mathbf{X}$
while $[\mathbf{x}]_{n}$ represents the $n$-th entry of the vector
$\mathbf{x}$. The symbols $[\cdot,\cdot]$ and $[\cdot;\cdot]$ represent the
column and row concentrations of two scalars/vectors/matrices. For a matrix
$\mathbf{X}$ of arbitrary size, $\mathbf{X}^{T}$, $\mathbf{X}^{*}$,
$\mathbf{X}^{H}$ and $\mathbf{X}^{\dagger}$ denote the transpose, conjugate,
conjugate transpose and pseudo-inverse of $\mathbf{X}$, respectively. For a
square matrix $\mathbf{Y}$, $\text{tr}(\mathbf{Y})$ represents the trace of
$\mathbf{Y}$, and $\mathbf{Y}^{-1}$ denotes the inverse of $\mathbf{Y}$ if
$\mathbf{Y}$ is full-rank. The symbol $||\cdot||$ denotes the Euclidean norm
of a complex vector, and $|\cdot|$ is the absolute value of a complex scalar.
The symbol $\text{diag}(x_{1},\cdots,x_{N})$ denotes a diagonal matrix whose
diagonal elements are set as $x_{1},\cdots,x_{N}$. ${\mathbb{C}}^{m\times n}$
denotes the space of $m\times n$ complex matrices. $\lceil x\rceil$ and
$\lfloor x\rfloor$ represent the minimum integer which is larger than $x$ and
the maximum integer that is smaller than $x$, respectively. The symbol $\%$
represents the remainder operation, i.e.,
$x_{1}\%x_{2}=x_{1}-x_{2}\lfloor\frac{x_{1}}{x_{2}}\rfloor$. $\mathbf{I}_{M}$,
$\mathbf{0}_{M\times N}$ and $\mathbf{1}_{M\times N}$ denote the $M\times M$
identity matrix, $M\times N$ all-zero matrix and $M\times N$ all-one matrix,
respectively. Finally, we define the complex normal distribution as
$\mathcal{C}\mathcal{N}(\mu,{\sigma}^{2})$ with mean $\mu$ and variance
${\sigma}^{2}$.

##  II System Model

Figure 1: An IRS-aided multiuser SIMO communication system.

We consider an IRS-aided uplink multiuser SIMO system as shown in Fig. 1,
where $K$ single-antenna users, denoted by $\\{U_{1},\cdots,U_{K}\\}$,
simultaneously communicate with a BS, and an IRS is deployed to enhance the
communication performance. The BS is equipped with a uniform linear array
(ULA) with $M$ antennas, and the IRS is equipped with an $N_{y}\times N_{z}$
uniform planar array (UPA) which contains $N=N_{y}N_{z}$ passive reflecting
elements. Let $\mathbf{h}_{d,k}\in\mathbb{C}^{M\times 1}$ and
$\mathbf{h}_{r,k}\in\mathbb{C}^{N\times 1}$ denote the baseband equivalent
channels from $U_{k}$ to the BS and the IRS, respectively. The symbol
${\mathbf{G}}\in\mathbb{C}^{M\times N}$ represents the channel from the IRS to
the BS with $\mathbf{g}_{n}$ as the $n$-th column of ${\mathbf{G}}$.

In this work, a quasi-static block-fading channel is assumed and all the
channels $\mathbf{h}_{d,k}$’s, $\mathbf{h}_{r,k}$’s, $\mathbf{G}$ remain
approximately constant over a coherence interval of length $T$ symbols. Since
both line-of-sight (LoS) and non-LoS (NLoS) components may exist in practical
channels due to the insufficient angular spread of the scattering environment
and closely spaced antennas/reflecting elements, $\mathbf{h}_{d,k}$ can be
modeled as the general spatially correlated Rician fading channel [23] as

${\mathbf{h}}_{d,k}=\sqrt{\frac{\beta^{\text{UB}}}{1+\beta^{\text{UB}}}}{\bar{\mathbf{z}}}_{d,k}+\sqrt{\frac{1}{1+\beta^{\text{UB}}}}{\bm{\Phi}}_{d}^{\frac{1}{2}}{{\mathbf{z}}}_{d,k},$
(1)

where $\beta^{\text{UB}}$ denotes the Rician factor;
${\bar{\mathbf{z}}}_{d,k}\in{\mathbb{C}}^{M\times 1}$ denotes LoS component in
${\mathbf{h}}_{d,k}$;
$[{{\mathbf{z}}}_{d,k}]_{m}\sim\mathcal{C}\mathcal{N}\left(0,l_{k,m}^{\text{UB}}\right)$
follows the independent and identically distributed (i.i.d.) Rayleigh fading
channel model with $\l_{k,m}^{\text{UB}}$ denoting the distance-dependent path
loss of the link from $U_{k}$ to the $m$-th BS antenna; ${\bm{\Phi}}_{d}$ is
the BS receive correlation matrix. Similarly, the user-IRS and IRS-BS channels
can be modeled as follows:

${\mathbf{h}}_{r,k}=\sqrt{\frac{\beta^{\text{UI}}}{1+\beta^{\text{UI}}}}{\bar{\mathbf{z}}}_{r,k}+\sqrt{\frac{1}{1+\beta^{\text{UI}}}}{\bm{\Phi}}_{r,k}^{\frac{1}{2}}{{\mathbf{z}}}_{r,k},\;\;{\mathbf{G}}=\sqrt{\frac{\beta^{\text{IB}}}{1+\beta^{\text{IB}}}}\bar{\mathbf{F}}+\sqrt{\frac{1}{1+\beta^{\text{IB}}}}{\bm{\Phi}}_{d}^{\frac{1}{2}}{\mathbf{F}}{\bm{\Phi}}_{r}^{\frac{1}{2}},$
(2)

where $\beta^{\text{UI}}$ and $\beta^{\text{IB}}$ denote the corresponding
Rician factors; ${\bm{\Phi}}_{r,k}$ represents the correlation matrix between
the IRS and $U_{k}$, and ${\bm{\Phi}}_{r}$ is the IRS transmit correlation
matrix;
$[{{\mathbf{z}}}_{r,k}]_{n}\sim\mathcal{C}\mathcal{N}\left(\mathbf{0},l^{\text{UI}}_{k,n}\right)$
and
$[{\mathbf{F}}]_{m,n}\sim\mathcal{C}\mathcal{N}\left(\mathbf{0},l^{\text{IB}}_{m,n}\right)$
denote the NLoS Rayleigh fading components with $l^{\text{UI}}_{k,n}$ and
$l^{\text{IB}}_{m,n}$ denoting the corresponding path losses, while
${\bar{\mathbf{z}}}_{r,k}$ and $\bar{\mathbf{F}}$ represent the LoS
components.

More specifically, the LoS components of the user-BS, user-IRS and IRS-BS
channels, i.e., ${\bar{\mathbf{z}}}_{d,k}$, ${\bar{\mathbf{z}}}_{r,k}$ and
$\bar{\mathbf{F}}$, are respectively modeled as:

$[{\bar{\mathbf{z}}}_{d,k}]_{m}=\sqrt{l^{\text{UB}}_{k,m}}{f}_{L}({\vartheta}_{k}^{\text{UB}}),\;\;[{\bar{\mathbf{z}}}_{r,k}]_{n}=\sqrt{l^{\text{UI}}_{k,n}}{f}_{P}(\varphi^{\text{UI}}_{k},\psi^{\text{UI}}_{k}),\;\;[{\bar{\mathbf{F}}}]_{m,n}=\sqrt{l^{\text{IB}}_{m,n}}{f}_{L}(\vartheta^{\text{IB}}){f}_{P}(\varphi^{\text{IB}},\psi^{\text{IB}}),$
(3)

where ${f}_{P}(\varphi,\psi)=e^{{{j2\pi
d_{\text{I}}\big{(}(n-\lfloor{n}/{N_{x}}\rfloor
N_{x})\sin({\varphi})\cos({\psi})+\lfloor{n}/{N_{x}}\rfloor\sin({\varphi})\sin({\psi})\big{)}}/{\lambda}}}$,
${f}_{L}(\vartheta)=e^{{{j2\pi(m-1)d_{\text{B}}\sin({\vartheta})}/{\lambda}}}$,
$\lambda$ is the wavelength, $d_{B}$ denotes the antenna spacing at the BS,
$d_{I}$ denotes the reflecting element spacing at the IRS;
${\vartheta}_{k}^{\text{UB}}$ and ${\vartheta}^{\text{IB}}$ denote the angles-
of-arrival (AoAs) from $U_{k}$ to the BS and that from the IRS to the BS,
respectively; ${\psi}^{\text{IB}}$ and ${\varphi}^{\text{IB}}$ are the azimuth
and elevation angles-of-departure (AoDs) from the IRS to the BS, respectively;
${\psi}^{\text{UI}}_{k}$ and ${\varphi}^{\text{UI}}_{k}$ represent the azimuth
and elevation AoAs from $U_{k}$ to the IRS, respectively. Furthermore,
according to [24, 25], the BS receive correlation matrix can be modeled as

${\bm{\Phi}}_{d}(i,j)=\left\\{\begin{array}[]{ll}(r_{d})^{j-i},&\text{ if
}i\leq j,\\\ {\bm{\Phi}}_{d}^{*}(j,i),&\text{ otherwise },\end{array}\right.$
(4)

where $r_{d}\in[0,1]$ is the correlation coefficient. Besides,
${\bm{\Phi}}_{r}$ (${\bm{\Phi}}_{r,k}$) can be modeled as
${\bm{\Phi}}_{r}={\bm{\Phi}}_{r}^{h}\otimes{\bm{\Phi}}_{r}^{v}$
(${\bm{\Phi}}_{r,k}={\bm{\Phi}}_{r,k}^{h}\otimes{\bm{\Phi}}_{r,k}^{v}$)
according to [26], where ${\bm{\Phi}}_{r}^{h}$ (${\bm{\Phi}}_{r,k}^{h}$) and
${\bm{\Phi}}_{r}^{v}$ (${\bm{\Phi}}_{r,k}^{v}$) denote the spatial correlation
matrices of the horizontal and vertical domains, respectively, and they are
similarly defined as in (4) with $r_{r}$ ($r_{r,k}$) denoting the correlation
coefficient.

In the considered IRS-aided uplink multiuser SIMO system, each reflecting
element at the IRS is able to induce an independent phase-shift change of the
incident signal sent by the users. Thus, by properly adjusting the IRS
reflection coefficients at different reflecting elements, a preferable
wireless signal propagation environment between the transmitters and receiver
can be created to enhance the communication performance. Let $\theta_{n,i}$
denote the reflection coefficient of the $n$-th reflecting element in time
slot $i$, which satisfies $|\theta_{n,i}|=1$ if the $n$-th reflecting element
is on and $|\theta_{n,i}|=0$ otherwise. It is shown that the phase of the
incident signal on the $n$-th reflecting element can only be adjusted when the
reflecting element is switched on. The received signal at the BS in time slot
$i,i\in[1,T]$, is given by

$\mathbf{y}^{(i)}=\sum_{k=1}^{K}\mathbf{h}_{d,k}\sqrt{p_{k}}s_{k}^{(i)}+\sum_{k=1}^{K}\sum_{n=1}^{N}\theta_{n,i}[\mathbf{h}_{r,k}]_{n}\mathbf{g}_{n}\sqrt{p_{k}}s_{k}^{(i)}+\mathbf{n}^{(i)},$
(5)

where $s_{k}$ and $p_{k}$ denote the transmit signal and transmit power of
$U_{k}$, respectively,
$\mathbf{n}^{(i)}\sim\mathcal{C}\mathcal{N}\left(\mathbf{0},\sigma^{2}\mathbf{I}\right)$
denotes the additive white Gaussian noise (AWGN) at the BS in time slot $i$.
Define the reflected channel from $U_{k}$ to the BS via the IRS as

${\mathbf{H}}_{k}\triangleq\mathbf{G}\text{diag}({\mathbf{h}}_{r,k})\in{\mathbb{C}}^{M\times
N},$ (6)

and let $\mathbf{h}_{k,n}$ denote the $n$-th column of $\mathbf{H}_{k}$, then,
the received signal in (5) can be equivalently rewritten as111In this work, we
assume for simplicity that the transmit power of each user is set to the
maximum value $p$ to maximize the channel training SNR, i.e., $p_{k}=p,\forall
k\in[1,K]$.

$\mathbf{y}^{(i)}=\sum_{k=1}^{K}\big{(}\mathbf{h}_{d,k}+{\mathbf{H}}_{k}\bm{\theta}_{i})\sqrt{p}s_{k}^{(i)}+\mathbf{n}^{(i)},$
(7)

where $\bm{\theta_{i}}=[\theta_{1,i};\cdots;\theta_{N,i}]$. We can observe
from (7) that in order to recover the transmit signal
$\\{s^{(i)}_{k}\\}_{k=1}^{K}$ from the received signal $\mathbf{y}^{(i)}$, all
the direct channels $\\{\mathbf{h}_{d,k}\\}_{k=1}^{K}$ and reflected channels
$\\{\mathbf{H}_{k}\\}_{k=1}^{K}$ need to be obtained.

Generally, each coherence block of length $T$ time slots is divided into the
channel estimation stage (consists of $\tau<T$ time slots) and data
transmission stage (consists of $T-\tau$ time slots). Let ${a}_{k,i}$ denote
the pilot symbol of $U_{k}$ in time slot $i$ which satisfies
$|a_{k,i}|\in\\{0,1\\},i\in[1,T]$. Since the IRS is not capable of
transmitting or receiving pilot symbols, the BS is expected to estimate
$KM+KMN$ unknown coefficients included in
$\\{\mathbf{h}_{d,k},\mathbf{H}_{k}\\}_{k=1}^{K}$ based on the overall
received signal
$\mathbf{Y}\triangleq[\mathbf{y}^{(1)},\cdots,\mathbf{y}^{(\tau)}]\in\mathbb{C}^{M\times\tau}$,
the known user pilot symbol $\\{a_{k,i}\\}_{k=1}^{K},i\in[1,\tau]$ and the IRS
reflection coefficients $\\{\theta_{n,i}\\}_{n=1}^{N},i\in[1,\tau]$. When the
numbers of BS antennas, IRS reflecting elements and/or users are large, a
great amount of channel training overhead is required to estimate such a large
number of channel coefficients, which will reduce the sum rate of the
considered IRS-aided communication system. Motivated by the 3PCE strategy in
[22], we propose a new channel estimation strategy to exploit the strong
correlation among $\mathbf{H}_{k}$’s for channel training overhead reduction
and in the meantime, alleviate the error propagation issue by using fewer
channel estimation phases. The robust beamforming design for data transmission
with given estimated channels and channel error distributions will be left to
our future work.

##  III Proposed 2PCE Strategy

In this section, a novel 2PCE strategy for the considered IRS-aided multiuser
SIMO system is proposed. The main idea is to alleviate the error propagation
issue by reducing the channel estimation phases, which is essentially achieved
by estimating the direct and reflected channels associated with each user
simultaneously. Specifically, the 3PCE strategy in [22] works as follows: 1)
in the first phase, all the direct channels are estimated by switching off all
the reflecting elements; 2) in the second phase, the IRS is switched on and
only one typical user is allowed to transmit non-zero pilot symbols, then the
reflected channel associated with this typical user is estimated based on its
estimated direct channel in the first phase; 3) in the last phase, due to the
fact that the reflected channels associated with the other users are scaled
versions of the reflected channel associated with the typical user, the
scaling factors (rather than the whole reflected channels) are estimated based
on the estimated reflected channel associated with the typical user and the
direct channels. We can observe from the above 3PCE strategy that the
estimated direct channels obtained in the first phase are required in the next
two phases, which means that the estimation errors of the direct channels will
affect the estimation accuracy of all the reflected channels, especially when
the channel estimation in the first phase is not very accurate.222Note that
similar error propagation issues also exist in other channel estimation
strategies, e.g., [27, 12, 13], etc. Besides, estimation errors of the
reflected channel associated with the typical user in the second phase will
also deteriorate the estimation performance of the reflected channels
associated with other users. Therefore, the performance of the state-of-the-
art 3PCE strategy in [22] is limited by this error propagations issue and it
can be seen that the channel estimation performance of the second and third
phases are critically dependent on those in the first and second phases.

Figure 2: Proposed 2PCE strategy.

In the following, we present the proposed 2PCE strategy to alleviate the
abovementioned error propagation issue. As shown in Fig. 2, the proposed
strategy is divided into two phases:

  * •

In Phase I which consists of $\tau_{1}$ time slots, the direct and reflected
channels associated with a typical user (denoted by $U_{1}$ for convenience),
i.e., $\mathbf{h}_{d,1}$ and $\mathbf{H}_{1}$, are estimated simultaneously;

  * •

In Phase II which consists of $\tau_{2}$ time slots, the direct and reflected
channels associated with $U_{k}$, $k\in[2,K]$, are estimated simultaneously
based on the estimated reflected channel associated with the typical user by
exploiting the correlation among $\mathbf{H}_{k}$’s and properly designing the
pilot symbols of the users and the reflection patterns.

Note that in the proposed 2PCE strategy, the reflected channel associated with
the typical user is estimated in the first phase and it does not depend on the
estimation accuracy of the direct channels, therefore the estimation accuracy
of $\mathbf{H}_{1}$ can be improved without increasing the channel training
overhead. Moreover, the direct and reflected channels associated with each
user are estimated jointly, hence the estimation errors of the direct channels
will not propagate to the estimation process of the reflected channels.

For a better illustration of the overall channel estimation process of the
2PCE strategy and the required channel training overhead (i.e., the specific
numbers of $\tau_{1}$ and $\tau_{2}$), we consider the ideal case without
receive noise at the BS in the rest of this section. How to implement our
proposed 2PCE framework in the practical case with noise and the detailed
analysis of the asymptotic MSE for estimating the channels will be presented
in the next section.

###  III-A Phase I: Direct and Reflected Channels Estimation for Typical User

In Phase I, all the reflecting elements are switched on and only the typical
user $U_{1}$ transmits non-zero pilots to the BS, such that the direct and
reflected channels associated with $U_{1}$ can be estimated. Accordingly, the
received signal at the BS in time slot $i$ is given by

$\mathbf{y}^{(i)}={\mathbf{H}}_{1}\bm{\theta}_{i}\sqrt{p}a_{1,i}+\mathbf{h}_{d,1}\sqrt{p}a_{1,i},\;\;i\in[1,\tau_{1}],$
(8)

which can be viewed as a linear system of $M$ equations. Since there are
$M(N+1)$ unknown channel coefficients included in ${\mathbf{h}}_{d,1}$ and
${\mathbf{H}}_{1}$, at least $M(N+1)$ linear equations of these unknowns are
required to estimate $\\{\mathbf{h}_{d,1},\mathbf{H}_{1}\\}$, thus at least
$N+1$ time slots are needed in this phase. By defining ${\tau}_{1}=N+1$ and
$a_{k,i}=1,i\in[1,\tau_{1}]$, the overall received signal at the BS in Phase I
can be written as

${\mathbf{Y}}^{\text{I}}\triangleq[{\mathbf{y}}^{(1)},\cdots,{\mathbf{y}}^{(N+1)}]=\sqrt{p}{\mathbf{X}}^{\text{I}}{\mathbf{V}}^{\text{I}},$
(9)

where
${\mathbf{X}}^{\text{I}}\triangleq[{\mathbf{h}}_{d,1},{\mathbf{H}}_{1}]\in{\mathbb{C}}^{M\times(N+1)}$
denotes the composite channel matrix associated with $U_{1}$, and
${\mathbf{V}}^{\text{I}}\in{\mathbb{C}}^{(N+1)\times(N+1)}$ represents the
reflection pattern matrix the $i$-th column of which is given by
$[1;\bm{\theta}_{i}]$. To minimize the variance of the channel estimation
error, $\mathbf{V}^{\text{I}}$ should be an orthogonal matrix and each of its
entry should satisfy the uni-modulus constraint [14]. In this paper, similar
to [14, 15, 7], we set $\mathbf{V}^{\text{I}}$ as an $(N+1)$-dimension DFT
matrix which satisfies
$({\mathbf{V}}^{\text{I}})^{-1}=\frac{({\mathbf{V}}^{\text{I}})^{H}}{N+1}$.
Hence, the estimated composite channel can be expressed as

$\hat{{\mathbf{X}}^{\text{I}}}\triangleq[\hat{{\mathbf{h}}}_{d,1},\hat{{\mathbf{H}}}_{1}]=\frac{{\mathbf{Y}}^{\text{I}}({\mathbf{V}}^{\text{I}})^{H}}{\sqrt{p}(N+1)}.$
(10)

###  III-B Phase II: Direct and Reflected Channels Estimation for Other Users

In Phase II, the other users, i.e., $\\{U_{k}\\}_{k=2}^{K}$, transmit non-zero
pilot symbols to the BS to estimate the direct and reflected channels, i.e.,
$\\{\mathbf{h}_{d,k},\mathbf{H}_{k}\\}_{k=2}^{K}$. The received signal at the
BS in time slot $i$ can be written as

$\mathbf{y}^{(i)}=\sum_{k=2}^{K}{\mathbf{H}}_{k}\bm{\theta}_{i}\sqrt{p}a_{k,i}+\sum_{k=2}^{K}\mathbf{h}_{d,k}\sqrt{p}a_{k,i},\;\;i\in[\tau_{1}+1,\tau_{1}+\tau_{2}].$
(11)

Intuitively, if the correlations among $\mathbf{H}_{k}$’s is not exploited, at
least $(K-1)(N+1)$ time slots are required to estimate the $(K-1)(MN+M)$
channel coefficients in $\\{\mathbf{h}_{d,k},\mathbf{H}_{k}\\}_{k=2}^{K}$.
When the numbers of BS antennas, reflecting elements and users are large, the
channel training overhead required in Phase II would be overwhelming, which
leads to low user transmission rate since the time left for data transmission
is limited in this case. To overcome this difficulty, we adopt the strategy in
[22] and propose to express the reflected channel associated with $U_{k}$ as a
scaled version of $\mathbf{H}_{1}$, which is shown as follows:

$\mathbf{H}_{k}=\mathbf{H}_{1}\text{diag}(\bm{\mu}_{k}),$ (12)

where $\bm{\mu}_{k}\in\mathbb{C}^{N\times 1}$ represents the scaling vector of
$\mathbf{H}_{k}$ and the $n$-th entry of $\bm{\mu}_{k}$, i.e.,
$\mu_{k,n}={[\mathbf{h}_{r,k}]_{n}}/{[\mathbf{h}_{r,1}]_{n}}$, represents the
scaling factor associated with the $n$-th reflecting element. Utilizing the
relationship in (12), $\\{\mathbf{H}_{k}\\}_{k=2}^{K}$ can be estimated based
on the knowledge of $\\{\bm{\mu}_{k}\\}_{k=2}^{K}$ and $\mathbf{H}_{1}$. Thus,
we can estimate
$\big{\\{}\\{\mathbf{h}_{d,k}\\}_{k=1}^{K},\mathbf{H}_{1},\\{\bm{\mu}_{k}\\}_{k=2}^{K}\big{\\}}$
instead of
$\big{\\{}\\{\mathbf{h}_{d,k}\\}_{k=1}^{K},\\{\mathbf{H}_{k}\\}_{k=1}^{K}\big{\\}}$,
and the number of unknown channel coefficients to be estimated can be
significantly reduced. As a result, the received signal in (11) can be
equivalently rewritten as

${\mathbf{y}}^{(i)}=\sqrt{p}\sum_{k=2}^{K}a_{k,i}({\mathbf{h}}_{d,k}+{\mathbf{H}}_{1}\text{diag}({\bm{\theta}}_{i}){\bm{\mu}}_{k}),\;\;i\in[{\tau}_{1}+1,{\tau}_{1}+{\tau}_{2}].$
(13)

Let us define
${\mathbf{x}}^{\text{II}}\triangleq[{\bm{\mu}}_{2};\cdots;{\bm{\mu}}_{K};{\mathbf{h}}_{d,2};\cdots;{\mathbf{h}}_{d,K}]$,
${\mathbf{y}}^{\text{II}}\triangleq[{\mathbf{y}}^{(\tau_{1}+1)};\cdots;{\mathbf{y}}^{(\tau_{1}+\tau_{2})}]$
and

${\mathbf{W}}^{\text{II}}\triangleq\left[\begin{array}[]{cccccc}a_{2,\tau_{1}+1}{\mathbf{H}}_{1}\text{diag}({\bm{\theta}}_{\tau_{1}+1})&\cdots&a_{K,\tau_{1}+1}{\mathbf{H}}_{1}\text{diag}({\bm{\theta}}_{\tau_{1}+1})&a_{2,\tau_{1}+1}{{\mathbf{I}}_{M}}&\cdots&a_{K,\tau_{1}+1}{{\mathbf{I}}_{M}}\\\
\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
a_{2,\tau_{1}+\tau_{2}}{\mathbf{H}}_{1}\text{diag}({\bm{\theta}}_{\tau_{1}+\tau_{2}})&\cdots&a_{K,\tau_{1}+\tau_{2}}{\mathbf{H}}_{1}\text{diag}({\bm{\theta}}_{\tau_{1}+\tau_{2}})&a_{2,\tau_{1}+\tau_{2}}{{\mathbf{I}}_{M}}&\cdots&a_{K,\tau_{1}+\tau_{2}}{{\mathbf{I}}_{M}}\end{array}\right],$
(14)

then the overall received signal at the BS in Phase II is given by

${\mathbf{y}}^{\text{II}}=\sqrt{p}{\mathbf{W}}^{\text{II}}{\mathbf{x}}^{\text{II}},$
(15)

which is a linear system with $\tau_{2}M$ equations and $(K-1)(M+N)$ unknowns.

To show how to design the pilot symbols $a_{k,i}$’s and the reflection
patterns $\bm{\theta}_{i}$’s such that (14) can be solved efficiently, we
resort to the following lemma which indicates the necessary and sufficient
condition for a linear system of equations to have a unique solution.

######  Lemma 1.

[28] For any linear system of equations $\mathbf{Wx}=\mathbf{y}$ with
$\mathbf{W}\in{\mathbb{C}}^{d_{1}\times
d_{2}},\mathbf{x}\in{\mathbb{C}}^{d_{2}\times 1}$ and
${\mathbf{y}}\in{\mathbb{C}}^{d_{1}\times 1}$, there exists a unique solution
if and only if the matrix $\mathbf{W}$ satisfies: 1) $d_{1}\geq d_{2}$, 2)
$\text{rank}({\mathbf{W}})=d_{2}$.

Based on Lemma 1, we observe that $a_{k,i}$’s and $\bm{\theta}_{i}$’s should
be carefully designed to satisfy
$\text{rank}({\mathbf{W}}^{\text{II}})=(K-1)(M+N)$ and
$\tau_{2}M\geq(K-1)(M+N)$ such that $\mathbf{x}^{\text{II}}$ can be uniquely
determined by (15). Let us consider the following two cases for details.

####  III-B1 $M\geq N$

In this case, all the users except $U_{1}$ are allowed to transmit non-zero
pilots simultaneously, and orthogonal pilot sequences are employed to
eliminate the inter-user interference. In order to obtain the unique solution
of (15), i.e., $\\{\bm{\mu}_{k},\mathbf{h}_{d,k}\\}_{k=2}^{K}$, we set
$\tau_{2}=2(K-1)$ and design the following pilot sequence
${\mathbf{A}}^{\text{II}}$ and reflection patterns
$\\{{\bm{\theta}}_{i}\\}_{i=\tau_{1}+1}^{\tau_{1}+\tau_{2}}$ to meet the
requirements in Lemma 1:

${\mathbf{A}}^{\text{II}}\triangleq[{\mathbf{a}}_{2}^{\text{II}},\cdots,{\mathbf{a}}_{K}^{\text{II}}]=[\mathbf{D}_{K-1};\mathbf{D}_{K-1}],\;\;{\bm{\theta}}_{i}=\left\\{\begin{array}[]{ll}\;\mathbf{1}_{N\times
1},&\text{ if }i\in[\tau_{1}+1,\tau_{1}+K-1],\\\ \;-\mathbf{1}_{N\times
1},&\text{ otherwise,}\end{array}\right.$ (16)

where
$\mathbf{a}^{\text{II}}_{k}=[a_{k,\tau_{1}};\cdots;a_{k,\tau_{1}+\tau_{2}}]$
and $\mathbf{D}_{K-1}$ denotes the $(K-1)$-dimension DFT matrix. According to
(16), $\mathbf{x}^{\text{II}}$ can be obtained by multiplying
$\mathbf{y}^{\text{II}}$ on the right hand side (RHS) with the pseudo-inverse
of $\mathbf{W}^{\text{II}}$. However, the direct calculation of
$(\mathbf{W}^{\text{II}})^{\dagger}$ can be very computationally intensive due
to the large dimension of $\mathbf{W}^{\text{II}}$. To address this problem,
we first design an auxiliary matrix

${\bm{\Sigma}}_{k}\triangleq\left[\begin{array}[]{cccccc}a_{k,\tau_{1}+1}^{*}{\mathbf{I}}_{M}&\cdots&a_{k,\tau_{1}+K-1}^{*}{\mathbf{I}}_{M}&{\mathbf{0}}_{M\times
M}&\cdots&{\mathbf{0}}_{M\times M}\\\ {\mathbf{0}}_{M\times
M}&\cdots&{\mathbf{0}}_{M\times
M}&a_{k,\tau_{1}+K}^{*}{\mathbf{I}}_{M}&\cdots&a_{k,\tau_{1}+2(K-1)}^{*}{\mathbf{I}}_{M}\end{array}\right]\in{\mathbb{C}}^{2M\times(N+M)(K-1)},$
(17)

then by left-multiplying $\mathbf{y}^{\text{II}}$ with
${\bm{\Sigma}}_{k},k\in[2,K]$, (15) can be decomposed into $K-1$ smaller
linear systems of equations, which are shown as follows:

${\bm{\Sigma}}_{k}{\mathbf{y}}^{\text{II}}=\sqrt{p}(K-1){\mathbf{V}}^{\text{II}}[{\bm{\mu}}_{k};{\mathbf{h}}_{d,k}],\;\;k\in[2,K],$
(18)

where
${\mathbf{V}}^{\text{II}}\triangleq[{\mathbf{H}}_{1},{{\mathbf{I}}_{M}};-{\mathbf{H}}_{1},{{\mathbf{I}}_{M}}]$.
As can be observed from the channel model in (2) and (6), the columns of
$\mathbf{H}_{1}$ are linearly independent with probability one, thus the
matrix ${\mathbf{V}}^{\text{II}}$ is of full column rank. Besides, the number
of rows in ${\mathbf{V}}^{\text{II}}$ is larger than its number of columns
when $M\geq N$. Therefore, for each $k\in[2,K]$, there exists a unique
solution to (18), and $\\{{\bm{\mu}}_{k},{\mathbf{h}}_{d,k}\\}_{k=2}^{K}$ can
be perfectly estimated as

$[{\hat{\bm{\mu}}}_{k};{\hat{\mathbf{h}}}_{d,k}]=\frac{1}{\sqrt{p}(K-1)}({\mathbf{V}}^{\text{II}})^{\dagger}{\bm{\Sigma}}_{k}{\mathbf{y}}^{\text{II}},\;\;k\in[2,K].$
(19)

####  III-B2 $M\textless N$

In this case, in order for (15) to have a unique solution, at least
$\lceil\frac{(K-1)(M+N)}{M}\rceil=K-1+\lceil\frac{(K-1)N}{M}\rceil$ time slots
are needed to meet the first requirement in Lemma 1. However, compared with
the case of $M\geq N$, it is much more challenging in this case to design a
matrix $\mathbf{W}^{\text{II}}$ that satisfies the full-rank requirement when
all the users are active. For ease of pilot sequence and reflection pattern
design, we only allow one or a few users to transmit non-zero pilot symbols in
each time slot, as inspired by [22]. Our main idea to achieve the minimum
channel training overhead, i.e., $\tau_{2}=K-1+\lceil\frac{(K-1)N}{M}\rceil$,
is to allow the users to share some of the $\tau_{2}$ time slots with a
selected number of reflecting elements switched on.

More specifically, the proposed strategy in this case can be further divided
into two subphases, i.e., Phase II-A and Phase II-B. Let us define
$\gamma\triangleq\lfloor\frac{M+N}{M}\rfloor$ and $\delta\triangleq
N-(\gamma-1)M$, then Phase II-A (which consists of $(K-1)\gamma$ time slots)
is divided into $K-1$ separate non-overlapping time intervals each consists of
$\gamma$ time slots and each user is allocated a separate time interval. In
this subphase, each user transmits an all-one pilot sequence with some pre-
selected reflecting elements at the IRS switched on, such that the
corresponding direct channel and $M(\gamma-1)$ scaling factors can be
estimated. In Phase II-B which consists of $\lceil\frac{(K-1)\delta}{M}\rceil$
time slots, more than one users are activated and the rest reflecting elements
that are not selected in Phase II-A are switched on. Thereby, the rest
$(K-1)\delta$ scaling factors associated with each user can be obtained based
on the estimated $\mathbf{H}_{1}$ obtained in Phase I as well as the estimated
direct channels and $M(\gamma-1)(K-1)$ scaling factors obtained in Phase II-A.

Let $\chi_{k}^{\text{A}}$ and $\chi_{k}^{\text{B}}$ denote the index sets of
the scaling factors which are estimated in Phase II-A and Phase II-B,
respectively, we let

$\chi_{k}^{\text{B}}=\big{\\{}\big{(}(k-2)\delta+1\big{)}\;\%\;N,\cdots,\big{(}(k-1)\delta\big{)}\;\%\;N\big{\\}},\;\chi_{k}^{\text{A}}=\mathcal{W}-\chi_{k}^{\text{B}}.$
(20)

It can be seen that $\chi_{k}^{\text{A}}$ and $\chi_{k}^{\text{B}}$ contain
non-overlapping elements. In Phase II-A, only $U_{k}$ is permitted to transmit
pilot symbols in time slots $\xi_{k}+1\leq i\leq\xi_{k}+\gamma$ with
$\xi_{k}=\tau_{1}+(k-2)\gamma$, and the corresponding pilot sequence and
reflection patterns are set as

$a_{k,i}=\left\\{\begin{array}[]{ll}1,&\text{ if
}\lceil\frac{i-\tau_{1}}{\gamma}\rceil+1=k,\\\ 0,&\text{
otherwise,}\end{array}\right.\;\;{{\theta}}_{n,i}=\left\\{\begin{array}[]{ll}1,&\text{
if }n\in\mathcal{E}_{i}\text{ and }i\in[\xi_{k}+1,\xi_{k}+\gamma-1],\\\
-1,&\text{ if }n\in\chi_{k}^{\text{A}}\text{ and }i=\xi_{k}+\gamma,\\\
0,&\text{ otherwise,}\end{array}\right.$ (21)

where
$\mathcal{E}_{i}\triangleq\big{\\{}\chi_{k}^{\text{A}}(\lfloor\frac{i-\xi_{k}}{\gamma}\rfloor
M+1),\cdots,\chi_{k}^{\text{A}}(\lfloor\frac{i-\xi_{k}}{\gamma}\rfloor
M+M)\big{\\}}$. By defining

${\mathbf{V}}^{\text{II}}_{k}=\left[\begin{array}[]{cc}{\mathbf{H}}_{1}^{k}\text{diag}(\bar{{\bm{\theta}}}_{\xi_{k}+1})&{{\mathbf{I}}_{M}}\\\
\vdots&\vdots\\\
{\mathbf{H}}_{1}^{k}\text{diag}(\bar{{\bm{\theta}}}_{\xi_{k}+\gamma})&{{\mathbf{I}}_{M}}\end{array}\right],\;{\mathbf{H}}_{1}^{k}=[{\mathbf{h}}_{1,\chi_{k}^{\text{A}}(1)},\cdots,{\mathbf{h}}_{1,\chi_{k}^{\text{A}}(N-\delta)}],\;\bar{{\bm{\theta}}}_{i}=[\theta_{\chi_{k}^{\text{A}}(1),i};\cdots;\theta_{\chi_{k}^{\text{A}}(N-\delta),i}],$
(22)

the received signal at the BS can be written as

${\mathbf{y}}^{\text{II}}_{k}\triangleq[{\mathbf{y}}^{\xi_{k}+1};\cdots;{\mathbf{y}}^{\xi_{k}+\gamma}]=\sqrt{p}{\mathbf{V}}^{\text{II}}_{k}[{{\mu}}_{k,\chi_{k}^{\text{A}}(1)},\cdots,{{\mu}}_{k,\chi_{k}^{\text{A}}(N-\delta)},{\mathbf{h}}_{d,k}^{T}]^{T},\;\;k\in[2,K].$
(23)

As a result, the scaling factors in $\chi_{k}^{\text{A}}$ and the direct
channel associated with $U_{k}$ can be estimated simultaneously as

$[{\hat{\mu}}_{k,\chi_{k}^{\text{A}}(1)},\cdots,{\hat{\mu}}_{k,\chi_{k}^{\text{A}}(N-\delta)},{\hat{\mathbf{h}}}_{d,k}^{T}]^{T}=\frac{1}{\sqrt{p}}({\mathbf{V}}^{\text{II}}_{k})^{\dagger}{\mathbf{y}}^{\text{II}}_{k}.$
(24)

In Phase II-B, the pilot sequences of the users and reflection patterns at the
IRS in time slot $i\geq\tau_{1}+(K-1)\gamma+1$ are designed as

$a_{k,i}=\left\\{\begin{array}[]{ll}1,&\text{ if }k\in\mathcal{R}_{i},\\\
0,&\text{
otherwise},\end{array}\right.\;\;\theta_{n,i}=\left\\{\begin{array}[]{ll}1,&\text{
if }n\in\mathcal{W}_{i}\\\ 0,&\text{ otherwise},\end{array}\right.$ (25)

where
$\mathcal{R}_{i}\triangleq\\{\lceil\frac{(i-\tau_{1}-(K-1)\gamma-1)M+1}{\delta}\rceil+1,\cdots,\lceil\frac{(i-\tau_{1}-(K-1)\gamma-1)M+M_{i}}{\delta}\rceil+1\\}$,
$\mathcal{W}_{i}\triangleq\\{((i-\tau_{1}-(K-1)\gamma-1)M+m)\%N,\forall
m\in[1,M_{i}]\\}$ and $M_{i}=\min(M,(K-1)\delta-(i-\tau_{1}-(K-1)\gamma-1)M)$.
Then, the received signal can be written as

${\mathbf{y}}^{\left(i\right)}=\sqrt{p}\sum\limits_{k\in\mathcal{R}_{i}}\mathbf{h}_{d,k}+\sqrt{p}\sum\limits_{k\in\mathcal{R}_{i}}\sum\limits_{n\in\mathcal{W}_{i}}\mu_{k,n}\mathbf{h}_{1,n},\;\;i\in[\tau_{1}+(K-1)\gamma+1,\tau],$
(26)

which includes $|\mathcal{R}_{i}||\mathcal{W}_{i}|$ scaling factors, i.e.,
$\mu_{k,n},k\in\mathcal{R}_{i},n\in\mathcal{W}_{i}$. Note that with the
designed pilot sequences and reflection patterns in (21), the scaling factors
$\mu_{\mathcal{R}_{i}(j_{1}),\mathcal{W}_{i}(j_{2})},j_{1},j_{2}\in[1,M_{i}],j_{1}\neq
j_{2}$ and the direct channels $\mathbf{h}_{d,k},k\in\mathcal{R}_{i}$ have
been estimated in Phase II-A according to (24), hence the rest scaling factors
$\mu_{\mathcal{R}_{i}(j),\mathcal{W}_{i}(j)},j\in[1,M_{i}]$ can be easily
estimated as follows:

$[{\hat{\mu}}_{\mathcal{R}_{i}(1),\mathcal{W}_{i}(1)};\cdots;{\hat{\mu}}_{\mathcal{R}_{i}(M_{i}),\mathcal{W}_{i}(M_{i})}]=[{\hat{\mathbf{h}}}_{1,\mathcal{W}_{i}(1)},\cdots,{\hat{\mathbf{h}}}_{1,\mathcal{W}_{i}(M_{i})}]^{\dagger}\frac{{\bar{\mathbf{y}}}^{(i)}}{\sqrt{p}},\;\;i\in[\tau_{1}+(K-1)\gamma+1,\tau],$
(27)

where
${\bar{\mathbf{y}}}^{(i)}={\mathbf{y}}^{\left(i\right)}-\sum_{k\in\mathcal{R}_{i}}\sqrt{p}{\hat{\mathbf{h}}}_{d,k}-\sum_{k\in\mathcal{R}_{i}}\sum_{n\in\mathcal{W}_{i}\cap\chi_{k}^{\text{A}}}\sqrt{p}\hat{\mu}_{k,n}{\hat{\mathbf{h}}}_{1,n}$.

To sum up, according to (19) and (24), the direct channel and
$\min(N,(\gamma-1)M)$ scaling factors associated with $U_{k},k\neq 1$, are
estimated simultaneously, hence the estimation accuracy of the scaling factors
will not be deteriorated by the estimation error of the corresponding direct
channel.

###  III-C Overall Channel Training Overhead

According to the above two subsections, the overall channel training overhead
required by the proposed 2PCE strategy is

$\tau=\tau_{1}+\tau_{2}=\left\\{\begin{array}[]{ll}N+2K-1,&\text{ if }M\geq
N,\\\ N+K+\lceil\frac{(K-1)N}{M}\rceil,&\text{ otherwise, }\end{array}\right.$
(28)

which is the same with that of the 3PCE strategy in [22], i.e,
$K+N+\max(K-1,\lceil\frac{(K-1)N}{M}\rceil)$. Note that both the proposed 2PCE
strategy and the 3PCE strategy in [22] are feasible solutions for achieving
the minimum required channel training overhead in the considered IRS-aided
multiuser SIMO system, the difference is that the proposed 2PCE strategy
contains less estimation phases and thus the error propagation can be well-
controlled.

##  IV  Asymptotic MSE analysis

In the previous section, we have shown how to perfectly estimate all the
channels in detail for the ideal case without receive noise at the BS. In this
section, we consider the practical case with noise and present an LS estimator
for the proposed 2PCE strategy.333Note that when the estimation error of the
reflected channel associated with the typical user is not neglected, the
derivation of the linear minimum MSE (LMMSE) estimator can be very complex,
hence for simplicity we only consider the LS estimation in this work. Further
investigation into more complicated estimators are left for future work.
Moreover, we derive the asymptotic MSE of the estimated $\mathbf{H}_{1}$,
$\\{\mathbf{h}_{d,k}\\}_{k=1}^{K}$ and $\\{\bm{\mu}_{k}\\}_{k=2}^{K}$ when $M$
becomes large to investigate how the error propagation issue is alleviated by
the proposed 2PCE strategy. Besides, performance comparison between the
proposed 2PCE strategy and the 3PCE strategy [22] in terms of the asymptotic
MSE is provided to show the advantages of the proposed strategy. In the
sequel, for ease of analysis, we assume that both the IRS-BS and user-IRS
channels follow Rayleigh distribution, i.e.,
$\mathbf{g}_{n}\sim\mathcal{CN}(\mathbf{0},l^{\textrm{IB}}\mathbf{I})$ and
$\mathbf{h}_{r,k}\sim\mathcal{CN}(\mathbf{0},l^{\textrm{UI}}_{k}\mathbf{I})$,
where $l^{\textrm{IB}}$ ($l^{\textrm{IU}}_{k}$) represents the path loss of
the link from the IRS to the BS (from $U_{k}$ to the IRS). In addition, some
necessary propositions and lemmas used throughout this section are listed in
Appendix A.

###  IV-A MSE Analysis in Phase I

With receive noise at the BS, the received signal at the BS given in (9) can
be rewritten as

${\mathbf{Y}}^{\text{I}}=\sqrt{p}{\mathbf{X}}^{\text{I}}{\mathbf{V}}^{\text{I}}+{\mathbf{N}}^{\text{I}},$
(29)

where
${\mathbf{N}}^{\text{I}}=[{\mathbf{n}}^{(1)},\cdots,{\mathbf{n}}^{(N+1)}]\in{\mathbb{C}}^{M\times(N+1)}$.
By right-multiplying $\mathbf{Y}^{\text{I}}$ by the inverse of
${\mathbf{V}}^{\text{I}}$ (or equivalently,
$\frac{(\mathbf{V}^{\text{I}})^{H}}{N+1}$), we can obtain the LS estimation of
$\mathbf{h}_{d,1}$ and $\mathbf{H}_{1}$ as

$\hat{{\mathbf{X}}^{\text{I}}}\triangleq[\hat{\mathbf{h}}_{d,1},\hat{\mathbf{H}}_{1}]={\mathbf{X}}^{\text{I}}+\frac{{\mathbf{N}}^{\text{I}}({\mathbf{V}}^{\text{I}})^{H}}{\sqrt{p}(N+1)}.$
(30)

From (30), it can be observed that the estimation errors of the direct and
reflected channels associated with $U_{1}$ have complex Gaussian entries since
the channel noise is modeled as AWGN with i.i.d. entries and the adopted LS
estimator only involves linear operations to $\mathbf{Y}^{\text{I}}$. Besides,
we can see that the CSI error matrix can be expressed as

$\Delta{\mathbf{X}}^{\text{I}}\triangleq[\Delta\mathbf{h}_{d,1},\Delta\mathbf{H}_{1}]={\mathbf{X}}^{\text{I}}-\hat{{\mathbf{X}}^{\text{I}}}=\frac{{\mathbf{N}}^{\text{I}}({\mathbf{V}}^{\text{I}})^{H}}{\sqrt{p}(N+1)},$
(31)

which satisfies $\mathbb{E}\\{\Delta{\mathbf{X}}^{\text{I}}\\}={\mathbf{0}}$
and
$\mathbb{E}\\{\Delta{\mathbf{X}}^{\text{I}}(\Delta{\mathbf{X}}^{\text{I}})^{H}\\}=\frac{\sigma^{2}}{p}\mathbf{I}_{M}$.
As a result, the MSEs for estimating $\mathbf{h}_{d,1}$ and $\mathbf{H}_{1}$
are respectively given by

$\epsilon_{d,1,\text{2P}}\triangleq\mathbb{E}\big{\\{}\text{tr}(\Delta{\mathbf{h}}_{d,1}\Delta{\mathbf{h}}_{d,1}^{H})\big{\\}}=\frac{M\sigma^{2}}{p(N+1)},\;\;\epsilon_{r,1,\text{2P}}\triangleq\mathbb{E}\big{\\{}\text{tr}(\Delta{\mathbf{H}}_{1}\Delta{\mathbf{H}}_{1}^{H})\big{\\}}=\frac{MN\sigma^{2}}{p(N+1)}.$
(32)

From (32), we can see that $\epsilon_{d,1,\text{2P}}$ decreases with the
increasing of $N$ since the estimation of $\mathbf{h}_{d,1}$ is obtained based
on the received signals in a total number of $N+1$ time slots, hence the MSE
for estimating $\mathbf{h}_{d,1}$ in the 2PCE strategy is expected to be lower
than that in the 3PCE strategy which only uses one time slot. Besides, we can
observe that $\epsilon_{d,1,\text{2P}}$ and $\epsilon_{r,1,\text{2P}}$ are not
related, which means that the estimation of the reflected channel associated
with the typical user will not be affected by that of the corresponding direct
channel.

###  IV-B Asymptotic MSE Analysis in Phase II

In Phase II, the direct and reflected channels associated with the non-typical
users are estimated simultaneously to alleviate the negative effects caused by
error propagation. By replacing $\mathbf{H}_{1}$ with
$\hat{\mathbf{H}}_{1}+\Delta\mathbf{H}_{1}$, the effective received signal at
the BS given in (13) can be re-expressed as

${\mathbf{y}}^{(i)}=\sqrt{p}\sum\limits_{k=2}^{K}a_{k,i}({\mathbf{h}}_{d,k}+{\hat{\mathbf{H}}}_{1}\text{diag}({\bm{\theta}}_{i}){\bm{\mu}}_{k})+\sqrt{p}\sum\limits_{k=2}^{K}a_{k,i}\Delta{\mathbf{H}}_{1}\text{diag}({\bm{\theta}}_{i}){\bm{\mu}}_{k}+{\mathbf{n}}^{(i)},\;\;i\in[\tau_{1}+1,\tau_{1}+\tau_{2}].$
(33)

Then, the overall received signal at the BS with receive noise can be written
as

${\mathbf{y}}^{\text{II}}\triangleq[{\mathbf{y}}^{(\tau_{1})};\cdots;{\mathbf{y}}^{(\tau_{1}+\tau_{2})}]=\sqrt{p}{\hat{\mathbf{W}}}^{\text{II}}{\mathbf{x}}^{\text{II}}+\sqrt{p}{\mathbf{W}}_{\Delta{\mathbf{H}}_{1}}[{\bm{\mu}}_{2};\cdots;{\bm{\mu}}_{K}]+{\mathbf{n}}^{\text{II}},$
(34)

where
${\mathbf{n}}^{\text{II}}=[{\mathbf{n}}^{(\tau_{1}+1)};\cdots;{\mathbf{n}}^{(\tau_{1}+\tau_{2})}]$,

${\hat{\mathbf{W}}}^{\text{II}}=\left[\begin{array}[]{cccccc}a_{2,\tau_{1}+1}{\hat{\mathbf{H}}}_{1}\text{diag}({\bm{\theta}}_{\tau_{1}+1})&\cdots&a_{K,\tau_{1}+1}{\hat{\mathbf{H}}}_{1}\text{diag}({\bm{\theta}}_{\tau_{1}+1})&a_{2,\tau_{1}+1}{{\mathbf{I}}_{M}}&\cdots&a_{K,\tau_{1}+1}{{\mathbf{I}}_{M}}\\\
\vdots&\ddots&\vdots&\vdots&\ddots&\vdots\\\
a_{2,\tau_{1}+\tau_{2}}{\hat{\mathbf{H}}}_{1}\text{diag}({\bm{\theta}}_{\tau_{1}+\tau_{2}})&\cdots&a_{K,\tau_{1}+\tau_{2}}{\hat{\mathbf{H}}}_{1}\text{diag}({\bm{\theta}}_{\tau_{1}+\tau_{2}})&a_{2,\tau_{1}+\tau_{2}}{{\mathbf{I}}_{M}}&\cdots&a_{K,\tau_{1}+\tau_{2}}{{\mathbf{I}}_{M}}\end{array}\right],$
(35)
${\mathbf{W}}_{\Delta{\mathbf{H}}_{1}}=\left[\begin{array}[]{ccc}a_{2,\tau_{1}+1}\Delta{\mathbf{H}}_{1}\text{diag}({\bm{\theta}}_{\tau_{1}+1})&\cdots&a_{K,\tau_{1}+1}\Delta{\mathbf{H}}_{1}\text{diag}({\bm{\theta}}_{\tau_{1}+1})\\\
\vdots&\ddots&\vdots\\\
a_{2,\tau_{1}+\tau_{2}}\Delta{\mathbf{H}}_{1}\text{diag}({\bm{\theta}}_{\tau_{1}+\tau_{2}})&\cdots&a_{K,\tau_{1}+\tau_{2}}\Delta{\mathbf{H}}_{1}\text{diag}({\bm{\theta}}_{\tau_{1}+\tau_{2}})\end{array}\right].$
(36)

Similar to Section III-B, the asymptotic MSE analysis of the proposed 2PCE
strategy in Phase II are also divided in two cases: $M\geq N$ and $M\textless
N$.

####  IV-B1 $M\geq N$

In order to estimate the scaling vector $\bm{\mu}_{k}$ and the direct channel
${\mathbf{h}}_{d,k}$, we left-multiply ${\mathbf{y}}^{\text{II}}$ with the
auxiliary matrix ${\bm{\Sigma}}_{k}$ (given in (17)) and thereby obtain the
following equation:

${\bm{\Sigma}}_{k}{\mathbf{y}}^{\text{II}}=\sqrt{p}(K-1)({\hat{\mathbf{V}}}^{\text{II}}[{\bm{\mu}}_{k};{\mathbf{h}}_{d,k}]+\Delta{\mathbf{S}}{\bm{\mu}}_{k})+{\bm{\Sigma}}_{k}{\mathbf{n}}^{\text{II}},\;\;k\in[2,K],$
(37)

where
${\hat{\mathbf{V}}}^{\text{II}}=[{\hat{\mathbf{H}}}_{1},{{\mathbf{I}}_{M}};-{\hat{\mathbf{H}}}_{1},{{\mathbf{I}}_{M}}]$
and $\Delta{{\mathbf{S}}}=[\Delta{\mathbf{H}}_{1};-\Delta{\mathbf{H}}_{1}]$.
Therefore, the LS estimation of ${\bm{\mu}}_{k}$ and ${\mathbf{h}}_{d,k}$ can
be obtained as follows:

$[{\hat{\bm{\mu}}}_{k};{\hat{\mathbf{h}}}_{d,k}]=[{\bm{\mu}}_{k};{\mathbf{h}}_{d,k}]+\underbrace{\frac{1}{\sqrt{p}(K-1)}({\hat{\mathbf{V}}}^{\text{II}})^{\dagger}{\bm{\Sigma}}_{k}{\mathbf{n}}^{\text{II}}+({\hat{\mathbf{V}}}^{\text{II}})^{\dagger}\Delta{\mathbf{S}}{\bm{\mu}}_{k}}_{\Delta\mathbf{x}_{k}^{\text{II}}},\;\;k\in[2,K],$
(38)

where $\Delta\mathbf{x}_{k}^{\text{II}}$ denotes the composite estimation
error matrix in Phase II. It can be seen from (38) that the estimation
accuracy in this phase is closely related to two factors, i.e., the channel
noise $\mathbf{n}^{\text{II}}$ and the estimation error of $\mathbf{H}_{1}$.
Since each entry of $\mathbf{n}^{\text{II}}$ and $\Delta\mathbf{S}$ follows
independent complex normal distribution, the covariance matrix of
$\Delta{\mathbf{x}}^{\text{II}}_{k}$ can be decomposed into two parts, i.e.,

$\begin{array}[]{l}\mathbb{E}\big{\\{}\Delta{\mathbf{x}}^{\text{II}}_{k}(\Delta{\mathbf{x}}^{\text{II}}_{k})^{H}\big{\\}}=\mathbf{P}_{k,1}^{\text{a}}+\mathbf{P}_{k,2}^{\text{a}},\end{array}$
(39)

where
$\mathbf{P}_{k,1}^{\text{a}}\triangleq\frac{\mathbb{E}\\{({\hat{\mathbf{V}}}^{\text{II}})^{\dagger}{\bm{\Sigma}}_{k}{\mathbf{n}}^{\text{II}}({\mathbf{n}}^{\text{II}}{\bm{\Sigma}}_{k})^{H}(({\hat{\mathbf{V}}}^{\text{II}})^{H})^{\dagger}\\}}{p(K-1)^{2}}$
and
$\mathbf{P}_{k,2}^{\text{a}}\triangleq\mathbb{E}\big{\\{}({\hat{\mathbf{V}}}^{\text{II}})^{\dagger}\Delta{\mathbf{S}}{\bm{\mu}}_{k}{\bm{\mu}}_{k}^{H}\Delta{\mathbf{S}}^{H}(({\hat{\mathbf{V}}}^{\text{II}})^{\dagger})^{H}\big{\\}}$.
Note that the first part $\mathbf{P}_{k,1}^{\text{a}}$ in (39) is related to
both $\mathbf{n}^{\text{II}}$ and $\Delta\mathbf{H}_{1}$ (contained in
${\hat{\mathbf{V}}}^{\text{II}}$), while the second part
$\mathbf{P}_{k,2}^{\text{a}}$ is only related to $\Delta\mathbf{H}_{1}$
(contained in ${\hat{\mathbf{V}}}^{\text{II}}$ and $\Delta\mathbf{S}$). The
exact characterization of $\mathbf{P}_{k,1}^{\text{a}}$ and
$\mathbf{P}_{k,2}^{\text{a}}$ is very challenging since the pseudo-inverse
operation of the random matrix ${\hat{\mathbf{V}}}^{\text{II}}$ or
$({\hat{\mathbf{V}}}^{\text{II}})^{H}$ is involved.444Note that in [23], the
MSE performance of the LMMSE estimator is derived under the assumption that
the reflected channel associated with the typical user is perfectly known. In
this work, we focus on the more practical case, where no such assumption is
made.  In the following, we derive trackable forms for these two parts based
on random matrix theory and some necessary lemmas presented in Appendix A.

First, we focus on the derivation of $\mathbf{P}_{k,1}^{\text{a}}$. Since the
entries of $\mathbf{n}^{\text{II}}$ and
$(\hat{\mathbf{V}}_{\text{II}})^{\dagger}$ are independent with each other,
$\mathbf{P}_{k,1}^{\text{a}}$ can be equivalently transformed into
$\mathbf{P}_{k,1}^{\text{a}}=\frac{\sigma^{2}}{p(K-1)}\mathbb{E}\big{\\{}\big{(}({\hat{\mathbf{V}}}^{\text{II}})^{H}{\hat{\mathbf{V}}}^{\text{II}}\big{)}^{-1}\big{\\}}$.
Moreover, we observe that
$({\hat{\mathbf{V}}}^{\text{II}})^{H}{\hat{\mathbf{V}}}^{\text{II}}=2[\hat{{\mathbf{H}}}_{1}^{H}{\hat{\mathbf{H}}}_{1},\mathbf{0}_{N\times
M};\mathbf{0}_{M\times N},{\mathbf{I}}_{M}]$, which is a $2\times 2$ block
matrix, therefore the inverse of
$({\hat{\mathbf{V}}}^{\text{II}})^{H}{\hat{\mathbf{V}}}^{\text{II}}$ can be
obtained according to Lemma 2 and $\mathbf{P}_{k,1}^{\text{a}}$ can be further
simplified as follows:

$\displaystyle\mathbf{P}_{k,1}^{\text{a}}=\frac{\sigma^{2}}{2p(K-1)}\left[\begin{array}[]{cc}\mathbb{E}\big{\\{}(\hat{{\mathbf{H}}}_{1}^{H}{\hat{\mathbf{H}}}_{1})^{-1}\big{\\}}&{\mathbf{0}}_{N\times
M}\\\ {\mathbf{0}}_{M\times N}&{\mathbf{I}}_{M}\end{array}\right].$ (40)

As for $\mathbf{P}_{2}^{\text{a}}$, we first define
$\mathbf{M}_{1}\triangleq\Delta{\mathbf{S}}{\bm{\mu}}_{k}{\bm{\mu}}_{k}^{H}\Delta{\mathbf{S}}^{H}=[\mathbf{E}_{k},-\mathbf{E}_{k};-\mathbf{E}_{k},\mathbf{E}_{k}]$
where $\mathbf{E}_{k}=\Delta\mathbf{H}_{1}\bm{\mu}_{k}\bm{\mu}_{k}^{H}$
$\Delta\mathbf{H}_{1}^{H}$, then the pseudo-inverse of $\mathbf{M}_{1}$ can be
obtained by
$\frac{1}{4}[\mathbf{E}_{k}^{\dagger},-\mathbf{E}_{k}^{\dagger};-\mathbf{E}_{k}^{\dagger},\mathbf{E}_{k}^{\dagger}]$
according to Lemma 3. As a result, $\mathbf{P}_{k,2}^{\text{a}}$ can be
simplified as

$\displaystyle\mathbf{P}_{k,2}^{\text{a}}=\mathbb{E}\big{\\{}\big{(}({\hat{\mathbf{V}}}^{\text{II}})^{H}\mathbf{M}_{1}^{\dagger}{\hat{\mathbf{V}}}^{\text{II}}\big{)}^{\dagger}\big{\\}}=\left[\begin{array}[]{cc}\mathbb{E}\\{\mathbf{R}_{k}\\}&\mathbf{0}_{N\times
M}\\\ \mathbf{0}_{M\times N}&\mathbf{0}_{M\times M}\end{array}\right],$ (41)

where
$\mathbf{R}_{k}=\hat{\mathbf{H}}_{1}^{\dagger}\Delta\mathbf{H}_{1}\bm{\mu}_{k}\bm{\mu}_{k}^{H}\Delta\mathbf{H}_{1}^{H}(\hat{\mathbf{H}}_{1}^{H})^{\dagger}$.
Based on (40) and (41), the covariance matrix of
$\Delta{\mathbf{x}}^{\text{II}}_{k}$ in (39) can be rewritten as

$\mathbb{E}\big{\\{}\Delta{\mathbf{x}}^{\text{II}}_{k}(\Delta{\mathbf{x}}^{\text{II}}_{k})^{H}\big{\\}}=\frac{\sigma^{2}}{2p(K-1)}\left[\begin{array}[]{cc}\mathbb{E}\big{\\{}(\hat{{\mathbf{H}}}_{1}^{H}{\hat{\mathbf{H}}}_{1})^{-1}\big{\\}}&{\mathbf{0}}_{N\times
M}\\\ {\mathbf{0}}_{M\times
N}&{\mathbf{I}}_{M}\end{array}\right]+\left[\begin{array}[]{cc}\mathbb{E}\\{\mathbf{R}_{k}\\}&\mathbf{0}_{N\times
M}\\\ \mathbf{0}_{M\times N}&\mathbf{0}_{M\times M}\end{array}\right].$ (42)

Accordingly, the MSEs for estimating $\mathbf{h}_{d,k}$ and
$\bm{\mu}_{k},k\in[2,K]$, are given by

$\epsilon_{d,k,\text{2P}}^{\text{a}}\triangleq\mathbb{E}\big{\\{}\text{tr}(\Delta{\mathbf{h}}_{d,k}\Delta{\mathbf{h}}_{d,k}^{H})\big{\\}}=\frac{M\sigma^{2}}{2p(K-1)},$
(43a)
$\displaystyle\epsilon_{\bm{\mu},k,\text{2P}}^{\text{a}}\triangleq\mathbb{E}\big{\\{}\text{tr}(\Delta\bm{\mu}_{k}\Delta\bm{\mu}_{k}^{H})\big{\\}}=\frac{\sigma^{2}}{2p(K-1)}\mathbb{E}\big{\\{}\text{tr}\big{(}(\hat{{\mathbf{H}}}_{1}^{H}{\hat{\mathbf{H}}}_{1})^{-1}\big{)}\big{\\}}+\mathbb{E}\big{\\{}\text{tr}(\mathbf{R}_{k})\big{\\}}.$
(43b)

It can be seen from (43a) that $\epsilon_{d,k,\text{2P}}^{\text{a}}$ is only
affected by the channel noise in Phase II and it is inversely proportional to
$2(K-1)$ since the estimation of $\mathbf{h}_{d,k}$ is obtained based on the
received signals in a total number of $2(K-1)$ time slots. Besides,
$\epsilon_{\bm{\mu},k,\text{2P}}^{\text{a}}$ is further related to the CSI
error matrix of $\mathbf{H}_{1}$, which is still difficult to characterize
since it involves the inverse operation of
$\hat{{\mathbf{H}}}_{1}^{H}{\hat{\mathbf{H}}}_{1}$.

Let
$\mathbf{P}_{\bm{\mu}_{k},1}^{\text{a}}\triangleq\frac{\sigma^{2}}{2p(K-1)}\mathbb{E}\big{\\{}\text{tr}\big{(}(\hat{{\mathbf{H}}}_{1}^{H}{\hat{\mathbf{H}}}_{1})^{-1}\big{)}\big{\\}}$
and
$\mathbf{P}_{\bm{\mu}_{k},2}^{\text{a}}\triangleq\mathbb{E}\big{\\{}\text{tr}(\mathbf{R}_{k})\big{\\}}$
denote the two terms on the RHS of (43b), in the following, we focus on
further simplifying $\mathbf{P}_{\bm{\mu}_{k},1}^{\text{a}}$ and
$\mathbf{P}_{\bm{\mu}_{k},2}^{\text{a}}$ when $M$ becomes large to draw useful
insights and we have the following proposition.

######  Proposition 1.

$\mathbf{P}_{\bm{\mu}_{k},1}^{\text{a}}$ and
$\mathbf{P}_{\bm{\mu}_{k},2}^{\text{a}}$ can be approximated by
$\mathbf{P}_{\bm{\mu}_{k},1}^{\text{a}}\approx\frac{\sigma^{2}(N+1)N}{2(K-1)M\big{(}p(N+1)l^{\text{UI}}_{1}l^{\text{IB}}+\sigma^{2}\big{)}}$
and
$\mathbf{P}_{\bm{\mu}_{k},2}^{\text{a}}\approx\frac{M\sigma^{2}\mathbb{E}\\{\bm{\mu}_{k}^{H}\bm{\mu}_{k}\\}}{Nl^{\text{UI}}_{1}l^{\text{IB}}p(N+1)+M\sigma^{2}}$
when $M$ becomes asymptotically large.

###### Proof.

Please refer to Appendix B. ∎

According to Proposition 1, by substituting the asymptotic values of
$\mathbf{P}_{\bm{\mu}_{k},1}^{\text{a}}$ and
$\mathbf{P}_{\bm{\mu}_{k},2}^{\text{a}}$ into (43b), we can obtain the
asymptotic MSE for estimating $\bm{\mu}_{k}$ as follows:

$\displaystyle\epsilon_{\bm{\mu}_{k},\text{2P}}^{\text{a}}\approx\frac{\sigma^{2}(N+1)N}{2(K-1)M\big{(}p(N+1)l^{\text{UI}}_{1}l^{\text{IB}}+\sigma^{2}\big{)}}+\frac{M\sigma^{2}\mathbb{E}\\{\bm{\mu}_{k}^{H}\bm{\mu}_{k}\\}}{Nl^{\text{UI}}_{1}l^{\text{IB}}p(N+1)+M\sigma^{2}}.$
(44)

From (44), we can observe that $\epsilon_{\bm{\mu}_{k},\text{2P}}^{\text{a}}$
is inversely proportional to the transmit power $p$, the path losses of the
links from the typical user to the IRS and from the IRS to the BS, i.e.,
$l^{\text{UI}}_{1}$ and $l^{\text{IB}}$, while it is directly proportional to
$\mathbb{E}\\{\bm{\mu}_{k}^{H}\bm{\mu}_{k}\\}$. The first term on the RHS of
(44) will dominate the second term when $N$ becomes asymptotically large, and
$\epsilon_{\bm{\mu}_{k},\text{2P}}^{\text{a}}$ is directly proportional to $N$
in this case. Besides, since the first term on the RHS of (44) is affected by
the channel noise contained in the $2(K-1)$ received signals,
$\epsilon_{\bm{\mu}_{k},\text{2P}}^{\text{a}}$ decreases with the increasing
of $K$.

####  IV-B2 $M\textless N$

In this case, we assume for simplicity an orthogonal channel training
strategy, i.e., the users do not share the time slots and only one user is
allowed to transmit its pilot symbol to the BS in each time slot.555 Note that
in the general case when the users are allowed to share some of the time
slots, the asymptotic MSE analysis can be very complicated. In this paper, we
focus on this simplified case, and further investigation into the more general
case is left for future work.  Specifically, in time slot
$i\in[\xi^{\prime}_{k}+1,\xi^{\prime}_{k}+\gamma+1]$ with
$\xi^{\prime}_{k}=\tau_{1}+(\gamma+1)(k-1)$, only $U_{k}$ is allowed to
transmit pilot symbol $a_{k,i}=1$ and some selected reflecting elements are
switched on, then the direct channel and corresponding scaling factors
associated with $U_{k}$ can be estimated. The designed reflection patterns in
time slot $i\in[\xi^{\prime}_{k}+1,\xi^{\prime}_{k}+\gamma+1]$ are set as

$\theta_{n,i}=\left\\{\begin{array}[]{ll}1,&\text{ if
}i-\xi^{\prime}_{k}=\lceil\frac{n}{M}\rceil,\\\ -1,&\text{ if
}i-\xi^{\prime}_{k}=\gamma+1,\\\ 0,&\text{ otherwise. }\end{array}\right.$
(45)

Let us define
$\Delta{{\mathbf{S}}}_{k}=[\Delta{{\mathbf{H}}}_{1}\text{diag}({{\bm{\theta}}}_{\xi^{\prime}_{k}+1});\cdots;\Delta{{\mathbf{H}}}_{1}\text{diag}({{\bm{\theta}}}_{\xi^{\prime}_{k}+\gamma+1})]$,

${\mathbf{Q}}^{\text{II}}_{k}=\left[\begin{array}[]{cc}{{\mathbf{H}}}_{1}\text{diag}({{\bm{\theta}}}_{\xi^{\prime}_{k}+1})&{{\mathbf{I}}_{M}}\\\
\vdots&\vdots\\\
{{\mathbf{H}}}_{1}\text{diag}({{\bm{\theta}}}_{\xi^{\prime}_{k}+\gamma+1})&{{\mathbf{I}}_{M}}\end{array}\right]\;\text{and}\;{\hat{\mathbf{Q}}}^{\text{II}}_{k}=\left[\begin{array}[]{cc}{\hat{\mathbf{H}}}_{1}\text{diag}({{\bm{\theta}}}_{\xi^{\prime}_{k}+1})&{{\mathbf{I}}_{M}}\\\
\vdots&\vdots\\\
{\hat{\mathbf{H}}}_{1}\text{diag}({{\bm{\theta}}}_{\xi^{\prime}_{k}+\gamma+1})&{{\mathbf{I}}_{M}}\end{array}\right],$

then the overall received signals in time slots
$i\in[\xi^{\prime}_{k}+1,\xi^{\prime}_{k}+\gamma+1]$ can be expressed as

${\mathbf{y}}^{\text{II}}_{k}\triangleq[\mathbf{y}^{(\xi^{\prime}_{k}+1)};\cdots;\mathbf{y}^{(\xi^{\prime}_{k}+\gamma+1)}]=\sqrt{p}{\mathbf{Q}}^{\text{II}}_{k}[{\bm{\mu}}_{k};{\mathbf{h}}_{d,k}]+{\mathbf{n}}_{k}^{\text{II}},\;\;k\in[2,K],$
(46)

where
${\mathbf{n}}_{k}^{\text{II}}=[{\mathbf{n}}^{(\xi^{\prime}_{k}+1)};\cdots;{\mathbf{n}}^{(\xi^{\prime}_{k}+\gamma+1)}]$.
According to (46), the LS estimation of ${\bm{\mu}}_{k}$ and
${{\mathbf{h}}}_{d,k}$ associated with $U_{k}$ are given by

$[\hat{{\bm{\mu}}}_{k};{\hat{\mathbf{h}}}_{d,k}]=[{\bm{\mu}}_{k};{{\mathbf{h}}}_{d,k}]+({\hat{\mathbf{Q}}}^{\text{II}}_{k})^{\dagger}\Delta{\mathbf{S}}_{k}{\bm{\mu}}_{k}+\frac{1}{\sqrt{p}}({\hat{\mathbf{Q}}}^{\text{II}}_{k})^{\dagger}{\mathbf{n}}_{k}^{\text{II}}.$
(47)

Let
$\Delta{\mathbf{x}}_{k}\triangleq\frac{1}{\sqrt{p}}({\hat{\mathbf{Q}}}^{\text{II}}_{k})^{\dagger}{\mathbf{n}}_{k}^{\text{II}}+({\hat{\mathbf{Q}}}^{\text{II}}_{k})^{\dagger}\Delta{\mathbf{S}}_{k}{\bm{\mu}}_{k}$
denote the CSI error matrix in (47), where the first term
$\frac{1}{\sqrt{p}}({\hat{\mathbf{Q}}}^{\text{II}}_{k})^{\dagger}{\mathbf{n}}_{k}^{\text{II}}$
is caused by the channel noise and the estimated reflected channel associated
with the typical user $\hat{\mathbf{H}}_{1}$ obtained in Phase I, and the
second term
$({\hat{\mathbf{Q}}}^{\text{II}}_{k})^{\dagger}\Delta{\mathbf{S}}_{k}{\bm{\mu}}_{k}$
is only related to $\hat{\mathbf{H}}_{1}$. Then, similar to (39), we can
obtain the covariance matrix of $\Delta{\mathbf{x}}_{k}$ as follows:

$\begin{array}[]{l}\mathbb{E}\\{\Delta{\mathbf{x}}_{k}\Delta{\mathbf{x}}_{k}^{H}\\}=\mathbf{P}_{k,1}^{\text{b}}+\mathbf{P}_{k,2}^{\text{b}},\end{array}$
(48)

where
$\mathbf{P}_{k,1}^{\text{b}}\triangleq\frac{1}{p}\mathbb{E}\big{\\{}({\hat{\mathbf{Q}}}^{\text{II}}_{k})^{\dagger}{\mathbf{n}}^{\text{II}}_{k}({\mathbf{n}}^{\text{II}}_{k})^{H}\big{(}({\hat{\mathbf{Q}}}^{\text{II}}_{k})^{H}\big{)}^{\dagger}\big{\\}}$
and
$\mathbf{P}_{k,2}^{\text{b}}\triangleq\mathbb{E}\big{\\{}({\hat{\mathbf{Q}}}^{\text{II}}_{k})^{\dagger}\Delta{\mathbf{S}}_{k}{\bm{\mu}}_{k}{\bm{\mu}}_{k}^{H}\Delta{\mathbf{S}}_{k}^{H}(({\hat{\mathbf{Q}}}^{\text{II}}_{k})^{\dagger})^{H}\big{\\}}$.

Similar to the case of $M\geq N$, we will simplify
$\mathbf{P}_{k,1}^{\text{b}}$ and $\mathbf{P}_{k,2}^{\text{b}}$ in the
following. First, since
$\mathbb{E}\big{\\{}{\mathbf{n}}^{\text{II}}_{k}({\mathbf{n}}^{\text{II}}_{k})^{H}\big{\\}}=\sigma^{2}{\mathbf{I}}$,
we obtain
$\mathbf{P}_{k,1}^{\text{b}}=\frac{\sigma^{2}}{p}\mathbb{E}\big{\\{}\big{(}({\hat{\mathbf{Q}}}^{\text{II}}_{k})^{H}{\hat{\mathbf{Q}}}^{\text{II}}_{k}\big{)}^{-1}\big{\\}}$.
Besides, as
$({\hat{\mathbf{Q}}}^{\text{II}}_{k})^{H}{\hat{\mathbf{Q}}}^{\text{II}}_{k}$
is a $2\times 2$ block matrix which can be expressed as
$({\hat{\mathbf{Q}}}^{\text{II}}_{k})^{H}{\hat{\mathbf{Q}}}^{\text{II}}_{k}=[\hat{\mathbf{S}}_{k}^{H}\hat{\mathbf{S}}_{k},\mathbf{0}_{N\times
M};\mathbf{0}_{M\times N},(\gamma+1)\mathbf{I}_{M}]$ with
$\hat{\mathbf{S}}_{k}=[{\hat{\mathbf{H}}}_{1}\text{diag}({{\bm{\theta}}}_{\xi_{k}+1});\cdots;{\hat{\mathbf{H}}}_{1}\text{diag}({{\bm{\theta}}}_{\xi_{k}+\gamma+1})]$,
we can further transform $\mathbf{P}_{k,1}^{\text{b}}$ into

$\mathbf{P}_{k,1}^{\text{b}}{=}\frac{\sigma^{2}}{p}\left[\begin{array}[]{cc}\mathbb{E}\big{\\{}(\hat{\mathbf{S}}_{k}^{H}\hat{\mathbf{S}}_{k})^{-1}\big{\\}}&\mathbf{0}_{N\times
M}\\\ \mathbf{0}_{M\times
N}&\frac{1}{\gamma+1}\mathbf{I}_{M}\end{array}\right],$ (49)

according to Lemma 2. Moreover, by employing the result in Lemma 3,
$\mathbf{P}_{k,2}^{\text{b}}$ can be re-expressed as follows:

$\displaystyle\mathbf{P}_{k,2}^{\text{b}}=\mathbb{E}\big{\\{}\big{(}({\hat{\mathbf{Q}}}^{\text{II}}_{k})^{H}(\Delta{\mathbf{S}_{k}}{\bm{\mu}}_{k}{\bm{\mu}}_{k}^{H}\Delta{\mathbf{S}}^{H}_{k})^{\dagger}{\hat{\mathbf{Q}}}^{\text{II}}_{k}\big{)}^{\dagger}\big{\\}}=\left[\begin{array}[]{cc}\mathbb{E}\big{\\{}{\hat{\mathbf{S}}}_{k}^{\dagger}\Delta{\mathbf{S}}_{k}{\bm{\mu}}_{k}{\bm{\mu}}_{k}^{H}\Delta{\mathbf{S}}_{k}^{H}({\hat{\mathbf{S}}}_{k}^{\dagger})^{H}\big{\\}}&\mathbf{0}_{N\times
M}\\\ \mathbf{0}_{M\times N}&\mathbf{0}_{M\times M}\end{array}\right].$ (50)

Based on the simplified $\mathbf{P}_{k,1}^{\text{b}}$ and
$\mathbf{P}_{k,2}^{\text{b}}$ in (49) and (50), the MSEs for estimating
$\mathbf{h}_{d,k}$ and $\bm{\mu}_{k},k\in[2,K]$, are given by

$\epsilon_{d,k,\text{2P}}^{b}\triangleq\mathbb{E}\big{\\{}\text{tr}(\Delta\mathbf{h}_{d,k}\Delta\mathbf{h}_{d,k}^{H})\big{\\}}=\frac{M\sigma^{2}}{p(\gamma+1)},$
(51a)
$\epsilon_{\bm{\mu}_{k},\text{2P}}^{b}\triangleq\mathbb{E}\big{\\{}\text{tr}(\Delta{\mathbf{x}}_{k}\Delta{\mathbf{x}}_{k}^{H})\big{\\}}=\frac{\sigma^{2}}{p}\mathbb{E}\big{\\{}\text{tr}\big{(}(\hat{\mathbf{S}}^{H}_{k}\hat{\mathbf{S}}_{k})^{-1}\big{)}\big{\\}}+\mathbb{E}\big{\\{}\text{tr}({\hat{\mathbf{S}}_{k}}^{\dagger}\Delta{\mathbf{S}}_{k}{\bm{\mu}}_{k}{\bm{\mu}}_{k}^{H}\Delta{\mathbf{S}}^{H}_{k}({\hat{\mathbf{S}}}_{k}^{\dagger})^{H})\big{\\}}.$
(51b)

Note that different from the case of $M\geq N$ where the MSE for estimating
$\mathbf{h}_{d,k}$, i.e., $\epsilon_{d,k,\text{2P}}^{a}$, is related to $K$,
$\epsilon_{d,k,\text{2P}}^{b}$ is related to $\gamma$ here since
$\mathbf{h}_{d,k}$ is estimated based on the received signals in $\gamma+1$
time slots.

Let
$\mathbf{P}_{\bm{\mu}_{k},1}^{b}\triangleq\frac{\sigma^{2}}{p}\mathbb{E}\big{\\{}\text{tr}\big{(}(\hat{\mathbf{S}}^{H}_{k}\hat{\mathbf{S}}_{k})^{-1}\big{)}\big{\\}}$
and
$\mathbf{P}_{\bm{\mu}_{k},2}^{b}\triangleq\mathbb{E}\big{\\{}\text{tr}({\hat{\mathbf{S}}_{k}}^{\dagger}\Delta{\mathbf{S}}_{k}{\bm{\mu}}_{k}{\bm{\mu}}_{k}^{H}\Delta{\mathbf{S}}^{H}_{k}({\hat{\mathbf{S}}}_{k}^{\dagger})^{H})\big{\\}}$
represent the two terms on the RHS of (51b). In the sequel, we focus on
further simplifying $\mathbf{P}_{\bm{\mu}_{k},1}^{b}$ and
$\mathbf{P}_{\bm{\mu}_{k},2}^{b}$ when $M$ becomes large as the following
proposition.

######  Proposition 2.

$\mathbf{P}_{\bm{\mu}_{k},1}^{\text{b}}$ and
$\mathbf{P}_{\bm{\mu}_{k},2}^{\text{b}}$ can be approximated by
$\mathbf{P}_{\bm{\mu}_{k},1}^{b}\approx\frac{\sigma^{2}N(N+1)}{2pM(N+1)l^{\text{UI}}_{1}l^{\text{IB}}+2M\sigma^{2}}$
and
$\mathbf{P}_{\bm{\mu}_{k},2}^{b}\approx\frac{(\delta\gamma^{2}+2\gamma+2\delta-M)\sigma^{2}\sum_{n=1}^{(\gamma-1)M}\mathbb{E}\\{\mu_{k,n}\mu_{k,n}^{*}\\}}{p\delta(\gamma+1)^{2}(N+1)l^{\text{UI}}_{1}l^{\text{IB}}+(\delta\gamma^{2}+2\gamma+2\delta-M)\sigma^{2}}+\frac{(M\gamma^{2}+(3M-1)\gamma+\delta)\sigma^{2}\sum_{n=(\gamma-1)M+1}^{N}\mathbb{E}\\{\mu_{k,n}\mu_{k,n}^{*}\\}}{p\delta(\gamma+1)^{2}(N+1)l^{\text{UI}}_{1}l^{\text{IB}}+(M\gamma^{2}+(3M-1)\gamma+\delta)\sigma^{2}}$
when $M$ becomes asymptotically large.

###### Proof.

Please refer to Appendix C. ∎

By substituting the asymptotic values of $\mathbf{P}_{\bm{\mu}_{k},1}^{b}$ and
$\mathbf{P}_{\bm{\mu}_{k},2}^{b}$ into (51b), the asymptotic MSE for
estimating $\bm{\mu}_{k}$ can be expressed as follows:

$\displaystyle\epsilon_{\bm{\mu}_{k},\text{2P}}^{b}\approx$
$\displaystyle\frac{\sigma^{2}N(N+1)}{2pM(N+1)l^{\text{UI}}_{1}l^{\text{IB}}+2M\sigma^{2}}+\frac{(\delta\gamma^{2}+2\gamma+2\delta-M)\sigma^{2}\sum_{n=1}^{(\gamma-1)M}\mathbb{E}\\{\mu_{k,n}\mu_{k,n}^{*}\\}}{p\delta(\gamma+1)^{2}(N+1)l^{\text{UI}}_{1}l^{\text{IB}}+(\delta\gamma^{2}+2\gamma+2\delta-M)\sigma^{2}}$
(52)
$\displaystyle+\frac{(M\gamma^{2}+(3M-1)\gamma+\delta)\sigma^{2}\sum_{n=(\gamma-1)M+1}^{N}\mathbb{E}\\{\mu_{k,n}\mu_{k,n}^{*}\\}}{p\delta(\gamma+1)^{2}(N+1)l^{\text{UI}}_{1}l^{\text{IB}}+(M\gamma^{2}+(3M-1)\gamma+\delta)\sigma^{2}}.$

Similar to (44), we can see that as $p$, $l^{\text{UI}}_{1}$ and
$l^{\text{IB}}$ decrease, $\epsilon_{\bm{\mu}_{k},\text{2P}}^{b}$ increases,
and the first term on the RHS of (52) dominates the other terms when $N$ is
asymptotically large. Besides, different from the case of $M\geq N$,
$\epsilon_{\bm{\mu}_{k},\text{2P}}^{b}$ is related to $\gamma$ and $\delta$ in
this case, which are determined by the quantitative relationship between $M$
and $N$.

###  IV-C Asymptotic MSE Comparison between 2PCE and 3PCE Strategies

In this subsection, we compare the performance of the 2PCE and 3PCE strategies
in terms of the asymptotic MSE. For a fair comparison, the asymptotic MSE
performance of the 3PCE strategy is analyzed in a similar way as those in
Section IV-A and Section IV-B, by assuming an LS estimator. The derived
results are summarized in Table I,

TABLE I: MSE performance of the 3PCE strategy in [22] $\epsilon_{d,\text{3P}}$ | $\frac{M\sigma^{2}}{p}$
---|---
$\epsilon_{r,1,\text{3P}}$ | $\frac{(1+K)M\sigma^{2}}{pK}$
$\epsilon_{r,k,\text{3P}}^{\text{a}}$ | $\frac{N(N-1)(1+K)\sigma^{2}}{KM\sigma^{2}+pKMNl^{\text{UI}}_{1}l^{\text{IB}}}+\frac{K(1+K)N\sigma^{2}}{K(K+N)M\sigma^{2}+pMNK^{2}l^{\text{UI}}_{1}l^{\text{IB}}}$
|
$+\frac{(K+N)M\sigma^{2}\mathbb{E}\\{\mu_{k,1}\mu_{k,1}^{*}\\}}{pN^{2}Kl^{\text{UI}}_{1}l^{\text{IB}}+(K+N)M\sigma^{2}}+\frac{M\sigma^{2}\sum_{n=2}^{N}\mathbb{E}\\{\mu_{k,n}\mu_{k,n}^{*}\\}}{pN^{2}l^{\text{UI}}_{1}l^{\text{IB}}+M\sigma^{2}}$
$\epsilon_{\bm{\mu}_{k},\text{3P}}^{b}$ | $\frac{K(1+K)N\sigma^{2}}{pK^{2}NMl^{\text{UI}}_{1}l^{\text{IB}}+K(K+N)M\sigma^{2}}+\frac{N(N-1)(1+K)\sigma^{2}}{pNMKl^{\text{UI}}_{1}l^{\text{IB}}+KM\sigma^{2}}+\frac{M\sigma^{2}\sum_{n=2}^{(\gamma-1)M}\mathbb{E}\\{\mu_{k,n}\mu_{k,n}^{*}\\}}{pNMl^{\text{UI}}_{1}l^{\text{IB}}+M\sigma^{2}}$
|
$+\frac{(K+N)M\sigma^{2}\mathbb{E}\\{\mu_{k,1}\mu_{k,1}^{*}\\}}{pKNMl^{\text{UI}}_{1}l^{\text{IB}}+(K+N)M\sigma^{2}}+\frac{M\sigma^{2}\sum_{n=(\gamma-1)M+1}^{N}\mathbb{E}\\{\mu_{k,n}\mu_{k,n}^{*}\\}}{pN\delta
l^{\text{UI}}_{1}l^{\text{IB}}+M\sigma^{2}}$

where $\epsilon_{d,\text{3P}}$ is the MSE for estimating all the direct
channels, $\epsilon_{r,1,\text{3P}}$ represents the MSE for estimating the
reflected channel associated with $U_{1}$,
$\epsilon_{\bm{\mu}_{k},\text{3P}}^{\text{a}}$ and
$\epsilon_{\bm{\mu}_{k},\text{3P}}^{b}$ denote the asymptotic MSEs for
estimating all the scaling factors in the cases of $M\geq N$ and $M<N$,
respectively. Note that in the 3PCE strategy,
$\\{\mathbf{h}_{d,k}\\}_{k=1}^{K}$, $\mathbf{H}_{1}$ and
$\\{\bm{\mu}_{k}\\}_{k=2}^{K}$ are estimated successively, thus
$\epsilon_{d,\text{3P}}$ is mainly due to the channel noise,
$\epsilon_{r,1,\text{3P}}$ depends on the channel noise and the imperfect
estimation of $\mathbf{h}_{d,1}$, while
$\epsilon_{\bm{\mu}_{k},\text{3P}}^{\text{a}}$
($\epsilon_{\bm{\mu}_{k},\text{3P}}^{b}$) is affected by the channel noise,
the imperfect estimation of $\mathbf{h}_{d,1}$ and the imperfect estimation of
$\mathbf{H}_{1}$. In the sequel, we focus on the performance comparison
between these two strategies and the main results are given in the following
two theorems.

######  Theorem 1.

In the case of $M\geq N$, the asymptotic MSEs achieved by the 2PCE strategy
for estimating all the direct and reflected channels are lower than that
achieved by the 3PCE strategy.

###### Proof.

Please refer to Appendix D. ∎

######  Theorem 2.

In the case of $M<N$, the MSE achieved by the 2PCE strategy for estimating all
the direct channels is lower than that achieved by the 3PCE strategy, if
$\gamma<K+2$. For estimating the reflected channels, the proposed 2PCE
strategy always outperforms the 3PCE strategy in terms of asymptotic MSE.

###### Proof.

Please refer to Appendix E. ∎

##  V Simulation Results

In this section, we present numerical results to verify the effectiveness of
the proposed 2PCE strategy. In our simulations, the distance-dependent path
losses of ${\mathbf{h}}_{d,k}$’s, ${\mathbf{h}}_{r,k}$’s and $\mathbf{G}$ are
modeled as
$l_{k,m}^{\text{UB}}=l_{0}\left({d_{k,m}}/{d_{0}}\right)^{-\alpha^{\text{UB}}}$,
$l_{k,n}^{\text{UI}}=l_{0}\left({d_{k,n}}/{d_{0}}\right)^{-\alpha^{\text{UI}}}$
and
$l^{\text{IB}}_{m,n}=l_{0}\left({d_{m,n}}/{d_{0}}\right)^{-\alpha^{\text{IB}}}$,
respectively, where $l_{0}=-30$ dB denotes the path loss at the reference
distance $d_{0}=1$ meter (m), $d_{k,m}$, $d_{k,n}$ and $d_{m,n}$ represent the
link distances from $U_{k}$ to the $m$-th BS antenna, from $U_{k}$ to the
$n$-th reflecting element and from the $n$-th reflecting element to the $m$-th
BS antenna, respectively, and $\alpha^{\text{UB}}$, $\alpha^{\text{UI}}$ and
$\alpha^{\text{IB}}$ are path-loss exponents for the user-BS, user-IRS and
IRS-BS links, respectively. We assume that the IRS is deployed to serve the
users that suffer from severe signal attenuation in the user-BS link, thus the
path-loss exponents are set to $\alpha^{\text{UB}}=5$,
$\alpha^{\text{IB}}=2.2$ and $\alpha^{\text{UI}}=2.2$. We consider a three-
dimensional coordinate system where the BS/IRS are deployed on the $x$-axis
and $y$-$z$ plane, and we fix $N_{y}=4$. The reference antenna/element at the
BS/IRS are located at $(2\;\textrm{m},0,0)$ and
$(0,45\;\textrm{m},2\;\textrm{m})$, and the antenna/element spacing are set to
${d_{B}}={\mu}/{2}$ and ${d_{I}}={\mu}/{8}$. The power spectrum of the AWGN at
the BS and the channel bandwidth are -169 dBm/Hz and 1 MHz. Besides, we employ
the normalized MSEs, i.e.,
$\text{NMSE}_{d}=\frac{\sum_{k=1}^{K}\mathbb{E}\left\\{\|\hat{\mathbf{h}}_{d,k}-\mathbf{h}_{d,k}\|^{2}\right\\}}{\sum_{k=1}^{K}\mathbb{E}\left\\{\|\mathbf{h}_{d,k}\|^{2}\right\\}}$
and
$\text{NMSE}_{r}=\frac{\sum_{k=1}^{K}\mathbb{E}\left\\{\|\hat{\mathbf{H}}_{k}-\mathbf{H}_{k}\|^{2}\right\\}}{\sum_{k=1}^{K}\mathbb{E}\left\\{\|\mathbf{H}_{k}\|^{2}\right\\}}$,
as the performance metrics for estimating the direct and reflected channels,
respectively. Other system parameters are set as follows unless otherwise
specified: $K=4,N=32,M=40,\beta^{\text{IB}}=3$ dB,
$\beta^{\text{UB}}=\beta^{\text{UI}}=0$, $r_{d}=r_{r}=r_{r,k}=0$. All the
results are averaged over 2000 independent channel realizations.

####  V-1 Impact of the transmit power, $p$

Figure 3: NMSE performance versus the transmit power, $p$.

Figure 4:  NMSE performance versus the path-loss exponents,
$\alpha^{\text{UB}}$ and $\alpha^{\text{IB}}$.

First, we investigate in Fig. 4 the NMSE performance of the proposed 2PCE
strategy versus the transmit power $p$, and the performance of the 3PCE
strategy [22] is regarded as the performance benchmark. Let
$\text{NMSE}_{d,1}$ ($\text{NMSE}_{r,1}$) and $\text{NMSE}_{d,k}$
($\text{NMSE}_{r,k}$) represent the NMSEs for estimating $\mathbf{h}_{d,1}$
($\mathbf{H}_{1}$) and $\\{\mathbf{h}_{d,k}\\}_{k=2}^{K}$
($\\{\mathbf{H}_{k}\\}_{k=2}^{K}$), respectively. First, it is observed that
the NMSE performance of both strategies improves with the increasing of $p$,
which is expected since larger $p$ implies higher channel training SNR.
Besides, we can observe from Fig. 4 (a) that for the typical user, the 2PCE
strategy can achieve much better NMSE performance than the 3PCE strategy when
estimating the direct channel $\mathbf{h}_{d,1}$, which is mainly due to the
fact that $\mathbf{h}_{d,1}$ is obtained based on the received signals in a
larger total number of time slots as compared to the 3PCE strategy. Moreover,
the 2PCE strategy outperforms the 3PCE strategy when estimating
$\mathbf{H}_{1}$, which is because the estimation error of $\mathbf{h}_{d,1}$
will not affect the estimation of $\mathbf{H}_{1}$ in the proposed 2PCE
strategy. Due to similar reasons, the NMSE performance achieved by the
proposed 2PCE strategy for the estimated channels associated with the other
users is also better than that achieved by the 3PCE strategy, as can be
observed from Fig. 4 (b). Specifically, to achieve the same NMSE performance,
the transmit power required by the proposed 2PCE strategy is about 2 dB less
than that required by the 3PCE strategy. Finally, it is noteworthy that this
performance gain is available for all training SNR regimes.

####  V-2 Impact of the path-loss coefficients, $\alpha^{\text{UB}}$ and
$\alpha^{\text{IB}}$

Fig. 4 shows the NMSE performance versus the path-loss exponents of the user-
BS link $\alpha^{\text{UB}}$ and the IRS-BS link $\alpha^{\text{IB}}$, where
the user transmit power is set to $p=16$ dBm. Besides, we fix
$\alpha^{\text{IB}}=2.2$ and $\alpha^{\text{UI}}=2.2$ in Fig. 4 (a), while in
Fig. 4 (b), we fix $\alpha^{\text{UB}}=5$ and $\alpha^{\text{UI}}=2.2$. It is
observed that the proposed 2PCE strategy outperforms the 3PCE strategy,
similar to the results in Fig. 4. Besides, from Fig. 4 (a), we observe that
the NMSEs achieved by the 2PCE and 3PCE strategies for estimating the direct
channels decrease with the increasing of $\alpha^{\text{UB}}$, while those for
estimating the reflected channels remain unchanged with different values of
$\alpha^{\text{UB}}$. This is because when $\alpha^{\text{UB}}$ increases, the
path loss of each user-BS link becomes more severe and it becomes more
difficult to accurately estimate the direct channels. Similarly, from Fig. 4
(b), it can be seen that the NMSEs achieved by these two strategies for
estimating the reflected channels decrease as $\alpha^{\text{IB}}$ increases,
while those for estimating the direct channels are not affected by the value
of $\alpha^{\text{IB}}$. For the 3PCE strategy, this is because the estimation
of the direct channels is conducted when the IRS is switched off and thus not
affected by the estimation of the reflected channels. While for the proposed
2PCE strategy, the results in Fig. 4 (b) are consistent with the analysis in
Section IV (see (43a)), i.e., although the estimation accuracy of the direct
channels associated with the other users in Phase II depends on the knowledge
of $\mathbf{H}_{1}$, the MSE performance is not affected by the estimation
accuracy of $\mathbf{H}_{1}$. Note that by fixing $\alpha^{\text{UB}}=5$,
$\alpha^{\text{IB}}=2.2$ and investigating the impact of $\alpha^{\text{UI}}$
on the NMSE performance, similar results as those in Fig. 4 (b) can be
obtained and therefore the details are omitted for brevity.

Figure 5: NMSE performance versus the number of BS antennas, $M$.

Figure 6:  NMSE performance versus the number of reflecting elements, $N$.

####  V-3 Impact of the number of BS antennas $M$

In Fig. 6, we investigate the impact of the number of BS antennas $M$ on the
NMSE performance, where we fix $N=20$, $p=24$ dBm and $K=2$. First, we observe
that the proposed 2PCE strategy outperforms the 3PCE strategy, as in Figs. 4
and 4. Second, from Fig. 6 (a), we observe that the NMSE performance is almost
invariant with different values of $M$. This is because the MSEs for
estimating $\mathbf{h}_{d,1}$ and $\mathbf{H}_{1}$ in the 2PCE strategy, i.e.,
$\epsilon_{d,1,\text{2P}}$ and $\epsilon_{r,1,\text{2P}}$, increase with $M$
(according to (32)), thus the corresponding NMSEs, i.e., $\text{NMSE}_{d,1}$
and $\text{NMSE}_{r,1}$, are not related to $M$ due to the normalization.
Then, one can see from Fig. 6 (b) that in the case of $M<N$, the
$\text{NMSE}_{d,k}$ performance of the 2PCE strategy deteriorates with the
increasing $M$ because the asymptotic MSE for estimating
$\mathbf{h}_{d,k},k\neq 1$, i.e., $\epsilon_{d,k,\text{2P}}^{b}$, is inversely
proportional to $\gamma+1$ according to (51a) and the value of $\gamma$
decreases when $M$ becomes larger. In the case of $M\geq N$, since
$\epsilon_{d,k,\text{2P}}^{a}$ is linearly proportional to $M$ according to
(43a), the $\text{NMSE}_{d,k}$ achieved by the 2PCE strategy is not related to
$M$ because of the normalization, and it is similar to that achieved by the
3PCE strategy. Besides, for the estimation of the reflected channels
associated with the other users, the NMSE performance achieved by the 2PCE
strategy does not change much when $M<N$, while it improves with the
increasing of $M$ when $M\geq N$. The main reason is that in the proposed 2PCE
strategy, the inverses of the $2M\times(M+N)$ matrix
${\hat{\mathbf{V}}}^{\text{II}}$ and the $(\gamma+1)M\times(M+N)$ matrix
${\hat{\mathbf{Q}}}^{\text{II}}_{k}$ are required to estimate
$\mathbf{H}_{k},k\neq 1$, for the cases of $M\geq N$ and $M<N$, respectively.
Let $\varrho_{V}=\frac{2M}{M+N}$ and $\varrho_{Q}=\frac{(\gamma+1)M}{M+N}$
denote the ratios of the number of rows to the number of columns of the
matrices ${\hat{\mathbf{V}}}^{\text{II}}$ and
${\hat{\mathbf{Q}}}^{\text{II}}_{k}$, respectively. Then, we can see that
$\varrho_{V}$ increases with $M$ when $M\geq N$ and
$\hat{\mathbf{V}}^{\text{II}}$ will become less ill-conditioned as $M$
increases, which causes less performance loss when estimating
$\mathbf{H}_{k}$. However, when $M<N$, $\varrho_{M}$ satisfies
$\varrho_{Q}\geq 1$ and $\varrho_{Q}\approx 1$ with different values of $M$,
therefore $\text{NMSE}_{r,k}$ is almost invariant with $M$ in this case.

####  V-4 Impact of the number of reflecting elements $N$

In Fig. 6, we investigate the NMSE performance versus the number of reflecting
elements $N$ with fixed $M=20$, $p=24$ dBm and $K=2$. Similar to the results
presented above, we can observe that the NMSE performance of the 2PCE strategy
is superior to that of the 3PCE strategy. From Fig. 6 (a), it can be seen that
the NMSE performance achieved by the 2PCE strategy for estimating
$\mathbf{h}_{d,1}$ and ${\mathbf{H}}_{1}$ improves with $N$, which is mainly
because both $\mathbf{h}_{d,1}$ and $\mathbf{H}_{1}$ are estimated based on
the received signals in $N+1$ time slots. Fig. 6 (b) shows that
$\text{NMSE}_{d,k}$ remains unchanged with varying $N$ when $M\geq N$. This is
due to the fact that $\epsilon_{d,k,\text{2P}}^{a}$ is not related to $N$
according to (43a). However, when $M<N$, $\text{NMSE}_{d,k}$ decreases with
the increasing of $N$, which is consistent with analysis in (51a).
Furthermore, focusing on the estimation of $\mathbf{H}_{k},k\neq 1$, we can
see that the NMSE performance achieved by the 2PCE strategy deteriorates with
the increasing of $N$ due to the reduced value of $\varrho_{V}$ when $M\geq N$
(similar to that in Fig. 6 (b)), while the performance is nearly invariant
when $M<N$ since $\varrho_{Q}\approx 1$ holds under different values of $N$.

Figure 7: NMSE performance versus the Rician factors, $\beta^{\text{IB}}$ and
$\beta^{\text{UI}}$.

Figure 8: NMSE performance versus the correlation coefficient, $r_{d}$.

####  V-5 Impact of the Rician factors, $\beta^{\text{IB}}$ and
$\beta^{\text{UI}}$

Fig. 8 investigates the NMSE performance versus the Rician factors of the IRS-
BS link $\beta^{\text{IB}}$ and the user-IRS link $\beta^{\text{UI}}$ with
fixed $p=16$ dBm. First, in Fig. 8 (a), we fix
$\beta^{\text{UI}}=\beta^{\text{UB}}=0$ and it can be seen that the NMSEs
achieved by the 2PCE and 3PCE strategies for estimating the reflected channels
$\mathbf{H}_{k},k\neq 1$, increase with the increasing of $\beta^{\text{IB}}$,
while the $\text{NMSE}_{d,1}$, $\text{NMSE}_{d,k}$ and $\text{NMSE}_{r,1}$
performance do not change much with varying $\beta^{\text{IB}}$. The main
reason is that as $\beta^{\text{IB}}$ increases, the IRS-BS channel becomes
more deterministic and the rank of the channel decreases in general (according
to the channel model in (2)), thus the channel matrix $\mathbf{H}_{1}$ becomes
more ill-conditioned, which will affect the estimation of $\mathbf{H}_{k}$
when calculating the pseudo-inverse of $\hat{\mathbf{V}}^{\text{II}}$ in the
proposed 2PCE strategy (or $\hat{\mathbf{H}}_{1}$ in the 3PCE strategy).
Second, in Fig. 8 (b), we fix $\beta^{\text{IB}}=0$ and $\beta^{\text{UB}}=0$
and it can be observed that the NMSE performance achieved by these two
strategies improves with the increasing of $\beta^{\text{UI}}$, which is
expected since the channel matrix $\mathbf{H}_{1}$ becomes less ill-
conditioned as $\beta^{\text{UI}}$ increases.666 According to the user-IRS
channel model in (2), as $\beta^{\text{UI}}$ increases, the user-IRS channel
associated with $U_{k}$, i.e., $\mathbf{h}_{r,k}$, becomes more deterministic
and more like its LoS element ${\bar{\mathbf{z}}}_{r,k}$ whose entries have
the same amplitude according to (2) and (3). As a result, the reflected
channel $\mathbf{H}_{k}$, whose $n$-th column is the scalar-vector
multiplication of the $n$-th entry of $\mathbf{h}_{r,k}$ and the $n$-th column
of the IRS-BS channel $\mathbf{G}$ (see the reflected channel model given in
(6)), is less ill-conditioned if $\beta^{\text{UI}}$ is larger.  Furthermore,
one can observe that the benefit brought by increasing $\beta^{\text{UI}}$ is
more pronounced on the proposed 2PCE strategy than that on the 3PCE strategy.

####  V-6 Impact of the correlation coefficient, $r_{d}$

Finally, we plot in Fig. 8 the NMSE performance achieved by the 2PCE and 3PCE
strategies versus the correlation coefficients $r_{d}$, where we assume
$r_{r}=r_{r,k}=0$ and $p=16$ dBm. As can be observed, only the NMSE
performance for estimating $\\{\mathbf{H}_{k}\\}_{k=2}^{K}$ is affected by the
increasing of $r_{d}$, while the impact of $r_{d}$ on $\text{NMSE}_{d,1}$,
$\text{NMSE}_{d,k}$ and $\text{NMSE}_{r,1}$ is not that significant. This is
due to the fact that as $r_{d}$ increases, the IRS-BS channel $\mathbf{G}$
becomes more correlated and thus the reflected channel $\mathbf{H}_{1}$
becomes more ill-conditioned, which could result in NMSE performance loss.
Besides, it is noteworthy that the proposed 2PCE strategy outperforms the 3PCE
strategy for all values of $r_{d}$. Similar result as in Fig. 8 can be
observed when fixing $r_{d}=r_{r,k}=0$ and investigating the impact of $r_{r}$
on the NMSE performance, while the NMSE performance is shown to be invariant
with $r_{r,k}$ by assuming $r_{d}=r_{r}=0$ and studying the impact of
$r_{r,k}$, hence their details are not shown here.

##  VI Conclusion

This work studied the uplink channel estimation problem for an IRS-aided
multiuser SIMO system and proposed a novel 2PCE strategy. By reducing the
channel estimation phases and estimating the direct and reflected channels
associated with each user simultaneously, the proposed 2PCE strategy is able
to alleviate the negative effects caused by error propagation and enhance the
channel estimation performance without increasing the channel training
overhead. In addition, the asymptotic MSE of the proposed 2PCE strategy is
analyzed when the LS channel estimation method is employed. Both theoretical
analysis and simulation results demonstrated that the proposed 2PCE strategy
can achieve better performance than the state-of-the-art 3PCE strategy [22]
under various system parameters.

Appendix A Necessary propositions and lemmas

######  Proposition 3.

For any random matrix $\mathbf{X}\in\mathbf{C}^{p\times n}$ with $p\leq n$, if
$\mathbb{E}\\{\mathbf{X}{\mathbf{X}}^{H}\\}$ is a diagonal matrix, i.e., the
columns of $\mathbf{X}$ are independent with each other, then
$\mathbb{E}\big{\\{}(\mathbf{X}\mathbf{X}^{H})^{-1}\big{\\}}=\big{(}\mathbb{E}\\{\mathbf{X}{\mathbf{X}}^{H}\\}\big{)}^{-1}$
holds as $n\rightarrow\infty$.

###### Proof.

First, it is easy to see that if $\mathbb{E}\\{\mathbf{X}{\mathbf{X}}^{H}\\}$
is a diagonal matrix, then each column of $\mathbf{X}$ is independent with
each other. By defining
$\mathbf{C}_{\mathbf{X}}=\mathbb{E}\\{\mathbf{X}{\mathbf{X}}^{H}\\}$ and
$\mathbf{Y}=\mathbf{C}_{\mathbf{X}}^{-\frac{1}{2}}\mathbf{X}$, we can see that
the random matrix $\mathbf{Y}=[\mathbf{y}_{1},\cdots,\mathbf{y}_{n}]$ consists
of $n$ independent random vectors $\\{\mathbf{y}_{i}\\}_{i=1}^{n}$, which are
from the same distribution with covariance matrix
$\mathbb{E}\\{\mathbf{y}\mathbf{y}^{H}\\}=\frac{1}{n}\mathbf{I}$. Furthermore,
defining $\hat{\mathbf{C}}_{\mathbf{y},n}$ as the sample covariance matrix of
$\mathbf{y}$ defined over $n$ samples, i.e.,
$\hat{\mathbf{C}}_{\mathbf{y},n}=\frac{1}{n}\sum\limits_{i=1}^{n}\mathbf{y}_{i}\mathbf{y}_{i}^{H}=\frac{1}{n}\mathbf{Y}\mathbf{Y}^{H}$,
then $\hat{\mathbf{C}}_{\mathbf{y},n}$ approximates the true covariance matrix
$\mathbb{E}\\{\mathbf{y}\mathbf{y}^{H}\\}$ as $n\rightarrow\infty$, and
equivalently we have $\mathbf{Y}\mathbf{Y}^{H}\rightarrow\mathbf{I}$ as
$n\rightarrow\infty$ [29, 30, 31]. Therefore, according to the relationship
between $\mathbf{X}$ and $\mathbf{Y}$, i.e.,
$\mathbf{Y}=\mathbf{C}_{\mathbf{X}}^{-\frac{1}{2}}\mathbf{X}$, it follows that
$\mathbf{X}\mathbf{X}^{H}=\mathbf{C}_{\mathbf{X}}.$ Since
$\mathbf{C}_{\mathbf{X}}$ is a constant diagonal matrix, we have
$(\mathbf{X}\mathbf{X}^{H})^{-1}\rightarrow\mathbf{C}_{\mathbf{X}}^{-1}$ as
$n\rightarrow\infty$ and hence
$\mathbb{E}\big{\\{}(\mathbf{X}\mathbf{X}^{H})^{-1}\big{\\}}=\mathbf{C}_{\mathbf{X}}^{-1}$
holds when $n\rightarrow\infty$. This completes the proof.

∎

######  Lemma 2.

[32] Given a $2\times 2$ block matrix
$\mathbf{M}=[\mathbf{A},\mathbf{B};\mathbf{C},\mathbf{D}]$, if both
$\mathbf{A}$ and $\mathbf{D}$ are square and both $\mathbf{A}$ and
$\mathbf{D}-\mathbf{CA}^{-1}\mathbf{B}$ are invertible, then the inverse of
$\mathbf{M}$ can be expressed as follows:

$\left[\begin{array}[]{cc}\mathbf{A}&\mathbf{B}\\\
\mathbf{C}&\mathbf{D}\end{array}\right]^{-1}=\left[\begin{array}[]{cc}\mathbf{A}^{-1}+\mathbf{A}^{-1}\mathbf{B}\left(\mathbf{D}-\mathbf{CA}^{-1}\mathbf{B}\right)^{-1}\mathbf{CA}^{-1}&-\mathbf{A}^{-1}\mathbf{B}\left(\mathbf{D}-\mathbf{CA}^{-1}\mathbf{B}\right)^{-1}\\\
-\left(\mathbf{D}-\mathbf{CA}^{-1}\mathbf{B}\right)^{-1}\mathbf{CA}^{-1}&\left(\mathbf{D}-\mathbf{CA}^{-1}\mathbf{B}\right)^{-1}\end{array}\right].$
(53)

######  Lemma 3.

[33] Given a $2\times 2$ block matrix
$\mathbf{M}=[\mathbf{A},\mathbf{B};\mathbf{C},\mathbf{D}]$, if
$(\mathbf{I}-\mathbf{A}\mathbf{A}^{\dagger})\mathbf{B}=\mathbf{0},\mathbf{C}(\mathbf{I}-\mathbf{A}^{\dagger}\mathbf{A})=\mathbf{0}$
and $\mathbf{D}-\mathbf{C}\mathbf{A}^{\dagger}\mathbf{B}=\mathbf{0}$, then the
pseudo-inverse of $\mathbf{M}$ is given by

$\mathbf{M}^{\dagger}=\left[\begin{array}[]{cc}\left(\mathbf{I}-\mathbf{K}_{B}\tilde{\mathbf{K}}_{B}^{-1}\mathbf{K}_{B}^{H}\right)\mathbf{A}^{\dagger}\left(\mathbf{I}-\mathbf{K}_{C}^{H}\tilde{\mathbf{K}}_{C}^{-1}\mathbf{K}_{C}\right)&\left(\mathbf{I}-\mathbf{K}_{B}\tilde{\mathbf{K}}_{B}^{-1}\mathbf{K}_{B}^{H}\right)\mathbf{A}^{\dagger}\mathbf{K}_{C}^{H}\tilde{\mathbf{K}}_{C}^{-1}\\\
\tilde{\mathbf{K}}_{B}^{-1}\mathbf{K}_{B}^{H}\mathbf{A}^{\dagger}\left(\mathbf{I}-\mathbf{K}_{C}^{H}\tilde{\mathbf{K}}_{C}^{-1}\mathbf{K}_{C}\right)&\tilde{\mathbf{K}}_{B}^{-1}\mathbf{K}_{B}^{H}\mathbf{A}^{\dagger}\mathbf{K}_{C}^{H}\tilde{\mathbf{K}}_{C}^{-1}\end{array}\right],$
(54)

where $\mathbf{K}_{B}=\mathbf{A}^{\dagger}\mathbf{B}$,
$\mathbf{K}_{C}=\mathbf{C}\mathbf{A}^{\dagger}$,
$\tilde{\mathbf{K}}_{B}=\mathbf{I}+\mathbf{K}_{B}^{H}\mathbf{K}_{B}$ and
$\tilde{\mathbf{K}}_{C}=\mathbf{I}+\mathbf{K}_{C}\mathbf{K}_{C}^{H}$.

######  Lemma 4.

For any square matrix $\mathbf{X}\in{\mathbb{C}}^{N\times N}$, we have
$\mathbb{E}\big{\\{}\text{tr}(\mathbf{X})\big{\\}}=\text{tr}\big{(}\mathbb{E}\\{\mathbf{X}\\}\big{)}$.

###### Proof.

Since
$\mathbb{E}\big{\\{}\text{tr}(\mathbf{X})\big{\\}}=\mathbb{E}\\{\sum_{l=1}^{L}\sum_{n=1}^{N}\mathbf{X}_{n,n}\\},$
and
$\text{tr}\big{(}\mathbb{E}\\{\mathbf{X}\\}\big{)}=\sum_{n=1}^{N}\mathbb{E}\\{\mathbf{X}_{n,n}\\}$,
we can easily prove
$\mathbb{E}\big{\\{}\text{tr}(\mathbf{X})\big{\\}}=\text{tr}\big{(}\mathbb{E}\\{\mathbf{X}\\}\big{)}$
as the expected value of the sum of some random variables is equal to the sum
of their individual expected values regardless of whether they are
independent. ∎

######  Lemma 5.

[34] Any product of conformable matrices can be permuted cyclically without
altering the trace of the product, i.e.,
$\text{tr}(\mathbf{ABC})=\text{tr}(\mathbf{BCA})=\text{tr}(\mathbf{CBA}).$

Appendix B Proof of Proposition 1

First, by resorting to Lemma 4, $\mathbf{P}_{\bm{\mu}_{k},1}^{\text{a}}$ can
be transformed into
$\mathbf{P}_{\bm{\mu}_{k},1}^{\text{a}}{=}\frac{\sigma^{2}}{2p(K-1)}\text{tr}\big{(}\mathbb{E}\big{\\{}(\hat{{\mathbf{H}}}_{1}^{H}{\hat{\mathbf{H}}}_{1})^{-1}\big{\\}}\big{)}$.
Then, we have
$\mathbb{E}\\{\hat{{\mathbf{H}}}_{1}^{H}{\hat{\mathbf{H}}}\\}\overset{\text{(a)}}{=}\mathbb{E}\\{\mathbf{H}_{1}^{H}\mathbf{H}_{1}\\}+\mathbb{E}\\{\Delta\mathbf{H}_{1}^{H}\Delta\mathbf{H}_{1}\\}\overset{\text{(b)}}{=}\frac{M(p(N+1)l^{\text{UI}}_{1}l^{\text{IB}}+\sigma^{2})}{p(N+1)}\mathbf{I}_{N}$,
where $(\text{a})$ is due to the fact that the elements in $\mathbf{H}_{1}$
and $\Delta\mathbf{H}_{1}$ are independent with each other, and $(\text{b})$
is obtained based on (31) and the reflected channel model given in (6), thus
$\mathbb{E}\\{\hat{{\mathbf{H}}}_{1}^{H}{\hat{\mathbf{H}}}\\}$ is a diagonal
matrix. Besides, when $M$ becomes asymptotically large, we have
$\mathbb{E}\big{\\{}(\hat{{\mathbf{H}}}_{1}^{H}{\hat{\mathbf{H}}}_{1})^{-1}\big{\\}}\approx\big{(}\mathbb{E}\\{\hat{{\mathbf{H}}}_{1}^{H}{\hat{\mathbf{H}}}\\}\big{)}^{-1}=\frac{p(N+1)}{M(p(N+1)l^{\text{UI}}_{1}l^{\text{IB}}+\sigma^{2})}\mathbf{I}_{N}$
according to Proposition 3. Therefore,
$\mathbf{P}_{\bm{\mu}_{k},1}^{\text{a}}$ can be approximated by
$\mathbf{P}_{\bm{\mu}_{k},1}^{\text{a}}\approx\frac{\sigma^{2}(N+1)N}{2(K-1)M(p(N+1)l^{\text{UI}}_{1}l^{\text{IB}}+\sigma^{2})}$.
Second, $\mathbf{P}_{\bm{\mu}_{k},2}^{\text{a}}$ can be equivalently rewritten
as

$\mathbf{P}_{\bm{\mu}_{k},2}^{\text{a}}\overset{\text{(a)}}{=}\mathbb{E}\Big{\\{}\text{tr}\big{(}\bm{\mu}_{k}\bm{\mu}_{k}^{H}\Delta\mathbf{H}_{1}^{H}(\hat{\mathbf{H}}_{1}^{H})^{\dagger}\hat{\mathbf{H}}_{1}^{\dagger}\Delta\mathbf{H}_{1}\big{)}\Big{\\}}\overset{\text{(b)}}{=}\text{tr}\Big{(}\mathbb{E}\big{\\{}\bm{\mu}_{k}\bm{\mu}_{k}^{H}\big{\\}}\mathbb{E}\big{\\{}(\Delta\mathbf{H}_{1}^{\dagger}\hat{\mathbf{H}}_{1}(\Delta\mathbf{H}_{1}^{\dagger}\hat{\mathbf{H}}_{1})^{H})^{-1}\big{\\}}\Big{)},$
(55)

where $(\text{a})$ is due to Lemma 5, and $(\text{b})$ is due to Lemma 4 and
the fact that the elements in $\bm{\mu}_{k}\bm{\mu}_{k}^{H}$ are independent
to those in
$\Delta\mathbf{H}_{1}^{H}(\hat{\mathbf{H}}_{1}^{H})^{\dagger}\hat{\mathbf{H}}_{1}^{\dagger}\Delta\mathbf{H}_{1}$.
Furthermore, since the elements in $\mathbf{H}_{1}$ and $\Delta\mathbf{H}_{1}$
are independent, it follows that
$\mathbb{E}\big{\\{}\Delta\mathbf{H}_{1}^{\dagger}\hat{\mathbf{H}}_{1}(\Delta\mathbf{H}_{1}^{\dagger}\hat{\mathbf{H}}_{1})^{H}\big{\\}}=\mathbb{E}\big{\\{}\Delta\mathbf{H}_{1}^{\dagger}\mathbb{E}\\{{\mathbf{H}}_{1}{\mathbf{H}}_{1}^{H}\\}(\Delta\mathbf{H}_{1}^{\dagger})^{H}\big{\\}}+\mathbf{I}_{N}\overset{\text{(a)}}{=}Nl^{\text{UB}}_{1}l^{\text{IB}}\mathbb{E}\big{\\{}(\Delta\mathbf{H}_{1}^{H}\Delta\mathbf{H}_{1})^{-1}\big{\\}}+\mathbf{I}_{N}\overset{\text{(b)}}{\approx}\frac{Nl^{\text{UI}}_{1}l^{\text{IB}}p(N+1)+M\sigma^{2}}{M\sigma^{2}}\mathbf{I}_{N}$,
where $(\text{a})$ is because
$\mathbb{E}\\{\mathbf{H}_{1}\mathbf{H}_{1}^{H}\\}=Nl^{\text{UI}}_{1}l^{\text{IB}}\mathbf{I}_{M}$
which is obtained based on the reflected channel model in (6), and (b) is due
to the fact that
$\mathbb{E}\\{\Delta\mathbf{H}_{1}^{H}\Delta\mathbf{H}_{1}\\}=\frac{M\sigma^{2}}{p(N+1)}\mathbf{I}_{N}$
is a diagonal matrix. Since
$\mathbb{E}\big{\\{}\Delta\mathbf{H}_{1}^{\dagger}\hat{\mathbf{H}}_{1}(\Delta\mathbf{H}_{1}^{\dagger}\hat{\mathbf{H}}_{1})^{H}\big{\\}}$
is a diagonal matrix, we obtain
$\mathbb{E}\big{\\{}\big{(}\Delta\mathbf{H}_{1}^{\dagger}\hat{\mathbf{H}}_{1}(\Delta\mathbf{H}_{1}^{\dagger}\hat{\mathbf{H}}_{1})^{H}\big{)}^{-1}\big{\\}}\approx\big{(}\mathbb{E}\big{\\{}\Delta\mathbf{H}_{1}^{\dagger}\hat{\mathbf{H}}_{1}(\Delta\mathbf{H}_{1}^{\dagger}\hat{\mathbf{H}}_{1})^{H}\big{\\}}\big{)}^{-1}$
according to Proposition 3. As a result,
$\mathbf{P}_{\bm{\mu}_{k},2}^{\text{a}}$ can be approximated by
$\mathbf{P}_{\bm{\mu}_{k},2}^{\text{a}}\approx\frac{M\sigma^{2}\mathbb{E}\\{\bm{\mu}_{k}^{H}\bm{\mu}_{k}\\}}{Nl^{\text{UI}}_{1}l^{\text{IB}}p(N+1)+M\sigma^{2}}$.
This thus completes the proof.

Appendix C Proof of Proposition 2

First, by resorting to Lemma 4, $\mathbf{P}_{\bm{\mu}_{k},1}^{b}$ can be
transformed into
$\mathbf{P}_{\bm{\mu}_{k},1}^{b}=\frac{\sigma^{2}}{p}\text{tr}\big{(}\mathbb{E}\big{\\{}(\hat{\mathbf{S}}^{H}_{k}\hat{\mathbf{S}}_{k})^{-1}\big{\\}}\big{)}$.
Then, we define
${\mathbf{S}}_{k}\triangleq[\text{diag}({{\bm{\theta}}}_{\xi_{k}+1}){{\mathbf{H}}}_{1}^{H},\cdots,\text{diag}({{\bm{\theta}}}_{\xi_{k}+\gamma+1}){{\mathbf{H}}}_{1}^{H}]^{H}$
and it is readily seen that the elements of $\mathbf{S}_{k}$ are independent
with those of $\Delta\mathbf{S}_{k}$ due to the fact that $\mathbf{H}_{1}$ and
$\Delta\mathbf{H}_{1}$ are uncorrelated. By replacing $\hat{\mathbf{S}}_{k}$
with ${\mathbf{S}}_{k}-\Delta{\mathbf{S}}_{k}$, we can obtain
$\mathbb{E}\\{\hat{\mathbf{S}}^{H}_{k}\hat{\mathbf{S}}_{k}\\}=\mathbb{E}\\{{{\mathbf{S}}}^{H}_{k}{{\mathbf{S}}}_{k}\\}+\mathbb{E}\\{\Delta{\mathbf{S}}^{H}_{k}\Delta{\mathbf{S}}_{k}\\}=\frac{2pM(N+1)l^{\text{UI}}_{1}l^{\text{IB}}+2M\sigma^{2}}{p(N+1)}\mathbf{I}_{N}$,
which is a diagonal matrix. Besides, when $M$ is asymptotically large, we have
$\mathbb{E}\\{(\hat{\mathbf{S}}^{H}_{k}\hat{\mathbf{S}}_{k})^{-1}\\}\approx(\mathbb{E}\\{\hat{\mathbf{S}}^{H}_{k}\hat{\mathbf{S}}_{k}\\})^{-1}$
according to Proposition 3, then $\mathbf{P}_{\bm{\mu}_{k},1}^{b}$ can be
approximated by
$\mathbf{P}_{\bm{\mu}_{k},1}^{b}\approx\frac{\sigma^{2}N(N+1)}{2pM(N+1)l^{\text{UI}}_{1}l^{\text{IB}}+2M\sigma^{2}}.$
Second, $\mathbf{P}_{\bm{\mu}_{k},2}^{b}$ can be equivalently rewritten as
$\mathbf{P}_{\bm{\mu}_{k},2}^{b}=\text{tr}\big{(}\mathbb{E}\\{\bm{\mu}_{k}\bm{\mu}_{k}^{H}\\}\mathbb{E}\\{\tilde{\mathbf{S}}_{k}^{-1}\\}\big{)}$
with
$\tilde{\mathbf{S}}_{k}=\Delta\mathbf{S}_{k}^{\dagger}\hat{\mathbf{S}}_{k}(\Delta\mathbf{S}_{k}^{\dagger}\hat{\mathbf{S}}_{k})^{H}$
similar to the derivation of (55). Furthermore, since the elements of
$\mathbf{S}_{k}$ and $\Delta\mathbf{S}_{k}$ are independent with each other,
it follows that
$\mathbb{E}\\{\tilde{\mathbf{S}}_{k}\\}=\mathbb{E}\big{\\{}\big{(}\Delta\mathbf{S}^{H}_{k}(\mathbb{E}\\{{{\mathbf{S}}}_{k}{{\mathbf{S}}}^{H}_{k}\\})^{\dagger}\Delta\mathbf{S}_{k}\big{)}^{-1}\big{\\}}+\mathbf{I}_{N}\overset{\text{(a)}}{\approx}\big{(}\mathbb{E}\big{\\{}\Delta\mathbf{S}_{k}^{H}(\mathbb{E}\\{{{\mathbf{S}_{k}}}{{\mathbf{S}_{k}}}^{H}\\})^{\dagger}\Delta\mathbf{S}_{k}\big{\\}}\big{)}^{-1}+\mathbf{I}_{N}$,
where $(\text{a})$ is because
$\mathbb{E}\big{\\{}\Delta\mathbf{S}_{k}^{H}(\mathbb{E}\\{{{\mathbf{S}_{k}}}{{\mathbf{S}_{k}}}^{H}\\})^{\dagger}\Delta\mathbf{S}_{k}\big{\\}}$
is a diagonal matrix and Proposition 3 is then satisfied. More specifically,
according to Lemma 2, the detailed expression of
$\mathbb{E}\big{\\{}\Delta\mathbf{S}_{k}^{H}(\mathbb{E}\\{{{\mathbf{S}_{k}}}{{\mathbf{S}_{k}}}^{H}\\})^{\dagger}\Delta\mathbf{S}_{k}\big{\\}}$
can be obtained as

$\begin{array}[]{c}\mathbb{E}\big{\\{}\Delta\mathbf{S}_{k}^{H}(\mathbb{E}\\{{{\mathbf{S}_{k}}}{{\mathbf{S}_{k}}}^{H}\\})^{\dagger}\Delta\mathbf{S}_{k}\big{\\}}=\frac{M\sigma^{2}}{p(N+1)l^{\text{UI}}_{1}l^{\text{IB}}}\left[\begin{array}[]{cccc}E_{Q,1}\mathbf{I}_{M}&\mathbf{0}_{M\times
M}&\cdots&\mathbf{0}_{M\times M}\\\ \mathbf{0}_{M\times
M}&E_{Q,2}\mathbf{I}_{M}&\cdots&\mathbf{0}_{M\times M}\\\
\vdots&\vdots&\ddots&\vdots\\\ \mathbf{0}_{M\times M}&\mathbf{0}_{M\times
M}&\cdots&E_{Q,\gamma}\mathbf{I}_{\delta}\end{array}\right],\end{array}$ (56)

where
$E_{Q,n}=\frac{\gamma^{2}+2\gamma+2}{M(\gamma+1)^{2}}-\frac{1}{\delta(\gamma+1)^{2}}$
if $n<\gamma$, and
$E_{Q,n}=-\frac{\gamma-1}{M(\gamma+1)^{2}}+\frac{\gamma^{2}+3\gamma}{\delta(\gamma+1)^{2}}$
otherwise. As a result, $\mathbb{E}\\{\tilde{\mathbf{S}}_{k}\\}$ (which is the
sum of two diagonal matrices) is also a diagonal matrix, thus
$\mathbb{E}\big{\\{}(\tilde{\mathbf{S}}_{k})^{-1}\big{\\}}\approx\big{(}\mathbb{E}\\{\tilde{\mathbf{S}}_{k}\\}\big{)}^{-1}$
holds based on Proposition 3. Accordingly, $\mathbf{P}_{\bm{\mu}_{k},2}^{b}$
can be approximated by
$\mathbf{P}_{\bm{\mu}_{k},2}^{b}\approx\frac{(\delta\gamma^{2}+2\gamma+2\delta-M)\sigma^{2}\sum_{n=1}^{(\gamma-1)M}\mathbb{E}\\{\mu_{k,n}\mu_{k,n}^{*}\\}}{p\delta(\gamma+1)^{2}(N+1)l^{\text{UI}}_{1}l^{\text{IB}}+(\delta\gamma^{2}+2\gamma+2\delta-M)\sigma^{2}}+\frac{(M\gamma^{2}+(3M-1)\gamma+\delta)\sigma^{2}\sum_{n=(\gamma-1)M+1}^{N}\mathbb{E}\\{\mu_{k,n}\mu_{k,n}^{*}\\}}{p\delta(\gamma+1)^{2}(N+1)l^{\text{UI}}_{1}l^{\text{IB}}+(M\gamma^{2}+(3M-1)\gamma+\delta)\sigma^{2}}$.
This thus completes the proof.

Appendix D Proof of Theorem 1

In the case of $M\geq N$, we first let
$\epsilon_{d,\text{2P}}^{\text{a}}=\epsilon_{d,1,\text{2P}}+\sum_{k=2}^{K}\epsilon_{d,k,\text{2P}}^{\text{a}}$
denote the overall MSE for estimating all the direct channels under the
proposed strategy, then we have
$\epsilon_{d,\text{3P}}-\epsilon_{d,\text{2P}}^{\text{a}}=\frac{M\sigma^{2}}{p}\big{(}\frac{1}{2}-\frac{1}{N+1}\big{)}\geq
0,$ which demonstrates the 2PCE strategy is better than the 3PCE strategy when
estimating the direct channels, and this performance gain holds regardless of
the numbers of BS antennas, reflecting elements and users. For the reflected
channel associated with the typical user, we have
$\epsilon_{r,1,\text{3P}}-\epsilon_{r,1,\text{2P}}=\frac{(1+K)M\sigma^{2}}{pK}-\frac{NM\sigma^{2}}{p(N+1)}>0$,
i.e., the proposed channel estimation strategy can also achieve a better
performance in this case. Besides, focusing on the estimation of
$\bm{\mu}_{k}$, it follows that
$\epsilon_{\bm{\mu}_{k},\text{3P}}^{\text{a}}-\epsilon_{\bm{\mu}_{k},\text{2P}}^{\text{a}}=A_{k,a}+B_{k,a}+C_{k,a}+D_{k,a}$,
where
$A_{k,a}=\frac{M\sigma^{2}\mathbb{E}\\{\mu_{k,1}\mu_{k,1}^{*}\\}}{pNl^{\text{BU}}_{1}l^{\text{BI}}\frac{NK}{K+N}+M\sigma^{2}}-\frac{M\sigma^{2}\mathbb{E}\\{\mu_{k,1}\mu_{k,1}^{*}\\}}{pNl^{\text{BU}}_{1}l^{\text{BI}}(N+1)+M\sigma^{2}}$,
$B_{k,a}=\frac{M\sigma^{2}\sum_{n=2}^{N}\mathbb{E}\\{\mu_{k,n}\mu_{k,n}^{*}\\}}{pN^{2}l^{\text{UI}}_{1}l^{\text{IB}}+M\sigma^{2}}-\frac{M\sigma^{2}\sum_{n=2}^{N}\mathbb{E}\\{\mu_{k,n}\mu_{k,n}^{*}\\}}{pN(N+1)l^{\text{UI}}_{1}l^{\text{IB}}+M\sigma^{2}}$,
$C_{k,a}=\frac{N\sigma^{2}}{pMl^{\text{UI}}_{1}l^{\text{IB}}\frac{NK}{(N-1)(K+1)}+M\sigma^{2}\frac{K}{(N-1)(K+1)}}-\frac{N\sigma^{2}}{2p(K-1)Ml^{\text{UI}}_{1}l^{\text{IB}}+2M\sigma^{2}\frac{K-1}{N+1}}$
and
$D_{k,a}=\frac{(1+K)N\sigma^{2}}{pMNKl^{\text{UI}}_{1}l^{\text{IB}}+(K+N)M\sigma^{2}}$.
As can be observed, $A_{k,a}$ is the difference of two fractions who have the
same numerator but different denominators, since the denominator of the first
fraction, i.e., $pNl^{\text{BU}}_{1}l^{\text{BI}}\frac{NK}{K+N}+M\sigma^{2}$,
is smaller than that of the second fraction, i.e.,
$pNl^{\text{BU}}_{1}l^{\text{BI}}(N+1)+M\sigma^{2}$, we obtain $A_{k,a}>0$.
Similarly, we have $B_{k,a}>0$ and $C_{k,a}>0$, together with the fact that
$D_{k,a}$ is also larger than 0, we obtain
$\epsilon_{\bm{\mu}_{k},\text{3P}}^{\text{a}}>\epsilon_{\bm{\mu}_{k},\text{2P}}^{\text{a}}$.
To summarize, we can conclude that in the case of $M\geq N$, the proposed 2PCE
strategy outperforms the 3PCE strategy when M is asymptotically large.

Appendix E Proof of Theorem 2

In the case of $M<N$, we first consider the estimation of all the direct
channels. Let
$\epsilon_{d,\text{2P}}^{b}=\epsilon_{d,1,\text{2P}}+\sum_{k=2}^{K}\epsilon_{d,k,\text{2P}}^{b}$
denote the overall MSE achieved by the proposed 2PCE strategy for estimating
the direct channels, then the corresponding MSE difference between the two
considered strategies is given by
$\epsilon_{d,\text{3P}}-\epsilon_{d,\text{2P}}^{b}=\frac{M\sigma^{2}}{p}\big{(}\frac{N}{N+1}-\frac{K-1}{\gamma+1}\big{)}.$
It can be observed that the sign of this MSE difference depends on the values
of $K$ and $\gamma$, i.e., if $K<\gamma+2$ is satisfied, the performance of
the 2PCE strategy is better, and vice versa. As for the estimation of the
scaling vectors, the MSE difference can be expressed as
$\epsilon_{\bm{\mu}_{k},\text{3P}}^{b}-\epsilon_{\bm{\mu}_{k},\text{2P}}^{b}=A_{k,b}+B_{k,b}+C_{k,b}+D_{k,b}+E_{k,b},$
where
$A_{k,b}=\frac{(1+K)N\sigma^{2}}{pKNMl^{\text{UI}}_{1}l^{\text{IB}}+(K+N)M\sigma^{2}}$,
$B_{k,b}=\frac{N\sigma^{2}}{pMl^{\text{UI}}_{1}l^{\text{IB}}\frac{NK}{(N-1)(K+1)}+M\sigma^{2}\frac{K}{(N-1)(K+1)}}-\frac{N\sigma^{2}}{2pMl^{\text{UI}}_{1}l^{\text{IB}}+\frac{2M\sigma^{2}}{N+1}}$,
$C_{k,b}=\frac{\sigma^{2}\mathbb{E}\\{\mu_{k,1}\mu_{k,1}^{*}\\}}{pl^{\text{UI}}_{1}l^{\text{IB}}\frac{KN}{K+N}+\sigma^{2}}-\frac{\sigma^{2}\mathbb{E}\\{\mu_{k,1}\mu_{k,1}^{*}\\}}{pl^{\text{UI}}_{1}l^{\text{IB}}\frac{(\delta\gamma^{2}+2\delta\gamma+\delta)(N+1)}{\delta\gamma^{2}+2\gamma+2\delta-M}+\sigma^{2}}$,
$D_{k,b}=\frac{\sigma^{2}\sum_{n=2}^{(\gamma-1)M}\mathbb{E}\\{\mu_{k,n}\mu_{k,n}^{*}\\}}{pNl^{\text{UI}}_{1}l^{\text{IB}}+\sigma^{2}}-\frac{\sigma^{2}\sum_{n=2}^{(\gamma-1)M}\mathbb{E}\\{\mu_{k,n}\mu_{k,n}^{*}\\}}{pl^{\text{UI}}_{1}l^{\text{IB}}\frac{(\delta\gamma^{2}+2\delta\gamma+\delta)(N+1)}{\delta\gamma^{2}+2\gamma+2\delta-M}+\sigma^{2}}$
and
$E_{k,b}=\frac{\sigma^{2}\sum_{n=(\gamma-1)M+1}^{N}\mathbb{E}\\{\mu_{k,n}\mu_{k,n}^{*}\\}}{pl^{\text{UI}}_{1}l^{\text{IB}}\frac{N}{M}+\sigma^{2}}-\frac{\sigma^{2}\sum_{n=(\gamma-1)M+1}^{N}\mathbb{E}\\{\mu_{k,n}\mu_{k,n}^{*}\\}}{pl^{\text{UI}}_{1}l^{\text{IB}}\frac{\delta\gamma^{2}+2\delta\gamma+\delta}{M\gamma^{2}+(3M-1)\gamma+\delta}+\sigma^{2}}$.
Notice that the numerators of the two fractions in $B_{k,b}$ are the same
while the denominator of the first fraction, i.e.,
$\frac{pNKMl^{\text{UI}}_{1}l^{\text{IB}}}{(N-1)(K+1)}+\frac{KM\sigma^{2}}{(N-1)(K+1)}$,
is smaller than that of the second fraction, i.e.,
$2pMl^{\text{UI}}_{1}l^{\text{IB}}+\frac{2M\sigma^{2}}{N+1}$, hence we have
$B_{k,b}>0$. Similarly, we can obtain $C_{k,b}>0$, $D_{k,b}>0$ and
$E_{k,b}>0$, therefore together with the fact that $A_{k,b}>0$, we have
$\epsilon_{\bm{\mu}_{k},\text{3P}}^{b}>\epsilon_{\bm{\mu}_{k},\text{2P}}^{b}$.
Besides, since $\epsilon_{r,1,\text{3P}}>\epsilon_{r,1,\text{2P}}$ holds
according to Appendix D, we can infer that the proposed 2PCE strategy can
outperform the 3PCE strategy when estimating the reflected channels. The proof
is thus completed.

## References

  * [1] F. Boccardi, R. W. Heath, A. Lozano, T. L. Marzetta, and P. Popovski, “Five disruptive technology directions for 5G,” _IEEE Commun. Mag._ , vol. 52, no. 2, pp. 74–80, Feb. 2014.
  * [2] M. Shafi, A. F. Molisch, P. J. Smith, T. Haustein, P. Zhu, P. D. Silva, F. Tufvesson, A. Benjebbour, and G. Wunder, “5G: A tutorial overview of standards, trials, challenges, deployment, and practice,” _IEEE J. Sel. Areas Commun._ , vol. 35, no. 6, pp. 1201–1221, Jun. 2017.
  * [3] Q. Wu and R. Zhang, “Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network,” _IEEE Commun. Mag._ , vol. 58, no. 1, pp. 106–112, Jun. 2020.
  * [4] E. Basar, M. D. Renzo, J. D. Rosny, M. Debbah, M. Alouini, and R. Zhang, “Wireless communications through reconfigurable intelligent surfaces,” _IEEE Access_ , vol. 7, pp. 116 753–116 773, Aug. 2019.
  * [5] Q. Wu, S. Zhang, B. Zheng, C. You, and R. Zhang, “Intelligent reflecting surface aided wireless communications: A tutorial,” _arXiv: 2007.02759v2_ , 2020.
  * [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, pp. 270–278.
  * [7] M. M. Zhao, Q. Wu, M. J. Zhao, and R. Zhang, “Exploiting amplitude control in intelligent reflecting surface aided wireless communication with imperfect CSI,” _arXiv: 2005.07002v2_ , 2020.
  * [8] H. Guo, Y. Liang, J. Chen, and E. G. Larsson, “Weighted sum-rate maximization for reconfigurable intelligent surface aided wireless networks,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 5, pp. 3064–3076, Feb. 2020.
  * [9] M. M. Zhao, A. Liu, Y. Wan, and R. Zhang, “Two-timescale beamforming optimization for intelligent reflecting surface aided multiuser communication with QoS constraints,” _arXiv:2011.02237v1_ , 2020.
  * [10] Q. Wu and R. Zhang, “Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 11, pp. 5394–5409, Nov. 2019.
  * [11] M. M. Zhao, Q. Wu, M. J. Zhao, and R. Zhang, “Intelligent reflecting surface enhanced wireless network: Two-timescale beamforming optimization,” _IEEE Trans. Wireless Commun._ , DOI: 10.1109/TWC.2020.3022297, 2020.
  * [12] D. Mishra and H. Johansson, “Channel estimation and low-complexity beamforming design for passive intelligent surface assisted MISO wireless energy transfer,” in _IEEE ICASSP_ , 2019, pp. 4659–4663.
  * [13] Y. Yang, B. Zheng, S. Zhang, and R. Zhang, “Intelligent reflecting surface meets OFDM: Protocol design and rate maximization,” _IEEE Trans. Commun._ , vol. 68, no. 7, pp. 4522–4535, Jul. 2020.
  * [14] T. L. Jensen and E. D. Carvalho, “An optimal channel estimation scheme for intelligent reflecting surfaces based on a minimum variance unbiased estimator,” in _IEEE ICASSP_ , May 2020, pp. 5000–5004.
  * [15] B. Zheng and R. Zhang, “Intelligent reflecting surface-enhanced OFDM: Channel estimation and reflection optimization,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 4, pp. 512–522, Apr. 2020.
  * [16] J. Lin, G. Wang, R. Fan, T. A. Tsiftsis, and C. Tellambura, “Channel estimation for wireless communication systems assisted by large intelligent surfaces,” _arXiv:1911.02158v1_ , 2019.
  * [17] C. You, B. Zheng, and R. Zhang, “Channel estimation and passive beamforming for intelligent reflecting surface: Discrete phase shift and progressive refinement,” _IEEE J. Select. Areas Commun._ , vol. 38, no. 11, pp. 2604–2620, Nov. 2020.
  * [18] Z. He and X. Yuan, “Cascaded channel estimation for large intelligent metasurface assisted massive MIMO,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 2, pp. 210–214, Oct. 2020.
  * [19] J. Chen, Y. Liang, H. V. Cheng, and W. Yu, “Channel estimation for reconfigurable intelligent surface aided multi-user MIMO systems,” _arXiv:1912.03619v1_ , 2019.
  * [20] Z. Wan, Z. Gao, and M. Alouini, “Broadband channel estimation for intelligent reflecting surface aided mmwave massive MIMO systems,” in _IEEE ICC_ , Jun. 2020, pp. 1–6.
  * [21] W. Zhang, J. Xu, W. Xu, D. W. K. Ng, and H. Sun, “Cascaded channel estimation for IRS-assisted mmwave multi-antenna with quantized beamforming,” _IEEE Commun. Lett._ , DOI: 10.1109/LCOMM.2020.3028878, 2020\.
  * [22] Z. Wang, L. Liu, and S. Cui, “Channel estimation for intelligent reflecting surface assisted multiuser communications: Framework, algorithms, and analysis,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 10, pp. 6607–6620, Jun. 2020.
  * [23] M. R. McKay and I. B. Collings, “General capacity bounds for spatially correlated Rician MIMO channels,” _IEEE Trans. Inf. Theory_ , vol. 51, no. 9, pp. 3121–3145, Sept. 2005.
  * [24] S. L. Loyka, “Channel capacity of MIMO architecture using the exponential correlation matrix,” _IEEE Commun. Lett._ , vol. 5, no. 9, pp. 369–371, Sept. 2001.
  * [25] Y. Wei, M. M. Zhao, M. J. Zhao, and M. Lei, “Learned conjugate gradient descent network for massive MIMO detection,” _IEEE Trans. Signal Process._ , vol. 68, pp. 6336–6349, Nov. 2020.
  * [26] J. Choi and D. J. Love, “Bounds on eigenvalues of a spatial correlation matrix,” _IEEE Commun. Lett._ , vol. 18, no. 8, pp. 1391–1394, Aug. 2014\.
  * [27] Q. Nadeem, H. Alwazani, A. Kammoun, A. Chaaban, M. Debbah, and M. Alouini, “Intelligent reflecting surface assisted multi-user MISO communication: Channel estimation and beamforming design,” _IEEE OJ-COMS_ , vol. 1, pp. 661–680, May 2020.
  * [28] R. L. Burden and J. D. Faires, _Numerical Analysis (5th ed.)_. Boston: Prindle, 1993.
  * [29] G. Livan, M. Novaes, and P. Vivo, _Introduction to Random Matrices Theory and Practice_. SpringerBriefs in Mathematical Physics, 2018.
  * [30] G. Cipolloni and L. Erdos, “Fluctuations for linear eigenvalue statistics of sample covariance random matrices,” _arXiv:1806.08751v2_ , 2018.
  * [31] T. Cai, C. Zhang, and H. Zhou, “Optimal rates of convergence for covariance matrix estimation,” _The Annals of Statistics_ , vol. 38, no. 4, pp. 2118–2144, 2010.
  * [32] D. Bernstein, _Matrix Mathematics_. UK: Princeton University Press, 2005.
  * [33] R. E. Hartwing, “Block generalized inverses,” _Arch. Rational Mech. Anal._ , no. 61, pp. 197–251, 1976.
  * [34] P. Babington, _Matrix Analysis and Applied Linear Algebra_. USA: Society for Industrial and Applied Mathematics, 2000\.