# Hybrid Vector Message Passing for Cascaded Channel Estimation in RIS-Aided Multi-User MIMO-OFDM Systems

Wenjun Jiang, Xiaojun Yuan, Marco Di Renzo

► To cite this version:

Wenjun Jiang, Xiaojun Yuan, Marco Di Renzo. Hybrid Vector Message Passing for Cascaded Channel Estimation in RIS-Aided Multi-User MIMO-OFDM Systems. IEEE Transactions on Wireless Communications, 2025, 24 (5), pp.4174–4189. 10.1109/TWC.2025.3536283. hal-05378144

HAL Id: hal-05378144

https://hal.science/hal-05378144v1

Submitted on 28 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Hybrid Vector Message Passing for Cascaded Channel Estimation in RIS-Aided Multi-User MIMO-OFDM Systems

Wenjun Jiang<sup>D</sup>, Graduate Student Member, IEEE, Xiaojun Yuan<sup>D</sup>, Senior Member, IEEE, and Marco Di Renzo, Fellow, IEEE

**Abstract**— This paper investigates the fundamental problem of cascaded channel estimation in reconfigurable intelligent surface (RIS) aided multi-user (MU) multiple-input multiple-output (MIMO) orthogonal frequency division multiplexing (OFDM) systems. To address the high pilot overhead required by existing methods, we introduce appropriate auxiliary variables to decompose the received signal model into a subcarrier-wise bilinear sub-model and two linear sub-models with respect to the MU-to-RIS and RIS-to-base-station (BS) cascaded channels. The proposed model decomposition maintains the matrix factorization structure of the two cascaded channels and avoids the problem of parameter expansion in exiting methods. In the two linear sub-models, in addition, the cascaded channels are transformed into the delay-angle domain to leverage the joint sparsity. Based on the established models, we consider the problem of simultaneously estimating the MU-to-RIS and RIS-to-BS channel matrices as a bilinear estimation problem. To tackle this problem, we formulate a Bayesian inference framework and develop a hybrid message-passing (HVMP) algorithm to achieve approximate Bayesian inference by leveraging the Bethe method. Notably, the HVMP algorithm infers the two cascaded channels iteratively, where the covariance of each cascaded channel matrix is estimated to characterize the correlation of the matrix elements. Simulation results show that the proposed algorithm achieves accurate channel estimation with low pilot overhead while state-of-the-art baseline schemes exhibit poor performance. Furthermore, the proposed algorithm can approach the estimation oracle bound

The work of Xiaojun Yuan was

supported in part by Sichuan Science and Technology Program under Grant 2024ZYD0036 and in part by the General Program of the National Natural Science Foundation of China under Grant 62071090. The work of Marco Di Renzo was supported in part by European Union through the Horizon Europe Project COVER under Grant 101086228, in part by the Horizon Europe Project UNITE under Grant 101129618, in part by the Horizon Europe Project INSTINCT under Grant 101139161, in part by the Horizon Europe Project TWIN6G under Grant 101182794, in part by the Agence Nationale de la Recherche (ANR) through the France 2030 Project ANR-PEPR Networks of the Future under Grant NF-PERSEUS 22-PEFT-004, and in part by the CHIST-ERA Project PASSIONATE under Grant CHIST-ERA-22-WAI-04 and Grant ANR-23-CHR4-0003-01. The associate editor coordinating the review of this article and approving it for publication was I. Bergel. (Corresponding author: Xiaojun Yuan.)

Wenjun Jiang and Xiaojun Yuan are with the National Key Laboratory of Wireless Communications, University of Electronic Science and Technology of China, Chengdu 611731, China (e-mail: wjjiang@std.uestc.edu.cn; xjyuan@uestc.edu.cn).

Marco Di Renzo is with Université Paris-Saclay, CNRS, Centrale-Supélec, Laboratoire des Signaux et Systèmes, 91192 Gif-sur-Yvette, France, and also with King's College London, Centre for Telecommunications Research—Department of Engineering, WC2R 2LS London, U.K. (e-mail: marco.di-renzo@universite-paris-saclay.fr; marco.di_renzo@kcl.ac.uk).

of the MU-to-RIS (or RIS-to-BS) channel which assumes perfect knowledge of the RIS-to-BS (or MU-to-RIS) channel.

**Index Terms**— Reconfigurable intelligent surface, multiple-input multiple-output, orthogonal frequency-division multiplexing, cascaded channel estimation, message passing.

## I. INTRODUCTION

RECONFIGURABLE intelligent surface (RIS) has emerged as a promising technology to empower sixth-generation (6G) wireless communications [1], [2], [3]. A typical RIS is a planar meta-surface comprising a large number of nearly-passive elements, where each element induces an independent phase shift to the incident electromagnetic wave consuming ultra-low power. By jointly optimizing the RIS phase shifts, RIS provides a highly reliable communication link between transceivers, a.k.a., nearly-passive beamforming [4], [5], [6]. Nearly-passive beamforming significantly enhances communication coverage, spectral efficiency, connection density, etc. Furthermore, RIS phase shifts can be adjusted to implicitly convey information that is embedded in the incident electromagnetic wave, a.k.a., nearly-passive information transfer [7], [8], [9]. The nearly-passive information transfer provides a new dimension for information delivery, enabling, for instance, wireless signaling interaction between an RIS and a base station (BS), green information transfer for connected IoT devices without radio-frequency chains, etc.

Efficient nearly-passive beamforming and/or nearly-passive information transfer require the knowledge of channel state information (CSI). However, the acquisition of CSI in an RIS-aided system is a challenging task, since the RIS comprises nearly-passive elements and is not capable of baseband signal processing, and it is not endowed with power amplifiers. More importantly, the RIS-related channels have a cascaded structure that is made of the link from the users to the RIS and the link from the RIS to the BS [10]. With the increase of the number of RIS elements, there is a large number of cascaded channel coefficients to be estimated, incurring an unprecedented increase in system overhead. Hence, it is crucial to develop an efficient channel estimation algorithm achieving a low pilot overhead and a high estimation accuracy.

Recently, cascaded channel estimation methods for narrow-band multiple-input multiple-output (MIMO) systems have
became a cutting-edge and actively researched topic [11], [12], [13], [14], [15]. Similar to a conventional MIMO channel, the cascaded channels exhibit angle domain sparsity in both the user-to-RIS link and the RIS-to-BS link [13], [14], [15]. By exploiting the angle domain sparsity, the number of non-zero channel coefficients is significantly reduced for pilot saving. When targeting broadband orthogonal frequency division multiplexing (OFDM) systems, a new challenge arises due to the channel frequency selectivity [16], [17], [18], [19], [20], [21], [22]. The authors of [17] and [18] propose to exploit the limited channel delay spread and introduce a new delay domain linear representation of the user-to-RIS and RIS-to-BS cascaded channels. With this linear representation, the problem of the cascaded channel estimation is converted into a linear inverse problem, which is simpler to solve but requires the estimation of additional parameters. The authors of [19], [20], and [21] leverage sporadic multi-paths in a millimeter wave system and propose to estimate the channel model parameters, e.g., the path delays/angles, based on a ray-tracing channel model. In a ray-tracing model, each RIS element with a size comparable to a wavelength acts as a diffuse scatterer [1]. Then, a path of the overall user-RIS-BS channel corresponds to the combination of a path in the user-to-RIS link and a path in the RIS-to-BS link. Therefore, the total number of parameters required in the channel representation used in [19], [20], and [21] is proportional to the product of the numbers of paths in the two cascaded links. Another research direction considered a hybrid nearly passive/active RIS architecture [22]. Based on the received signal from the active RIS elements, a deep learning-assisted compressive sensing method was proposed in [22] to estimate the user-to-RIS channel.

In RIS-aided OFDM systems, existing channel estimation methods are mostly developed for single-user systems. However, fifth generation (5G) BSs operate using MIMO-OFDM with MU transmission on the same time-frequency resource block for attaining high spectral efficiency. It is crucial to study RIS-assisted MIMO-OFDM systems in the MU scenario. To the best of our knowledge, the problem of cascaded channel estimation has not been fully addressed yet in an RIS-aided MU MIMO-OFDM system. Extending the methods in [17], [18], [19], [20], and [21] to the MU scenario is possible, but the corresponding pilot overhead will become prohibitive, since an increase of the number of users leads to a sharp increase of the number of cascaded channel coefficients. A recent work [23] exploits the property that all the users share a common RIS-to-BS channel and proposes a two-phase scheme. However, the two-phase scheme in [23] relies on a known line-of-sight (LOS) path in the RIS-to-BS link (as a prior) to assist the cascaded channel estimation.

In this paper, we aim to develop a new cascaded channel estimation method that can achieve a low pilot overhead and a high estimation accuracy for application to RIS-aided MU MIMO-OFDM systems. To this end, we first introduce appropriate auxiliary variables to decompose the received signal model into three sub-models, i.e., a subcarrier-wise bilinear sub-model and two linear sub-models with respect to the MU-to-RIS and RIS-to-BS channels. The proposed model

<table><caption>TABLE I
COMPARISON OF EXISTING METHODS FOR CASCADED CHANNEL ESTIMATION</caption><thead><tr><td>References</td><td>Consider multi-antenna at BS?</td><td>Consider OFDM?</td><td>Consider multi-user?</td></tr></thead><tbody><tr><td>[12]–[15]</td><td>✓</td><td>×</td><td>✓</td></tr><tr><td>[16], [17]</td><td>×</td><td>✓</td><td>×</td></tr><tr><td>[19]–[22]</td><td>✓</td><td>✓</td><td>×</td></tr><tr><td>[18], [23]</td><td>✓</td><td>✓</td><td>✓</td></tr></tbody></table>

decomposition maintains the matrix factorization structure of the MU-to-RIS and RIS-to-BS cascaded channels, which avoids the parameter expansion problem that characterizes the linearization approach in [17] and [18] and the path combination problem in [19], [20], and [21]. Furthermore, in the two linear sub-models, we transform the two cascaded channel matrices into the delay-angle domain to leverage the joint channel sparsity. In the proposed method, the number of non-zero channel coefficients is proportional to the sum of the numbers of paths in the two cascaded links, and is significantly smaller than the product of the numbers of paths in [19], [20], and [21].

With the established signal model, the main challenge is to simultaneously estimate the MU-to-RIS and RIS-to-BS channel matrices from the received signal function as a bilinear function of the two cascaded channels. To address this issue, we formulate the bilinear estimation problem into a Bayesian inference problem. In the Bayesian inference framework, we employ the Bethe method [24] to develop a hybrid vector message passing (HVMP) algorithm. Different from a sequential estimation of the two cascaded channels in [23], the HVMP algorithm iteratively estimates the MU-to-RIS and RIS-to-BS channel matrices. In each iteration, the HVMP algorithm is designed to output an estimate of each channel covariance matrix to characterize the correlation of the matrix entries, by borrowing the idea of vector-form message-passing in [25]. This correlation is usually not considered in existing message-passing designs for cascaded channel estimation [12], [13]. We conduct extensive simulations to evaluate the performance of the proposed algorithm, which show that: (a) the HVMP algorithm converges within a few iterations; (b) the HVMP algorithm achieves accurate channel estimation with a low pilot overhead while various state-of-the-art baseline schemes exhibit poor performance or even fail to work properly; (c) the HVMP algorithm can approach the estimation oracle bound, with the oracle-bound for the MU-to-RIS channel (or the RIS-to-BS channel) assuming perfect knowledge of the RIS-to-BS channel (or the MU-to-RIS channel).

*Organization:* In Sec. II, we introduce the signal model for the considered RIS-aided multi-user MIMO-OFDM system. In Sec. III, we decompose the signal model into three sub-models. Based on the sub-models, we formulate the problem of cascaded channel estimation into a Bayesian inference problem. In Sec. IV, we develop a message-passing algorithm based on the proposed Bayesian inference model. In Sec. V, we present numerical results. In Sec. VI, we conclude this paper.
<figure><img src="image_2.png" /><figcaption>Fig. 1. An RIS-aided multi-user MIMO-OFDM system, where the RIS is placed on the yz plane.</figcaption></figure>

Notation: We use bold capital letters (e.g., **X**) for matrices and bold lowercase letters (e.g., **x**) for vectors. (·)T, (·)*, and (·)H denote the transpose, conjugate, and conjugate transpose, respectively. δ(·) denotes the Dirac delta function. The matrix I_N denotes the identity matrix with size N × N. ⊗ and ⊙ denote the Kronecker product and the point-wise product, respectively. We use diag(**X**(n)) to denote a block diagonal matrix with the n-th diagonal block being **X**(n), n = 1, ..., N. We use diag(**x**) to denote a diagonal matrix with the i-th diagonal entry being the i-th element of **x**. tr(·) and vec(·) denote the trace and vectorization operators, respectively. (·)⁻¹ denotes the Moore-Penrose inverse. For a random vector **x** (or matrix **X**), we denote its probability density function (pdf) by p(**x**) (or p(**X**)). E[·] and cov(·) denote the expectation and covariance operators, respectively. The pdf of a complex Gaussian random variable x ∈ C with mean m and variance v is denoted by CN(x; m, v) = exp(-|x - m|^2/v)/(πv). The probability density function (pdf) of a complex Gaussian random vector **x** ∈ C^N with mean **m** and precision matrix **Σ** is denoted by CN(**x**; **m**, **Σ**) = |**Σ**|exp(-(**x** - **m**)ᵀ**Σ**(**x** - **m**)))/π^N. The pdf of a complex Gaussian random matrix **X** ∈ C^{N₁×N₂} with mean **M** and precision matrix **Σ**^{N₁×N₁} is denoted by CN(**X**; **M**, **Σ**) = |**Σ**|^{N₁}exp(-tr(**X**-**M**)ᵀ**Σ**(**X**-**M**))/π^{N₁N₂}.

## II. SYSTEM MODEL

### A. RIS-Aided Multi-User MIMO-OFDM System

Consider the RIS-aided MIMO-OFDM system illustrated in Fig. 1, where the BS is equipped with a uniform linear array (ULA) of M antennas, and K single-antenna users communicate simultaneously with the BS over J continuous subcarriers. The communication between the users and the BS is assisted by an RIS consisting of N nearly-passive reflecting elements arranged in the form of an N₁×N₂ uniform rectangular array (URA) with N₁N₂ = N. Each RIS element induces an independent phase shift on the incident signal.

We ignore mutual coupling and multiple reflections at the RIS [26]. As for the MU-to-RIS channel, we assume a block-fading channel model, i.e., the channel remains constant in a transmission block and varies independently across blocks. Denote by g_k^UR(j) ∈ C^N the channel from the k-th user to the RIS at the j-th subcarrier in a transmission block. The channel vector g_k^UR(j) is given by [27] and [28]

$$g_k^\text{UR}(j) = \sum_{\ell=1}^{P_k^\text{UR}} \beta_{k,\ell}^\text{UR} \mathbf{a}_R (\phi_{k,\ell}^\text{UR}, \theta_{k,\ell}^\text{UR}) e^{-i2\pi \Delta f \tau_{k,\ell}^\text{UR} j}, \quad (1)$$

where P_k^UR is the number of multi-paths of user k; β_{k,ℓ}^UR is the complex-valued channel coefficient of path l; Δf is the subcarrier spacing; τ_{k,ℓ}^UR is the l-th path delay; φ_{k,ℓ}^UR is the ℓ-th zenith angle-of-arrival (AOA) and θ_{k,ℓ}^UR is the ℓ-th azimuth AOA; a_R(φ, θ) ∈ C^N is the array response of the URA given by [29]

$$\mathbf{a}_R(\phi, \theta) = \mathbf{f}_{N_1}(\sin \theta \cos \phi) \otimes \mathbf{f}_{N_2}(\sin \phi),$$

with

$$\mathbf{f}_N(x) \triangleq \frac{1}{\sqrt{N}} [1, e^{-i2\pi \frac{d}{\lambda} x}, \dots, e^{-i2\pi \frac{d}{\lambda} (N-1)x}]^T,$$

where λ is the wavelength and d is the inter-element spacing. We assume d/λ = 1/2.

Next, we model the RIS-to-BS channel G^RB(j). The locations of the RIS and the BS are generally static after deployment [13]. Then, the channel links originating from static scattering clusters (e.g., the building and the road in Fig. 1) remain constant over a duration significantly longer than a transmission block. Meanwhile, a small number of channel links originating from moving scattering clusters (e.g., the bus in Fig. 1) vary across transmission blocks. Thus, we divide the RIS-to-BS channel into static and dynamic components as

$$G^{\text{RB}}(j) = \bar{G}^{\text{RB}}(j) + \tilde{G}^{\text{RB}}(j), \quad (2)$$

with

$$\tilde{G}^{\text{RB}}(j) = \sum_{\ell=1}^{\tilde{P}^{\text{RB}}} \beta_{\ell}^{\text{RB}} \mathbf{a}_B(\tilde{\psi}_{\ell}^{\text{RB}}) \mathbf{a}_R^H(\tilde{\phi}_{\ell}^{\text{RB}}, \tilde{\theta}_{\ell}^{\text{RB}}) e^{-i2\pi \Delta f \tau_{\ell}^{\text{RB}} j}, \quad (3)$$

$$\tilde{G}^{\text{RB}}(j) = \sum_{\ell=1}^{\tilde{P}^{\text{RB}}} \tilde{\beta}_{\ell}^{\text{RB}} \mathbf{a}_B(\tilde{\psi}_{\ell}^{\text{RB}}) \mathbf{a}_R^H(\tilde{\phi}_{\ell}^{\text{RB}}, \tilde{\theta}_{\ell}^{\text{RB}}) e^{-i2\pi \Delta f \tilde{\tau}_{\ell}^{\text{RB}} j}, \quad (4)$$

where P^RB and P^RB are the numbers of multi-path components of the static and dynamic channels, respectively; φ^RB, θ^RB, and ψ^RB are the ℓ-th azimuth angle-of-departure (AOD), the ℓ-th zenith AOD, and the ℓ-th azimuth AOA of G^RB(j), respectively; φ̃^RB, θ̃^RB, and ψ̃^RB have a similar meaning for G^RB(j); a_B(θ) ∈ C^N is the array response of the ULA given by a_B(θ) = f_M(cos θ). {β^RB, ψ^RB, φ^RB, θ^RB, τ^RB}l remain constants for a relatively long time, and {β̃^RB, ψ̃^RB, φ̃^RB, θ̃^RB, τ̃^RB}l change between transmission blocks. Then, the static channel G^RB(j), ∀j, can be easily obtained by averaging over historical estimates of
$G^{\text{RB}}(j), \forall j$, at the initial stage when the communication link is established. Thus, we assume that $\bar{G}^{\text{RB}}(j), \forall j$, is known.<fnref n="1" />

We assume that the MU-to-BS channel is accurately estimated in advance<fnref n="2" /> and focus on the estimation of $g_k^{\text{UR}}(j)$ and $\bar{G}^{\text{RB}}(j)$. At the beginning of each transmission block, $T$ OFDM symbols are used for pilot transmission. On the $t$-th symbol and the $j$-th subcarrier, we denote by $x_{k,t}(j) \in \mathbb{C}$ the pilot symbol of user $k$. Following the 5G standard TS 38.211 [32], we adopt an orthogonal pilot assignment, i.e., if $x_{k,t}(j) \neq 0, \forall t, j$ for user $k$, then $x_{i,t}(j) = 0$ for any user $i \neq k$. The average power of the pilot symbols is denoted by $P_x$. As for the RIS configuration, we denote by $\Theta \in \mathbb{C}^{N \times T}$ the phase of the reflecting coefficients whose magnitude is assumed to be one. Also, we assume that the length of the cyclic prefix (CP) is longer than the maximum channel delay spread. After removing the CP and applying the discrete Fourier transform (DFT), the received signal at the $j$-th subcarrier is expressed as [5]

$$Y(j) = G^{\text{RB}}(j) (\Theta \odot (G^{\text{UR}}(j)X(j))) + N(j), \quad \forall j, \quad (5)$$

where $G^{\text{UR}}(j) = [g_1^{\text{UR}}(j), \dots, g_K^{\text{UR}}(j)] \in \mathbb{C}^{N \times K}$ is the MU-to-RIS channel; $X(j) \in \mathbb{C}^{K \times T}$ is the pilot matrix of $K$ users with the $(k, t)$-th entry being $x_{k,t}(j)$; $N(j) \in \mathbb{C}^{M \times T}$ is an additive white Gaussian noise (AWGN) matrix with elements independently drawn from $\mathcal{CN}(0, \sigma_w^2)$.

## B. Delay-Angle Domain Representation

In (1), the path delays $\tau_{k,l}^{\text{UR}}, \forall k, l$ are typically smaller than the symbol duration $1/\Delta_f$. We introduce the equivalent delay-space domain channel $h_{\text{ds},k}^{\text{UR}} = [h_{\text{ds},k}^{\text{UR}}(1)^T, \dots, h_{\text{ds},k}^{\text{UR}}(L_1)^T] \in \mathbb{C}^{L_1 N}$ with $h_{\text{ds},k}^{\text{UR}}(l) \in \mathbb{C}^N$ and delay tap number $L_1 \ll J$. Then, we represent the frequency-space domain channel $g_k^{\text{UR}} = [g_k^{\text{UR}}(1)^T, \dots, g_k^{\text{UR}}(J)^T] \in \mathbb{C}^{JN}$ with

$$g_k^{\text{UR}} = \mathbf{F}_1 \mathbf{h}_{\text{ds},k}^{\text{UR}}, \quad (6)$$

where $\mathbf{F}_1 = \mathbf{D}_{L_1} \otimes \mathbf{I}_N$, and $\mathbf{D}_{L_1}$ consists of $L_1$ selected columns of the $J \times J$ DFT matrix. In practice, the selected $L_1$ columns comprise the $L_1^{\text{first}}$ first columns and the $L_1^{\text{last}}$ last columns of the DFT matrix with $L_1 = L_1^{\text{first}} + L_1^{\text{last}}$, where $L_1^{\text{first}}$ and $L_1^{\text{last}}$ are determined according to the maximum delay spread and the energy leakage, respectively. Due to the limited scattering, the delays $\{\tau_{k,l}^{\text{UR}}\}$ and angles $\{\phi_{k,l}^{\text{UR}}, \theta_{k,l}^{\text{UR}}\}$ are clustered. In other words, $h_{\text{ds},k}^{\text{UR}}$ exhibits sparsity when transformed into the delay-angle domain. We define the delay-angle domain channel $h_k^{\text{UR}} = [h_k^{\text{UR}}(1)^T, \dots, h_k^{\text{UR}}(L_1)^T]^T \in \mathbb{C}^{L_1 N}$ and express $h_{\text{ds},k}^{\text{UR}}$ as [33]

$$\mathbf{h}_{\text{ds},k}^{\text{UR}} = \mathbf{U}_R \mathbf{h}_k^{\text{UR}}, \quad (7)$$

with $\mathbf{U}_R = \mathbf{I}_{L_1} \otimes \mathbf{A}_R$, where $\mathbf{A}_R = [\mathbf{f}_{N_1}(-1), \dots, \mathbf{f}_{N_1}(\frac{N_1-2}{N_1})] \otimes [\mathbf{f}_{N_2}(-1), \dots, \mathbf{f}_{N_2}(\frac{N_2-2}{N_2})]$

<figure><img src="image_3.png" /><figcaption>Fig. 2. An example for the delay-angle domain representation of the user-to-RIS channel. $N = 10$, $J = 10$ and $L_1 = 5$.</figcaption></figure>

is an angle domain transformation matrix for the URA obtained by uniformly sampling $f_{N_1}(\cdot)$ and $f_{N_2}(\cdot)$ with spacing $\frac{2}{N_1}$ and $\frac{2}{N_2}$, respectively. Note that the $n$-th entry of $h_k^{\text{UR}}(l)$ represents the channel coefficient on the $l$-th delay tap and the $n$-th AOA grid point. Fig. 2 gives an example for the user-to-RIS channel representation in the delay-angle domain, where we set $N = 10$, $J = 10$, and $L_1 = 5$. The dashed squares represent the subgrid to be removed, which involves $J - L_1$ delay taps in the middle of the OFDM symbol. In the selected $L_1 \times N$ delay-angle grid, the empty and shaded squares correspond to the zero and non-zero channel coefficients, respectively.

Similar to (6), we represent the *dynamic* channel component $\tilde{G}^{\text{RB}} = [\tilde{G}^{\text{RB}}(1), \dots, \tilde{G}^{\text{RB}}(J)]^T \in \mathbb{C}^{JN \times M}$ with

$$\tilde{G}^{\text{RB}} = \mathbf{F}_2 \mathbf{H}_{\text{ds}}^{\text{RB}}, \quad (8)$$

where $\mathbf{F}_2 = \mathbf{D}_{L_2} \otimes \mathbf{I}_N$ with $\mathbf{D}_{L_2}$ consisting of $L_2$ selected columns of the $J \times J$ DFT matrix, and $\mathbf{H}_{\text{ds}}^{\text{RB}} \in \mathbb{C}^{L_2 N \times M}$ represents the delay-space domain channel. Then, similar to (7), the delay-angle domain representation of $\mathbf{H}_{\text{ds}}^{\text{RB}}$ is given by

$$\mathbf{H}_{\text{ds}}^{\text{RB}} = \mathbf{U}_B \mathbf{H}^{\text{RB}} \mathbf{A}_B^H, \quad (9)$$

where $\mathbf{A}_B = [\mathbf{f}_M(-1), \dots, \mathbf{f}_M(\frac{M-2}{M})]$ is an angle domain transformation matrix for the ULA, and $\mathbf{U}_B = \mathbf{I}_{L_2} \otimes \mathbf{A}_R$. Note that the $(n, m)$-th entry of $\mathbf{H}^{\text{RB}}(l)$ represents the channel coefficient on the $l$-th delay tap, the $n$-th AOD grid point, and the $m$-th AOA grid point.

## III. PROBLEM FORMULATION

Given the noisy observation $\mathbf{Y}$ in (5), we aim to estimate the cascaded channels $\{\mathbf{h}_k^{\text{UR}}\}_k$ and $\mathbf{H}^{\text{RB}}$ by exploiting the delay-angle domain representations in (6)-(9). By introducing auxiliary variables, we establish an equivalent signal model as

$$\mathbf{Y}(j) = \mathbf{G}^{\text{RB}}(j)\mathbf{R}(j) + \mathbf{N}(j), \quad \forall j, \quad (10a)$$

$$\mathbf{r}_t = \sum_{k=1}^{K} \text{diag}(\mathbf{x}_{k,t} \otimes \theta_t) \mathbf{F}_1 \mathbf{U}_R \mathbf{h}_k^{\text{UR}}, \quad \forall t, \quad (10b)$$

$$\mathbf{G}^{\text{RB}} = \bar{\mathbf{G}}^{\text{RB}} + \mathbf{F}_2 \mathbf{U}_B \mathbf{H}^{\text{RB}} \mathbf{A}_B^H, \quad (10c)$$

<footnote n="1">¹The assumption of knowing $\bar{G}^{\text{RB}}(j)$ is made without loss of generality, since any estimation error that originates from channel averaging can be absorbed into the dynamic channel $G^{\text{RB}}(j)$.</footnote>

<footnote n="2">²To estimate the users-to-BS channel, we can adjust the RIS to absorb incident electromagnetic waves [30], [31] and adopt conventional channel estimation methods.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 3. A factor graph representation of the probability factorization in (13).</figcaption></figure>

and $p(\mathbf{H}^{\mathrm{UR}}|\mathbf{Y})$, respectively. As for each factor node and the neighboring variable nodes, we introduce $b_{f_Y}(\mathbf{R}, \mathbf{G}^{\mathrm{RB}})$, $b_{f_R}(\mathbf{R}, \mathbf{H}^{\mathrm{UR}})$, $b_{f_G^{\mathrm{RB}}}(\mathbf{G}^{\mathrm{RB}}, \mathbf{H}^{\mathrm{RB}})$, $b_{f_H^{\mathrm{UR}}}(\mathbf{H}^{\mathrm{UR}})$, and $b_{f_H^{\mathrm{RB}}}(\mathbf{H}^{\mathrm{RB}})$ to approximate $p(\mathbf{R}, \mathbf{G}^{\mathrm{RB}}|\mathbf{Y})$, $p(\mathbf{R}, \mathbf{H}^{\mathrm{UR}}|\mathbf{Y})$, $p(\mathbf{G}^{\mathrm{RB}}, \mathbf{H}^{\mathrm{RB}}|\mathbf{Y})$, $p(\mathbf{H}^{\mathrm{UR}}|\mathbf{Y})$, and $p(\mathbf{H}^{\mathrm{RB}}|\mathbf{Y})$, respectively. These two types of density functions, termed as *beliefs*, fulfill the marginalization consistency as

$$b_v(\mathbf{v}) = \int b_f(\mathbf{u})d\mathbf{u}\setminus\mathbf{v}, \qquad (15)$$

where $f \in \mathcal{F} = \{f_Y, f_R, f_G^{\mathrm{RB}}, f_H^{\mathrm{UR}}, f_H^{\mathrm{RB}}\}$; $\mathbf{v} \in \mathcal{V} = \{\mathbf{R}, \mathbf{G}^{\mathrm{RB}}, \mathbf{H}^{\mathrm{UR}}, \mathbf{H}^{\mathrm{RB}}\}$; $\mathbf{u}$ denotes all the variables associated with the factor node $f$; $\mathbf{u} \setminus \mathbf{v}$ denotes $\mathbf{u}$ except $\mathbf{v}$. Given $b_v(\mathbf{v})$ and $b_f(\mathbf{u})$ along with the tree structure of the factor graph, the optimal $b(\mathbf{H}^{\mathrm{UR}}, \mathbf{H}^{\mathrm{RB}}, \mathbf{G}^{\mathrm{RB}}, \mathbf{R})$ (which minimizes (14)) fulfills [24]

$$b(\mathbf{H}^{\mathrm{UR}}, \mathbf{H}^{\mathrm{RB}}, \mathbf{G}^{\mathrm{RB}}, \mathbf{R}) = \frac{\prod_{f \in \mathcal{F}} b_f(\mathbf{u})}{\prod_{v \in \mathcal{V}} b_v(\mathbf{v})}. \qquad (16)$$

Substituting (13) and (16) into (14) and applying some manipulations, we express the variational free energy as a function of the beliefs $\{b_v(\mathbf{v})\}$ and $\{b_f(\mathbf{u})\}$, i.e.,

$$F(\{b_v(\mathbf{v})\}, \{b_f(\mathbf{u})\}) = \sum_{f \in \mathcal{F}} \int b_f(\mathbf{u}) \ln \frac{b_f(\mathbf{u})}{f(\mathbf{u})} d\mathbf{u} \\ - \sum_{v \in \mathcal{V}} \int b_v(\mathbf{v}) \ln b_v(\mathbf{v}) d\mathbf{v}. \qquad (17)$$

We aim to optimize $\{b_v(\mathbf{v})\}$ and $\{b_f(\mathbf{u})\}$ in (17) with respect to the marginalization-consistency constraint in (15), which, however, involves a high-dimensional integral. To avoid the high-dimensional integral, we relax the marginalization-consistency constraint by imposing a moment-matching constraint as

$$\mathbb{E}[v|b_v(v)] = \mathbb{E}\left[v \mid \int b_f(\mathbf{u})d\mathbf{u}\setminus\mathbf{v}\right], \qquad (18a)$$

$$\mathbb{E}[vv^H|b_v(v)] = \mathbb{E}\left[vv^H \mid \int b_f(\mathbf{u})d\mathbf{u}\setminus\mathbf{v}\right], \qquad (18b)$$

which is typically used in EP [37]. We then adopt a factorization relaxation for $b_{f_Y}(\mathbf{R}, \mathbf{G}^{\mathrm{RB}})$:

$$b_{f_Y}(\mathbf{R}, \mathbf{G}^{\mathrm{RB}}) = b_{f_Y,R}(\mathbf{R})b_{f_Y,G}(\mathbf{G}^{\mathrm{RB}}), \qquad (19)$$

as commonly used in VMP [36], where $b_{f_Y,R}(\mathbf{R})$ and $b_{f_Y,G}(\mathbf{G}^{\mathrm{RB}})$ denote two factors for the factorization of $b_{f_Y}(\mathbf{R}, \mathbf{G}^{\mathrm{RB}})$. Then, the stationary points of (17) s.t. (18)-(19) are given by the following theorem:

**Theorem 1: Consider the variational optimization problem**

$$\min_{b(\cdot)} F(\{b_v(\mathbf{v})\}, \{b_f(\mathbf{u})\}) \quad \text{s.t.} \quad (18)-(19). \qquad (20)$$

The stationary points fulfill the fixed-point equations as

$$b_f(\mathbf{u}) \propto f(\mathbf{u}) \prod_{v \in \mathcal{N}(f)} \mathcal{M}_{v \to f}(v), \quad \forall f \in \mathcal{F} \setminus f_Y, \qquad (21a)$$

$$b_{f_Y, v}(\mathbf{v}) \propto e^{\left( \int \prod_{\tilde{v} \in \mathcal{N}(f_Y) \setminus v} b_{f_Y, \tilde{v}}(\tilde{v}) \ln f_Y(\mathbf{u}) d\mathbf{u} \right)} \times \mathcal{M}_{v \to f_Y}(\mathbf{v}), \mathbf{v} \in \{\mathbf{R}, \mathbf{G}^{\mathrm{RB}}\}, \qquad (21b)$$

$$b_v(\mathbf{v}) \propto \prod_{f \in \mathcal{N}(v)} \mathcal{M}_{f \to v}(\mathbf{v}), \quad \forall \mathbf{v} \in \mathcal{V}, \qquad (21c)$$

with

$$\mathcal{M}_{v \to f}(\mathbf{v}) \propto \prod_{\tilde{f} \in \mathcal{N}(v) \setminus f} \mathcal{M}_{\tilde{f} \to v}(\mathbf{v}), \qquad (22a)$$

$$\mathcal{M}_{f \to v}(\mathbf{v}) \propto \frac{\mathcal{P}\left(\int \prod_{\tilde{v} \in \mathcal{N}(f)} \mathcal{M}_{\tilde{v} \to f}(\tilde{v}) f(\mathbf{u}) d\mathbf{u} \setminus \mathbf{v}\right)}{\mathcal{M}_{v \to f}(\mathbf{v})}, \quad f \in \mathcal{F} \setminus f_Y, \qquad (22b)$$

$$\mathcal{M}_{f_Y \to v}(\mathbf{v}) \propto \frac{\mathcal{P}\left(e^{\int \prod_{\tilde{v} \in \mathcal{N}(f_Y) \setminus v} b_{f_Y, \tilde{v}}(\tilde{v}) \ln f_Y(\mathbf{u}) d\mathbf{u} \setminus \mathbf{v}} \mathcal{M}_{v \to f_Y}(\mathbf{v})\right)}{\mathcal{M}_{v \to f_Y}(\mathbf{v})}, \qquad (22c)$$

where $\mathcal{N}(f)$ is the set of variables associated with the factor node $f$; $\mathcal{N}(v)$ is the set of factor nodes associated with the variable node $v$; $\mathcal{P}(b(\cdot))$ denotes the projection of the pdf $b(\cdot)$ to a Gaussian pdf whose mean and precision match those of $b(\cdot)$.

**Proof:** The main idea of the proof follows [25, Appendix A]. For completeness, we present the proof in Appendix A. $\square$

In Theorem 1, $\mathcal{M}_{v \to f}(\mathbf{v})$ (or $\mathcal{M}_{f \to v}(\mathbf{v})$) is termed as the message from the variable node $v$ to the factor node $f$ (or from the factor node $f$ to the variable node $v$). To solve the stationary point equations in Theorem 1, we adopt an iterative method, i.e., iteratively compute the messages in (22a)-(22c). By substituting the messages into (21c), we obtain the estimates of $v$.

## B. Message-Passing Derivations

It is noteworthy that each variable node in Fig. 3 is connected to two factor nodes. Then, the message-passing in (22a) is simplified to $\mathcal{M}_{v \to f}(\mathbf{v}) = \mathcal{M}_{\tilde{f} \to v}(\mathbf{v})$. Therefore, it suffices to update $\mathcal{M}_{f \to v}(\mathbf{v})$, $b_f(\mathbf{u})$, and $b_v(\mathbf{v})$. Moreover, we emphasize that the structure of the factor graph is symmetrical for $(\mathbf{R}, \mathbf{H}^{\mathrm{UR}})$ and $(\mathbf{G}^{\mathrm{RB}}, \mathbf{H}^{\mathrm{RB}})$, which means that the message-passing equations for $(\mathbf{R}, \mathbf{H}^{\mathrm{UR}})$ and $(\mathbf{G}^{\mathrm{RB}}, \mathbf{H}^{\mathrm{RB}})$ share similar forms. In what follows, we derive the message-passing equations for $(\mathbf{R}, \mathbf{H}^{\mathrm{UR}})$ in detail and those for $(\mathbf{G}^{\mathrm{RB}}, \mathbf{H}^{\mathrm{RB}})$ are described briefly. We assume that the large-scale fading coefficient of each user is known to the BS. Then, we extract the large-scale fading coefficient and combine it with the pilot symbols, resulting in $\mathbb{E}[||\mathbf{G}^{\mathrm{RB}}(j)||_F^2]/NM = 1$ and $\mathbb{E}[||\mathbf{G}^{\mathrm{UR}}(j)||_F^2]/KN = 1$.
In practice, the large-scale fading coefficient of each user can be estimated from periodic signaling [32].

1) Message From the Factor Node $f_Y$ to the Variable Node $\mathbf{R}$: To start with, we derive the message $\mathcal{M}_{f_Y \to \mathbf{R}}(\mathbf{R})$. From (22c), we obtain

$$ \begin{aligned} \mathcal{M}_{f_Y \to \mathbf{R}}(\mathbf{R}) & \propto \frac{\mathcal{P}\left(e^{\int b_{f_Y,G}(\mathbf{G}^{\mathrm{RB}}) \ln f_Y(\mathbf{G}^{\mathrm{RB}}, \mathbf{R})} d\mathbf{G}^{\mathrm{RB}}\right) \mathcal{M}_{\mathbf{R} \to f_Y(\mathbf{R})}}{\mathcal{M}_{\mathbf{R} \to f_Y(\mathbf{R})}} \\ & \propto e^{\int b_{f_Y,G}(\mathbf{G}^{\mathrm{RB}}) \ln f_Y(\mathbf{G}^{\mathrm{RB}}, \mathbf{R})} d\mathbf{G}^{\mathrm{RB}}. \end{aligned} \quad (23) $$

(23) is obtained by noting that $f_Y(\mathbf{R}, \mathbf{G}^{\mathrm{RB}})$ is Gaussian from (10a), which makes $e^{\int b_{f_Y,G}(\mathbf{G}^{\mathrm{RB}}) \ln f_Y(\mathbf{G}^{\mathrm{RB}}, \mathbf{R})} d\mathbf{G}^{\mathrm{RB}}$ proportional to a Gaussian pdf of $\mathbf{R}$ (see Appendix B). Then, the Gaussian projector $\mathcal{P}(\cdot)$ can be removed, obtaining (23). Since $f_Y(\mathbf{R}, \mathbf{G}^{\mathrm{RB}})$ is Gaussian, we further express $b_{f_Y,G}(\mathbf{G}^{\mathrm{RB}})$ in (23) as

$$ b_{f_Y,G}(\mathbf{G}^{\mathrm{RB}}) \propto e^{\int b_{f_Y,G} (\mathbf{G}^{\mathrm{RB}}) \ln f_Y(\mathbf{G}^{\mathrm{RB}}, \mathbf{R}) d\mathbf{R}} \mathcal{M}_{\mathbf{G}^{\mathrm{RB}} \to f_Y(\mathbf{G}^{\mathrm{RB}})}, \quad (24a) $$

$$ \propto \mathcal{M}_{f_Y \to \mathbf{G}^{\mathrm{RB}}}(\mathbf{G}^{\mathrm{RB}}) \mathcal{M}_{\mathbf{G}^{\mathrm{RB}} \to f_Y}(\mathbf{G}^{\mathrm{RB}}), \quad (24b) $$

$$ \propto b_{\mathbf{G}^{\mathrm{RB}}}(\mathbf{G}^{\mathrm{RB}}), \quad (24c) $$

where (24a) follows from (21b), (24b) from (22b), and (24c) from combining (21c) and (22a). For the belief $b_{\mathbf{G}^{\mathrm{RB}}}(\mathbf{G}^{\mathrm{RB}})$ in (24c), we initialize $b_{\mathbf{G}^{\mathrm{RB}}}(\mathbf{G}^{\mathrm{RB}}) = \mathcal{CN}(\mathbf{G}^{\mathrm{RB}}; \hat{\mathbf{G}}^{\mathrm{RB}}, \Sigma_{\mathbf{G}}^{\mathrm{RB}})$ with $\hat{\mathbf{G}}^{\mathrm{RB}} = \mathbf{G}^{\mathrm{RB}}$ and $\Sigma_{\mathbf{G}}^{\mathrm{RB}} = \mathbf{I}$ in the first iteration. The update of $b_{\mathbf{G}^{\mathrm{RB}}}(\mathbf{G}^{\mathrm{RB}})$ in the subsequent iteration is shown in (54). By substituting (24c) into (23) and applying some algebraic manipulations, we obtain

$$ \mathcal{M}_{f_Y \to \mathbf{R}}(\mathbf{R}) = \mathcal{CN}(\mathbf{R}; \mathbf{R}_{f_Y \to \mathbf{R}}, \Sigma_{f_Y \to \mathbf{R}}), \quad (25) $$

where

$$ \mathbf{R}_{f_Y \to \mathbf{R}} = \beta_1^{-1} \Sigma_{f_Y \to \mathbf{R}}^{-1} \mathrm{diag}(\hat{\mathbf{G}}^{\mathrm{RB}}(j)^H \mathbf{Y}), \quad (26a) $$

$$ \Sigma_{f_Y \to \mathbf{R}} = \beta_1^{-1} \mathrm{diag} \left( (\hat{\mathbf{G}}^{\mathrm{RB}}(j)^H \hat{\mathbf{G}}^{\mathrm{RB}}(j) + \mathrm{cov}_{b_{\mathbf{G}^{\mathrm{RB}}}[\mathbf{G}^{\mathrm{RB}}(j)]}) \right), \quad (26b) $$

with $\beta_1 = \sigma_w^2$. The derivation of (25)-(26) is provided in Appendix B. Note that $\mathbf{R}_{f_Y \to \mathbf{R}}$ and $\mathrm{tr}(\Sigma_{f_Y \to \mathbf{R}})/JN$ are treated as the mean and variance of $\mathbf{R}$, respectively. Given $\mathbf{R}_{f_Y \to \mathbf{R}}$ as the mean of $\mathbf{R}$, we find from simulations that $\mathrm{tr}(\Sigma_{f_Y \to \mathbf{R}})/JN$ significantly deviates from $E[\|\mathbf{R} - \mathbf{R}_{f_Y \to \mathbf{R}}\|_F^2]/(JNT)$. This is because, with the variational approximation in (19), $\mathrm{tr}(\Sigma_{f_Y \to \mathbf{R}})/JN$ is a biased estimator of the variance [38]. Then, the algorithm may converge to a local minimum. To avoid this issue, many methods have been introduced such as increasing the initial variances of the messages and re-starting the algorithm several times [39]. Here, instead of using these ad-hoc methods, we propose a method based on a model approximation. Specifically, we rearrange the bilinear model in (10a) as a linear model of $\mathbf{R}(j)$, i.e.,

$$ \mathbf{Y}(j) = \hat{\mathbf{G}}^{\mathrm{RB}}(j)\mathbf{R}(j) + \mathbf{Z}(j), \quad \forall j, \quad (27) $$

<figure><img src="image_5.png" /><figcaption>Fig. 4. A factorization of variable node $\mathbf{H}^{\mathrm{UR}}$ and factor node $f_{\mathbf{H}}^{\mathrm{UR}}$.</figcaption></figure>

where $\mathbf{Z}(j) = (\mathbf{G}^{\mathrm{RB}}(j) - \hat{\mathbf{G}}^{\mathrm{RB}}(j))\mathbf{R}(j) + \mathbf{W}(j)$ with average power $\beta_1$. Given the message $b_{\mathbf{G}^{\mathrm{RB}}}(\mathbf{G}^{\mathrm{RB}}) = \mathcal{CN}(\mathbf{G}^{\mathrm{RB}}; \hat{\mathbf{G}}^{\mathrm{RB}}, \Sigma_{\mathbf{G}}^{\mathrm{RB}})$ with $\mathbf{G}^{\mathrm{RB}} = [\mathbf{G}^{\mathrm{RB}}(1), \dots, \mathbf{G}^{\mathrm{RB}}(J)]^T$ and the covariance $(\Sigma_{\mathbf{G}}^{\mathrm{RB}})^{-1} = \bar{v}_h^{\mathrm{RB}} \mathbf{F}_2 \mathbf{F}_2^H = \bar{v}_h^{\mathrm{RB}} (\mathbf{D}_{L_2} \mathbf{D}_{L_2}^H) \otimes \mathbf{I}_N$, we obtain $\mathbf{G}^{\mathrm{RB}}(j) - \hat{\mathbf{G}}^{\mathrm{RB}}(j)$ on subcarrier $j$ as a Gaussian random matrix with uncorrelated entries. Although the entries of $\mathbf{G}^{\mathrm{RB}}(j) - \hat{\mathbf{G}}^{\mathrm{RB}}(j)$ on different subcarriers are correlated, this correlation is ignored as we consider a subcarrier-wise estimation of $\mathbf{R}(j), \forall j$, from (27). Thus, $\mathbf{Z}(j)$ is treated as an AWGN, and the variance is $\beta_1 = P_x K \mathrm{tr}((\Sigma_{\mathbf{G}}^{\mathrm{RB}})^{-1})/J + \sigma_w^2$. Then, (27) is an AWGN model of $\mathbf{R}(j)$ with a linear mixing matrix $\hat{\mathbf{G}}^{\mathrm{RB}}(j)$, which corresponds to the pdf

$$ \begin{aligned} & e^{\frac{1}{\beta_1} \mathrm{Tr}((\mathbf{Y}(j) - \hat{\mathbf{G}}^{\mathrm{RB}}(j)\mathbf{R}(j))^H ((\mathbf{Y}(j) - \hat{\mathbf{G}}^{\mathrm{RB}}(j))\mathbf{R}(j)))}, \\ & \propto e^{\frac{1}{\beta_1} \mathrm{Tr}(\mathbf{R}(j)^H \hat{\mathbf{G}}^{\mathrm{RB}}(j)^H \hat{\mathbf{G}}^{\mathrm{RB}}(j)\mathbf{R}(j) - 2\Re(\mathbf{R}(j)^H \hat{\mathbf{G}}^{\mathrm{RB}}(j)^H \mathbf{Y}(j)))}, \\ & \propto \mathcal{CN}(\mathbf{R}(j); \mathbf{R}_{f_Y(j) \to \mathbf{R}(j)}, \Sigma_{f_Y(j) \to \mathbf{R}(j)}), \\ & \triangleq \mathcal{M}_{f_Y(j) \to \mathbf{R}(j)}(\mathbf{R}(j)), \end{aligned} \quad (28) $$

where $\Sigma_{f_Y(j) \to \mathbf{R}(j)} = \beta_1^{-1} \hat{\mathbf{G}}^{\mathrm{RB}}(j)^H \hat{\mathbf{G}}^{\mathrm{RB}}(j)$ and $\mathbf{R}_{f_Y(j) \to \mathbf{R}(j)} = \beta_1^{-1} \Sigma_{f_Y(j) \to \mathbf{R}(j)}^{-1} \hat{\mathbf{G}}^{\mathrm{RB}}(j)^H \mathbf{Y}(j)$. The compact form on $J$ subcarriers is expressed as $\mathcal{M}_{f_Y \to \mathbf{R}}(\mathbf{R}) = \mathcal{CN}(\mathbf{R}; \mathbf{R}_{f_Y \to \mathbf{R}}, \Sigma_{f_Y \to \mathbf{R}})$ with

$$ \Sigma_{f_Y \to \mathbf{R}} \mathbf{R}_{f_Y \to \mathbf{R}} = \beta_1^{-1} \mathrm{diag} \left( \hat{\mathbf{G}}^{\mathrm{RB}}(j)^H \mathbf{Y} \right), \quad (29a) $$

$$ \Sigma_{f_Y \to \mathbf{R}} = \beta_1^{-1} \mathrm{diag} \left( \hat{\mathbf{G}}^{\mathrm{RB}}(j)^H \hat{\mathbf{G}}^{\mathrm{RB}}(j) \right), \quad (29b) $$

In (29a), $\Sigma_{f_Y \to \mathbf{R}} \mathbf{R}_{f_Y \to \mathbf{R}}$ (instead of $\mathbf{R}_{f_Y \to \mathbf{R}}$) is given to avoid the matrix inversion $\Sigma_{f_Y \to \mathbf{R}}^{-1}$, and we show later in (32a) that we only require $\Sigma_{f_Y \to \mathbf{R}} \mathbf{R}_{f_Y \to \mathbf{R}}$ for message-passing.

2) Message From the Factor Node $f_R$ to the Variable Node $\mathbf{H}^{\mathrm{UR}}$: We obtain $\mathcal{M}_{f_Y \to \mathbf{R}}(\mathbf{R}) = \mathcal{M}_{\mathbf{R} \to f_R}(\mathbf{R})$ from (22a). We then compute the message $\mathcal{M}_{f_R \to \mathbf{H}^{\mathrm{UR}}}(\mathbf{H}^{\mathrm{UR}})$. However, combining $\mathcal{M}_{f_R \to \mathbf{H}^{\mathrm{UR}}}(\mathbf{H}^{\mathrm{UR}})$ and $p(\mathbf{H}^{\mathrm{UR}})$ involves $2^{MN L_2}$ terms due to the Bernoulli (0/1) component of $p(\mathbf{H}^{\mathrm{UR}})$ in (12). To reduce the computational complexity, we construct a further factorization for the variable node $\mathbf{H}^{\mathrm{UR}}$ and the factor node $p(\mathbf{H}^{\mathrm{UR}})$, as illustrated in Fig. 4. Specifically, the factor node $f_{\mathbf{H}}^{\mathrm{UR}}(\mathbf{H}^{\mathrm{UR}})$ in factorized to $NKL_1$ factor nodes $\{f_{h_{k,n}^{\mathrm{UR}}(l)}\}$ based on $p(\mathbf{H}^{\mathrm{UR}}) = \prod_{k,n,l} p(h_{h_{k,n}^{\mathrm{UR}}(l)})$, and the variable node $\mathbf{H}^{\mathrm{UR}}$ is factorized to $NKL_1$ variable nodes $\{h_{h_{k,n}^{\mathrm{UR}}(l)}\}$. Correspondingly, the factor node $f_{\mathbf{H}}^{\mathrm{UR}}(\mathbf{H}^{\mathrm{UR}})$ and the variable node $\mathbf{H}^{\mathrm{UR}}$ in Fig. 3 need to be replaced with $\{f_{h_{h_{k,n}^{\mathrm{UR}}(l)}}\}$ and $\{h_{h_{k,n}^{\mathrm{UR}}(l)}\}$ in Fig. 4. With this replacement,
we use (22b) to compute

$$
\begin{align}
& \mathcal{M}_{f_R \to h_{k,n}^{\text{UR}}(l)}(h_{k,n}^{\text{UR}}(l)) \nonumber \\
& = \left( \iint p(\mathbf{R}|\mathbf{H}^{\text{UR}}) \mathcal{M}_{\mathbf{R} \to f_R(\mathbf{R})} \prod_k \mathcal{M}_{f_{h_k^{\text{UR}} \to h_k^{\text{UR}}}(\mathbf{h}_k^{\text{UR}})} d\mathbf{R} \right. \nonumber \\
& \qquad \left. d\mathbf{H}^{\text{UR}} \setminus h_{k,n}^{\text{UR}}(l) \right) / \mathcal{M}_{f_{h_{k,n}^{\text{UR}}(l)} \to h_{k,n}^{\text{UR}}(l)}(h_{k,n}^{\text{UR}}(l)), \tag{30}
\end{align}
$$

with $\mathcal{M}_{f_{h_k^{\text{UR}} \to h_k^{\text{UR}}}(\mathbf{h}_k^{\text{UR}})} \triangleq \prod_{k,n,l} \mathcal{M}_{f_{h_{k,n}^{\text{UR}}}(l) \to h_{k,n}^{\text{UR}}(l)}$
$(h_{k,n}^{\text{UR}}(l)) = \mathcal{CN}(\mathbf{h}_k^{\text{UR}}; \mathbf{h}_{f_{h_k^{\text{UR}} \to h_k^{\text{UR}}}, \Sigma_{f_{h_k^{\text{UR}} \to h_k^{\text{UR}}}}(l))$, where
$\mathbf{h}_{f_{h_k^{\text{UR}} \to h_k^{\text{UR}}}} = [h_{f_{h_{k,1}^{\text{UR}}(1)} \to h_{k,1}^{\text{UR}}(1)}, \dots, h_{f_{h_{k,N}^{\text{UR}}(L_1)} \to h_{k,N}^{\text{UR}}(L_1)}]^T$
and $\Sigma_{f_{h_k^{\text{UR}} \to h_k^{\text{UR}}}}$ is a diagonal matrix with diagonal elements
$[v_{f_{h_{k,1}^{\text{UR}}(1)} \to h_{k,1}^{\text{UR}}(1)}^{\text{U}} \dots v_{f_{h_{k,N}^{\text{UR}}(L_1)} \to h_{k,N}^{\text{UR}}(L_1)}]$. We initialize
$h_{f_{h_{k,n}^{\text{UR}}(l)} \to h_{k,n}^{\text{UR}}(l)} = 0$ and $v_{f_{h_{k,n}^{\text{UR}}(l)} \to h_{k,n}^{\text{UR}}(l)} = 1$, and the
following update is given by (41). Based on the orthogonal
pilot assignment, the integral with respect to $\mathbf{R}$ in (30) can
be expressed as a message product of $\mathbf{h}_k^{\text{UR}}$, $k = 1, \dots, K$,
i.e.,

$$
\int_{\mathbf{R}} p(\mathbf{R} | \mathbf{H}^{\text{UR}}) \mathcal{M}_{\mathbf{R} \to f_{\mathbf{R}}(\mathbf{R})} \prod_k \mathcal{M}_{f_{h_k^{\text{UR}} \to h_k^{\text{UR}}}(\mathbf{h}_k^{\text{UR}})} \\
= \prod_{k=1}^{K} \mathcal{CN}(\mathbf{h}_k^{\text{UR}}; \mathbf{h}_{f_{\mathbf{R}}, \mathbf{h}_k^{\text{UR}}, \Sigma_{f_{\mathbf{R}}, k}}), \quad (31)
$$

where

$$
\mathbf{h}_{f_R, \mathbf{h}_k^{\text{UR}}} = \mathbf{\Sigma}_{f_R, k}^{-1} \left( \sum_{t=1}^{T} \mathbf{S}_{k,t}^H \mathbf{\Sigma}_{f_Y \to \mathbf{R} \mathbf{r}_{f_Y \to \mathbf{r}_t}} \right. \\
\left. + \mathbf{\Sigma}_{f_{\mathbf{h}_k^{\text{UR}}} \to \mathbf{h}_k^{\text{UR}}} \mathbf{h}_{f_{\mathbf{h}_k^{\text{UR}}} \to \mathbf{h}_k^{\text{UR}}} \right), \quad (32a)
$$

$$
\Sigma_{f_R, k} = \sum_{t=1}^{T} S_{k,t}^H \Sigma_{f_Y \to R} S_{k,t} + \Sigma_{f_k^{\text{UR}} \to h_k^{\text{UR}}}, \quad (32b)
$$

$$
S_{k,t} = \mathrm{diag}(\mathbf{x}_{k,t} \otimes \boldsymbol{\theta}_t) \mathbf{F}_1 \mathbf{U}_R. \quad (32c)
$$

The derivation of (31)-(32) is given in Appendix C. Note that without the pilot orthogonality, (31) becomes a joint message of $\mathbf{h}_k^{\text{UR}}$, $k = 1, \dots, K$, and the corresponding computational complexity increases due to performing higher-dimensional matrix computations. Let the diagonal elements of $\mathbf{\Sigma}_{f_R,k}^{-1}$ be $[\sigma_{f_R,h_{k,1}^{\text{UR}}(1)}, \dots, \sigma_{f_R,h_{k,N}^{\text{UR}}(L_1)}]$. Substituting (31) into (30) yields

$$
\begin{align}
& \mathcal{M}_{f_R \to h_{k,n}^{\text{UR}}(l)}(h_{k,n}^{\text{UR}}(l)) \nonumber \\
& = \mathcal{CN}(h_{k,n}^{\text{UR}}(l); h_{f_R \to h_{k,n}^{\text{UR}}(l)}, v_{f_R \to h_{k,n}^{\text{UR}}(l)}), \tag{33}
\end{align}
$$

with

$$
h_{f_R \to h_{k,n}^{\text{UR}}(l)} = v_{f_R \to h_{k,n}^{\text{UR}}(l)} \left( \frac{h_{f_R, h_{k,n}^{\text{UR}}(l)}}{\sigma_{f_R, h_{k,n}^{\text{UR}}(l)}^{-1}} - \frac{h_{f_{h_{k,n}^{\text{UR}}(l)} \to h_{k,n}^{\text{UR}}(l)}}{v_{f_{h_{k,n}^{\text{UR}}(l)} \to h_{k,n}^{\text{UR}}(l)}} \right), \quad (34a)
$$

$$
v_{f_R \to h_{k,n}^{\text{UR}}(l)} = \left( \sigma_{f_R, h_{k,n}^{\text{UR}}(l)} - v_{f_{h_{k,n}^{\text{UR}}(l)} \to h_{k,n}^{\text{UR}}(l)}^{-1} \right)^{-1}, \quad (34b)
$$

by using the Gaussian combining property.³

3) Message From the Factor Node $f_{\mathbf{H}}^{\text{UR}}$ to the Variable Node $\mathbf{H}^{\text{UR}}$ and Beliefs at the Variable Nodes $\mathbf{H}^{\text{UR}}$ and $\mathbf{R}$: Given $\mathcal{M}_{f_R \to h_{k,n}^{\text{UR}}(l)}(h_{k,n}^{\text{UR}}(l))$, we apply (21a) to compute

$$
b_{f_{h_{k,n}^{\text{UR}}(l)}(h_{k,n}^{\text{UR}}(l))} = p(h_{k,n}^{\text{UR}}(l)) \mathcal{M}_{f_R \to h_{k,n}^{\text{UR}}(l)}(h_{k,n}^{\text{UR}}(l)). \quad (35)
$$

By noting that the product of a Gaussian pdf and a Bernoulli Gaussian pdf gives a Bernoulli Gaussian pdf, we obtain

$$
b_{f_{h_{k,n}^{\text{UR}}}(l)}(h_{k,n}^{\text{UR}}(l)) = (1 - \hat{\lambda}_{k,n}^{\text{UR}}(l))\delta(h_{k,n}^{\text{UR}}(l)) \\
+ \hat{\lambda}_{k,n}^{\text{UR}}(l)\mathcal{CN}(h_{k,n}^{\text{UR}}(l); \mu_{k,n}^{\text{UR}}(l); \phi_{k,n}^{\text{UR}}(l)), \tag{36}
$$

where

$$
\hat{\lambda}_{k,n}^{\text{UR}}(l) = \left( 1 + \frac{(1 - \lambda_1) \mathcal{CN}(h_{f_R \to h_{k,n}^{\text{UR}}(l)}, v_{f_R \to h_{k,n}^{\text{UR}}(l)})}{\lambda_1 \mathcal{CN}(h_{f_R \to h_{k,n}^{\text{UR}}(l)}, v_{f_R \to h_{k,n}^{\text{UR}}(l)} + \vartheta^{\text{UR}}(l))} \right)^{-1} \quad (37a)
$$

$$
\phi_{k,n}^{\text{UR}}(l) = (\vartheta^{\text{UR}}(l)^{-1} + v_{f_R \to h_{k,n}^{\text{UR}}(l)}^{-1})^{-1}, \quad (37b)
$$

$$
\mu_{k,n}^{\text{UR}}(l) = \phi_{k,n}^{\text{UR}}(l) v_{f_R \to h_{k,n}^{\text{UR}}(l)}^{-1} h_{f_R \to h_{k,n}^{\text{UR}}(l)}. \quad (37c)
$$

With $b_{f_{h_{k,n}^{\text{UR}}(l)}(h_{k,n}^{\text{UR}}(l))}$, we apply (22b) to compute

$$
\mathcal{M}_{f_{h_{k,n}^{\text{UR}}(l)} \to h_{k,n}^{\text{UR}}(l)}(h_{k,n}^{\text{UR}}(l)) = \frac{\mathcal{P}\left(b_{f_{h_{k,n}^{\text{UR}}(l)}}(h_{k,n}^{\text{UR}}(l))\right)}{\mathcal{M}_{f_{h_{k,n}^{\text{UR}}(l)}}(h_{k,n}^{\text{UR}}(l))}. \quad (38)
$$

Note that $\mathcal{P}(b_{f_{h_{k,n}^{\text{UR}}(l)}}(h_{k,n}^{\text{UR}}(l)))$ in the numerator of the
right hand-side of (38) is exactly $b_{h_{k,n}^{\text{UR}}(l)}(h_{k,n}^{\text{UR}}(l))$ according
to (21c). The corresponding mean and variance are respec-
tively expressed as

$$
\hat{h}_{k,n}^{\text{UR}}(l) = \hat{\lambda}_{k,n}^{\text{UR}}(l) \mu_{k,n}^{\text{UR}}(l), \quad (39a)
$$

$$
v_{h_{k,n}^{\text{UR}}(l)} = \hat{\lambda}_{k,n}^{\text{UR}}(l) (|\mu_{k,n}^{\text{UR}}(l)|^2 + \phi_{k,n}^{\text{UR}}(l)) - |\hat{h}_{k,n}^{\text{UR}}(l)|^2. \quad (39b)
$$

We approximate $v_{h_{k,n}^{\text{UR}}(l)}, \forall k, n, l$ by their average $\bar{v}_h^{\text{UR}} = \sum_{k,n,l} v_{h_{k,n}^{\text{UR}}(l)}/KN\bar{L}_1$, which is a widely used simplification in message-passing algorithms [40], [41]. Substituting (39) into (38) yields

$$
\begin{align}
& \mathcal{M}_{f_{h_{k,n}^{\text{UR}}(l)} \to h_{k,n}^{\text{UR}}(l)}(h_{k,n}^{\text{UR}}(l)) \nonumber \\
& = \mathcal{CN}(h_{k,n}^{\text{UR}}(l); h_{f_{h_{k,n}^{\text{UR}}(l)} \to h_{k,n}^{\text{UR}}(l)}, v_{f_{h_{k,n}^{\text{UR}}(l)} \to h_{k,n}^{\text{UR}}(l)}), \tag{40}
\end{align}
$$

with

$$
h_{f_{h_{k,n}^{\text{UR}}(l)} \to h_{k,n}^{\text{UR}}(l)} = v_{f_{h_{k,n}^{\text{UR}}(l)} \to h_{k,n}^{\text{UR}}(l)} \left( \frac{\hat{h}_{k,n}^{\text{UR}}(l)}{\bar{v}_h^{\text{UR}}} - \frac{h_{f_R \to h_{k,n}^{\text{UR}}(l)}}{v_{f_R \to h_{k,n}^{\text{UR}}(l)}} \right), \quad (41a)
$$

$$
v_{f_{h_{k,n}^{\text{UR}}(l)} \to h_{k,n}^{\text{UR}}(l)} = \left( (\bar{v}_h^{\text{UR}})^{-1} - v_{f_R \to h_{k,n}^{\text{UR}}(l)}^{-1} \right)^{-1} . \quad (41b)
$$

$^3\mathcal{C}\mathcal{N}(x; a, A)\mathcal{C}\mathcal{N}(x; b, B) = d\mathcal{C}\mathcal{N}(x; c, C)$ with $d = \mathcal{C}\mathcal{N}(0; a - b, A + B)$, $c = C(A^{-1}a + B^{-1}b)$, and $C = (A^{-1} + B^{-1})^{-1}$.
Based on $b_{h_{k,n}^{\text{UR}}}(h_{k,n}^{\text{UR}}(l))$ in (36) and $\mathbf{r}_t = \sum_{k=1}^{K} \mathbf{S}_{k,t} \mathbf{h}_k^{\text{UR}}, \forall t$ in (10b), we obtain $b_{\mathbf{R}}(\mathbf{R}) = \prod_{t=1}^{T} b_{\mathbf{r}_t}(\mathbf{r}_t) = \prod_{t=1}^{T} \mathcal{C}\mathcal{N}(\mathbf{r}_t; \hat{\mathbf{r}}_t, \mathbf{\Sigma}_{\mathbf{r}_t})$ with

$$ \hat{\mathbf{r}}_t = \sum_{k=1}^{K} \mathbf{S}_{k,t} \hat{\mathbf{h}}_k^{\text{UR}}, \quad (42a) $$

$$ \mathbf{\Sigma}_{\mathbf{r}_t}^{-1} = \sum_{k=1}^{K} \bar{v}_h^{\text{UR}} \mathbf{S}_{k,t} \mathbf{S}_{k,t}^{\text{H}}, \quad (42b) $$

where $\hat{\mathbf{h}}_k^{\text{UR}} = [\hat{h}_{1,1}^{\text{UR}}(1), \dots, \hat{h}_{K,N}^{\text{UR}}(L_1)]^T$.

4) Message From the Factor Node $f_Y$ to the Variable Node $\mathbf{G}^{\text{RB}}$: The derivation of $\mathcal{M}_{f_Y \to \mathbf{G}^{\text{RB}}}(\mathbf{G}^{\text{RB}})$ is similar to that of $\mathcal{M}_{f_Y \to \mathbf{R}}(\mathbf{R})$ in Sec. IV-B.1. We obtain the message $\mathcal{M}_{f_Y \to \mathbf{G}^{\text{RB}}}(\mathbf{G}^{\text{RB}})$ as

$$ \mathcal{M}_{f_Y \to \mathbf{G}^{\text{RB}}}(\mathbf{G}^{\text{RB}}) = \mathcal{C}\mathcal{N}(\mathbf{G}^{\text{RB}}; \mathbf{G}_{f_Y \to \mathbf{G}^{\text{RB}}}, \mathbf{\Sigma}_{f_Y \to \mathbf{G}^{\text{RB}}}), \quad (43) $$

with

$$ \mathbf{\Sigma}_{f_Y \to \mathbf{G}^{\text{RB}}} \mathbf{G}_{f_Y \to \mathbf{G}^{\text{RB}}} = \beta_2^{-1} \text{diag}(\hat{\mathbf{R}}(j)^*) \tilde{\mathbf{Y}}, \quad (44a) $$

$$ \mathbf{\Sigma}_{f_Y \to \mathbf{G}^{\text{RB}}} = \beta_2^{-1} \text{diag} \left( \hat{\mathbf{R}}^*(j) \hat{\mathbf{R}}(j)^T \right), \quad (44b) $$

where $\tilde{\mathbf{Y}} = [\mathbf{Y}(1), \dots, \mathbf{Y}(J)]^T$, $[\hat{\mathbf{R}}(1)^T, \dots, \hat{\mathbf{R}}(J)^T]^T = [\hat{\mathbf{r}}_1, \dots, \hat{\mathbf{r}}_T]$, and $\beta_2 = \sum_t \text{tr}(\mathbf{\Sigma}_{\mathbf{r}_t}^{-1})/(JT) + \sigma_w^2$.

5) Message From the Factor Node $f_{\mathbf{G}}^{\text{RB}}$ to the Variable Node $\mathbf{H}^{\text{RB}}$: Similar to the approach in Sec. IV-B.2, we factorize $f_{\mathbf{H}}^{\text{RB}}(\mathbf{H}^{\text{RB}})$ into the $\text{NML}_2$ factor nodes $\{f_{h_{n,m}^{\text{RB}}}(l)\}$ based on $p(\mathbf{H}^{\text{RB}}) = \prod_{n,m,l} p(h_{n,m}^{\text{RB}}(l))$, and factorize $\mathbf{H}^{\text{RB}}$ into the $\text{NML}_2$ variable nodes $\{h_{n,m}^{\text{RB}}(l)\}$. Similarly to (30)-(32), we obtain the quantity

$$ \begin{aligned} & q_{f_{\mathbf{G}}^{\text{RB}}}(\mathbf{H}^{\text{RB}}) \\ & \triangleq \int_{\mathbf{G}^{\text{RB}}} p(\mathbf{G}^{\text{RB}}|\mathbf{H}^{\text{RB}}) \mathcal{M}_{f_Y \to \mathbf{G}^{\text{RB}}}(\mathbf{G}^{\text{RB}}) \mathcal{M}_{f_{\mathbf{H}}^{\text{RB}} \to \mathbf{H}^{\text{RB}}}(\mathbf{H}^{\text{RB}}) \\ & = \mathcal{C}\mathcal{N}(\mathbf{H}^{\text{RB}}; \mathbf{H}_{f_{\mathbf{G}}^{\text{RB}}}, \mathbf{\Sigma}_{f_{\mathbf{G}}^{\text{RB}}}), \end{aligned} \quad (45) $$

where

$$ \mathbf{H}_{f_{\mathbf{G}}^{\text{RB}}} = \mathbf{\Sigma}_{f_{\mathbf{G}}^{\text{RB}}}^{-1} \left( \mathbf{\Sigma}_{f_{\mathbf{G}}^{\text{RB}}} \to \mathbf{H}^{\text{RB}} \mathbf{H}_{f_{\mathbf{G}}^{\text{RB}}} \to \mathbf{H}^{\text{RB}} + \mathbf{V}_{f_{\mathbf{H}}^{\text{RB}}} \to \mathbf{H}^{\text{RB}} \mathbf{H}_{f_{\mathbf{H}}^{\text{RB}}} \to \mathbf{H}^{\text{RB}} \right), \quad (46a) $$

$$ \mathbf{\Sigma}_{f_{\mathbf{G}}^{\text{RB}}} = \mathbf{\Sigma}_{f_{\mathbf{G}}^{\text{RB}}} \to \mathbf{H}^{\text{RB}} + \mathbf{V}_{f_{\mathbf{H}}^{\text{RB}}} \to \mathbf{H}^{\text{RB}}, \quad (46b) $$

$$ \mathbf{\Sigma}_{f_{\mathbf{G}}^{\text{RB}} \to \mathbf{H}^{\text{RB}}} \mathbf{H}_{f_{\mathbf{G}}^{\text{RB}}} \to \mathbf{H}^{\text{RB}} = \mathbf{U}_{\mathbf{B}}^{\text{H}} \mathbf{F}_{2}^{\text{H}} \mathbf{\Sigma}_{f_{Y} \to \mathbf{G}^{\text{RB}}} \times (\mathbf{G}_{f_{Y} \to \mathbf{G}^{\text{RB}}} - \bar{\mathbf{G}}^{\text{RB}}) \mathbf{A}_{\mathbf{B}}^{\text{H}}, \quad (46c) $$

$$ \mathbf{\Sigma}_{f_{\mathbf{G}}^{\text{RB}} \to \mathbf{H}^{\text{RB}}} = \mathbf{U}_{\mathbf{B}}^{\text{H}} (\mathbf{F}_{2}^{\text{H}} \mathbf{\Sigma}_{f_{Y} \to \mathbf{G}^{\text{RB}}} \mathbf{F}_{2}) \mathbf{U}_{\mathbf{B}}. \quad (46d) $$

Let the diagonal elements of $\mathbf{\Sigma}_{f_{\mathbf{G}}^{\text{RB}}}^{-1}$ be $[\sigma_{f_{\mathbf{G}}^{\text{RB}}, h_{1,m}^{\text{RB}}(1)}, \dots, \sigma_{f_{\mathbf{G}}^{\text{RB}}, h_{n,m}^{\text{RB}}(L_2)}^{-1}]$, $\forall m$. We then apply (45) to compute

$$ \begin{aligned} & \mathcal{M}_{f_{\mathbf{G}}^{\text{RB}} \to h_{n,m}^{\text{RB}}(l)}(h_{n,m}^{\text{RB}}(l)) \\ & = \frac{\int q_{f_{\mathbf{G}}^{\text{RB}}}(\mathbf{H}^{\text{RB}}) d\mathbf{H}^{\text{RB}} \setminus h_{n,m}^{\text{RB}}(l)}{\mathcal{M}_{f_{h_{n,m}^{\text{RB}}(l)} \to h_{n,m}^{\text{RB}}(l)}(h_{n,m}^{\text{RB}}(l))}, \\ & = \mathcal{C}\mathcal{N}(h_{n,m}^{\text{RB}}(l); h_{f_{\mathbf{G}}^{\text{RB}} \to h_{n,m}^{\text{RB}}(l)}, v_{f_{\mathbf{G}}^{\text{RB}} \to h_{n,m}^{\text{RB}}(l)}), \end{aligned} \quad (47) $$

with

$$ \begin{aligned} & h_{f_{\mathbf{G}}^{\text{RB}} \to h_{n,m}^{\text{RB}}(l)} = v_{f_{\mathbf{G}}^{\text{RB}} \to h_{n,m}^{\text{RB}}(l)} (h_{f_{\mathbf{G}}^{\text{RB}}, h_{n,m}^{\text{RB}}(l)} / \sigma_{f_{\mathbf{G}}^{\text{RB}}, h_{n,m}^{\text{RB}}(l)}^{-1} \\ & \qquad - h_{f_{h_{n,m}^{\text{RB}}(l)} \to h_{n,m}^{\text{RB}}(l)} / v_{f_{h_{n,m}^{\text{RB}}(l)} \to h_{n,m}^{\text{RB}}(l)}), \end{aligned} \quad (48a) $$

$$ v_{f_{\mathbf{G}}^{\text{RB}} \to h_{n,m}^{\text{RB}}(l)} = \left( \sigma_{f_{\mathbf{G}}^{\text{RB}}, h_{n,m}^{\text{RB}}(l)} - v_{f_{h_{n,m}^{\text{RB}}(l)} \to h_{n,m}^{\text{RB}}(l)}^{-1} \right)^{-1}. \quad (48b) $$

6) Message From the Factor Node $f_{\mathbf{H}}^{\text{RB}}$ to the Variable Node $\mathbf{H}^{\text{RB}}$ and Beliefs at the Variable Nodes $\mathbf{H}^{\text{RB}}$ and $\mathbf{G}^{\text{RB}}$: Similar to (36), we obtain the belief $b_{f_{h_{n,m}^{\text{RB}}(l)}}(h_{n,m}^{\text{RB}}(l))$ as

$$ \begin{aligned} b_{f_{h_{n,m}^{\text{RB}}(l)}}(h_{n,m}^{\text{RB}}(l)) = & (1 - \hat{\lambda}_{n,m}^{\text{RB}}(l)) \delta(h_{n,m}^{\text{RB}}(l)) \\ & + \hat{\lambda}_{n,m}^{\text{RB}}(l) \mathcal{C}\mathcal{N}(h_{n,m}^{\text{RB}}(l); \mu_{n,m}^{\text{RB}}(l); \phi_{n,m}^{\text{RB}}(l)), \end{aligned} \quad (49) $$

where

$$ \begin{aligned} & \hat{\lambda}_{n,m}^{\text{RB}}(l) \\ & = \left( 1 + \frac{(1 - \lambda_2) \mathcal{C}\mathcal{N}(h_{f_{\mathbf{G}}^{\text{RB}} \to h_{n,m}^{\text{RB}}(l)}, v_{f_{\mathbf{G}}^{\text{RB}} \to h_{n,m}^{\text{RB}}(l)})}{\lambda_2 \mathcal{C}\mathcal{N}(h_{f_{\mathbf{G}}^{\text{RB}} \to h_{n,m}^{\text{RB}}(l)}, v_{f_{\mathbf{G}}^{\text{RB}} \to h_{n,m}^{\text{RB}}(l)} + \vartheta^{\text{RB}}(l))} \right)^{-1}, \end{aligned} \quad (50a) $$

$$ \phi_{n,m}^{\text{RB}}(l) = (\vartheta^{\text{RB}}(l)^{-1} + v_{f_{\mathbf{G}}^{\text{RB}} \to h_{n,m}^{\text{RB}}(l)}^{-1})^{-1}, \quad (50b) $$

$$ \mu_{n,m}^{\text{RB}}(l) = \phi_{n,m}^{\text{RB}}(l) v_{f_{\mathbf{G}}^{\text{RB}} \to h_{n,m}^{\text{RB}}(l)}^{-1} h_{f_{\mathbf{G}}^{\text{RB}} \to h_{n,m}^{\text{RB}}(l)}. \quad (50c) $$

Then, we obtain $b_{f_{h_{n,m}^{\text{RB}}(l)}}(h_{n,m}^{\text{RB}}(l)) = \mathcal{P}(b_{f_{h_{n,m}^{\text{RB}}(l)}}(h_{n,m}^{\text{RB}}(l)))$ with mean and variance respectively given by

$$ \hat{h}_{n,m}^{\text{RB}}(l) = \hat{\lambda}_{n,m}^{\text{RB}}(l) \mu_{n,m}^{\text{RB}}(l), \quad (51a) $$

$$ v_{h_{n,m}^{\text{RB}}(l)} = \hat{\lambda}_{n,m}^{\text{RB}}(l) (|\mu_{n,m}^{\text{RB}}(l)|^2 + \phi_{n,m}^{\text{RB}}(l)) - |\hat{h}_{n,m}^{\text{RB}}(l)|^2. \quad (51b) $$

Similar to the approximation of $v_{h_{n,m}^{\text{UR}}(l)}$ in Sec. IV-B.3, we approximate $v_{h_{n,m}^{\text{RB}}(l)}, \forall n, m, l$ by $\bar{v}_h^{\text{RB}} = \sum_{n,m,l} v_{h_{n,m}^{\text{RB}}(l)} / (\text{NML}_2)$. Then, we compute

$$ \begin{aligned} & \mathcal{M}_{f_{h_{n,m}^{\text{RB}}(l)} \to h_{n,m}^{\text{RB}}(l)}(h_{n,m}^{\text{RB}}(l)) \\ & = \mathcal{C}\mathcal{N}(h_{n,m}^{\text{RB}}(l); h_{f_{h_{n,m}^{\text{RB}}(l)} \to h_{n,m}^{\text{RB}}(l)}, v_{f_{h_{n,m}^{\text{RB}}(l)} \to h_{n,m}^{\text{RB}}(l)}), \end{aligned} \quad (52) $$

with

$$ h_{f_{h_{n,m}^{\text{RB}}(l)} \to h_{n,m}^{\text{RB}}(l)} = v_{f_{h_{n,m}^{\text{RB}}(l)} \to h_{n,m}^{\text{RB}}(l)} (\hat{h}_{n,m}^{\text{RB}}(l) / \bar{v}_h^{\text{RB}} - h_{f_{h_{n,m}^{\text{RB}}(l)} \to h_{n,m}^{\text{RB}}(l)} / v_{f_{h_{n,m}^{\text{RB}}(l)} \to h_{n,m}^{\text{RB}}(l)}), \quad (53a) $$

$$ v_{f_{h_{n,m}^{\text{RB}}(l)} \to h_{n,m}^{\text{RB}}(l)} = \left( (\bar{v}_h^{\text{RB}})^{-1} - v_{f_{h_{n,m}^{\text{RB}}(l)} \to h_{n,m}^{\text{RB}}(l)}^{-1} \right)^{-1}. \quad (53b) $$

Using $b_{f_{h_{n,m}^{\text{RB}}(l)}}(h_{n,m}^{\text{RB}}(l))$ in (49) and $\mathbf{G}^{\text{RB}} = \bar{\mathbf{G}}^{\text{RB}} + \mathbf{F}_2 \mathbf{U}_B \mathbf{H}_{\text{RB}}^{\text{H}} \mathbf{A}_B^{\text{H}}$ in (10c), we obtain $b_{\mathbf{G}^{\text{RB}}}(\mathbf{G}^{\text{RB}}) = \mathcal{C}\mathcal{N}(\mathbf{G}^{\text{RB}}; \hat{\mathbf{G}}^{\text{RB}}, \mathbf{\Sigma}_{\mathbf{G}^{\text{RB}}})$ with

$$ \hat{\mathbf{G}}^{\text{RB}} = \mathbf{G}^{\text{RB}} + \mathbf{F}_2 \mathbf{U}_B \hat{\mathbf{H}}^{\text{RB}} \mathbf{A}_B^{\text{H}}, \quad (54a) $$

$$ (\mathbf{\Sigma}_{\mathbf{G}}^{\text{RB}})^{-1} = \bar{v}_h^{\text{RB}} \mathbf{F}_2 \mathbf{F}_2^{\text{H}}, \quad (54b) $$

which completes the message derivations.
### C. EM Learning

In practice, the prior parameters $\omega = \{\lambda_1, \lambda_2, \vartheta^{\text{UR}}(l), \vartheta^{\text{RB}}(l)\}$ in the Bernoulli-Gaussian models (11)-(12) are typically unknown. We adopt the expectation maximization (EM) method [42, Chap. 9] to learn $\omega$ from the received signal. The EM learning is formulated as

$$ \omega^{(i+1)} = \arg\max_{\omega} \mathbb{E} \left[ \ln p(\mathbf{H}^{\text{UR}}, \mathbf{H}^{\text{RB}}, \mathbf{G}^{\text{RB}}, \mathbf{R}, \mathbf{Y}; \omega) \mid \mathbf{Y}; \omega^{(i)} \right] \quad (55) $$

where $\omega^{(i)}$ is the estimate of $\omega$ in the $i$-th iteration. $\mathbb{E}[\cdot|\mathbf{Y}; \omega^{(i)}]$ represents the expectation over the posterior distribution $p(\mathbf{H}^{\text{UR}}, \mathbf{H}^{\text{RB}}, \mathbf{G}^{\text{RB}}, \mathbf{R}|\mathbf{Y})$ with $\omega = \omega^{(i)}$. In the HVMP framework, we use the message product $b_{\mathbf{H}^{\text{UR}}}(\mathbf{H}^{\text{UR}})b_{\mathbf{H}^{\text{RB}}}(\mathbf{H}^{\text{RB}})b_{\mathbf{G}^{\text{RB}}}(\mathbf{G}^{\text{RB}})b_{\mathbf{R}}(\mathbf{R})$ (updated in the $i$-th iteration) to approximate the posterior distribution for computing $\mathbb{E}[\cdot|\mathbf{Y}; \omega^{(i)}]$. Then, we set the derivative of $\mathbb{E}[\cdot|\mathbf{Y}; \omega^{(i)}]$ w.r.t. $\lambda_1$ to be zero, which yields

$$ \lambda_1^{(i+1)} = \frac{1}{KNL_1} \sum_{k,n,l} \hat{\lambda}_{k,n}^{\text{UR}}(l). \quad (56) $$

Similar to (56), the EM estimate of $\lambda_2$ is given by

$$ \lambda_2^{(i+1)} = \frac{1}{MNL_2} \sum_{m,n,l} \hat{\lambda}_{m,n}^{\text{RB}}(l). \quad (57) $$

We set the derivative of $\mathbb{E}[\cdot|\mathbf{Y}; \omega^{(i)}]$ w.r.t. $\vartheta^{\text{UR}}(l)$ to be zero, which yields

$$ \vartheta^{\text{UR},(i+1)}(l) = \frac{\sum_{k,n} \hat{\lambda}_{k,n}^{\text{UR}}(l) (|\hat{h}_{k,n}^{\text{UR}}(l)|^2 + v_{h_{k,n}^{\text{UR}}(l)})}{\sum_{k,n} \hat{\lambda}_{k,n}^{\text{UR}}(l)}. \quad (58) $$

Similar to (58), the EM estimate of $\vartheta^{\text{RB}}(l)$ is given by

$$ \vartheta^{\text{RB},(i+1)}(l) = \frac{\sum_{n,m} \hat{\lambda}_{n,m}^{\text{RB}}(l) (|\hat{h}_{n,b}^{\text{RB}}(l)|^2 + v_{h_{n,b}^{\text{RB}}(l)})}{\sum_{n,b} \hat{\lambda}_{n,b}^{\text{RB}}(l)}. \quad (59) $$

### D. Complete Algorithm

We summarize the proposed HVMP algorithm in Algorithm 1. In each HVMP iteration, the message of $\mathbf{H}^{\text{UR}}$ in steps 2-3 is employed to update the message of $\mathbf{H}^{\text{RB}}$ in steps 5-6, and vice versa. In other words, the messages of $\mathbf{H}^{\text{UR}}$ and $\mathbf{H}^{\text{RB}}$ alternatively serve as soft pilots. To improve the convergence of the algorithm, we adopt the damping strategy given in [43, Sec. II] for the update of the mean and variance of $h_{k,n}^{\text{UR}}(l)$ in (39) and those of $h_{n,m}^{\text{RB}}(l)$ in (51). When the HVMP iteration converges, we output the means of $b_{\mathbf{H}^{\text{UR}}}(\mathbf{H}^{\text{UR}})$ and $b_{\mathbf{H}^{\text{RB}}}(\mathbf{H}^{\text{RB}})$ as the estimates of $\mathbf{H}^{\text{UR}}$ and $\mathbf{H}^{\text{RB}}$, respectively.

One advantage of Algorithm 1 is that the variable nodes and the factor nodes are matrices and joint distributions of matrices, respectively, which result in the joint Gaussian message $M_{f_Y \to R}(\mathbf{R})$ in (25) (or $M_{f_Y \to G^{\text{RB}}}(\mathbf{G}^{\text{RB}})$ in (43)). With $R_{f_Y \to R}$ (or $G_{f_Y \to G^{\text{RB}}}$) treated as the mean of $\mathbf{R}$ (or $\mathbf{G}^{\text{RB}}$) in the joint Gaussian message, the corresponding covariance characterizes the correlation of the entries of $\mathbf{R}$ (or $\mathbf{G}^{\text{RB}}$),

### Algorithm 1 Proposed HVMP Algorithm

$$ \text{Input: } \mathbf{Y}, \theta_t, x_{k,t}(j), \sigma_w^2 \\ \text{Initialization: } \mathbf{G}^{\text{RB}} = \mathbf{G}^{\text{RB}}, \Sigma_{\mathbf{G}}^{\text{RB}} = \mathbf{I}, h_{f_{h_{n,m}^{\text{RB}}}(l)} \to h_{n,m}^{\text{RB}}(l) = 0, \\ h_{f_{h_{n,m}^{\text{UR}}}(l)} \to h_{n,m}^{\text{UR}}(l) = 0, v_{f_{h_{n,m}^{\text{RB}}}(l)} \to h_{n,m}^{\text{RB}}(l) = 1, v_{f_{h_{n,m}^{\text{UR}}}(l)} \to h_{n,m}^{\text{UR}}(l) = 1. $$

while the stopping criterion is not met do

% Messages and beliefs for the estimation of $\mathbf{H}^{\text{UR}}$.

1: Update $M_{f_Y \to R}(\mathbf{R})$ by (25) and (29);

2: Update $M_{f_R \to H^{\text{UR}}}(\mathbf{H}^{\text{UR}})$ by (32)-(34);

3: Update $b_{\mathbf{H}^{\text{UR}}}(\mathbf{H}^{\text{UR}})$ by (37)-(39), $M_{f_{\mathbf{H}^{\text{UR}}} \to \mathbf{H}^{\text{UR}}}(\mathbf{H}^{\text{UR}})$ by (40)-(41), and $b_{\mathbf{R}}(\mathbf{R})$ by (42).

% Messages and beliefs for the estimation of $\mathbf{H}^{\text{RB}}$.

4: Update $M_{f_Y \to G^{\text{RB}}}(\mathbf{G}^{\text{RB}})$ by (43)-(44);

5: Update $M_{f_R \to H^{\text{RB}}}(\mathbf{H}^{\text{RB}})$ by (46)-(48);

6: Update $b_{\mathbf{H}^{\text{RB}}}(\mathbf{H}^{\text{RB}})$ by (50)-(51), $M_{f_{\mathbf{H}^{\text{RB}}} \to \mathbf{H}^{\text{RB}}}(\mathbf{H}^{\text{RB}})$ by (52)-(53), and $b_{\mathbf{G}^{\text{RB}}}(\mathbf{G}^{\text{RB}})$ by (54).

% EM learning.

7: Update $\{\lambda_1, \lambda_2, \vartheta^{\text{UR}}(l), \vartheta^{\text{RB}}(l)\}$ by (56)-(59).

end

Output: $\hat{\mathbf{H}}^{\text{UR}}$ and $\hat{\mathbf{H}}^{\text{RB}}$.

which is then considered in the computation of the belief of $\mathbf{H}^{\text{UR}}$ (or $\mathbf{H}^{\text{RB}}$), leading to an improvement of the algorithm convergence. For comparison, most existing message-passing methods, e.g., the bilinear generalized approximate message passing (BIG-AMP) [41], are designed based on a factor graph which consists of scalar variable nodes. As such, the corresponding messages are estimates of marginal pdfs that fail to characterize the correlation between variables.

We next analyze the computational complexity of HVMP. Specifically, in steps 1-2 of Algorithm 1, the matrix inversion dominates the complexity. Specifically, the computation of $\Sigma_{f_R, k}^{-1}, \forall k$ in (32) is $O(KL_1^3N^3)$. The main computations in step 3 are point-wise operations, where the computational complexity is linear with the parameters $K, N$, and $L_1$. Besides, the multiplication of a vector and a transformation matrix given in (6)-(9) can be fulfilled by using the fast Fourier transform (FFT). Steps 4-6 are similar to steps 1-3, and the corresponding computational complexity is dominated again by the matrix inversion in (46), which is $O(L_2^3N^3)$. Note that the matrix inversion can be efficiently calculated by using the Cholesky decomposition, whose computational complexity is comparable to a matrix multiplication, since the matrix to be computed is positive-definite. In summary, the computational complexity of the HVMP algorithm is $O(KL_1^3N^3) + O(L_2^3N^3)$ per iteration. To reduce the computational complexity, the sampling intervals in (7) and (9) can be increased, resulting in channel representations with smaller sizes. To avoid the missed detection of some path angles that do not match the predetermined sampling grid, we can employ a super-resolution method, similar to the approach in [44], to fine-tune the sampling grid.

## V. NUMERICAL RESULTS

### A. Numerical Results Under Channel Models (11) and (12)

In this subsection, we evaluate the performance of the proposed algorithm by assuming that the cascaded channels are generated by the probability models in (11) and (12) with
known prior parameters. We consider that the RIS is deployed
near the users [5]. We set $L_1 = 3$ for the MU-to-RIS link and
$L_2 = 12$ for the RIS-to-BS link. The angle domain sparsities
for $\mathbf{H}^{\text{UR}}$ and $\mathbf{H}^{\text{RB}}$ are set to $\lambda_1 = 0.5$ and $\lambda_2 = 0.1$,
respectively. The channel gain of the *dynamic channel* $\mathbf{H}^{\text{RB}}$
is 5 dB lower than that of the *static channel*. The *static*
channel has $L_2$ delay-taps, and the delay-angle domain
sparsity is 0.5 with the non-zero coefficients being i.i.d.
Gaussian random variables. We adopt a frequency-division
pilot assignment. More specifically, the pilot symbol $x_{k,t}(j)$ of
user $k$ is equal to $\exp(i\alpha_{k,t}(j))$ with $\alpha_{k,t}(j)$ uniformly taken
in $[0, 2\pi]$ if $j \in \{k, k+K, k+2K, \dots, k+(\frac{J}{K}-1)K\}$ and
is equal to 0 if $j \notin \{k, k+K, k+2K, \dots, k+(\frac{J}{K}-1)K\}$.
As for the RIS, the phase shift of each RIS element
is uniformly taken in $[0, 2\pi]$. We define SNR =
$\sum_j \frac{\mathbb{E}\|\mathbf{G}^{\text{RB}}(j)(\mathbf{\Theta} \odot (\mathbf{G}^{\text{UR}}(j)\mathbf{X}(j)))\|_F^2}{N \sum_j \mathbb{E}\|\mathbf{N}(j)\|_F^2}$.
Unless otherwise specified, we set $J = 288$, $\Delta f = 30$ KHz,
$M = 48$, $K = 18$, $T = 8$, and $N = 72$ (with $N_x = 12$ and
$N_y = 6$). The baseline schemes are as follows:

*   **LS:** Following the model in [17, eq. (4)], $\mathbf{h}_{\text{ds},k}^{\text{UR}}$ of user $k$ in (6) and $\mathbf{H}_{\text{ds}}^{\text{RB}}$ in (8) are combined into an effective channel of user $k$, which is estimated using the LS method based on the received signal $\mathbf{Y}(j)$, $j \in \{k, k+K, k+2K, \dots, k+(\frac{J}{K}-1)K\}$.

*   **MCMP:** Based on the received signal model in (5), the channel matrices $\tilde{\mathbf{G}}^{\text{RB}}(j)$, $\forall j$ and $\mathbf{G}^{\text{UR}}(j)$, $\forall j$ are estimated by using the matrix-calibration-based message-passing (MCMP) algorithm in [13]. Note that the MCMP algorithm is designed based on the BIG-AMP framework and exploits the angle domain sparsities of the cascaded channels. The pilot matrix $\mathbf{X}(j)$, $\forall j$ is a random matrix with i.i.d. Gaussian entries.

*   **Modified MCMP:** Each iteration of MCMP can be decomposed into two stages. Stage one gives coarse estimates of the cascaded channel variables from noisy observations; stage two refines the estimates of $\tilde{\mathbf{G}}^{\text{RB}}(j)$, $\forall j$ and $\mathbf{G}^{\text{UR}}(j)$, $\forall j$ using [13, eq. (55)] and [13, eq. (64)], respectively. We modify stage two by transforming the channel estimates (from stage 1) into the delay-angle domain based on (6)-(9). Then, we modify [13, eq. (55)] and [13, eq. (64)] by replacing the original prior probabilities with the prior probabilities in (11) and (12), respectively.

*   **Two-phase scheme:** Following [23], we select a typical user $k$ with the maximum channel power in the MU-to-RIS link. According to the two-phase scheme in [23], we use one path of the *static* channel (with maximum power) to assist the estimation of $\{\beta_{k,l}^{\text{UR}}, \phi_{k,l}^{\text{UR}}, \theta_{k,l}^{\text{UR}}, \tau_{k,l}^{\text{UR}}\}_{l}$ in (1) for the typical user $k$, and then the estimated parameters are used to estimate $\{\tilde{\beta}_{l}^{\text{RB}}, \tilde{\phi}_{l}^{\text{RB}}, \tilde{\theta}_{l}^{\text{RB}}, \tilde{\psi}_{l}^{\text{RB}}, \tilde{\tau}_{l}^{\text{RB}}\}_{l}$ in (4). Finally, we estimate $\{\beta_{i,l}^{\text{UR}}, \phi_{i,l}^{\text{UR}}, \theta_{i,l}^{\text{UR}}, \tau_{i,l}^{\text{UR}}\}_{l}$ for user $i \neq k$. Given these parameter estimates, we reconstruct $\mathbf{G}^{\text{UR}}(j)$, $\forall j$ and $\mathbf{G}^{\text{RB}}(j)$, $\forall j$ based on (1)-(4).

*   **Oracle bounds:** The oracle bound of $\mathbf{H}^{\text{UR}}$ represents the MSE of $\mathbf{H}^{\text{UR}}$ by executing Algorithm 1 with perfectly known $\mathbf{H}^{\text{RB}}$; the oracle bound of $\mathbf{H}^{\text{RB}}$ represents the MSE of $\mathbf{H}^{\text{RB}}$ by executing Algorithm 1 with perfectly known $\mathbf{H}^{\text{UR}}$.

<figure><img src="image_6.png" /><figcaption>Fig. 5. MSE of channel estimation versus the iteration number for $\mathbf{G}^{\text{RB}}(j)$ and $\mathbf{G}^{\text{UR}}(j)$ in the left and right sub-plots, respectively. SNR=-2 dB.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. MSE of channel estimation versus the SNR in the left and right sub-plots, respectively.</figcaption></figure>

Fig. 5 shows the MSE of the channel estimation against the iteration number, where the left and right sub-plots provide the MSE performance of estimating $\mathbf{G}^{\text{RB}}(j)$ and $\mathbf{G}^{\text{UR}}(j)$, respectively. We see that the proposed algorithm converges within a few iterations, and converges to the oracle bounds. On the contrary, the baseline schemes exhibit poor performance. Particularly, the poor performance of the MCMP method is due to the ignorance of the delay domain sparsity. The modified MCMP method improves the performance by exploiting the delay domain sparsity but still performs much worse than the HVMP scheme. This is because the BIG-AMP algorithm used in the MCMP method can only handle scalar messages instead of matrix messages as in HVMP. The two-phase scheme fails to work since it assumes a known LOS path in the RIS-to-BS link to assist the channel estimation, which is not applicable to the setting of multiple non-line-of-sight (NLOS) paths with comparable powers.

In Fig. 6, we present the MSE of channel estimation against the SNR, where the left and right sub-plots show the MSE of estimating $\mathbf{G}^{\text{RB}}(j)$ and $\mathbf{G}^{\text{UR}}(j)$, respectively. As the SNR increases, the performance gap between the proposed algorithm and the baseline schemes becomes larger. Besides, the proposed algorithm consistently approaches the oracle bounds within the SNR range of [-15, 5] dB.

Among the considered baseline schemes, the LS method
estimates the equivalent channel coefficients from a linear
<figure><img src="image_8.png" /><figcaption>Fig. 7. MSE of channel estimation versus the SNR for $G^{\text{RB}}(j) G^{\text{UR}}(j)$. $\lambda_1 = 0.5$ and $\lambda_2 = 0.1$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. MSE of channel estimation versus the SNR for $G^{\text{RB}}(j) G^{\text{UR}}(j)$. $\lambda_1 = 0.25$ and $\lambda_2 = 0.05$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. MSE of channel estimation versus the number of OFDM pilot symbols for $G^{\text{RB}}(j) G^{\text{UR}}(j)$. SNR=0 dB.</figcaption></figure>

representation of the two cascaded channels. To compare the LS method with the HVMP method, we report the MSE of the estimate of $G^{\text{UR}}(j)G^{\text{RB}}(j)$ in Fig. 7. The trend is similar to that in Fig. 6, where the proposed algorithm significantly outperforms the baseline schemes. In Fig. 8, we reduce $\lambda_1$ and $\lambda_2$ to 0.25 and 0.05, respectively. Compared to the results in Fig. 7, the performance of the proposed method improves by nearly 3 dB due to the increase of channel sparsity.

Fig. 9 shows the estimation performance of the cascaded channels by varying the number of pilot symbols. Particularly, the proposed algorithm achieves MSE = -20 dB using $T = 4$ pilot symbols while the performance of the baseline schemes is above -5 dB when $T < 14$. These observations imply that the proposed scheme can significantly reduce the required pilot overhead while maintaining reasonably accurate channel estimates.

## B. Numerical Results Under a More Practical Channel Model

In this subsection, we adopt a more practical channel model based on 3GPP TR 38.901 [45]. The carrier frequency is set

to 6 GHz. The distance of the RIS-to-BS link is 300 m, and the distance of each user-to-RIS link is randomly generated within the region of [10, 50] m. For the MU-to-RIS channel $H^{\text{UR}}$, we set a LOS path and multiple scattering paths due to a short distance between the users and the RIS. The Rice factor follows [45, Table 7.5-6] for the UMa LOS scenario. For the scattering paths, we adopt the TDL-B channel model in [45, Table 7.7.2-2] with 23 channel paths. The root mean square (r.m.s.) delay spread of these 23 paths is 30 ns which corresponds to a relatively short delay spread [45, Table 7.3.3-1]. Each path consists of 40 sub-paths, has a central azimuth/zenith angle randomly taken in $(-\pi/2, \pi/2)$, and has a 15-degree angle spread. For the RIS-to-BS channel, we consider multiple scattering paths due to the longer distance between the RIS and the BS. Specifically, the channel model of the *static* channel component is the TDL-B model with 23 channel paths and a 150 ns r.m.s. delay spread which corresponds to a relatively nominal delay spread [45]. The number of paths of the *dynamic* channel is 2, where the path delays are randomly taken within the range of the delay spread of the *static* channel. Each path comprises 40 sub-paths, and the sub-path angles are generated following those of the MU-to-RIS channel. The EM method is adopted to learn the prior parameters with initialization $\lambda_1 = 0.5, \lambda_2 = 0.5, \vartheta^{\text{UR}}(l) = 1, \forall l$, and $\vartheta^{\text{RB}}(l) = 0.1, \forall l$. The parameters $J, \Delta f, M, K, T$, and $N$ are set as in Subsection A.

Fig. 10 shows the MSE of channel estimation against the iteration number, where the left and right sub-plots provide the MSE of estimating $G^{\text{RB}}(j)$ and $G^{\text{UR}}(j)$, respectively. It is seen that the MCMP scheme and the modified MCMP scheme perform much worse than the HVMP algorithm. It is also seen that the proposed scheme converges within a few iterations. This demonstrates the advantage of considering matrix messages in HVMP.

Fig. 11 shows the MSE of the estimation of $G^{\text{RB}}(j)G^{\text{UR}}(j)$. As the SNR increases, the performance of the proposed scheme approaches -18 dB while the performance of the baseline schemes is above -7 dB. Fig. 12
<figure><img src="image_11.png" /><figcaption>Fig. 10. MSE of channel estimation versus the iteration number for $G^{\text{RB}}(j)$ and $G^{\text{UR}}(j)$ in the left and right sub-plots, respectively. SNR=-4 dB.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 11. MSE of channel estimation versus the SNR for $G^{\text{RB}}(j)$ and $G^{\text{UR}}(j)$.</figcaption></figure>

shows the MSE of estimating $G^{\text{RB}}(j)G^{\text{UR}}(j)$ for a varying number of OFDM pilot symbols. We see that the proposed scheme achieves MSE ≈ -15 dB using $T = 6$ pilot symbols, providing almost 10 dB gain compared with the baseline schemes. The obtained results indicate that the proposed algorithm is able to effectively reduce the pilot overhead while maintaining excellent channel estimation performance.

## VI. CONCLUSION

This paper studied the problem of cascaded channel estimation in RIS-aided MU MIMO-OFDM systems. We established an equivalent signal model which comprises one subcarrier-wise bilinear model and two linear models for the MU-to-RIS and RIS-to-BS channels. The proposed equivalent signal model maintains the matrix factorization structure of the two cascaded channels in the delay-angle domain, where the number of non-zero channel coefficients is proportional to the sum of the numbers of paths in the two cascaded links. Based on the established model, we leveraged the Bethe method to design the HVMP algorithm for the cascaded channel estimation. In each iteration of HVMP, the covariance of each channel matrix is estimated to characterize the correlation of the matrix elements. Simulation results showed that the HVMP algorithm can significantly reduce the

<figure><img src="image_13.png" /><figcaption>Fig. 12. MSE of channel estimation versus the number of OFDM pilot symbols for $G^{\text{RB}}(j)$ and $G^{\text{UR}}(j)$. SNR=-5 dB.</figcaption></figure>

required pilot overhead compared to various baseline schemes, and can approach the oracle estimation bounds for the MU-to-RIS and RIS-to-BS channels.

## APPENDIX A
### PROOF OF THEOREM 1

We substitute (19) into (17) and construct the Lagrange function as

$$ \begin{aligned} \mathcal{L} = & \sum_{f \in \mathcal{F} \setminus \mathcal{F}_Y} \int_{\mathcal{U}} b_f(\mathbf{u}) \ln \frac{b_f(\mathbf{u})}{f(\mathbf{u})} - \sum_{\mathbf{v} \in \mathcal{V}} \int_{\mathcal{V}} b_v(\mathbf{v}) \ln b_v(\mathbf{v}) \\ & + \int_{\mathbf{R}, \mathbf{G}^{\text{RB}}} b_{f_{\mathbf{Y}, \mathbf{R}}}(\mathbf{R}) b_{f_{\mathbf{Y}, \mathbf{G}}}(\mathbf{G}^{\text{RB}}) \ln \frac{b_{f_{\mathbf{Y}, \mathbf{R}}}(\mathbf{R}) b_{f_{\mathbf{Y}, \mathbf{G}}}(\mathbf{G}^{\text{RB}})}{p(\mathbf{Y} | \mathbf{R}, \mathbf{G}^{\text{RB}})} \\ & + \sum_{\mathbf{v} \in \mathcal{V}} \sum_{f \in \mathcal{N}(\mathbf{v})} (\lambda_{fv}^{(1)})^H \int_{\mathcal{V}} \mathbf{v} \left( b_v(\mathbf{v}) - \int_{\mathcal{U} \setminus \mathcal{V}} b_f(\mathbf{u}) \right) \\ & + \sum_{\mathbf{v} \in \mathcal{V}} \sum_{f \in \mathcal{N}(\mathbf{v})} (\lambda_{fv}^{(2)})^H \text{vec} \left( \int_{\mathcal{V}} \mathbf{v} \mathbf{v}^H (b_v(\mathbf{v}) - \int_{\mathcal{U} \setminus \mathcal{V}} b_f(\mathbf{u})) \right), \end{aligned} $$

where $\lambda_{fv}^{(1)}$ and $\lambda_{fv}^{(2)}$ are the Lagrange multipliers for the moment-matching constraint in (18).

Setting the derivative of $\mathcal{L}$ with respect to $b_f(\mathbf{u})$, $\forall f \in \mathcal{F} \setminus \mathcal{F}_Y$ equal to zero, we obtain

$$ b_f(\mathbf{u}) \propto f(\mathbf{u}) \prod_{\mathbf{v} \in \mathcal{N}(f)} e^{(\lambda_{fv}^{(1)})^H \mathbf{v} + (\lambda_{fv}^{(2)})^H \text{vec}(\mathbf{v}\mathbf{v}^H)}. \quad (60) $$

Setting the derivative of $\mathcal{L}$ with respect to $b_{f_{\mathbf{Y}, \mathbf{v}}}(\mathbf{v})$, $\forall \mathbf{v} \in \{\mathbf{R}, \mathbf{G}^{\text{RB}}\}$ equal to zero, we obtain

$$ b_{f_{\mathbf{Y}, \mathbf{v}}}(\mathbf{v}) \propto e^{\int \prod_{\tilde{\mathbf{v}} \in \mathcal{N}(\mathbf{f}_{\mathbf{Y}}) \setminus \mathbf{v}} b_{f_{\mathbf{Y}, \tilde{\mathbf{v}}}}(\tilde{\mathbf{v}}) \ln f_{\mathbf{Y}}(\mathbf{u}) du \setminus \mathbf{v}} \\ \times e^{(\lambda_{fv}^{(1)})^H \mathbf{v} + (\lambda_{fv}^{(2)})^H \text{vec}(\mathbf{v}\mathbf{v}^H)}. \quad (61) $$

Setting the derivative of $\mathcal{L}$ with respect to $b_v(\mathbf{v})$, $\forall \mathbf{v} \in \mathcal{V}$ equal to zero, we obtain

$$ b_v(\mathbf{v}) \propto \prod_{f \in \mathcal{N}(\mathbf{v})} e^{(\lambda_{fv}^{(1)})^H \mathbf{v} + (\lambda_{fv}^{(2)})^H \text{vec}(\mathbf{v}\mathbf{v}^H)}. \quad (62) $$
To simplify the expressions in (60)-(62), we define the quantities

$$
\begin{align}
\mathcal{M}_{v \to f}(v) &= \prod_{\tilde{f} \in \mathcal{N}(v) \setminus f} \mathcal{M}_{\tilde{f} \to v}(v) \nonumber \\
&\triangleq e^{(\lambda_{f\tilde{v}}^{(1)})^H v + (\lambda_{f\tilde{v}}^{(2)})^H \text{vec}(v\tilde{v}^H)}, \tag{63}
\end{align}
$$

where $\mathcal{M}_{v \to f}(v)$ and $\mathcal{M}_{f \to v}(v)$ are Gaussian pdfs. Substituting (63) into (60)-(62) yields

$$
b_f(u) \propto f(u) \prod_{v \in \mathcal{N}(f)} \mathcal{M}_{v \to f}(v), \quad f \in \mathcal{F} \setminus f_Y, \quad (64a)
$$

$$
b_{f_{\tilde{Y}},v}(v) \propto e^{\int \prod_{\tilde{v} \in \mathcal{N}(f_Y) \setminus v} b_{f_{\tilde{Y}},\tilde{v}}(\tilde{v}) \ln f_Y(u) du \setminus v} \mathcal{M}_{v \to f_Y}(v), \quad (64b)
$$

$$
b_v(v) \propto \prod_{f \in \mathcal{N}(v)} \mathcal{M}_{v \to f}(v) = \prod_{f \in \mathcal{N}(v)} \mathcal{M}_{f \to v}(v). \quad (64c)
$$

where (64c) follows from the fact that each factor
node $f$ is connected to two variable nodes yielding
$\prod_{f \in \mathcal{N}(v)} \mathcal{M}_{v \to f}(v) = \prod_{f \in \mathcal{N}(v)} \mathcal{M}_{\tilde{f} \to v}(\tilde{v}) = \prod_{f \in \mathcal{N}(v)} \mathcal{M}_{f \to v}(v)$.

We next consider the moment-matching constraint between
the beliefs $b_f(\mathbf{u})$, $\forall \mathbf{u} \in \mathcal{F} \setminus f_Y$ and $b_v(\mathbf{v})$. Based on (63),
we rewrite (64c) as

$$
\begin{align}
b_v(v) \propto \prod_{f \in \mathcal{N}(v)} \mathcal{M}_{f \to v}(v) &= \mathcal{M}_{f \to v}(v) \prod_{\tilde{f} \in \mathcal{N}(v) \setminus f} \mathcal{M}_{\tilde{f} \to v}(\tilde{v}) \nonumber \\
&= \mathcal{M}_{f \to v}(v) \mathcal{M}_{v \to f}(v). \tag{65}
\end{align}
$$

Based on (64a), we have

$$
\int b_f(\mathbf{u}) du \setminus v \propto \int f(\mathbf{u}) \prod_{\tilde{v} \in \mathcal{N}(f)} \mathcal{M}_{\tilde{v} \to f}(\tilde{v}) du \setminus v, \forall f \in \mathcal{F} \setminus f_Y. \quad (66)
$$

Recall that $\mathcal{M}_{v \to f}(v)$ and $\mathcal{M}_{f \to v}(v)$ are Gaussian pdfs, leading to a Gaussian $b_v(v)$ in (65). Given (65)-(66), the moment-matching constraint between $b_v(v)$ and $b_f(u)$, $\forall f \in \mathcal{F} \setminus f_Y$ is equivalently described as

$$
\begin{align}
\mathcal{M}_{f \to v}(v) \mathcal{M}_{v \to f}(v) \propto & \mathcal{P} \left( \int f(u) \prod_{\tilde{v} \in \mathcal{N}(f)} \mathcal{M}_{\tilde{v} \to f}(\tilde{v}) du \setminus v \right), \nonumber \\
& \forall f \in \mathcal{F} \setminus f_Y. \tag{67}
\end{align}
$$

Similar to (65)-(67), the moment-matching constraint
between $b_v(v)$ and $b_{f_Y, v}(v)$ is equivalently described as

$$
\begin{align*}
\mathcal{M}_{f_Y \to v}(v) \mathcal{M}_{v \to f_Y}(v) \\
&\propto \mathcal{P}(b_{f_Y, v}(v)), \\
&\propto \mathcal{P}\left(e^{\int \prod_{\tilde{v} \in \mathcal{N}(f_Y) \setminus v} b_{f_Y, \tilde{v}}(\tilde{v}) \ln f_Y(u) du \setminus v} \mathcal{M}_{v \to f_Y}(v)\right),
\end{align*}
$$

which completes the proof.

APPENDIX B

DERIVATION OF (25)-(26)

Given $f_Y(\mathbf{R}, \mathbf{G}^{\text{RB}}) = \mathcal{CN}(\mathbf{Y}; \text{diag}(\mathbf{G}^{\text{RB}}(j))\mathbf{R}, \beta_1\mathbf{I})$ in (10a) and $b_{\text{GRB}}(\mathbf{G}^{\text{RB}}) = \mathcal{CN}(\mathbf{G}^{\text{RB}}; \hat{\mathbf{G}}^{\text{RB}}, \boldsymbol{\Sigma}_{\mathbf{G}}^{\text{RB}})$ in (54), $\mathcal{M}_{f_Y \to \mathbf{R}}(\mathbf{R})$ in (23) is expressed as

$$
\mathcal{M}_{f_Y \to \mathbf{R}}(\mathbf{R})
$$

$$
\begin{align*}
&\propto \exp\left(\int b_{\mathbf{G}^{\text{RB}}}(\mathbf{G}^{\text{RB}}) \ln f_Y(\mathbf{G}^{\text{RB}}, \mathbf{R}) d\mathbf{G}^{\text{RB}}\right), \\
&\propto \exp\left(\int -b_{\mathbf{G}^{\text{RB}}}(\mathbf{G}^{\text{RB}}) \frac{1}{\beta_1} \text{Tr}\left((\mathbf{Y} - \text{diag}(\mathbf{G}^{\text{RB}}(j))\mathbf{R})^H \right. \right. \\
&\qquad \left. \left. \times (\mathbf{Y} - \text{diag}(\mathbf{G}^{\text{RB}}(j))\mathbf{R})\right) d\mathbf{G}^{\text{RB}}\right), \\
&\propto \exp\left(\frac{1}{\beta_1} \text{Tr}\left(2\mathcal{R}(\mathbf{R}^H \text{diag}(\hat{\mathbf{G}}^{\text{RB}}(j)^H)\mathbf{Y}) - \mathbf{R}^H \right. \right. \\
&\qquad \left. \left. \times \text{diag}(\mathbb{E}_{b_{\mathbf{G}^{\text{RB}}}}[\mathbf{G}^{\text{RB}}(j)^H \mathbf{G}^{\text{RB}}(j)])\mathbf{R})\right), \\
&\propto \mathcal{CN}(\mathbf{R}; \mathbf{R}_{f_Y \to \mathbf{R}}, \boldsymbol{\Sigma}_{f_Y \to \mathbf{R}}).
\end{align*}
$$

APPENDIX C
DERIVATION OF (31)-(32)

Recall that the relationship between **R**(j), ∀j and **r**<sub>t</sub>, ∀t, is
[**R**(1)<sup>T</sup>, . . . , **R**(J)<sup>T</sup>]<sup>T</sup> = [**r**<sub>1</sub>, . . . , **r**<sub>T</sub>]. Given this relationship
together with (10b) and (32c), we obtain p(**R**|**H**<sup>UR</sup>) =
Π<sub>t</sub> δ(**r**<sub>t</sub> − Σ<sub>k</sub> **S**<sub>k,t</sub>**h**<sub>k</sub><sup>UR</sup>). Plugging p(**R**|**H**<sup>UR</sup>) = Π<sub>t</sub> δ(**r**<sub>t</sub> −
Σ<sub>k</sub> **S**<sub>k,t</sub>**h**<sub>k</sub><sup>UR</sup>), (25) and (40) into the left side of (31) gives

$$
\prod_t \mathcal{CN} \left( \sum_k \mathbf{S}_{k,t} \mathbf{h}_k^{\text{UR}} ; \mathbf{r}_{f_Y \to r_t}, \boldsymbol{\Sigma}_{f_Y \to \mathbf{R}} \right) \\
\times \prod_k \mathcal{CN} \left( \mathbf{h}_k^{\text{UR}} ; \mathbf{h}_{\mathbf{h}_k^{\text{UR}} \to \mathbf{h}_k^{\text{UR}}}^{\text{UR}}, \boldsymbol{\Sigma}_{\mathbf{h}_k^{\text{UR}} \to \mathbf{h}_k^{\text{UR}}}^{\text{UR}} \right). \quad (68)
$$

Based on the pilot orthogonality in Sec II-A, i.e.,
$x_{k_1,t}(j)x_{k_2,t}(j) = 0, \forall k_1 \neq k_2$, the first term in (68) is
simplified as

$$
\begin{align}
& \prod_t \mathcal{CN} \left( \sum_k \mathbf{S}_{k,t} \mathbf{h}_k^{\text{UR}} ; \mathbf{r}_{f_Y \to r_t}, \boldsymbol{\Sigma}_{f_Y \to \mathbf{R}} \right) \nonumber \\
& \propto \exp \left( - \sum_t \left( \mathbf{r}_{f_Y \to r_t} - \sum_k \mathbf{S}_{k,t} \mathbf{h}_k^{\text{UR}} \right)^H \boldsymbol{\Sigma}_{f_Y \to \mathbf{R}} \left( \mathbf{r}_{f_Y \to r_t} - \sum_k \mathbf{S}_{k,t} \mathbf{h}_k^{\text{UR}} \right) \right), \tag{69a}
\end{align}
$$

$$
\begin{align}
& \propto \exp\left(-\sum_k \left( (\mathbf{h}_k^{\text{UR}})^H \left(\sum_t \mathbf{S}_{k,t}^H \boldsymbol{\Sigma}_{f_Y \to \mathbf{R}} \mathbf{S}_{k,t}\right) \mathbf{h}_k^{\text{UR}} - 2\text{Re}\left(\sum_t \mathbf{r}_{f_Y \to r_t}^H \boldsymbol{\Sigma}_{f_Y \to \mathbf{R}} \mathbf{S}_{k,t} \mathbf{h}_k^{\text{UR}}\right) \right)\right), \tag{69b} \\
& \propto \prod_k \mathcal{CN} \left( \mathbf{h}_k^{\text{UR}} ; \boldsymbol{\Sigma}_S^{-1} \sum_t \mathbf{S}_{k,t}^H \boldsymbol{\Sigma}_{f_Y \to \mathbf{R}} \mathbf{r}_{f_Y \to r_t} ; \boldsymbol{\Sigma}_S \right), \tag{69c}
\end{align}
$$

where (69b) follows from the property $\mathbf{S}_{k_1,t}^H \boldsymbol{\Sigma}_{f_Y \to \mathbf{R}} \mathbf{S}_{k_2,t} = 0, \forall k_1 \neq k_2$ which is obtained by combining (29b), (32c), and $x_{k_1,t}(j)x_{k_2,t}(j) = 0, \forall k_1 \neq k_2; \boldsymbol{\Sigma}_S = \sum_t \mathbf{S}_{k,t}^H \boldsymbol{\Sigma}_{f_Y \to \mathbf{R}} \mathbf{S}_{k,t}$ in (69c). Then, by substituting (69c) into (68) and applying the Gaussian combining property, we obtain (31) and (32).

REFERENCES

[1] M. Di Renzo et al., “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, Nov. 2020.

[2] X. Yuan, Y.-J. A. Zhang, Y. Shi, W. Yan, and H. Liu, “Reconfigurable-intelligent-surface empowered wireless communications: Challenges and opportunities,” IEEE Wireless Commun., vol. 28, no. 2, pp. 136–143, Apr. 2021.
[3] C. Pan et al., "An overview of signal processing techniques for RIS/IRS-aided wireless systems," *IEEE J. Sel. Topics Signal Process.*, vol. 16, no. 5, pp. 883-917, Aug. 2022.

[4] Q. Wu and R. Zhang, "Intelligent reflecting surface enhanced wireless network via joint active and passive beamforming," *IEEE Trans. Wireless Commun.*, vol. 18, no. 11, pp. 5394-5409, Nov. 2019.

[5] S. Zhang and R. Zhang, "Capacity characterization for intelligent reflecting surface aided MIMO communication," *IEEE J. Sel. Areas Commun.*, vol. 38, no. 8, pp. 1823-1838, Aug. 2020.

[6] A. Abrardo, D. Dardari, and M. Di Renzo, "Intelligent reflecting surfaces: Sum-rate optimization based on statistical position information," *IEEE Trans. Commun.*, vol. 69, no. 10, pp. 7121-7136, Oct. 2021.

[7] W. Yan, X. Yuan, and X. Kuai, "Passive beamforming and information transfer via large intelligent surface," *IEEE Wireless Commun. Lett.*, vol. 9, no. 4, pp. 533-537, Apr. 2020.

[8] S. Guo, S. Lv, H. Zhang, J. Ye, and P. Zhang, "Reflecting modulation," *IEEE J. Sel. Areas Commun.*, vol. 38, no. 11, pp. 2548-2561, Nov. 2020.

[9] W. Jiang and X. Yuan, "Simultaneous active and passive information transfer for RIS-aided MIMO systems: Iterative decoding and evolution analysis," *IEEE Trans. Signal Process.*, vol. 71, pp. 1683-1698, 2023.

[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] 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.

[12] Y. Guo et al., "Efficient channel estimation for RIS-aided MIMO communications with unitary approximate message passing," *IEEE Trans. Wireless Commun.*, vol. 22, no. 2, pp. 1403-1416, Feb. 2023.

[13] H. Liu, X. Yuan, and Y.-J. A. Zhang, "Matrix-calibration-based cascaded channel estimation for reconfigurable intelligent surface assisted multiuser MIMO," *IEEE J. Sel. Areas Commun.*, vol. 38, no. 11, pp. 2621-2636, Nov. 2020.

[14] Z. Peng et al., "Channel estimation for RIS-aided multi-user mmWave systems with uniform planar arrays," *IEEE Trans. Commun.*, vol. 70, no. 12, pp. 8105-8122, Dec. 2022.

[15] Z. Peng et al., "Two-stage channel estimation for RIS-aided multiuser mmWave systems with reduced error propagation and pilot overhead," *IEEE Trans. Signal Process.*, vol. 71, pp. 3607-3622, 2023.

[16] 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.

[17] B. Zheng, C. You, and R. Zhang, "Fast channel estimation for IRS-assisted OFDM," *IEEE Wireless Commun. Lett.*, vol. 10, no. 3, pp. 580-584, Mar. 2021.

[18] H. Dong, C. Ji, L. Zhou, J. Dai, and Z. Ye, "Sparse channel estimation with surface clustering for IRS-assisted OFDM systems," *IEEE Trans. Commun.*, vol. 71, no. 2, pp. 1083-1095, Feb. 2023.

[19] X. Xu, S. Zhang, F. Gao, and J. Wang, "Sparse Bayesian learning based channel extrapolation for RIS assisted MIMO-OFDM," *IEEE Trans. Commun.*, vol. 70, no. 8, pp. 5498-5513, Aug. 2022.

[20] 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, Jun. 2022.

[21] Y. Lin, S. Jin, M. Matthaiou, and X. You, "Channel estimation and user localization for IRS-assisted MIMO-OFDM systems," *IEEE Trans. Wireless Commun.*, vol. 21, no. 4, pp. 2320-2335, Apr. 2022.

[22] S. Liu, Z. Gao, J. Zhang, M. D. Renzo, and M. Alouini, "Deep denoising neural network assisted compressive channel estimation for mmWave intelligent reflecting surfaces," *IEEE Trans. Veh. Technol.*, vol. 69, no. 8, pp. 9223-9228, Aug. 2020.

[23] T. Ling, Y. Han, S. Jin, and M. Matthaiou, "Two-phase parameter-based separate channel estimation in RIS-aided MIMO OFDM systems," in *Proc. IEEE Int. Conf. Commun.*, May 2023, pp. 4329-4334.

[24] J. S. Yedidia, W. T. Freeman, and Y. Weiss, "Constructing free-energy approximations and generalized belief propagation algorithms," *IEEE Trans. Inf. Theory*, vol. 51, no. 7, pp. 2282-2312, Jul. 2005.

[25] H. Jiang, X. Yuan, and Q. Guo, "Generalized bilinear factorization via hybrid vector message passing," *IEEE Trans. Signal Process.*, vol. 72, pp. 5675-5690, 2024.

[26] G. Gradoni and M. Di Renzo, "End-to-end mutual coupling aware communication model for reconfigurable intelligent surfaces: An electromagnetic-compliant approach based on mutual impedances," *IEEE Wireless Commun. Lett.*, vol. 10, no. 5, pp. 938-942, May 2021.

[27] E. Basar and I. Yildirim, "Reconfigurable intelligent surfaces for future wireless networks: A channel modeling perspective," *IEEE Wireless Commun.*, vol. 28, no. 3, pp. 108-114, Jun. 2021.

[28] B. Xiong et al., "A statistical MIMO channel model for reconfigurable intelligent surface assisted wireless communications," *IEEE Trans. Commun.*, vol. 70, no. 2, pp. 1360-1375, Feb. 2022.

[29] C. A. Balanis, *Antenna Theory: Analysis and Design*. Hoboken, NJ, USA: Wiley, 2016.

[30] A. Albanese, F. Devoti, V. Sciancalepore, M. Di Renzo, and X. Costa-Pérez, "MARISA: A self-configuring metasurfaces absorption and reflection solution towards 6G," in *Proc. IEEE Conf. Comput. Commun.*, May 2022, pp. 250-259.

[31] A. Albanese, F. Devoti, V. Sciancalepore, M. D. Renzo, A. Banchs, and X. Costa-Pérez, "ARES: Autonomous RIS solution with energy harvesting and self-configuration towards 6G," in *Proc. IEEE Trans. Mobile Comput.*, May 2024, pp. 1-14.

[32] *Physical Channels and Modulation*, Standard TS 38.211, 3GPP, Apr. 2017.

[33] D. Tse and P. Viswanath, *Fundamentals of Wireless Communication*. Cambridge, U.K.: Cambridge Univ. Press, 2005.

[34] E. Riegler, G. E. Kirkelund, C. N. Manchon, M.-A. Badiu, and B. H. Fleury, "Merging belief propagation and the mean field approximation: A free energy approach," *IEEE Trans. Inf. Theory*, vol. 59, no. 1, pp. 588-602, Jan. 2013.

[35] D. Zhang, X. Song, W. Wang, G. Fettweis, and X. Gao, "Unifying message passing algorithms under the framework of constrained Bethe free energy minimization," *IEEE Trans. Wireless Commun.*, vol. 20, no. 7, pp. 4144-4158, Jul. 2021.

[36] J. Winn and C. Bishop, and T. Jaakkola, "Variational message passing," *J. Mach. Learn. Res.*, vol. 6, no. 23, pp. 661-694, Apr. 2005.

[37] T. P. Minka, "Expectation propagation for approximate Bayesian inference," in *Proc. 7th Conf. Uncertainty Artif. Intell.*, San Francisco, CA, USA, Aug. 2001, pp. 362-369.

[38] T. Minka, "Divergence measures and message passing," *Citeseer*, Princeton, NY, USA, Tech. Rep. MSR-TR-2005-173, 2005.

[39] J. T. Parker, P. Schniter, and V. Cevher, "Bilinear generalized approximate message passing—Part II: Applications," *IEEE Trans. Signal Process.*, vol. 62, no. 22, pp. 5854-5867, Nov. 2014.

[40] S. Rangan, "Generalized approximate message passing for estimation with random linear mixing," in *Proc. IEEE Int. Symp. Inf. Theory*, Jul. 2011, pp. 2168-2172.

[41] J. T. Parker, P. Schniter, and V. Cevher, "Bilinear generalized approximate message passing—Part I: Derivation," *IEEE Trans. Signal Process.*, vol. 62, no. 22, pp. 5839-5853, Nov. 2014.

[42] C. M. Bishop, *Pattern Recognition and Machine Learning*, vol. 2. Cham, Switzerland: Springer, 2006, pp. 1122-1128.

[43] S. Rangan, P. Schniter, A. K. Fletcher, and S. Sarkar, "On the convergence of approximate message passing with arbitrary matrices," *IEEE Trans. Inf. Theory*, vol. 65, no. 9, pp. 5339-5351, Sep. 2019.

[44] H. Liu, X. Yuan, and Y. J. Zhang, "Super-resolution blind channel-and-signal estimation for massive MIMO with one-dimensional antenna array," *IEEE Trans. Signal Process.*, vol. 67, no. 17, pp. 4433-4448, Sep. 2019.

[45] *5G: Study on Channel Model for Frequencies From 0.5 To 100 GHz*, Standard TR 38.901, 3GPP, Nov. 2020.

<figure><img src="image_14.png" /></figure>

Wenjun Jiang (Graduate Student Member, IEEE) received the B.S. degree in communication engineering from the University of Electronic Science and Technology of China in 2019, where he is currently pursuing the Ph.D. degree in electrical and information engineering with the National Key Laboratory of Wireless Communications. His research interests include statistical signal processing, wireless communications, and machine learning.
<figure><img src="image_15.png" /></figure>

**Xiaojun Yuan** (Senior Member, IEEE) received the Ph.D. degree in electrical engineering from the City University of Hong Kong in 2009. From 2009 to 2011, he was a Research Fellow with the Department of Electronic Engineering, City University of Hong Kong. He was also a Visiting Scholar with the Department of Electrical Engineering, University of Hawaii, Manoa, in Spring 2009, Summer 2009, and 2010. From 2011 to 2014, he was a Research Assistant Professor with the Institute of Network Coding, The Chinese University

of Hong Kong. From 2014 to 2017, he was an Assistant Professor with the School of Information Science and Technology, ShanghaiTech University. He is currently a state specially recruited Professor with the National Key Laboratory of Wireless Communications, University of Electronic Science and Technology of China. He has served on a number of technical programs for international conferences. His research interests include signal processing, machine learning, and wireless communications, such as intelligent communications, structured signal reconstruction, Bayesian approximate inference, and distributed learning. He has published over 200 peer-reviewed research papers in leading international journals and conferences in related areas. He was a co-recipient of the Best Paper Award of the IEEE International Conference on Communications (ICC) in 2014 and also a co-recipient of the Best Journal Paper Award of the IEEE Technical Committee on Green Communications and Computing (TCGCC) in 2017. He is an Editor of IEEE TRANSACTIONS ON WIRELESS COMMUNICATIONS.

<figure><img src="image_16.png" /></figure>

**Marco Di Renzo** (Fellow, IEEE) received the Laurea (cum laude) and Ph.D. degrees in electrical engineering from the University of L'Aquila, Italy, in 2003 and 2007, respectively, and the Habilitation à Diriger des Recherches (Doctor of Science) degree from University Paris-Sud (now Paris-Saclay University), Paris, France, in 2013. He was a France-Nokia Chair of Excellence in ICT with the University of Oulu, Finland, a Tan Chin Tuan Exchange Fellow in engineering with Nanyang Technological University, Singapore, a Fulbright

Fellow with The City University of New York, USA, a Nokia Foundation Visiting Professor with Aalto University, Finland, and a Royal Academy of Engineering Distinguished Visiting Fellow with Queen's University Belfast, U.K. Currently, he is a CNRS Research Director (Professor) and the Head of the Intelligent Physical Communications Group, Laboratory of Signals and Systems (L2S), CNRS, CentraleSupélec, Paris-Saclay University, and a Chair Professor in telecommunications engineering with the Centre for Telecommunications Research, Department of Engineering, King's College London, London, U.K. He is the Fellow of IET, EURASIP, and AAIA; an Academician of AIIA; an Ordinary Member of European Academy of Sciences and Arts, an Ordinary Member of the Academia Europaea; an Ambassador of European Association on Antennas and Propagation; and a Highly Cited Researcher. His recent research awards include the Michel Monpetit Prize conferred by the French Academy of Sciences and the IEEE Communications Society Marconi Prize Paper Award in wireless communications. His current main roles within the IEEE Communications Society include serving as a Voting Member of the Fellow Evaluation Standing Committee, the Chair of the Publications Misconduct Ad Hoc Committee, and the Director of Journals. He served as the Editor-in-Chief for IEEE COMMUNICATIONS LETTERS from 2019 to 2023.