#  Power Measurement Enabled Channel Autocorrelation Matrix Estimation for
IRS-Assisted Wireless Communication

Ge Yan, Lipeng Zhu, Rui Zhang, Part of this work was presented at the IEEE
Global Communications Conference Workshops 2023, Kuala Lumpur, Malaysia [1].G.
Yan is with the NUS Graduate School, National University of Singapore,
Singapore 119077, and also with the Department of Electrical and Computer
Engineering, National University of Singapore, Singapore 117583 (e-mail:
geyan@u.nus.edu). L. Zhu is with the Department of Electrical and Computer
Engineering, National University of Singapore, Singapore 117583
(zhulp@nus.edu.sg). R. Zhang is with The Chinese University of Hong Kong,
Shenzhen, and Shenzhen Research Institute of Big Data, Shenzhen, China 518172
(e-mail: rzhang@cuhk.edu.cn). He is also with the Department of Electrical and
Computer Engineering, National University of Singapore, Singapore 117583
(e-mail: elezhang@nus.edu.sg).

###### Abstract

By reconfiguring wireless channels via passive signal reflection, intelligent
reflecting surface (IRS) can bring significant performance enhancement for
wireless communication systems. However, such performance improvement
generally relies on the knowledge of channel state information (CSI) for IRS-
involved links. Prior works on IRS CSI acquisition mainly estimate IRS-
cascaded channels based on the extra pilot signals received at the users/base
station (BS) with time-varying IRS reflections, which, however, needs to
modify the existing channel training/estimation protocols of wireless systems.
To address this issue, we propose in this paper a new channel estimation
scheme for IRS-assisted communication systems based on the received signal
power measured at the user terminal, which is practically attainable without
the need of changing the current protocol. Due to the lack of signal phase
information in measured power, the autocorrelation matrix of the BS-IRS-user
cascaded channel is estimated by solving an equivalent rank-minimization
problem. To this end, a low-rank-approaching (LRA) algorithm is proposed by
employing the fractional programming and alternating optimization techniques.
To reduce computational complexity, an approximate LRA (ALRA) algorithm is
also developed. Furthermore, these two algorithms are extended to be robust
against the receiver noise and quantization error in power measurement.
Simulation results are provided to verify the effectiveness of the proposed
channel estimation algorithms as well as the IRS passive reflection design
based on the estimated channel autocorrelation matrix.

###### Index Terms:

Intelligent reflecting surface (IRS), channel estimation, channel
autocorrelation matrix, passive reflection design.

##  I INTRODUCTION

In recent years, intelligent reflecting surface (IRS) has received great
attention due to its appealing capability of reconfiguring wireless channels.
By applying tunable phase shifts to incident wireless signals, IRS can
effectively control their propagation channels and thereby significantly
enhance the wireless communication performance, such as spectral/energy
efficiency and transmission reliability [2, 3, 4, 5]. Given such benefits as
well as its high deployment flexibility, low hardware cost, and low power
consumption, IRS has been identified as a key enabling technology for future
wireless networks such as 6G [2, 3]. However, to reap the high performance
gain by IRS, it is essential to acquire the channel state information (CSI)
for the IRS channels with its assisting base station (BS) and users, which is
practically difficult due to the following reasons. On one hand, the passive
IRS is not equipped with wireless transceivers at its reflecting elements,
making it impossible to estimate the BS-IRS and IRS-user channels separately.
Instead, only the cascaded BS-IRS-user/user-IRS-BS channel can be estimated at
the user/BS [6]. On the other hand, to compensate for the significant product-
distance path loss of the IRS-cascaded link, the number of IRS reflecting
elements needs to be sufficiently large in practice, e.g., several tens or
even hundreds [7, 8]. This results in high-dimensional IRS channel
vectors/matrices that incur prohibitive overhead and computational complexity
to estimate.

To tackle the above challenges, extensive studies have been devoted to the
cascaded channel estimation for IRS-assisted wireless communication systems,
while these works mainly adopted conventional pilot-based channel estimation
approaches by sending additional pilot signals with concurrent time-varying
IRS reflections [9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20]. For example,
by switching on only one reflecting element at one time, the IRS-cascaded
channel for each element was estimated based on the received pilot signals at
the user in [9, 10]. To exploit the array gain of IRS for channel estimation,
reflection codebooks for channel estimation were designed in [12, 13, 14] with
all reflecting elements switched on. In particular, the discrete Fourier
transform (DFT)-based IRS reflection codebook was shown to yield the minimum-
mean-square-error (MMSE) estimation of the cascaded channel in IRS-assisted
multiple-input single-output (MISO) systems [11] and IRS-enhanced orthogonal
frequency division multiplexing (OFDM) systems [12]. In [14], the DFT-based
codebook was employed for reflection training in an IRS-assisted multi-user
orthogonal frequency division multiple access (OFDMA) system. Furthermore,
under the practical discrete-phase-shift constraint on IRS reflection
coefficients, a Hadamard matrix-based IRS reflection pattern was proposed in
[13], while more sophisticated codebooks were designed in [7, 21, 22] to
achieve higher training efficiency. Besides, the authors in [12, 13, 14, 15]
proposed to group adjacent IRS reflecting elements into sub-surfaces so that
only the effective cascaded channel for each sub-surface needs to be
estimated, thus reducing the number of training pilots. In addition, various
compressed sensing algorithms were developed for IRS channel estimation by
utilizing the sparsity of the channel paths in the angular domain [16, 17, 18,
19]. Moreover, the deep residual network was employed in [20] to refine the
least-sqaure (LS) estimation of IRS channels.

In the aforementioned works, IRS-involved CSI is estimated based on the
received complex-valued pilot signals at the users/BS. However, in the
protocol of existing wireless communication systems, such as 4G/5G [23], the
pilot signals are dedicated to estimating the BS-user direct channels only. As
such, substantially additional pilots are required to estimate the new IRS-
cascaded CSI, which thus requires significant modifications of the existing
channel estimation/training protocols. To address this issue, IRS reflection
designs based on the received signal power measurement at the users have been
proposed [24, 25], which do not require additional pilot signals for explicit
IRS CSI estimation. As user power measurement is commonly adopted and easy to
obtain in existing wireless systems, such as reference signal received power
(RSRP) [26], this approach can be practically implemented without any change
of the current protocols. For example, for the conditional sample mean (CSM)
method proposed in [24], the received signal power at the user was modelled as
a random variable and a large number of IRS reflection patterns are randomly
generated for user power measurement. After the user measures the received
signal power values for all IRS reflection patterns, each IRS reflecting
element sets the reflection coefficient as the one that achieves the maximum
received power expectation conditioned on its value. The majority voting
algorithm proposed in [25] employed a similar idea but the received signal
strength indicator (RSSI) was used instead of RSRP. Both methods in [24, 25]
directly design the IRS reflection coefficients based on the power measurement
without estimating the channel. However, to obtain the optimal IRS reflection
performance, an excessively large number of IRS training reflections/power
measurement are generally required (in the quadratic order of the number of
IRS reflecting elements). Therefore, these methods may still be time-consuming
for practical implementation.

It is worth noting that the methods in [25, 24] did not fully exploit the
power measurement to obtain explicit CSI of IRS-cascaded channels, thus
resulting in their high overhead and low efficiency. To improve the existing
IRS channel estimation/reflection designs based on user power measurement,
this paper proposes a new channel autocorrelation matrix estimation scheme.
Specifically, since there is no signal phase information in received power
measurement, the IRS-cascaded channel vector cannot be completely recovered,
while its autocorrelation matrix can be (uniquely) estimated. By equivalently
transforming the channel autocorrelation estimation problem into a rank-
minimization problem, a low-rank-approaching (LRA) algorithm is proposed to
recover the channel autocorrelation matrix given the received signal power
measurement with time-varying random IRS reflections. In particular, the LRA
algorithm converts the rank-minimization problem into a fractional programming
problem and alternating optimization is applied to obtain its solution. To
reduce the computational complexity of the LRA algorithm, an approximate LRA
(ALRA) algorithm is also proposed, where the fractional programming is
approximated by low-complexity quadratic programming. Furthermore, the robust
extensions of both LRA and ALRA algorithms are designed by considering the
effects of receiver noise and quantization error in power measurement. The
convergence and estimation accuracy of the proposed algorithms are verified
via simulations. Besides, it is validated that with the IRS reflection designs
based on the estimated channel autocorrelation matrix, significantly higher
passive reflection gains can be achieved compared to other benchmark schemes,
which demonstrates the effectiveness of the proposed channel estimation
methods for IRS-assisted wireless communication systems based on user power
measurement.

Notations: Boldface letters refer to vectors (lower case) or matrices (upper
case). For square matrix $\bm{A}$, $\text{tr}(\bm{A})$ denotes its trace and
$\bm{A}^{-1}$ denotes its inverse matrix. For matrix $\bm{B}$, let
$\bm{B}^{T}$, $\bm{B}^{H}$, $\text{rank}(\bm{B})$, $\bm{B}^{\dagger}$,
$\|\bm{B}\|_{F}$, and $\text{vec}(\bm{B})$ denote the transpose, conjugate
transpose, rank, pseudo inverse, Frobenius norm, and vectorization of
$\bm{B}$, respectively. $\bm{I}_{N}$ denotes an $N\times N$ identity matrix.
$\bm{0}_{N\times M}$ denotes an $N\times M$ all-zero matrix. $\bm{x}^{T}$,
$\bm{x}^{H}$, and $\|\bm{x}\|_{2}$ denote the transpose, conjugate transpose,
and Euclidean norm of vector $\bm{x}$, respectively. Vector $\bm{1}_{K}$
denotes an all-one vector of size $K\times 1$. For vector $\bm{x}$,
$\text{diag}(\bm{x})$ denotes the diagonal matrix whose main diagonal elements
are extracted from $\bm{x}$. For matrix $\bm{A}$, $\text{diag}(\bm{A})$
denotes the vector whose elements are extracted from the main diagonal
elements of $\bm{A}$. For complex number $c$, let $\text{Re}(c)$,
$\text{Im}(c)$, and $|c|$ denote its real part, imaginary part, and magnitude,
respectively. $\mathbb{C}^{a\times b}$ and $\mathbb{R}^{a\times b}$ denote the
space of $a\times b$-dimensional complex and real matrices, respectively.
$\mathbb{E}[\cdot]$ denotes the statistical expectation. Symbol $j$ denotes
the imaginary unit $\sqrt{-1}$.

##  II System Model

Figure 1: An IRS-aided wireless communication system.

As shown in Fig. 1, we consider an IRS-aided downlink communication system
with a multi-antenna BS serving a single-antenna user111The solutions proposed
in this paper are also applicable to multiuser systems, by performing received
signal power measurement at different users simultaneously and employing the
proposed algorithms for each user’s channel estimation. , where an IRS is
deployed near the user to establish a reflected link to assist in their
communication. To focus on the IRS passive reflection design, we assume that
the BS fixes its precoding vector to guarantee the coverage of the user’s
residing area, and thus it can be considered having a single directional
antenna equivalently. The IRS is composed of $N_{irs}=N_{x}N_{z}$ reflecting
elements, where $N_{x}$ and $N_{z}$ are the number of reflecting elements in
the horizontal and vertical dimensions, respectively. Each element of the IRS
introduces a phase shift to the reflected signal. Let $u_{n}$ denote the
reflection coefficient of the $n$-th element, $n=1,2,\ldots,N_{irs}$, while
$\bm{u}=[u_{1},\ldots,u_{N_{irs}}]^{T}\in\mathbb{C}^{N_{irs}\times 1}$ and
$\bm{\Theta}=\text{diag}(\bm{u})\in\mathbb{C}^{N_{irs}\times N_{irs}}$ denote
the IRS reflection coefficient vector and matrix, respectively. Due to the
unit amplitude constraint on the reflecting coefficients, we set $|u_{n}|=1$,
$\forall n$. Furthermore, denote the number of bits for controlling the
discrete phase shift of each element as $b$. Then, the reflection coefficient
$u_{n}$ should be selected from a discrete set
$\Phi_{b}=\\{e^{j\Delta\theta},\ldots,e^{j2^{b}\Delta\theta}\\}$, with
$\Delta\theta=2\pi/2^{b}$. Denoting
$\Phi_{b}^{M}=\\{\bm{x}\in\mathbb{C}^{M\times
1}|x_{n}\in\Phi_{b},n=1,\ldots,M\\}$ as the set of $M$-dimensional vectors
whose elements are selected from $\Phi_{b}$, we thus have
$\bm{u}\in\Phi_{b}^{N_{irs}}$.

The baseband equivalent channels of the BS-IRS link, BS-user link, and IRS-
user link are denoted as $\bm{g}\in\mathbb{C}^{N_{irs}\times 1}$,
$h_{d}\in\mathbb{C}$, and $\bm{h}_{r}\in\mathbb{C}^{N_{irs}\times 1}$,
respectively. The received signal corrupted by noise at the user is given by

$y=\left(\bm{h}_{r}^{H}\bm{\Theta}\bm{g}+h_{d}^{*}\right)x+z,$ (1)

where $x$ is the transmitted signal with zero mean and power $p_{0}$ and
$z\sim\mathcal{CN}(0,\sigma^{2})$ is the independent and identically
distributed (i.i.d.) additive complex Guassian noise at the user receiver with
mean zero and power $\sigma^{2}$. Due to
$\bm{h}_{r}^{H}\bm{\Theta}=\bm{h}_{r}^{H}\text{diag}(\bm{u})=\bm{u}^{H}\text{diag}(\bm{h}_{r}^{H})$,
we can rewrite (1) as

$y=\left(\bm{u}^{H}\text{diag}(\bm{h}_{r}^{H})\bm{g}+h_{d}^{*}\right)x+z.$ (2)

Let $N=N_{irs}+1$ and define
$\bar{\bm{h}}=[\bm{g}^{H}\text{diag}(\bm{h}_{r}),h_{d}]^{H}\in\mathbb{C}^{N\times
1}$ as the equivalent channel and $\bm{v}=[\bm{u}^{T},1]^{T}\in\Phi_{b}^{N}$
as the equivalent IRS reflection vector. Then, the received noisy signal can
be simplified as $y=(\bm{v}^{H}\bar{\bm{h}})x+z$ and the power of the signal
(ignoring the noise for the time being) is given by

$p=\mathbb{E}\left[\left|\left(\bm{h}_{r}^{H}\bm{\Theta}\bm{g}+h_{d}^{*}\right)x\right|^{2}\right]=p_{0}\left|\bm{v}^{H}\bar{\bm{h}}\right|^{2}=p_{0}\text{tr}(\bar{\bm{H}}\bm{V}),$
(3)

where $\bar{\bm{H}}=\bar{\bm{h}}\bar{\bm{h}}^{H}$ and
$\bm{V}=\bm{v}\bm{v}^{H}$ are the autocorrelation matrices of the equivalent
channel $\bar{\bm{h}}$ and the equivalent IRS reflection vector $\bm{v}$,
respectively.

##  III Channel Autocorrelation Matrix Estimation

We assume that the user is quasi-static and its channels (or predominant
deterministic channel components) with the BS and IRS do not change within a
sufficiently long channel coherence time, during which the IRS first changes
its reflection coefficients $T_{p}$ times in total. In the meanwhile, the user
measures the corresponding received signal power values and feed them back to
a central processing unit (e.g., BS or IRS controller) that can design
accordingly and then set the IRS reflection coefficients for data
transmission. Specifically, for the $t$-th power measurement,
$t=1,\ldots,T_{p}$, a random IRS reflection vector $\bm{u}_{t}$ is applied and
$p_{t}$ denotes the received signal power at the user. Define
$\bm{v}_{t}=[\bm{u}_{t},1]^{T}$ and $\bm{V}_{t}=\bm{v}_{t}\bm{v}_{t}^{H}$, and
then we have $p_{t}=p_{0}\text{tr}(\bar{\bm{H}}\bm{V}_{t})$. With all the
received signals’ power obtained, the channel autocorrelation matrix is
estimated based on $p_{t}$ and $\bm{V}_{t}$, $t=1,\ldots,T_{p}$.

Given power measurement
$\bm{p}=[p_{1},\ldots,p_{T_{p}}]^{T}\in\mathbb{R}^{T_{p}\times 1}$, the
estimation problem can be formulated as finding a semidefinite rank-one matrix
$\bm{H}$ that satisfies $p_{0}\text{tr}(\bm{H}\bm{V}_{t})=p_{t}$,
$t=1,\ldots,T_{p}$, i.e.,

$\displaystyle\text{find}\ {\bm{H}}\in\mathbb{S}_{+}^{N}$ (4)
$\displaystyle~{}\mathrm{s.t.}\
p_{0}\text{tr}(\bm{H}\bm{V}_{t})=p_{t},t=1,\ldots,T_{p},$ (4a)
$\displaystyle~{}~{}~{}~{}~{}~{}\text{rank}(\bm{H})=1,$ (4b)

where $\mathbb{S}_{+}^{N}$ denotes the set of all positive semidefinite
hermitian matrices of dimension $N\times N$. The form of problem (4) is the
same as the PhaseLift problem studied in [27, 28], where the trace-
minimization relaxation was applied to find an approximate solution. However,
the performance of the approximate solution relies on the assumption that
vectors $\\{\bm{v}_{t},t=1,\ldots,T_{p}\\}$ are i.i.d. Guassian random
vectors, which is not applicable to the considered IRS reflection vector due
to its unit-amplitude elements with discrete phase shifts. Thus, we propose to
design more efficient methods customized to the practical constraint on
$\\{\bm{v}_{t},t=1,\ldots,T_{p}\\}$ to achieve better performance but with
lower computational complexity.

Due to the discrete phase shift constraint on IRS reflecting elements, the
solution for problem (4) may not be unique. Specifically, constraint (4a)
forms a system of linear equations and the $N\times N$ hermitian matrices form
a linear space of dimension $N^{2}$ [29]. Thus, the hermitian autocorrelation
matrix $\bm{H}$ can be uniquely determined if the maximum number of linearly
independent matrices in $\\{\bm{V}_{t},t=1,\ldots,T_{p}\\}$, denoted as
$D_{V}$ ($D_{V}\leq T_{p}$), is no less than $N^{2}$.  However, the unit-
amplitude and discrete phase of the entries in $\bm{V}_{t}$ result in
$\\{\bm{V}_{t},t=1,\ldots,T_{p}\\}$ always being confined in a subspace of
$\mathbb{C}^{N\times N}$ with its dimension smaller than $N^{2}$, as shown in
the following lemma.

######  Lemma 1.

For any $T_{p}$, $D_{V}\leq\mathcal{D}_{N}^{(b)}<N^{2}$ always holds, where
$\mathcal{D}_{N}^{(b)}=(N^{2}-N)/2+1$ for $b=1$ and
$\mathcal{D}_{N}^{(b)}=N^{2}-N+1$ for $b\geq 2$.

###### Proof:.

See Appendix A.  ∎

According to Lemma 1, constraint (4a) forms an underdetermined system of
equations due to the discrete phase shifts of IRS reflection vectors. As a
result, the uniqueness of the solution for problem (4) depends on the value of
$b$, as illustrated in the following proposition.

######  Proposition 1.

For $b=1$, problem (4) has only two solutions, i.e., $\bar{\bm{H}}$ and its
conjugate matrix $\bar{\bm{H}}^{*}$, if $N$ and $D_{V}$ are sufficiently large
(i.e., $N\geq 6$, $D_{V}=\mathcal{D}_{N}^{(1)}$). For $b\geq 2$, problem (4)
has one unique solution $\bar{\bm{H}}$ if $N$ and $D_{V}$ are sufficiently
large (i.e., $N\geq 3$, $D_{V}=\mathcal{D}_{N}^{(b)}$).

###### Proof:.

See Appendix B.  ∎

As the solutions for $b=1$ and $b\geq 2$ are different, we solve problem (4)
for these two cases separately. For the case of $b\geq 2$, we derive the
autocorrelation matrix $\bar{\bm{H}}$ by directly solving problem (4). For the
case of $b=1$, the autocorrelation matrix cannot be uniquely determined, while
the real part of the matrix, denoted as
$\bar{\bm{H}}_{r}=\text{Re}(\bar{\bm{H}})=\text{Re}(\bar{\bm{H}}^{*})$, is
unique according to Proposition 1. Since the equivalent IRS reflection vector
$\bm{v}$ is always a real vector for $b=1$, the received signal power at the
user always satisfies
$p=p_{0}\text{tr}(\bar{\bm{H}}\bm{V})=p_{0}\text{tr}(\bar{\bm{H}}_{r}\bm{V})$.
Therefore, we only need to estimate $\bar{\bm{H}}_{r}$ for optimizing the IRS
reflection vector for data transmission. Note that
$\bar{\bm{H}}_{r}=(\bar{\bm{H}}+\bar{\bm{H}}^{*})/2$ is semidefinite with its
rank no more than two because $\bar{\bm{H}}$ and $\bar{\bm{H}}^{*}$ are
semidefinite rank-one matrices. As such, for $b=1$, we consider to estimate
$\bar{\bm{H}}_{r}$ instead of $\bar{\bm{H}}$ by solving the following problem:

$\displaystyle\text{find}\ {\bm{H}_{r}}\in\mathbb{M}_{+}^{N}$ (5)
$\displaystyle~{}\mathrm{s.t.}\ p_{0}\text{tr}(\bm{H}_{r}\bm{V}_{t})=p_{t},\
t=1,\ldots,T_{p},$ (5a)
$\displaystyle~{}~{}~{}~{}~{}~{}\text{rank}(\bm{H}_{r})\leq 2,$ (5b)

where $\mathbb{M}_{+}^{N}$ denotes the set of all positive semidefinite real
symmetric matrices of dimension $N\times N$. The following proposition ensures
the existence and uniqueness of the solution for problem (5)222Note that
Propositions 1 and 2 only provide sufficient conditions on the uniqueness of
the solution for problem (4) for the two cases of $b=1$ and $b\geq 2$,
respectively. However, as revealed by simulations in Section VII, the unique
solution $\bar{\bm{H}}$ for $b\geq 2$ and $\bar{\bm{H}}_{r}$ for $b=1$ may
also be found for smaller values of $D_{V}$ (or $T_{p}$) than that given in
the two propositions. .

######  Proposition 2.

For $b=1$, problem (5) has one unique solution $\bar{\bm{H}}_{r}$ if $N$ and
$D_{V}$ are sufficiently large (i.e., $N\geq 6$,
$D_{V}=\mathcal{D}_{N}^{(1)}$).

###### Proof:.

See Appendix C.  ∎

In the following sections, the LRA and ALRA algorithms are proposed for $b=1$
and $b\geq 2$ by solving problems (4) and (5), respectively, based on received
signal power.

##  IV LRA Algorithm

In this section, the channel autocorrelation matrix estimation problems (4)
and (5) are solved by transforming them into equivalent rank-minimization
problems. Instead of directly finding a rank-one/two matrix solution, the
proposed LRA algorithm iteratively approaches a low-rank matrix solution via
alternating optimization.

###  IV-A Solution for Problem (4)

For $b\geq 2$, problem (4) can be equivalently transformed into the following
rank-minimization problem:

$\displaystyle\mathop{\min_{\bm{H}}}\ \text{rank}(\bm{H})$ (6)
$\displaystyle~{}\mathrm{s.t.}\ p_{0}\text{tr}(\bm{H}\bm{V}_{t})=p_{t},\
t=1,\ldots,T_{p},$ (6a)
$\displaystyle~{}~{}~{}~{}~{}~{}{\bm{H}}\in\mathbb{S}_{+}^{N}.$ (6b)

The equivalence between problems (4) and (6) is analyzed as follows. As the
channel autocorrelation matrix $\bar{\bm{H}}$ is feasible to problem (6), any
optimal solution for this problem, denoted by $\hat{\bm{H}}$, should satisfy
$\text{rank}(\hat{\bm{H}})\leq\text{rank}(\bar{\bm{H}})=1$, leading to
$\text{rank}(\hat{\bm{H}})=1$ and thus $\hat{\bm{H}}$ is also a solution for
problem (4). Reversely, any solution $\hat{\bm{H}}^{\prime}$ for problem (4)
is feasible to the rank-minimization problem (6) and satisfies
$\text{rank}(\hat{\bm{H}}^{\prime})=1$, which indicates that
$\hat{\bm{H}}^{\prime}$ is an optimal solution for problem (6).

Since $\bm{H}\in\mathbb{S}_{+}^{N}$ is nonzero, all the eigenvalues of
$\bm{H}$ are real and non-negative and $\text{tr}(\bm{H})>0$. Define the
eigenvalue-ratio function for matrix $\bm{H}$ as

$g(\bm{H})=\frac{\lambda_{1}(\bm{H})}{\text{tr}(\bm{H})},\
\bm{H}\in\mathbb{S}_{+}^{N},$ (7)

where $\lambda_{1}(\bm{H})$ is the largest eigenvalue of $\bm{H}$. Obviously,
$0<g(\bm{H})\leq 1$ holds for any nonzero $\bm{H}\in\mathbb{S}_{+}^{N}$, and
it is easy to verify that $\text{rank}(\bm{H})=1$ if and only if
$g(\bm{H})=1$. As we have mentioned above, any solution $\hat{\bm{H}}$ for
problem (6) satisfies $\text{rank}(\hat{\bm{H}})=1$. Thus, we have
$g(\hat{\bm{H}})=1$, which means that $\hat{\bm{H}}$ maximizes the eigenvalue-
ratio function $g(\bm{H})$. On the other hand, any matrix
$\hat{\bm{H}}^{\prime}$ that maximizes $g(\bm{H})$ subject to constraints (6a)
and (6b) also minimizes $\text{rank}(\bm{H})$. Therefore, the solutions for
problem (6) are the same as the solutions that maximize the eigenvalue-ratio
function $g(\bm{H})$ subject to constraints (6a) and (6b). Note that
$\lambda_{1}(\bm{H})=\mathop{\max_{\|\bm{x}\|_{2}\leq
1}}{\bm{x}^{H}\bm{H}\bm{x}}$. Thus, problem (6) can be written as

$\mathop{\max_{\bm{H}}\max_{\|\bm{x}\|_{2}\leq 1}}\
f(\bm{H},\bm{x})=\frac{\bm{x}^{H}\bm{H}\bm{x}}{\text{tr}(\bm{H})},~{}~{}\mathrm{s.t.}~{}\eqref{prob:cov-
est-rank-min-power},\eqref{prob:cov-est-rank-min-semidefinite}.$ (8)

This optimization problem is non-convex, while alternating optimization can be
employed to obtain a suboptimal solution for it. Given $\bm{H}$, an optimal
$\bm{x}$ can be obtained as the normalized eigenvector of $\bm{H}$
corresponding to the largest eigenvalue. Given $\bm{x}$, the optimization of
$\bm{H}$ is simplified as

$\mathop{\max_{\bm{H}}}\
\frac{\text{tr}(\bm{H}\bm{X})}{\text{tr}(\bm{H})},~{}~{}\mathrm{s.t.}~{}\eqref{prob:cov-
est-rank-min-power},\eqref{prob:cov-est-rank-min-semidefinite},$ (9)

with $\bm{X}=\bm{x}\bm{x}^{H}$. This is a fractional programming problem and
can be transformed into a convex optimization problem. Specifically, define
$\bm{G}=\bm{H}/\text{tr}(\bm{H})$ and $\gamma=1/\text{tr}(\bm{H})$. Then,
problem (9) can be transformed into

$\displaystyle\mathop{\max_{\bm{G},\gamma>0}}\ \text{tr}(\bm{G}\bm{X})$ (10)
$\displaystyle~{}\mathrm{s.t.}\ p_{0}\text{tr}(\bm{G}\bm{V}_{t})=p_{t}\gamma,\
t=1,\ldots,T_{p},$ (10a)
$\displaystyle~{}~{}~{}~{}~{}~{}{\bm{G}}\in\mathbb{S}_{+}^{N},\
\text{tr}(\bm{G})=1,$ (10b)

which is a convex semidefinite programming (SDP) problem and thus can be
solved by CVX [30]. Denoting the optimal solution for problem (10) as
$\hat{\bm{G}}$ and $\hat{\gamma}$, then the optimal solution for problem (9)
is obtained as $\hat{\bm{H}}=\hat{\bm{G}}/\hat{\gamma}$.

The proposed LRA algorithm to solve problem (4) for $b\geq 2$ is summarized in
Algorithm 1. The computational complexity of the algorithm is dominated by the
SDP problem (10), which is given by $\mathcal{O}(N^{4.5})$ [31]. Note that in
line $3$, only the eigenvector corresponding to the largest eigenvalue is
needed. The well-known power method [32] with the complexity of
$\mathcal{O}(N^{2})$ can be applied. Thus, the computational complexity of
Algorithm 1 is given by $\mathcal{O}(N^{4.5}I_{1})$, which is in the same
order as the trace-minimization method in [27], where $I_{1}$ is the total
number of iterations.

The convergence of the LRA algorithm is analyzed as follows. In Algorithm 1,
variables $\bm{H}$ and $\bm{x}$ are updated as $\bm{H}^{(i)}$ and
$\bm{x}^{(i)}$ in the $i$-th iteration. It is worth noting that $\bm{x}^{(i)}$
maximizes $f(\bm{H},\bm{x})$ given $\bm{H}=\bm{H}^{(i-1)}$ and thus

$g(\bm{H}^{(i-1)})=\max_{\|\bm{x}\|_{2}\leq
1}f(\bm{H}^{(i-1)},\bm{x})=f(\bm{H}^{(i-1)},\bm{x}^{(i)})$ (11)

holds. Moreover, $\bm{H}^{(i)}$ is the optimal solution for problem (9) with
$\bm{X}=\bm{X}^{(i)}$, which means

$\displaystyle f(\bm{H}^{(i)},\bm{x}^{(i)})$
$\displaystyle=\frac{\text{tr}(\bm{H}^{(i)}\bm{X}^{(i)})}{\text{tr}(\bm{H}^{(i)})}\geq\frac{\text{tr}(\bm{H}^{(i-1)}\bm{X}^{(i)})}{\text{tr}(\bm{H}^{(i-1)})}$
(12) $\displaystyle=f(\bm{H}^{(i-1)},\bm{x}^{(i)}),\ i=1,\ldots,I_{1}.$

Thus, we have

$\displaystyle g(\bm{H}^{(i)})$ $\displaystyle=\max_{\|\bm{x}\|_{2}\leq
1}f(\bm{H}^{(i)},\bm{x})\geq f(\bm{H}^{(i)},\bm{x}^{(i)})$ (13a)
$\displaystyle\geq f(\bm{H}^{(i-1)},\bm{x}^{(i)})=g(\bm{H}^{(i-1)}),\ \forall
i.$ (13b)

Therefore, the eigenvalue-ratio function $g(\bm{H}^{(i)})$ is non-decreasing
during the iterations for $i=0,1,\ldots,I_{1}$. Since $g(\bm{H})$ is upper-
bounded by $1$, the convergence of Algorithm 1 is guaranteed. Meanwhile, the
monotonic increase of the eigenvalue-ratio function indicates that matrix
$\bm{H}^{(i)}$ gradually approaches a rank-one matrix during the iterations.

Algorithm 1 LRA algorithm for $b\geq 2$.

0: $\\{\bm{v}_{t}$, $t=1,\ldots,T\\}$, $\bm{p}\in\mathbb{R}^{T\times 1}$,
threshold $\epsilon$.

1: Initialization: Solve $\bm{H}^{(0)}$ via trace-minimization relaxation
[27]; iteration index $i\leftarrow 1$.

2: while $g(\bm{H}^{(i-1)})\leq\epsilon$ do

3: Let $\bm{x}^{(i)}$ be the normalized eigenvector of $\bm{H}^{(i-1)}$
corresponding to the largest eigenvalue, and
$\bm{X}^{(i)}\leftarrow\bm{x}^{(i)}{\bm{x}^{(i)}}^{H}$.

4: Solve $\bm{G}^{(i)}$ and $\gamma^{(i)}$ from problem (10) with
$\bm{X}=\bm{X}^{(i)}$, and obtain
$\bm{H}^{(i)}\leftarrow\bm{G}^{(i)}/\gamma^{(i)}$.

5: $i\leftarrow i+1$.

6: end while

7: return The estimated matrix $\hat{\bm{H}}\leftarrow\bm{H}^{(i-1)}$.

###  IV-B Solution for Problem (5)

Next, we consider $b=1$. To solve matrix $\bar{\bm{H}}_{r}$ from problem (5),
a similar method to the case of $b\geq 2$ can be applied. Consider the
following rank-minimization problem:

$\displaystyle\mathop{\min_{\bm{H}_{r}}}\ \text{rank}(\bm{H}_{r})$ (14)
$\displaystyle~{}\mathrm{s.t.}\ p_{0}\text{tr}(\bm{H}_{r}\bm{V}_{t})=p_{t},\
t=1,\ldots,T_{p},$ (14a)
$\displaystyle~{}~{}~{}~{}~{}~{}{\bm{H}_{r}}\in\mathbb{M}_{+}^{N}.$ (14b)

Obviously, matrix $\bar{\bm{H}}_{r}=\text{Re}(\bar{\bm{H}})$ is feasible to
problem (14). Any optimal solution for problem (14), denoted by
$\hat{\bm{H}}_{r}$, satisfies
$\text{rank}(\hat{\bm{H}}_{r})\leq\text{rank}(\bar{\bm{H}}_{r})\leq 2$, which
means that $\hat{\bm{H}}_{r}$ is also a solution for problem (5). Thus, the
solution for problem (5) can be obtained by solving problem (14).

Define the generalized eigenvalue-ratio function as

$g_{r}(\bm{H}_{r})=\frac{\lambda_{1}(\bm{H}_{r})+\lambda_{2}(\bm{H}_{r})}{\text{tr}(\bm{H}_{r})},\
\bm{H}_{r}\in\mathbb{S}_{+}^{N},$ (15)

where $\lambda_{1}(\bm{H}_{r})$ and $\lambda_{2}(\bm{H}_{r})$ are the first
and second largest eigenvalues of $\bm{H}_{r}$, respectively. For nonzero
$\bm{H}_{r}\in\mathbb{S}_{+}^{N}$, we have $\text{tr}(\bm{H}_{r})>0$,
$0<g_{r}(\bm{H}_{r})\leq 1$, and $g_{r}(\bm{H}_{r})=1$ holds if and only if
$\text{rank}(\bm{H}_{r})\leq 2$. Therefore, solving problem (14) is equivalent
to maximizing $g_{r}(\bm{H}_{r})$ subject to constraints (14a) and (14b).
Furthermore, since we have

$\displaystyle\lambda_{1}(\bm{H}_{r})+$
$\displaystyle\lambda_{2}(\bm{H}_{r})=\mathop{\max_{\bm{x}_{1},\bm{x}_{2}}}\
\bm{x}_{1}^{T}\bm{H}_{r}\bm{x}_{1}+\bm{x}_{2}^{T}\bm{H}_{r}\bm{x}_{2},$ (16)
$\displaystyle~{}\mathrm{s.t.}\
\|\bm{x}_{1}\|_{2}=1,~{}\|\bm{x}_{2}\|_{2}=1,~{}\bm{x}_{1}^{T}\bm{x}_{2}=0,$
(16a)

problem (14) can be equivalently written as

$\displaystyle\mathop{\max_{\bm{H}_{r}}\max_{\bm{x}_{1},\bm{x}_{2}}}\
f_{r}(\bm{H}_{r},\bm{x}_{1},\bm{x}_{2})=\frac{\bm{x}_{1}^{T}\bm{H}_{r}\bm{x}_{1}+\bm{x}_{2}^{T}\bm{H}_{r}\bm{x}_{2}}{\text{tr}(\bm{H}_{r})}$
(17) $\displaystyle~{}\mathrm{s.t.}\ \eqref{prob:largest-two-eigenvalues-
unitary},\eqref{prob:cov-real-est-rank-min-power},\eqref{prob:cov-real-est-
rank-min-semidefinite}.$

Alternating optimization can be applied to solve problem (17) sub-optimally.
Given $\bm{H}_{r}$, the optimal $\bm{x}_{1}$ and $\bm{x}_{2}$ can be obtained
as the normalized eigenvectors corresponding to the first and second largest
eigenvalues of $\bm{H}_{r}$, respectively. Given $\bm{x}_{1}$ and
$\bm{x}_{2}$, $\bm{H}_{r}$ can be optimized via

$\mathop{\max_{\bm{H}_{r}}}\
\frac{\text{tr}(\bm{H}_{r}\bm{X})}{\text{tr}(\bm{H}_{r})},~{}~{}\mathrm{s.t.}~{}\eqref{prob:cov-
real-est-rank-min-power},\eqref{prob:cov-real-est-rank-min-semidefinite},$
(18)

with $\bm{X}=\bm{x}_{1}\bm{x}_{1}^{T}+\bm{x}_{2}\bm{x}_{2}^{T}$. This problem
is also a fractional programming problem and can be solved in the same way as
problem (9). The proposed LRA algorithm to solve problem (5) for the case of
$b=1$ is summarized in Algorithm 2, with a complexity of
$\mathcal{O}(N^{4.5}I_{2})$ which is also in the same order as the trace-
minimization method in [27], where $I_{2}$ is the total number of iterations.
The convergence can be guaranteed by the monotonic increase of the generalized
eigenvalue-ratio function $g_{r}(\bm{H}_{r}^{(i)})$, similar to Algorithm 1.

Algorithm 2 LRA algorithm for $b=1$.

0: $\\{\bm{v}_{t}$, $t=1,\ldots,T\\}$, $\bm{p}\in\mathbb{R}^{T\times 1}$,
threshold $\epsilon$.

1: Initialization: Solve $\bm{H}_{r}^{(0)}$ via trace-minimization relaxation
[27]; iteration index $i\leftarrow 1$.

2: while $g_{r}(\bm{H}_{r}^{(i-1)})\leq\epsilon$ do

3: Let $\bm{x}_{1}^{(i)}$ and $\bm{x}_{2}^{(i)}$ be the normlized eigenvectors
of $\bm{H}_{r}^{(i-1)}$ corresponding to the first and second largest
eigenvalues, and
$\bm{X}^{(i)}\leftarrow\bm{x}_{1}^{(i)}{\bm{x}_{1}^{(i)}}^{T}+\bm{x}_{2}^{(i)}{\bm{x}_{2}^{(i)}}^{T}$.

4: Solve $\bm{H}_{r}^{(i)}$ from problem (18) with $\bm{X}=\bm{X}^{(i)}$.

5: $i\leftarrow i+1$.

6: end while

7: return The estimated matrix $\hat{\bm{H}_{r}}\leftarrow\bm{H}_{r}^{(i-1)}$.

##  V Approximate LRA Algorithm

The proposed LRA algorithms can approach a low-rank matrix solution via
alternating maximization of the eigenvalue-ratio functions, but their
computational complexity may be high in practice because SDP is involved for
each iteration. In this section, proper approximations are implemented in the
optimization of $\bm{H}$ and $\bm{H}_{r}$ for $b\geq 2$ and $b=1$,
respectively, and the ALRA algorithm is proposed to reduce the computational
complexity. The basic principle of the LRA algorithms is retained, while the
eigenvalue-ratio functions are replaced by quadratic functions. By exploiting
the vector representations of hermitian matrices, closed-form solutions are
derived for the ALRA algorithm during the iterations and thus the
computational complexity is significantly reduced.

###  V-A Distance-Minimization Approximation

For each iteration of Algorithm 1 with $b\geq 2$, a normalized eigenvector
$\bm{x}$ is obtained and then $\bm{H}$ is optimized via problem (9). The
objective of problem (9) can be written as

$\mathop{\max_{\bm{H}}}\frac{\text{tr}(\bm{H}\bm{X})}{\text{tr}(\bm{H})}=\text{tr}(\bm{X})\mathop{\max_{\bm{H}}}\frac{\text{tr}(\bm{H}\bm{X})}{\text{tr}(\bm{H})\text{tr}(\bm{X})},$
(19)

where $\bm{X}=\bm{x}\bm{x}^{H}$. Note that
$\text{tr}(\bm{H}\bm{X})=\text{vec}(\bm{H})^{H}\text{vec}(\bm{X})$ and for
$\bm{H}\in\mathbb{S}_{+}^{N}$, we have
$\text{tr}(\bm{H})\geq\sqrt{\text{tr}(\bm{H}^{H}\bm{H})}=\|\text{vec}(\bm{H})\|_{2}$.
Thus, the objective function on the right-hand side of (19) is upper-bounded
by

$\displaystyle\frac{\text{tr}(\bm{H}\bm{X})}{\text{tr}(\bm{H})\text{tr}(\bm{X})}$
$\displaystyle\leq\frac{\text{tr}(\bm{H}\bm{X})}{\sqrt{\text{tr}(\bm{H}^{H}\bm{H})}\sqrt{\text{tr}(\bm{X}^{H}\bm{X})}}$
(20a)
$\displaystyle=\frac{\text{vec}(\bm{H})^{H}\text{vec}(\bm{X})}{\|\text{vec}(\bm{H})\|_{2}\|\text{vec}(\bm{X})\|_{2}},$
(20b)

which is the inner product of two normalized vectors
$\text{vec}(\bm{H})/\|\text{vec}(\bm{H})\|_{2}$ and
$\text{vec}(\bm{X})/\|\text{vec}(\bm{X})\|_{2}$. As a result, the optimization
of $\bm{H}$ can be approximately seen as aligning the direction of
$\text{vec}(\bm{H})$ to that of $\text{vec}(\bm{X})$, i.e., minimizing the
cosine distance between them, which is defined as

$d_{c}(\text{vec}(\bm{H}),\text{vec}(\bm{X}))=1-\frac{\text{vec}(\bm{H})^{H}\text{vec}(\bm{X})}{\|\text{vec}(\bm{H})\|_{2}\|\text{vec}(\bm{X})\|_{2}}.$
(21)

To reduce the computational complexity, we consider to minimize the Euclidean
distance between $\text{vec}(\bm{H})$ and $\mu\text{vec}(\bm{X})$, where $\mu$
is a scaling factor to be determined. Specifically, the distance-minimization
problem is given by

$\displaystyle\mathop{\min_{\bm{H},\mu}}\ d_{e}(\bm{H},\mu\bm{X})$ (22)
$\displaystyle~{}\mathrm{s.t.}\ p_{0}\text{tr}(\bm{H}\bm{V}_{t})=p_{t},\
t=1,\ldots,T_{p},$ (22a)
$\displaystyle~{}~{}~{}~{}~{}~{}\bm{H}\in\mathbb{S}_{+}^{N},$ (22b)

where $d_{e}(\bm{H},\mu\bm{X})$ denotes the Euclidean distance between vectors
$\text{vec}(\bm{H})$ and $\text{vec}(\mu\bm{X})$ and it is defined as

$d_{e}(\bm{H},\mu\bm{X})=\|\text{vec}(\bm{H})-\text{vec}(\mu\bm{X})\|_{2}=\|\bm{H}-\mu\bm{X}\|_{F}.$
(23)

To further simplify the optimization of $\bm{H}$, we relax the semidefinite
constraint (22b) to a hermitian constraint and solve the following problem:

$\displaystyle\mathop{\min_{\bm{H},\mu}}\ \|\bm{H}-\mu\bm{X}\|_{F}^{2}$ (24)
$\displaystyle~{}\mathrm{s.t.}\ p_{0}\text{tr}(\bm{H}\bm{V}_{t})=p_{t},\
t=1,\ldots,T_{p},$ (24a) $\displaystyle~{}~{}~{}~{}~{}~{}\bm{H}^{H}=\bm{H}.$
(24b)

Similarly, for Algorithm 2 with $b=1$, the optimization of $\bm{H}_{r}$ in
problem (18) can be approximated by minimizing the Euclidean distance
function. For $b=1$,
$\bm{X}=\bm{x}_{1}\bm{x}_{1}^{T}+\bm{x}_{2}\bm{x}_{2}^{T}$ is rank-two. To
improve the recovery accuracy, we consider a vector
$\bm{\mu}=[\mu_{1},\mu_{2}]^{T}\in\mathbb{R}^{2\times 1}$ such that
$\bm{H}_{r}$ is optimized via

$\displaystyle\mathop{\min_{\bm{H}_{r},\bm{\mu}}}\
\|\bm{H}_{r}-\left(\mu_{1}\bm{X}_{1}+\mu_{2}\bm{X}_{2}\right)\|_{F}^{2}$ (25)
$\displaystyle~{}\mathrm{s.t.}\ p_{0}\text{tr}(\bm{H}_{r}\bm{V}_{t})=p_{t},\
t=1,\ldots,T_{p},$ (25a)
$\displaystyle~{}~{}~{}~{}~{}~{}\bm{H}_{r}^{T}=\bm{H}_{r}\in\mathbb{R}^{N\times
N},$ (25b)

where $\bm{X}_{1}=\bm{x}_{1}\bm{x}_{1}^{T}$ and
$\bm{X}_{2}=\bm{x}_{2}\bm{x}_{2}^{T}$ are obtained in the previous iteration
in Algorithm 2. Since problems (24) and (25) are quadratic optimization
problems with linear constraints, their optimal solutions can be given in
closed form, which are derived in the following context.

###  V-B Solution for Problem (24)

According to [29], all hermitian matrices of size $N\times N$ form a linear
space of dimension $N^{2}$ with real combination coefficients. In particular,
each $N\times N$ hermitian matrix $\bm{A}$ can be represented by an
$N^{2}$-dimensional real vector $\bm{w}_{a}$, which is the coordinate of
$\bm{A}$ in the hermitian matrix space. Following the definition in [29],
denote the bijective mapping from $\bm{A}$ to $\bm{w}_{a}$ as function
$\mathcal{M}:\mathbb{C}^{N\times N}\to\mathbb{R}^{N^{2}\times 1}$, such that
$\bm{w}_{a}=\mathcal{M}(\bm{A})$ and $\bm{A}=\mathcal{M}^{-1}(\bm{w}_{a})$. It
has been shown in [29] that for any hermitian matrices $\bm{A}$ and $\bm{B}$,
we have $\textup{tr}(\bm{A}\bm{B})=\bm{w}_{a}^{T}\bm{w}_{b}$, where
$\bm{w}_{a}=\mathcal{M}(\bm{A})$ and $\bm{w}_{b}=\mathcal{M}(\bm{B})$. Based
on this property, problem (24) can be equivalently written as a quadratic
optimization problem with vector variables. Denote
$\bm{w}=\mathcal{M}(\bm{H})$, $\bm{w}_{X}=\mathcal{M}(\bm{X})$, and
$\bm{w}_{t}=\mathcal{M}(\bm{V}_{t})$, $t=1,\ldots,T_{p}$. Then, we have
$p_{t}=p_{0}\text{tr}(\bm{H}\bm{V}_{t})=p_{0}\bm{w}^{T}\bm{w}_{t}$, $\forall
t$, and

$\displaystyle\|\bm{H}-\mu\bm{X}\|_{F}^{2}$
$\displaystyle=\text{tr}\left((\bm{H}-\mu\bm{X})^{H}(\bm{H}-\mu\bm{X})\right)$
(26) $\displaystyle=\|\bm{w}-\mu\bm{w}_{X}\|_{2}^{2}.$

Define
$\bm{C}=p_{0}[\bm{w}_{1},\ldots,\bm{w}_{T_{p}}]\in\mathbb{R}^{N^{2}\times
T_{p}}$, and then the power constraints
$p_{t}=p_{0}\text{tr}(\bm{H}\bm{V}_{t})$ can be expressed as
$\bm{C}^{T}\bm{w}=\bm{p}$. Thus, problem (24) can be equivalently written as

$\mathop{\min_{\bm{w},\mu}}\
\|\bm{w}-\mu\bm{w}_{X}\|_{2}^{2},~{}~{}\mathrm{s.t.}\
\bm{C}^{T}\bm{w}=\bm{p},$ (27)

where the hermitian constraint (22b) is dropped because it is guaranteed by
the transformation from $\bm{H}$ to $\bm{w}$. For any given $\mu$, the
Lagrange function for variable $\bm{w}$ is given by

$L(\bm{w},\bm{\lambda})=\|\bm{w}-\mu\bm{w}_{X}\|_{2}^{2}+\bm{\lambda}^{T}(\bm{p}-\bm{C}^{T}\bm{w}),$
(28)

where $\bm{\lambda}\in\mathbb{R}^{T_{p}\times 1}$ is the Lagrange multiplier.
Then, $\bm{w}$ and $\bm{\lambda}$ can be solved via the Karush-Kuhn-Tucker
(KKT) conditions, i.e.,

$\displaystyle\frac{\partial L}{\partial\bm{\lambda}}$
$\displaystyle=\bm{p}-\bm{C}^{T}\bm{w}=\bm{0}_{T_{p}\times 1},$ (29a)
$\displaystyle\frac{\partial L}{\partial\bm{w}}$
$\displaystyle=2\bm{w}-2\mu\bm{w}_{X}-\bm{C}\bm{\lambda}=\bm{0}_{N^{2}\times
1}.$ (29b)

By assuming that $T_{p}\ll N^{2}$ and $\text{rank}(\bm{C})=T_{p}$ hold333In
practical systems, the number of power measurement $T_{p}$ is generally much
smaller than $N^{2}$ to avoid excessively large training overhead. Meanwhile,
we can properly design the IRS reflection coefficients for training such that
$\text{rank}(\bm{C})=D_{V}=T_{p}$ for $T_{p}\ll\mathcal{D}_{N}^{(b)}<N^{2}$. .
, these two equations can be solved as

$\displaystyle\bm{\lambda}$
$\displaystyle=2(\bm{C}^{T}\bm{C})^{-1}(\bm{p}-\mu\bm{C}^{T}\bm{w}_{X}),$
(30a) $\displaystyle\bm{w}$
$\displaystyle=\bm{D}\bm{p}+\mu(\bm{I}_{N^{2}}-\bm{D}\bm{C}^{T})\bm{w}_{X},$
(30b)

where $\bm{D}=\bm{C}(\bm{C}^{T}\bm{C})^{-1}$. Thus, the optimal $\mu$, denoted
as $\mu^{\star}$, can be solved by substituting (30b) to the objective
function in (27), i.e.,

$\mu^{\star}=\mathop{\arg\min_{\mu}}{\|\bm{D}\bm{p}-\mu\bm{D}\bm{C}^{T}\bm{w}_{X}\|_{2}^{2}}.$
(31)

If $\bm{D}\bm{C}^{T}\bm{w}_{X}\neq\bm{0}_{N^{2}\times 1}$, $\mu^{\star}$ is
obtained as

$\mu^{\star}=\frac{\bm{w}_{X}^{T}\bm{C}\bm{D}^{T}\bm{D}\bm{p}}{\|\bm{D}\bm{C}^{T}\bm{w}_{X}\|_{2}^{2}}=\frac{\bm{w}_{X}^{T}\bm{D}\bm{p}}{\bm{w}_{X}^{T}\bm{D}\bm{C}^{T}\bm{w}_{X}}.$
(32)

If $\bm{D}\bm{C}^{T}\bm{w}_{X}=\bm{0}_{N^{2}\times 1}$, $\mu^{\star}$ can be
any real value. For simplicity, we choose $\mu^{\star}=0$. Then, the optimal
solution $\bm{w}^{\star}$ can be obtained via (30b) with $\mu=\mu^{\star}$,
which yields the optimal solution for problem (24) as
$\bm{H}^{\star}=\mathcal{M}^{-1}(\bm{w}^{\star})$.

###  V-C Solution for Problem (25)

For $b=1$, $\bm{H}_{r}$ is a real symmetric matrix, which can also be
transformed into a real vector, similar to the case of $b\geq 2$. By defining
${\bm{w}_{r}}=\mathcal{M}(\bm{H}_{r})$,
${\bm{w}}_{X1}=\mathcal{M}(\bm{X_{1}})$,
${\bm{w}}_{X2}=\mathcal{M}(\bm{X}_{2})$ and
${\bm{w}_{rt}}=\mathcal{M}(\bm{V}_{t})$, $t=1,\ldots,T_{p}$, problem (25) can
be equivalently written as

$\mathop{\min_{{\bm{w}_{r}},\bm{\mu}}}\
\|{\bm{w}_{r}}-{\bm{W}}_{X}\bm{\mu}\|_{2}^{2},~{}~{}\mathrm{s.t.}\
{\bm{C}}^{T}{\bm{w}_{r}}=\bm{p},$ (33)

where
${\bm{C}}=p_{0}[{\bm{w}}_{r1},\ldots,{\bm{w}}_{rT_{p}}]\in\mathbb{R}^{N^{2}\times
T_{p}}$ and
${\bm{W}}_{X}=[{\bm{w}}_{X1},{\bm{w}}_{X2}]\in\mathbb{R}^{N^{2}\times 2}$.
Similar to the case of $b\geq 2$, for any given $\bm{\mu}$, the optimal vector
${\bm{w}_{r}}$ is given by

${\bm{w}_{r}}={\bm{D}}\bm{p}+(\bm{I}_{N^{2}}-{\bm{D}}{\bm{C}}^{T}){\bm{W}}_{X}\bm{\mu},$
(34)

where ${\bm{D}}={\bm{C}}({\bm{C}}^{T}{\bm{C}})^{-1}$. If
${\bm{D}}{\bm{C}}^{T}{\bm{W}}_{X}\neq\bm{0}_{N^{2}\times 2}$, the optimal
$\bm{\mu}$ for problem (33) is given by

$\bm{\mu}^{\star}=\left({\bm{W}}_{X}^{T}{\bm{D}}{\bm{C}}^{T}{\bm{W}}_{X}\right)^{\dagger}{\bm{W}}_{X}^{T}{\bm{D}}\bm{p}.$
(35)

If ${\bm{D}}{\bm{C}}^{T}{\bm{W}}_{X}=\bm{0}_{N^{2}\times 2}$,
$\bm{\mu}^{\star}=\bm{0}_{2\times 1}$ is chosen. Then, the optimal
${\bm{w}_{r}^{\star}}$ for problem (33) is obtained via (34) with
$\bm{\mu}=\bm{\mu}^{\star}$, and the optimal solution for problem (25) is
given by $\bm{H}_{r}^{\star}=\mathcal{M}^{-1}({\bm{w}_{r}^{\star}})$.

###  V-D Complexity and Convergence Analysis

For $b\geq 2$, the ALRA algorithm is implemented by replacing line $4$ in
Algorithm 1 with $\bm{H}^{\star}$ given by Section V-B. For $b=1$, the ALRA
algorithm is implemented by replacing line $4$ in Algorithm 2 with
$\bm{H}_{r}^{\star}$ given by Section V-C. It is worth noting that in the
original LRA algorithm, the trace-minimization method used for initialization
also requires SDP and may cause high computational complexity. Thus, we also
replace the trace-minimization initializations with solutions for problem (24)
and (25) as $\bm{X}=\bm{0}_{N^{2}\times N^{2}}$ and
$\bm{X}_{1}=\bm{X}_{2}=\bm{0}_{N^{2}\times N^{2}}$, respectively. The
computational complexity for the solutions in (30b),(32) and (34),(35) is
given by $\mathcal{O}(N^{2}T_{p})$. Meanwhile, the optimization of $\bm{X}$
can be implemented by the power method for both $b\geq 2$ and $b=1$ with the
complexity of $\mathcal{O}(N^{2})$. By denoting $I_{1}^{\prime}$ and
$I_{2}^{\prime}$ as the total number of iterations for $b\geq 2$ and $b=1$,
respectively, the total computational complexity of the ALRA algorithm is
given by $\mathcal{O}(N^{2}T_{p}I_{1}^{\prime})$ for $b\geq 2$ and
$\mathcal{O}(N^{2}T_{p}I_{2}^{\prime})$ for $b=1$, which is much lower than
those of the corresponding LRA algorithms and the trace-minimization method in
[27].

Moreover, the convergence of the ALRA algorithm is guaranteed, which is
analyzed as follows. For $b\geq 2$, $\bm{x}^{(i)}$, $\mu^{(i)}$, and
$\bm{H}^{(i)}$ are obtained for the $i$-th iteration. Specifically,
$\bm{x}^{(i)}$ is the eigenvector of $\bm{H}^{(i-1)}$ corresponding to the
largest eigenvalue, denoted as $\lambda_{1}^{(i-1)}$. Thus, matrix
$\lambda_{1}^{(i-1)}\bm{X}^{(i)}=\lambda_{1}^{(i-1)}\bm{x}^{(i)}{\bm{x}^{(i)}}^{H}$
is the rank-one matrix closest to $\bm{H}^{(i-1)}$, i.e.,

$\lambda_{1}^{(i-1)}\bm{X}^{(i)}=\mathop{\arg\min_{\text{rank}(\bm{A})\leq
1}}{\|\bm{H}^{(i-1)}-\bm{A}\|_{F}^{2}}.$ (36)

On the other hand, $\bm{H}^{(i)}$ and $\mu^{(i)}$ are jointly optimized to
minimize $\|\bm{H}-\mu\bm{X}^{(i)}\|_{F}^{2}$ under constraints (24a) and
(24b). Therefore, for $i\geq 1$, we have

$\displaystyle\left\|\bm{H}^{(i)}-\mu^{(i)}\bm{X}^{(i)}\right\|_{F}^{2}$
$\displaystyle\leq\left\|\bm{H}^{(i-1)}-\lambda_{1}^{(i-1)}\bm{X}^{(i)}\right\|_{F}^{2}$
(37a)
$\displaystyle\leq\left\|\bm{H}^{(i-1)}-\mu^{(i-1)}\bm{X}^{(i-1)}\right\|_{F}^{2},$
(37b)

which guarantees the convergence of the ALRA algorithm for $b\geq 2$.
Similarly, the ALRA algorithm for $b=1$ converges with a non-increasing
sequence
$\|\bm{H}_{r}^{(i)}-(\mu_{1}^{(i)}\bm{X}_{1}^{(i)}+\mu_{2}^{(i)}\bm{X}_{2}^{(i)})\|_{F}$,
indicating that the estimated channel autocorrelation matrix approaches a
rank-two matrix iteratively.

##  VI Robust Estimation with Practical Power measurement

In the above sections, the channel autocorrelation matrix was estimated based
on the received signal power. However, the exact value of the received signal
power cannot be obtained in practical communication systems. This is because
the received signals are corrupted by noise. Besides, the power measurement at
the user can only be fed back to the central processing unit approximately
after quantization, which thus incurs quantization error. In this section, the
robust extensions of the proposed LRA and ALRA algorithms are designed based
on practical power measurement, where the effects of receiver noise and
quantization error are considered.

###  VI-A Power Measurement

According to (1), the noisy signal power received at the user is given by

$p_{r}=\left|\left(\bm{v}^{H}\bar{\bm{h}}\right)x+z\right|^{2}=p_{0}\text{tr}(\bar{\bm{H}}\bm{V})+Z,$
(38)

where $Z=|z|^{2}+2\text{Re}((\bm{v}^{H}\bar{\bm{h}})xz^{*})$ is the noise
power satisfying $\mathbb{E}[Z]=\sigma^{2}$. By averaging $N_{0}$ noisy power
values, each practical power measurement is given by

$q=\frac{1}{N_{0}}\sum_{i=1}^{N_{0}}{\left|\left(\bm{v}^{H}\bar{\bm{h}}\right)x+z_{i}\right|^{2}}=p_{0}\text{tr}(\bar{\bm{H}}\bm{V})+\frac{1}{N_{0}}\sum_{i=1}^{N_{0}}{Z_{i}},$
(39)

where $z_{i}\sim\mathcal{CN}(0,\sigma^{2})$, $i=1,\ldots,N_{0}$, are i.i.d.
random noise variables and
$Z_{i}=|z_{i}|^{2}+2\text{Re}((\bm{v}^{H}\bar{\bm{h}})xz_{i}^{*})$. Define
$\bar{Z}=\frac{1}{N_{0}}\sum_{i=1}^{N_{0}}{Z_{i}}$ as the average noise power
in each power measurement. If $N_{0}$ is sufficiently large, $\bar{Z}$
converges to a Gaussian random variable with expectation
$\mathbb{E}[\bar{Z}]=\sigma^{2}$, according to the central limit theorem.

Despite the above model of practical power measurement, its value cannot be
exactly acquired because of the quantization. For example, in the current
protocol [33], the RSRP value is mapped into equally quantized levels within
the range from $-156$ dBm to $-44$ dBm. Thus, to investigate the effect of
power quantization mapping on the performance of channel autocorrelation
matrix estimation, we split the range $[-156,-44]$ dBm into $M$ equal levels
and define $D=112/M$ dB as the interval for each level. As such, the set of
quantization power levels is given by
$\Omega_{D}=\\{-44-(l+1/2)D~{}$dBm$,l=0,\ldots,M-1\\}$. For the $l$-th level,
the index $l$ is reported to the central processing unit if the power value
$q$ lies between $-44-(l+1)D$ dBm and $-44-lD$ dBm.

###  VI-B Reformulation of Channel Estimation Problem

Denote the actual power value for the $t$-th measurement as
$q_{t}=p_{0}\text{tr}(\bar{\bm{H}}\bm{V}_{t})+\bar{Z}_{t}$, where
$\bar{Z}_{t}$ is the noise power. Given the quantized power level index
$l_{t}$, we have $\zeta_{t}\leq q_{t}\leq\xi_{t}$, where
$\zeta_{t}=10^{-7.4-(l_{t}+1)D/10}$ and $\xi_{t}=10^{-7.4-l_{t}D/10}$ are the
lower and upper bounds of the $l_{t}$-th power level. For $b\geq 2$, the
channel autocorrelation matrix estimation problem is reformulated as

$\displaystyle\mathop{\min_{\bm{H},\bm{\delta}\geq\bm{0}_{T_{p}\times 1}}}\
{\|\bm{\delta}\|_{2}^{2}}$ (40) $\displaystyle~{}\mathrm{s.t.}\
\zeta_{t}-\delta_{t}\leq
p_{0}\text{tr}(\bm{H}\bm{V}_{t})+\sigma^{2}\leq\xi_{t}+\delta_{t},\forall t,$
(40a)
$\displaystyle~{}~{}~{}~{}~{}~{}\text{rank}(\bm{H})=1,~{}\bm{H}\in\mathbb{S}_{+}^{N},$
(40b)

where $\bm{\delta}=[\delta_{1},\ldots,\delta_{T_{p}\times
1}]\in\mathbb{R}^{T_{p}\times 1}$ is a slack vector which is introduced to
tackle the uncertainty of noise power $\bar{Z}_{t}$. For $b=1$, the problem is
given by

$\displaystyle\mathop{\min_{\bm{H}_{r},\bm{\delta}\geq\bm{0}_{T_{p}\times
1}}}\ {\|\bm{\delta}\|_{2}^{2}}$ (41) $\displaystyle~{}\mathrm{s.t.}\
\zeta_{t}-\delta_{t}\leq
p_{0}\text{tr}(\bm{H}_{r}\bm{V}_{t})+\sigma^{2}\leq\xi_{t}+\delta_{t},\forall
t,$ (41a) $\displaystyle~{}~{}~{}~{}~{}~{}\text{rank}(\bm{H}_{r})\leq
2,~{}\bm{H}_{r}\in\mathbb{M}_{+}^{N}.$ (41b)

It should be noticed that problems (40) and (41) are an non-convex
optimization problems because of the rank constraints, whose optimal solutions
cannot be solved by existing optimization tools in polynomial time. In the
following, we generalize the proposed LRA and ALRA algorithms to solve
problems (40) and (41) for achieving robust estimation of the channel
autocorrelation matrix.

Algorithm 3 Robust LRA algorithm for $b\geq 2$.

0: $\\{\bm{v}_{t}$, $t=1,\ldots,T\\}$, $\bm{l}\in\mathbb{R}^{T\times 1}$,
noise variance $\sigma^{2}$, penalty parameter $\rho$, and threshold
$\epsilon$.

1: Initialization: Compute $\bm{\zeta}$ and $\bm{\xi}$; solve $\bm{H}^{(0)}$
and $\bm{\delta}^{(0)}$ from problem (45); iteration index $i\leftarrow 1$.

2: while $g(\bm{H}^{(i-1)})\leq\epsilon$ do

3: Let $\bm{x}^{(i)}$ be the normalized eigenvector of $\bm{H}^{(i-1)}$
corresponding to the largest eigenvalue, and
$\bm{X}^{(i)}\leftarrow\bm{x}^{(i)}{\bm{x}^{(i)}}^{H}$.

4: Given $\bm{X}=\bm{X}^{(i)}$, solve $\bm{G}^{(i)}$ and $\gamma_{1}$ from
problem (44) with $\bm{\delta}=\bm{\delta}^{(i-1)}$ fixed, and then solve
$\bm{\delta}^{(i)}$ and $\gamma^{(i)}$ with $\bm{G}=\bm{G}^{(i)}$ fixed.
Obtain $\bm{H}^{(i)}\leftarrow\bm{G}^{(i)}/\gamma^{(i)}$.

5: $i\leftarrow i+1$.

6: end while

7: return The estimated matrix $\hat{\bm{H}}\leftarrow\bm{H}^{(i-1)}$.

###  VI-C Robust LRA Algorithm

In this subsection, the robust extension of the LRA algorithm is proposed. The
eigenvalue-ratio function defined in (7) and the generalized eigenvalue-ratio
function defined in (15) are added to the objective functions of problems (40)
and (41), respectively, as relaxations for the rank constraint. Specifically,
for $b\geq 2$, problem (40) is relaxed to

$\displaystyle\mathop{\max_{\bm{H},\bm{\delta}\geq\bm{0}_{T_{p}\times
1}}\max_{\|\bm{x}\|_{2}\leq 1}}\ f(\bm{H},\bm{x})-\rho\|\bm{\delta}\|_{2}^{2}$
(42) $\displaystyle~{}\mathrm{s.t.}\ \zeta_{t}-\delta_{t}\leq
p_{0}\text{tr}(\bm{H}\bm{V}_{t})+\sigma^{2}\leq\xi_{t}+\delta_{t},\ \forall
t,$ (42a) $\displaystyle~{}~{}~{}~{}~{}~{}{\bm{H}}\in\mathbb{S}_{+}^{N},$
(42b)

where $\rho>0$ is a predefined parameter. Compared to problem (8), problem
(42) approaches a low-rank solution and minimizes the penalty term
$\|\bm{\delta}\|_{2}^{2}$ simultaneously, thus leading to a robust estimation
of the channel autocorrelation matrix. The alternating optimization method can
be applied similar to that for solving problem (8), while the optimization for
$\bm{H}$ should be modified to decrease $\|\bm{\delta}\|_{2}^{2}$. Given
matrix $\bm{X}=\bm{x}\bm{x}^{H}$, problem (42) is simplified as

$\mathop{\max_{\bm{H},\bm{\delta}\geq\bm{0}_{T_{p}\times 1}}}\
\frac{\text{tr}(\bm{H}\bm{X})}{\text{tr}(\bm{H})}-\rho\|\bm{\delta}\|_{2}^{2},~{}~{}\mathrm{s.t.}\
\eqref{prob:cov-est-maxmax-robust-multiary-power},\eqref{prob:cov-est-maxmax-
robust-multiary-semidefinite}.$ (43)

With $\bm{G}=\bm{H}/\text{tr}(\bm{H})$ and $\gamma=1/\text{tr}(\bm{H})$,
problem (43) can be transformed into

$\displaystyle\mathop{\max_{\bm{G},\bm{\delta}\geq\bm{0}_{T_{p}\times
1},\gamma>0}}\ \text{tr}(\bm{G}\bm{X})-\rho\|\bm{\delta}\|_{2}^{2}$ (44)
$\displaystyle~{}\mathrm{s.t.}\ \zeta_{t}-\delta_{t}\leq
p_{0}\text{tr}(\bm{G}\bm{V}_{t})/\gamma+\sigma^{2}\leq\xi_{t}+\delta_{t},\
\forall t,$ (44a)
$\displaystyle~{}~{}~{}~{}~{}~{}{\bm{G}}\in\mathbb{S}_{+}^{N},\
\text{tr}(\bm{G})=1.$ (44b)

Although problem (44) is nonconvex, it can be transformed into a convex one
with respect to any two of the variables $\bm{G},\bm{\delta}$, and $\gamma$,
with the third variable given and fixed. To obtain an approximate solution, we
first solve $\bm{G}$ and $\gamma$ in problem (44) with $\bm{\delta}$ being
fixed to maximize the eigenvalue-ratio function, and then solve $\gamma$ and
$\bm{\delta}$ with $\bm{G}$ being fixed to minimize the penalty term
$\|\bm{\delta}\|_{2}^{2}$. For initialization, the robust trace-minimization
problem is employed as

$\mathop{\min_{\bm{H},\bm{\delta}\geq\bm{0}_{T_{p}\times 1}}}\
\text{tr}(\bm{H})+\rho\|\bm{\delta}\|_{2}^{2},~{}~{}\mathrm{s.t.}\
\eqref{prob:cov-est-maxmax-robust-multiary-power},\eqref{prob:cov-est-maxmax-
robust-multiary-semidefinite}.$ (45)

The robust LRA algorithm for $b\geq 2$ is summarized in Algorithm 3. It is
easy to verify that its computational complexity is
$\mathcal{O}(N^{4.5}I_{3})$, which is the same as that of Algorithm 1, and the
convergence of the algorithm is guaranteed, where $I_{3}$ is the total number
of iterations for Algorithm 3. For $i\geq 1$, $(\bm{G}^{(i)},\gamma_{1})$ is
solved from problem (44) with $\bm{\delta}=\bm{\delta}^{(i-1)}$ being fixed,
i.e., $\text{tr}(\bm{G}\bm{X}^{(i)})$ is maximized subject to constraints
(44a) and (44b) with $\bm{\delta}=\bm{\delta}^{(i-1)}$. Meanwhile, since
$(\bm{G}^{(i-1)},\gamma^{(i-1)},\bm{\delta}^{(i-1)})$ is feasible to problem
(44), we have
$\text{tr}(\bm{G}^{(i)}\bm{X}^{(i)})\geq\text{tr}(\bm{G}^{(i-1)}\bm{X}^{(i)})$.
Additionally, it is easy to verify that
$f(\bm{H}^{(i)},\bm{x})=\text{tr}(\bm{G}^{(i)}\bm{X})$ for any $\bm{x}$ and
$\bm{X}=\bm{x}\bm{x}^{H}$. Therefore, we have

$\displaystyle g(\bm{H}^{(i)})$ $\displaystyle\geq
f(\bm{H}^{(i)},\bm{x}^{(i)})=\text{tr}(\bm{G}^{(i)}\bm{X}^{(i)})$ (46a)
$\displaystyle\geq\text{tr}(\bm{G}^{(i-1)}\bm{X}^{(i)})=f(\bm{H}^{(i-1)},\bm{x}^{(i)})$
(46b) $\displaystyle=g(\bm{H}^{(i-1)}),~{}i=1,\ldots,I_{3},$ (46c)

which guarantees the convergence of Algorithm 3.

Similarly, for $b=1$, problem (41) is relaxed to

$\displaystyle\mathop{\max_{\bm{H}_{r},\bm{\delta}\geq\bm{0}_{T_{p}\times
1}}\max_{\bm{x}_{1},\bm{x}_{2}}}\
f_{r}(\bm{H}_{r},\bm{x}_{1},\bm{x}_{2})-\rho\|\bm{\delta}\|_{2}^{2}$ (47)
$\displaystyle~{}\mathrm{s.t.}\
\|\bm{x}_{1}\|_{2}=1,~{}\|\bm{x}_{2}\|_{2}=1,~{}\bm{x}_{1}^{T}\bm{x}_{2}=0,$
(47a) $\displaystyle~{}~{}~{}~{}~{}~{}\zeta_{t}-\delta_{t}\leq
p_{0}\text{tr}(\bm{H}_{r}\bm{V}_{t})+\sigma^{2}\leq\xi_{t}+\delta_{t},\
\forall t,$ (47b)
$\displaystyle~{}~{}~{}~{}~{}~{}{\bm{H}_{r}}\in\mathbb{M}_{+}^{N}.$ (47c)

By applying alternating optimization, the optimization of $\bm{H}_{r}$ given
matrix $\bm{X}=\bm{x}_{1}\bm{x}_{1}^{T}+\bm{x}_{2}\bm{x}_{2}^{T}$ is expressed
as

$\mathop{\max_{\bm{H}_{r},\bm{\delta}\geq\bm{0}_{T_{p}\times 1}}}\
\frac{\text{tr}(\bm{H}_{r}\bm{X})}{\text{tr}(\bm{H}_{r})}-\rho\|\bm{\delta}\|_{2}^{2},~{}~{}\mathrm{s.t.}\
\eqref{prob:cov-est-maxmax-robust-binary-power},~{}\eqref{prob:cov-est-maxmax-
robust-binary-semidefinite}.$ (48)

Problem (48) has the same form as problem (43) and thus it can be solved
similarly. For initialization, the following problem is considered:

$\mathop{\max_{\bm{H}_{r},\bm{\delta}\geq\bm{0}_{T_{p}\times 1}}}\
\text{tr}(\bm{H}_{r})+\rho\|\bm{\delta}\|_{2}^{2},~{}~{}\mathrm{s.t.}\
\eqref{prob:cov-est-maxmax-robust-binary-power},~{}\eqref{prob:cov-est-maxmax-
robust-binary-semidefinite}.$ (49)

The robust LRA algorithm for $b=1$ is summarized in Algorithm 4. By denoting
the total number of iterations as $I_{4}$, its computational complexity can be
obtained as $\mathcal{O}(N^{4.5}I_{4})$, and the convergence of Algorithm 4 is
ensured by the non-decreasing property of the sequence
$g_{r}(\bm{H}_{r}^{(i)})$, $i=1,\ldots,I_{4}$.

Algorithm 4 Robust LRA algorithm for $b=1$.

0: $\\{\bm{v}_{t}$, $t=1,\ldots,T\\}$, $\bm{l}\in\mathbb{R}^{T\times 1}$,
noise variance $\sigma^{2}$, penalty parameter $\rho$, and threshold
$\epsilon$.

1: Initialization: Compute $\bm{\zeta}$ and $\bm{\xi}$; solve
$\bm{H}_{r}^{(0)}$ and $\bm{\delta}^{(0)}$ from problem (49); iteration index
$i\leftarrow 1$.

2: while $g(\bm{H}^{(i-1)})\leq\epsilon$ do

3: Let $\bm{x}_{1}^{(i)}$ and $\bm{x}_{2}^{(i)}$ be the normlized eigenvectors
of $\bm{H}_{r}^{(i-1)}$ corresponding to the first and second largest
eigenvalues, and
$\bm{X}^{(i)}\leftarrow\bm{x}_{1}^{(i)}{\bm{x}_{1}^{(i)}}^{T}+\bm{x}_{2}^{(i)}{\bm{x}_{2}^{(i)}}^{T}$.

4: Given $\bm{X}=\bm{X}^{(i)}$, solve $\bm{H}_{r}^{(i)}$ from problem (48).

5: $i\leftarrow i+1$.

6: end while

7: return The estimated matrix $\hat{\bm{H}}_{r}\leftarrow\bm{H}_{r}^{(i-1)}$.

###  VI-D Robust ALRA Algorithm

In this subsection, the robust ALRA algorithm is proposed with a lower
computational complexity. Similar to that in Section V, the Euclidean distance
functions are employed as approximations of the eigenvalue-ratio functions of
the robust LRA algorithm. Specifically, given $\bm{X}$ for $b\geq 2$, the
optimization of $\bm{H}$ (i.e., problem (43)) is approximated by

$\displaystyle\mathop{\min_{\bm{H},\bm{\delta}\geq\bm{0}_{T_{p}\times
1},\mu}}\ \|\bm{H}-\mu\bm{X}\|_{F}^{2}+\rho\|\bm{\delta}\|_{2}^{2}$ (50)
$\displaystyle~{}\mathrm{s.t.}\ \zeta_{t}-\delta_{t}\leq
p_{0}\text{tr}(\bm{H}\bm{V}_{t})+\sigma^{2}\leq\xi_{t}+\delta_{t},\forall t,$
(50a) $\displaystyle~{}~{}~{}~{}~{}~{}\bm{H}^{H}=\bm{H}.$ (50b)

Note that problem (50) is a quadratic minimization problem with inequality
constraints, which is convex but the optimal solution cannot be obtained in
closed form. To derive an explicit solution, we define
$\hat{q}_{t}=(\zeta_{t}+\xi_{t})/2$, and approximate constraint (50a) with
$\hat{q}_{t}-\delta_{t}\leq
p_{0}\text{tr}(\bm{H}\bm{V}_{t})+\sigma^{2}\leq\hat{q}_{t}+\delta_{t}$,
$\forall t$. Then, problem (50) can be further approximated by

$\displaystyle\mathop{\min_{\bm{H},\mu}}\
\|\bm{H}-\mu\bm{X}\|_{F}^{2}+\rho\sum_{t=1}^{T_{p}}{\left|p_{0}\text{tr}(\bm{H}\bm{V}_{t})+\sigma^{2}-\hat{q}_{t}\right|^{2}}$
(51) $\displaystyle~{}\mathrm{s.t.}\ \bm{H}^{H}=\bm{H},$ (51a)

which is a quadratic minimization problem with linear constraints. By
leveraging the transformation $\mathcal{M}$, problem (51) can be equivalently
written as

$\mathop{\min_{\bm{w},\mu}}\
\|\bm{w}-\mu\bm{w}_{X}\|_{2}^{2}+\rho\|\bm{C}^{T}\bm{w}-(\hat{\bm{q}}-\sigma^{2}\bm{1}_{T_{p}})\|_{2}^{2},$
(52)

with $\hat{\bm{q}}=[\hat{q}_{1},\ldots,\hat{q}_{T_{p}}]^{T}$. By taking the
first derivative of the objective function, the optimal solutions can be
obtained in closed form as

$\displaystyle\mu^{\star}$
$\displaystyle=\frac{\rho}{\|\bm{\psi}\|_{2}^{2}}\bm{\psi}^{T}\bm{\Upsilon}\bm{C}(\hat{\bm{q}}-\sigma^{2}\bm{1}_{T_{p}}),$
(53a) $\displaystyle\bm{w}^{\star}$
$\displaystyle=\rho\bm{\Upsilon}\bm{C}(\hat{\bm{q}}-\sigma^{2}\bm{1}_{T_{p}})+\mu^{\star}\bm{\Upsilon}\bm{w}_{X},$
(53b) $\displaystyle\bm{\psi}$
$\displaystyle=\left(\bm{I}_{N^{2}}-\bm{\Upsilon}\right)\bm{w}_{X},$ (53c)
$\displaystyle\bm{\Upsilon}$
$\displaystyle=\left(\bm{I}_{N^{2}}+\rho\bm{C}\bm{C}^{T}\right)^{-1},$ (53d)

and thus the optimal solution for problem (51) is given by
$\bm{H}^{\star}=\mathcal{M}^{-1}(\bm{w}^{\star})$. Therefore, the robust ALRA
algorithm for $b\geq 2$ can be implemented by replacing line $4$ in Algorithm
3 with the optimal solution for problem (51) given $\bm{X}$ from the previous
iteration, i.e., $\bm{H}^{\star}=\mathcal{M}^{-1}(\bm{w}^{\star})$, where
$\bm{w}^{\star}$ is given by (53). Besides, the optimal solution
$\bm{H}^{\star}$ for problem (51) given $\bm{X}=\bm{0}_{N^{2}\times N^{2}}$ is
employed as the initialization.

Similarly, for $b=1$, the robust estimation of $\bm{H}_{r}$ can be realized by
converting problem (25) to the following problem:

$\displaystyle\begin{aligned} \mathop{\min_{\bm{H}_{r},\bm{\mu}}}\
&\|\bm{H}_{r}-(\mu_{1}\bm{X}_{1}+\mu_{2}\bm{X}_{2})\|_{F}^{2}\\\
&~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\rho\sum_{t=1}^{T_{p}}{\left|p_{0}\text{tr}(\bm{H}_{r}\bm{V}_{t})+\sigma^{2}-\hat{q}_{t}\right|^{2}}\end{aligned}$
(54) $\displaystyle~{}\mathrm{s.t.}\ \bm{H}_{r}^{T}=\bm{H}_{r}.$ (54a)

According to the derivations for solving problem (51), the optimal solution
for problem (54) can be obtained as
$\bm{H}_{r}^{\star}=\mathcal{M}^{-1}(\bm{w}_{r}^{\star})$, where
$\bm{w}_{r}^{\star}$ is solved from

$\mathop{\min_{{\bm{w}_{r}},\bm{\mu}}}\
\|{\bm{w}_{r}}-{\bm{W}}_{X}\bm{\mu}\|_{2}^{2}+\rho\|{\bm{C}}^{T}{\bm{w}_{r}}-(\hat{\bm{q}}-\sigma^{2}\bm{1}_{T_{p}})\|_{2}^{2}.$
(55)

Specifically, the optimal solutions ${\bm{w}}_{r}^{\star}$ and
$\bm{\mu}^{\star}$ are given by

$\displaystyle\bm{\mu}^{\star}$
$\displaystyle=\rho\left({\bm{\Psi}}^{T}{\bm{\Psi}}\right)^{\dagger}{\bm{\Psi}}^{T}{\bm{\Upsilon}}{\bm{C}}(\hat{\bm{q}}-\sigma^{2}\bm{1}_{T_{p}}),$
(56a) $\displaystyle{\bm{w}}_{r}^{\star}$
$\displaystyle=\rho{\bm{\Upsilon}}{\bm{C}}(\hat{\bm{q}}-\sigma^{2}\bm{1}_{T_{p}})+{\bm{\Upsilon}}{\bm{W}}_{X}\bm{\mu}^{\star},$
(56b) $\displaystyle{\bm{\Psi}}$
$\displaystyle=\left(\bm{I}_{N^{2}}-{\bm{\Upsilon}}\right){\bm{W}}_{X},$ (56c)
$\displaystyle{\bm{\Upsilon}}$
$\displaystyle=\left(\bm{I}_{N^{2}}+\rho{\bm{C}}{\bm{C}}^{T}\right)^{-1}.$
(56d)

Thus, the robust ALRA algorithm for $b=1$ is implemented by replacing line $4$
in Algorithm 4 with the optimal solution for problem (54) given $\bm{X}_{1}$
and $\bm{X}_{2}$ from the previous iteration, i.e.,
$\bm{H}_{r}^{\star}=\mathcal{M}^{-1}(\bm{w}_{r}^{\star})$, where
$\bm{w}_{r}^{\star}$ is given by (56). The initialization is obtained as
$\bm{H}_{r}^{\star}$ with $\bm{X}_{1}=\bm{X}_{2}=\bm{0}_{N^{2}\times N^{2}}$.

Note that $\bm{\Upsilon}\in\mathbb{R}^{N^{2}\times N^{2}}$ can be calculated
offline since $\rho$ and $\bm{C}$ are known. Thus, the online complexity of
the robust estimation for the channel autocorrelation matrix is dominated by
that of the multiplication by $\bm{\Upsilon}$ to vectors, which is
$\mathcal{O}(N^{4})$. To further reduce the complexity, we apply the Woodbury
matrix identity [34] to $\bm{\Upsilon}$ such that we have

$\displaystyle\bm{\Upsilon}$
$\displaystyle=\left(\bm{I}_{N^{2}}+\rho\bm{C}\bm{C}^{T}\right)^{-1}$ (57a)
$\displaystyle=\bm{I}_{N^{2}}-\bm{C}\left({\rho}^{-1}\bm{I}_{T_{p}^{2}}+\bm{C}^{T}\bm{C}\right)^{-1}\bm{C}^{T}.$
(57b)

Given $\bm{C}\in\mathbb{R}^{N^{2}\times T_{p}}$, matrix
$\bm{S}=({\rho}^{-1}\bm{I}_{T_{p}^{2}}+\bm{C}^{T}\bm{C})^{-1}\in\mathbb{R}^{T_{p}\times
T_{p}}$ can be calculated offline. Then, the multiplication to a vector by
$\bm{\Upsilon}=\bm{I}_{N^{2}}-\bm{C}\bm{S}\bm{C}^{T}$ can be implemented with
complexity $\mathcal{O}(N^{2}T_{p})$. Such implementation can be applied to
both cases of $b\geq 2$ and $b=1$. Therefore, the computational complexity of
the robust ALRA algorithm is $\mathcal{O}(N^{2}T_{p}I_{3}^{\prime})$ for
$b\geq 2$ and $\mathcal{O}(N^{2}T_{p}I_{4}^{\prime})$ for $b=1$, where
$I_{3}^{\prime}$ and $I_{4}^{\prime}$ are the total number of iterations for
the robust ALRA algorithm for $b\geq 2$ and $b=1$, respectively. The
complexity is thus significantly reduced compared to the robust LRA algorithms
and trace-minimization method [27].

Next, the convergence of the proposed robust ALRA algorithm is analyzed. For
$b\geq 2$, it can be shown that the objective function of problem (51), which
is denoted as $\varphi(\bm{H},\mu,\bm{X})$ and can be seen as a penaltized
version of the squared distance function $d_{e}^{2}(\bm{H},\mu\bm{X})$, is
non-increasing during the iterations. In the $i$-th iteration for $b\geq 2$,
$\bm{x}^{(i)}$ is the eigenvector of $\bm{H}^{(i-1)}$ corresponding to the
largest eigenvalue $\lambda_{1}^{(i-1)}$, and thus (36) still holds, i.e.,
$\lambda_{1}^{(i-1)}\bm{X}^{(i)}$ is the rank-one matrix closest to
$\bm{H}^{(i-1)}$. Moreover, $\bm{H}^{(i)}$ and $\mu^{(i)}$ are solved
optimally from problem (51) given $\bm{X}^{(i)}$. Thus, it is easy to verify
that

$\displaystyle\varphi(\bm{H}^{(i)},\mu^{(i)},\bm{X}^{(i)})$
$\displaystyle\leq\varphi(\bm{H}^{(i-1)},\lambda_{1}^{(i-1)},\bm{X}^{(i)})$
(58a) $\displaystyle\leq\varphi(\bm{H}^{(i-1)},\mu^{(i-1)},\bm{X}^{(i-1)})$
(58b)

always holds for $i=1,\cdots,I_{3}^{\prime}$. Since
$\varphi(\bm{H},\mu,\bm{X})\geq 0$, the robust ALRA algorithm converges for
$b\geq 2$. The similar analysis applies to the case of $b=1$, where the
objective function of problem (54), denoted as
$\varphi_{r}(\bm{H}_{r},\bm{\mu},\bm{X}_{1},\bm{X}_{2})$, is also non-
increasing during the iterations.

##  VII Performance Evaluation

###  VII-A Simulation Setup

In the simulation, the BS and the IRS are located at $(50,-150,20)$ and
$(-2,-1,0)$ in meters (m) in a three-dimensional coordinate system,
respectively. The location of the user is randomly generated as
$(x_{u},y_{u},0)$ with $0\leq x_{u},y_{u}\leq 10$. The size of IRS is set as
$N_{x}\times N_{z}=8\times 8=64$, and thus we have $N=64+1=65$. The path loss
model for all channels is given by $\eta=C_{0}d^{-\alpha}$, where $d$ is the
signal propagation distance, while $C_{0}$ and $\alpha$ are the channel gain
at the reference distance of $1$ m and the path-loss exponent, which are
denoted for the BS-user, BS-IRS, and IRS-user links as $C_{0,BU}$ and
$\alpha_{BU}$, $C_{0,BI}$ and $\alpha_{BI}$, and $C_{0,IU}$ and $\alpha_{IU}$,
respectively. The corresponding path loss coefficients are denoted as
$\eta_{BU}$, $\eta_{BI}$ and $\eta_{IU}$, respectively. For all channels,
Rician fading is assumed with the Rician factor denoted as $\beta_{BU}$,
$\beta_{BI}$ and $\beta_{IU}$ for the BS-user, BS-IRS and IRS-user links,
respectively. Specifically, the expression of the BS-IRS channel vector
$\bm{g}$ is given below as an example:

$\bm{g}=\sqrt{\frac{\beta_{BI}}{1+\beta_{BI}}}\bm{g}^{\text{LoS}}+\sqrt{\frac{1}{1+\beta_{BI}}}\bm{g}^{\text{NLoS}},$
(59)

where the vector $\bm{g}^{\text{NLoS}}=\sqrt{\eta_{BI}}\tilde{\bm{g}}$ is the
Gaussian non-line-of-sight (NLoS) component with
$\tilde{\bm{g}}\sim\mathcal{CN}(\bm{0}_{N_{irs}\times 1},\bm{I}_{N_{irs}})$
and $\bm{g}^{\text{LoS}}$ is the deterministic line-of-sight (LoS) component
given by

$\bm{g}^{\text{LoS}}=\sqrt{\eta_{BI}}\bm{b}_{N}(\omega,\psi).$ (60)

Vector $\bm{b}_{N}(\omega,\psi)\in\mathbb{C}^{N\times 1}$ is the steering
vector of the LoS path from the BS to the IRS, where $\omega\in[0,\pi)$ and
$\psi\in[0,\pi)$ are the physical azimuth and elevation angles of arrival
(AoAs) at the IRS, respectively. Specifically, define
$\bm{a}_{N}(\phi)=[e^{j\pi
0\phi},\ldots,e^{j\pi(N-1)\phi}]^{T}\in\mathbb{C}^{N\times 1}$ as the
$N$-dimensional steering vector. Then, $\bm{b}_{N}(\omega,\psi)$ is defined as
$\bm{b}_{N}(\omega,\psi)=\bm{a}_{N_{x}}(\cos{(\omega)}\sin{(\psi)})\otimes\bm{a}_{N_{z}}(\cos{(\psi)})$,
where $\otimes$ is the Kronecker product. In addition, we set $C_{0,BU}=-33$
dB, $C_{0,BI}=C_{0,IU}=-30$ dB, $\alpha_{BU}=3.7$,
$\alpha_{BI}=\alpha_{IU}=2$, $\beta_{BU}=0$, $\beta_{BI}=10$, $\beta_{IU}=1$
and $p_{0}=30$ dBm, $\sigma^{2}=-90$ dBm, $\epsilon=0.95$ and $\rho=10$,
unless specified otherwise. The IRS reflection vectors
$\\{\bm{v}_{t},t=1,\ldots,T_{p}\\}$ are randomly generated subject to the
discrete phase shift constraint.  Each point in simulation results is averaged
over $1000$ random user locations and channel realizations.

(a) LRA convergence.

(b) ALRA convergence.

Figure 2: Convergence of the LRA and ALRA algorithms assuming perfect power
measurement ($T_{p}=65$).

###  VII-B Performance with Perfect Power Measurement

In this subsection, the performance of both LRA and ALRA algorithms is
evaluated assuming perfect received signal power measurement as considered in
Sections IV and V.

####  VII-B1 Algorithm Convergence

As discussed in Sections IV and V, the convergence of the proposed algorithms
is guaranteed with the alternating optimization process. In particular, for
the LRA algorithm, the eigenvalue-ratio function and the generalized
eigenvalue-ratio function, i.e., $g(\bm{H}^{(i)})$ in (7) for $b\geq 2$ and
$g_{r}(\bm{H}_{r}^{(i)})$ in (15) for $b=1$, are non-decreasing with $i$ and
upper-bounded by $1$ for the LRA algorithm. On the other hand, the Euclidean
distance functions for the ALRA algorithm, i.e.,
$d_{e}(\bm{H}^{(i)},\mu^{(i)}\bm{X}^{(i)})$ for $b\geq 2$ and
$d_{e}(\bm{H}_{r}^{(i)},\mu_{1}^{(i)}\bm{X}_{1}^{(i)}+\mu_{2}^{(i)}\bm{X}_{2}^{(i)})$
for $b=1$ as defined in (23), are positive and non-increasing with $i$. In
Figs. 2(a) and 2(b), the values of eigenvalue-ratio functions and the
normalized distance functions over the iterations for both $b=1$ and $b=2$ are
shown for LRA and ALRA algorithms, respectively. The total number of power
measurement is set as $T_{p}=65$. It is observed from Fig. 2(a) that the
eigenvalue-ratio function values start from small values with the
initializations, and converge to $1$ within $6$ iterations, while the ALRA
algorithm converges after $50$ iterations, as is shown in Fig. 2(b). This
verifies that the estimated channel autocorrelation matrices obtained by the
proposed LRA and ALRA algorithms iteratively approach a rank-one matrix for
$b=2$ and a rank-two matrix for $b=1$, respectively.

Figure 3: Normalized error by different channel estimation schemes assuming
perfect power measurement.

(a) $b=1$.

(b) $b=2$.

Figure 4: Effective channel gain with different channel estimation schemes
assuming perfect power measurement.

####  VII-B2 Channel Autocorrelation Matrix Estimation Error

For the case of $b\geq 2$, the channel autocorrelation matrix $\bar{\bm{H}}$
is estimated as $\hat{\bm{H}}$, and the normalized estimation error is defined
as
$\mathcal{E}_{b}=\|\hat{\bm{H}}-\bar{\bm{H}}\|_{F}^{2}/\|\bar{\bm{H}}\|_{F}^{2}$.
For the case of $b=1$, the matrix $\bar{\bm{H}}_{r}$ is estimated as
$\hat{\bm{H}}_{r}$, and thus the normalized estimation error is defined as
$\mathcal{E}_{b}=\|\hat{\bm{H}}_{r}-\bar{\bm{H}}_{r}\|_{F}^{2}/\|\bar{\bm{H}}_{r}\|_{F}^{2}$.

Fig. 3 shows the normalized estimation error under different numbers of power
measurement, i.e., $T_{p}$, for both LRA and ALRA algorithms, where “TM”
represents the error of the estimated channel autocorrelation matrix by
employing the trace-minimization relaxation method given in [27]. It is
observed that the estimation error of the LRA algorithms decreases rapidly
with $T_{p}$ and is much smaller than that of the TM benchmark. Moreover, the
normalized estimation error vanishes for LRA and TM algorithms when $T_{p}\geq
185$ and $320$, respectively, which means that the unique solution for the
channel autocorrelation matrix is recovered successfully and the proposed LRA
algorithm outperforms TM significantly in terms of estimation accuracy. On the
other hand, the normalized estimation error of the ALRA algorithm decreases
slowly with $T_{p}$, which is caused by the distance-minimization
approximation, while its performance is good for sufficiently large $T_{p}$.
Additionally, it is shown that the estimation error for $b=1$ is smaller than
that for $b=2$ for all algorithms. This is because only the real part of
$\bar{\bm{H}}$ needs to be estimated for $b=1$.

####  VII-B3 IRS Reflection Design with Estimated Channel

After the channel autocorrelation matrix is estimated, the IRS reflection
vector $\bm{v}$ can be optimized to maximize the effective channel gain
between the BS and user, denoted as
$\bar{\gamma}=\text{tr}(\bar{\bm{H}}\bm{V})$, for data transmission. For
$b\geq 2$, we apply the eigenvalue decomposition to the estimated matrix
$\hat{\bm{H}}$ and define $\hat{\lambda}_{1}$ as the largest eigenvalue of
$\hat{\bm{H}}$ and $\hat{\bm{x}}_{1}$ as the corresponding normalized
eigenvector. Since $\hat{\bm{H}}\approx\bar{\bm{H}}$ is nearly rank-one, the
effective channel gain can be approximated as
$\bar{\gamma}\approx\text{tr}(\hat{\bm{H}}\bm{V})/p_{0}\approx\hat{\lambda}_{1}|\hat{\bm{x}}_{1}^{H}\bm{v}|^{2}/p_{0}$.
Then, the IRS beamforming vector $\bm{v}$ is optimized to maximize
$|\hat{\bm{x}}_{1}^{H}\bm{v}|^{2}$ subject to the discrete phase shift
constraint $\bm{v}\in\Phi_{b}^{N}$, which can be solved optimally by using the
method proposed in [35]. For $b=1$, similarly, eigenvalue decomposition is
applied to the estimated matrix $\hat{\bm{H}}_{r}$, where $\hat{\lambda}_{1}$
and $\hat{\lambda}_{2}$ denote the first and second largest eigenvalues of
$\hat{\bm{H}}_{r}$, with $\hat{\bm{x}}_{1}$ and $\hat{\bm{x}}_{2}$ denoting
the corresponding eigenvectors, respectively. Since the IRS beamforming vector
$\bm{v}$ is always a real vector for $b=1$, the effective channel gain can be
approximated by

$\displaystyle\bar{\gamma}\approx\frac{1}{p_{0}}\text{tr}(\hat{\bm{H}}_{r}\bm{V})$
$\displaystyle\approx\frac{1}{p_{0}}\left|\hat{\lambda}_{1}^{\frac{1}{2}}\hat{\bm{x}}_{1}^{T}\bm{v}\right|^{2}+\frac{1}{p_{0}}\left|\hat{\lambda}_{2}^{\frac{1}{2}}\hat{\bm{x}}_{2}^{T}\bm{v}\right|^{2}$
(61)
$\displaystyle=\frac{1}{p_{0}}\left|\left(\hat{\lambda}_{1}^{\frac{1}{2}}\hat{\bm{x}}_{1}+j\hat{\lambda}_{2}^{\frac{1}{2}}\hat{\bm{x}}_{2}\right)^{H}\bm{v}\right|^{2},$

and the optimal vector $\bm{v}$ can be obtained according to the method
proposed in [35].

(a) R-LRA convergence.

(b) R-ALRA convergence.

Figure 5: Convergence of robust LRA and robust ALRA algorithms based on noisy
power measurement for $b=1$ ($\sigma^{2}=-85$ dBm and $N_{0}=1$).

(a) Normalized estimation error versus $T_{p}$ for different values of $N_{0}$
and $\sigma^{2}$ ($b=1$).

(b) Effective channel gain versus $T_{p}$ for different values of $N_{0}$ and
$\sigma^{2}$ ($b=1$).

Figure 6: Performance of different schemes versus $T_{p}$ with various noise
levels for $b=1$.

For comparison, three benchmark schemes for IRS reflection design based on
power measurement are listed as follows: i) TM (trace-minimization): The IRS
reflection vector is optimized as illustrated above based on the channel
autocorrelation matrix estimated via the trace-minimization method in [27].
ii) UB (upper bound): The upper bound on the effective channel gain is
obtained by using the optimal IRS reflection vector based on the perfect CSI
$\bm{h}$; iii) RMS (random-max sampling): A large number of random IRS
reflection vectors are applied with $u_{n}$ uniformly distributed in
$\Phi_{b}$, $\forall n$, and the one achieving the largest received signal
power is used; iv) CSM (conditional sample mean): This is the method proposed
in [24], where a large number of random IRS reflection vectors are applied,
and the empirical expectation of the received signal power is calculated
conditioned on $u_{n}$ fixed at every possible value, $\forall n$. For each
element, CSM selects the phase shift that maximizes the empirical expectation
conditioned on its value.

The effective channel gain obtained based on the estimated channels by the
proposed algorithms as well as other benchmark schemes for $b=1$ and $b=2$ are
shown in Fig. 4(a) and 4(b), respectively. Due to the discrete phase shift for
IRS, the upper bound for the effective channel gain for $b=2$ is higher than
that for $b=1$. As can be observed, the effective channel gain achieved by the
proposed LRA scheme increases rapidly with the number of power measurement,
and reaches the upper bound when $T\geq 180$ for both $b=1$ and $b\geq 2$
cases. In comparison, the ALRA scheme achieve lower effective channel gain and
approach the upper bound when $T_{p}\geq 350$. Similar to the case of
estimation error, the effective channel gain achieved by TM also falls in
between those of LRA and ALRA, while reaching the upper bound when $T_{p}\geq
270$. Moreover, both LRA and ALRA schemes outperform RMS and CSM schemes
significantly for almost all values of $T_{p}$. These results validate the
effectiveness of the proposed channel estimation algorithms in improving the
effective channel gain between the BS and user with optimized IRS reflections.

###  VII-C Performance with Noisy Power Measurement

In this subsection, the robust LRA and robust ALRA algorithms proposed in
Section VI are evaluated based on noisy power measurement, which are labeled
as R-LRA and R-ALRA, respectively. The impact of noise received at the user on
the performance is analyzed while the quantization effect is ignored for the
time being. Due to space limitation, we only consider the case of $b=1$ in
this subsection and the effective channel gain for TM is omitted for brevity.

####  VII-C1 Algorithm Convergence

To verify the convergence of the R-LRA and R-ALRA algorithms based on noisy
power measurement, Fig. 5(a) shows the values of the generalized eigenvalue-
ratio function $g_{r}(\bm{H}_{r}^{(i)})$ defined in (15), while Fig. 5(b)
shows the normalized values of the penaltized distance function
$\varphi_{r}(\bm{H}_{r}^{(i)},\bm{\mu}^{(i)},\bm{X}_{1}^{(i)},\bm{X}_{2}^{(i)})$
defined in Section VI-D during the iterations. We set $\sigma^{2}=-85$ dBm and
$N_{0}=1$. It is observed that the convergence of the proposed robust
algorithms is still guaranteed. However, the generalized eigenvalue-ratio
function of R-LRA may not converge to $1$ due to the impact of random noise.
Especially for large $T_{p}$, it becomes increasingly difficult to find a low-
rank matrix $\bm{H}_{r}$ that can keep the penalty term small enough. The
penaltized distance function for R-ALRA also increases with $T_{p}$. This
indicates that both the proposed R-LRA and R-ALRA algorithms approach low-rank
solutions and reduce the penalty terms at the same time, yielding robust
estimations for channel autocorrelation matrices.

####  VII-C2 Estimation Error and Effective Channel Gain

The normalized estimation error and effective channel gain achieved by the
proposed and benchmark schemes for $b=1$ are shown versus the number of power
measurement in Figs. 6(a) and 6(b), respectively. For three subfigures in Fig.
6(a), the number of reference signals for each power measurement, $N_{0}$, is
set to be $1$, $8$ and $80$, respectively, and the estimation errors with
$\sigma^{2}=-90$ dBm and $\sigma^{2}=-85$ dBm are shown in each subfigure for
comparison. Their corresponding effective channel gains are shown in the three
subfigures in Fig. 6(b).

In Fig. 6(a), it is observed that the normalized estimation error decreases
with $T_{p}$ and $N_{0}$ but increases with $\sigma^{2}$ for all schemes, and
both the R-LRA and R-ALRA algorithms always achieve a lower error than TM.
Notably, the estimation error of the R-ALRA algorithm is even lower than R-LRA
when $T_{p}$ is large. This is because, with larger noise, the effectiveness
of the trace-minimization relaxation is undermined for TM because of the
greater uncertainty of the measured power values. In comparison, the R-LRA
algorithm may be more likely to converge to a locally optimal solution since
the estimation problem is non-convex, especially when $T_{p}$ is large. For
R-ALRA, however, the optimizations of $\bm{H}$ and $\bm{H}_{r}$ for $b\geq 2$
and $b=1$, are both convex and always have only one globally optimal solution,
leading to better estimation robustness with higher noise power and larger
$T_{p}$. In contrast, when $T_{p}$ is small, e.g., $T_{p}\leq 300$, the R-LRA
algorithm achieves much lower estimation error than R-ALRA.

Despite that the estimation error varies greatly for different values of
$\sigma^{2}$ and $N_{0}$, the effective channel gain achieved by IRS
reflection design based on estimated channel autocorrelation matrices does not
change significantly, as shown in Fig. 6(b). For $N_{0}=1$, the effective
channel gains achieved by all schemes for $\sigma^{2}=-90$ dBm and $-85$ dBm
are close to each other. Besides, the performance difference between $N_{0}=8$
and $N_{0}=80$ is small. For the proposed schemes based on the R-LRA and
R-ALRA algorithms, the optimized IRS reflection vector $\bm{v}$ based on the
estimated channel autocorrelation matrix $\hat{\bm{H}}$ or $\hat{\bm{H}}_{r}$
is still near-optimal even when the estimation error is large. By comparing
Fig. 6(a) to Fig. 6(b), it can be observed that as long as the normalized
estimation error is no larger than $10^{-1}$, the optimized IRS reflection
vector based on the estimated channel autocorrelation matrix is nearly
optimal, and thus the effective channel gain almost achieves the upper bound.
As such, significant improvements can still be obtained by the R-LRA and
R-ALRA schemes compared to other benchmark schemes.

###  VII-D Performance with Quantized Power Measurement

Finally, the impact of quantization errors on the performance of R-LRA and
R-ALRA algorithms is evaluated. We fix $T_{p}=200$, $\sigma^{2}=-90$ dBm and
$N_{0}=8$, while the estimation error and effective channel gain are evaluated
with different values of the quantization gap $D$.

(a) Estimation error.

(b) Effective channel gain.

Figure 7: Performance with quantized power measurement for $b=1$
($\sigma^{2}=-90$ dBm, $N_{0}=8$, and $T_{p}=200$).

In Fig. 7(a), the normalized estimation error of different schems versus $D$
is shown for $b=1$. Correspondingly, the effective channel gain under the same
condition is given in Fig. 7(b). It can be observed that the estimation error
increases and the effective channel gain decreases with the quantization gap
$D$. The estimation error of R-LRA is always lower than that of TM or R-ALRA,
while the error of R-ALRA becomes higher than TM for $T_{p}\geq 8$ dB. This is
caused by introducing the variable $\hat{\bm{q}}$ and approximating problem
(50) with problem (51) to obtain closed-form solutions. With larger $D$, the
accuracy of the approximation decreases and thus the estimation error
increases. Nevertheless, the effective channel gains of both R-LRA and R-ALRA
are close to each other and they are higher compared to all benchmark schemes
for all values of $D$, which verifies the robustness of the proposed schemes
against quantization error. It is worth mentioning that in the current
protocol [33], $D=1$ dB is employed, for which the performance loss caused by
quantization for both R-LRA and R-ALRA is quite small based on Fig. 7(b).

##  VIII Conclusion

This paper studied the IRS-cascaded channel estimation problem based on
received signal power measurement at the user. Since the signal phase
information is not available in power measurement, the channel autocorrelation
matrix was estimated by solving a rank-minimization problem. Based on the
analysis of existence and uniqueness of the solution, two algorithms of
different computational complexity were proposed to solve the channel
autocorrelation estimation problem. By relaxing the non-convex low-rank
constrained problem to a fractional programming problem, the proposed LRA
algorithm applied the alternating optimization method to iteratively approach
a low-rank matrix solution. To reduce computational complexity, the ALRA
algroithm was also developed to obtain an approximate solution in closed-form
during each iteration. Moreover, robust extensions of the LRA and ALRA
algorithms were proposed based on practical power measurement, where the
effects of receiver noise and quantization error in power measurement were
considered. The convergence and efficiency of the proposed algorithms were
validated via simulations, demonstrating that small estimation errors of the
channel autocorrelation matrix can be achieved. It was also shown that by
applying the IRS reflection designs based on the channel autocorrelation
matrices estimated by the proposed LRA and ALRA algorithms, significant
improvement in the effective channel gain between BS and user can be achieved
compared to other benchmark schemes. Furthermore, the performance of the
robust LRA and robust ALRA algorithms was evaluated using practical power
measurement with noise and quantization errors, which demonstrated their
robustness in improving the effective channel gain for IRS-assisted
communication systems, even with imperfect power measurement.

##  Appendix A Proof of Lemma 1

Note that the set of $(N\times N)$-dimensional hermitian matrices forms an
$N^{2}$-dimensional linear space over $\mathbb{R}$ [29]. Define
$\mathcal{S}_{V}^{b}(N)=\\{\bm{V}\in\mathbb{C}^{N\times
N}|\bm{V}=\bm{v}\bm{v}^{H},\bm{v}\in\Phi_{b}^{N}\\}$ as the set of all
possible autocorrelation matrices of the IRS reflection vector and
$\mathcal{S}_{X}^{b}(N)$ as the linear space spanned by elements from
$\mathcal{S}_{V}^{b}(N)$ with real coefficients. Due to the unit-amplitude of
the entries of $\bm{v}$, the diagonal entries in matrix $\bm{V}$ are always
equal to one, which results in a dimensional defficiency for
$\mathcal{S}_{X}^{b}(N)$ in the $N^{2}$-dimensional hermitian matrix space.
Specifically, it can be proved that
$\text{dim}(\mathcal{S}_{X}^{b}(N))=\mathcal{D}_{N}^{(b)}=N^{2}-N+1$ for
$b\geq 2$ and $(N^{2}-N)/2+1$ for $b=1$ (see Appendix D). Since we have
$D_{V}\leq\text{dim}(\mathcal{S}_{X}^{b}(N))=\mathcal{D}_{N}^{(b)}$, Lemma 1
is proved.

##  Appendix B Proof of Proposition 1

The existence of $\bar{\bm{H}}$ is straightforward because problem (4) follows
(3), where the actual channel autocorrelation matrix is always a feasible
solution.

Next, the uniqueness is analyzed. For $b\geq 2$,
$\\{\bm{v}_{t}\in\Phi_{b}^{N}$, $t=1,\ldots,T_{p}\\}$ are complex vectors.
When $T_{p}\geq D_{V}=N^{2}-N+1$, matrices $\\{\bm{V}_{t},t=1,\ldots,T_{p}\\}$
span the space $\mathcal{S}_{X}^{b}(N)$ and thus the projection of $\bm{H}$ in
$\mathcal{S}_{X}^{b}(N)$ is determined. By analyzing the orthonormal basis of
$\mathcal{S}_{X}^{b}(N)$, it can be shown that the values of
$\text{tr}(\bm{H})$ and all non-diagonal entries in $\bm{H}$ can be uniquely
determined by (4a), i.e., $\text{tr}(\bm{H})=\text{tr}(\bar{\bm{H}})$ and
$H_{nm}=\bar{H}_{nm}$, $\forall n\neq m$. Under such conditions, it can be
proved that $\bm{H}=\bar{\bm{H}}$ is the only rank-one solution for problem
(4) when $N\geq 3$ (see Appendix E).

For $b=1$, $\bm{v}_{t}\in\Phi_{1}^{N}=\\{\pm 1\\}^{N}$ is always a real vector
for any $t$, indicating that $\bm{V}_{t}=\bm{V}_{t}^{T}$, $\forall t$. Then,
we have

$\text{tr}(\bm{H}\bm{V}_{t})=\text{tr}((\bm{H}\bm{V}_{t})^{T})=\text{tr}(\bm{V}_{t}^{T}\bm{H}^{*})=\text{tr}(\bm{H}^{*}\bm{V}_{t}),$
(62)

which means that the imaginary part of $\bm{H}$ does not influence the
received signal power at the user for $b=1$. Thus, for any solution
$\hat{\bm{H}}$ for problem (4), its conjugate $\hat{\bm{H}}^{*}$ is also a
solution. Similar to the case of $b\geq 2$, by analyzing the orthonormal basis
of $\mathcal{S}_{X}^{b}(N)$ for $b=1$, it can be shown that the values of
$\text{tr}(\bm{H})$ and real parts of all non-diagonal entries in $\bm{H}$ can
be uniquely determined by (4a) when $T_{p}\geq D_{V}=\mathcal{D}_{N}^{(1)}$,
i.e., $\text{tr}(\bm{H})=\text{tr}(\bar{\bm{H}})$ and
$\text{real}(H_{nm})=\text{real}(\bar{H}_{nm})$, $\forall n\neq m$. Based on
that, it can be proved that $\bar{\bm{H}}$ and $\bar{\bm{H}}^{*}$ are the only
two rank-one solutions for problem (4) when $N\geq 6$ (see Appendix E).

##  Appendix C Proof of Proposition 2

It is easy to verify that $\bar{\bm{H}}_{r}$ is a solution for problem (5),
which guarantees the existence of the solution. For the uniqueness, any
solution for problem (5), denoted by $\hat{\bm{H}}_{r}$, is symmetric and
semidefinite with $\text{rank}(\hat{\bm{H}}_{r})\leq 2$. With eigenvalue
decomposition, we have
$\hat{\bm{H}}_{r}=\bm{U}_{r}\text{diag}([\alpha_{1},\alpha_{2}])\bm{U}_{r}^{T}$,
where $\bm{U}_{r}=[\bm{q}_{r1},\bm{q}_{r2}]\in\mathbb{R}^{N\times 2}$ and
$\alpha_{1},\alpha_{2}\geq 0$. Denote
$\hat{\bm{h}}_{1}=\sqrt{\alpha_{1}}\bm{q}_{r1},\hat{\bm{h}}_{2}=\sqrt{\alpha_{2}}\bm{q}_{r2}\in\mathbb{R}^{N\times
1}$, then
$\hat{\bm{H}}_{r}=\hat{\bm{h}}_{1}\hat{\bm{h}}_{1}^{T}+\hat{\bm{h}}_{2}\hat{\bm{h}}_{2}^{T}$.
Consider $\hat{\bm{h}}=\hat{\bm{h}}_{1}+j\hat{\bm{h}}_{2}$ and
$\hat{\bm{H}}=\hat{\bm{h}}\hat{\bm{h}}^{H}$. Obviously, we have
$\hat{\bm{H}}_{r}=\text{Re}(\hat{\bm{H}})$, $\text{rank}(\hat{\bm{H}})=1$ and
$p_{0}\text{tr}(\hat{\bm{H}}\bm{V}_{t})=p_{0}\text{tr}(\hat{\bm{H}}_{r}\bm{V}_{t})=p_{t}$,
$t=1,\ldots,T_{p}$. Thus, $\hat{\bm{H}}$ is a solution for problem (4), which
means $\hat{\bm{H}}=\bar{\bm{H}}$ or $\bar{\bm{H}}^{*}$ according to
Proposition 1. Then, we have
$\hat{\bm{H}}_{r}=\text{Re}(\hat{\bm{H}})=\bar{\bm{H}}_{r}$, which is uniquely
determined.

##  Appendix D Complementary Proof of Lemma 1 in Appendix A

In this section, the detailed proof for the dimension of space
$\mathcal{S}_{X}^{b}(N)$ defined in Appendix A is given to complete the proof
of Lemma 1, where it is claimed that
$\text{dim}(\mathcal{S}_{X}^{b}(N))=\mathcal{D}_{N}^{(b)}$. The function
$\mathcal{D}_{N}^{(b)}$ is defined in Lemma 1 as

$\mathcal{D}_{N}^{(b)}=\left\\{\begin{aligned}
&\frac{N^{2}-N}{2}+1,~{}~{}~{}b=1\\\ &N^{2}-N+1,~{}~{}~{}~{}b\geq
2\end{aligned}\right..$ (63)

To prove $\text{dim}(\mathcal{S}_{X}^{b}(N))=\mathcal{D}_{N}^{(b)}$, we first
denote the $N^{2}$-dimensional hermitian matrices space as
$\mathcal{S}_{H}(N)$ and explicitly define an orthonormal basis for it. Then,
we show that $\mathcal{S}_{X}^{b}(N)$ can be spanned by
$\mathcal{D}_{N}^{(b)}$ matrices in the basis.

###  D-A Orthonormal Basis of Space $\mathcal{S}_{H}(N)$

Consider vector $\bm{b}_{1}=\bm{1}_{N}/\sqrt{N}\in\mathbb{R}^{N\times 1}$, and
denote $\bm{b}_{2},\ldots,\bm{b}_{N}$ as the vectors that form an orthonormal
basis along with $\bm{b}_{1}$ for the $N$-dimensional real vector space
$\mathbb{R}^{N\times 1}$. Then, an orthonormal basis for space
$\mathcal{S}_{H}(N)$ defined in [29] can be given by the following lemma.

######  Lemma 2.

Matrices
$\\{\bm{B}_{1},\ldots,\bm{B}_{N};\bm{E}_{nl}^{(R)},\bm{E}_{nl}^{(I)},1\leq
n<l\leq N\\}$ form an orthonormal basis for $\mathcal{S}_{H}(N)$ defined as

$\displaystyle\bm{B}_{n}$ $\displaystyle=\textup{diag}(\bm{b}_{n}),~{}1\leq
n\leq N,$ (64a) $\displaystyle\bm{E}_{nl}^{(R)}$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\bm{e}_{n}\bm{e}_{l}^{T}+\bm{e}_{l}\bm{e}_{n}^{T}\right),~{}1\leq
n<l\leq N,$ (64b) $\displaystyle\bm{E}_{nl}^{(I)}$
$\displaystyle=\frac{j}{\sqrt{2}}\left(\bm{e}_{n}\bm{e}_{l}^{T}-\bm{e}_{l}\bm{e}_{n}^{T}\right),~{}1\leq
n<l\leq N,$ (64c)

where $\bm{e}_{n}\in\mathbb{R}^{N\times 1}$ has a $1$ for the $n$-th element
and all other elements equal to zero.

###### Proof:.

It is easy to verify that these matrices are hermitian and orthonormal to each
other under the Frobenius product $\text{tr}(\bm{A}^{H}\bm{B})$ for arbitrary
matrices $\bm{A}$ and $\bm{B}$. Meanwhile, on one hand, for any $N\times N$
hermitian matrix $\bm{A}$, its main diagonal entries form an $N$-dimensional
real vector $\bm{d}_{\bm{A}}=\text{diag}(\bm{A})$. Thus, there exists a real
vector
$\bm{c}_{\bm{A}}=[c_{\bm{A},1},\ldots,c_{\bm{A},N}]^{T}\in\mathbb{R}^{N\times
1}$ such that $\bm{d}_{\bm{A}}=\sum_{n=1}^{N}{c_{\bm{A},n}\bm{b}_{n}}$, i.e.,

$\text{diag}(\bm{d}_{\bm{A}})=\sum_{n=1}^{N}{c_{\bm{A},n}\bm{B}_{n}}.$ (65)

On the other hand, entries of matrices $\bm{E}_{nl}^{(R)}$ and
$\bm{E}_{nl}^{(I)}$ are zero except for two on the $n$-th row and $l$-th
column, and the $l$-th row and $n$-th column, $\forall n<l$. Thus, it can be
verified that $\bm{A}$ can be written as

$\bm{A}=\sqrt{2}\sum_{n<l}{\text{Re}(A_{nl})\bm{E}_{nl}^{(R)}+\text{Im}(A_{nl})\bm{E}_{nl}^{(I)}}+\sum_{n=1}^{N}{c_{\bm{A},n}\bm{B}_{n}},$
(66)

which is the linear combination of $N(N-1)+N=N^{2}$ matrices with $N^{2}$ real
coefficients. Therefore, the $N^{2}$ matrices
$\\{\bm{B}_{1},\ldots,\bm{B}_{N};\bm{E}_{nl}^{(R)},\bm{E}_{nl}^{(I)}$, $1\leq
n<l\leq N\\}$ form an orthonormal basis for the space of hermitian matrices. ∎

###  D-B Dimension of Space $\mathcal{S}_{X}^{b}(N)$

To show $\text{dim}(\mathcal{S}_{X}^{b}(N))=\mathcal{D}_{N}^{(b)}$, we provide
a two-fold proof, where we first prove
$\text{dim}(\mathcal{S}_{X}^{b}(N))\leq\mathcal{D}_{N}^{(b)}$ and then
$\text{dim}(\mathcal{S}_{X}^{b}(N))\geq\mathcal{D}_{N}^{(b)}$.

For the first part, we show
$\text{dim}(\mathcal{S}_{X}^{b}(N))\leq\mathcal{D}_{N}^{(b)}$. Since
$\mathcal{S}_{X}^{b}(N)$ is spanned by $\mathcal{S}_{V}^{b}(N)$ with real
coefficients, we have $\mathcal{S}_{X}^{b}(N)\in\mathcal{S}_{H}(N)$ and thus
$\text{dim}(\mathcal{S}_{X}^{b}(N))\leq N^{2}$. Due to the unit-amplitude of
the IRS reflecting coefficients, however, the diagonal entries in any matrix
$\bm{V}\in\mathcal{S}_{V}^{b}(N)$ are always equal to 1, i.e.,
$\text{diag}(\bm{V})=\bm{1}_{N}=\sqrt{N}\bm{b}_{1}$. Then, it is easy to
verify that $\bm{V}$ is orthogonal to $\bm{B}_{2},\ldots,\bm{B}_{N}$ for any
values of $b$, i.e.,

$\text{tr}(\bm{V}^{H}\bm{B}_{n})=\sqrt{N}\bm{b}_{1}^{T}\bm{b}_{n}=0,~{}2\leq
n\leq N.$ (67)

Moreover, for $b=1$, $\bm{V}$ is always a real matrix and thus we have
$\text{tr}(\bm{V}^{H}\bm{E}_{nl}^{(I)})=\text{Im}(V_{nl})=0$, $\forall n\neq
l$. Thus, for $b\geq 2$, $\mathcal{S}_{X}^{b}(N)$ should be a subspace of the
space spanned by $\\{\bm{B}_{1};\bm{E}_{nl}^{(R)},\bm{E}_{nl}^{(I)},1\leq
n<l\leq N\\}$, leading to $\text{dim}(\mathcal{S}_{X}^{b}(N))\leq
N^{2}-N+1=\mathcal{D}_{N}^{(b)}$, while for $b=1$, $\mathcal{S}_{X}^{1}(N)$
should be a subspace of the space spanned by
$\\{\bm{B}_{1};\bm{E}_{nl}^{(R)},1\leq n<l\leq N\\}$, leading to
$\text{dim}(\mathcal{S}_{X}^{1}(N))\leq N(N-1)/2+1=\mathcal{D}_{N}^{(1)}$.
Therefore, $\text{dim}(\mathcal{S}_{X}^{b}(N))\leq\mathcal{D}_{N}^{(b)}$ holds
for any $b$.

For the second part of the proof, we prove by induction. To prove
$\text{dim}(\mathcal{S}_{X}^{b}(N))\geq\mathcal{D}_{N}^{(b)}$, it suffices to
show that there exist at least $\mathcal{D}_{N}^{(b)}$ matrices in
$\mathcal{S}_{V}^{b}(N)$ such that they are linearly independent with real
coefficients. For the following, we first prove the case of $b=1$, and then
consider the case of $b\geq 2$.

For $b=1$, we show that there always exists
$\mathcal{D}_{N}^{(1)}=(N^{2}-N)/2+1$ matrices in $\mathcal{S}_{V}^{1}(N)$
that are linearly independent with real coefficients. It is easy to verify
that this is true when $N\leq 3$. Assume that this is true for $N=L\geq 3$,
and let $K_{1}=\mathcal{D}_{L}^{(1)}$ for simplicity, then we have to show
that this is true for $N=L+1$. Suppose $\bm{V}_{k}=\bm{v}_{k}\bm{v}_{k}^{H}$,
$k=1,\ldots,{K_{1}}$ are linearly independent matrices in
$\mathcal{S}_{V}^{1}(L)$, where $\bm{v}_{k}\in{\\{\pm 1\\}}^{L}$,
$k=1,\ldots,K_{1}$. Define matrix
$\bm{W}=[\bm{v}_{1},\ldots,\bm{v}_{K_{1}}]\in\mathbb{R}^{L\times{K_{1}}}$,
which has full row rank $L$ according to the following lemma.

######  Lemma 3.

The matrix $\bm{W}$ has full row rank $L$ for $L\geq 3$.

###### Proof:.

As $K_{1}>L$ for $L\geq 3$, we have $\text{rank}(\bm{W})\leq L<{K_{1}}$. If
$\text{rank}(\bm{W})<L$, then the space spanned by the columns of $\bm{W}$,
denoted as $\text{span}(\bm{W})$, has a dimension lower than $L$. Since all
vectors in $\\{\pm 1\\}^{L}$ span the $L$-dimensional space $\mathbb{R}^{L}$,
there exists a vector $\bm{x}\in\\{\pm 1\\}^{L}$ such that
$\bm{x}\notin\text{span}(\bm{W})$. Let $\bm{x}=\bm{v}_{x}+\bm{y}$, where
$\bm{v}_{x}\in\text{span}(\bm{W})$ while $\bm{y}\perp\text{span}(\bm{W})$.
Note that $\bm{X}=\bm{x}\bm{x}^{H}\in\mathcal{S}_{V}^{1}(N)$ and thus matrices
$\\{\bm{V}_{1},\ldots,\bm{V}_{K_{1}},\bm{X}\\}$ must be linearly dependent
with real coefficients according to the first part of the proof, because there
are ${K_{1}}+1>\mathcal{D}_{L}^{(1)}$ hermitian matrices from
$\mathcal{S}_{V}^{1}(N)$ in total. Hence, there exists
$\bm{\mu}\in\mathbb{R}^{{K_{1}}\times 1}$ and $\mu_{x}\in\mathbb{R}$, such
that they are not all zeros and

$\mu_{x}\bm{X}+\sum_{k=1}^{{K_{1}}}{\mu_{k}\bm{V}_{k}}=\bm{0}_{L\times L}.$
(68)

Consider the matrix $\bm{Y}=\bm{y}\bm{y}^{H}$ and we have

$\displaystyle 0$
$\displaystyle=\text{tr}\left(\bm{Y}^{H}\left(\mu_{x}\bm{X}+\sum_{k=1}^{{K_{1}}}{\mu_{k}\bm{V}_{k}}\right)\right)$
(69a)
$\displaystyle=\mu_{x}\text{tr}(\bm{Y}^{H}\bm{X})+\sum_{k=1}^{{K_{1}}}{\mu_{k}\text{tr}(\bm{Y}^{H}\bm{V}_{k})}.$
(69b)

Note that
$\text{tr}(\bm{Y}^{H}\bm{X})=|\bm{y}^{H}(\bm{v}_{x}+\bm{y})|^{2}=\|\bm{y}\|_{2}^{4}$
and $\text{tr}(\bm{Y}^{H}\bm{V}_{k})=|\bm{y}^{H}\bm{v}_{k}|^{2}=0$, $\forall
k$. Thus, equation (69) is equivalent to $\mu_{x}\|\bm{y}\|_{2}^{4}=0$, which
leads to $\mu_{x}=0$ as $\bm{y}$ is nonzero. Then, equation (68) becomes
$\sum_{k=1}^{{K_{1}}}{\mu_{k}\bm{V}_{k}}=\bm{0}$ for some
$\bm{\mu}\neq\bm{0}$. This is contradictory to the assumption that
$\\{\bm{V}_{k},k=1,\ldots,{K_{1}}\\}$ are linearly independent with real
coefficients. Therefore, $\bm{W}$ should have full row rank $L$. ∎

Given Lemma 3, we implicitly construct $\mathcal{D}_{L+1}^{(1)}$ linearly
independent matrices in $S_{V}^{1}(L+1)$ for the case of $N=L+1$. Define
$(L+1)$-dimensional vectors

$\bm{u}_{k}=\left[\begin{array}[]{c}\bm{v}_{k}\\\
1\end{array}\right],~{}\bm{y}_{k}=\left[\begin{array}[]{c}\bm{v}_{k}\\\
-1\end{array}\right],~{}k=1,\ldots,{K_{1}},$ (70)

and matrices $\bm{U}_{k}=\bm{u}_{k}\bm{u}_{k}^{H}$ and
$\bm{Y}_{k}=\bm{y}_{k}\bm{y}_{k}^{H}$, $\forall k$. Obviously,
$\bm{U}_{k},\bm{Y}_{k}\in S_{V}^{1}(L+1)$, $\forall k$, and there are
$2{K_{1}}=2\mathcal{D}_{L}^{(1)}$ matrices in total. Note that for $L\geq 3$,
we have $2\mathcal{D}_{L}^{(1)}>\mathcal{D}_{L+1}^{(1)}$, so matrices
$\\{\bm{U}_{k},\bm{Y}_{k},k=1,\ldots,{K_{1}}\\}$ are linearly dependent with
real coefficients according to the first part of the proof, i.e., there exists
real vectors
$\bm{\alpha}=[\alpha_{1},\ldots,\alpha_{K_{1}}]^{T},\bm{\beta}=[\beta_{1},\ldots,\beta_{K_{1}}]^{T}\in\mathbb{R}^{{K_{1}}\times
1}$ such that

$\displaystyle\bm{0}_{(L+1)\times(L+1)}$
$\displaystyle=\sum_{k=1}^{{K_{1}}}{\alpha_{k}\bm{U}_{k}+\beta_{k}\bm{Y}_{k}}$
(71a)
$\displaystyle=\sum_{k=1}^{{K_{1}}}{\left[\begin{array}[]{cc}(\alpha_{k}+\beta_{k})\bm{V}_{k}&(\alpha_{k}-\beta_{k})\bm{v}_{k}\\\
(\alpha_{k}-\beta_{k})\bm{v}_{k}^{H}&\alpha_{k}+\beta_{k}\end{array}\right]}.$
(71d)

Thus, we have $\sum_{k=1}^{{K_{1}}}{(\alpha_{k}+\beta_{k})\bm{V}_{k}}=\bm{0}$.
However, $\\{\bm{V}_{k},k=1,\ldots,{K_{1}}\\}$ are linearly independent with
real coefficients by assumption, which leads to
$\bm{\alpha}+\bm{\beta}=\bm{0}$. Then,
$\sum_{k=1}^{{K_{1}}}(\alpha_{k}-\beta_{k})\bm{v}_{k}=2\sum_{k=1}^{{K_{1}}}{\alpha_{k}\bm{v}_{k}}=2\bm{W}\bm{\alpha}=\bm{0}$.
Reversely, it is easy to verify that (71) holds if
$\bm{\alpha}+\bm{\beta}=\bm{0}$ and $\bm{W}\bm{\alpha}=\bm{0}$. Thus, (71)
holds if and only if

$\bm{F}\left[\begin{array}[]{c}\bm{\alpha}\\\
\bm{\beta}\end{array}\right]=\bm{0}_{(K_{1}+L)\times
1},~{}\bm{F}=\left[\begin{array}[]{cc}\bm{I}_{K_{1}}&\bm{I}_{K_{1}}\\\
\bm{W}&\bm{0}_{L\times K_{1}}\end{array}\right].$ (72)

Therefore, the maximum number of linearly independent matrices among
$\\{\bm{U}_{k},\bm{Y}_{k},k=1,\ldots,{K_{1}}\\}$ with real coefficients,
denoted as $M_{1}$, should be equal to the number of linearly independent
columns of matrix $\bm{F}\in\mathbb{R}^{(K_{1}+L)\times(2K_{1})}$, i.e.,
$M_{1}=\text{rank}(\bm{F})$, which can be easily verified to be $K_{1}+L$.
Since ${K_{1}}=\mathcal{D}_{L}^{(1)}=(L^{2}-L)/2+1$, we have
$M_{1}=(L^{2}+L)/2+1=\mathcal{D}_{L+1}^{(1)}$. Thus, there exists
$\mathcal{D}_{L+1}^{(1)}$ linearly independent matrices with real coefficients
in the set $S_{V}^{1}(L+1)$, which completes the proof of
$\text{dim}(\mathcal{S}_{X}^{1}(N))\geq\mathcal{D}_{N}^{(1)}$.

For $b\geq 2$, we have $\Phi_{2}=\\{\pm 1,\pm j\\}\subseteq\Phi_{b}$ and
$\mathcal{D}_{N}^{(b)}=\mathcal{D}_{N}^{(2)}=N^{2}-N+1$. Therefore, it
suffices to check the case of $b=2$ only, i.e., showing that there always
exists $\mathcal{D}_{N}^{(2)}$ matrices in $S_{V}^{2}(N)$ that are linearly
independent with real coefficients. Then, for any $b\geq 2$, the same matrices
can be found in $\mathcal{S}_{V}^{b}(N)$ such that they are linearly
independent with real coefficients.

For $b=2$, it is easy to verify the cases where $N\leq 3$. Assume that for
$N=L\geq 3$, there exists $K_{2}=\mathcal{D}_{L}^{(2)}=L^{2}-L+1$ matrices in
$S_{V}^{2}(L)$ that are linearly independent with real coefficients, which are
denoted as $\bm{V}_{k}=\bm{v}_{k}\bm{v}_{k}^{H}$, $k=1,\ldots,K_{2}$, with
$\bm{v}_{k}\in{\\{\pm 1,\pm j\\}}^{L}$, $\forall k$. Similar to the case of
$b=1$, we implicitly construct $\mathcal{D}_{L+1}^{(2)}$ linearly independent
matrices in $S_{V}^{2}(L+1)$ to prove the case of $N=L+1$. Specifically,
define $\bm{d}_{k}=(1-j)\bm{v}_{k}\in\\{1+j,1-j,-1+j,-1-j\\}^{L}$,
$k=1,\ldots,K$, then we have $\bm{V}_{k}=\bm{d}_{k}\bm{d}_{k}^{H}/2,\forall
k$. Note that the real and imaginary parts of vector $\bm{d}_{k}$ can be
considered as two binary variables independent of each other. Thus, it can be
decomposed into two parts, i.e., $\bm{d}_{k}=\bm{d}_{k,r}+j\bm{d}_{k,i}$,
where $\bm{d}_{k,r},\bm{d}_{k,i}\in\\{\pm 1\\}^{L}=\Phi_{1}^{L}$. Then, define
vector $\bm{q}_{k}$ by concatenating $\bm{d}_{k,r}$ and $\bm{d}_{k,i}$
together, i.e.,

$\bm{q}_{k}=\left[\begin{array}[]{c}\bm{d}_{k,r}\\\
\bm{d}_{k,i}\end{array}\right]\in\\{\pm 1\\}^{2L},~{}k=1,\ldots,K_{2},$ (73)

and matrix
$\bm{Q}=[\bm{q}_{1},\ldots,\bm{q}_{K_{2}}]\in\mathbb{R}^{2L\times{K_{2}}}$.
Let $R=2L-\text{rank}(\bm{Q})\geq 0$. If $R>0$, i.e.,
$\text{rank}(\bm{Q})<2L$, there exists vectors
$\bm{q}_{{K_{2}}+1},\ldots,\bm{q}_{{K_{2}}+R}\in\\{\pm 1\\}^{2L}$ such that
$\tilde{\bm{Q}}=[\bm{Q},\bm{q}_{{K_{2}}+1},\ldots,\bm{q}_{{K_{2}}+R}]\in\mathbb{R}^{2L\times({K_{2}}+R)}$
has full row rank $2L$. Correspondingly, for $k={K_{2}}+1,\ldots,{K_{2}}+R$,
we decompose $\bm{q}_{k}$ as
$\bm{q}_{k}=[\bm{d}_{k,r}^{T},\bm{d}_{k,i}^{T}]^{T}$ with
$\bm{d}_{k,r},\bm{d}_{k,i}\in\mathbb{R}^{L\times 1}$ and define
$\bm{d}_{k}=\bm{d}_{k,r}+j\bm{d}_{k,i}$ and
$\bm{v}_{k}=\bm{d}_{k}/(1-j)\in\Phi_{2}^{L}$. If $R=0$, however, let
$\tilde{\bm{Q}}=\bm{Q}$.

Next, for $N=L+1$, consider $(L+1)$-dimensional vectors

$\bm{u}_{k}=\left[\begin{array}[]{c}\bm{v}_{k}\\\
1\end{array}\right],~{}\bm{y}_{k}=\left[\begin{array}[]{c}\bm{v}_{k}\\\
-1\end{array}\right],~{}k=1,\ldots,{K_{2}}+R,$ (74)

and matrices $\bm{U}_{k}=\bm{u}_{k}\bm{u}_{k}^{H}$ and
$\bm{Y}_{k}=\bm{y}_{k}\bm{y}_{k}^{H}$, $\forall k$. Matrices $\bm{U}_{k}$ and
$\bm{Y}_{k}$ are included in set $S_{V}^{2}(L+1)$ and there are $2({K_{2}}+R)$
matrices in total. Since $2({K_{2}}+R)\geq 2{K_{2}}>\mathcal{D}_{L+1}^{(2)}$
for $L\geq 3$, they must be linearly dependent with real coefficients
according to the first part of the proof. Thus, there exists real vectors
$\bm{\alpha}=[\alpha_{1},\ldots,\alpha_{{K_{2}}+R}]^{T},\bm{\beta}=[\beta_{1},\ldots,\beta_{{K_{2}}+R}]^{T}\in\mathbb{R}^{({K_{2}}+R)\times
1}$ such that

$\displaystyle\bm{0}_{{(L+1)}\times{(L+1)}}$
$\displaystyle=\sum_{k=1}^{{K_{2}}+R}{\alpha_{k}\bm{U}_{k}+\beta_{k}\bm{Y}_{k}}$
(75a)
$\displaystyle=\sum_{k=1}^{{K_{2}}+R}{\left[\begin{array}[]{cc}(\alpha_{k}+\beta_{k})\bm{V}_{k}&(\alpha_{k}-\beta_{k})\bm{v}_{k}\\\
(\alpha_{k}-\beta_{k})\bm{v}_{k}^{H}&\alpha_{k}+\beta_{k}\end{array}\right]}.$
(75d)

Define $\mathcal{G}_{s}$ as the set of all vectors
$\bm{s}\in\mathbb{R}^{(K_{2}+R)\times 1}$ such that
$\sum_{k=1}^{{K_{2}}+R}{s_{k}\bm{V}_{k}}=\bm{0}$ holds, then
$\bm{\alpha}+\bm{\beta}\in\mathcal{G}_{s}$. Obviously, $\mathcal{G}_{s}$ is a
linear space whose dimension is bounded by
$\text{dim}(\mathcal{G}_{s})\leq{K_{2}}+R-{K_{2}}=R$ because
$\bm{V}_{1},\ldots,\bm{V}_{K_{2}}$ are linearly independent with real
coefficients. Meanwhile, by defining $\bm{z}=\bm{\alpha}-\bm{\beta}$, we have
$\sum_{k=1}^{{K_{2}}+R}{z_{k}\bm{v}_{k}}=\bm{0}_{L\times 1}$ and thus
$\sum_{k=1}^{{K_{2}}+R}{z_{k}\bm{d}_{k}}=(1-j)\sum_{{K_{2}}=1}^{K_{2}+R}{z_{k}\bm{v}_{k}}=\bm{0}_{L\times
1}$. As $\bm{z}$ is a real vector, we have

$\tilde{\bm{Q}}\bm{z}=\sum_{k=1}^{{K_{2}}+R}{z_{k}\left[\begin{array}[]{c}\bm{d}_{k,r}\\\
\bm{d}_{k,i}\end{array}\right]}=\sum_{k=1}^{{K_{2}}+R}\left[\begin{array}[]{c}{z_{k}\text{Re}(\bm{d}_{k})}\\\
{z_{k}\text{Im}(\bm{d}_{k})}\end{array}\right]=\bm{0}_{2L\times 1},$ (76)

which means $\bm{z}$ is in the null space of $\tilde{\bm{Q}}$, i.e.,
$\mathcal{N}(\tilde{\bm{Q}})$. The dimension of $\mathcal{N}(\tilde{\bm{Q}})$
is ${K_{2}}+R-2L$, because $\tilde{\bm{Q}}\in\mathbb{R}^{2L\times({K_{2}}+R)}$
has full row rank $2L$ by definition. Reversely, one can easily verify that
equation (75) is true as long as $\bm{\alpha}+\bm{\beta}\in{\mathcal{G}_{s}}$
and $\bm{\alpha}-\bm{\beta}\in\mathcal{N}(\tilde{\bm{Q}})$. Therefore,
$\sum_{k=1}^{{K_{2}}+R}{\alpha_{k}\bm{U}_{k}+\beta_{k}\bm{Y}_{k}}=\bm{0}_{{K_{2}}\times{K_{2}}}$
if and only if $\bm{\alpha}+\bm{\beta}\in\mathcal{G}_{s}$ and
$\bm{\alpha}-\bm{\beta}\in\mathcal{N}(\tilde{\bm{Q}})$, i.e.,

$\displaystyle\left[\begin{array}[]{c}\bm{\alpha}+\bm{\beta}\\\
\bm{\alpha}-\bm{\beta}\end{array}\right]\in\mathcal{G}$
$\displaystyle=\bigg{\\{}\bm{x}\in\mathbb{R}^{2({K_{2}}+R)\times
1}\bigg{|}\bm{x}=\left[\begin{array}[]{c}\bm{s}\\\ \bm{z}\end{array}\right],$
(77)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\bm{s}\in\mathcal{G}_{s},\
\bm{z}\in\mathcal{N}(\tilde{\bm{Q}})\bigg{\\}}.$

where $\mathcal{G}$ is a linear space with dimension
$\text{dim}(\mathcal{G})=\text{dim}(\mathcal{G}_{s})+\text{dim}(\mathcal{N}(\tilde{\bm{Q}}))\leq
K_{2}+2R-2L$. Note that

$\left[\begin{array}[]{c}\bm{\alpha}\\\
\bm{\beta}\end{array}\right]=\bm{F}_{s}\left[\begin{array}[]{c}\bm{\alpha}+\bm{\beta}\\\
\bm{\alpha}-\bm{\beta}\end{array}\right],~{}\bm{F}_{s}=\frac{1}{2}\left[\begin{array}[]{cc}\bm{I}&\bm{I}\\\
\bm{I}&-\bm{I}\end{array}\right],$ (78)

where the transformation matrix $\bm{F}_{s}$ is invertible, so the dimension
of the solution space of $[\bm{\alpha}^{T},\bm{\beta}^{T}]^{T}$ for (75)
equals to the dimension of $\mathcal{G}$. Thus, the maximum number $M_{2}$ of
linearly independent matrices among
$\\{\bm{U}_{k},\bm{Y}_{k},k=1,\ldots,{K_{2}}+R\\}$ with real coefficients can
be obtained as

$\displaystyle M_{2}$ $\displaystyle=2({K_{2}}+R)-\text{dim}(\mathcal{G})$
(79a) $\displaystyle\geq 2({K_{2}}+R)-({K_{2}}+2R-2L)$ (79b)
$\displaystyle={K_{2}}+2L.$ (79c)

As ${K_{2}}=\mathcal{D}_{L}^{(2)}=L^{2}-L+1$, we have
$M_{2}\geq{K_{2}}+2L=L^{2}+L+1=\mathcal{D}_{L+1}^{(2)}$. Therefore, there
exists at least $\mathcal{D}_{L+1}^{(2)}$ linearly independent matrices in the
set $S_{V}^{2}(L+1)$, which completes the proof of
$\text{dim}(\mathcal{S}_{X}^{b}(N))\geq\mathcal{D}_{N}^{(b)}$ for the case of
$b=2$, and thus the case of $b\geq 2$.

Combining two parts of the proof, i.e.,
$\text{dim}(\mathcal{S}_{X}^{b}(N))\leq\mathcal{D}_{N}^{(b)}$ and
$\text{dim}(\mathcal{S}_{X}^{b}(N))\geq\mathcal{D}_{N}^{(b)}$, we conclude
that $\text{dim}(\mathcal{S}_{X}^{b}(N))=\mathcal{D}_{N}^{(b)}$, as stated in
Appendix A. Additionally, for $b\geq 2$, it has been shown that
$\mathcal{S}_{X}^{b}(N)$ is a subspace of the space spanned by
$\\{\bm{B}_{1};\bm{E}_{nl}^{(R)},\bm{E}_{nl}^{(I)},1\leq n<l\leq N\\}$ with
real coefficients, whose dimension is $N^{2}-N+1$, which is the same as the
dimension of $\mathcal{S}_{X}^{b}(N)$ given above. Hence,
$\mathcal{S}_{X}^{b}(N)$ equals to this spanned space and
$\\{\bm{B}_{1};\bm{E}_{nl}^{(R)},\bm{E}_{nl}^{(I)},1\leq n<l\leq N\\}$ form an
orthonormal basis of $\mathcal{S}_{X}^{b}(N)$ for $b\geq 2$. Similarly, an
orthonormal basis of $\mathcal{S}_{X}^{b}(N)$ for $b=1$ is given by
$\\{\bm{B}_{1};\bm{E}_{nl}^{(R)},1\leq n<l\leq N\\}$.

##  Appendix E Complementary Proof of Proposition 1 in Appendix B

In this section, the uniqueness of the solution for problem (4) demonstrated
in Proposition 1 is proved in details. First we consider the case of $b\geq 2$
and then $b=1$.

For $b\geq 2$, suppose there exists a matrix $\hat{\bm{H}}\neq\bar{\bm{H}}$
that is also a solution to Problem (4). Then, $\hat{\bm{H}}$ is a rank-one
hermitian matrix and we have
$\text{tr}(\hat{\bm{H}}\bm{V}_{t})=p_{t}=\text{tr}(\bar{\bm{H}}\bm{V}_{t})$,
$t=1,\ldots,T_{p}$. By defining $\hat{\bm{H}}=\hat{\bm{h}}\hat{\bm{h}}^{H}$
and $\mathcal{\bm{E}}=\bar{\bm{H}}-\hat{\bm{H}}\in\mathbb{C}^{N\times N}$, we
have $\text{rank}(\mathcal{\bm{E}})\leq 2$ and
$\text{tr}(\mathcal{\bm{E}}\bm{V}_{t})=0$ for $\forall t$. If
$D_{V}=\mathcal{D}_{N}^{(b)}=\text{dim}(\mathcal{S}_{X}^{b}(N))$, the matrix
$\mathcal{\bm{E}}$ is perpendicular to the space $\mathcal{S}_{X}^{b}(N)$,
which results in
$\text{tr}(\mathcal{\bm{E}}\bm{E}_{nl}^{(R)})=\text{tr}(\mathcal{\bm{E}}\bm{E}_{nl}^{(I)})=0$,
$1\leq n<l\leq N$, and also
$\text{tr}(\mathcal{\bm{E}}\bm{B}_{1})=\text{tr}(\mathcal{\bm{E}})/\sqrt{N}=0$.
According to the definitions of matrices
$\bm{E}_{nl}^{(R)},\bm{E}_{nl}^{(I)},1\leq n<l\leq N$, all non-diagonal
entries of $\mathcal{\bm{E}}$ are zero and thus we have

$\displaystyle\mathcal{\bm{E}}$
$\displaystyle=\left[\begin{array}[]{cccc}\mathcal{E}_{11}&0&\ldots&0\\\
0&\mathcal{E}_{22}&\ldots&0\\\ \vdots&\vdots&\ddots&0\\\
0&0&0&\mathcal{E}_{NN}\end{array}\right],$ (80)

where $\mathcal{E}_{nn}=|\bar{h}_{n}|^{2}-|\hat{h}_{n}|^{2}$, $\forall n$. Due
to $\text{rank}(\mathcal{\bm{E}})\leq 2$, we assume $\mathcal{E}_{nn}=0$ for
$n\geq 3$ without loss of generality, and thus $|\bar{h}_{n}|=|\hat{h}_{n}|$
for $n\geq 3$. By assuming $N\geq 3$, we have
$\mathcal{E}_{n1}=\bar{h}_{n}\bar{h}_{1}^{*}-\hat{h}_{n}\hat{h}_{1}^{*}=0$,
$n=3,\ldots,N$, leading to $|\bar{h}_{n}\bar{h}_{1}|=|\hat{h}_{n}\hat{h}_{1}|$
for $\forall n\geq 3$. Since $\bar{h}_{n}\neq 0$ holds with probability $1$
for any $n$, we find that $|\bar{h}_{1}|=|\hat{h}_{1}|$. Similarly,
$|\bar{h}_{2}|=|\hat{h}_{2}|$ can be obtained. Then,
$\mathcal{E}_{11}=\mathcal{E}_{22}=0$ and thus $\mathcal{\bm{E}}=\bm{0}$,
which indicates $\bar{\bm{H}}=\hat{\bm{H}}$. Hence, $\bar{\bm{H}}$ is the
unique solution to Problem (4) with probability $1$ for $b\geq 2$ when $N\geq
3$ and $D_{V}=\mathcal{D}_{N}^{(b)}$.

For $b=1$, the IRS reflection vectors are all real vectors and it is easy to
verify that $\bar{\bm{H}}$ and $\bar{\bm{H}}^{*}$ are both solutions for
problem (4). Suppose there exists a matrix $\hat{\bm{H}}$ that does not equal
to $\bar{\bm{H}}$ or $\bar{\bm{H}}^{*}$ but is also a solution for problem
(4), and its real part is denoted as
$\hat{\bm{H}}_{r}=\text{Re}(\hat{\bm{H}})$. Note that the ranks of both
$\bar{\bm{H}}_{r}=\text{Re}(\bar{\bm{H}})$ and $\hat{\bm{H}}_{r}$ are no more
than two and
$\text{tr}(\bar{\bm{H}}_{r}\bm{V}_{t})=\text{tr}(\bar{\bm{H}}\bm{V}_{t})=p_{t}=\text{tr}(\hat{\bm{H}}\bm{V}_{t})=\text{tr}(\hat{\bm{H}}_{r}\bm{V}_{t})$.
Similar to the case where $b\geq 2$, we define
$\mathcal{\bm{E}}_{r}=\bar{\bm{H}}_{r}-\hat{\bm{H}}_{r}\in\mathbb{R}^{N\times
N}$, which is a real symmetric matrix satisfying
$\text{rank}(\mathcal{\bm{E}}_{r})\leq 4$ and
$\text{tr}(\mathcal{\bm{E}}_{r}\bm{V}_{t})=0$, $t=1,\ldots,T_{p}$. When
$D_{V}=\mathcal{D}_{N}^{(1)}=\text{dim}(\mathcal{S}_{X}^{1}(N))$, the matrix
$\mathcal{\bm{E}}_{r}$ is perpendicular to the space $\mathcal{S}_{X}^{b}(N)$
and thus $\text{tr}(\mathcal{\bm{E}}_{r}\bm{E}_{nl}^{(R)})=0$, $1\leq n<l\leq
N$. As a result, $\mathcal{\bm{E}}_{r}$ can be written as

$\displaystyle\mathcal{\bm{E}}_{r}$
$\displaystyle=\left[\begin{array}[]{cccc}\mathcal{E}_{r,11}&0&\ldots&0\\\
0&\mathcal{E}_{r,22}&\ldots&0\\\ \vdots&\vdots&\ddots&0\\\
0&0&0&\mathcal{E}_{r,NN}\end{array}\right],$ (81)

where $\mathcal{E}_{r,nn}$ is the $n$-th diagonal element of
$\mathcal{\bm{E}}_{r}$. Without loss of generality, we have
$\mathcal{E}_{r,nn}=0$ for $n\geq 5$ due to
$\text{rank}(\mathcal{\bm{E}}_{r})\leq 4$. Since
$\bar{\bm{H}}=\bar{\bm{h}}\bar{\bm{h}}^{H}$, it is easy to verify that
$\bar{\bm{H}}_{r}=\bar{\bm{h}}_{r}\bar{\bm{h}}_{r}^{T}+\bar{\bm{h}}_{m}\bar{\bm{h}}_{m}^{T}$,
where $\bar{\bm{h}}_{r}=\text{Re}(\bar{\bm{h}})$ and
$\bar{\bm{h}}_{m}=\text{Im}(\bar{\bm{h}})$. Similarly, by defining
$\hat{\bm{H}}=\hat{\bm{h}}\hat{\bm{h}}^{H}$, we have
$\hat{\bm{H}}_{r}=\hat{\bm{h}}_{r}\hat{\bm{h}}_{r}^{T}+\hat{\bm{h}}_{m}\hat{\bm{h}}_{m}^{T}$,
where $\hat{\bm{h}}_{r}=\text{Re}(\hat{\bm{h}})$ and
$\hat{\bm{h}}_{m}=\text{Im}(\hat{\bm{h}})$. Consider real matrices

$\displaystyle\bar{\bm{A}}$
$\displaystyle=\left[\begin{array}[]{cc}\bar{h}_{r,1}&\bar{h}_{m,1}\\\
\vdots&\vdots\\\ \bar{h}_{r,4}&\bar{h}_{m,4}\end{array}\right],$
$\displaystyle\hat{\bm{A}}$
$\displaystyle=\left[\begin{array}[]{cc}\hat{h}_{r,1}&\hat{h}_{m,1}\\\
\vdots&\vdots\\\ \hat{h}_{r,4}&\hat{h}_{m,4}\end{array}\right],$ (82g)
$\displaystyle\bar{\bm{B}}$
$\displaystyle=\left[\begin{array}[]{cc}\bar{h}_{r,5}&\bar{h}_{m,5}\\\
\vdots&\vdots\\\ \bar{h}_{r,N}&\bar{h}_{m,N}\end{array}\right],$
$\displaystyle\hat{\bm{B}}$
$\displaystyle=\left[\begin{array}[]{cc}\hat{h}_{r,5}&\hat{h}_{m,5}\\\
\vdots&\vdots\\\ \hat{h}_{r,N}&\hat{h}_{m,N}\end{array}\right],$ (82n)

where $\bar{h}_{r,n}$, $\bar{h}_{m,n}$, $\hat{h}_{r,n}$ and $\hat{h}_{m,n}$
are the $n$-th elements of vectors $\bar{\bm{h}}_{r}$, $\bar{\bm{h}}_{m}$,
$\hat{\bm{h}}_{r}$ and $\hat{\bm{h}}_{m}$, respectively. Then, the matrix
$\mathcal{\bm{E}}_{r}$ can be equivalently written as

$\displaystyle\mathcal{\bm{E}}_{r}$
$\displaystyle=\bar{\bm{H}}_{r}-\hat{\bm{H}}_{r}$ (83a)
$\displaystyle=\left[\begin{array}[]{c}\bar{\bm{A}}\\\
\bar{\bm{B}}\end{array}\right]\left[\bar{\bm{A}}^{T},\bar{\bm{B}}^{T}\right]-\left[\begin{array}[]{c}\hat{\bm{A}}\\\
\hat{\bm{B}}\end{array}\right]\left[\hat{\bm{A}}^{T},\hat{\bm{B}}^{T}\right]$
(83f)
$\displaystyle=\left[\begin{array}[]{cc}\bar{\bm{A}}\bar{\bm{A}}^{T}-\hat{\bm{A}}\hat{\bm{A}}^{T}&\bar{\bm{A}}\bar{\bm{B}}^{T}-\hat{\bm{A}}\hat{\bm{B}}^{T}\\\
\bar{\bm{B}}\bar{\bm{A}}^{T}-\hat{\bm{B}}\hat{\bm{A}}^{T}&\bar{\bm{B}}\bar{\bm{B}}^{T}-\hat{\bm{B}}\hat{\bm{B}}^{T}\end{array}\right],$
(83i)

where $\bar{\bm{A}}\bar{\bm{A}}^{T}-\hat{\bm{A}}\hat{\bm{A}}^{T}$ is a
$4\times 4$ diagonal matrix and all other three blocks in (83i) are zeros,
according to equation (81) and the analysis above. Thus, we have
$\bar{\bm{A}}\bar{\bm{B}}^{T}=\hat{\bm{A}}\hat{\bm{B}}^{T}$ and
$\bar{\bm{B}}\bar{\bm{B}}^{T}=\hat{\bm{B}}\hat{\bm{B}}^{T}$. Obviously, the
space spanned by the columns of $\hat{\bm{B}}$ is the same as that of
$\bar{\bm{B}}$, which means the columns of $\hat{\bm{B}}$ can be represented
by $\bar{\bm{B}}$ with a coefficient matrix $\bm{Y}\in\mathbb{R}^{2\times 2}$,
i.e., $\hat{\bm{B}}=\bar{\bm{B}}\bm{Y}$. Thus,
$\bar{\bm{B}}\bar{\bm{B}}^{T}=\hat{\bm{B}}\hat{\bm{B}}^{T}=\bar{\bm{B}}\bm{Y}\bm{Y}^{T}\bar{\bm{B}}^{T}$.
Since $\bar{\bm{h}}$ represents the equivalent channel from the BS to user,
the matrix $\bar{\bm{B}}$ has full column rank $2$ with probability $1$ for
$N\geq 6$, and thus $\bar{\bm{B}}^{T}\bar{\bm{B}}\in\mathbb{R}^{2\times 2}$ is
invertible. Hence, $\bm{Y}\bm{Y}^{T}$ can be solved as

$\bm{Y}\bm{Y}^{T}=\left(\bar{\bm{B}}^{T}\bar{\bm{B}}\right)^{-1}\bar{\bm{B}}^{T}\left(\bar{\bm{B}}\bar{\bm{B}}^{T}\right)\bar{\bm{B}}\left(\bar{\bm{B}}^{T}\bar{\bm{B}}\right)^{-1}=\bm{I},$
(84)

which means that $\bm{Y}\in\mathbb{R}^{2\times 2}$ is an orthogonal matrix. On
the other hand, we have
$\bar{\bm{A}}\bar{\bm{B}}^{T}=\hat{\bm{A}}\hat{\bm{B}}^{T}=\hat{\bm{A}}\bm{Y}^{T}\bar{\bm{B}}^{T}$.
By applying $\bar{\bm{B}}$ on the left to both sides, it can be simplified as
$\bar{\bm{A}}=\hat{\bm{A}}\bm{Y}^{T}$, or equivalently,
$\hat{\bm{A}}=\bar{\bm{A}}\bm{Y}$. Then,
$\bar{\bm{A}}\bar{\bm{A}}^{T}=\hat{\bm{A}}\bm{Y}^{T}\bm{Y}\hat{\bm{A}}^{T}=\hat{\bm{A}}\hat{\bm{A}}^{T}$
holds and $\mathcal{\bm{E}}_{r}=\bm{0}$ follows, which leads to
$\bar{\bm{H}}_{r}=\hat{\bm{H}}_{r}$.

Meanwhile, it is worth noting that

$\bar{\bm{h}}=\left[\bar{\bm{h}}_{r},\bar{\bm{h}}_{m}\right]\left[\begin{array}[]{c}1\\\
j\end{array}\right]=\left[\begin{array}[]{c}\bar{\bm{A}}\\\
\bar{\bm{B}}\end{array}\right]\left[\begin{array}[]{c}1\\\
j\end{array}\right].$ (85)

Therefore, the matrix $\bar{\bm{H}}$ can be written as

$\bar{\bm{H}}=\bar{\bm{h}}\bar{\bm{h}}^{H}=\left[\begin{array}[]{c}\bar{\bm{A}}\\\
\bar{\bm{B}}\end{array}\right]\left[\begin{array}[]{cc}1&-j\\\
j&1\end{array}\right]\left[\bar{\bm{A}}^{T},\bar{\bm{B}}^{T}\right].$ (86)

Similarly, $\hat{\bm{H}}$ can be equivalently written as

$\displaystyle\hat{\bm{H}}$
$\displaystyle=\left[\begin{array}[]{c}\hat{\bm{A}}\\\
\hat{\bm{B}}\end{array}\right]\left[\begin{array}[]{cc}1&-j\\\
j&1\end{array}\right]\left[\hat{\bm{A}}^{T},\hat{\bm{B}}^{T}\right]$ (87e)
$\displaystyle=\left[\begin{array}[]{c}\bar{\bm{A}}\\\
\bar{\bm{B}}\end{array}\right]\bm{Y}\left[\begin{array}[]{cc}1&-j\\\
j&1\end{array}\right]\bm{Y}^{T}\left[\bar{\bm{A}}^{T},\bar{\bm{B}}^{T}\right].$
(87j)

Define matrix $\bm{J}\in\mathbb{R}^{2\times 2}$ as

$\bm{J}=\bm{Y}\left[\begin{array}[]{cc}1&-j\\\
j&1\end{array}\right]\bm{Y}^{T}.$ (88)

Let $\bm{Y}=[\bm{y}_{1},\bm{y}_{2}]$ with $\bm{y}_{1}=[y_{11},y_{12}]^{T}$ and
$\bm{y}_{2}=[y_{21},y_{22}]^{T}$. Then, we have
$\bm{y}_{1}\bm{y}_{1}^{T}+\bm{y}_{2}\bm{y}_{2}^{T}=\bm{Y}\bm{Y}^{T}=\bm{I}$
and

$\displaystyle\bm{J}$
$\displaystyle=\left(\bm{y}_{1}+j\bm{y}_{2}\right)\left(\bm{y}_{1}^{T}-j\bm{y}_{2}^{T}\right)$
(89a)
$\displaystyle=\left(\bm{y}_{1}\bm{y}_{1}^{T}+\bm{y}_{2}\bm{y}_{2}^{T}\right)+j\left(\bm{y}_{2}\bm{y}_{1}^{T}-\bm{y}_{1}\bm{y}_{2}^{T}\right)$
(89b) $\displaystyle=\bm{I}+j\left[\begin{array}[]{cc}0&-\eta\\\
\eta&0\end{array}\right],$ (89e)

where $\eta=y_{11}y_{22}-y_{12}y_{21}=\text{det}(\bm{Y})=\pm 1$ because
$\bm{Y}$ is an orthogonal matrix. If $\text{det}(\bm{Y})=1$, we have

$\bm{J}=\left[\begin{array}[]{cc}1&-j\\\ j&1\end{array}\right],\
\hat{\bm{H}}=\left[\begin{array}[]{c}\bar{\bm{A}}\\\
\bar{\bm{B}}\end{array}\right]\bm{J}\left[\bar{\bm{A}}^{T},\bar{\bm{B}}^{T}\right]=\bar{\bm{H}}.$
(90)

If $\text{det}(\bm{Y})=-1$, however, $\bm{J}$ and $\hat{\bm{H}}$ are given as

$\bm{J}=\left[\begin{array}[]{cc}1&j\\\ -j&1\end{array}\right],\
\hat{\bm{H}}=\left[\begin{array}[]{c}\bar{\bm{A}}\\\
\bar{\bm{B}}\end{array}\right]\bm{J}\left[\bar{\bm{A}}^{T},\bar{\bm{B}}^{T}\right]=\bar{\bm{H}}^{*}.$
(91)

Therefore, $\bar{\bm{H}}$ and $\bar{\bm{H}}^{*}$ are the only two solutions to
Problem (4) with probability $1$ for $b=1$ when $N\geq 6$ and
$D_{V}=\mathcal{D}_{N}^{(1)}$. Hereto, the proof for the existence and
uniqueness of the solutions to Problem (4) claimed in Proposition 1 have been
completed.

## References

  * [1] G. Yan, L. Zhu, and R. Zhang, “Channel autocorrelation estimation for IRS-aided wireless communications based on power measurements,” in _2023 IEEE Globecom Workshops (GC Wkshps)_ , Kuala Lumpur, Malaysia, Dec. 2023, pp. 1457–1462.
  * [2] Q. Wu, B. Zheng, C. You, L. Zhu, K. Shen, X. Shao, W. Mei, B. Di, H. Zhang, E. Basar, L. Song, M. D. Renzo, Z.-Q. Luo, and R. Zhang, “Intelligent surfaces empowered wireless network: Recent advances and the road to 6G,” arXiv:2312.16918, 2023. [Online]. Available: https://arxiv.org/abs/2312.16918.
  * [3] Q. Wu, S. Zhang, B. Zheng, C. You, and R. Zhang, “Intelligent reflecting surface-aided wireless communications: A tutorial,” _IEEE Trans. Commun._ , vol. 69, no. 5, pp. 3313–3351, May 2021.
  * [4] G. Yan, L. Zhu, and R. Zhang, “Passive reflection optimization for IRS-aided multicast beamforming with discrete phase shifts,” _IEEE Wireless Commun. Lett._ , vol. 12, no. 8, pp. 1424–1428, Aug. 2023.
  * [5] W. Ma, L. Zhu, and R. Zhang, “Passive beamforming for 3-D coverage in IRS-assisted communications,” _IEEE Wireless Commun. Lett._ , vol. 11, no. 8, pp. 1763–1767, Aug. 2022.
  * [6] B. Zheng, C. You, W. Mei, and R. Zhang, “A survey on channel estimation and practical passive beamforming design for intelligent reflecting surface aided wireless communications,” _IEEE Commun. Surveys Tuts._ , vol. 24, no. 2, pp. 1035–1071, 2nd quart. 2022.
  * [7] M. Najafi, V. Jamali, R. Schober, and H. V. Poor, “Physics-based modeling and scalable optimization of large intelligent reflecting surfaces,” _IEEE Trans. Commun._ , vol. 69, no. 4, pp. 2673–2691, Apr. 2021.
  * [8] E. Björnson and L. Sanguinetti, “Power scaling laws and near-field behaviors of massive MIMO and intelligent reflecting surfaces,” _IEEE Open J. Commun. Soc._ , vol. 1, pp. 1306–1324, Sep. 2020.
  * [9] D. Mishra and H. Johansson, “Channel estimation and low-complexity beamforming design for passive intelligent surface assisted MISO wireless energy transfer,” in _Proc. IEEE Int. Conf. Acoust. Speech Signal Process. (ICASSP)_ , Brighton, UK, May 2019, pp. 4659–4663.
  * [10] Z.-Q. 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, Feb. 2020.
  * [11] T. L. Jensen and E. De Carvalho, “An optimal channel estimation scheme for intelligent reflecting surfaces based on a minimum variance unbiased estimator,” in _Proc. IEEE Int. Conf. Acoust. Speech Signal Process. (ICASSP)_ , Barcelona, Spain, May 2020, pp. 5000–5004.
  * [12] B. Zheng and R. Zhang, “Intelligent reflecting surface-enhanced OFDM: Channel estimation and reflection optimization,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 4, pp. 518–522, Apr. 2020.
  * [13] C. You, B. Zheng, and R. Zhang, “Channel estimation and passive beamforming for intelligent reflecting surface: Discrete phase shift and progressive refinement,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 11, pp. 2604–2620, Nov. 2020.
  * [14] B. Zheng, C. You, and R. Zhang, “Intelligent reflecting surface assisted multi-user OFDMA: Channel estimation and training design,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 12, pp. 8315–8329, Dec. 2020.
  * [15] 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, July 2020.
  * [16] P. Wang, J. Fang, H. Duan, and H. Li, “Compressed channel estimation for intelligent reflecting surface-assisted millimeter wave systems,” _IEEE Signal Process. Lett._ , vol. 27, pp. 905–909, May 2020.
  * [17] X. Zheng, P. Wang, J. Fang, and H. Li, “Compressed channel estimation for IRS-assisted millimeter wave OFDM systems: A low-rank tensor decomposition-based approach,” _IEEE Wireless Commun. Lett._ , vol. 11, no. 6, pp. 1258–1262, June 2022.
  * [18] J. He, H. Wymeersch, and M. Juntti, “Channel estimation for RIS-aided mmwave MIMO systems via atomic norm minimization,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 9, pp. 5786–5797, Sep. 2021.
  * [19] X. Wei, D. Shen, and L. Dai, “Channel estimation for RIS assisted wireless communications—part II: An improved solution based on double-structured sparsity,” _IEEE Commun. Lett._ , vol. 25, no. 5, pp. 1403–1407, May 2021\.
  * [20] C. Liu, X. Liu, D. W. K. Ng, and J. Yuan, “Deep residual learning for channel estimation in intelligent reflecting surface-assisted multi-user communications,” _IEEE Trans. Wireless Commun._ , vol. 21, no. 2, pp. 898–912, Aug. 2022.
  * [21] V. Jamali, M. Najafi, R. Schober, and H. V. Poor, “Power efficiency, overhead, and complexity tradeoff of IRS codebook design—quadratic phase-shift profile,” _IEEE Commun. Lett._ , vol. 25, no. 6, pp. 2048–2052, June 2021\.
  * [22] W. R. Ghanem, V. Jamali, M. Schellmann, H. Cao, J. Eichinger, and R. Schober, “Optimization-based phase-shift codebook design for large IRSs,” _IEEE Commun. Lett._ , vol. 27, no. 2, pp. 635–639, Feb. 2023.
  * [23] 3GPP, “NR - physical channels and modulation,” 3rd Generation Partnership Project (3GPP), Technical Specification (TS) 38.211, 2018. [Online]. Available: https://www.3gpp.org.
  * [24] S. Ren, K. Shen, Y. Zhang, X. Li, X. Chen, and Z.-Q. Luo, “Configuring intelligent reflecting surface with performance guarantees: Blind beamforming,” _IEEE Trans. Wireless Commun._ , vol. 22, no. 5, pp. 3355–3370, May 2023.
  * [25] V. Arun and H. Balakrishnan, “RFocus: Beamforming using thousands of passive antennas,” _USENIX Symp. Netw. Sys. Design Implementation (NSDI)_ , pp. 1047–1061, Feb. 2020.
  * [26] 3GPP, “Evolved universal terrestrial radio access (E-UTRA); physical layer; measurements,” 3rd Generation Partnership Project (3GPP), Technical Specification (TS) 36.214, 2022. [Online]. Available: https://www.3gpp.org.
  * [27] E. Candes, T. Strohmer, and V. Voroninski, “Phaselift: Exact and stable signal recovery from magnitude measurements via convex programming,” _Commun. Pure Appl. Math._ , vol. 66, no. 8, pp. 1241–1274, Aug. 2013.
  * [28] L. Demanet and P. Hand, “Stable optimizationless recovery from phaseless linear measurements,” _J. Fourier Anal. Appl._ , vol. 20, pp. 199–221, Aug. 2012. [Online]. Available: https://api.semanticscholar.org/CorpusID:17332752.
  * [29] K. W. Choi, D. I. Kim, and M. Y. Chung, “Received power-based channel estimation for energy beamforming in multiple-antenna rf energy transfer system,” _IEEE Trans. Signal Process._ , vol. 65, no. 6, pp. 1461–1476, March 2017.
  * [30] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1,” http://cvxr.com/cvx, Mar. 2014.
  * [31] Z.-Q. Luo, W.-K. Ma, A. M.-C. So, Y. Ye, and S. Zhang, “Semidefinite relaxation of quadratic optimization problems,” _IEEE Signal Process. Mag._ , vol. 27, no. 3, pp. 20–34, May 2010.
  * [32] S. Andrilli and D. Hecker, “Chapter 9 - numerical techniques,” in _Elementary Linear Algebra (Sixth Edition)_ , sixth edition ed., S. Andrilli and D. Hecker, Eds. Academic Press, 2023, pp. 401–439. [Online]. Available: https://www.sciencedirect.com/science/article/pii/B9780128229781000195.
  * [33] 3GPP, “Evolved universal terrestrial radio access (E-UTRA); requirements for support of radio resource management,” 3rd Generation Partnership Project (3GPP), Technical Specification (TS) 36.133, 2023. [Online]. Available: https://www.3gpp.org.
  * [34] M. A. Woodbury, _Inverting modified matrices_. Princeton University, Princeton, NJ, 1950, statistical Research Group, Memo. Rep. no. 42,.
  * [35] S. Ren, K. Shen, X. Li, X. Chen, and Z.-Q. Luo, “A linear time algorithm for the optimal discrete IRS beamforming,” _IEEE Wireless Commun. Lett._ , vol. 12, no. 3, pp. 496–500, Mar. 2023.