Beamspace Channel Estimation for Millimeter Wave Massive MIMO System with Hybrid Precoding and Combining

In this paper, a framework of beamspace channel estimation in millimeter wave (mmWave) massive MIMO system is proposed. The framework includes the design of hybrid precoding and combining matrix as well as the search method for the largest entry of over-sampled beamspace receiving matrix. Then based on the framework, three channel estimation schemes including identity matrix approximation (IA)-based scheme, scattered zero off-diagonal (SZO)-based scheme and concentrated zero off-diagonal (CZO)-based scheme are proposed. These schemes together with the existing channel estimation schemes are compared in terms of computational complexity, estimation error and total time slots for channel training. Simulation results show that the proposed schemes outperform the existing schemes and can approach the performance of the ideal case. In particular, total time slots for channel training can be substantially reduced.

Abstract-In this paper, a framework of beamspace channel estimation in millimeter wave (mmWave) massive MIMO system is proposed. The framework includes the design of hybrid precoding and combining matrix as well as the search method for the largest entry of over-sampled beamspace receiving matrix. Then based on the framework, three channel estimation schemes including identity matrix approximation (IA)-based scheme, scattered zero off-diagonal (SZO)-based scheme and concentrated zero offdiagonal (CZO)-based scheme are proposed. These schemes together with the existing channel estimation schemes are compared in terms of computational complexity, estimation error and total time slots for channel training. Simulation results show that the proposed schemes outperform the existing schemes and can approach the performance of the ideal case. In particular, total time slots for channel training can be substantially reduced.

I. INTRODUCTION
Millimeter wave (mmWave) communication is a promising technology for next generation wireless communication owing to its abundant frequency spectrum resource [1], [2]. However, realizing mmWave massive MIMO in practice is not a trivial task, which faces the problem of high propagation loss due to the high carrier frequency [3]. To compensate for the propagation loss, antenna arrays are usually used to form directional beamforming. Fortunately, thanks to the short wavelength of the mmWave frequency, large antenna arrays are possible to be packed into small form factors.
In order to exploit the spatial degree of freedom, the hybrid analog and digital precoding is usually employed [4]- [6]. A small number of RF chains are tied to a large antenna array. This structure enables parallel transmission, and thus provides the potential to approach the capacity bound that can be achieved by digital precoding. On the other hand, the large antenna arrays challenge the low-complexity design of hybrid precoding and channel estimation [7]. In particular, the hybrid precoding may require matrix operations with a scale of antenna size, which is generally large in mmWave communication [4]. Moreover, the channel estimation is also rather time consuming due to the large number of antennas at both transmitting and receiving sides [8]. To reduce the complexity of channel estimation in mmWave massive MIMO system, some advanced schemes based on the beamspace channel have been proposed very recently [8]- [14]. A hybrid and multiresolutional codebook (HMC)based channel estimation method and a JOINT-based channel estimation method are proposed in [8] and [9], respectively. The common idea of [8] and [9] is to use the hierarchical codebook, where the precision of channel estimation relies on the number of the layers in the hierarchical codebook.
The key ideas of [11]- [14] are to efficiently explore the sparsity of beamspace channel by sparse signal processing techniques. An enhanced compressive sensing (ECS)-based channel estimation scheme is proposed in [11], which explores the channel sparsity in beamspace. With higher resolution of phase shifters, improved channel estimation accuracy can be achieved since higher freedom is available for the design of measurement matrix for the sparse recovery [4], [5].
However, considering the limited beamspace resolution, the sparsity of beamspace channel may be impaired by power leakage [15], indicating that the beamspace channel is not ideally sparse and there are many small nonzero entries [16]. Therefore, it brings extra challenge for the sparse recovery [17], [18]. To solve this problem, an adaptive support detection (ASD)-based channel estimation scheme is proposed to iteratively detect and adjust the channel support to find one with the largest channel power [14]. Note that the beamspace resolution is proportional to the reciprocal of the antenna number. Therefore, ASD-based scheme cannot obtain high precision channel estimation regarding the channel power leakage [8]. A discrete compressive sensing (DCS)-based channel estimation scheme is proposed in [15], where the beamspace resolution can be set arbitrarily and can be higher than the reciprocal of the antenna number, leading to better channel estimation performance than ASD-based scheme. An over-sampling compressive sensing (OCS)-based channel estimation scheme is proposed in [17], where the measurement matrix is consisted of a large number of over-sampling steering vectors and is capable of estimating the angle of arrival (AoA) and angle of departure (AoD) more accurate than conventional measurement matrix. However, both DCS-based and OCSbased schemes use random measurement matrices, which can not always achieve the optimal performance. In particular, the estimated AoA and AoD can not be always within the scale of quantization error even without any noise. Therefore, it is better to use deterministic measurement matrix.
In this paper, we first propose a framework of beamspace channel estimation, which is divided into two subproblems, including the design of the hybrid precoding and combining matrix, and the search method for the largest entry of oversampled beamspace receiving matrix. Note that the design of the hybrid precoding and combining matrix as well as the search method for the largest entry is original. Then based on the framework, we propose three channel estimation schemes. 1) We propose an identity matrix approximation (IA)-based channel estimation scheme. We formulate the design of hybrid combining and hybrid precoding as two optimization problems with the constraint of total power and constraint of the constant envelope required by phase shifters. Due to the non-convexity of the problems, we decouple the design of analog combining and digital combining, and then obtain closed-form solutions. We propose an algorithm for the design of hybrid combining matrix. The algorithm repeatedly fixes the analog combining matrix to obtain digital combining matrix, and then fixes the digital combining matrix to obtain the analog combining matrix in turn, until the stop condition is satisfied. Detailed steps summarized in an algorithm table are provided. Since the design of the hybrid precoding matrix is similar, we briefly describe the design of the hybrid precoding matrix. Note that the hybrid combining matrix and the hybrid precoding matrix can be designed off-line before the channel training. After that, we propose an algorithm to search the largest entry of the over-sampled beamspace receiving matrix. The algorithm is based on trichotomy search and includes two stages, where we find the main lobe in the first stage, and we find the largest entry corresponding to the peak within the main lobe in the second stage.
2) We design hybrid precoding matrix and hybrid combining matrix so that the coordinates of the largest entry of oversampled beamspace receiving matrix are the AoA and AoD within the quantization error. Then we convert the discrete problem into a continuous problem to obtain the derivative. By setting the derivative zero, we get two solutions where the first solution is found to be meaningless. Based on the second solution, we propose two zero off-diagonal (ZO) beamspace channel estimation schemes, namely, scattered zero off-diagonal (SZO)-based scheme and concentrated zero offdiagonal (CZO)-based scheme.
• In the SZO-based channel estimation scheme, the nonzero diagonal entries are uniformly distributed with the same interval. Since the main lobe and the side lobes have the same envelope, we first make beam training based on codebook to find the main lobe and then use complementary channel estimation to further estimate the channel AoA and AoD within the main lobe. An integrationbased codebook design method which results in closedform expression of codewords is proposed and compared with the existing sparse-based codebook design method. Note that the proposed integration-based codebook design method is original. • In the CZO-based channel estimation scheme, the nonzero diagonal entries are concentrated on the upper left corner of the matrix. Since the envelope of the main lobe and the side lobe is different, we can directly employ the algorithm proposed in the IA-based scheme to search the largest entry of over-sampled beamspace receiving matrix corresponding to the channel AoA and AoD within the main lobe. Additionally, we also compare the above three schemes together with the existing channel estimation schemes in terms of computational complexity, estimation error and total time slots for channel training.
The rest of the paper is organized as follows. In Section II, the system model and problem formulation of beamspace channel estimation with hybrid precoding and combining are provided. In Section III, we propose three beamspace channel estimation schemes. The simulation results are provided in Section IV. Finally, Section V concludes the paper.
The notations are defined as follows. Symbols for matrices (upper case) and vectors (lower case) are in boldface.
, · , ⌊·⌋, Tr(·), Z, R{·} and CN , denote the transpose, conjugate transpose (Hermitian), inverse, identity matrix of size L, vector of size L with all entries being 1, the set of M × N complex-valued matrices, kronecker product, vectorization, expectation, order of complexity, zero vector of size M , M × N zero matrix, l 0norm, l 2 -norm, Frobenius norm, entry of A at the pth row and qth column, round function, floor function, trace, set of integer, real part and complex Gaussian distribution, respectively.

A. System Model
We consider a time division duplexing (TDD) multi-user mmWave massive MIMO system comprising a base station (BS) and U users. We focus on the uplink transmission. Both the BS and users are equipped with an uniform linear array (ULA) [19]. Let N A , M A , N R and M R denote the number of antennas at the BS, number of antennas at each user, number of RF chains at the BS and number of RF chains at each user. In practical mmWave massive MIMO with hybrid precoding and combining, the number of RF chains is much smaller than that of antennas, i.e., N R ≪ N A and M R ≪ M A . This is because we use large antenna array to form directional beamforming, which can compensate for the propagation loss caused by the high carrier frequency [20].
For uplink transmission, each user performs analog precoding in RF and digital precoding in the baseband, while the BS performs analog combining in RF and digital combining in the baseband. The received signal vector at the BS can be represented as where F B,u ∈ C MR×MR , F R,u ∈ C MA×MR , W B ∈ C NR×NR , and W R ∈ C NR×NA are the digital precoding matrix, analog precoding matrix, digital combining matrix, and analog combining matrix for the u(u = 1, 2, . . . , U )th user, respectively. s u ∈ C MR denotes the signal vector satisfying E{s u s H u } = I MR . n ∈ C NA denotes additive white Gaussian noise (AWGN) vector satisfying n ∼ CN (0, σ 2 I NA ).
H u ∈ C NA×MA denotes the channel matrix between the BS and the uth user and can be expressed according to the widely used Saleh-Valenzuela channel model [1] as where L u and g u,i denote the total number of resolvable paths and the channel fading coefficient of the ith path for the uth user, respectively. Usually, there is a strong line-of-sight (LOS) path (i = 1) and L u − 1 much weaker non-line-ofsight (NLOS) paths (2 ≤ i ≤ L u ). The mmWave transmission essentially relies on the strong LOS path. The steering vector α(N, θ) is defined as Define the AoA and AoD of the ith path of the uth user as Θ u,i and Φ u,i , respectively. Further define θ u,i 2dBS λ sin Θ u,i and ϕ u,i 2dUE λ sin Φ u,i , where d BS and d UE denote the antenna interval of the BS and users, respectively. We usually set d BS = d UE = λ/2, where λ is the wavelength of mmWave signal. In practice, both Θ u,i and Φ u,i obey the uniform distribution [−π, π] [8], [21].

B. Problem Formulation
Note that y in (1) is a combination of signal from different users. We use T 1 different digital precoding matrices and analog precoding matrices, denoted as F t1 B,u ∈ C MR×MR and F t1 R,u ∈ C MA×MR , respectively, t 1 = 1, 2, . . . , T 1 , at the u(u = 1, 2, . . . , U )th user. We use T 2 different digital combining matrices and analog combining matrices, denoted as W t2 B ∈ C NR×NR and W t2 R ∈ C NR×NA , respectively, t 2 = 1, 2, . . . , T 2 , at the BS. To distinguish different user signal at the BS, each user repeatedly transmits an orthogonal pilot sequence p u ∈ C U for T 1 T 2 times. For simplicity, we suppose each user transmit the same pilot sequence for all M R RF chains, where the pilot matrix for the uth user can be defined as P u [p u , p u , . . . , p u ] H = 1 MR p H u ∈ C MR×U . The channel keeps constant during T T 1 T 2 U time slots [16]. During the T 1 repetitive transmission of pilot sequence from the ((t 2 −1)T 1 +1)th transmission to (t 2 T 1 )th transmission, we use T 1 different F t1 B,u and F t1 R,u for hybrid precoding while using the same W t2 B and W t2 R for hybrid combining, where the received pilot matrix Y t1,t2 ∈ C NR×U can be denoted as with N t1,t2 ∈ C NA×U representing the AWGN matrix. Each entry of N t1,t2 independently obeys complex Gaussian distribution with zero mean and variance of σ 2 . To ease the notation, we define N t1,t2 W t1 B W t1 R N t1,t2 . Due to the orthogonality of p u , i.e., p H u p u =1 and p H u p i = 0, ∀u, i ∈ {1, 2, . . . , U }, i = u [14], we can obtain the measurement vector r t1,t2 u ∈ C NR for the uth user by multiplying Y t1,t2 with p u as where Note that we can distinguish different user signal utilizing the orthogonality of pilot sequences. In this way, the spatial distribution of users has no impact on the performance of our proposed schemes. Define T 3 T 2 N R . We stack the T 2 received pilot sequences together and have where Further define We have Considering the LOS channel path concentrates most channel power in mmWave massive MIMO systems, we usually use LOS channel path to transmit data for the uth user [15]. Therefore, it is important to design W and F u to estimate the AoA and the AoD of the LOS path of H u in (10) for the uth user, which will be discussed in the following sections.

III. BEAMSPACE CHANNEL ESTIMATION
In this section, we first propose a framework of beamspace channel estimation. Then based on this framework, three channel estimation schemes are proposed. Finally, the comparisons of these three schemes together with the existing HMCbased [8], JOINT-based [9], ECS-based [11], DCS-based [15] and OCS-based [17] channel estimation schemes are also presented.

A. Framework of Beamspace Channel Estimation
The beamspace channel matrix for the uth userH v u ∈ C NA×MA can be represented as [14] H where D(N, K) ∈ C N ×K is the sampling matrix, which is defined as In fact, D(N, K) samples the beamspace [-1,1] in an interval of 2/K by K steering vectors. ForH v u , the AoA and AoD is sampled in an interval of 2/N A and 2/M A , respectively. Therefore the quantization error for the estimated AoA and AoD is 2/N A and 2/M A , respectively. In order to decrease the quantization error, we introduce the over-sampled beamspace channel matrix for the uth user H v u ∈ C K×K as where K is the number of steering vectors with K > N A and K > M A . Then the coordinates of the largest entry of H v u are the AoA and AoD of the LOS path with the quantization error of 2/K. To reduce the quantization error and improve channel estimation, we can use a large K by finding the largest entry of H v u . However, we cannot directly obtain H v u based on R u in (10), due to the hybrid precoding and combining operations of F u and W , respectively. Note that the dimension of In order to obtain the over-sampled beamspace channel matrix as described in (14), we multiply R u with W H on the left and F H u on the right, which can make the dimension of W H W H u F u F H u the same as that of H u . Now we can obtain an over-sampled beamspace receiving matrix However, in this case, it requires that T 3 ≥ N A and T 1 ≥ M A , leading to huge pilot overhead. For example, in an mmWave massive MIMO system with N A = 64, M A = 16, N R = 4 and M R = 1, the pilot sequence should be repetitively transmitted for T 1 T 2 = 256 times. Therefore, it is required to reduce the pilot overhead in practice, where can be satisfied. Now we have the following two subproblems.
Subproblem 1 (Hybrid Precoding and Combining Matrix Design): W and F u should be well designed so that the coordinates of the largest entry of R v u are the AoA and AoD of the LOS path with the quantization error of 2/K. Therefore, the AoA and AoD of the LOS path can be estimated by finding the largest entry of R v u . Subproblem 2 (Search the Largest Entry): Due to the large dimension of R v u , finding the largest entry of R v u is computationally expensive. Therefore, it is better to design a low-complexity search algorithm fully regarding the structure of R v u . (15) is the over-sampled beamspace channel matrix with noise if (16) is satisfied. However, due to the fact that T 3 < N A , T 1 < M A , the rank of W H W and F u F H u is less than N A and M A , respectively. So we cannot find a proper W and F u satisfying (16).

B. IA-based Beamspace Channel Estimation Scheme 1) Hybrid Precoding and Combining Matrix Design:
The optimization problem for hybrid combining matrix design is where W R is the set of all feasible analog combining matrix and P W = 1 to normalize the hybrid combining matrix. Note that the design of hybrid precoder and hybrid combiner is independent. We can well design the hybrid precoding and combining matrices before the transmission of pilot sequences. Note that D( We have (11). Then (17) can be further rewritten as It is seen that W defined in (8) is a flat matrix where the columns are more than the rows. Therefore, it is infeasible that W H W equals γ N I NA . Instead, it is important to design W so that W H W /γ N approximates the identity matrix, a.k.a, identity matrix approximation (IA). To minimize [22]. For example, we obtain U by singular value decomposition (SVD) of a N A × N A random matrix A, i.e., A = U ΣV H , where each entry of A obeys the uniform distribution [0,1]. In this way, we obtain W . According to (8), we can obtain W t2 , t 2 = 1, 2, . . . , T 2 , which is essentially dividing W into T 2 submatrices. Then (18) is converted into T 2 subproblems, where each subproblem can be expressed as We cannot directly obtain solutions for (19) due to the nonconvexity of the constraints. Note that [23] has proved that temporarily neglecting the second power constraint of (19) during the optimization of W t2 B and W t2 R will have little impact on the optimality of the hybrid precoding problem. After W t2 B and W t2 R are obtained, we may set W t2 B as to satisfy the second constraint of (19). To mitigate the interference among different data streams, we impose a common constraint that the columns of the digital combining matrix are mutually orthogonal, i.e., W t2 Therefore (19) can be expressed as It shows that W t2 R and W t2 D are decoupled. Given W t2 D , the solution of W t2 R can be expressed as where arg(A, R) first normalizes each entry of A and then quantizes the normalized matrix in terms of R. Note that the quantization is required since the resolution of phase shifters is limited in practice, e.g., the resolution is 2/64 if the phase shifter is 6 bits and the range of the angle is [-1,1]. Similarly, given W t2 R , the optimization of W t2 D based on (22) can be expressed as It is seen that (24) is similar to the orthogonal Procrustes problem [24]. Then the solution to (24) can be obtained as Note that both (23) and (25) are optimal closed-form solutions to the problems expressed in (22) and (24), respectively.
As shown in Algorithm 1, we propose an algorithm of hybrid combining matrix design for the IA-based channel estimation. We repeatedly fix W t2 R to obtain W t2 B via (25), and then fix W t2 R to obtain W t2 B via (23) in turn. Define the normalized iteration error ǫ as Obtain W t2 via (6). 12: end for 13: Obtain W based on W t2 , t 2 = 1, 2, . . . , T 2 via (8). 14: Output: W .
the stop condition is that the iterative update of both W t2 R and W t2 B is stable, i.e., ǫ < δ, where δ is the threshold. The optimization problem for hybrid precoding matrix design is where F R is the set of all feasible analog precoding matrix and P F is the given power for the hybrid precoding. Similar to (18), (27) can be rewritten as We can obtain a solution as F u by selecting the first T 1 columns of According to (9), we can obtain f t1 u as the t 1 -th column of Then we have Similar to (19), (28) can be converted to T 1 subproblems, where each subproblem is expressed as Similar to (19), the second constraint of (31) can also be temporarily neglected. Therefore, we may replace W t2 In order to run Algorithm 1 to obtain F u , we have to further replace N A , T 2 , W and W with M A , T 1 , F u and The routine of the algorithm is exactly the same, except that an additional operation to obtain F t1 B,u as is required after finishing step 9; and W t2 B should be replaced by (F t1 B,u ) H at step 10. In summary, in the first half of the IA-based channel estimation, we design of hybrid combining matrix and hybrid precoding matrix; while in the other half to be discussed, we will search the largest entry of the over-sampled beamspace channel matrix.
2) Search the Largest Entry: Instead of exhaustively search the largest entry from R v u defined in (15), we can improve the efficiency of the search algorithm by analyzing the structure of R v u . Neglecting the term from the additive noise and assuming there is single path, we rewrite R v u as Since r N is a column vector and r H M is a row vector, the largest entry of R v u essentially depends on the largest entry of r N and r M . Therefore we will analyze the structure of r N and r M .
In the ideal case, i.e., F u F H u D(M A , K) = γ M D(M A , K), r M is an over-sampled transmit steering vector of α(M A , ϕ u,1 ) with an interval of 2/K. As shown in Fig. 1, we illustrate the envelope of r M , where the position corresponding to the peak of r M is the channel AoD with quantization error of 2/K. In order to apply fast algorithms such as Obtain Q 1 , Q 2 , Q 3 and Q 4 via (41). 10: Obtain Q min via (43). 11: Update Γ and Υ. 12: end while 13:θ u,1 = Γ 1 ,φ u,1 = Υ 1 14: Output:θ u,1 ,φ u,1 . trichotomy search to find the peak of the curve, we have to first find the main lobe with the width of 4/M A ; otherwise these algorithms may stop the search at the peak of side lobes. Now we propose Algorithm 2 to fast search the largest entry of R v u , considering the structure of the steering vectors. Algorithm 2 includes two stages. In the first stage, we find the main lobe of r M and r N without oversampling. In the second stage with oversampling, we apply the trichotomy search to find the peak of the main lobe.
In the first stage from step 3 to step 6, we search the main lobe, which is formed by two adjacent columns and two adjacent rows ofH v u defined in (12) [25]. We obtain the beamspace receiving matrixR v u ∈ C NA×MA at step 3 as R is a noise term. We first find two adjacent columns indexed by {s q , s q + 1} with the largest channel power fromH whereR v u,q ∈ C NA×2 represents a submatrix consisted of two consecutive columns ofR whereR v u,p ∈ C 2×MA represents a submatrix consisted of two consecutive rows ofR v u , with row indices denoted as p and p + 1. By finding out the largest two adjacent columns {s q , s q + 1} and rows {s p , s p + 1} of beamspace channel matrixH v u , the search of AoA and AoD can be limited to the range of respectively, where In this way, we narrow down the search space of the AoA and AoD from [−1, 1] to Γ and Υ, respectively. As shown in Fig. 2, the two adjacent rows and two adjacent columns found in the first stage are illustrated in a grey square area.
In the second stage from step 8 to step 14, we find the coordinates of the largest entry of R v u corresponding to the AoA in Γ and AoD in Υ. Note that the quantization error of the AoD and AoA is reduced from 2/M A and 2/N A to both 2/K by oversampling. We apply the trichotomy search. The two points that divide Γ into three equal parts are 2Γ 1 /3 + Γ 2 /3 and Γ 1 /3 + 2Γ 2 /3. Similarly, the two points that divide Υ into three equal parts are 2Υ 1 /3+Υ 2 /3 and Υ 1 /3+2Υ 2 /3. These four points marked in light blue in Fig. 2 can divide the area of R v u corresponding to Γ and Υ into nine smaller areas. The entries corresponding to these four points are where quan() is the quantization function to quantize the consecutive θ into K discrete points, which is defined as Then we compare the amplitude of Q 1 , Q 2 , Q 3 and Q 4 to find the smallest one, which is expressed as We delete the parts that include Q min and update Γ and Υ. As shown in Fig. 2, Q min is marked in dark blue, where the entries within the area of the updated Γ and Υ are marked in dark grey. For example, if Q min = Q 1 , the updated Γ 1 , Γ 2 , Υ 1 and Υ 2 , denoted asΓ 1 ,Γ 2 ,Ῡ 1 andῩ 2 , respectively, can be represented as We repeat the procedures until Γ 2 − Γ 1 ≤ 2/K and Υ 2 − Υ 1 ≤ 2/K, which means the resolution 2/K of over-sampling search is reached for both the AoA and AoD. Finally we output the estimated AoA and AoD of LOS path at step 15.

C. ZO-based Beamspace Channel Estimation Schemes
In the IA-based beamspace channel estimation, we first solve two optimization problems (17) and (27) to design the hybrid combining matrix and hybrid precoding matrix. However, considering that W H W does not exactly equal γ N I NA and F u F H u does not exactly equal γ M I MA , there are possible errors for the search of the largest entry even without noise. Also note that the off-diagonal entries of W H W and F u F H u are not restricted to be zero, which may introduce some interference for beamspace channel estimation.
In this subsection, we design hybrid precoding matrix F u and hybrid combining matrix W so that the coordinates of the largest entry of R v u are the AoA and AoD of the LOS path with the quantization error of 2/K for single path. Since the largest entry of R v u essentially depends on the largest entry of r N in (34) and r M in (35), respectively, we will first analyze r N and r M .
The kth entry of r M can be represented as r M [k] = α H (M A , ϕ u,1 )F u F H u α(M A , −1 + 2(k − 1)/K). Therefore, F u should be designed so that the largest entry of r M is r M [quan(ϕ u,1 )], which indicates the quantization error is 2/K. However, it is difficult to solve this problem for discrete K points. Now we convert this discrete problem into continuous problem to obtain the derivative. To ease the notation, we define α M (ϕ) α(M A , ϕ) and ϕ g ϕ u,1 . We further define a function of ϕ as R(ϕ g , ϕ) which can be regarded as the continuous version of r M . The problem can be formulated as where ϕ g is the genuine channel parameter defined in (2) but is unknown to the receiver. It is expected that one of the solutions to (45) is ϕ g .
which is a diagonal matrix with M A − 1 nonzero entries. Then (46) can be further written as Since ϕ g is one of the solutions to (45), we have Then we have (51) There are two solutions of F + u to (51), denoted as Solution 1 and Solution 2. while In terms of Solution 1, if (52) is satisfied, we set ϕ g = −1, −1 + 2/M A , . . . , −1 + 2(M A − 1)/M A and obtain M A equations. We combine these M A equations together, having where A 1 ∈ C MA×MA is a matrix with zero diagonal entries, leading to Tr(A 1 ) = 0.
But the trace of the expression on the left side of (54) is where (a) is true because F + u is positive semi-definite, and the equality of (a) holds only when F + u = 0 MA . Simultaneously satisfying (55) and (56) leads to However, in practice F + u can not be zero matrix due to the power constraint in (27). Therefore Solution 1 is meaningless.
In terms of Solution 2, we first rewrite the expression on the left side of (53) as Therefore (53) can be further written as where ξ is a constant. Similar to (54), we have where A 2 ∈ C MA×MA is a matrix whose diagonal entries all equal ξ, leading to The trace of the expression on the left of (60) is Tr F + u /M A according to (56). Combining (61) and (56), we have Now we solve (59) to obtain the solution of F + . The expression on the left side of (59) is essentially in a quadratic form, which can be written in detail as Based on (59), (62) and (63), we have Since in practice ϕ g can be any value in [-1,1], F + u should satisfy which means the summation of off-diagonal entries on the same line parallel to the main diagonal line of F + u is zero. To simplify the problem, we set all off-diagonal entries of F + u to be zero, a.k.a., zero off-diagonal (ZO). In this way, we have to shut down some antenna ports for ZO-based schemes, while the number of active RF chains keeps the same. In fact, according to the existing literature [26]- [28], some antenna ports are shut down in mmWave massive MIMO systems. In [26] and [27], some antenna ports are shut down to achieve the antenna selection aiming at the sum-rate maximization. In [28], some antenna ports are shut down to achieve wide mainlobe of beams for mmWave beam training based on hierarchical codebooks. Since F u ∈ C MA×T1 is a tall matrix with M A > T 1 and F + u = F u F H u , the rank of F + u is no more than T 1 , indicating that there are only T 1 nonzero entries of the diagonal matrix F + u . We set the T 1 nonzero entries to be the same γ M M A /T 1 and the left M A − T 1 diagonal entries to be zero.
Given F + u , we can obtain F u by making SVD of F + u , i.e., is an unitary matrix and Σ M ∈ C T1×T1 is a real diagonal matrix. Then we can obtain F u by We set W + to be a diagonal matrix, where T 3 nonzero diagonal entries are set to be the same γ N N A /T 3 and the left N A − T 3 diagonal entries are zero. Given W + , we can obtain W by making SVD of W + , i.e., W + = U N Σ N U H N , where U N ∈ C NA×T3 is an unitary matrix and Σ N ∈ C T3×T3 is a diagonal matrix. Then we can obtain W by According to different layout of nonzero diagonal entries, now we propose a scattered zero off-diagonal (SZO) and a concentrated zero off-diagonal (CZO) based beamspace channel estimation scheme. In the SZO scheme, the nonzero diagonal entries of F + u and W + are uniformly distributed with the same interval. In the CZO scheme, the nonzero diagonal entries of F + u and W + are concentrated on the upper left corner of the matrix.
1) SZO-based Beamspace Channel Estimation Scheme: We design F + u as where Z N denotes an N × N matrix where only the entry on the upper-left corner is one and all the other entries are zero. We design W + as (70) Then we analyze the largest entry of r M . The amplitude of r M [k] can be derived by (71). It is verified by (71) that the off-diagonal entries of F + u are all zero. We observe that r M [k] is essentially the inner product between a steering vector α(T 1 , M A (−1+2(k −1)/K)/T 1 ) and a steering vector α(T 1 , ϕ g M A /T 1 ), indicating that r M [k] is maximized when the angles of these two steering vectors are the closest, i.e., (73) Note that the angles of these two steering vectors may be out of [−1, 1] because M A /T 1 > 1. Therefore we need to add the term of 2l in (73) to guarantee that (−1 + 2 It is seen that the number of solutions to (74) is ⌊M A /T 1 ⌋. Define ϕ g −1 + 2( k − 1)/K as an estimation of ϕ g . We have Suppose T 1 K/M A to be an integer. Then (75) can be rewritten as where 2 (ϕ g + 1)K/2 /K − 1 is essentially the quantization of ϕ g with resolution of 2/K. It seen from (76) that ϕ g is periodic with l, indicating that the main lobe and the side lobes of |r M | have the same envelope, which is illustrated in the third sub-figure of Fig. 1.
Similarly, we define θ g θ u,1 to ease the notation, and further define θ g as an estimation of θ g , we have It seen from (77) that θ g is periodic with i.
In order to eliminate the uncertainty of l in (76) and i in (77), we resort to beam training based on codebook [29] to find the main lobes of |r M | and |r N |, and then design F u and W to estimate the AoA and AoD of the LOS path within the main lobe.
(i) Codebook Design and Beam Training: As shown in Fig. 3, a typical hierarchical codebook has S layers, which satisfy N = 2 S , where N is the number of antennas. In the s(s = 1, 2, . . . , S)th layer, there are 2 s codewords with the same beam width but different steering angles. The union of beam angle of all the codewords in each layer is [−1, 1].
Denote c(s, n) as the n(n = 1, 2, . . . , 2 s )th codeword in the sth layer, covering the beam angle where ω 1 (s, n) −1 + (n − 1)/2 s−1 and ω 2 (s, n) −1+n/2 s−1 . Now we propose an integration-based codebook design method. The codeword c(s, n) is the integration of steering vectors α(N, θ) with θ from ω 1 (s, n) to ω 2 (s, n), which is expressed in (79). In Fig. 4, we compare the integration-based codebook design method with the existing sparse-based codebook design method [8], with respect to the codewords in the first and second layer of the codebook. It shows the codewords designed by the integration-based method have the same beam pattern as those designed by the sparse-based method. However, our method has a closed-form expression, which is much easier for beam generation than the sparse-based method.
After transmitting the beams formed by the codewords of the Sth layer, the AoA and AoD can be coarsely estimated (ii) Complementary Channel Estimation: After the codebook training in (i), the main lobes of |r M | and |r N | are identified. However, the estimation precision of the AoA and AoD is not satisfied. Therefore, we make complementary channel estimation after the codebook training. During the complementary channel estimation, limited pilot transmission with small pilot overhead is needed.
As illustrated in the third sub-figure of Fig. 1, the width of the main lobe of |r M | is 4/M A , which requires the period of 2T 1 /M A no smaller than 4/M A , resulting in T 1 ≥ 2.
Similarly, we have T 3 = T 2 N R ≥ 2, resulting in T 2 ≥ 1. We set T 1 = 2 and T 2 = 1, meaning that we only need two different hybrid precoding matrices F u and one hybrid combining matrix W , where the overhead of pilot training is substantially reduced compared to the IA-based scheme.
Finally, we run the second stage of Algorithm 2 with Γ 1 −1 + 2(n W − 3/2)/N A , Γ 2 −1 + 2(n W − 1/2)/N A , Υ 1 −1+2(n F −3/2)/M A and Υ 2 −1+2(n F −1/2)/M A . Note that we cannot run the first stage of Algorithm 2 to find the main lobe of |r M | and |r N |, since the main lobe and the side lobes have the same envelope, as shown in the third sub-figure of Fig. 1. We can only use the beam training based on codebook to find the main lobe.
2) CZO-based Beamspace Channel Estimation Scheme: In this scheme, the nonzero diagonal entries of F + u and W + are concentrated on the upper left corner of the matrix. We design F + u and W + as (81) Then we analyze the largest entry of r M . Similar to (71), the amplitude of r M [k] can be derived by (72). We also observe that r M [k] is the inner product between a steering vector α(T 1 , −1 + 2(k − 1)/K)) and a steering vector α(T 1 , ϕ g ), indicating that r M [k] is maximized when the angles of these two steering vectors are the closest, i.e., Define ϕ g −1 + 2( k − 1)/K as an estimation of ϕ g . We have ϕ g = 2 ϕ g + 1 K/2 /K − 1.
where 2 (ϕ g + 1)K/2 /K − 1 is essentially the quantization of ϕ g with resolution of 2/K. Unlike (75), ϕ g in (84) is not periodic. Therefore, the envelope of the main lobe and the side lobe is different, as shown in the last sub-figure of Fig. 1, where the beam training based on codebook to identify the main lobe is not necessary. We can directly employ Algorithm 2 to find the largest entry, which corresponds to the AoD and AoA of the LOS path.

D. Comparisons
Now we compare the proposed three schemes together with the HMC-based, JOINT-based, ECS-based, DCS-based and OCS-based channel estimation scheme in terms of computational complexity, estimation error and total time slots for channel training, which are summarized in Table I. 1) Computational Complexity: As shown in (65), (67) and Algorithm 1, the proposed hybrid precoding and combining schemes are independent of the channel matrix, so we can design the hybrid precoding and combining matrix well before the transmission of pilot sequences. In this way, we do not need to consider the channel coherence time when designing the hybrid precoding and combing matrix. In other words, the computational complexity mainly comes from the search of the largest entry in Algorithm 2. For both the IA-based scheme and CZO-based scheme, in the first stage, we need to compute Frobenius norm of the N A × 2 matrixR For the SZO-based scheme, where only the second stage of Algorithm 2 is needed, we use the trichotomy search to find the largest entry among K/N A × K/M A entries, resulting in the total computational complexity to be For the DCS-based channel estimation scheme [15], where the main lobe is first searched and then the largest entry within the main lobe is further searched by the exhaustive search method, the computational complexity is  For the ECS-based [11] and OCS-based [17] channel estimation schemes, the largest entry is directly searched by the exhaustive search method, the computational complexity is Since K > N A and K > M A , the computational complexity of the IA-based scheme, SZO-based scheme and CZO-based scheme is much lower than that of the ECS-based, OCSbased and DCS-based channel estimation scheme. Moreover, for the HMC-based and JOINT-based schemes, since they use hierarchical codebook for beam training instead of computing multiplication in (15) to obtain R v u , the computational complexity is almost zero.
2) Estimation Error: For four schemes including the IAbased scheme, ECS-based scheme, DCS-based scheme and OCS-based scheme, the estimation error of the AoA and AoD is larger than 2/K without noise. However, as shown in (76) and (84), the estimation error of the AoA and AoD is 2/K without noise for both SZO-based and CZO-based scheme. For the HMC-based scheme and JOINT-based scheme, the precision of channel estimation relies on the number of the layers in the hierarchical codebook. In fact, the number of time slots for training is proportional to the layers in the codebook. However, for the IA-based scheme, ECS-based scheme, DCSbased scheme, OCS-based scheme, SZO-based scheme and CZO-based scheme, K can be set independent of the number of time slots for training.
As shown in the last sub-figure of Fig. 1 for the CZObased scheme, the width of main lobe of r M and r N is 4/T 1 and 4/T 3 , respectively, which are larger than that of the IAbased scheme and SZO-based scheme. Therefore, r N and r M in the CZO-based scheme are flatter than those of the IAbased scheme and SZO-based scheme, meaning that the CZObased scheme has the largest estimation error among the three schemes with the same noise power.

3) Total Time Slots for Channel Training:
For the SZObased scheme where the beam training based on codebook is used, the number of time slots consumed by the beam training is U (5 log 2 M A +2 log 2 (N A /M A )). Additionally, we need 2U time slots to transmit pilot sequences. Therefore, the number of total times slots for channel training is U (5 log 2 M A + 2 log 2 (N A /M A )+2) for the SZO-based scheme. However, for the IA-based scheme, CZO-based scheme, ECS-based scheme, DCS-based scheme, OCS-based scheme, HMC-based scheme and JOINT-based scheme, the total time slots are flexible, depending on T 1 and T 2 .

IV. SIMULATION RESULTS
Now we evaluate the performance of three proposed schemes. We consider uplink transmission of a multi-user mmWave massive MIMO system. The BS serving U = 4 users has N A = 64 antennas and N R = 4 RF chains, while each user has M A = 16 antennas and M R = 1 RF chain. Suppose we use 6 bit and 4 bit digital phase shifters at the BS and users, respectively. The number of resolvable paths in mmWave channel is set to be L u = 3, while g u,1 ∼ CN (0, 1) and g u,i ∼ CN (0, 0.01) for i = 2, 3. We use K = 1024 over-sampling steering vectors. For SZO-based scheme, the number of total time slots for channel training is fixed to be U (5 log 2 M A + 2 log 2 (N A /M A ) + 2) = 144. For the IAbased scheme, SZO-based scheme, ECS-based scheme, HMCbased scheme and JOINT-based scheme, the total time slots for channel training are flexible.
As shown in Fig. 5, we compare the channel estimation performance in terms of normalized mean square error (NMSE) for the proposed three schemes, the DCS-based scheme [15] and the OCS-based scheme [17] with different SNR. The NMSE is defined (θ u,1 − θ u,1 ) 2 + (φ u,1 − ϕ u,1 ) 2 , which reflects the estimation accuracy of AoA and AoD. We set T 1 = 4 and T 2 = 4. Then T 3 = T 2 N R = 16. In order to make   based, HMC-based and JOINT-based schemes, respectively. The reason for the performance of the DCS-based and OCSbased schemes is that both the DCS-based and OCS-based schemes employ random precoding and random combining matrix, which are not optimal. Although both the DCSbased and OCS-based schemes outperform the CZO-based scheme in low SNR region, at high SNR region such as SNR of 20dB, the CZO-based scheme has 86.8% and 34.5% improvement compared with the DCS-based and OCS-based schemes, respectively. The reason is that the width of the main lobe in CZO-based scheme is larger than DCS-based and OCS-based scheme. Besides, both the SZO-based scheme and CZO-based scheme can approach the ideal performance as SNR increases, which shows the superiority of SZO-based scheme and CZO-based scheme in their small estimation error. The reason for the unsatisfactory performance of the ECSbased scheme is that the ECS-based scheme does not consider the power leakage due to the limited beamspace resolution. Therefore the ideal sparse property of beamspace channel is impaired. The reason for the unsatisfactory performance of the HMC-based and JOINT-based schemes is the number of codewords in the last layer of hierarchical codebook is only 8 and much smaller than K = 1024, where the total time slots for training are set to be the same with other schemes for fair comparisons. In addition, compared with the CZO-based scheme, the IA-based scheme has smaller estimation error in the low SNR region. The IA-based scheme can achieve better NMSE performance than the CZO-based scheme when SNR < 20 dB, which means we can use the IA-based scheme for better NMSE performance when SNR < 20dB. As shown in Fig. 6, we compare the sum-rate for the proposed three schemes, the DCS-based, OCS-based, ECSbased, HMC-based and JOINT-based schemes. It is seen that the proposed three schemes achieve better performance than the other schemes, especially in low SNR region. At SNR of 10dB, the performance gap between the proposed three schemes and the ideal case is around 1dB.
As shown in Fig. 7, we compare the channel estimation performance in terms of NMSE for different schemes with different number of total time slots for channel training, which is T = U T 1 T 2 . Since the number of total time slots for channel training is fixed to be 144 for the SZO-based scheme, SZObased scheme is not included for comparison. We fix SNR to be 15dB. Given T , we set T 1 = T 2 = T /U . It is seen that when T is small, the CZO-based scheme performs the best; when T is large, the IA-based scheme performs the best. When T = 16, the CZO-based scheme has 78.2%, 77.8%, 83.5%, 76.0% and 76.1% performance improvement compared with the DCS-based, OCS-based, ECS-based, HMC-based and JOINT-based schemes, respectively. When T = 100, the IA-based scheme has 76.3%, 11.2%, 81.4%, 91.4% and 91.5% improvement compared with the DCS-based, OCSbased, ECS-based, HMC-based and JOINT-based schemes, respectively. In the ideal case that both W H W = γ N I NA and F u F H u = γ M I MA can be achieved, the number of total time slots for channel training is U M A N A = 4096. It is seen from Fig. 7 that with substantially reduced training overhead, i.e., T ≪ 4096, satisfactory performance can almost be achieved, e.g., 0.0027 of NMSE with T = 100 for the IA-based scheme. The number of total time slots for channel training is flexible and fixed for the IA-based scheme and the SZO-based scheme, respectively, which indicates the flexibility of the IA-based scheme. To achieve the NMSE of 0.0027, T = 100 is enough for the IA-based scheme, while T = 144 is fixed for the SZObased scheme.
As shown in Fig. 8, we compare the sum-rate for different number of total time slots for channel training. We fix SNR to be 15dB. We observe that the CZO-based scheme can achieve better performance than the DCS-based, OCS-based, ECS-based, HMC-based and JOINT-based schemes when T is small. When T is close to 100, the sum-rate of the IA-based scheme and CZO-based scheme is almost invariant, indicating that T = 100 is enough to achieve the maximal sum-rate.
As shown in Fig. 9, we verify the convergence of Algorithm 1. Suppose we use 5, 6 and 7 bit digital phase shifters at the BS, respectively. It is seen that ǫ in (26) decreases rapidly as the number of iterations grows, and ǫ decreases faster with higher resolution of digital phase shifters. Using 6 bit digital phase shifters at the BS, ǫ is smaller than 10 −3 when the number of iterations is larger than 25, which means 25 iterations is enough for δ = 10 −3 .

V. CONCLUSIONS
This paper has investigated beamspace channel estimation for multi-user mmWave massive MIMO system. A framework of beamspace channel estimation has been proposed. Then based on the this framework, three channel estimation schemes have been proposed. These schemes together with the existing channel estimation schemes have been compared in terms of computational complexity, estimation error and total time slots for channel training. Simulation results have shown that the proposed schemes outperform the existing schemes and can approach the performance of the ideal case with substantially reduced training overhead. Future work will focus on the design of hybrid precoding and hybrid combining for multiuser data transmission regarding the energy efficiency and the theoretical proof of the convergence of iterative hybrid precoding and combining design.