Structured Channel Estimation for RIS-Assisted THz Communications

IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 74, NO. 3, MARCH 2025 5175 Structured Channel Estimation for RIS-Assisted THz Communications Fazal-E-Asim , Senior Member, IEEE, Bruno Sokal , Member, IEEE,AndréL.F.de Almeida , Senior Member, IEEE, Behrooz Makki , Senior Member, IEEE, and Gábor Fodor , Senior Member, IEEE Abstract—This paper proposes tensor-based channel estimation for reconfigurable intelligent surface (RIS)-assisted communication networks. We exploit the inherent geometrical structure of the Terahertz propagation channel, including the antenna array geometries at the base station, the RIS, and the user equipment to design a tensor-based channel estimator, referred to as the higher-dimensional rank-one approximations (HDR) method. By exploiting the geometrical structure of the combined base station-RISuser equipment channel, the proposed HDR estimator recasts parametric channel estimation as a single sixth-order rank-one tensor approximation problem, which can be efficiently solved using higher-order singular value decomposition to deliver parallel estimates of each channel component vector. Numerical results show that the proposed method provides significantly more accurate parameter estimates than competing state-of-the-art tensorbased RIS channel estimation, Khatri-Rao factorization, and least squares methods. For higher-rank channels, the HDR method shows similar spectral efficiency compared to its competitors while having similar computational complexity to the classical least squares estimator. Index Terms—Rank-one approximation, Terahertz (THz) communications, channel parameter estimation, and tensor modeling. I. INTRODUCTION Reconfigurable intelligent surface (RIS) is a promising candidate for the enhancement of future communication networks. Moreover, the introduction of millimeter-wave (mmWave) and Terahertz (THz)/(micrometer) bands facilitates the deployment of a massive number of antennas at the base station (BS) and the user equipment (UE). To deal with these problems, RISs that help to concentrate energy towards the desired location and minimize the interference elsewhere is a promising solution to boost the signal-to-interference-plus-noise ratio in an energy-efficient fashion. On the other hand, estimating Received 7 August 2024; accepted 2 November 2024. Date of publication 7 November 2024; date of current version 5 March 2025. The work of Gábor Fodor was supported by the EU Horizon Europe Program under Grant 101139176 6G-MUSICAL. This work was supported in part by Ericsson Research, Sweden, in part by Ericsson Innovation Center, Brazil, under Grant UFC.51, in part by the Fundação Cearense de Apoio ao Desenvolvimento Científico e Tecnológico (FUNCAP) under Grant FC3-00198-00056.01.00/22 and Grant ITR-0214-00041.01.00/23, in part by the National Institute of Science and Technology (INCT-Signals) through Brazil’s National Council for Scientific and Technological Development (CNPq) under Grant 406517/2022-3, and in part by CNPq under Grant 312491/2020-4 and Grant 443272/2023-9. The review of this article was coordinated by Prof. Mugen Peng. (Corresponding author: Fazal-E-Asim.) Fazal-E-Asim, Bruno Sokal, and André L. F. de Almeida are with the Wireless Telecommunications Research Group (GTEL), Department of Teleinformatics Engineering, Federal University of Ceara, Fortaleza 60455-970, Brazil (e-mail: [email protected]; [email protected]; [email protected]). Behrooz Makki is with the Ericsson Research, Ericsson, 417 56 Gothenburg, Sweden (e-mail: [email protected]). Gábor Fodor is with the Ericsson Research, 16480 Stockholm, Sweden, and also with the Division of Decision and Control, KTH Royal Institute of Technology, 11428 Stockholm, Sweden (e-mail: gabor[email protected]). Digital Object Identifier 10.1109/TVT.2024.3492998 wireless channels is challenging when dealing with large RIS panels in combination with massive multiple input multiple output (MIMO) infrastructure nodes. Acquiring channel state information (CSI) is especially challenging when passive RISs are deployed, where the burden of pilot transmission, reception, and associated signal processing must be managed by the end nodes of the communication network [1]. Furthermore, with the use of higher frequencies, the communication channel is dominated by a few usable propagation paths, which are more susceptible to blockages that decrease the coverage, throughput, and quality of services provided by the network. Here, RISs can effectively avoid the blockage effects by establishing a virtual path. To summarize, a passive RIS cannot process pilot signals, which indirectly means that channel estimation/parameter estimation must be performed at the end nodes of the wireless network, potentially overloading the chosen node. Secondly, a large panel size usually implies many channel coefficients to be estimated. This may require many pilot resources, affecting the overall system spectral efficiency (SE). Finally, the accuracy of the CSI acquisition is crucial because the gains the RIS provides depend on a sufficiently accurate knowledge of the channel. The paper [2] briefly explains various channel estimation methods in RIS-assisted communications, considering both unstructured and structured channel models. Algorithms introduced for estimating unstructured channels are simple and easy to implement but require extensively huge training overhead in RIS-based systems, compromising their achievable rates due to the estimation of many channel coefficients. On the other hand, geometric channel models lead to an estimation of fewer channel parameters, and therefore, much smaller training is required for improved estimation performance. These benefits, however, come at the cost of increased algorithmic complexity, model order estimation, and inevitable modeling errors. The authors in [3] use the inductive matrix completion followed by a root multiple signal classification (MUSIC) algorithm to estimate the respective angle of departure and arrival in RIS networks. However, the method is only shown for a uniform linear array (ULA), including RIS, with no clear method for angle pairing. Also, [4] exploits the sparsity of the mmWave propagation channel to jointly design beamforming and estimating the cascaded BS-RIS-UE channel. However, the proposed method is designed for a single-antenna user with a ULA BS. The work in [5] designs a joint channel estimation and data detection algorithm for hybrid reconfigurable intelligent surface (HRIS) using orthogonal time frequency space modulation. A new transmission structure is formulated so that the partial HRIS elements are alternatively activated. The channel and unknown data symbols are estimated by iteratively executing the message passing and expectation-maximization algorithms. Therein, a ULA BS and a single-antenna user are assumed, which impedes full exploitation of the channel structure. The works in [6], [7] used distributed semi-passive IRS for sensing and passive IRS for communications in an integrated sensing and communication (ISAC) system for single and multi-user cases. The total least squares (LS) estimation of signal parameters via the rotational invariance technique is used to estimate the angles of arrival effectively. Then, as a second step, the MUSIC algorithm is used to obtain effective pairing. However, single antenna users are considered in these works. Tensor decompositions have been successfully applied to formulate channel estimation methods for mmWave/massive MIMO systems in different scenarios. In particular, the widespread use of the parallel factor (PARAFAC) tensor decomposition in wireless communications 0018-9545 © 2024 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See https://www.ieee.org/publications/rights/index.html for more information. Authorized licensed use limited to: TEL AVIV UNIVERSITY. Downloaded on August 10,2025 at 10:58:26 UTC from IEEE Xplore. Restrictions apply. 5176 IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 74, NO. 3, MARCH 2025 comes from its capability to exploit the intrinsic multilinear structure of signals/channels and its powerful uniqueness properties [8]. The authors in [9] combine tensor decomposition and compressive sensing to formulate a sparse channel estimation method for massive MIMO-OFDM (orthogonal frequency division multiplexing) systems. The work [10] proposes a joint hybrid precoder and combiner design for maximizing the achievable sum rate of mmWave MIMO-OFDM systems, where the analog precoder and combiner are designed as a Tucker2 tensor decomposition problem. Tensors have recently been exploited to derive channel estimation methods for RIS-assisted communications [11],[12],[13],[14],[15], [16]. The authors in [11] capitalize on the PARAFAC decomposition to formulate an iterative algorithm based on alternating least squares (ALS) [8] to solve the individual channel estimation problem for RIS-assisted MIMO communications. The work [17] jointly exploits the sparsity and tensor decomposition structures to improve channel estimation performance for RIS in a multiuser scenario. The work [13] introduces sparsity-structured-based PARAFAC tensor factorization. A two-stage sparse recovery problem is proposed for mmWave channels, referred to as two-stage RIS-aided channel estimation (TRICE). Therein, the first stage provides the directions of departure and arrival estimates, while in the second stage, the cascaded channel is estimated from the previously extracted angular parameters. The work [14] exploits low-rank channels and designs a PARAFAC decomposition, i.e., tensor-based RIS channel estimation (TenRICE), where the channel matrices are estimated using an ALS method. However, the latter TenRICE method outperforms the TRICE method. In [15], tensor-based channel estimation schemes are proposed for IRS-MIMO systems considering hardware imperfections affecting the phase shifts. The work [16] formulates a tensor-based approach to channel estimation and data tracking in RIS-assisted MIMO systems. The reader is referred to [18] and [19] for overviews on tensor decompositions applied to wireless communications. Reference [20] proposes PARAFAC-based unified tensor modeling for wireless communication and analyses various applications subject to frequency-selective multipath fading, such as multiuser equalization. Paper [21] provides an overview of constrained PARAFAC models having linear dependencies among columns of the factor matrices of the tensor decomposition. The work in [11] jointly exploits the sparsity and tensor decomposition structures to improve channel estimation performance for RIS-assisted communication in a multiuser scenario having a single antenna each. That work introduces sparsity-structured-based PARAFAC tensor factorization. Finally, [22] proposes an adaptive grid-matching pursuit estimation algorithm by transforming the estimation problem into a sparse recovery one, while in [23] a variational inference-sparse Bayesian learning estimator is proposed for RIS-assisted networks. Most of the previous works either use single antenna users or deploy ULA both at the UE and the BS. Please note that although the proposed work in this paper, assumes a uniform rectangular array (URA) at both the BS and the UE, the proposed method still applies to ULA deployed at both ends. In this paper, we exploit the geometrical structure of the channel estimation problem and recast it as a rank-one higher-order tensor approximation problem to estimate the channel parameters associated with the dominant communication links. More specifically, we reformulate the BS-RIS-UE channel as an approximated sixth-order rank-one channel tensor, which allows us to recast parametric channel estimation as a single rank-one tensor approximation problem, which can be efficiently solved using the higher-order singular value decomposition (HOSVD) algorithm [24]. Our proposed higher-dimensional rank-one approximations (HDR) method offers remarkable gains over its competitors in terms of normalized mean square error (NMSE) in the low signal-to-noise-ratio (SNR) regime. Moreover, our proposed HDR Fig. 1. System model. method outperforms the Cramér-Rao lower bound (CRLB) derived in [11], due to its efficient noise rejection by exploiting the rank-one tensor channel structure. The contributions of the present paper are summarized as follows: rExploiting the geometrical structure of the propagation channels, we resort to a tensor modeling formalism to formulate the received pilot signal as a higher-order tensor. rUsing a data reshuffling scheme, we recast the structured channel estimation problem as an HDR problem to estimate the channel parameters associated with the dominant communication links. rAnalyzing the performance of the proposed HDR method for several scenarios and comparing it with the classical LS method [11] and the theoretical CRLB [11] for unstructured channels. rComparing our solution with state-of-the-art competing methods such as Khatri-Rao factorization (KRF) [11] and TenRICE [14]. Notation: Scalars are denoted by lower-case italic letters a, vectors by bold lower-case italic letters a, matrices by bold upper-case italic letters a, tensors are defined by calligraphic upper-case letters A. AT,A∗,AHstand for transpose, conjugate and Hermitian of a.The operators ⊗,,◦,anddefine the Kronecker, the Khatri-Rao, the outer product, and the Hadamard (element-wise) product, respectively. diag{a}represents the square diagonal matrix with elements of vector aacross the main diagonal. vec{.}vectorizes an I×Jmatrix argument, while unvec(I×J){.}does the opposite operation. E[·]is expectation operator. For an Nth order tensor Y∈CL1×···×LN,then-mode unfolding of Yis the matrix [Y](n)=CLn×L1...Ln−1Ln+1...LN. II. SYSTEM AND CHANNEL MODEL Let us consider RIS-assisted MIMO communications for a THz scenario, in which we consider the deployment of a URA [25] placed in the y-z plane installed at the BS having M=MyMznumber of transmit antennas, where Myis the number of antenna elements along the y-axis and Mzis the number of antenna elements along the z-axis. The direct path between the BS and the UE is unavailable due to blockage or deep fade. Similarly, a URA is deployed at the UE having Q=QyQznumber of receive antenna elements, where Qyand Qz denote the antenna elements along the y and z axes, respectively. Moreover, N=NyNzis the number of reflecting elements deployed at the RIS, with Nyand Nzbeing the number of reflecting elements along the y and z axes, respectively (see Fig. 1). The geometrical channel between the BS and the RIS is represented as, H= R  r=1 αrb(φr risA,θr risA)aT(φr bs,θr bs)∈CN×M,(1) where αris the rth complex path gain, a(φr bs,θr bs)∈CM×1is the rth two-dimensional channel steering vector at the BS, φr bs is the rth azimuth of departure (AoD) and θr bs as rth elevation of departure (EoD). Similarly, b(φr risA,θr risA)∈CN×1is the rth two-dimensional channel steering vector at the RIS, with φr risAbeing the azimuth of arrival (AoA) and θr risAthe elevation of arrival (EoA). Analogously, the channel Authorized licensed use limited to: TEL AVIV UNIVERSITY. Downloaded on August 10,2025 at 10:58:26 UTC from IEEE Xplore. Restrictions apply. IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 74, NO. 3, MARCH 2025 5177 between the RIS and the UE is represented as G= L  l=1 βlq(φl ue,θl ue)pT(φl risD,θl risD)∈CQ×N,(2) where βlis the lth complex path gain, p(φl risD,θl risD)∈CN×1is the lth two-dimensional channel steering vector at the RIS, with φl risDand θl risDbeing the corresponding AoD and EoD, respectively. Similarly, q(φl ue,θl ue)∈CQ×1is the lth two-dimensional channel steering vector at the UE, with φl ue and θl ue being the associated directional angular parameters. The BS array response in y-z plane can be written as [a(φr bs,θr bs)]m=e−jπ[(my−1)sinθr bs sin φr bs+(mz−1)cosθr bs],(3) where m=mz+(my−1)Mz,my∈{1,...,M y},andmz∈ {1,...,M z}, as explained in [26]. Defining the spatial frequencies as μr bs =πsin θr bs sin φr bs and ψr bs =πcos θr bs, the channel steering vector can be expressed as the Kronecker product between two channel steering vectors, as a(μr bs,ψr bs)=ay(μr bs)⊗az(ψr bs)∈CM×1,(4) where ay(μr bs)=1,e −jμr bs ,...,e −j(My−1)μr bs T∈CMy×1 az(ψr bs)=1,e −jψr bs ,...,e −j(Mz−1)ψr bs T∈CMz×1. In a similar fashion, the channel steering vector b(φr risA,θr risA)can also be factorized as the Kronecker product of by(μr risA)and bz(ψr risA), respectively. The involved channel components Hand Gmentioned implicitly in (1) and (2), respectively, can be expanded in terms of their respective spatial frequencies and Kronecker products as H= R  r αrby(μr risA)⊗bz(ψr risA)[ay(μr bs)⊗az(ψr bs)]T,(5) G= L  l βlqy(μl ue)⊗qz(ψl ue)py(μl risD)⊗pz(ψl risD)T.(6) Using the properties (A⊗B)T=AT⊗BTand (A⊗B)(C⊗ D)=(AC ⊗BD), the rank-one channel matrices Hand Gdefined in (5)–(6) can be reformulated as a Kronecker product of their corresponding y and z domain components as H= R  r=1 αrby(μr risA)aT y(μr bs)   Hr y ⊗bz(ψr risA)aT z(ψr bs)   Hr z ,(7) G= L  l=1 βlqy(μl ue)pT y(μl risD)   Gl y ⊗qz(ψl ue)pT z(ψl risD)   Gl z .(8) Note that decomposing the channels into their yth and zth components allows us to exploit their intrinsic rank-one components, which is how we formulate our rank-one approximation problem. Under the assumption of THz propagation [27], the channels consist of a dominant line-of-sight (LOS) component and weaker non-LOS (NLOS) components. When the NLOS terms are negligible, the involved MIMO channels are approximately rank-one matrices, given as H≈Hy⊗Hz, where Hy∈CNy×My,andHz∈CNz×Mzare horizontal and vertical domain channels. Similarly, G≈Gy⊗Gz,whereGy∈CQy×Nyand Gz∈CQz×Nzare the associated horizontal and vertical component matrices. The received pilot sequence at the UE via RIS by assigning the k-th phase-shift pattern at the RIS is represented as Xk=Gdiag {ωk}HS +Vk∈CQ×T,(9) where k=1,...,K, is the associated received pilot block, while S∈CM×Tdenotes the pilot symbol matrix holding the length-Tpilot sequences. The diagonal matrix diag{ωk}∈CN×Nholds the RIS reflection coefficients (k-th column of the discrete Fourier transform (DFT) matrix) used at the k-th block for probing the propagation channel. It is further assumed that Hand Gremain constant during the channel probing stage. III. HIGHER-DIMENSIONAL RANK-ONE APPROXIMATIONS (HDR) METHOD Our idea is to exploit the Kronecker decomposition of BS-RIS and RIS-UE channels given in (7) and (8) to solve the parametric channel estimation as a whole, using a higher-order rank-one tensor approximation problem. For the convenience of presentation, in the following derivation steps, we focus on the signal part of (9) by excluding the noise term (which will be added later). Let Xk=Gdiag{ωk}HS,sothatXk= Xk+Vk. Using the properties vec{ABC}=(CT⊗A)vec{B}, vec{Adiag{d}B}=(BTA)d, and applying the vec{.}operator to Xkleads to xk . =vec{Xk}=ST⊗IQHTGωk∈CQT ×1. Defining the matrix X=[ x1,...,xK]∈CQT ×Kthat collects the resulting signal over Kreceived blocks yields X=ST⊗IQHTGΩ,(10) where Ω=[ω1,...,ωK]∈CN×K. Recasting (10) as tensor X=TG,H×1IQ×2ST×3ΩT,(11) where X∈CQ×T×Kis the pilot signal tensor, and TG,H∈CQ×M×N is the composite channel tensor expressed as TG,H=I3,N ×1G×2HT×3IN.(12) Using the structure in TG,H, the 1-mode unfolding of the noiseless pilot signal tensor Xdefined in (11) can be factorized [X](1)=[TG,H](1)(Ω⊗S)=G(INHT)(Ω⊗S).(13) A PARAFAC tensor model is proposed in [11] for unstructured channels. Let Ψ. =Ω⊗S∈CMN×TK be the effective joint training matrix that combines the pilot and RIS phase shifts matrices. We design the pilot sequence length Tand the RIS training block number Ksuch that Ψis row-wise orthonormal satisfying ΨΨH=IMN. Several choices ensuring this property are possible, such as a (possibly truncated) TK ×MDFT or Hadamard matrix. Note that such a joint design of the pilot and phase shift matrices requires TK ≥MN.1Now, adding 1The joint design of Ψ=Ω⊗Sgives us the flexibility to trade off the pilot sequence length and the number of RIS training blocks to meet the desired orthogonality. For instance, increasing Tand decreasing Kaccordingly yields a finer time delay resolution and a “wider” spatial probing. Alternatively, we could also split the joint filtering into two sequential filtering steps, which corresponds to E=˜ X× 2S∗×3Ω∗. Although this approach is a bit less complex, it would require satisfying both K≥Nand T≥M, which is a more restrictive condition than the joint design. Authorized licensed use limited to: TEL AVIV UNIVERSITY. Downloaded on August 10,2025 at 10:58:26 UTC from IEEE Xplore. Restrictions apply. 5178 IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 74, NO. 3, MARCH 2025 Algorithm 1: HDR Method. the noise tensor V∈CQ×T×K,wehave [˜ X](1)=G(INHT)Ψ+[V](1),(14) Applying right-filtering on (14) using the known joint pilot-RIS training matrix leads to [E](1). =[˜ X](1)ΨH≈G(INHT),(15) where E≈I3,N ×1G×2HT×3IN∈CQ×M×N.Let E. =[E]T (3)≈HTG∈CQM×N(16) be an LS approximation to the combined BS-RIS-UE “Khatri-Rao channel”. It was shown in [11] that unstructured estimates of the individual channel matrices Hand Gcan be obtained from the filtered signal in (15) by solving the Khatri-Rao factorization (KRF) problem ˆ H,ˆ G=argmin E−HTG 2 F.(17) Herein, we extend this idea by exploiting the channel’s geometric structure employing the Kronecker decomposition of BS-RIS and RIS-UE channels into their respective y and z components given in (7) and (8), which leads to the following mixed Khatri-Rao-Kronecker factorization problem arg min   E−(Hy⊗Hz)T(Gy⊗Gz)   2 F.(18) Interestingly, we arrive at an equivalent problem by resorting to the property (A⊗B)(C⊗D)=P1[(AC)⊗(BD)], where P1∈RQzMzQyMy×QzQyMzMyis a block permutation matrix the structure of which is derived in Appendix. Applying this property and defining Z. =P1Eallows us to rewrite (18) as a Kronecker factorization problem in terms of y and z Khatri-Rao channels arg min  Z−(HT yGy)⊗(HT zGz) 2 F.(19) Let z. =vec{Z}∈CQzMzQyMyNzNy×1and consider the property vec{A}⊗vec{B}=P2(vec{A⊗B}),where P2∈RQzMzNzQyMyNy×QzMzQyMyNzNyis a block permutation matrix, the structure of which is provided in the Appendix. Defining z. =P2z∈CQzMzNzQyMyNy×1allows to recast problem (19) in a vector form as arg min  z−vec HT yGy⊗vec HT zGz 2 F.(20) Note that vec{HT yGy}=ny(μy)⊗ay(μbs)⊗qy(μue)and vec {HT zGz}=nz(ψz)⊗az(ψbs)⊗qz(ψue), and define the sixthorder tensor Z∈CQz×Mz×Nz×Qy×My×Nyobtained by reshaping the vector z∈CQzMzNzQyMyNy×1accordingly. Such vector-tensor data mapping is represented by the following “tensorization” operator Z= tens(z)that rearranges the entries of the one-dimensional vector as Fig. 2. HDR receiver. those of a sixth-order tensor. Hence, problem (20) is equivalent to a rank-one sixth-order tensor approximation problem ˆ ay,ˆ az,ˆ qy,ˆ qz,ˆ ny,ˆ nz=argmin  Z−qz(ψue)◦az(ψbs)◦nz(ψz)◦qy(μue)◦ay(μbs)◦ny(μy) 2 F, (21) where ny(μy)=by(μrisA)py(μrisD)∈CNy×1and nz(ψz)=bz (ψrisA)pz(ψrisD)∈CNz×1are the effective y-th and z-th domains RIS steering vectors. The solution to the problem (21) can be obtained by resorting to the HOSVD algorithm, which in our case consists of six independent rank-one matrix approximations to each matrix unfolding of the tensor Z( due to limited space, we refer the reader to [24] for further details on the HOSVD steps). This problem can also be solved in an iterative way using the higher-order power method [28], which has also been exploited in [29] to decode multi-linear Kronecker-structured constellations. It is worth noting that the proposed method enjoys parallel processing since all the six rank-one approximation steps of the HOSVD algorithm can be executed in parallel, thus delivering fast estimates of the involved steering vectors, which is important in scenarios with low latency requirements. Finally, the combined complex path gain γ=αβ, i.e., ˆγ=ˆ E× 1ˆ qH×2ˆ aH×3ˆ nH,(22) where ˆ n=1 N(ˆ ny⊗ˆ nz),ˆ a=1 M(ˆ ay⊗ˆ az),andˆ q=1 Q(ˆ qy⊗ˆ qz).By taking into account the filtering step in (15), followed by the tensor reshaping E=reshape QM×N{[˜ X](1)ΨH}leading to (16), as well as the two block permutation steps Z. =P1Eand z. =P2vec{Z}used in the derivations from (18)–(19) and from (19)–(20), respectively, we can arrive at a simpler sequence of steps representing the proposed receiver processing linking the received pilot signals to the HOSVD input block, which is translated by the input-output relationship z= P2vec{P1E}=P2(IN⊗P1)vec{E}, or equivalently, by filtering, reshaping, and block permutation steps modeled as z=Pvec reshapeQM×N{[˜ X](1)ΨH}, where P. =P2(IN⊗P1)represents an effective QzMzNzQyMy Ny×QzQyMzMyNzNyblock permutation matrix that properly shuffles the filtered signal, which is then tensorized to feed the HOSVD block. Fig. 2illustrates the block diagram of the receiver processing steps, and a summary of the proposed HDR method is provided in Algorithm 1. Upon completion of the HOSVD stage, an estimate of the KhatriRao channel ˆ E=ˆ HTˆ Gis obtained from the estimated steering vectors by rebuilding the y and z channels. The complexity of HDR is dominated by the filtering and HOSVD steps and is given by OQ2MNTK +QMN (Qz+Qy+Mz+My+Nz+Ny). Authorized licensed use limited to: TEL AVIV UNIVERSITY. Downloaded on August 10,2025 at 10:58:26 UTC from IEEE Xplore. Restrictions apply. IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 74, NO. 3, MARCH 2025 5179 Fig. 3. NMSE-based performance comparisons. (a) NMSE performance for rank-one channels. (b) NMSE performance for higher-rank channels. Fig. 4. SE performance and complexity comparisons. (a) SE performance under channel estimation errors assuming multi-rank channels. (b) Computational Complexity comparison. IV. SIMULATION RESULTS We adopt a URA at the BS and the UE, and the channel parameters are modeled and chosen according to [27],[30].TheAoDφr bs,φl risDand AoA φr risA,φl ue are generated with uniform distribution by assuming one sector of a cell as φr bs,φ r risA,φ l risD,φ l ue ∼U(−60◦,+60◦).EoDθr bs,θl risD and EoA θr risA,θl ue are generated by assuming uniform distribution as θr bs,θr risA,θl risD,θl ue ∼U(90◦,130◦)[27],[30]. The complex path gains are generated as αr,β l∼CN(0,1).TheLOSpathisassumedtobe stronger than the NLOS paths. Note that, for higher-frequency signals, such as mmWave and THz bands, the channel measurement campaigns reveal that the signal power of the LOS component is about 13 dB higher than the sum of the power of NLOS components [31]. Moreover, the number of antennas deployed at the BS is M=16, where the number of antennas along the horizontal axis is My=4 and the number of antennas along the vertical axis is Mz=4. Similarly, the number of antennas deployed at the UE is Q=16 having Qy=4andQz=4. Finally, the number of reflecting elements deployed at the RIS is N= 64, where Ny=8, and Nz=8 are fixed for Figs. 3(a),to4(a). The total transmit power is assumed to be PT=1Wdefining SNR = PT/σ2 n. The NMSE for the reconstructed channel is given as NMSE = E  E−ˆ E   2 F/E2 F.(23) Fig. 3(a) evaluates the NMSE-based performance for rank-one channels of our proposed HDR method in comparison with the TenRICE [14], KRF and LS methods [11]. The Khatri-Rao channel given in (16) is reconstructed using the respective estimations of the channel steering vectors and combined complex path gain. The HDR method outperforms the TenRICE [14], the KRF, and the LS methods [11]. This is because of exploiting the inherent factorization structure and modeling this as a higher-order tensor to obtain a significant tensor gain. The TenRICE method [14] partially exploits the geometrical structure while both the LS and KRF methods [11] are unable to exploit the geometrical structure of the problem. The LS method for unstructured channel estimation satisfies the normalized CRLB [11], which is derived for unstructured channels. The HDR method and the TenRICE [14] method outperform the normalized CRLB due to noise rejection by fully and partially exploiting the geometrical structure of the channels, respectively. On the other hand, the KRF method [11] shows improvement compared with the LS method and the normalized CRLB due to noise rejection involved in the LS Khatri-Rao factorization step and due to the harmless influence of noise on the reconstruction of the factor matrices, respectively. The HDR method shows almost 100 times better performance compared to the KRF method [11], and 10 times compared to the TenRICE method [14] at 0 dB SNR. To further analyze the performance of our proposed HDR method in higher-rank channels, Fig. 3(b) shows that the HDR method still presents the best performance in the low SNR regime (−5 dB to 0 dB) but does not show further improvement after 5 dB SNR due to the approximation error caused by the assumption of a THz propagation scenario. The performance of the proposed HDR algorithm is further evaluated in Fig. 4(a), in terms of SE as defined in [32] assuming multi-rank channels. The active and passive beamforming vectors are calculated using the solution of [33], where an additional noise rejection is achieved. Here, the design of the active and passive beamforming vectors is based on the estimated channels. We also show the SE results for perfectly known channels as a benchmark. The HDR-based channel estimation shows almost similar performance compared to the TenRICE [14] and the KRF methods [11] for higher-rank channels. The reason behind the minor improvement in the case of the TenRICE method is that the SVD-based precoder and the combiner design are based on estimated multi-rank channels while the HDR method only estimates the dominant path of the channel, and, therefore, the SVD-based precoder and the combiner are based on the approximated rank-one channel, which shows a little degradation in the SE performance. Still, HDR shows a slightly better SE performance than KRF at SNR=−15 dB but a similar performance compared to TenRICE due to the approximation error of higher-rank channels as rank-one channels. Fig. 4(b) depicts a complexity comparison of our proposed HDR algorithm with the TenRICE [14], KRF, and LS methods [11]. Exploiting the inherent factorization structure and modeling it as a higher-order tensor help to reduce the overall complexity. The HDR method has almost the same computational complexity as the LS method [11] for different numbers of RIS elements but outperforms the TenRICE [14], and the KRF [11] methods. V. CONCLUSION We proposed a tensor-based channel estimation method for THz channels, which exploits the inherent geometrical structures of the involved channel matrices to recast the channel parameter estimation as a single sixth-order tensor rank-one approximation problem that can be efficiently solved using the HOSVD. The proposed HDR method outperforms the competing TenRICE, KRF, and LS methods in terms of estimation accuracy with a similar spectral efficiency. Additionally, HDR has a lower complexity than its competitors. APPENDIX BLOCK PERMUTATION MATRICES P1AND P2 The two block-permutation matrices used in step 4 of Algorithm 1 can be defined in a unified and general form using Kronecker products Authorized licensed use limited to: TEL AVIV UNIVERSITY. Downloaded on August 10,2025 at 10:58:26 UTC from IEEE Xplore. Restrictions apply. 5180 IEEE TRANSACTIONS ON VEHICULAR TECHNOLOGY, VOL. 74, NO. 3, MARCH 2025 of canonical unit vectors as follows Pr . = I  i=1 J  j=1 K  k=1 L  l=1 p(i,k,j,l) b(r)p(i,j,k,l)T a(r)∈RIKJL×IJKL,(24) r=1,2, where pa(r)=e(r) l⊗e(r) k⊗e(r) j⊗e(r) i(25) pb(r)=e(r) l⊗e(r) j⊗e(r) k⊗e(r) i.(26) For r=1, we have the following correspondences (I,J,K,L)↔ (Qz,Q y,M z,M y), while for r=2, we have (I,J,K,L)↔ (QzMz,Q yMy,N z,N y). 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, doi: 10.1109/JSAC.2020.3007211. [2] A. L. Swindlehurst, G. Zhou, R. Liu, C. Pan, and M. Li, “Channel estimation with reconfigurable intelligent surfaces—a general framework,” Proc. IEEE, vol. 110, no. 9, pp. 1312–1338, Sep. 2022, doi: 10.1109/JPROC.2022.3170358. [3] K. F. Masood, J. Tong, J. Xi, J. Yuan, and Y. Yu, “Inductive matrix completion and root-MUSIC-based channel estimation for intelligent reflecting surface (IRS)-aided hybrid MIMO systems,” IEEE Trans. Wirel. Commun., vol. 22, no. 11, pp. 7917–7931, Nov. 2023, doi: 10.1109/TWC.2023.3257138. [4] P. Wang, J. Fang, H. Duan, and H. Li, “Compressed channel estimation for intelligent reflecting surface-assisted millimeter wave systems,” IEEE Sig. Proc. Lett, vol. 27, pp. 905–909, 2020, doi: 10.1109/LSP.2020.2998357. [5] M. Li, S. Zhang, Y. Ge, F. Gao, and P. Fan, “Joint channel estimation and data detection for hybrid RIS aided millimeter wave OTFS systems,” IEEE Trans. Commun., vol. 70, no. 10, pp. 6832–6848, Oct. 2022, doi: 10.1109/TCOMM.2022.3199019. [6] X. Hu, C. Liu, M. Peng, and C. Zhong, “IRS-based integrated location sensing and communication for mmwave SIMO systems,” IEEE Trans. Wireless Commun., vol. 22, no. 6, pp. 4132–4145, Jun. 2023, doi: 10.1109/TWC.2022.3223428. [7] Z. Yu, X. Hu, C. Liu, M. Peng, and C. Zhong, “Location sensing and beamforming design for IRS-enabled multi-user ISAC systems,” IEEE Trans. Signal Process., vol. 70, pp. 5178–5193, 2022, doi: 10.1109/ TSP.2022.3217353. [8] P. Comon, X. Luciani, and A. L. F. de Almeida, “Tensor decompositions, alternating least squares and other tales,” J. Chemometrics, vol. 23, no. 7-8, pp. 393–405, 2009. [9] D. C. Araújo, A. L. F. de Almeida, J. P. C. L. Da Costa, and R. T. de Sousa, “Tensor-based channel estimation for massive MIMO-OFDM systems,” IEEE Access, vol. 7, pp. 42133–42147, 2019, doi: 10.1109/ACCESS.2019.2908207. [10] G. M. Zilli and W.-P. Zhu, “Constrained tensor decomposition-based hybrid beamforming for mmwave massive MIMO-OFDM communication systems,” IEEE Trans. Veh. Technol., vol. 70, no. 6, pp. 5775–5788, Jun. 2021, doi: 10.1109/TVT.2021.3076691. [11] G. T. de Araújo, A. L. F. de Almeida, and R. Boyer, “Channel estimation for intelligent reflecting surface assisted MIMO systems: A tensor modeling approach,” IEEE J. Sel Top Signal Process, vol. 15, no. 3, pp. 789–802, Apr. 2021, doi: 10.1109/JSTSP.2021.3061274. [12] G. T. de Araújo, A. L. F. de Almeida, R. Boyer, and G. Fodor, “Semi-blind joint channel and symbol estimation for IRS-assisted MIMO systems,” IEEE Trans. Signal Process., vol. 71, pp. 1184–1199, 2023, doi: 10.1109/TSP.2023.3263257. [13] K. Ardah, S. Gherekhloo, A. L. F. de Almeida, and M. Haardt, “TRICE: A channel estimation framework for RIS-aided millimeter-wave MIMO systems,” IEEE Signal Process. Lett., vol. 28, pp. 513–517, 2021, doi: 10.1109/LSP.2021.3059363. [14] S. Gherekhloo, K. Ardah, A. L. F. de Almeida, and M. Haardt, “Tensorbased channel estimation and reflection design for RIS-aided millimeterwave MIMO communication systems,” in Proc. IEEE 55th Asilomar Conf. Signals, Syst., Comput., 2021, pp. 1683–1689, doi: 10.1109/IEEE CONF53345.2021.9723362. [15] P. R. B. Gomes, G. T. de Araújo, B. Sokal, A. L. F. d. Almeida, B. Makki, and G. Fodor, “Channel estimation in RIS-assisted MIMO systems operating under imperfections,” IEEE Trans. Veh. Technol., vol. 72, no. 11, pp. 14200–14213, Nov. 2023, doi: 10.1109/TVT.2023.3279805. [16] K. B. A. Benício, A. L. F. de Almeida, B. Sokal, Fazal-EAsim, B. Makki, and G. Fodor, “Tensor-based channel estimation and data-aided tracking in IRS-assisted MIMO systems,” IEEE Wireless Commun. Lett., vol. 13, no. 2, pp. 333–337, Feb. 2024, doi: 10.1109/LWC.2023.3328838. [17] X. Zhang, X. Shao, Y. Guo, Y. Lu, and L. Cheng, “Sparsity-structured tensor-aided channel estimation for RIS-assisted MIMO communications,” IEEE Commun. Lett., vol. 26, no. 10, pp. 2460–2464, Oct. 2022, doi: 10.1109/LCOMM.2022.3194687. [18] A. L. F. de Almeida, G. Favier, J. P. C. L. da Costa, and J. C. M. Mota, “Overview of tensor decompositions with applications to communications,” in Signals and Images: Advances and Results in Speech, Estimation, Compression, Recognition, Filtering, and Processing, R. Coelho, V. Nascimento, R. de Queiroz, J. Romano, and C. Cavalcante, Eds. Boca Raton, FL, USA: CRC-Press, 2016, no. Chapter 12, pp. 325–356. [19] H. Chen, F. Ahmad, S. Vorobyov, and F. Porikli, “Tensor decompositions in wireless communications and MIMO radar,” IEEE J. Sel. Topics Signal Process., vol. 15, no. 3, pp. 438–453, Apr. 2021, doi: 10.1109/JSTSP.2021.3061937. [20] A. L. F. de Almeida, G. Favier, and J. C. M. Mota, “PARAFAC-based unified tensor modeling for wireless communication systems with application to blind multiuser equalization,” Signal Process., vol. 87, no. 2, pp. 337–351, 2007, doi: 10.1016/j.sigpro.2005.12.014. [21] G. Favier and A. L. F. de Almeida, “Overview of constrained PARAFAC models,” EURASIP J. Adv. Signal Process., vol. 2014, pp. 1–25, 2014, doi: 10.1186/1687-6180-2014-142. [22] C. Jia, J. Cheng, H. Gao, and W. Xu, “High-resolution channel estimation for intelligent reflecting surface-assisted mmwave communications,” in in Proc. IEEE 31st Annu. Int. Symp. Personal, Indoor Mobile Radio Commun., London, U.K., 2020, pp. 1–6, doi: 10.1109/ PIMRC48278.2020.9217312. [23] I.-S. Kim, M. Bennis, J. Oh, J. Chung, and J. Choi, “Bayesian channel estimation for intelligent reflecting surface-aided mmwave massive MIMO systems with semi-passive elements,” 2022. [Online]. Available: https: //arxiv.org/abs/2206.06605 [24] L. De Lathauwer, B. De Moor, and J. Vandewalle, “A multilinear singular value decomposition,” SIAM J. Matrix Anal. Appl., vol. 21, no. 4, pp. 1253–1278, 2000. [25] R. W. Heath, N. González-Prelcic, S. Rangan, W. Roh, and A. M. Sayeed, “An overview of signal processing techniques for millimeter wave MIMO systems,” IEEE J. Sel. Topics Signal Process., vol. 10, no. 3, pp. 436–453, Apr. 2016, doi: 10.1109/JSTSP.2016.2523924. [26] Fazal-EAsim, F. Antreich, C. C. Cavalcante, A. L. F. de Almeida, and J. A. Nossek, “Two-dimensional channel parameter estimation for millimeterwave systems using butler matrices,” IEEE Trans. Wirel. Commun., vol. 20, no. 4, pp. 2670–2684, Apr. 2021, doi: 10.1109/TWC.2020.3043958. [27] C. Han et al., “Terahertz wireless channels: A holistic survey on measurement, modeling, and analysis,” IEEE Commun. Surv. Tuts., vol. 24, no. 3, pp. 1670–1707, thirdquarter 2022, doi: 10.1109/COMST.2022.3182539. [28] P. Regalia and E. Kofidis, “The higher-order power method revisited: Convergence proofs and effective initialization,” in Proc. 2000 IEEE Int. Conf. Acoust., Speech, Signal Process. Proc., 2000, pp. 2709–2712, doi: 10.1109/ICASSP.2000.861047. [29] F. E-Asim, A. L. F. de Almeida, M. Haardt, C. C. Cavalcante, and J. A. Nossek, “Rank-one detector for Kronecker-structured constant modulus constellations,” IEEE Signal. Process. Lett, vol. 27, pp. 1420–1424, Jul. 2020, doi: 10.1109/LSP.2020.3010133. [30] M. K. Samimi and T. S. Rappaport, “3-D millimeter-wave statistical channel model for 5G wireless system design,” IEEE Trans. Microw. Theory Techn., vol. 64, no. 7, pp. 2207–2225, Jul. 2016, doi: 10.1109/ TMTT.2016.2574851. [31] Z. Muhi-Eldeen, L. Ivrissimtzis, and M. Al-Nuaimi, “Modelling and measurements of millimeter wavelength propagation in urban environments,” IET Microw., Antennas & Propag., vol. 4, no. 9, pp. 1300–1309, 2010. [32] F.-E. Asim, C. C. Cavalcante, F. Antreich, A. L. F. de Almeida, and J. A. Nossek, “Efficient hybrid A/D beamforming for millimeter-wave systems using butler matrices,” IEEE Trans. Wirel. Commun., vol. 22, no. 2, pp. 1001–1013, Feb. 2023, doi: 10.1109/TWC.2022.3200298. [33] A. Zappone, M. Di Renzo, F. Shams, X. Qian, and M. Debbah, “Overheadaware design of reconfigurable intelligent surfaces in smart radio environments,” IEEE Trans. Wirel. Commun., vol. 20, no. 1, pp. 126–141, Jan. 2021, doi: 10.1109/TWC.2020.3023578. Authorized licensed use limited to: TEL AVIV UNIVERSITY. Downloaded on August 10,2025 at 10:58:26 UTC from IEEE Xplore. Restrictions apply.