Joint Channel Estimation and Data Rate Maximization for Intelligent Reflecting Surface Assisted Terahertz MIMO Communication Systems

Terahertz (THz) communications recently attract significant attention and become an emerging technology pillar for sixth generation (6G) wireless systems. Due to the serious path attenuation of THz signals, THz communication is applicable for the short-distance indoor scenarios. However, the THz waves are easily blocked by obstacles, leading to a communication interruption. To this end, an intelligent reflecting surface (IRS), which interacts with incident THz waves in a controlled manner by adjusting the discrete phase shifts of the IRS elements, is considered as a promising technology to mitigate blockage vulnerability and enhance coverage capability for indoor scenarios. In light of graphene-enabled hardware structure of an IRS, the IRS-assisted THz multiple-input multiple-output (MIMO) system model is developed. Moreover, an iterative atom pruning based subspace pursuit (IAP-SP) scheme is developed for channel estimation. Compared to the classical subspace pursuit (SP) scheme, the proposed IAP-SP algorithm can substantially reduce the computational complexity while maintaining accurate channel recovery. With the estimated channel, a data rate maximization problem is formulated, which can be converted to a discrete phase shift search problem. The exhaustive search method is firstly proposed to obtain the optimal transmission rate but endure extremely high computational burden. Then, a local search method is proposed to decrease the number of possible discrete phase candidates of IRS while undergoes obvious performance loss. Interestingly, a novel feedforward fully connected structure based deep neural network (DNN) scheme is put forward, which has the ability to learn how to output the optimal phase shift configurations by inputting the features of estimated channel. Simulation results demonstrate that, in contrast with the exhaustive search scheme and the local search scheme, the proposed DNN-based scheme achieves a near-optimal communication rate performance. Meanwhile, the DNN-based scheme enormously alleviates the computational complexity and allows for dynamic parameter adaption in rapid-varying channel conditions.

[3]- [6], terahertz (THz) band (from 0.1 THz to 10 THz) has drawn more and more attention for 6G wireless networks. Proverbially, THz frequency band bridges the gap between millimeter wave (mmWave) band and optical frequency band. Compared with mmWave communication, THz communication possesses more abundant spectrum resources and is able to provide much higher transmission rates from hundreds of Gigabit-per-second (Gbps) to several Terabit-per-second (Tbps). In terms of the optical communication, THz communication is more favorable to adapt the inconvenient climate conditions including fog, dust, snow and so on. Besides, THz frequency band has the possibility to employ the reflected paths to enhance the communication performance in the complex propagation environments. Stated thus, the THz frequency band is envisioned as the most appropriate candidate to enable the 6G communication systems.
Except for the advantages mentioned above, THz communication still faces some challenges when it is applied to the indoor application scenario. Due to the serious path attenuation and strong directivity of THz waves, the line-ofsight (LOS) THz path can be easily blocked by the indoor obstacles (e.g., walls, doors, wardrobes, human bodies) [7], [8]. Thus, the non-line-of-sight (NLOS) paths play a main role in such a communication scenario. To tackle this issue, an intelligent reflecting surface (IRS) is newly proposed to improve the spectral efficiency and enhance the coverage performance recently [9]- [16]. The IRS is a physical meta-surface consisting of a large number of small-unit reflectors. In general, the IRS is equipped with a number of passive elements and is controlled by a central processor. Each element of the IRS is able to reflect incident THz waves independently with an adjustable phase-shift. Different from the existing candidate techniques, such as conventional reflect-arrays [17], [18], amplify-and-forward relaying [19]- [21] and backscatter communication [22]- [24], IRS owns some distinct advantages. In contrast to reflect-arrays and backscatter communication that require the radio frequency (RF) chains to realize the phase shifts, the hardware complexity of IRS is much lower due to the lack of RF chains, and thus much more reflecting elements can be scaled up in practice. Additionally, since IRS is designed with passive reflecting elements, IRS is also an energy efficient technique. Instead, the amplify-and-forward relaying system usually works with high power consumption. Based on the inherent characteristics, IRS has become a popular research priority extensively.
To the best of our knowledge, the initial concept about IRS is an intelligent wall, which is regarded as an autonomous part of a smart indoor environment to control radio coverage by make use of active frequency-selective surfaces [25]- [27]. In [28], the reconfigurable reflect-array which has the ability to enhance the strength of received signal and eliminate interference simultaneously, is proposed to increase the spectrum sharing capacity for indoor environments. Recently, the IRS-enhanced multiple-input single-output (MISO) wireless systems are explore to maximize the total received power at the user side by jointly optimizing the active beamforming at the base station (BS) and passive beamforming at the IRS [29]- [31]. However, these research contents consider that the phase shifts of the IRS elements are continuous, and the IRS with continuous phase shift results in high fabrication cost and complex hardware design. Later, in [13], [32], the IRS with only a limited number of phase shifts is adopted optimize the energy efficiency and sum-rate performance for both the discrete phase shifts of IRS elements and transmit power allocation. It is worth noting that all the prior work assumes that the channel status information (CSI) is available at the BS side. In other words, we need to solve the channel estimation problem firstly. The main challenge of the channel estimation in IRS-assisted communication system is that the IRS elements are passive and are unable to process the incident signals. Given this, the authors in [33] propose a novel IRS structure that combines the passive elements with a few active elements. These active elements are connected to the baseband processor, and then are used to obtain the whole phase shifts of IRS elements by training a small number of active elements. Except for the low channel estimation accuracy, these active elements increase the power consumption of the IRS-assisted communication system. Instead, the authors present a cascaded channel estimation method with only passive reflecting elements [34]. Nevertheless, this method does not consider the channel characteristics (e.g., sparsity, path loss) in THz band, and it is unable to estimate the BS-IRS channel and IRS-user channel simultaneously. Through the above discussion, the IRS-assisted THz communication system is still worthy of further exploration.
In this paper, we present an IRS-assisted THz multipleinput multiple-output (MIMO) communication system for 6G indoor application scenarios. Considering the hardware design of IRS in THz band, the tunable graphene, which is a two-dimensional carbon material, is selected to achieve a wide phase response range of IRS elements by controlling the applied voltage. With the given hardware structure of graphene-based IRS, the IRS-assisted THz communication system is modeled firstly. Next, before the data transmission, it is necessary to acquire the CSI according to the system model where the geometric channel model is adopted that contains few limited scattering paths [35]- [37]. Then, by exploiting the sparsity feature of the THz channel, the compressed sensing (CS) technique [38] can be utilized to estimate the channel parameters, including path gain, angle of arrival (AoA) and angle of departure (AoD) [39]. With the estimated CSI, we then optimize the communication rate maximization problem by searching the optimal phase shift combination of the IRS elements. Furthermore, three different phase shift search methods are investigated, such as exhaustive search method, local search method and the deep neural network (DNN) based method. Finally, the communication rate performance and the complexity analysis of diverse phase shift search methods are also provided in this paper. The major contributions can be summarized as follows: • To begin with, to obtain a wide phase response range of reflecting elements, the electric properties of graphene at THz frequency band is investigated by using the CST Microwave Studio which confirms that the tunable graphene is an effective material to design IRS. Based on the electric characteristics of the graphene, the hardware structure of graphene-based IRS is creatively proposed in this paper, where the phase shift of each reflecting element belongs to a discrete phase set. With the designed IRS architecture, the IRS-assisted THz communication model is derived in detail.
• To realize the CSI acquisition of the IRS assisted THz system, an iterative atom pruning based subspace pursuit (IAP-SP) scheme with low complexity is proposed. Compared with the classical subspace pursuit (SP) method [40], the proposed IAP-SP is able to reduce the dimension of the sensing matrix without any recovery performance loss, where the columns (atoms) of the sensing matrix that least correlate with the signal residual tends to be eliminated in the iterative process.
• With the estimated channel information, we note that the communication rate maximization problem of IRS-assisted THz system is a non-convex problem, and then we convert this NP-hard problem to the phase shift search problem. Although the exhaustive search method achieves the optimal performance, it inherits the really high complexity. Fortunately, the local search method is proposed to decrease the computational burden by searching the partial phase shift combinations.
• Moreover, in order to fully explore the relationship between the channel features and the phase shift selection of IRS elements, the DNN-based phase shift search method is employed in this paper to optimize the communication rate. Concretely, the proposed DNN-based search method has the ability to unveil the mapping between the channel features and the phase shift selection. Thus, the DNN-based search method can realize the near-optimal communication rate performance with extremely low complexity. The rest of this paper is organized as follows. In Section II, we introduce the hardware structure of tunable graphenebased IRS, and the simulation analysis of tunable graphene properties are also presented. In Section III, the IRSassisted THz communication system and channel model are presented. Also, the channel estimation problem is settled. In Section IV, the communication rate maximization problem is formulated and three different solutions are proposed to search the optimal phase shift combination. Then the simulation results are presented in Section V to verify the proposed algorithms. Finally, we conclude this paper in Section VI.
In this paper, we use the following notation: A is a matrix, a is a vector, a is a scalar. A F is the Frobenius norm, whereas A H , A * , A T , A −1 , A † and rank (A) are conjugate transpose, conjugate, transpose, inverse, pseudo-inverse and the rank of A, respectively. A ⊗ B is the Kronecker product of A and B. diag (a) is a diagonal matrix with elements of a on its diagonal. Tr (A) is the trace of matrix A. E [·] is used to denote the expectation. {·} is the operation to extract the real part of a vector, and {·} is the operation to extract the image part of a vector. supp (a) denotes the set of non-zero elements in vector a. vec (A) is the column-ordered vectorization of matrix A, and vec −1 (A) is the reverse operation of vec (A).
A [ ] represents the sub-matrix of A that only contains those columns contained in the index set . = supp (a) denotes the index set corresponding to non-zero elements of vector a. M (a, K ) is a index set that corresponds to the largest K elements of a. Z (a) is a index set that corresponds to the zero elements of a, and Z −1 (a) is the reverse operation of Z (a), respectively.

II. HARDWARE PRELIMINARIES OF IRS
In this section, the design theory and working principle of IRS will be explained, and an achievable structure of IRS is proposed, which is suitable for THz communication systems. From the perspective of hardware structure, IRS employed in THz communication scenarios is a kind of tunable metasurface, which refers to a series of artificial structures with special physical properties, such as negative refractive index, near-zero index and so on [41]. The IRS are composed of periodically arranged IRS elements (also namely meta-atoms), and each element can control the amplitude and phase of the reflected waves. It is worth noting that the size of each IRS element is at a subwavelength level, and thus the propagation state of the reflected waves can be controlled through the arrangement of IRS elements. In general, the structure of each IRS element is designed first to attain the amplitude and phase control of reflected waves, and then these reflecting elements will be arranged properly to construct the reflective wavefront.

A. ELECTRIC PROPERTIES OF GRAPHENE IN THz BAND
As previously discussed, we know that IRS is a kind of tunable metasurface. In order to actively control the reflected waves, the existing devices (e.g., varactor [42], transistor [43]) can be integrated with the structure design. However, due to the short wavelength of the THz waves, the size of each IRS element becomes much smaller and it is difficult to combine those semiconductor devices with IRS elements. In this case, adopting the novel tunable materials, such as graphene, is a feasible choice.
Graphene is a monolayer of carbon atoms arranged in a honeycomb lattice, which has attracted tremendous attention for its extraordinary properties. To verify that the conductivity of graphene can be effectively tuned by applied voltage [44], the complex conductivity of graphene needs to be considered firstly. According to the [45], the complex conductivity of graphene consists of interband and intraband transition contributions, where the interband transition contribution can be neglected and the intraband transition contribution plays a primary role. Then, the intraband transition contribution of VOLUME 8, 2020 the complex conductivity of graphene can be given as where ω is the angular frequency, e is the elementary charge, h is the reduced Planck constant, k B is the Boltzmann constant, T is the temperature, E F is the Fermi level and τ = 1 × 10 −12 s is the relaxation time, respectively. In terms of the THz frequency band, we know that the conductivity of graphene is mainly dependent on intraband transition, which is related to Fermi energy [46]. Furthermore, the Fermi energy can be formulated as where the ν F = 1 × 10 6 m · s −1 is the Fermi velocity, ε is the dielectric constant of the substrate, t s refers to its thickness and V g is the applied voltage, respectively. From the formulation above, it is clear that the conductivity of graphene can be controlled by applying different voltages.  Fig. 1 illustrates a typical hardware structure of graphenebased IRS consisting of a array of reflecting elements. These subwavelength elements are closely arranged and only a single element is unable to realize the function of phase control due to the scattering of THz waves. For example, the size of IRS is 20mm × 20mm while such an array contains 100×100 reflecting elements. From Fig. 1 (a) we can see that, each reflecting element of IRS comprises the graphene patch, the silicon dioxide substrate and the gold ground plane. This hardware structure of an IRS can be considered as a Fabry-Perot resonance cavity [47], and the reflecting phase can be written as

B. HARDWARE DESIGN OF GRAPHENE-BASED IRS
where w is the width of graphene patch, k 0 = 2π λ 0 is the wave number of free space, ϕ is the phase shift of reflection wave, m is an integer, n sp is the effective refraction index of the surface plasmon, which is the function of graphene's effective permittivity ε eff and the angular frequency ω. The permittivity of the graphene can be calculated as where t g is the thickness of graphene. According to (3) and (4), we can obtain the relationship between the phase response and graphene's conductivity, that is, the phase shift of each IRS element can be altered by changing graphene's conductivity through applied voltage. Fig. 1 (b) shows a schematic diagram of controlling each reflecting element. The attached material, termed ion-gel, is used to connect graphene and the electrode, and the voltage V g is applied between the electrode and the ground plane. To verify the effectiveness of the graphene-based hardware structure, simulation results of each IRS element by using CST Microwave Studio is shown in Fig. 2. Fig. 2 depicts the phase response of the IRS elements along with the frequency and the diverse Fermi level values. One may note that the proposed IRS just work well at a narrow frequency band (about 3 GHz) once the hardware structure of IRS is fabricated. For a normally incident wave at 0.22 THz, the controllable phase shift range of each IRS element is able to reach the phase range of 351.25 • with E F = 0.3 eV producing a phase of 88.71 • and 0.4 eV generating a phase of −262.54 • . More importantly, by properly adjusting the parameters of hardware structure, like the length of the substrate L and the width of the graphene patch w, the available phase shift range can be consistent with the interval [0, 2π] [48]. Remarkably, we can take any desired phase shift from the phase range by applying the voltages continuously, which results in high complexity hardware circuits and expensive system overhead in practice.
Based on the above analysis about the hardware characteristics of IRS, we indicate that the phase-shift ϕ of each IRS element can be adjusted within the range [ϕ min , ϕ max ], where ϕ min is the controllable minimum phase shift and ϕ max is the controllable maximum phase shift for a given IRS structure. Since the power consumption and hardware cost of IRS with discrete phase-shifts is much more practical than the IRS with continuous phase-shifts, we therefore consider discrete structure of the phase-shifts in our work. The phase shift of each reflecting element ϕ belongs to a finite phase set F = ϕ min , ϕ min + ϕ · · · , ϕ min + ϕ 2 b − 1 , where b is a bit quantization number and ϕ = (ϕ max − ϕ min ) 2 b is the phase interval. For example, when b = 2 bits, ϕ min = 0, and ϕ max = 2π, F = 0, π 2, π, (3π ) 2 .

III. SYSTEM MODEL AND CHANNEL ESTIMATION
In this section, the IRS-assisted THz communication model is introduced from a MIMO perspective. Then, the channel estimation problem and the corresponding solution are formulated in the following.

A. SYSTEM MODEL
Consider a downlink MIMO THz communication system and the communication link is assisted with IRS that consists of N IRS passive reflecting elements, as shown in Fig. 3. The hybrid beamforming architecture is adopted, since we take the hardware implementation and system cost into consideration. Assume that a BS with N BS antennas and M BS RF chains serves a single mobile station (MS) with N MS antennas and M MS RF chains. In general, the number of antennas is larger than the number of the RF chains due to the serious power consumption of the RF circuits, i.e. N BS > M BS , N MS > M MS . There are N s data streams that are transmitted between BS and MS. Due to the complex environment of indoor communication scenarios, the LOS path between BS and MS is usually blocked by the obstacles and NLOS paths play a major role. Also, the NLOS paths that are just reflected by an IRS for the first time are considered here. Furthermore, a time-division duplexing (TDD) scheme can be used for the downlink CSI acquisition by utilizing the characteristic of the channel reciprocity.
For the channel training scheme in such a downlink MIMO system, we suppose that the BS uses P different precoding vectors at P successive time frame. Each time frame contains Q time slots, and a combining vector is employed by the MS to detect the transmitted signal at each time slot. At pth time frame, the transmitted beamforming vector at the BS side can be expressed as where s p ∈ C N s is the transmitted pilot symbol vector and F p is the N BS × N s BS precoding matrix that combines the RF precoding vector F RF p ∈ C N BS ×M BS and the baseband precoding vector F BB p ∈ C M BS ×N s . Note that the pilot vector satisfies E[sps H p ] = P s N s I N s and P s is the total transmitted power of N s pilot symbols.
At the receiver side, the MS successively adopts Q combining vectors to detect the transmitted signal. For a given time frame p, the MS employs a combining vector w q to combine the received signal at time slot q, then the resulting signal formulation is written as where w q ∈ C N MS is the combining vector at qth time slot, z q,p is the additive Gaussian noise and H ∈ C N MS ×N BS is the cascaded channel matrix, respectively. More importantly, the cascaded channel H is composed of the BS-IRS channel which differs from conventional channel model without IRS and will be later introduced in this paper. After collecting the received signals processed by Q combining vectors, the received signal vector at pth time frame can be formulated as where With regard to the total P time frames, the MS utilizes the identical combiner W ∈ C N MS ×Q to detect the received signals, and thus, the received matrix at the MS side can be written as where Y = [y 1 , y 2 , · · · , y P ] is a Q × P received matrix, is a Q × P noise matrix by concatenating the P noise vectors, respectively.

B. CHANNEL MODEL
Before formulating the channel estimation problem, the channel model for the IRS-assisted MIMO THz communication system should be described firstly. Different from the conventional channel model [39], [49], [50], the channel model of the IRS-based THz system includes the BS-IRS channel H 1 , the IRS-MS channel H 2 , and the phase shift matrix diagonal matrix to indicate the phase shifts of IRS elements. Then, the cascaded channel of the IRS assisted system can be expressed as Without loss of generality, in THz frequency band, we assume that both H 1 and H 2 are able to be presented by the geometric channel model with few scattering paths. Furthermore, each scattering path is supposed to decide a single propagation path. Based on the geometric channel model, the BS-IRS channel H 1 can be formulated as where L 1 denotes the scattering paths of BS-IRS link, α l 1 is the complex gain of l 1 th path, γ AOA l 1 ∈ [0, 2π] is the AoA and φ l 1 ∈ [0, 2π] is the AoD of l 1 th path corresponding to BS and IRS, respectively. The complex gains of the propagation paths are mainly influenced by the path loss and molecular absorption, and more details are recommended in [51]. For simplification, the uniform linear arrays (ULAs) are considered here. Hence, we can get the array response vectors a BS (φ l 1 ) ∈ C N BS and a IRS (γ AOA l 1 ) ∈ C N IRS at the BS and IRS. Specifically, the array response vectors a BS (φ l 1 ) and a IRS (γ AOA l 1 ) can be written as follow where λ is the wavelength of the THz signals and d is the distance between adjacent antenna elements or IRS elements, which is usually defined as d = λ/2. Apart from the normal expression (11), a more compact form of BS-IRS channel can be formulated as where α= N BS N IRS L 1 α 1 , α 2 , · · · , α L 1 T , and the array response matrices can be respectively given as Similar to the BS-IRS channel H 1 , the IRS-MS channel H 2 can be formulated as where By substituting (14) and (17) into (10), the entire channel model can be written as

C. PROPOSED CHANNEL ESTIMATION SCHEME
In this subsection, the channel estimation problem and the corresponding solution of IRS-aided THz MIMO system is to introduced in detail. The main challenge is to estimate the BS-IRS channel and the IRS-MS channel simultaneously.
With the geometric channel model mentioned above, we discover that the channel estimation problem can be converted to the sparse recovery problem by utilizing the poor scattering characteristics of THz channel. Then, the CS technique is an efficient tool to handle such a sparse problem.
To be specific, before formulating the spare problem, we need to vectorize the received matrix Y in (9), which can be expressed as For the sake of forming the CS problem, the grid quantization procedure is achieved necessarily [39]. Here, we consider that the AoAs and AoDs of H 1 and H 2 are selected from a uniform grid of N 1 (N 1 >> L 1 ) points and a uniform grid of N 2 (N 2 >> L 2 ) points, respectively. Once the grid quantization procedure of azimuth angles is completed, the quantized AoAs and AoDs of H 1 and H 2 can be described, i.e.,φ l 1 ,γ AOA . According to (15) 99570 VOLUME 8, 2020 and (16), the quantized array response matrices can be defined asĀ Subsequently, the cascaded channel model (19) after the grid quantization operation can be rewritten as where G 1 ∈ C N 1 ×N 1 is the quantized path gain matrix of H 1 and G 2 ∈ C N 2 ×N 2 is the quantized path gain matrix of H 2 , respectively. It is worth noting that the non-zero elements of G 1 are equivalent to the elements of α due to the grid quantization principle [52]. Also, such a mapping relationship is available for G 2 and β.
Combining (20) and (22), the sparse formulation of CS problem can be expressed as where is a PQ × 1 vectorized noise vector, respectively. Based on structure properties of (23), it can be treated as a sparse recovery problem, and thus different CS techniques are able to be leveraged to estimate the channel parameters, such as SP, matching pursuit (MP) [53], orthogonal matching pursuit (OMP) [54], [55], and iterative hard thresholding (IHT) [56].
Proof: The detailed proof for the above (23) is provided in Appendix A.
In terms of the sparse reconstruction problem (23), the CS technique involves the unacceptably computational complexity due to the high dimension of sensing matrix A v . For instance, we suppose P = 6, Q = 12, N 1 = 30, N 2 = 30, then the dimension of A v is 72 × 30 4 . It is worth noting that the high-dimensional sensing matrix is caused by the double grid quantization operation, where the double grid quantization means that the grid quantization of H 1 and the grid quantization of H 2 are operated synchronously. Thus, the critical objective is to reduce the dimension of A v as well as the complexity.
In this paper, the SP algorithm is employed to perform the better recovery performance, which is able to both remove and add the elements in an active set. To efficiently strike the computational difficulty, we propose a low complexity IAP-SP algorithm. In recovery process of classical SP, we find that plenty of columns of the sensing matrix A v have no correlation with the signal residual. In other words, partial elements of the resulting vector that combines the sensing matrix and the residual are zero. Consequently, during iterative process, it is redundant to make the columns of A v with no contributions to the signal residual participate in calculations. According to this phenomenon, the proposed IAP-SP algorithm attains an attractive complexity reduction by iteratively eliminating the redundant columns in sensing matrix. The proposed low complexity IAP-SP algorithm is described in Algorithm 1.

Algorithm 1 Proposed IAP-SP Algorithm
, then t = t + 1, return to step 2, end if 9: Compressed recovery vector: g 10: Output estimated vector: Now, some important steps of Algorithm 1 are to be explained in the following. Firstly, as for step 1-2, the algorithm initialization and the support set estimation are presented, respectively. Secondly, in step 3, the index set of the redundant columns in A v is updated, where these redundant columns are not involved in subsequent calculations. It is worth mentioning that step 3 is a very pivotal operation in our proposed low complexity IAP-SP algorithm. Then, for steps 4-9, the operation principle of the proposed low complexity IAP-SP scheme is the same as classical SP algorithm VOLUME 8, 2020 in order to find the K -sparse solution to (23). Finally, for step 10, we need to fill the zero-elements to the corresponding positions that is caused by step 3. According to the analysis of the algorithm steps, the computational complexity of classical SP algorithm and our proposed IAP-SP in terms of complex multiplications can be expressed as respectively, where T represents the assumed maximum iterations.
Although the sparse vectorḡ v is reconstructed, estimating theḠ 1 andḠ 2 separately is still a tricky problem, where the recoveredḡ v ,Ḡ 1 , andḠ 2 correspond to actual g v , G 1 , and G 2 , respectively. Considering the practical communication scenarios, the positions of BS and IRS are fixed once the communication environment is determined. Therefore, we assume that the sum of the total path gains between BS and IRS is a constant C 0 for a specific environment, where C 0 is decided by the transmit power at BS side. Then we can get the following expression as whereᾱ l 1 ∈ supp vec Ḡ 1 and l 1 ∈ [1, · · · , L 1 ]. Here, we defineḠ v = vec −1 ḡ v from a block-matrix perspective, and thenḠ v can be rewritten as where Combining (26) and (27), the matrixḠ 2 can be derived as N 2 ]. Subsequently, the desired matrixḠ 1 can be formulated as Proof: The detailed proofs for (28) and (29) are provided in Appendix B and C, respectively.
Once the path gain matricesḠ 1 andḠ 2 are completed, the AoAs and AoDs can be determined as well as the BS-IRS channelH 1 and IRS-MS channelH 2 by utilizing the regulation of grid quantization accordingly. For evaluating the channel estimation performance, the normalized mean squared error (NMSE) is defined to compare the validity of diverse sparse recovery methods, which can be formulated as whereH is the estimated channel andH =H 2 H 1 .  Taking the computational burden into consideration, the number of grid points of the grid quantization process is set as N 1 = N 2 = 30. From Fig. 4 we note that the classical SP achieves the best sparse reconstruction performance compared to the other considered algorithms, in which the selected atoms are possible to be removed or added to form the estimated support set during the iterative reconstruction process, just as mentioned in [40]. Interestingly, our proposed IAP-SP algorithm is able to realize approximately the same NMSE performance as classical SP, while it outstrips the rest CS techniques, i.e., MP [53], OMP [54] and IHT [56]. Especially, when SNR = 20, the NMSE values of the IAP-SP (basically same as classical SP), OMP, MP and IHT are about 3 × 10 −2 , 7 × 10 −2 , 2.9 × 10 −1 and 3.2 × 10 −1 , respectively. Additionally, by substituting these parameter settings of the THz communication system into (24) and (25), the proposed IAP-SP algorithm can achieve around 88.79% complexity reduction when the number of iterations is T = 9. Through the analysis of computational complexity and recovery performance, the proposed IAP-SP scheme can make a better trade-off between computational calculations and recovery accuracy compared with classical SP scheme. With the number of the grid quantization points increasing, more complexity reduction is to be attained correspondingly. In addition, according to [57], further improved performance can be achieved by utilizing the redundant dictionary. Therefore, the proposed IAP-SP algorithm is an efficient tool to solve the channel estimation problem of the IRS-assisted THz MIMO communication system.

IV. DATA-RATE MAXIMIZATION FOR IRS-ASSISTED THz COMMUNICATION SYSTEM
With the estimated CSI analyzed above, the major objective is to handle the data-rate maximization problem of the IRS-assisted THz MIMO communication system by jointly designing the IRS matrix and the hybrid combiner matrix W = W RF W BB , which is an intractable non-convex optimization problem. Given this, we reconsider this NP-hard problem as the phase shift search problem of IRS elements. In terms of the phase shift search issue, we propose three different methods are proposed to optimize the IRS phaseshift matrix, including the exhaustive search scheme, local search scheme and deep learning based scheme. After that, we leverage the classical singular value decomposition (SVD) operation to design the optimal combining matrix W as well as the maximum communication rate.
Firstly, we present the achievable communication rate of the IRS-aided THz MIMO system as where ρ is the total transmission power, δ 2 is the noise power, and F = F RF F BB is the hybrid precoding matrix, respectively. We note that the BS-IRS channelH 1 has been estimated in previous section. Considering the practical communication scenarios, the locations of the BS and the IRS are fixed andH 1 also maintains unchanged. Given this, the precoder F at the BS side is a constant matrix and can be obtained by the SVD operation ofH 1 . Specifically, sincē H 1 is able to be decomposed intoH 1 = U 1 1 V H 1 , then we have the relation as F =Ṽ 1 , whereṼ 1 is a sub-matrix with the first N s column vectors of V 1 [49]. Once the optimal hybrid precoder F or combiner W is obtained, the OMPbased solution is adopted to decompose F or W into RF part and baseband (BB) part separately [58]. For simplicity, we assume that the IRS satisfies a perfect reflecting mode where each IRS element is set as {µ i } N IRS i=1 = 1. Then, the datarate maximization problem can be written as Based on the above discussion, the optimization objective (32) is a a non-convex and discrete optimization problem with two matrix variables, and W . For the sake of solving the optimization problem (32), existing optimization techniques are hard to operate directly. Nevertheless, it is worth noting that there are still some useful properties existing in problem (32). Firstly, the number of possible candidates of the phase shift matrix is finite, because each entry {ϕ n } N IRS n=1 is discrete according to the hardware structure of IRS. Secondly, the combining matrix W is an unconstrained combiner matrix. Moreover, the phase shift matrix and the combining matrix W are mutually independent in the IRS assisted THz communication system.
Motivating by these distinguishing features, one available way of settling the above optimization problem is to first optimize the phase shift matrix and later design the combining matrix W . To achieve this, we suppose that the effective channelH e is defined asH e =H 2 H 1 . By utilizing the SVD operation, the channelH e can be decom- is a rank H e × rank H e diagonal matrix of singular values arranged in a decreasing order, and V is a N BS × rank H e unitary matrix, respectively. Furthermore, we also define the unitary matrix U as U = [U 1 , U 2 ], and then the optimal combiner is simply given by W opt = U 1 [49], where U 1 is of dimension N MS × N s . To realize the better communication performance, we calculate the combiner matrix W by using the SVD operation of the cascaded channelH e rather than the IRS-MS chan-nelH 2 , since the cascaded channelH e contains the CSI of joint W and . Thus, the optimization problem (32) can be rewritten as A. EXHAUSTIVE SEARCH SCHEME In terms of the optimization problem (33), the number of alternative matrix is finite since the phase shift of each reflecting element is discrete and is selected from the predefined phase shift set F, where ϕ n ∈ F, ∀n = 1, · · · , N IRS . Under such a condition, the non-convex optimization problem (33) is able to be converted to a phase shift search problem whose optimal solution can be obtained by an exhaustive search method. Specifically, the exhaustive search method demands to successively generate the phase shift value {ϕ n } N n=1 from F, and then the phase shift matrix is able to be determined. Then, with a given , the cascaded channelH e is definite and the corresponding combining matrix W can be obtained by utilizing the SVD operation of the cascaded channelH e . Once the phase shift and the combining matrix is determined, the optimal data rate for the IRS-aided THz communication system can be got by computing all the possible candidates of and W . After traversing VOLUME 8, 2020 all the alternative matrices and the corresponding matrices W , we are able to get the optimal communication rate and the expected phase shift matrix, denoted as opt . During this procedure, there are total |F| N IRS alternative phase shift matrices to participate in calculations. Therefore, the exhaustive search method suffers from severe computational complexity that involves a huge search operations for possible matrices . To decrease the computational complexity, it is pressing to put forward the low complexity methods in practical application scenarios.

B. LOCAL SEARCH SCHEME
To begin with, the main complexity of the exhaustive search method comes from the exponential growth of the possible phase shift combinations when the number of the IRS elements increases. Attractively, different from the exhaustive search method, the proposed local search method can greatly reduce the complexity by just selecting the partially better combinations of discrete phase shifts. Similar to the exhaustive search method, the local search method transforms the optimization problem (33) into the phase shift search problem. The local search method just achieves the suboptimal data rate performance while it can greatly decrease the computational complexity. According to [59], the elements of the diagonal matrix { (i, i)} N s i=1 are regarded as the virtual path gains. and the phase shift matrix is influenced by the virtual path gains. Instead, the combiner W is related with the unitary matrix U 1 rather than the virtual path gains. Subsequently, we turn to maximize the sum of these virtual path gains rather than directly solving the non-convex optimization problem (33). Therefore, the problem (33) can be further transformed as where the definition (a) comes from the linear algebra that for any matrix A ∈ C N s ×N s and Tr From the problem (35) we note that this optimization problem is a phase shift search problem, in which just contains one matrix variable . Compared with exhaustive search method, the proposed local search method owns some obvious advantages. On the one hand, by analyzing the structure features of the matrixH e , the matrix variable W can be eliminated and the obtained problem (35) is not involved with SVD operations. On the other hand, the phase shift search strategy is optimized in the local search method. Instead of calculating all the phase shift combinations, the proposed local search method only takes advantage of the relatively better combinations of the discrete phase shifts. Thus, based on these superiorities, the proposed local search method is capable of lessening the computational complexity tremendously by avoiding the exhaustive search of the phase shift combinations.

Algorithm 2 Local Search Algorithm
Require: BS-IRS channel:H 1 , IRS-MS channel:H 2 , the bit quantization number: b, the controllable phase shift range of each reflecting element: ϕ. 1: Initialize: Phase shift matrix: = 0 N IRS ×N IRS , the phase shift set: = diag e jϕ 1 , · · · , e jϕ N IRS , 6: Calculate the term γ with a given : if γ > γ max do 8: γ max = γ , opt = , 9: end if 10: end for 11 The detailed procedure of the proposed local search algorithm for the data-rate maximization problem is illustrated in Algorithm 2. Here, the important steps of the proposed algorithm are introduced in the following. Firstly, we initialize some key parameters that are utilized in the iterative process, such as the phase shift matrix, the phase shift set for each IRS element and so on. Next, the iterative process is to be described in detail. During the ith iteration process, we fix the phase shifts of reflecting elements from i + 1 to N IRS , and select the optimal phase shift from F for the ith reflecting element by maximizing the term Tr H eH H e . Continue such iteration process for N IRS times, we are able to find the optimal phase shift matrix opt , the optimal combiner W opt and the corresponding data-rate R max in the end.

C. DEEP NEURAL NETWORK BASED SCHEME
Through the analysis of optimization problem as well as the hardware constraints, intriguingly, we note that the phase shift matrix is related to inherent channel features, such as path gain, AoAs and AoDs. In other words, there is a nonlinear functional relationship between the phase shifts of IRS elements and the estimated CSI. However, such an implicit relationship is hard to be described mathematically by the IRS mapping function. Thus, the DNN based method, which maps from the estimated CSI to the phase shift matrix, is an appropriate choice to fit the desired mapping function [60]. The local search method discussed above not only endures obvious performance loss, but also suffers from high computational complexity due to the enormous SVD operations and the matrix multiplication. Therefore, another advantage of our developed method is that the DNN is able to further mitigate the computational burden of the phase shift search problem by fitting the IRS mapping function instead of searching the phase shift combinations of reflecting elements. To this end, the DNN based method is developed in this subsection. From [61] we know that DNN based method well fits the functional relationship when the training samples of the data sets is sufficient enough. Since the estimated CSI in previous section can provide enough data sets, the proposed DNN based method is able to tackle the optimization problem (33) with the low complexity and slight performance penalty compared with the exhaustive search method. Once the DNN based search model is established, the optimal phase shift matrix as well as the corresponding communication rate R can be obtained quickly for the new CSI input. In the following, the neural network architecture of the proposed DNN base method is introduced, as shown in Fig. 5. The feedforward fully connected network is designed here, which contains an input layer, S − 2 hidden layers and an output layer. Also, any two neural units in different layers are interconnected with each other.
In the following, the proposed DNN based phase shift search algorithm is presented in detail. Generally, the proposed DNN based method consists of two distinct phases, namely (I) the offline (training) phase and (II) the online phase.

1) OFFLINE PHASE
In the training phase, the database can be generated by the channel estimation method mentioned above and the exhaustive search method, where the exhaustive search procedure operated online involves tremendous latency overhead. The estimated CSI provides the channel features while the exhaustive search method can find the mapping relationship between the optimal phase shift matrix (optimal data rate) and the channel features. Thus, the input vectorD contains the path gain, AoAs and AoDs. Correspondingly, the output vector O = ϕ 1 , · · · , ϕ N IRS , R is the optimal phase shifts of IRS elements and the desired data rate. By aggregating the input data and the output data together, a complete database is created to train the proposed DNN model. Once the DNN model is trained well, the optimal phase shift matrix can be gained with much lower computational complexity compared with the local search method and the exhaustive search method. It is worth noting that the collected channel samples are generated dynamically during the offline training process. Although the channel varies in practice, the trained DNN model still possesses the extremely high adaptability.

2) ONLINE PHASE
In the online training phase, we suppose that both the database and DNN model are placed on the IRS side. Before the DNN model comes into use, we firstly estimate a new channel sample. With the given CSI, the IRS feeds the preprocessing channel information to the trained DNN model. Finally, the optimal R max and opt can be gained through the online processing of the proposed DNN model. We note that the online prediction only involves simple calculation operations in the neural network, such as the multiplication operation, and thus the proposed DNN based method has the better ability to reduced the computational complexity significantly.
As discussed above, the DNN, as one of the branches of deep learning [60]- [62], has an advantage in fitting the mapping functions between the phase shift matrix and the CSI. With this motivation, the S-layer DNN model is designed to seek the appropriate mapping functions. The critical parts of the DNN model are introduced concretely in the following.

a: FEATURES SELECTION
Through the analysis of optimization objective and the training samples, the array response vectors and the path gains of channel information are finally selected as the input features in the adopted DNN model. In the IRS-assisted THz MIMO communication model, we assume that the BS-IRS channel and IRS-MS channel contain L 1 and L 2 propagation paths, respectively. Letā BS (φ l 1 ),ā IRS (γ AOA l 1 ),ā IRS (γ AOD l 2 ),ā MS (θ l 2 ) represent the estimated array response vectors. Meanwhile, theᾱ l 1 andβ l 2 are the estimated path gains, where l 1 ∈ (1, 2, · · · , L 1 ) and l 2 ∈ (1, 2, · · · , L 2 ).

b: DATA PREPROCESSING
Before feeding the channel samples to the proposed DNN model, the input vectorD needs to be normalized to conclude VOLUME 8, 2020 the uniform statistical distribution of the training samples. In practice, since the THz wave endures serious path loss and molecular absorption loss, the propagation path gains of THz signals changes widely. Thus, the normalization process is indispensable to narrow the gap among the different values of the data samples and accelerate the convergence of the DNN model. After the data preprocessing of the original data sets, the processed data samples can be regarded as the effective input of the DNN model. Next, the formulation of the processed input vectorD can be represented aŝ where D max and D min are the maximum and minimum elements in D .

c: INPUT
For a specific communication scenario, once the number of input features is determined, the dimension of the input layer is also definite. Although the input features of DNN have been selected according to the database, another challenge is that the input vectorD of the proposed DNN has to be the real numbers. To face this difficulty, the complex input vector is essential to be divided into the real part and the imaginary part. Then, the real part and the imaginary part are integrated together to form the final input vector. For the sake of simplifying this explanation, we select the array responsē a BS (φ l 1 ) as an example. Specifically, the vectorā BS (φ l 1 ) of path l 1 is split into the real parts (ā BS (φ l 1 )) and the imaginary parts (ā BS (φ l 1 )). Similarly, all the feature vectors can be split into real and imaginary parts. Therefor, the dimension of the expected input vector is (N BS L 1 + N IRS L 1 + N IRS L 2 + N MS L 2 + L 1 + L 2 ) in the proposed DNN model, and the input vector can be further expressed as It is worth mentioning that the data-rate (R) is correlated with the optimal phase shift matrix , so the proposed DNN model can output the R and the diagonal elements ϕ = (ϕ 1 , · · · , ϕ N IRS ) of the matrix simultaneously. Thus, the final output vector is represented asŌ = ϕ 1 , · · · , ϕ N IRS , R .

e: TRAINING PROCESS OF DNN
On the basis of the working principle of the neural network, each current neural unit sums all the input values which are related to the neural units from the previous layer and the corresponding weights. More specifically, the input value for each current neural unit sums the product of the output values from the previous layer and the corresponding weights. In the hidden layers, the summation values need to be processed by rectified linear units (ReLUs) to improve nonlinear relations and alleviate gradient dissipation problems. Then, the output values of current neural units can be further utilized for the next layer. To avoid the overfitting issue, the dropout units are also designed in the proposed neural network architecture, where the dropout units indicate that a random neural unit is selected as an inactive unit during the learning process.
In order to obtain the optimal and R with small batches of training samples as well as the low computational complexity, the stochastic gradient descent (SGD) algorithm is employed in the proposed neural network architecture. According to the SGD algorithm, the parameters of the DNN model at kth layer can be expressed as where η is the learning rate of the algorithm, is the set of all the neural network parameters, and ∇L (·) is the gradient of the batch's loss function. After that all the parameters are determined, the DNN model is established accordingly. According to the trained DNN model, when a new data sample is input, the predicted outputÔ can be obtained swiftly. To demonstrate the effectiveness of the proposed DNN based method, the minimum mean square error (MMSE) criterion is chosen to minimize the loss function L( ). Also, the MMSE can evaluate the prediction accuracy of the proposed DNN model by calculating the difference between the predicted outputÔ and the sampled outputŌ. During the training process of the DNN model, we need to minimize the regression loss function until the proposed DNN converges. So, the loss function L( ) is expressed as

D. COMPLEXITY ANALYSIS OF DIVERSE SCHEMES
In this subsection, the detailed complexity of different phase shift search schemes are analyzed, including the exhaustive search scheme, local search scheme, and DNN based scheme for the IRS aided THz communication scenarios.
To begin with, the computational complexity of the exhaustive search method mainly caused by three different aspects. Firstly, calculating the multiplication of two matrix Then, the complexity of searching all the possible phase shift matrices is O(N 2 b IRS ). In addition, the complexity of the SVD operation about the matrixH e is O (N MS N BS N s ). Therefore, the total complexity of the exhaustive search scheme is

V. NUMERICAL RESULTS
In this section, numerical results are presented to demonstrate the performance of the IRS-aided THz MIMO system and the conventional THz communication system without IRS in terms of the data rate performance and the computational complexity. Here, the conventional scheme without IRS means that all the system settings are the same as the proposed IRS assisted system except for the existence of IRS, and thus we can see the benefits brought by the IRS. For a given downlink communication scenario, the BS employs N BS = 32 antennas and the MS uses N MS = 32 antennas. As for the IRS, we adopt N IRS = 16 reflecting elements with 2 bits quantization precision and the working frequency of the IRS elements is at 0.22 THz. Particularly, the distance between adjacent antenna elements or the neighbouring reflecting elements is set as the half the wavelength of the THz wave. Due to the extremely sparse nature of the THz channelH 1 andH 2 , we denote that the number of the propagation paths is L 1 =L 3 = 3. Additionally, the complex gains are selected from the circularly-symmetric Gaussian distribution and satisfiesᾱ l 1 ,β l 2 ∈ CN (0, 1), where l 1 , l 2 ∈ [1, 2, 3]. On the basis of the geometric channel model, the physical AoAs and AoDs are randomly generated from [0, 2π] with an equal probability. Ultimately, it worth noting that the channel parameters are estimated in previous section and thus are known here. A. DATA RATE PERFORMANCE COMPARISONS Fig. 6 depicts the data-rate performance comparisons of diverse phase shift search methods along with the increase of the SNR. In the practical THz communication scenarios, it is obvious that all the proposed phase shift search schemes with IRS achieves higher data rate performance than the conventional scheme without IRS. Among these proposed phase shift search schemes, the exhaustive search method performs the optimal performance while the proposed local search method suffers from the apparent performance penalty. Specifically, when the SNR is 10 dB, the exhaustive search method is about 7.59 bps/Hz higher than the local search method. Interestingly, the simulation curve of the proposed DNN based method can be basically consistent with the curve of the optimal method (around 1.14 bps/Hz performance gap at SNR = 10 dB). One may also note that the DNN based method is able to realize the better communication performance compared with the proposed local search method in indoor THz communication systems. Fig. 7 illustrates the data rate performance versus the number of reflecting elements N IRS among different schemes under the condition that SNR = 10 dB. From Fig. 7 we can see that, the data rate performance of the conventional scheme without IRS remains unchanged when the number of the reflecting elements increases. Fortunately, the number of IRS elements has a great impact on the communication performance of our proposed schemes, including the exhaustive search scheme, the local search scheme and the DNN based scheme. In other words, by increasing the number of the IRS elements, all the proposed schemes is able to enhance the data rate performance significantly. Specifically, in the case that the number of IRS elements is configured as 60, the data rate of the proposed DNN based scheme is 38.36 bps/Hz, which is greatly better than that of local search VOLUME 8, 2020  method with 12.24 bps/Hz. Moreover, it is worth mentioning that the data rate performance of the proposed DNN based scheme approaches the simulation curve of the exhaustive search scheme. Fig. 8 presents the data rate performance of different scheme versus the number of training samples. It is obvious that the proposed local search method and the exhaustive search method are not influenced by the number of training samples. Nevertheless, the data rate of the proposed DNN based method becomes continuously better along with the increasing number of the training samples, and has a tendency to approach the the exhaustive search method. Notably, when the number of the training samples is 8000, the DNN based method has almost converged in terms of SNR = 5 dB and SNR = 10 dB. Under the case that the training samples are sufficient at SRN = 5 dB, the DNN based method inherits about 1.17 bps/Hz performance loss compared with the exhaustive search method, and possesses about 3.12 bps/Hz performance enhancement compared with the local search method. In short, as long as the database is large enough, the DNN based method can attain the near optimal data rate performance for the IRS aided THz communication systems.  Fig. 9 investigates the computational complexity comparisons of the proposed three different schemes with the increasing of the IRS elements. For the sake of a fair comparison condition, the proposed diverse phase shift search methods are based on the same parameter settings just as mentioned before. Especially, the proposed DNN model contains four lays where n 2 = 16 and n 3 = 16. Firstly, from Fig. 9 we find that the proposed local search method can greatly reduce the computational overhead in comparison with the exhaustive search method. Moreover, the further complexity reduction is able to be achieved dramatically by the proposed DNN based method compared with the local search method. In more concrete terms, when the number of the IRS elements is 16, the complexity of the exhaustive search, the local search method and the DNN based method is 2.15 × 10 9 , 2.10 × 10 6 and 5.42 × 10 3 , respectively. Generally, the complexity of offline phase is not taken into account here. Additionally, one may point out that, as the number of the reflecting elements increases, the computational complexity of both the local search method and the exhaustive search method grows to a great extent, while the proposed DNN based method varies slowly. Thus, in other words, the advantage of the proposed DNN based method method in terms of the complexity tends to be more distinct. Ultimately, combining the Fig. 6 and Fig. 9, the proposed DNN based scheme is capable of making a better trade-off between the data rate performance and the computational complexity, which will be widely deployed in future IRS assisted THz communication scenarios.

VI. CONCLUSION
THz frequency band is regarded as a promising alternative to provide large bandwidth and high data rates for 6G wireless communication systems. To settle the poor coverage capability and high path loss of THz waves existing in indoor THz communication systems, the IRS has been proposed recently to enhance the communication rates of the terminal equipment as well as the coverage capacity. By utilizing the CS technique, we put forward a low complexity IAP-SP scheme to solve the channel estimation problem that involves the double grid quantization issue. Fortunately, the proposed IAP-SP algorithm achieves the basically same recovery performance compared with classical SP algorithm. With the estimated CSI, we then propose three diverse data rate maximization methods by searching the phase shifts of IRS elements. Remarkably, the proposed DNN based scheme makes a better tradeoff between computational complexity and data rate performance. Specifically, as long as the training samples are large enough,the proposed DNN-based scheme is able to achieve the near-optimal performance and debases the computational burden sufficiently. It is also worth noting that these proposed frameworks in this paper are able to be extended to multiuser MIMO communication scenarios. Ultimately, the IRS-assisted THz communication systems will act an important role in indoor 6G application scenarios.

B. PROOF OF (28)
Considering the (27), we note that . . . . . . where Based on the fact thatḡ 1 i 1 ,j 1 is a complex number, so we can get Thus, the proof for (28) is completed.
C. PROOF OF (29) On the basic of (47), the entry ofḠ 1 can be formulated as whereḡ 2 i 2 ,j 2 ∈Ḡ 2 , a i 1 ,j 1 (i 2 , j 2 ) ∈ A i 1 ,j 1 , ∀ i 1 , j 1 ∈ [1, N 1 ], ∀ i 2 , j 2 ∈ [1, N 2 ]. By resorting to (49), the proof for (29) can be completed directly. In 1984, he joined UESTC as an Academic Member, where he has been a Professor of information and communication systems, since 1997, and also been a Ph.D. Supervisor, since 2000. He is currently the Director of the National Key Laboratory of Science and Technology on Communications, UESTC. He holds over 60 granted and filed patents. His general interests include wireless and mobile communications, anti-jamming technologies, and signal processing for communications subjects, on which he has authored or coauthored over 100 journal articles, 100 conference papers, and two edited books. His current research interests include multiple-antenna signal processing technologies for mobile communications, cognitive radios, and coding and modulation for next-generation mobile broadband communications systems. He has been a member of the  VOLUME 8, 2020