{'fraction_non_alphanumeric': 0.03985956841499955, 'fraction_numerical': 0.006809721255409945, 'mean_word_length': 4.375305026842362, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Compressed CPD-Based Channel Estimation and Joint Beamforming for RIS-Assisted Millimeter Wave Communications

We consider the problem of channel estimation and joint active and passive beamforming for reconfigurable intelligent surface (RIS) assisted millimeter wave (mmWave) multiple-input multiple-output (MIMO) orthogonal frequency division multiplexing (OFDM) systems. We show that, with a well-designed frame-based training protocol, the received pilot signal can be organized into a low-rank third-order tensor that admits a canonical polyadic decomposition (CPD). Based on this observation, we propose two CPD-based methods for estimating the cascade channels associated with different subcarriers. The proposed methods exploit the intrinsic low-rankness of the CPD formulation, which is a result of the sparse scattering characteristics of mmWave channels, and thus have the potential to achieve a significant training overhead reduction. Specifically, our analysis shows that the proposed methods have a sample complexity that scales quadratically with the sparsity of the cascade channel. Also, by utilizing the singular value decomposition-like structure of the effective channel, this paper develops a joint active and passive beamforming method based on the estimated cascade channels. Simulation results show that the proposed CPD-based channel estimation methods attain mean square errors that are close to the Cramer-Rao bound (CRB) and present a clear advantage over the compressed sensing-based method. In addition, the proposed joint beamforming method can effectively utilize the estimated channel parameters to achieve superior beamforming performance.


I. INTRODUCTION
Millimeter wave (mmWave) and terahertz (THz) communications are able to support extremely high data rate transmissions [1], [2]. Nevertheless, due to the reduced diffraction effect and high penetration loss, mmWave/THz systems require more antennas and active nodes such as access points (APs) and relays to improve the signal coverage [3], [4]. Deploying active nodes incurs an additional energy consumption and meanwhile presents a serious network interference issue. It is thus of practical significance to develop innovative technologies to address the coverage issue of future mmWave/THz wireless networks with a low cost and complexity.
Recently, reconfigurable intelligent surface (RIS) has emerged as an energy-efficient and cost-effective solution to tackle the above challenges. Generally speaking, RIS intelligently adjusts the signal reflection through a large number of low-cost passive reflection elements, which can dynamically reshape the wireless propagation environment and thereby improve the system performance [5]- [8]. An important advantage of RIS is that it does not require any active circuits such as ratio frequency (RF) chains for signal transmission/reception, which reduces hardware complexity as well as energy consumptions compared to traditional active transceivers/relays. Furthermore, RIS can be easily attached to different objects (such as walls and ceilings), thus showing great flexibility and compatibility in practical deployment [9].
Channel state information (CSI) acquisition is a pre-requisite to realize the full potential of RIS-assisted mmWave systems. Nevertheless, since RIS is usually composed of a large number of passive elements, CSI acquisition for RIS-assisted mmWave systems faces the difficulty of requiring a large amount of training overhead. In [10]- [13], by utilizing the limited scattering nature of mmWave channels, fast beam training and alignment methods were proposed for RIS-assisted mmWave systems, where the objective is to simultaneously identify the best beam alignment for both the base station (BS)-RIS link and the RISuser link. Analyses and experimental results [13] show that such methods only require a modest amount of training overhead to establish an effective virtual line-of-sight (LOS) path for data transmission. Beam training methods, however, can only acquire the angular parameters associated with the dominant path, which prevents from utilizing the spatial diversity and achieving a higher spectral efficiency. Different from beam training methods, other studies, e.g. [14]- [23], aim to obtain the full CSI that suffices for active and passive optimization. Specifically, in order to reduce the training overhead, the inherent sparse structure of the cascade BS-RIS-user mmWave channel was exploited and the cascade channel estimation is cast into a compressed sensing framework [16]- [19]. In addition to compressed sensing-based methods, tensor decomposition-based methods were developed by exploiting the intrinsic multi-dimensional structure of the channel [20]- [23]. Most of these tensor decomposition-based methods [20]- [22], however, did not utilize the sparse scattering characteristics of mmWave channels and have a CP rank that is equal to the number of reflecting elements. As a result, these methods require a training overhead proportional to the number of reflecting elements, which is usually large in practice.
In this paper, we study the problem of channel estimation and joint beamforming for RIS-assisted mmWave MIMO-OFDM systems. We show that by exploring the sparse scattering characteristics and the intrinsic multi-dimensional structure of the mmWave cascade channel, the received signal can be formulated into a low-rank tensor that admits a canonical polyadic decomposition (CPD). Based on this formulation, an alternating least squares (ALS) method and a Vandemonde structured-based method are developed for channel estimation. Theoretical analysis shows that the proposed methods have a sample complexity of O(U 2 ). Here U denotes the sparsity of the cascade channel. Since U is usually small relative to the dimension of the cascade channel, the proposed methods can achieve a substantial training overhead reduction. The proposed methods, unlike compressed sensing-based techniques, are essentially a gridless approach and therefore are free of the grid discretization errors.
In addition to channel estimation, this paper also considers the problem of joint active and passive beamforming design. For point-to-point RIS-assisted MIMO systems, existing beamforming methods, e.g. [24], [25], usually require the global CSI knowledge, i.e. the channels of both the BS-RIS link and the RIS-user link, for joint beamforming. In this paper, by exploring the inherent structure of the effective channel, we show that the knowledge of the cascade channel alone suffices for the purpose of joint active and passive beamforming, and we develop a manifold optimization-based scheme for active and passive optimization based on the estimated cascade channel.
The current work is an extension of our previous work [23] which developed a CPD-based channel estimation method for RIS-assisted mmWave MISO systems. The extension of the current work consists of two aspects. First, we extend the CPD-based method to the MIMO scenarios. Note that such an extension is nontrivial. In fact, for the MISO case, as pointed out in [23], the Kruskal's condition which is essential to the uniqueness of the CPD does not hold and hence the classical ALS method cannot be applied. We show that due to the diversity brought by multi-antenna at the receiver, the Kruskal's condition can be satisfied for MIMO scenarios. Such a fact enables us to develop an ALS-based channel estimation method that is more robust against noise. Second, besides channel estimation, the current work also considers how to optimize the active and passive beamforming coefficients based on the estimated cascade channel, which is a challenging problem (particularly for MIMO systems) and was not studied in our previous work.
The rest of the paper is organized as follows. In Section II, the system model and the formulation of the channel estimation problem are discussed. CPD-based channel estimation methods are developed in Section III. In Section IV, the problem of joint active and passive beamforming design is studied. Finally, simulation results are presented in Section V.

II. SYSTEM MODEL AND PROBLEM FORMULATION
Consider a point-to-point RIS-aided mmWave MIMO-OFDM system, where an RIS is deployed to assist the downlink data transmission of N s data streams from the BS to the user (see Fig. 1). For simplicity, we assume that the direct link between the BS and the user is blocked due to poor propagation conditions. The BS is equipped with a uniform linear array (ULA) with N t antennas and R t radio frequency (RF) chains, and the user is equipped with a ULA with N r antennas and R r RF chains,

A. Channel Model
In this paper, we adopt a geometric wideband mmWave channel model [26] to characterize the channel between the BS (RIS) and the IRS (user). Specifically, the BS-RIS channel in the delay domain can be expressed as where L is the total number of paths between the BS and the RIS, α l is the complex gain associated with the lth path, φ l represents the spatial angle of departure (AoD), {ϑ r l , χ r l } denote the spatial azimuth and elevation angles of arrival (AoAs), τ l denotes the time delay, δ (τ ) denotes the Dirac-delta function, a IRS (ϑ, χ) and a BS (φ) denote the normalized array response vectors associated with the IRS and the BS, respectively. Similarly, the IRS-user channel in the delay domain is modeled as where L r is the number of paths between the RIS and the user, β l denotes the associated complex path gain, θ l represents the spatial AoA, {ϑ t l , χ t l } denote the spatial azimuth and elevation AoDs, κ l is the time delay, and a UE (θ) denotes the normalized array response vectors associated with the user. The normalized array response vectors a BS (φ), a UE (θ) and a IRS (ϑ, χ) are respectively defined as where ⊗ denotes the Kronecker product, φ Accordingly, the frequency-domain BS-RIS and RIS-user channels associated with the pth subcarrier can be respectively expressed as where f s is the sample frequency, and P 0 denotes the total number of OFDM tones.

B. Downlink Training and Signal Model
To facilitate the algorithm development, we employ a frame-based downlink training protocol (see where s t,p ∈ C Ns denotes the pth subcarrier's pilot symbol vector, F BB,t,p ∈ C Rt×Ns is the baseband precoding matrix associated with the pth subcarrier, and F RF,t ∈ C Nt×Rt is a radio frequency (RF) precoder common to all subcarriers. The transmitted signal arrives at the user via propagation through the BS-RIS-user channel. The user employs the hybrid combiner W t,p = W RF,t W BB,t,p ∈ C Nr×Ns to combine the received signal, where W BB,t,p ∈ C Rr×Ns denotes the baseband combiner associated with the pth subcarrier and W RF,t ∈ C Nr×Rr is a RF combiner common to all subcarriers. Hence, the received signal y q,t,p ∈ C Ns associated with the pth sub-carrier at the qth time slot of the tth time frame can be written as In the channel estimation stage, for simplicity, we assume F BB,t,p = F BB,t and s t,p = s t , in which case we have f t,p = f t F RF,t F BB,t s t . Similarly, let W BB,t,p = W BB,t . We have W t,p = W t W RF,t W BB,t . The received signal can thus be expressed as where we define X t f t ⊗ W * t ∈ C NtNr×Ns , ⊙ denotes the Khatri-Rao product, and H p G T p ⊙ R p ∈ C NtNr×M denotes the cascade channel associated with the pth subcarrier. Recalling (8) and (9), the cascade channel H p can be further expressed as where a S (φ m , θ n ) a * BS (φ m ) ⊗ a UE (θ n ), and the mapping process (a) is defined as where ⌈x⌉ denotes the ceiling function which gives the least integer greater than or equal to x.
Our objective of this paper is two-fold. First, we wish to develop a method to estimate the cascade channel matrices {H p } from the received measurements {y q,t,p }. In particular, the proposed method is expected to provide a reliable channel estimate by using as few pilot symbols as possible. After the cascade channel matrices are obtained, another purpose of this work is to develop a joint active and passive beamforming method which aims to maximize the spectral efficiency by exploiting the knowledge of the cascade channels.

A. Low-Rank Tensor Representation
Substituting (13) into (12), we have Define Y t,p y 1,t,p · · · y Q,t,p T ∈ C Q×Ns . The received signal at the tth time frame can be written where Collecting the received signals from all time frames and defining Y p Since signals associated with multiple subcarriers are available, the entire received signal {Y p } p can be expressed as a third-order tensor Y ∈ C Q×T Ns×P whose three modes respectively represent the subframe, the time frame and the subcarrier. Each slice of the tensor Y, Y p , is a weighted sum of a common set of rank-one outer products. Hence the third-order tensor Y admits a CPD which decomposes a tensor into a sum of rank-one component tensors, i.e.
where U LL r , N ∈ C Q×T Ns×P is the noise tensor, and Define Here {A, B, C} are the factor matrices of the tensor Y. We see that the channel parameters can be readily estimated from these factor matrices. Therefore our objective is to first obtain the factor matrices from the observed tensor Y, and then estimate the associated channel parameters from the estimated factor matrices.

B. Uniqueness Condition
Before proceeding to the CPD, we first study the condition that ensures the uniqueness of the CPD since the uniqueness of the CPD is essential to the success of our proposed method. This condition also sheds light on the sample complexity of the proposed method, i.e. the amount of training overhead required to reliably estimate the channel.
A well-known sufficient condition for the uniqueness of the CPD is known as Kruskal's condition summarized as follows [27].
Theorem 1: Let X ∈ C I×J×K be a third-order tensor decomposed of three factor matrices A (1) ∈ is satisfied, then the CPD of X is unique up to scaling and permutation ambiguities.
Here k A denotes the k-rank of A, which is defined as the largest value of k A such that every subset of k A columns of A is linearly independent. From the above theorem, we know that if then the CP decomposition of Y is essentially unique.
We first examine the k-rank of A. Recall that Assume that each entry of V ∈ C M ×Q is randomly generated from the unit circle, i.e. v m,n = e jσm,n , Recall that a IRS (ζ i , ξ i ) is a Kronecker product of two steering vectors. Hence we have a H [28]. In reality, due to the random nature of the channel parameters, the angles {ζ u } U u=1 are mutually distinct with probability one, and so are the angles {ξ u } U u=1 . We thus have E a H m,i a n,j ≈ 0 even for the case m = n. On the other hand, according to the central limit theorem, a m,i approximately follows a Gaussian distribution. Therefore entries of A can be considered as i.i.d.
Gaussian variables with zero mean and unit variance. As a result, we have The factor matrix B can be written as Similarly, for a randomly generated F whose entries are uniformly chosen from a unit circle, we can arrive at the conclusion that the k-rank of B is equal to Note that due to the mapping process, the set {φ u } only contains L distinct elements, and {θ u } only pair of (φ u , θ u ) is unique according to the mapping rule defined in (14).
As for the factor matrix C, it is a Vandermonde matrix with distinct generators, i.e. ι i = ι j , ∀i = j.
Thus we have Based on the above results, we know that Kruskal's condition is equivalent to Note that the total number of pilot signals for downlink training is QT P . To meet condition (35), we can set Q ≥ U , T N s ≥ U and P ≥ 2, in which case the amount of training overhead is in the order of We see that the sample complexity of the proposed method only depends on the sparsity of the cascade channel U . As U is usually small relative to the dimension of the cascade channel, a substantial training overhead reduction can be achieved.

C. CPD: An ALS-Based Approach
For generic CPD problems, an alternating least squares (ALS) method is usually employed to search for the factor matrices. Specifically, assume that the CP rank, U , is known a priori. The CP decomposition of Y can be accomplished by solving: The above optimization can be solved by an ALS procedure which alternatively solves the following least squares problemŝ where Y (n) denotes the mode-n unfolding of the tensor Y.
If the CP rank U is unknown, a sparsity-promoting regularizer can be added to the objective function such that one can automatically determine the CP rank to find a low-rank representation of the observed tensor. LetÛ > U denote an overestimated CP rank. The problem can be formulated as [29], [30] min where µ is a regularization parameter to control the tradeoff between low-rankness and the data fitting Again, the above optimization problem can be solved via the ALS method.

D. CPD: A Vandermonde Structured-Based Approach
The ALS-based CPD method involves a high computational complexity as it needs to solve several large-scale least squares problems at each iteration. Notice that, for our CPD problem, one of the factor matrices is a Vandermonde matrix. Such a Vandermonde structure along with its linear algebra properties can be utilized to devise an efficient CPD method which yields a closed-form solution of the factor matrices [31], [32]. Details of the Vandermonde structured-inspired CPD method are provided below.
To utilize the Vandermonde structure of the factor matrix C, we consider the mode-1 unfolding of the received tensor Y, which is a matrix constructed by concatenating the mode-1 fibers of Y and can be expressed as: We perform the truncated singular value decomposition (SVD) Y T To facilitate the exposition, we ignore the noise. Since (C ⊙ B) is full column rank, we know that there exists a nonsingular matrix M ∈ C U ×U such that Consequently, we have where A denotes a submatrix of A obtained by removing the top row of A, A denotes a submatrix of A obtained by removing the bottom row of A, and On the other hand, by utilizing the Vandermonde structure of C, we have where Z diag (z 1 , · · · , z U ), and z u e −j2π fs P 0 ιu is the generator of the factor matrix C. Hence we arrive at Since C ⊙ B is full column rank, both U 1 and U 2 are full column rank. Therefore (48) can be further rewritten as The above equation implies that the generators {z u } U u=1 and M can be estimated from the eigenvalue decomposition (EVD) of U † 1 U 2 . Based on the estimated generators {ẑ u }, each column of the factor matrix C can be estimated asĉ u = ẑ uẑ Finally, givenB andĈ, the factor matrix A can be estimated aŝ

E. Channel Estimation
After obtainingÂ,B andĈ, we now proceed to estimate the channel parameters. From the CPD theory, we know that the estimated {Â,B,Ĉ} and the true factor matrices {A, B, C} are related aŝ where are estimation errors, and Γ is an unknown permutation matrix. This permutation matrix Γ is common to all factor matrices, and thus can be ignored.
Recall that each column of the factor matrix A is characterized by spatial angle parameters {ζ u , ξ u }.
Therefore these two spatial angle parameters can be estimated through a correlation-based estimator: whereâ u denotes the uth column ofÂ.
Similarly, the spatial AoD associated with the BS φ u and the spatial AoA associated with the user θ u can be estimated as whereb u denotes the uth column ofB.
Note that the uth column of C is given by g (ι u ). Therefore the time delay ι u can be estimated viâ whereĉ u denotes the uth column ofĈ.

IV. JOINT ACTIVE AND PASSIVE BEAMFORMING DESIGN
In this section, we consider the problem of joint active and passive beamforming design based on the estimated channel parameters. Specifically, we aim to maximize the spectral efficiency by jointly optimizing the precoding matrices at the transmitter, the combining matrices at the receiver, and the reflection coefficients at the RIS. In the data transmission stage, the received signal associated with the pth subcarrier can be expressed as whereH p R p ΦG p ∈ C Nr×Nt is the equivalent (i.e. effective) channel associated with the pth subcarrier, and n p ∼ CN (0, σ 2 I) denotes the additive white Gaussian noise.
Assuming that the transmitted signal obeys a Gaussian distribution, the achievable spectral efficiency can be calculated as where † denotes the Moore-Penrose inverse.

A. Problem Formulation
With the objective of maximizing the spectral efficiency, the joint beamforming problem can be formulated as where S F RF , {F BB,p } P p=1 , W RF {W BB,p } P p=1 , Φ is the set of optimization variables, and ρ is the transmit power constraint. To simplify the problem, we first ignore the constraint introduced by the hybrid analog/digital structure and consider a fully digital precoder/combiner. Let F p ∈ C Nt×Ns and W p ∈ C Nr×Ns denote a fully digital precoder and a fully digital combiner, respectively.
Once an optimal fully digital precoder/combiner is found, we can use the manifold optimization method [33] to search for a hybrid precoder/combiner to approximate the optimal fully digital precoder/combiner.
Due to the sparse scattering nature of mmWave channels, such a strategy has been proven effective in previous studies, e.g. [34], [35], which showed that hybrid beamforming/combining with a small number of RF chains can asymptotically approach the performance of fully digital beamforming/combining.
Given the reflection matrix Φ, we first study how to devise the fully digital precoder/combiner. Let r p = rank H p . Define the truncated SVD of the effective channelH p as where U p ∈ C Nr×rp , V p ∈ C Nt×rp , and Σ p is an r p ×r p diagonal matrix. Also, we have U 1,p ∈ C Nr×Ns , Σ 1,p ∈ C Ns×Ns and V 1,p ∈ C Nt×Ns . Given a specified Φ, the optimal fully digital precoder/combiner is given as where and denotes the optimal amount of power allocated to the ith data stream, 1/λ p is the water level satisfying Ns i=1 ρ i,p = ρ. Thanks to the massive array gain provided by the RIS, the effective signal-to-noise ratio (SNR) is usually large, in which case an equal power allocation scheme is near-optimal. Therefore we can approximate F opt,p as: Substituting the optimal precoder/combiner {F opt,p } and {W opt,p } into (66), we arrive at a problem which concerns only the optimization of the passive beamforming matrix Φ: (72)

B. Passive Beamforming Design
We assume |α 1 | ≥ |α 2 | ≥ . . . ≥ |α L |, |β 1 | ≥ |β 2 | ≥ . . . ≥ |β Lr |. To gain insight into the passive beamforming design, we write the equivalent channel as follows: where we define D p (m, n) α m β n e −j2πfs(τm+κn) p P 0 d mn , and In the above equation, α m β n is the complex gain associated with the (m, n)th BS-RIS-user composite path which is composed of the mth path from the BS to the RIS and the nth path from the RIS to the user, {ϑ r m − ϑ t n , χ r m − χ t n } are the composite RIS angles associated with the (m, n)th composite path, and d mn is referred to as the passive beamforming gain associated with the (m, n)th composite path.
Clearly, there are a total number of U = LL r composite paths.
When the numbers of antennas N t and N r are sufficiently large, A UE and A BS can be considered as orthonormal matrices with unit-norm and mutually orthogonal columns. If the phase shift vector v is properly devised such that the off-diagonal elements of D p are small relative to entries on the main diagonal, thenH p = A UE D p A H BS can be approximated as a truncated SVD ofH p , in which case the optimization problem (72) turns into where δ is a small positive value and the constraint |d ij | < δ is imposed to ensure thatH p = A UE D p A H BS is a good approximation of the truncated SVD ofH p . As analyzed in [36], due to the asymptotic orthogonality of RIS's array response vectors characterized with different angular parameters, the constraint |d ij | < δ can be neglected and the solution to the simplified problem can automatically guarantee that off-diagonal entries of D p are small relative to entries on its main diagonal. Specifically, by ignoring the constraint |d ij | < δ, the optimization (77) can be simplified as Such an optimization can be efficiently solved via a manifold optimization technique, whose details can be found in [36].
From (78), we see that to optimize the reflection coefficients, we only need the knowledge of the composite gains {α i β i } Ns i=1 and the composite RIS angles Recall that, in the channel estimation stage, the following channel parameters {ζ u ,ξ u ,φ u ,θ u ,ι u ,̺ u } U u=1 are obtained, in which we have We see that our proposed channel estimator can provide an estimate of the composite gains as well as the composite RIS angles associated with all U composite paths. The problem now is how to appropriately choose N s composite paths from these U composite paths. Randomly choosing N s composite paths certainly does not work. In fact, from (73), it is easy to know that the composite paths corresponding to the diagonal entries of D p must have mutually distinct AoDs at the BS and mutually distinct AoAs at the user. Also, to improve the spectral efficiency, clearly we should choose those composite paths whose composite gains are as large as possible. Based on the above considerations, the N s composite paths can be selected based on the following criterion: where the last two constraints are imposed to ensure that the selected composite paths have mutually distinct AoDs at the BS and mutually distinct AoAs at the user, in which δ BS and δ UE are small positive parameters of user's choice. Based on the selected N s composite paths, the optimization problem (78) can be further written as max v i∈I The above optimization problem can be efficiently solved via the manifold optimization method proposed in [33].

C. Active Beamforming Design
After the passive beamforming vector v is determined, according to (73), the pth subcarrier's equivalent channelH p can be estimated aŝ where the second equality is a truncated SVD ofĤ p . Based on the previous discussion, the optimal fully digital precoder/combiner can be obtained as After the optimal fully digital precoder/combiner is obtained, we search for a common analog precoding (combining) matrix F RF (W RF ) and a set of baseband precoding (combining) matrices {F BB,p } ({W BB,p }) to approximate the optimal precoder (combiner) {F opt,p } ({W opt,p }). The problem can be formulated as The above optimization problem can be solved via the manifold optimization technique introduced in [37].

V. SIMULATION RESULTS
We present simulation results to evaluate the performance of the proposed CPD-based channel estimation methods and the joint beamforming scheme. In Section III, two different approaches are introduced to perform the CPD, namely, the ALS method and the Vandemonde structure-based method. The corresponding channel estimation methods are referred to as CPD-ALS and CPD-VS, respectively.
In our simulations, we assume that the BS employs a ULA with N t = 32 antennas and R t = 4 RF chains, the IRS is equipped with M = 16 × 16 passive reflecting elements, and the user employs a ULA with N r = 32 and R r = 4 RF chains. The angular parameters {ϑ r l , χ r l , φ l } L l=1 , ϑ t l , χ t l , θ l Lr l=1 are randomly generated from [0, 2π], where we set L = 3 and L r = 3. The delay spreads {τ l } L l=1 , {κ l } Lr l=1 are drawn from a uniform distribution U (0, 100ns). The complex gains {α l } L l=1 and {β l } Lr l=1 follow a circularly symmetric Gaussian distribution CN (0, 1). The number of data streams is set to N s = 2. The total number of subcarriers is set to P 0 = 128, among which P subcarriers are used for training. The sampling rate is set to f s = 0.32GHz. The signal-to-noise ratio (SNR) is defined as We first examine the estimation accuracy of the channel parameters {ζ u , ξ u , φ u , θ u , ι u , ̺ u } U u=1 and the overall channel estimation performance. Note that the channel estimation problem being considered in this work can be cast as a multi-measurement vector (MMV) compressed sensing problem, and the simultaneous-OMP method (SOMP) [38] can be used to estimate the cascade channel. For the SOMP  an estimation accuracy close to the theoretical lower bound. Also, the CPD-ALS method is superior to the CPD-VS method in the low SNR regime. This is probably because the CPD-VS which depends critically on the structure of the factor matrix is more sensitive to noise. Also, it can be observed that both CPD-based methods present a substantial advantage over the SOMP method.
In Fig. 4(a), we plot the estimation performance of respective methods as a function of the SNR. The performance is evaluated via the normalized mean squared error (NMSE) of the cascaded channel, which is defined as P p=1 F . Again, we see that our proposed methods present a significant performance improvement over the SOMP method. In particular, the CPD-ALS method achieves a decent estimation performance even in a low SNR regime, say SNR = 0dB. Note that in mmWave communications, due to the severe path loss, the SNR for channel estimation is usually low, with the range of interest from 0dB to 10dB at most. Thus the ability of delivering an accurate channel estimate in the low SNR regime is highly desirable and has important practical implications.
In Fig. 4(b), we plot the estimation performance of respective methods as a function of the number of subcarriers P , where we set Q = 10 and T = 5. It can be seen that the proposed methods provide a reliable channel estimate when P ≥ 8, which corresponds to a total number of 400 measurements for training. As a comparison, note that the cascade channel H p to be estimated has a size of N t N r × M = 1024 × 256, which has more than 2.6 × 10 5 parameters. This result indicates that the proposed methods can achieve a substantial training overhead reduction. Fig. 4(c) plots the estimation performance versus the number of time slots T , where we set Q = 10 and P = 10. This result, again, demonstrates the superiority of the proposed methods over the compressed sensing-based method.
Next, we examine the beamforming performance attained by the joint beamforming scheme proposed in Section IV. To illustrate the effectiveness of the proposed channel estimator, we include the beamforming performance attained by assuming the perfect knowledge of the CSI, which serves as an upper bound on the beamforming performance attained by using the estimated CSI. Fig. 5 plots the spectral efficiency of the proposed joint beamforming scheme as a function of SNR and the number of time frames, respectively. We see that our proposed CPD-based estimators achieve performance close to that attained by assuming perfect CSI knowledge even in the low SNR regime, which verifies the effectiveness of the proposed estimation method. Also, our proposed methods present a clear performance improvement over the compressed sensing-based method. Particularly, when the training overhead is low, say, P = Q = 16 and T = 2, the proposed CPD-based methods can still achieve decent beamforming performance, whereas the compressed sensing-based method incurs a significant performance loss.

VI. CONCLUSIONS
In this paper, by exploiting the intrinsic multi-dimensional structure as well as the sparse scattering characteristics of the mmWave channels, we developed two CPD-aided channel estimation methods, namely, an ALS-based CPD method and a Vandemonde structure-based CPD method, for RIS-assisted mmWave MIMO-OFDM systems. The proposed methods effectively utilize the low-rankness of the CPD formulation and can achieve a substantial training overhead reduction. We also developed a joint beamforming scheme that utilizes the estimated cascade channel parameters for optimizing the system's active and passive variables. Simulation results show that our proposed methods present a significant performance advantage over the compressed sensing method, and can achieve superior channel estimation and beamforming performance with a low training overhead.

A. Partial
in whichã a,u jV T D 1 a IRS (ζ u , ξ u ) and where Re {·} represents the real part of a complex number, e u is a unit vector whose uth entry equals to one and all other entries equal to zeros, andÃ a [ã a,1 · · ·ã a,U ]. Similarly, we have whereÃ e [ã e,1 · · ·ã e,U ],ã e,u = jV T D 2 a IRS (ζ u , ξ u ), and D 2 diag(0, · · · , M z − 1

B. Calculation of Fisher Information Matrix Ω (p)
We first calculate the entries in the principal minors of Ω (p).
where N (1) is the mode-1 unfolding of N , vec(N H (1) ) ∼ CN (0, σ 2 I). Since w 1 is a linear transformation of vec(N H (1) ), it also follows a circularly symmetric complex Gaussian distribution. Its covariance matrix C w1 ∈ C U 2 ×U 2 and second-order moments M w1 ∈ C U 2 ×U 2 are respectively given by where C w2 = 1 σ 2 (Ã where Other entries in the off-principal minors of Ω (p) can be similarly calculated.

C. Cramér-Rao Bound
After obtaining the FIM, the CRB for the parameters p can be calculated as