# Channel Estimation in RIS-Assisted MIMO Systems Operating Under
Imperfections

Paulo R. B. Gomes, Gilderlan T. de Araújo, Bruno Sokal, André L. F. de
Almeida, Behrooz Makki, and Gábor Fodor Paulo R. B. Gomes, Gilderlan T. de
Araújo, Bruno Sokal, and André L. F. de Almeida are with the Wireless Telecom
Research Group (GTEL), Department of Teleinformatics Engineering, Federal
University of Ceará, Fortaleza-CE. E-mails:
<EMAIL_ADDRESS>Makki is with Ericsson
Research, G$\ddot{\text{o}}$teborg, Sweden. E-mail:
<EMAIL_ADDRESS>Fodor is with Ericsson Research and KTH Royal
Institute of Technology, Stockholm, Sweden. E-mail:
<EMAIL_ADDRESS>work was supported by the Ericsson Research,
Sweden, and Ericsson Innovation Center, Brazil, under UFC.48 Technical
Cooperation Contract Ericsson/UFC. This study was financed in part by the
Coordenação de Aperfeiçoamento de Pessoal de Nível Superior - Brasil
(CAPES)-Finance Code 001, and CAPES/PRINT Proc. 88887.311965/2018-00. André L.
F. de Almeida acknowledges CNPq for its financial support under the grant
312491/2020-4. Gábor Fodor was partially supported by the Digital Futures
project PERCy.Part of this work has been submitted for possible presentation
in IEEE GLOBECOM 2022 [1].

###### Abstract

Reconfigurable intelligent surface (RIS) is a potential technology component
of future wireless networks due to its capability of shaping the wireless
environment. The promising gains of RIS-assisted multiple-input multiple-
output (MIMO) systems in terms of extended coverage and enhanced capacity are,
however, critically dependent on the accuracy of the channel state
information. However, traditional channel estimation schemes are not
applicable in RIS-assisted MIMO networks, since passive RISs typically lack
the signal processing capabilities that are assumed by channel estimation
algorithms. This becomes most problematic when physical imperfections or
electronic impairments affect the RIS due to its exposition to different
environmental effects or caused by hardware limitations from the circuitry.
While these real-world effects are typically ignored in the literature, in
this paper we propose efficient channel estimation schemes for RIS-assisted
MIMO systems taking different imperfections into account. Specifically, we
propose two sets of tensor-based algorithms, based on the parallel factor
analysis decomposition schemes. First, by assuming a long-term model – in
which the RIS imperfections, modeled as unknown phase shifts, are static
within the channel coherence time – we formulate an iterative alternating
least squares (ALS)-based algorithm for the joint estimation of the
communication channels and the unknown phase deviations. Next, we develop the
short-term imperfection model, which allows both amplitude and phase RIS
imperfections to be non-static with respect to the channel coherence time. We
propose two iterative ALS-based and closed-form higher order singular value
decomposition-based algorithms for the joint estimation of the channels and
the unknown impairments. Moreover, we analyze the identifiability and
computational complexity of the proposed algorithms and study the effects of
various imperfections on the channel estimation quality. Simulation results
demonstrate the effectiveness of the proposed tensor-based algorithms in terms
of estimation accuracy and computational complexity.

###### Index Terms:

Alternating least squares (ALS), channel estimation, hardware impairments,
higher-order singular value decomposition (HOSVD), imperfection detection,
multiple input multiple output (MIMO), parallel factor analysis (PARAFAC),
reconfigurable intelligent surface (RIS).

##  I Introduction

Wireless communications have become a necessity in our daily lives providing
significant benefits to individuals, businesses and the society at large. The
steadily increasing demands for ubiquitous wireless services drive the efforts
by the research and standardization communities to improve coverage, system
capacity as well as the reliability and quality of a growing number of
applications [2, 3]. The continuous growth of the number of mobile
subscriptions, devices and traffic increases the number of deployed
infrastructure nodes, which makes capital and operational expenditures as well
as energy consumption challenging for mobile network operators. [4].
Therefore, to guarantee green sustainable wireless networks, one needs to
carefully take the energy consumption and the hardware costs into account [5,
6].

Thanks to the recent development of meta-materials, reconfigurable intelligent
surface (RIS) has emerged as a potential technology applicable in future
wireless networks through the novel smart and programming environment paradigm
[7]. Unlike conventional networks, RIS-assisted networks enable the system to
shape the wireless environment to be more suitable for wireless communication
[8, 9]. In conventional networks, the wireless transmission medium is seen as
an uncontrollable element in the system due to the randomness in the radio
environment. Therefore, the propagation of electromagnetic waves through the
wireless channel cannot be controlled after they are emitted from the
transmitters and before they reach the receivers [10]. An RIS is a planar
meta-surface that has a large number of low-cost passive reflecting elements
with adjustable capability of the reflected signal parameters such as
amplitude, phase, frequency and polarization. With a well-configured RIS, the
propagation conditions are improved by controlling the scattering
characteristics to create passive and active beamforming at the RIS and the
wireless transceivers, respectively, to achieve high beamforming gains [11].
In addition, a passive RIS does not require active radio frequency chains for
signal transmission, reception and processing, since it simply relies on
passive signal reflection, which makes it a power-efficient and less-
complexity technology [12, 13, 14, 15].

Despite its potential benefits, the performance gains achieved in RIS-assisted
systems are strongly dependent on, among others factors, the quality of the
channel state information (CSI), where the channel between the wireless
transceivers is effectively split into two-hop RIS-assisted channels. This is
because CSI is required to jointly design the passive beamforming at the RIS
and the active beamforming at the transceivers [16]. However, due to the
passive nature of the RIS, channel estimation is not performed at the RIS but
only on gNB or UE side.

Recognizing this issue, recent works propose different strategies to tackle
the channel estimation problem in RIS-assisted wireless communications under
various configurations (see, e.g., [17, 18, 19] and references therein). The
most popular approach in the literature, known as the on/off method, consists
of switching groups of RIS elements (or each individual element) to estimate
the associated channel coefficients sequentially [20]. More specifically, each
element of the cascaded channel is estimated at each time slot from pilots
sent by the transmitter and reflected by the corresponding RIS element e.g.,
[21].

In the context of RIS-assisted multiple-input multiple-output (MIMO) systems,
[22] proposes a two-stage algorithm including sparse matrix factorization and
matrix completion to estimate the channels in massive array setups. In [23],
the authors leverage compressed sensing (CS) and deep learning tools to reduce
the training overhead. In a similar way, [24] also uses CS techniques and
proposes a two-step iterative solution to solve the problem in a multi-user
scenario. Then, [25] proposes an iterative method to estimate the channel
parameters (angle of departure, angle of arrival and propagation path gains)
in RIS-assisted millimeter wave (mmWave) MIMO systems, resulting in limited
training overhead.

Over the past two decades, tensor algebra and higher-order tensor
decompositions have been successfully applied to signal modeling and
processing in wireless communications. Tensor-based signal processing exploits
the intrinsic multidimensional structure of wireless channels and signals to
achieve an improved estimation accuracy under more relaxed and flexible system
parameter settings, compared to conventional matrix-based approaches [26, 27,
28].

Recently, tensor approaches have been proposed in the context of RIS-assisted
wireless communications. Particularly, [29] capitalizes on the parallel factor
(PARAFAC) decomposition to formulate an efficient iterative algorithm based on
the alternating least squares (ALS) concept to solve the channel estimation
problem in the downlink of a multi-user multiple input single output (MISO)
network. Also, in [30] and [31], the authors develop simple iterative and
closed-form channel estimation algorithms based on a PARAFAC modeling of the
single-user MIMO scenario.

Figure 1: An illustration of a RIS-assisted MIMO wireless communication system
operating under imperfections.

The main motivation for RIS is complexity/energy reduction, otherwise
integrated access and backhaul (IAB) nodes, relays, and/or repeaters may be
alternatives. Then, complexity reduction may affect the hardware capability
and, consequently, its performance. Therefore, to have a realistic view on the
performance and usefulness of RIS, one needs to carefully take the hardware
impairments as well as the environmental effects, such as water
precipitations, flying debris, air particles, snowflakes, freezing rain,
sleet, dry/damp sand particles and dirt into account, as shown in Fig. 1.

In practice, such blocking objects as well as the hardware impairments induce
unwanted attenuation and phase shift on the reflected signals by the RIS
introducing static or time-varying distortions in the received signal, which
directly affect the channel estimation accuracy and, consequently, the system
performance. However, the effects of hardware impairments and environmental
imperfections have been rarely studied in the literature. For instance, [32,
33, 34, 35], consider RISs operating under finite resolution of the phase
shifts or phase estimation errors from imperfect channel estimation. Then,
[36] and [37] consider different environmental effects on the RIS and propose
methods to jointly estimate the channel and array blockage parameters in
mmWave RIS-assisted systems.

As we explain in the following, typical channel estimation methods may not be
able to deal with different imperfections and may fail to properly estimate
the channel. Therefore, it is necessary to continuously monitor the channel
and compensate for imperfections in order to maintain robust system operation.

In this paper, we propose tensor-based algorithms for the joint estimation of
the involved channels and imperfections in RIS-assisted MIMO systems. We take
both the long- and short-term imperfections into account. First, we show that
the received signal under the long- and short-term imperfection models can be
recast as tensors following trilinear and quadrilinear PARAFAC models,
respectively. Exploiting the multi-linear structure of these models, we derive
two sets of tensor-based algorithms. For the long-term imperfection (LTI)
model, where the RIS imperfections, modeled as unknown phase shifts, are
static within the channel coherence time, we formulate an iterative trilinear
ALS-based algorithm, named TALS-LTI, for the joint estimation of the involved
channels and the unknown RIS phase deviations. Next, we generalize the
imperfections behavior to be non-static with respect to the channel coherence
time, referred to as the short-term imperfection (STI) model. For such a more
challenging scenario, we propose iterative and closed-form tensor
decomposition-based algorithms named TALS-STI and HOSVD-STI, respectively, to
solve the joint channel and RIS imperfections estimation. We also study the
identifiability of the proposed estimators, discuss their computational
complexity and investigate the effect of imperfections on the network
performance. The key features of the proposed tensor-based algorithms are
their ability to properly estimate the channel and their robustness to
different kinds of real-world imperfections at the RIS.

The simulation results show that, compared to the state-of-the-art methods,
the proposed algorithms properly estimate the involved channels when different
kinds of imperfections are takes into account. As a example, in the high
signal-to-noise ratio (SNR) regime, our proposed algorithms improve the
channel estimation by approximately $100$x compared to the method of [19],
while present performance close to the lower-bound least squares (LS)
estimator. Also, our proposed algorithms reduce considerably the overall
computational complexity, compared to the related state-of-the-art method.
Finally, the proposed TALS-LTI and TALS-STI algorithms are more flexible for
the choices of training parameters compared to the proposed HOSVD-STI
algorithm. Thus, the TALS-LTI and TALS-STI algorithms are preferable when more
flexible choices for training parameters are required, while the HOSVD-STI is
preferred when low processing delay is desired.

The rest of this paper is organized as follows. In Section II, the signal
model of the RIS-assisted MIMO communication system operating under
imperfections is introduced. We distinguish between two different types of
imperfections under the RIS operation, namely, LTI and STI. Then, the channel
estimation problem is discussed for the LTI and STI scenarios. These signal
models are reformulated as higher-order tensors that follow trilinear and
quadrilinear PARAFAC models, from which two sets of iterative and closed-form
tensor-based algorithms for the joint estimation of the involved channels and
the RIS imperfections are developed in Section III. A detailed identifiability
analysis and its link to the system design recommendations as well as the
computational complexity of the proposed tensor-based algorithms are provided
in Section IV. Simulation results are presented in Section V. Finally,
conclusions are drawn in Section VI.

###  I-A Notations and Properties

The notation conventions and the properties that will be used throughout this
paper are defined in the following. Scalars are denoted by lower-case letters
($a,b,\ldots$), column vectors by bold lower-case letters
($\mathbf{a},\mathbf{b},\ldots$), matrices by bold upper-case letters
($\mathbf{A},\mathbf{B},\ldots$) and tensors are represented by upper-case
calligraphic letters ($\mathcal{A},\mathcal{B},\ldots$). Then,
$\mathbf{A}^{\text{T}}$ and $\mathbf{A}^{{\dagger}}$ stand for the transpose
and Moore-Penrose pseudo-inverse of $\mathbf{A}$, respectively. The operator
$\text{vec}(\cdot)$ vectorizes its matrix argument by stacking its columns on
top of each other, while $\text{vecd}(\cdot)$ forms a vector out of the
diagonal of its matrix argument. Also, $\|\cdot\|_{\text{F}}$ represents the
Frobenius norm of a matrix or a tensor, which is defined as the square root of
the sum of the squared of its elements. $\lceil x\rceil$ is equal to the
smallest integer that is greater than or equal to $x$. Moreover,
$\mathbf{I}_{M}$ is the $M\times M$ identity matrix and $j=\sqrt{-1}$ is the
imaginary unit. The operator $\mathbf{D}_{i}\left(\mathbf{A}\right)$ forms a
diagonal matrix from the $i$-th row of its matrix argument $\mathbf{A}$, while
the operator $\text{diag}(\mathbf{a})$ forms a diagonal matrix out of its
vector argument $\mathbf{a}$. We define the Kronecker, Hadamard (element-wise
product) and the outer product operators by $\otimes$, $\odot$ and $\circ$,
respectively. The Khatri-Rao product (column-wise Kronecker product) between
two matrices is defined as

$\mathbf{A}\diamond\mathbf{B}=\left[\mathbf{a}_{1}\otimes\mathbf{b}_{1},\ldots,\mathbf{a}_{Q}\otimes\mathbf{b}_{Q}\right]\in\mathbb{C}^{IJ\times
Q},$ (1)

or, equivalently,

$\mathbf{A}\diamond\mathbf{B}=\left[\mathbf{D}_{1}\left(\mathbf{A}\right)\mathbf{B}^{\text{T}},\ldots,\mathbf{D}_{Q}\left(\mathbf{A}\right)\mathbf{B}^{\text{T}}\right]^{\text{T}},$
(2)

where
$\mathbf{A}=\left[\mathbf{a}_{1},\ldots,\mathbf{a}_{Q}\right]\in\mathbb{C}^{I\times
Q}$ and
$\mathbf{B}=\left[\mathbf{b}_{1},\ldots,\mathbf{b}_{Q}\right]\in\mathbb{C}^{J\times
Q}$.

We shall make use of the following properties of the Khatri-Rao and Kronecker
products

$\displaystyle\text{vec}\left(\mathbf{A}\text{diag}\left(\mathbf{c}\right)\mathbf{B}^{\text{T}}\right)$
$\displaystyle=$
$\displaystyle\left(\mathbf{B}\diamond\mathbf{A}\right)\mathbf{c},\forall\mathbf{A},\mathbf{B},\mathbf{c},$
(3) $\displaystyle\mathbf{a}\otimes\mathbf{b}\otimes\mathbf{c}$
$\displaystyle=$
$\displaystyle\text{vec}\left(\mathbf{c}\circ\mathbf{b}\circ\mathbf{a}\right),\forall\mathbf{a},\mathbf{b},\mathbf{c}.$
(4)

Furthermore, the definitions and operations involving tensors are in
accordance with [38] and [39]. The $n$-mode unfolding matrix of $\mathcal{A}$
along its $n$-th mode (or dimension) is represented by
$\left[\mathcal{A}\right]_{(n)}$. The $n$-mode product between $\mathcal{A}$
and $\mathbf{B}$, returns a tensor
$\mathcal{C}=\mathcal{A}\times_{n}\mathbf{B}$ such that
$\left[\mathcal{C}\right]_{(n)}=\mathbf{B}\left[\mathcal{A}\right]_{(n)}$.

Figure 2: An illustration of the considered channel estimation protocol. The
training time is divided into $K$ time-blocks of duration $T$ symbol periods
each. The RIS activation pattern $\mathbf{s}[k]$ is fixed during the $k$-th
time-block and varies between different time-blocks while the pilot symbols
$\mathbf{x}[1],\ldots,\mathbf{x}[T]$ are reused from block-to-block.

##  II Signal Model and Problem Description

In this section, we first introduce the signal model and describe in detail
the two considered LTI and STI models.

###  II-A RIS Operating under LTI

We consider a single-user narrowband RIS-assisted MIMO communication system,
in which the transmitter (Tx) and the receiver (Rx) are equipped with arrays
composed of $M$ and $L$ antennas, respectively. Without loss of generality,
although this work assumes a downlink communication, our signal models also
apply to the uplink case by inverting the roles of the transmitter and the
receiver. To improve the communication performance, an RIS with $N$
individually adjustable passive reflecting elements is deployed in a proper
place that creates an alternative Tx-RIS-Rx link. The direct Tx-Rx link is
assumed to be too weak or unavailable due to unfavorable propagation
conditions. The considered system setup is illustrated in Fig. 1. We also
assume a block-fading channel where the Tx-RIS and RIS-Rx channels remain
constant for at least $k=1,\ldots,K,$ time-blocks, each with duration of $T$
symbols, representing a total duration of $KT$ symbol periods dedicated for
channel estimation within the channel coherence time.

Let us define the RIS activation pattern
$\mathbf{s}\left[k\right]\in\mathbb{C}^{N\times 1}$ configured at the $k$-th
time-block as

$\mathbf{s}\left[k\right]=\left[\beta_{1,k}e^{j\phi_{1,k}},\ldots,\beta_{N,k}e^{j\phi_{N,k}}\right]^{\text{T}}\in\mathbb{C}^{N\times
1},$ (5)

where $0\leq\phi_{n,k}\leq 2\pi$ and $0\leq\beta_{n,k}\leq 1$ denote the phase
shift and the amplitude reflection coefficient of the $n$-th RIS element tuned
at the $k$-th time-block $\forall n=1,\ldots,N$, and $\forall k=1,\ldots,K$,
respectively. As a protocol for channel estimation we assume the following
(see Fig. 2):

1) At each time-block $k=1,\ldots,K$, of duration $T$, the elements of
$\mathbf{s}\left[k\right]$ are dynamically tuned in a passive way via the
smart controller;

2) The activation pattern $\mathbf{s}\left[k\right]$ remains constant within
the $k$-th time-block but may vary between different time-blocks, yielding a
total of $\mathbf{s}\left[1\right],\ldots,\mathbf{s}\left[K\right]$ different
adjustable patterns to the RIS during the channel estimation stage;

3) The pilot symbol $\mathbf{x}\left[t,k\right]\in\mathbb{C}^{M\times 1}$
transmitted at the $t$-th symbol period within the $k$-th time-block is reused
for each $k=1,\ldots,K$, i.e.,
$\mathbf{x}\left[t,k\right]=\mathbf{x}\left[t\right]$ $\forall k=1,\ldots,K$.

The baseband received pilot signal
$\mathbf{y}\left[t,k\right]\in\mathbb{C}^{L\times 1}$ associated with the
$t$-th symbol period at the $k$-th time-block can be expressed as

$\mathbf{y}\left[t,k\right]=\mathbf{G}\text{diag}\left(\mathbf{s}\left[k\right]\right)\mathbf{H}^{\text{T}}\mathbf{x}\left[t\right]+\mathbf{v}\left[t,k\right].$
(6)

Collecting the received signals during the $T$ symbol periods at the $k$-th
time-block, the model in (6) can be rewritten as

$\mathbf{Y}\left[k\right]=\mathbf{G}\text{diag}\left(\mathbf{s}\left[k\right]\right)\mathbf{H}^{\text{T}}\mathbf{X}+\mathbf{V}\left[k\right]\in\mathbb{C}^{L\times
T},$ (7)

where
$\mathbf{Y}\left[k\right]=\left[\mathbf{y}\left[1,k\right],\ldots,\mathbf{y}\left[T,k\right]\right]\in\mathbb{C}^{L\times
T}$. The matrices $\mathbf{H}\in\mathbb{C}^{M\times N}$ and
$\mathbf{G}\in\mathbb{C}^{L\times N}$ denote the Tx-RIS and RIS-Rx channels,
respectively, while
$\mathbf{X}=\left[\mathbf{x}\left[1\right],\ldots,\mathbf{x}\left[T\right]\right]\in\mathbb{C}^{M\times
T}$ collects the pilot signals transmitted within the $k$-th time-block, and
$\mathbf{V}\left[k\right]=\left[\mathbf{v}\left[1,k\right],\ldots,\mathbf{v}\left[T,k\right]\right]\in\mathbb{C}^{L\times
T}$ is the additive white Gaussian noise (AWGN) matrix with zero mean and unit
variance elements. In order to simplify our formulation and analysis, without
loss of generality, we assume the transmission of the pilot signal
$\mathbf{X}=\mathbf{I}_{M}$.

In this work, a special attention is given to the structure of the RIS
activation pattern. In practice, some imperfections at the RIS elements are
common to occur. Initially, we assume the case in which such imperfections
induce long-term static phase shift perturbations at the RIS response. Such
imperfections may come from, e.g., phase noise due to the finite resolution of
the phase shifts or by phase estimation errors from imperfect channel
estimation. In the presence of these imperfections, the structure of the RIS
activation pattern in (5) is modified in an undesired manner leading to the
following resulting RIS reflection pattern that incorporates the imperfection
contributions:

$\bar{\mathbf{s}}\left[k\right]=\left[\beta_{1,k}e^{j(\phi_{1,k}+\theta_{1})},\ldots,\beta_{N,k}e^{j(\phi_{N,k}+\theta_{N})}\right]^{\text{T}}.$
(8)

Here, $0\leq\theta_{n}\leq 2\pi$ $\forall n=1,\ldots,N$ denotes the phase
shift perturbation that affects the $n$-th RIS element. Also, (8) can be
alternatively represented in a more attractive form for our formulation as

$\bar{\mathbf{s}}\left[k\right]=\mathbf{e}\odot\mathbf{s}\left[k\right]\in\mathbb{C}^{N\times
1},$ (9)

where the entries of the random vector $\mathbf{e}\in\mathbb{C}^{N\times 1}$
that collects all unknown existing phase perturbations are defined as

$e_{n}=\left\\{\begin{array}[]{ll}1,&\quad\text{{\color[rgb]{0,0,0}non-
impaired case}}\\\
e^{j\theta_{n}},&\quad\text{{\color[rgb]{0,0,0}otherwise}},\end{array}\right.$
(10)

for $n=1,\ldots,N$. Making use of definitions in (9) and (10), the impaired
version of the received signal at the Rx node in (7) can be expressed as

$\mathbf{Y}\left[k\right]=\mathbf{G}\text{diag}\left(\underbrace{\mathbf{e}\odot\mathbf{s}\left[k\right]}_{\bar{\mathbf{s}}\left[k\right]}\right)\mathbf{H}^{\text{T}}+\mathbf{V}\left[k\right]\in\mathbb{C}^{L\times
M},$ (11)

$\forall k=1,\ldots,K$. Equivalently, in matrix form we have

$\mathbf{Y}\left[k\right]=\mathbf{G}\text{diag}\left(\mathbf{e}\right)\mathbf{D}_{k}\left(\mathbf{S}\right)\mathbf{H}^{\text{T}}+\mathbf{V}\left[k\right],$
(12)

where
$\mathbf{S}=\left[\mathbf{s}\left[1\right],\ldots,\mathbf{s}\left[K\right]\right]^{\text{T}}\in\mathbb{C}^{K\times
N}$ collects in its rows the RIS activation patterns used accross $K$ time-
blocks specifically configured for the channel estimation.

Particularly, the impaired received signal in (12) considers an RIS operating
under LTI that induces phase shift perturbations in the reflected signals. In
other words, in this model we assume that the vector
$\mathbf{e}\in\mathbb{C}^{N\times 1}$ is formed only by phase components that
remain static within the $KT$ symbol periods. Figure 3(a) illustrates the
considered LTI model. This occurs, for instance, when the behavior of the
imperfections at the RIS elements are static compared to the channel coherence
time.

(a) An illustration of the LTI model. The vector $\mathbf{e}$ with phase
perturbations is static during the $KT$ symbol periods.

(b) An illustration of the generalized STI model with time-varying RIS
fluctuations during the channel coherence time.

Figure 3: Illustration of the LTI and STI models considered in the paper.

###  II-B RIS Operating under STI

Here, we assume that the channel estimation occurs after the receiver collects
$p=1,\ldots,P$, frames composed of $KT$ symbol periods each. The division of
the reception time into $P$ frames is motivated by possible short-term
variations caused by the RIS imperfections, i.e., the imperfections changing
more frequently. In this approach, the behavior of the imperfections induces
both amplitude and phase perturbations in the RIS elements and have a non-
static nature with respect to channel coherence time but they present
stationary characteristics within each frame. Figure 3(b) illustrates the
considered STI model. In contrast to the received signal model formulated in
(12), this assumption implies that such imperfection model takes into account
rapid amplitude and phase fluctuations at the RIS elements during the channel
coherence time. Note that the STI induces a block-fading effect in the signals
reflected by the RIS during the $PKT$ symbol periods for channel estimation.
Therefore, by considering an RIS operating under this STI model the resulting
reflection pattern related to the $k$-th time-block at the $p$-th frame can be
written as

$\bar{\mathbf{s}}\left[p,k\right]=\left[e_{n,p}\beta_{1,k}e^{j\phi_{1,k}},\ldots,e_{N,p}\beta_{N,k}e^{j\phi_{1,k}}\right]^{\text{T}},$
(13)

or, equivalently,

$\bar{\mathbf{s}}\left[p,k\right]=\mathbf{e}\left[p\right]\odot\mathbf{s}\left[k\right]\mathbb{C}^{N\times
1},$ (14)

where
$\mathbf{e}\left[p\right]=\left[e_{1,p},\ldots,e_{N,p}\right]^{\text{T}}\in\mathbb{C}^{N\times
1}$ $\forall p=1,\ldots,P$, depends on the $p$-th received frame and models
the unknown non-static amplitude and phase fluctuations along the training
time. The entries of the random vector $\mathbf{e}\left[p\right]$ associated
with the $p$-th frame are defined as

$e_{n,p}=\left\\{\begin{array}[]{ll}1,&\quad{\color[rgb]{0,0,0}\text{non-
impaired case}}\\\ \alpha_{n,p}\cdot
e^{j\theta_{n,p}},&\quad{\color[rgb]{0,0,0}\text{otherwise}},\end{array}\right.$
(15)

where $0\leq\alpha_{n,p}\leq 1$ and $0\leq\theta_{n,p}\leq 2\pi$ $\forall
n=1,\ldots,N$, and $\forall p=1,\ldots,P$, denote the unwanted amplitude
attenuation and phase shift perturbations that affect the $n$-th RIS element
at the $p$-th frame, respectively. It is important to note that the model in
(15) captures different kinds of real-world imperfections at the RIS. For
example, we can note the follow situations:

1) $\alpha_{n,p}\neq 0$ and $\theta_{n,p}\neq 0$ represent the amplitude
absortion and phase shift caused by an object suspended on the $n$-th RIS
element [36], or caused by hardware impairments in the electronic circuits
that make up the RIS [40].

2) $\alpha_{n,p}=1$ and $\theta_{n,p}\neq 0$ represent the phase noise
perturbations from low-resolution phase shifts or phase errors from imperfect
channel estimation [32]-[35].

3) $\alpha_{n,p}=0$ represents the maximum absorption i.e., the $n$-th RIS
element is completely blocked [36].

4) $\alpha_{n,p}=1$ and $\theta_{n,p}=0$ represents the non-impaired RIS in
which no imperfection affects its $n$-th element. Note that in this ideal case
$\bar{\mathbf{s}}\left[p,k\right]=\mathbf{s}\left[p,k\right]$ holds since
$e_{n,p}=1$ $\forall n=1,\ldots,N$, and $\forall p=1,\ldots,P$.

We also observe that the LTI model presented in Section II-A is a particular
case of the generalized STI model when $\alpha_{n,p}=1$ for $P=1$ and
$n=1,\ldots,N$. However, we study these two cases separately for a clearer
description of the static and non-static imperfection scenarios.

By considering an RIS operating under the STI model, we can rewrite the
received signal in (11) as

$\mathbf{Y}\left[p,k\right]=\mathbf{G}\text{diag}\left(\underbrace{\mathbf{e}\left[p\right]\odot\mathbf{s}\left[k\right]}_{\bar{\mathbf{s}}\left[p,k\right]}\right)\mathbf{H}^{\text{T}}+\mathbf{V}\left[p,k\right],$
(16)

$\forall p=1,\ldots,P$, and $\forall k=1,\ldots,K$. In a more convenient form
for our formulation, the received signal (16) can be written in its complete
matrix and decoupled format as

$\mathbf{Y}\left[p,k\right]=\mathbf{G}\mathbf{D}_{p}\left(\mathbf{E}\right)\mathbf{D}_{k}\left(\mathbf{S}\right)\mathbf{H}^{\text{T}}+\mathbf{V}\left[p,k\right],$
(17)

where each row of the matrix
$\mathbf{E}=\left[\mathbf{e}\left[1\right],\ldots,\mathbf{e}\left[P\right]\right]^{\text{T}}\in\mathbb{C}^{P\times
N}$ collects the amplitude and phase parameters for the RIS elements impaired
at the $p$-th frame. Throughout this work, for the two approaches formulated
in Sections II-A and II-B, we assume that a number of $N_{B}=NR_{B}$ random
elements at the RIS are subject to imperfections, where $R_{B}\in[0,1]$
denotes its occurrence probability.

###  II-C Problem Description

In (9) and (14), the RIS reflection pattern corrupted with errors
$\bar{\mathbf{s}}\left[k\right]$ and $\bar{\mathbf{s}}\left[p,k\right]$ induce
unwanted amplitude and/or phase shift responses in the reflected signal by the
RIS creating a mismatch between the ideal reflection pattern and the one that
is actually applied by the RIS. To deal with these imperfections, in this work
we propose to jointly estimate the involved channels $\mathbf{G}$ and
$\mathbf{H}$, as well as the unknown imperfections that affect the $N_{B}$
impaired elements at the RIS. Decoupled estimations of the channels are
required, for instance, to optimize the phase shifts at the RIS, the transmit
precoder at the transmitter and the receive combiner at the receiver in order
to maximize the rate and energy efficiencies in the data transmission phase
[16, 41]. However, the issues of how to utilize the estimated channels to
jointly optimize the RIS phase shifts, the transmitter and receiver active
beamformers, as well as the control overhead problem of optimal phase shifts
[42, 43] are out of the scope of this work and will be addressed in a future
work.

In the following, we show that the received signal models in (12) and (17) can
be represented as third- and fourth-order tensor models, respectively. Then,
we show how such higher-order representations serve as reference models for
the development of efficient PARAFAC-based algorithms to solve the channel
estimation problem for more realistic scenarios where the RIS operates under
LTI and STI.

##  III Proposed Channel Estimation Algorithms

In this section, we propose three different channel estimation algorithms for
the cases with LTI and STI.

###  III-A TALS Algorithm for RIS Operating under LTI (TALS-LTI)

We initially consider the scenario with $N_{B}$ unknown passive elements of
the RIS being affected by LTI as presented in Section II-A. The values of the
phase perturbations as well as their positions are assumed to be unknown at
the receiver.

For simplicity of presentation and without loss of generality, we neglect the
noise term in our formulations. The noiseless part of the received signal in
(12) can be naturally identified as the $k$-th frontal slice of a third-order
tensor $\mathcal{Y}\in\mathbb{C}^{L\times M\times K}$ that admits the
following PARAFAC decomposition [38]

$\mathcal{Y}=\mathcal{I}_{3,N}\times_{1}\mathbf{G}\times_{2}\mathbf{H}\times_{3}\bar{\mathbf{S}}.$
(18)

Here, $\mathcal{I}_{3,N}$ represents a third-order identity tensor of size
$N\times N\times N$. According to (18), the factor matrices related to 1-mode,
2-mode and 3-mode of $\mathcal{Y}$ are respectively
$\mathbf{G}\in\mathbb{C}^{L\times N}$, $\mathbf{H}\in\mathbb{C}^{M\times N}$
and
$\bar{\mathbf{S}}=\left[\bar{\mathbf{s}}\left[1\right],\ldots,\bar{\mathbf{s}}\left[K\right]\right]^{\text{T}}\in\mathbb{C}^{K\times
N}$ where according to (9),
$\bar{\mathbf{s}}\left[k\right]=\mathbf{e}\odot\mathbf{s}\left[k\right]$ or,
equivalently, in matrix notation
$\mathbf{D}_{k}\left(\bar{\mathbf{S}}\right)=\text{diag}\left(\mathbf{e}\right)\mathbf{D}_{k}\left(\mathbf{S}\right)$,
$\forall k=1,\ldots,K$.

Resorting to the multilinear structure of the PARAFAC decomposition in (18),
higher degrees of freedom for signal processing can be achieved by exploiting
the dimensions of interest of the received signal tensor. In this sense,
$\mathcal{Y}$ can also be expressed with respect to its 1-mode and 2-mode
unfoldings, which can be expressed as

$\displaystyle\left[\mathcal{Y}\right]_{(1)}$ $\displaystyle=$
$\displaystyle\mathbf{G}\text{diag}\left(\mathbf{e}\right)\left(\mathbf{S}\diamond\mathbf{H}\right)^{\text{T}}\in\mathbb{C}^{L\times
MK},$ (19) $\displaystyle\left[\mathcal{Y}\right]_{(2)}$ $\displaystyle=$
$\displaystyle\mathbf{H}\text{diag}\left(\mathbf{e}\right)\left(\mathbf{S}\diamond\mathbf{G}\right)^{\text{T}}\in\mathbb{C}^{M\times
LK},$ (20)

where $\left[\mathcal{Y}\right]_{(1)}\in\mathbb{C}^{L\times MK}$ and
$\left[\mathcal{Y}\right]_{(2)}\in\mathbb{C}^{M\times LK}$ are obtained from
(12) by stacking the frontal slices into wide matrices, as defined by

$\displaystyle\left[\mathcal{Y}\right]_{(1)}$ $\displaystyle=$
$\displaystyle\left[\mathbf{Y}\left[1\right],\ldots,\mathbf{Y}\left[K\right]\right],$
(21) $\displaystyle\left[\mathcal{Y}\right]_{(2)}$ $\displaystyle=$
$\displaystyle\left[\mathbf{Y}^{\text{T}}\left[1\right],\ldots,\mathbf{Y}^{\text{T}}\left[K\right]\right].$
(22)

Additionally, (19) can also be represented in a covenient vectorized form. By
applying the property (3) to (19) we obtain

$\text{vec}\left(\left[\mathcal{Y}\right]_{(1)}\right)=\left(\mathbf{S}\diamond\mathbf{H}\diamond\mathbf{G}\right)\mathbf{e}\in\mathbb{C}^{LMK\times
1}.$ (23)

In the following, we describe an iterative way to estimate the channel
matrices $\mathbf{H}$ and $\mathbf{G}$ from the received signal tensor
$\mathcal{Y}$ that models the RIS-assisted MIMO system with RIS operating
under LTI. The estimation problem can be solved by computing a rank-$N$
approximation to the PARAFAC decomposition of $\mathcal{Y}$, i.e,

$\underset{\mathbf{G},\mathbf{H},\mathbf{e}}{\text{min}}\sum_{k=1}^{K}\left\|\mathbf{Y}\left[k\right]-\mathbf{G}\text{diag}\left(\mathbf{e}\right)\mathbf{D}_{k}\left(\mathbf{S}\right)\mathbf{H}^{\text{T}}\right\|_{\text{F}}^{2}.$
(24)

Problem (24) can be solved in an efficient form by an ALS algorithm [38, 44].
It is a well-known iterative method for estimating the factor matrices of a
tensor model thanks to its implementation simplicity and monotonic convergence
property in which the update of every given matrix at each iteration may
either improve or maintain but cannot worsen the current fit, leading usually
to global minimum solution [45, 46]. The decoupled estimates of $\mathbf{G}$,
$\mathbf{H}$ and $\mathbf{e}$ can be obtained by converting the trilinear
fitting problem in (24) into the following three simplest linear least squares
(LS) sub-problems formulated from (19), (20) and (23), respectively

$\displaystyle\hat{\mathbf{G}}$ $\displaystyle=$
$\displaystyle\underset{{\color[rgb]{0,0,0}\mathbf{G}\mid\mathbf{e},\mathbf{H}}}{\text{argmin}}\left\|\left[\mathcal{Y}\right]_{(1)}-\mathbf{G}\text{diag}\left(\mathbf{e}\right)\left(\mathbf{S}\diamond\mathbf{H}\right)^{\text{T}}\right\|_{\text{F}}^{2},$
(25) $\displaystyle\hat{\mathbf{H}}$ $\displaystyle=$
$\displaystyle\underset{{\color[rgb]{0,0,0}\mathbf{H}\mid\mathbf{e},\mathbf{G}}}{\text{argmin}}\left\|\left[\mathcal{Y}\right]_{(2)}-\mathbf{H}\text{diag}\left(\mathbf{e}\right)\left(\mathbf{S}\diamond\mathbf{G}\right)^{\text{T}}\right\|_{\text{F}}^{2},$
(26) $\displaystyle\hat{\mathbf{e}}$ $\displaystyle=$
$\displaystyle\underset{{\color[rgb]{0,0,0}\mathbf{e}\mid\textbf{H},\mathbf{G}}}{\text{argmin}}\left\|\text{vec}\left(\left[\mathcal{Y}\right]_{(1)}\right)-\left(\mathbf{S}\diamond\mathbf{H}\diamond\mathbf{G}\right)\mathbf{e}\right\|_{\text{F}}^{2}.$
(27)

According to (25), the conditional LS update for $\hat{\mathbf{G}}$ is given
by

$\hat{\mathbf{G}}=\left[\mathcal{Y}\right]_{(1)}\left[\text{diag}\left(\mathbf{e}\right)\left(\mathbf{S}\diamond\mathbf{H}\right)^{\text{T}}\right]^{{\dagger}}.$
(28)

Similarly, according to (26) and (27), the conditional LS updates for
$\hat{\mathbf{H}}$ and $\hat{\mathbf{e}}$ are respectively given by

$\displaystyle\hat{\mathbf{H}}$ $\displaystyle=$
$\displaystyle\left[\mathcal{Y}\right]_{(2)}\left[\text{diag}\left(\mathbf{e}\right)\left(\mathbf{S}\diamond\mathbf{G}\right)^{\text{T}}\right]^{{\dagger}},$
(29) $\displaystyle\hat{\mathbf{e}}$ $\displaystyle=$
$\displaystyle\left(\mathbf{S}\diamond\mathbf{H}\diamond\mathbf{G}\right)^{{\dagger}}\text{vec}\left(\left[\mathcal{Y}\right]_{(1)}\right).$
(30)

1. Set $i=0$;

1 Keep $\mathbf{S}$ fixed; Initialize randomly the matrix
$\hat{\mathbf{H}}_{(i=0)}$

and the imperfections vector $\hat{\mathbf{e}}_{(i=0)}$;

2. $i\leftarrow i+1$;

3. According to (28), obtain an LS estimate of $\hat{\mathbf{G}}_{(i)}$:

$\hat{\mathbf{G}}_{(i)}=\left[\mathcal{Y}\right]_{(1)}\left[\text{diag}\left(\hat{\mathbf{e}}_{(i-1)}\right)\left(\mathbf{S}\diamond\hat{\mathbf{H}}_{(i-1)}\right)^{\text{T}}\right]^{{\dagger}};$

4. According to (29), obtain an LS estimate of $\hat{\mathbf{H}}_{(i)}$:

$\hat{\mathbf{H}}_{(i)}=\left[\mathcal{Y}\right]_{(2)}\left[\text{diag}\left(\hat{\mathbf{e}}_{(i-1)}\right)\left(\mathbf{S}\diamond\hat{\mathbf{G}}_{(i)}\right)^{\text{T}}\right]^{{\dagger}};$

5. According to (30), obtain an LS estimate of $\hat{\mathbf{e}}_{(i)}$:

$\hat{\mathbf{e}}_{(i)}=\left(\mathbf{S}\diamond\hat{\mathbf{H}}_{(i)}\diamond\hat{\mathbf{G}}_{(i)}\right)^{{\dagger}}\text{vec}\left(\left[\mathcal{Y}\right]_{(1)}\right);$

2 6. Calculate the residual error
${\color[rgb]{0,0,0}\epsilon}_{(i)}=\|\left[\mathcal{Y}\right]_{(1)}-\hat{\left[\mathcal{Y}\right]}_{(1){\color[rgb]{0,0,0}(i)}}\|_{\text{F}}^{2}$

where

$\hat{\left[\mathcal{Y}\right]}_{(1){\color[rgb]{0,0,0}(i)}}=\hat{\mathbf{G}}_{(i)}\text{diag}\left(\hat{\mathbf{e}}_{(i)}\right)\left(\mathbf{S}\diamond\hat{\mathbf{H}}_{(i)}\right)^{\text{T}};$

7. Repeat Steps 2-6 until $|{\color[rgb]{0,0,0}\epsilon}_{(i)}-{\color[rgb]{0,0,0}\epsilon}_{(i-1)}|\leq{\color[rgb]{0,0,0}\delta}$.

Algorithm 1 TALS-LTI Algorithm

The proposed TALS-STI algorithm consists of three iterative and alternating
update steps formulated from the LS solutions in (28), (29) and (30). At each
step, the fitting error is minimized with respect to one given factor matrix
by fixing the other matrices to their values obtained at previous updating
steps. This procedure is repeated until the convergence of the algorithm at
the $i$-th iteration determined by the designer.

Define

$\hat{\left[\mathcal{Y}\right]}_{(1){\color[rgb]{0,0,0}(i)}}=\hat{\mathbf{G}}_{(i)}\text{diag}\left(\hat{\mathbf{e}}_{(i)}\right)\left(\mathbf{S}\diamond\hat{\mathbf{H}}_{(i)}\right)^{\text{T}}\in\mathbb{C}^{L\times
MK},$ (31)

as the reconstructed version of $\left[\mathcal{Y}\right]_{(1)}$ obtained from
the estimates of $\hat{\mathbf{G}}_{(i)}$, $\hat{\mathbf{H}}_{(i)}$, and
$\hat{\mathbf{e}}_{(i)}$, and the residual error as

${\color[rgb]{0,0,0}\epsilon_{(i)}}=\left\|\left[\mathcal{Y}\right]_{(1)}-\hat{\left[\mathcal{Y}\right]}_{(1){\color[rgb]{0,0,0}(i)}}\right\|_{\text{F}}^{2},$
(32)

computed at the end of the $i$-th iteration. The convergence of the algorithm
is declared when
$|{\color[rgb]{0,0,0}\epsilon}_{(i)}-{\color[rgb]{0,0,0}\epsilon}_{(i-1)}|\leq\delta$,
with $\delta$ being a constant considered by the designer, meaning that the
reconstruction error does not significantly change between two successive
iterations. In this work, we set $\delta=10^{-6}$ as a convergence threshold.
The implementation steps of the proposed iterative TALS-LTI algorithm are
summarized in the pseudocode shown in Algorithm 1. For the complexity analysis
of Algorithm 1, see Section IV.

###  III-B TALS Algorithm for RIS Operating under STI (TALS-STI)

In order to derive proposed channel estimators for a scenario with STI, let us
first establish a link between the received signal in (17) and the PARAFAC
decomposition. According to [38], the noiseless signal part of (17) expresses
the ($p,k$)-th frontal slice of a fourth-order tensor
$\mathcal{Y}\in\mathbb{C}^{L\times M\times K\times P}$ that follows the
PARAFAC decomposition

$\mathcal{Y}=\mathcal{I}_{4,N}\times_{1}\mathbf{G}\times_{2}\mathbf{H}\times_{3}\mathbf{S}\times_{4}\mathbf{E}.$
(33)

Here, $\mathcal{I}_{4,N}$ denotes the fourth-order identity tensor of size
$N\times N\times N\times N$, while $\mathbf{G}$, $\mathbf{H}$, $\mathbf{S}$
and $\mathbf{E}$ are the 1,2,3,4-mode factor matrices of the decomposition,
respectively.

By stacking column-wise the noiseless received signal in (17) for the $K$
time-blocks at frame $p$ as the matrix
$\mathbf{Y}_{p}=\left[\mathbf{Y}\left[p,1\right],\ldots,\mathbf{Y}\left[p,K\right]\right]\in\mathbb{C}^{L\times
MK}$, we have

$\mathbf{Y}_{p}=\mathbf{G}\mathbf{D}_{p}\left(\mathbf{E}\right)\left[\mathbf{D}_{1}\left(\mathbf{S}\right)\mathbf{H}^{\text{T}},\ldots,\mathbf{D}_{K}\left(\mathbf{S}\right)\mathbf{H}^{\text{T}}\right],$
(34)

$\forall p=1,\ldots,P$. Applying the property (2) to the right-hand side of
(34), a more compact form is obtained as

$\mathbf{Y}_{p}=\mathbf{G}\mathbf{D}_{p}\left(\mathbf{E}\right)\left(\mathbf{S}\diamond\mathbf{H}\right)^{\text{T}}\in\mathbb{C}^{L\times
MK}.$ (35)

From (35), we can define the new column-wise collection
$\left[\mathcal{Y}\right]_{(1)}=\left[\mathbf{Y}_{1},\ldots,\mathbf{Y}_{P}\right]\in\mathbb{C}^{L\times
MKP}$ as the 1-mode matrix unfolding of the received signal tensor
$\mathcal{Y}\in\mathbb{C}^{L\times M\times K\times P}$ in (33), which is given
by

$\left[\mathcal{Y}\right]_{(1)}=\mathbf{G}\left[\mathbf{D}_{1}\left(\mathbf{E}\right)\left(\mathbf{S}\diamond\mathbf{H}\right)^{\text{T}},\ldots,\mathbf{D}_{P}\left(\mathbf{E}\right)\left(\mathbf{S}\diamond\mathbf{H}\right)^{\text{T}}\right].$
(36)

By applying property (2) to the right-hand side of (36), we finally obtain

$\left[\mathcal{Y}\right]_{(1)}=\mathbf{G}\left(\mathbf{E}\diamond\mathbf{S}\diamond\mathbf{H}\right)^{\text{T}}\in\mathbb{C}^{L\times
MKP}.$ (37)

Additionally, for our purpose, we also need to define the 2-mode, 3-mode and
4-mode matrix unfoldings of the fourth-order received signal tensor
$\mathcal{Y}\in\mathbb{C}^{L\times M\times K\times P}$ since they will be
exploited to formulate our second set of channel estimation algorithms in the
sequel. The remaining unfoldings can be deduced using a similar procedure by
permuting the factor matrices in (17). This leads to the following
factorizations to the other unfoldings

$\displaystyle\left[\mathcal{Y}\right]_{(2)}$ $\displaystyle=$
$\displaystyle\mathbf{H}\left(\mathbf{E}\diamond\mathbf{S}\diamond\mathbf{G}\right)^{\text{T}}\in\mathbb{C}^{M\times
LKP},$ (38) $\displaystyle\left[\mathcal{Y}\right]_{(3)}$ $\displaystyle=$
$\displaystyle\mathbf{S}\left(\mathbf{E}\diamond\mathbf{H}\diamond\mathbf{G}\right)^{\text{T}}\in\mathbb{C}^{K\times
LMP},$ (39) $\displaystyle\left[\mathcal{Y}\right]_{(4)}$ $\displaystyle=$
$\displaystyle\mathbf{E}\left(\mathbf{S}\diamond\mathbf{H}\diamond\mathbf{G}\right)^{\text{T}}\in\mathbb{C}^{P\times
LMK}.$ (40)

In the following, we show that the channel matrices $\mathbf{H}$ and
$\mathbf{G}$ can be also estimated when STI are assumed. From the received
signal tensor described in (33), the estimates of $\mathbf{G}$ and
$\mathbf{H}$ can be obtained by minimizing the following quadrilinear LS
fitting problem

$\underset{\mathbf{G},\mathbf{H},\mathbf{E}}{\text{min}}\sum_{p=1}^{P}\sum_{k=1}^{K}\left\|\mathbf{Y}\left[p,k\right]-\mathbf{G}\mathbf{D}_{p}\left(\mathbf{E}\right)\mathbf{D}_{k}\left(\mathbf{S}\right)\mathbf{H}^{\text{T}}\right\|_{\text{F}}^{2}.$
(41)

Similar to the TALS-LTI algorithm, we also propose to solve this optimization
problem by means of the ALS algorithm. Since the matrix $\mathbf{S}$ is known
at the receiver, the quadrilinear fitting problem in (41) is simplified to a
trilinear fitting problem that reduces to iteratively minimize the following
linear LS sub-problems formulated from (37), (38) and (40), respectively

$\displaystyle\hat{\mathbf{G}}$ $\displaystyle=$
$\displaystyle\underset{{\color[rgb]{0,0,0}\mathbf{G}\mid\mathbf{E},\mathbf{H}}}{\text{argmin}}\left\|\left[\mathcal{Y}\right]_{(1)}-\mathbf{G}\left(\mathbf{E}\diamond\mathbf{S}\diamond\mathbf{H}\right)^{\text{T}}\right\|_{\text{F}}^{2},$
(42) $\displaystyle\hat{\mathbf{H}}$ $\displaystyle=$
$\displaystyle\underset{{\color[rgb]{0,0,0}\mathbf{H}\mid\mathbf{E},\mathbf{G}}}{\text{argmin}}\left\|\left[\mathcal{Y}\right]_{(2)}-\mathbf{H}\left(\mathbf{E}\diamond\mathbf{S}\diamond\mathbf{G}\right)^{\text{T}}\right\|_{\text{F}}^{2},$
(43) $\displaystyle\hat{\mathbf{E}}$ $\displaystyle=$
$\displaystyle\underset{{\color[rgb]{0,0,0}\mathbf{E}\mid\mathbf{H},\mathbf{G}}}{\text{argmin}}\left\|\left[\mathcal{Y}\right]_{(4)}-\mathbf{E}\left(\mathbf{S}\diamond\mathbf{H}\diamond\mathbf{G}\right)^{\text{T}}\right\|_{\text{F}}^{2}.$
(44)

It follows from (42), (43) and (44) that the conditional LS updates of
$\hat{\mathbf{G}}$, $\hat{\mathbf{H}}$ and $\hat{\mathbf{E}}$ are given by

$\displaystyle\hat{\mathbf{G}}$ $\displaystyle=$
$\displaystyle\left[\mathcal{Y}\right]_{(1)}\left[\left(\mathbf{E}\diamond\mathbf{S}\diamond\mathbf{H}\right)^{\text{T}}\right]^{{\dagger}},$
(45) $\displaystyle\hat{\mathbf{H}}$ $\displaystyle=$
$\displaystyle\left[\mathcal{Y}\right]_{(2)}\left[\left(\mathbf{E}\diamond\mathbf{S}\diamond\mathbf{G}\right)^{\text{T}}\right]^{{\dagger}},$
(46) $\displaystyle\hat{\mathbf{E}}$ $\displaystyle=$
$\displaystyle\left[\mathcal{Y}\right]_{(4)}\left[\left(\mathbf{S}\diamond\mathbf{H}\diamond\mathbf{G}\right)^{\text{T}}\right]^{{\dagger}},$
(47)

respectively. In the same way as in Algorithm 1, the updates of
$\hat{\mathbf{G}}$, $\hat{\mathbf{H}}$ and $\hat{\mathbf{E}}$ are obtained by
iteratively performing (45), (46) and (47) until the convergence. The proposed
iterative TALS-STI algorithm is detailed in the pseudocode shown in Algorithm
2. For the complexity analysis of Algorithm 2, see Section IV.

1. Set $i=0$;

Keep $\mathbf{S}$ fixed;

1 Initialize randomly the matrices $\hat{\mathbf{H}}_{(i=0)}$ and
$\hat{\mathbf{E}}_{(i=0)}$;

2. $i\leftarrow i+1$;

3. According to (45), obtain an LS estimate of $\hat{\mathbf{G}}_{(i)}$:

$\hat{\mathbf{G}}_{(i)}=\left[\mathcal{Y}\right]_{(1)}\left[\left(\hat{\mathbf{E}}_{(i-1)}\diamond\mathbf{S}\diamond\hat{\mathbf{H}}_{(i-1)}\right)^{\text{T}}\right]^{{\dagger}};$

4. According to (46), obtain an LS estimate of $\hat{\mathbf{H}}_{(i)}$:

$\hat{\mathbf{H}}_{(i)}=\left[\mathcal{Y}\right]_{(2)}\left[\left(\hat{\mathbf{E}}_{(i-1)}\diamond\mathbf{S}\diamond\hat{\mathbf{G}}_{(i)}\right)^{\text{T}}\right]^{{\dagger}};$

5. According to (47), obtain an LS estimate of $\hat{\mathbf{E}}_{(i)}$:

$\hat{\mathbf{E}}_{(i)}=\left[\mathcal{Y}\right]_{(4)}\left[\left(\mathbf{S}\diamond\hat{\mathbf{H}}_{(i)}\diamond\hat{\mathbf{G}}_{(i)}\right)^{\text{T}}\right]^{{\dagger}};$

2 6. Calculate the residual error
${\color[rgb]{0,0,0}\epsilon}_{(i)}=\|\left[\mathcal{Y}\right]_{(1)}-\hat{\left[\mathcal{Y}\right]}_{(1){\color[rgb]{0,0,0}(i)}}\|_{\text{F}}^{2}$

where

$\hat{\left[\mathcal{Y}\right]}_{(1){\color[rgb]{0,0,0}(i)}}=\hat{\mathbf{G}}_{(i)}\left(\hat{\mathbf{E}}_{(i)}\diamond\mathbf{S}\diamond\hat{\mathbf{H}}_{(i)}\right)^{\text{T}};$

7. Repeat Steps 2-6 until $|{\color[rgb]{0,0,0}\epsilon}_{(i)}-{\color[rgb]{0,0,0}\epsilon}_{(i-1)}|\leq{\color[rgb]{0,0,0}\delta}$.

Algorithm 2 TALS-STI Algorithm

Remark 1: In the application context of this work, the non-impaired RIS
activation pattern matrix $\mathbf{S}$ is assumed to be known at the receiver
as indicated in the first step of Algorithms 1-2, respectively. This is a
feasible assumption in accordance with the channel estimation protocol shown
in Fig. 2. Among different design possibilities, we set $\mathbf{S}$ as a
semi-unitary matrix satisfying
$\mathbf{S}^{\text{H}}\mathbf{S}=K\mathbf{I}_{N}$. According to [29], a good
choice is to consider $\mathbf{S}$ as a deterministic truncated discrete
Fourier transform (DFT) matrix. This choice guarantees a good performance of
the proposed algorithms since the correlation properties of the additive noise
are not affected during the estimation processing. More details on the optimal
design of $\mathbf{S}$ are found in [47]. Moreover, despite the iterative
nature of the proposed TALS-LTI and TALS-STI algorithms, the convergence to
the global minimum is always achieved within a few iterations (usually less
than $80$ iterations as verified in our simulation results) due to the
knowledge of $\mathbf{S}$ that remains fixed during the iterations.

###  III-C HOSVD Algorithm for RIS Operating under STI (HOSVD-STI)

We now derive a closed-form solution based on higher order singular value
decomposition (HOSVD) for channel estimation under the STI model. According to
(39), the transpose of the 3-mode unfolding of $\mathcal{Y}$ is denoted by

$\left[\mathcal{Y}\right]^{\text{T}}_{(3)}=\left(\mathbf{E}\diamond\mathbf{H}\diamond\mathbf{G}\right)\mathbf{S}^{\text{T}}.$
(48)

The first processing step at the receiver is to apply a bilinear time-domain
matched-filtering by multiplying both sides in (48) by the pseudo-inverse of
$\mathbf{S}^{\text{T}}$, resulting in

$\tilde{\mathbf{Y}}=\mathbf{E}\diamond\mathbf{H}\diamond\mathbf{G}\in\mathbb{C}^{LMP\times
N},$ (49)

where
$\tilde{\mathbf{Y}}=\left[\mathcal{Y}\right]^{\text{T}}_{(3)}\left(\mathbf{S}^{\text{T}}\right)^{{\dagger}}$.
From (49), decoupled estimates of the channel matrices and RIS imperfections
can be obtained by separating each factor matrix in the Khatri-Rao product. In
this sense, the estimates can be obtained by minimizing the following cost
function

$\underset{\mathbf{G},\mathbf{H},\mathbf{E}}{\text{min}}\left\|\tilde{\mathbf{Y}}-\mathbf{E}\diamond\mathbf{H}\diamond\mathbf{G}\right\|_{\text{F}}^{2}.$
(50)

Here, we propose to solve this problem by means of multiple rank-one tensor
approximations via the HOSVD [39]. To this end, let us define
$\tilde{\mathbf{Y}}=\left[\tilde{\mathbf{y}}_{1},\ldots,\tilde{\mathbf{y}}_{N}\right]\in\mathbb{C}^{LMP\times
N}$, and note that according to (1), the $n$-th column of $\tilde{\mathbf{Y}}$
can be written as

$\tilde{\mathbf{y}}_{n}=\mathbf{e}_{n}\otimes\mathbf{h}_{n}\otimes\mathbf{g}_{n}\in\mathbb{C}^{LMP\times
1},$ (51)

where $\mathbf{e}_{n}\in\mathbb{C}^{P\times 1}$,
$\mathbf{h}_{n}\in\mathbb{C}^{M\times 1}$ and
$\mathbf{g}_{n}\in\mathbb{C}^{L\times 1}$ denote the $n$-th column of
$\mathbf{E}$, $\mathbf{H}$ and $\mathbf{G}$, respectively. Using the
equivalence property in (4) that relates the Kronecker product to the outer
product, we can rewrite (51) as

$\tilde{\mathbf{y}}_{n}=\text{vec}\left(\mathbf{g}_{n}\circ\mathbf{h}_{n}\circ\mathbf{e}_{n}\right)\in\mathbb{C}^{LMP\times
1},$ (52)

that represents the vectorized form of the following third-order rank-one
tensor

$\tilde{\mathcal{Y}}_{n}=\mathbf{g}_{n}\circ\mathbf{h}_{n}\circ\mathbf{e}_{n}\in\mathbb{C}^{L\times
M\times P}.$ (53)

Thus, the optimization problem in (50) is equivalent to finding the estimates
of $\mathbf{H}$, $\mathbf{G}$ and $\mathbf{E}$ that minimize a set of $N$
rank-one tensor approximations, i.e,

$\left(\hat{\mathbf{G}},\hat{\mathbf{H}},\hat{\mathbf{E}}\right)=\underset{\mathbf{G},\mathbf{H},\mathbf{E}}{\text{argmin}}\sum_{n=1}^{N}\left\|\tilde{\mathcal{Y}}_{n}-\mathbf{g}_{n}\circ\mathbf{h}_{n}\circ\mathbf{e}_{n}\right\|_{\text{F}}^{2}.$
(54)

for $n=1,\ldots,N$

1 1. Rearrange the $n$-th column of $\tilde{\mathbf{Y}}$ in Equation (49)

as the rank-one tensor $\tilde{\mathcal{Y}}_{n}$ in Equation (53);

2. HOSVD procedure

2.1 Compute $\mathbf{U}_{n}^{(1)}$ as the $L$ left singular vectors of
$\left[\tilde{\mathcal{Y}}_{n}\right]_{(1)}$:
$\left[\tilde{\mathcal{Y}}_{n}\right]_{(1)}=\mathbf{U}_{n}^{(1)}\cdot\mathbf{\Sigma}_{n}^{(1)}\cdot\mathbf{V}_{n}^{(1)\text{H}};$

2.2 Compute $\mathbf{U}_{n}^{(2)}$ as the $M$ left singular vectors of
$\left[\tilde{\mathcal{Y}}_{n}\right]_{(2)}$:
$\left[\tilde{\mathcal{Y}}_{n}\right]_{(2)}=\mathbf{U}_{n}^{(2)}\cdot\mathbf{\Sigma}_{n}^{(2)}\cdot\mathbf{V}_{n}^{(2)\text{H}};$

2.3 Compute $\mathbf{U}_{n}^{(3)}$ as the $P$ left singular vectors of
$\left[\tilde{\mathcal{Y}}_{n}\right]_{(3)}$:
$\left[\tilde{\mathcal{Y}}_{n}\right]_{(3)}=\mathbf{U}_{n}^{(3)}\cdot\mathbf{\Sigma}_{n}^{(3)}\cdot\mathbf{V}_{n}^{(3)\text{H}};$

2.4 Compute the HOSVD core tensor $\mathcal{G}_{n}$ as:
$\mathcal{G}_{n}=\tilde{\mathcal{Y}}_{n}\times_{1}\mathbf{U}_{n}^{(1)\text{H}}\times_{2}\mathbf{U}_{n}^{(2)\text{H}}\times_{3}\mathbf{U}_{n}^{(3)\text{H}};$

end procedure

2 3. Obtain the estimates for $\hat{\mathbf{g}}_{n}$, $\hat{\mathbf{h}}_{n}$
and $\hat{\mathbf{e}}_{n}$ from Equations

(56), (57) and (58), respectively;

end

3 4. Return the matrices
$\hat{\mathbf{G}}=\left[\hat{\mathbf{g}}_{1},\ldots,\hat{\mathbf{g}}_{N}\right]$,
$\hat{\mathbf{H}}=\left[\hat{\mathbf{h}}_{1},\ldots,\hat{\mathbf{h}}_{N}\right]$

and
$\hat{\mathbf{E}}=\left[\hat{\mathbf{e}}_{1},\ldots,\hat{\mathbf{e}}_{N}\right]$.

Algorithm 3 HOSVD-STI Algorithm

Let us introduce the HOSVD of $\tilde{\mathcal{Y}}_{n}$ as

$\tilde{\mathcal{Y}}_{n}=\mathcal{G}_{n}\times_{1}\mathbf{U}_{n}^{(1)}\times_{2}\mathbf{U}_{n}^{(2)}\times_{3}\mathbf{U}_{n}^{(3)}\in\mathbb{C}^{L\times
M\times P},$ (55)

where $\mathbf{U}_{n}^{(1)}\in\mathbb{C}^{L\times L}$,
$\mathbf{U}_{n}^{(2)}\in\mathbb{C}^{M\times M}$ and
$\mathbf{U}_{n}^{(3)}\in\mathbb{C}^{P\times P}$ are unitary matrices, while
$\mathcal{G}_{n}\in\mathbb{C}^{L\times M\times P}$ denotes the HOSVD core
tensor. The estimates of the vectors $\mathbf{g}_{n}$, $\mathbf{h}_{n}$ and
$\mathbf{e}_{n}$ that solve the LS problem in (54) can be obtained by
truncating the HOSVD of $\tilde{\mathcal{Y}}_{n}$ to its dominant rank-one
component, yielding

$\displaystyle\hat{\mathbf{g}}_{n}$ $\displaystyle=$
$\displaystyle\sqrt[3]{\left(\mathcal{G}_{n}\right)_{1,1,1}}\cdot\mathbf{u}^{(1)}_{1,n},$
(56) $\displaystyle\hat{\mathbf{h}}_{n}$ $\displaystyle=$
$\displaystyle\sqrt[3]{\left(\mathcal{G}_{n}\right)_{1,1,1}}\cdot\mathbf{u}^{(2)}_{1,n},$
(57) $\displaystyle\hat{\mathbf{e}}_{n}$ $\displaystyle=$
$\displaystyle\sqrt[3]{\left(\mathcal{G}_{n}\right)_{1,1,1}}\cdot\mathbf{u}^{(3)}_{1,n},$
(58)

where $\mathbf{u}^{(1)}_{1,n}\in\mathbb{C}^{L\times 1}$,
$\mathbf{u}^{(2)}_{1,n}\in\mathbb{C}^{M\times 1}$ and
$\mathbf{u}^{(3)}_{1,n}\in\mathbb{C}^{P\times 1}$ are the first higher order
singular vectors, i.e., the first column of $\mathbf{U}_{n}^{(1)}$,
$\mathbf{U}_{n}^{(2)}$ and $\mathbf{U}_{n}^{(3)}$, respectively. Here,
$\left(\mathcal{G}_{n}\right)_{1,1,1}$ is the first element of the core tensor
$\mathcal{G}_{n}$. The estimates of $\hat{\mathbf{G}}$, $\hat{\mathbf{H}}$ and
$\hat{\mathbf{E}}$ are obtained by repeating the procedure of (53)-(58) for
the $N$ columns of $\tilde{\mathbf{Y}}$ in (49). In other words, a total of
$N$ rank-one tensor approximations via HOSVD are necessary to obtain the full
estimates of the matrices
$\hat{\mathbf{G}}=\left[\hat{\mathbf{g}}_{1},\ldots,\hat{\mathbf{g}}_{N}\right]$,
$\hat{\mathbf{H}}=\left[\hat{\mathbf{h}}_{1},\ldots,\hat{\mathbf{h}}_{N}\right]$
and
$\hat{\mathbf{E}}=\left[\hat{\mathbf{e}}_{1},\ldots,\hat{\mathbf{e}}_{N}\right]$
in a closed-form manner. The implementation steps of the proposed closed-form
HOSVD-STI algorithm are summarized in Algorithm 3. For the complexity analysis
of Algorithm 3, see Section IV.

##  IV Identifiability and Computational Complexity

In this section, we examine the identifiability aspects and the computational
complexity associated with the proposed tensor-based TALS-LTI, TALS-STI and
HOSVD-STI algorithms.

1) TALS-LTI algorithm: According to (28) and (29), the uniqueness of the LS
estimates of $\hat{\mathbf{G}}$ and $\hat{\mathbf{H}}$ requires that
$\text{diag}\left(\mathbf{e}\right)\left(\mathbf{S}\diamond\mathbf{H}\right)^{\text{T}}\in\mathbb{C}^{N\times
KM}$ and
$\text{diag}\left(\mathbf{e}\right)\left(\mathbf{S}\diamond\mathbf{G}\right)^{\text{T}}\in\mathbb{C}^{N\times
KL}$ are full row-rank to be right-invertible. Additionally, the uniqueness of
the LS estimate of $\hat{\mathbf{e}}$ requires that
$\left(\mathbf{S}\diamond\mathbf{H}\diamond\mathbf{G}\right)\in\mathbb{C}^{KLM\times
N}$ be full column-rank to be left-invertible in (30). This means that the
conditions $N\leq KM$, $N\leq KL$ and $N\leq KLM$ must be jointly satisfied.
By combining these three necessary and sufficient conditions, we obtain the
lower bound on the number of time-blocks $K$ necessary for channel estimation
so that Steps 3, 4 and 5 in Algorithm 1 yield unique solutions:

$K\geq\left\lceil\frac{N}{\text{min}\left(L,M\right)}\right\rceil.$ (59)

2) TALS-STI algorithm: We can note from (45), (46) and (47) that unique
estimates of $\hat{\mathbf{G}}$, $\hat{\mathbf{H}}$ and $\hat{\mathbf{E}}$ in
the LS sense requires that
$\left(\mathbf{E}\diamond\mathbf{S}\diamond\mathbf{H}\right)^{\text{T}}\in\mathbb{C}^{N\times
KMP}$,
$\left(\mathbf{E}\diamond\mathbf{S}\diamond\mathbf{G}\right)^{\text{T}}\in\mathbb{C}^{N\times
KLP}$ and
$\left(\mathbf{S}\diamond\mathbf{H}\diamond\mathbf{G}\right)^{\text{T}}\in\mathbb{C}^{N\times
KLM}$ are full row-rank to be right-invertible. This means that the conditions
$N\leq KMP$, $N\leq KLP$ and $N\leq KLM$ must be satisfied. The combination of
these inequalities leads to the following necessary and sufficient condition
to be satisfied

$K\geq\left\lceil\frac{N}{\text{min}\left(MP,LP,LM\right)}\right\rceil.$ (60)

The condition in (60) establishes the lower-bound on the required number of
time-blocks $K$ so that Steps 3, 4 and 5 in Algorithm 2 provide unique
solutions when STI are assumed.

3) HOSVD-STI algorithm: In contrast to the iterative TALS-LTI and TALS-STI
algorithms in which three LS conditions must be jointly satisfied, the
proposed HOSVD-STI algorithm is a closed-form solution requiring only that the
RIS activation pattern matrix $\mathbf{S}\in\mathbb{C}^{K\times N}$ has full
column-rank in order to guarantee the uniqueness in the LS sense when the
bilinear time-domain matched-filtering preprocessing is performed at the
receiver side, as indicated in (49). This leads to the following necessary and
sufficient condition to the use of the HOSVD-STI algorithm

$K\geq N.$ (61)

This condition is satisfied under the truncated DFT design discussed in Remark
1, where the RIS activation pattern matrix has orthonormal columns.

According to (59)-(61), independent of the imperfection model, the required
number of time-blocks $K$ scales with the number of RIS elements $N$ at least
linearly. However, the proposed TALS-LTI and TALS-STI algorithms present more
flexible operation conditions for channel estimation compared to our HOSVD-STI
algorithm. On the other hand, the HOSVD-STI is a closed-form approach that
presents lower computational complexity compared to the others one, as
discussed in the following.

4) Computational complexity: As can be observed in Algorithms 1-2, the
computational complexity of the proposed TALS-LTI and TALS-STI algorithms is
dominated by the cost associated with the computation of the matrix pseudo-
inverses in three LS update steps that calculate the estimates of the channels
and imperfection matrices in an iterative and alternating way. Therefore, the
computational complexity of the TALS-LTI and TALS-STI algorithms are
$\mathcal{O}\left(N^{2}K[M+L+ML]\right)$ and
$\mathcal{O}\left(N^{2}K[PM+PL+ML]\right)$ per iteration, respectively.
Regarding the proposed HOSVD-STI in Algorithm 3, its computational complexity
is dominated by the HOSVD of the third-order rank-one tensor in (53), which is
equivalent to compute the truncated SVDs of its 1-mode, 2-mode and 3-mode
unfolding matrices to rank-one. These truncated SVDs are repeated $N$ times.
Therefore, the HOSVD-STI algorithm has a complexity
$\mathcal{O}\left(NMLP\right)$, which is clearly less than that of the TALS-
LTI and TALS-STI algorithms. Specially, as opposed to Algorithms 1-2, where
the computational complexity is proportional to $N^{2},$ in Algorithm 3 the
complexity scales with the RIS size $N$ linearly. As a result, compared to
Algorithms 1-2, the relative gain of Algorithm 3, in terms of complexity
increases rapidly as the number of RIS elements increases.

Remark 2: By comparing the conditions in (60) and (61), we observe that the
TALS-STI algorithm presents a less restrictive requirement on the minimum
number $K$ of time-blocks necessary for the channel estimation compared to
HOSVD-STI algorithm that works well when $K\geq N$ is satisfied. On the other
hand, in the HOSVD-STI algorithm a single HOSVD is computed independently for
each column of the post-filtered signal in (49). This leads to a total of $N$
independent rank-one tensor approximations via HOSVD that can be parallelized
if more than one processor is available leading to a considerable reduction in
the total processing time for channel estimation. Based on this, we can
conclude that the TALS-STI algorithm may be preferred in situations where more
flexible choices for the number of time-blocks, $K$, are required while the
HOSVD-STI algorithm becomes more attractive when low processing delay is
needed and the receiver is equipped with multiple processors. Hence, there is
a tradeoff between computational complexity and operation conditions for the
proposed tensor-based algorithms.

Remark 3: Since $\mathbf{S}$ is known at the receiver, by satisfying the
design recommendations given by (59), (60) and (61) the estimated factor
matrices $\hat{\mathbf{G}}$, $\hat{\mathbf{H}}$ and $\hat{\mathbf{E}}$ do not
suffer from column permutation ambiguity. The scaling ambiguity affecting the
columns of the estimated matrices can be eliminated with a simple
normalization procedure, as performed in [22] and [29].

(a)

(b)

(c)

Figure 4: NMSE performance of the TALS-LTI algorithm versus the SNR (dB) for
$M=3$ transmit antennas, $L=2$ receive antennas, $K=20$ time-blocks and
different RIS elements $N$ and impairments occurrence probability $R_{B}$
assuming i.i.d. Rayleigh fading channels. Figure 5: Number of FLOPS required
by the TALS-LTI algorithm versus the number of the RIS elements $N$ for $M=3$
transmit antennas, $L=2$ receive antennas, impairments occurrence probability
$R_{B}=0.5$ and different time-blocks $K$.

##  V Simulation Results

In this section, we present the simulation results for performance evaluation
of our proposed tensor-based algorithms, in terms of computational complexity
and estimation accuracy of the channels and imperfections, while comparing to
benchmark approaches. The results presented here are averaged over
$\Omega=3\cdot 10^{3}$ independent Monte Carlo runs. Each run corresponds to a
different realization of the involved communication channels, RIS patterns,
impairment parameters and noise. We have designed the RIS pattern matrix
$\mathbf{S}$ as a DFT matrix. The amplitude and phase impairments parameters
embedded in the vector $\mathbf{e}$ and in the matrix $\mathbf{E}$ follow a
uniform distribution between $[0,1]$ and between $[0,2\pi]$, respectively. The
location of the unknown $N_{B}$ impaired elements are assumed to be random
with occurrence probability $R_{B}$, totalizing $N_{B}=NR_{B}$ impaired
elements at the RIS in each simulated setup. The metric used to evaluate the
estimation accuracy is the normalized mean square error (NMSE) between the
true and estimated matrices that provides a relative measure for the
estimation error of the proposed algorithms. For the estimated channel
$\hat{\mathbf{H}}$, we define

$\text{NMSE}(\hat{\mathbf{H}})=\frac{1}{\Omega}\sum_{\omega=1}^{\Omega}\frac{\|\mathbf{H}^{(\omega)}-\hat{\mathbf{H}}^{(\omega)}\|_{\text{F}}^{2}}{\|\mathbf{H}^{(\omega)}\|_{\text{F}}^{2}},$
(62)

where $\mathbf{H}^{(\omega)}$ and $\hat{\mathbf{H}}^{(\omega)}$ denote the
true channel and its estimate both related to the $\omega$-th run,
respectively. Similar definitions as in (62) apply to the estimates of
$\hat{\mathbf{G}}$, $\hat{\mathbf{E}}$ and $\hat{\mathbf{e}}$. Finally, note
that, while we present the simulation results for a set of parameter setting,
we have test the results for a broad range of parameter settings and observed
the same qualitative conclusions as those presented.

###  V-A TALS-LTI Performance

We first examine, in Figs. 4 and 5, the performance of the proposed TALS-LTI
algorithm. These plots show the NMSE and computational complexity in terms of
floating-point operations (FLOPS) assuming i.i.d. Rayleigh fading channels. As
a benchmark, we compare our TALS-LTI algorithm with the method proposed in
[19], which is also a PARAFAC-based algorithm but formulated to the ideal case
in which no impairments affect the RIS elements. Additionally, as a lower-
bound for comparison, we also plot the performance of the clairvoyant LS
estimators of $\hat{\mathbf{H}}$, $\hat{\mathbf{G}}$ and $\hat{\mathbf{e}}$ in
(28), (29) and (30) obtained when the true factor matrices in the right-hand
side of these equations are perfectly known.

From Figs. 4(a) and 4(b), we can observe that the method in [19] is not
suitable to tackle the channel estimation problem when RIS impairments are
present. In contrast, our proposed TALS-LTI algorithm provides accurate
estimates in which the NMSE of the estimated channels decreases linearly when
the SNR increases, and is not sensitive to the number of impaired elements at
the RIS. For instance, the TALS-LTI presents constant gaps when compared to
the lower-bound LS estimator equal to $4$ dB and $3$ dB providing satisfactory
performance in terms of the channel estimation for all simulated SNR range.
Also, in contrast to method in [19], the TALS-LTI accurately estimates the RIS
impairments that are treated as independent variables estimated beyond the
involved channels, as illustrated in Fig. 4(c).

In terms of computational complexity, we can see from Fig. 5 that the TALS-LTI
presents greater complexity than competitor method requiring approximately
$2.2\text{x}$ more FLOPS per iteration. This happens because the TALS-LTI has
one more update step (per iteration) to estimate the RIS imperfections
compared to the method in [19] that does not perform the estimation of
imperfections. Also, the complexity of both approaches increases with the
number of RIS elements $N$ and time-blocks $K$.

(a)

(b)

(c)

Figure 6: NMSE performance of the TALS-STI and HOSVD-STI algorithms versus the
SNR (dB) for $M=3$ transmit antennas, $L=2$ receive antennas, $N=50$ RIS
elements, $K=50$ time-blocks, $P=5$ frames and different impairments
occurrence probability $R_{B}$ assuming a typical mmWave propagation
environment. Figure 7: Number of FLOPS required by the TALS-STI and HOSVD-STI
algorithms versus the number of the RIS elements $N$ for $M=3$ transmit
antennas, $L=2$ receive antennas, $K=50$ time-blocks, impairments occurrence
probability $R_{B}=0.5$ and different frames $P$.

###  V-B TALS-STI and HOSVD-STI Performance

Here, we evaluate the performance of the proposed TALS-STI and HOSVD-STI
algorithms assuming i.i.d. Rayleigh fading channels and a typical mmWave
propagation environment. For the mmWave setup, we assume that the transmitter
and the receiver are equipped with half-wavelength spaced uniform linear
arrays while the RIS has half-wavelength spaced reflecting elements disposed
in a uniform rectangular grid. In this case, the channel matrices are
generated according to the widely used geometric channel model [7], in which
the angles of arrival and departure are randomly generated according to a
uniform distribution. The azimuth and elevation angles are drawn between
$[-\pi/2,\pi/2]$ and $[0,\pi/2]$, respectively. The number of propagation
paths is set to $1$ in the Tx-RIS link and to $2$ in the RIS-Rx link,
respectively. The complex channel gains follow uniform distributions. The
lower-bound LS estimators of $\hat{\mathbf{H}}$, $\hat{\mathbf{G}}$ and
$\hat{\mathbf{E}}$ are obtained similarly to Section V-A, but now from (45),
(46) and (47), respectively.

From Fig. 6, similar conclusions to the results in Fig. 4, in terms of channel
estimation performance, can be made by comparing the proposed TALS-STI and
HOSVD-STI algorithms with the method in [19]. However, more accurate estimates
for the channels and impairments are obtained when the mmWave propagation
scenario is considered. In this experiment, the TALS-STI and HOSVD-STI
present, respectively, constant gaps of approximately $1$ dB and $2$ dB,
compared to the lower-bound LS estimator for the considered range of SNR,
confirming the effectiveness of the proposed tensor-based algorithms in terms
of estimation accuracy. Similar results are obtained for the proposed TALS-LTI
algorithm, but omitted here due to space limitation. Also, we see from Fig.
6(c), that in the low/medium SNR regime, the TALS-STI algorithm outperforms
the HOSVD-STI one in terms of imperfections estimation, indicating that the
TALS-STI algorithm becomes preferable with considerable gain in low SNR when
complexity issues are not taken into account. The tradeoff between
computational complexity and estimation performance is analyzed in the sequel.
The presented result indicates that the proposed algorithms are able to work
with high accuracy under different propagation environments and kinds of
imperfections, i.e., our approaches are general techniques that provide
accurate estimates for different channel models.

Figure 8: Number of iterations required by the TALS-STI algorithm versus the
SNR (dB) for $M=3$ transmit antennas, $L=2$ receive antennas, $K=50$ time-
blocks, impairments occurrence probability $R_{B}=0.5$ and different RIS
elements $N$ and frames $P$.

In Fig. 7, we evaluate the computational complexity (in terms of FLOPS) of the
TALS-STI and HOSVD-STI algorithms. We can see that the HOSVD-STI is less
complex than TALS-STI for different values of RIS elements $N$ and frames $P$,
since it is a closed-form solution for joint channel and imperfections
estimation. The HOSVD-STI algorithm provides a remarkable gain of the order of
$NK$ FLOPS compared to the TALS-STI one. It can be seen that the complexity of
both methods grows when the number of RIS elements $N$ and frames $P$
increases, which is an expected result since the number of entries in
$\mathbf{H}$, $\mathbf{G}$ and $\mathbf{E}$ also increases with $N$ and $P$.
However, in contrast to the HOSVD-STI which is a closed-form solution, for a
complete analysis on the overall complexity of the proposed techniques, Fig. 8
shows the number of iterations necessary for the convergence of the TALS-STI
algorithm considering different values of $N$ and $P$. It can be seen that the
TALS-STI algorithm rapidly converges thanks to the knowledge of the non-
impaired matrix $\mathbf{S}$ that remains fixed during each iteration. Also,
the number of iterations required for the convergence decreases as a function
of the SNR. In the high SNR regime, the convergence of the proposed TALS-STI
algorithm is no more sensitive to $N$, $P$ and SNR values. When the SNR is
higher than $20$ dB, its convergence is quickly achieved within approximately
$5$ iterations.

Figure 9: Average runtime (in seconds) for the TALS-STI and HOSVD-STI
algorithms versus the SNR (dB) for $M=3$ transmit antennas, $L=2$ receive
antennas, $K=50$ time-blocks, $P=5$ frames, impairments occurrence probability
$R_{B}=0.5$ and different RIS elements $N$.

(a)

(b)

(c)

Figure 10: NMSE performance of the TALS-STI and HOSVD-STI algorithms versus
the number of the RIS elements $N$ for $M=3$ transmit antennas, $L=2$ receive
antennas, $K=50$ time-blocks, $P=5$ frames, SNR $=$ 20 dB and different
impairments occurrence probability $R_{B}$ assuming i.i.d. Rayleigh fading
channels.

(a)

(b)

(c)

Figure 11: NMSE performance of the TALS-STI and HOSVD-STI algorithms versus
the impairments occurrence probability $R_{B}$ for $M=3$ transmit antennas,
$L=2$ receive antennas, $N=50$ RIS elements, $K=50$ time-blocks, $P=5$ frames,
SNR $=$ 20 dB assuming i.i.d. Rayleigh fading channels.

In Fig. 9, we evaluate the overall computational complexity of the TALS-STI
and HOSVD-STI algorithms using as metric the average runtime (in seconds).
This metric takes into account the number of iterations for convergence,
making it possible to compare the complexity between iterative and closed-form
solutions. It can be seen that the runtime grows when the number of RIS
elements $N$ increases, confirming the results shown in Fig. 7. The runtime
required by HOSVD-STI is not sensitive to the SNR since it is a closed-form
algorithm. On the other hand, the runtime of the TALS-STI algorithm decreases
with the SNR as it is an iterative algorithm, achieving performance close to
HOSVD-STI algorithm in the high SNR regime thanks to its rapid convergence
(Fig. 8).

We can conclude that TALS-STI and HOSVD-STI algorithms outperform the method
in [19] in terms of both channel estimation accuracy and overall computational
complexity. This is because due to the RIS imperfections, the method in [19]
does not estimate channels properly and has very slow convergence which
directly increases its runtime. However, the TALS-STI can operate under more
flexible choices for the number of time-blocks $K$ for channel estimation
compared to HOSVD-STI (see Fig. 10). From this, we can note that when the
condition in (61) is not satisfied, the TALS-STI significantly outperforms the
HOSVD-STI in terms of estimation accuracy. Therefore, a tradeoff between
overall computational complexity, estimation performance and operation
conditions for the proposed solutions can be observed. Thus, the TALS-STI may
be attractive when more flexible choices for the number of time-blocks for
channel estimation $K$ are required, while the HOSVD-STI is preferred
specially when low processing delay is desired. We can also observe that the
estimation performance degrades as the number of RIS elements $N$ increases,
which is an expected result since the number of channel and imperfections
coefficients to be estimated also increases with $N$.

In our last experiment, shown in Fig. 11, we evaluate the NMSE as a function
of the impairments occurrence probability $R_{B}$. It can be seen that, in
terms of channel estimation, for a broad range of parameter settings, our
proposed approaches are not sensitive to the number of impaired elements at
the RIS. This is because the imperfections matrix is estimated as an
independent variable by both methods. However, in terms of imperfections
estimation, the TALS-STI proves to be more accurate for all considered values
of $R_{B}$. It is also important to note that the proposed TALS-STI and HOSVD-
STI algorithms outperform the method in [19], even when ideal non-impaired RIS
is assumed. However, the gains achieved by our methods become more evident in
the challenging scenario, when the number of impaired elements at the RIS
increases. In this case, the method of [19] is sensitive to $R_{B}$. Thus, if
the number of impaired elements changes, it does not show a stable/predictable
performance. This leads to the need for channel estimation correction and
retransmissions. In contrast, our proposed methods are quite robust to $R_{B}$
variation which simplifies the system design.

##  VI Conclusion

We have proposed different efficient tensor-based algorithms for channel
estimation in RIS-assisted MIMO systems, in which the RIS elements are
affected by real-world imperfections. We resort to the multidimensional
structure of the received signal to solve these non-idealized channel
estimation problems by means of trilinear and quadrilinear PARAFAC models. The
proposed TALS-LTI algorithm solves the problem when static imperfections are
assumed. In a generalized way, we have formulated the TALS-STI and HOSVD-STI
algorithms for the more challenging scenario in which the behavior of the RIS
imperfections is non-static with respect to the channel coherence time. The
TALS-LTI and TALS-STI algorithms are iterative solutions that relax the system
design requirements, operating under more flexible choices for the training
parameters. In contrast, the HOSVD-STI algorithm is a closed-form solution
that has a lower computational complexity compared to the competing ALS-based
solutions, while affording parallel processing. Simulation results illustrate
the high estimation performance of the proposed tensor-based algorithms for
different kinds of imperfections, channel models, and system configurations.
The TALS-STI and HOSVD-STI algorithms present similar channel estimation
performances. However, the TALS-STI is preferable for the imperfections
detection in the low SNR regime and when more flexible choices for training
parameters is required, while the HOSVD-STI is preferred when low processing
delay is desired.

## References

  * [1] P. R. B. Gomes, G. T. de Araújo, B. Sokal, A. L. F. de Almeida, B. Makki, and G. Fodor, “Tensor-based channel estimation for ris-assisted networks operating under imperfections,” 2022. [Online]. Available: https://arxiv.org/abs/2206.03557
  * [2] G. Forecast, “Cisco visual networking index: global mobile data traffic forecast update, 2017–2022,” _Update_ , vol. 2017, p. 2022, 2019.
  * [3] ITU, “IMT traffic estimates for the years 2020 to 2030 ,” ITU-R M, TR 2370-0, Jan. 2015, V.15.0.0.
  * [4] S. Zhang, Q. Wu, S. Xu, and G. Y. Li, “Fundamental green tradeoffs: Progresses, challenges, and impacts on 5G networks,” _IEEE Commun. Surveys Tuts._ , vol. 19, no. 1, pp. 33–56, Feb. 2016.
  * [5] 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, Jul. 2020.
  * [6] Q. Wu, G. Y. Li, W. Chen, D. W. K. Ng, and R. Schober, “An overview of sustainable green 5G networks,” _IEEE Wireless Communications_ , vol. 24, no. 4, pp. 72–80, Aug. 2017.
  * [7] G. Zhou, C. Pan, H. Ren, P. Popovski, and A. L. Swindlehurst, “Channel estimation for ris-aided multiuser millimeter-wave systems,” _IEEE Transactions on Signal Processing_ , vol. 70, pp. 1478–1492, 2022.
  * [8] M. Di Renzo, M. Debbah, D.-T. Phan-Huy, A. Zappone, M.-S. Alouini, C. Yuen, V. Sciancalepore, G. C. Alexandropoulos, J. Hoydis, H. Gacanin _et al._ , “Smart radio environments empowered by reconfigurable ai meta-surfaces: An idea whose time has come,” _EURASIP Journal on Wireless Communications and Networking_ , vol. 2019, no. 1, pp. 1–20, May. 2019.
  * [9] Y.-C. Liang, R. Long, Q. Zhang, J. Chen, H. V. Cheng, and H. Guo, “Large intelligent surface/antennas (LISA): Making reflective radios smart,” _Journal of Communications and Information Networks_ , vol. 4, no. 2, pp. 40–50, Jun 2019.
  * [10] M. Di Renzo, A. Zappone, M. Debbah, M.-S. Alouini, C. Yuen, J. de Rosny, and S. Tretyakov, “Smart radio environments empowered by reconfigurable intelligent surfaces: How it works, state of research, and the road ahead,” _IEEE J. Sel. Areas Commun._ , vol. 38, no. 11, pp. 2450–2525, Jul. 2020\.
  * [11] N. Rajatheva, I. Atzeni, S. Bicais, E. Bjornson, A. Bourdoux, S. Buzzi, C. D’Andrea, J.-B. Dore, S. Erkucuk, M. Fuentes, K. Guan, Y. Hu, X. Huang, J. Hulkkonen, J. M. Jornet, M. Katz, B. Makki, R. Nilsson, E. Panayirci, K. Rabie, N. Rajapaksha, M. Salehi, H. Sarieddeen, S. Shahabuddin, T. Svensson, O. Tervo, A. Tolli, Q. Wu, and W. Xu, “Scoring the terabit/s goal:broadband connectivity in 6G,” Feb. 2021, arXiv:2008.07220v2.
  * [12] E. Basar, M. Di Renzo, J. De Rosny, M. Debbah, M.-S. Alouini, and R. Zhang, “Wireless communications through reconfigurable intelligent surfaces,” _IEEE access_ , vol. 7, pp. 116 753–116 773, Jun. 2019.
  * [13] Q. Wu and R. Zhang, “Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network,” _IEEE Communications Magazine_ , vol. 58, no. 1, pp. 106–112, Nov. 2019.
  * [14] H. Guo, B. Makki, M. Åström, M.-S. Alouini, and T. Svensson, “Dynamic blockage pre-avoidance using reconfigurable intelligent surfaces,” 2022. [Online]. Available: https://arxiv.org/abs/2201.06659
  * [15] 3GPP, “RWS-210300 NR repeaters and reconfigurable intelligent surface,” in _3GPP TSG RAN Rel-18 workshop_. Electronic Meeting, June 28 - July 2, 2021.
  * [16] A. Zappone, M. Di Renzo, F. Shams, X. Qian, and M. Debbah, “Overhead-aware design of reconfigurable intelligent surfaces in smart radio environments,” _IEEE Trans. Wireless Commun._ , vol. 20, no. 1, pp. 126–141, Jan. 2021.
  * [17] A. L. Swindlehurst, G. Zhou, R. Liu, C. Pan, and M. Li, “Channel estimation with reconfigurable intelligent surfaces–a general framework,” _Proceedings of the IEEE_ , pp. 1–27, 2022.
  * [18] G. T. de Araújo, P. R. B. Gomes, A. L. F. de Almeida, G. Fodor, and B. Makki, “Semi-blind joint channel and symbol estimation in irs-assisted multi-user mimo networks,” _IEEE Wireless Commun. Lett._ , pp. 1–1, 2022\.
  * [19] G. T. de Araújo, A. L. de Almeida, and R. Boyer, “Channel estimation for intelligent reflecting surface assisted MIMO systems: A tensor modeling approach,” _IEEE J. Sel. Topics Signal Process._ , vol. 15, no. 3, pp. 789–802, Feb. 2021.
  * [20] 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, Mar. 2020.
  * [21] D. Mishra and H. Johansson, “Channel estimation and low-complexity beamforming design for passive intelligent surface assisted MISO wireless energy transfer,” in _proc. ICASSP “2019”_ , Brighton, UK, Apr.
  * [22] 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. 2019.
  * [23] A. Taha, M. Alrabeiah, and A. Alkhateeb, “Enabling large intelligent surfaces with compressive sensing and deep learning,” _IEEE Access_ , vol. 9, pp. 44 304–44 321, Mar. 2021.
  * [24] J. Chen, Y.-C. Liang, H. V. Cheng, and W. Yu, “Channel estimation for reconfigurable intelligent surface aided multi-user MIMO systems,” Dec. 2019, arXiv:1912.03619v1 [eess.SP].
  * [25] J. He, M. Leinonen, H. Wymeersch, and M. Juntti, “Channel estimation for RIS-aided mmwave MIMO systems,” in _proc. GLOBECOM “2020”_ , Taipei, Taiwan, Feb.
  * [26] N. Sidiropoulos, G. Giannakis, and R. Bro, “Blind PARAFAC receivers for DS-CDMA systems,” _IEEE Trans. Signal Process._ , vol. 48, no. 3, pp. 810–823, Mar. 2000.
  * [27] H. Chen, F. Ahmad, S. Vorobyov, and F. Porikli, “Tensor decompositions in wireless communications and mimo radar,” _IEEE J.Sel. Topics in Signal Process._ , vol. 15, no. 3, pp. 438–453, Feb. 2021.
  * [28] N. D. Sidiropoulos, L. De Lathauwer, X. Fu, K. Huang, E. E. Papalexakis, and C. Faloutsos, “Tensor decomposition for signal processing and machine learning,” _IEEE Trans. Signal Process._ , vol. 65, no. 13, pp. 3551–3582, Feb. 2017.
  * [29] L. Wei, C. Huang, G. C. Alexandropoulos, and C. Yuen, “Parallel factor decomposition channel estimation in RIS-assisted multi-user MISO communication,” in _proc. SAM “2020”_ , Hangzhou, China, Jun.
  * [30] G. T. de Araújo and A. L. de Almeida, “PARAFAC-based channel estimation for intelligent reflective surface assisted MIMO system,” in _proc. SAM “2020”_ , Hangzhou, China, Jun.
  * [31] G. T. de Araújo and A. L. F. de Almeida, “Channel estimation for MIMO system assisted by intelligent reflective surface,” in _proc. SBrT “2020”_ , Santa Catarina, Brazil, Nov.
  * [32] X. Qian, M. Di Renzo, J. Liu, A. Kammoun, and M.-S. Alouini, “Beamforming through reconfigurable intelligent surfaces in single-user MIMO systems: SNR distribution and scaling laws in the presence of channel fading and phase noise,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 1, pp. 77–81, Sep. 2021.
  * [33] M.-A. Badiu and J. P. Coon, “Communication through a large reflecting surface with phase errors,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 2, pp. 184–188, Feb. 2020.
  * [34] S. Zhou, W. Xu, K. Wang, M. Di Renzo, and M.-S. Alouini, “Spectral and energy efficiency of IRS-assisted MISO communication with hardware impairments,” _IEEE wireless commun. lett._ , vol. 9, no. 9, pp. 1366–1369, Sep. 2020.
  * [35] K. Zhi, C. Pan, H. Ren, and K. Wang, “Uplink achievable rate of intelligent reflecting surface-aided millimeter-wave communications with low-resolution ADC and phase noise,” _IEEE Wireless Commun. Lett._ , vol. 10, no. 3, pp. 654–658, Mar. 2021.
  * [36] B. Li, Z. Zhang, Z. Hu, and Y. Chen, “Joint array diagnosis and channel estimation for RIS-aided mmwave MIMO system,” _IEEE Access_ , vol. 8, pp. 193 992–194 006, Oct. 2020.
  * [37] B. Li, Z. Zhang, and Z. Hu, “Channel estimation for reconfigurable intelligent surface-assisted multiuser mmwave MIMO system in the presence of array blockage,” _Transactions on Emerging Telecommunications Technologies_ , vol. 32, no. 11, p. e4322, Jun. 2021.
  * [38] T. G. Kolda and B. W. Bader, “Tensor decompositions and applications,” _SIAM review_ , vol. 51, no. 3, pp. 455–500, Aug. 2009.
  * [39] L. De Lathauwer, B. De Moor, and J. Vandewalle, “A multilinear singular value decomposition,” _SIAM Journal on Matrix Analysis and Applications_ , vol. 21, no. 4, pp. 1253–1278, Apr. 2000.
  * [40] S. Abeywickrama, R. Zhang, Q. Wu, and C. Yuen, “Intelligent reflecting surface: Practical phase shift model and beamforming optimization,” _IEEE Trans. Commun._ , vol. 68, no. 9, pp. 5849–5863, Sep. 2020.
  * [41] Z. Yigit, E. Basar, and I. Altunbas, “Low complexity adaptation for reconfigurable intelligent surface-based MIMO systems,” _IEEE Commun. Lett._ , vol. 24, no. 12, pp. 2946–2950, Dec. 2020.
  * [42] B. Sokal, P. R. B. Gomes, A. L. F. de Almeida, B. Makki, and G. Fodor, “IRS phase-shift feedback overhead-aware model based on rank-one tensor approximation,” 2022. [Online]. Available: https://arxiv.org/abs/2205.12024
  * [43] ——, “Reducing the control overhead of intelligent reconfigurable surfaces via a tensor-based low-rank factorization approach,” 2022. [Online]. Available: https://arxiv.org/abs/2206.05341
  * [44] R. Bro, “Multi-way analysis in the food industry-models, algorithms, and applications,” in _MRI, EPG and EMA,” Proc ICSLP 2000_. Citeseer, 1998.
  * [45] X.-F. Gong and Q.-H. Lin, “Spatially constrained parallel factor analysis for semi-blind beamforming,” in _2011 Seventh International Conference on Natural Computation_ , vol. 1. IEEE, 2011, pp. 416–420.
  * [46] Z. Wang, C. Cai, F. Wen, and D. Huang, “A quadrilinear decomposition method for direction estimation in bistatic MIMO radar,” _IEEE Access_ , vol. 6, pp. 13 766–13 772, Mar. 2018.
  * [47] 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. ICASSP “2020”_ , Barcelona, Spain, May.