Tensor-Based Parameter Estimation of Double Directional Massive MIMO Channel with Dual-Polarized Antennas

TENSOR-B ASED P ARAMETER ESTIMA TION OF DOUBLE DIRECTIONAL MASSIVE MIMO CHANNEL WITH DU AL-POLARIZED ANTENNAS Cheng Qian † Xiao Fu ? Nicholas D. Sidir opoulos † Y e Y ang ‡ † Dept. of Electrical and Computer Engineering, Uni versity of V irginia, Charlottesville, V A 22904 ? School of Electrical Engineering and Computer Science, Oregon State Uni v ersity , Corvallis, OR 97331 ‡ Physical Layer & RRM IC Algorithm Dept., WN Huawei Co., Ltd, Shanghai, China ABSTRA CT The 3GPP suggests to combine dual polarized (DP) antenna arrays with the double directional (DD) channel model for downlink channel estimation. This combination strikes a good balance between high-capacity communications and parsimonious channel modeling, and also brings limited feed- back schemes for downlink channel estimation within reach. Howe v er , most e xisting channel estimation work under the DD model has not considered DP arrays, perhaps because of the complex array manifold and the resulting difficulty in algorithm design. In this paper , we first reveal that the DD channel with DP arrays at the transmitter and receiver can be naturally modeled as a lo w-rank four-way tensor , and thus the parameters can be effecti vely estimated via tensor decompo- sition algorithms. T o reduce computational complexity , we show that the problem can be recast as a four-snapshot three- dimensional harmonic retrie v al problem, which can be solv ed using computationally ef ficient subspace methods. On the theory side, we sho w that the DD channel with DP arrays is identifiable under very mild conditions, lev eraging identi- fiability of low-rank tensors. Numerical simulations are em- ployed to sho wcase the effecti v eness of our methods. Index T erms — Channel estimation, massi ve MIMO, dual-polarized array , tensor factorization, harmonic retriev al. 1. INTRODUCTION The dual-polarized (DP) antenna array has many appealing features and is thus considered a key technique for next gen- eration communications and massiv e MIMO [1–3]. For ex- ample, Foschini and Gans [4] sho wed that the capacity for systems with DP antennas at the transmitter can be increased up to 50% compared to systems without polarization. Besides the increased capacity , DP antennas ha ve other key advan- tages such as small size, easy installation, good interference mitigation performance, high link reliability , and high ability of interference filtering, just to name a few [1–4]. Author e-mails: alextoqc@gmail.com, xiao.fu@oregonstate.edu, nikos@virginia.edu, yangye@hua wei.com. In the recent releases of technical specifications suggested by 3GPP , the DP array and the double directional (DD) chan- nel model are considered ke y techniques [2]. The DD channel model is parsimonious for multipath channels with a small number of dominant paths, and parsimony is really essential for designing limited feedback schemes for downlink chan- nel estimation in massi v e MIMO [1–3, 6]. Specifically , 3GPP suggests that the mobile users estimate the DD channel pa- rameters such as directions-of-arriv al (DOAs), directions-of- departure (DODs), the path loss associated with each path and the polarization angles of the DP array , and then feed back these parameters to the base station (BS). This strat- egy is rather economical, as it is e xpected that the number of dominant paths will be small to moderate in practical de- ployments. On the other hand, to the best of our knowledge, there is very limited work related to the DD-DP parameter estimation problem. Most of the existing channel estimation algorithms such as [5–8] do not take polarization into consid- eration, and thus cannot be applied to this particular system. The early algorithm proposed in [9] in the context of array processing can only handle a small number of paths, since the maximum number of identifiable paths in [9] is restricted by the size of the receiv e array . In this work, we focus on the parameter estimation prob- lem under the DD channel model with DP arrays. Specifi- cally , we first sho w that the DD channel with DP arrays at the transmitter and recei ver can be naturally modeled as a low- rank four -way tensor . Leveraging this structure, we recast the associated parameter estimation problem as a Parallel Factor Analysis (P ARAF A C) decomposition problem [10] and han- dle it using effecti ve tensor decomposition algorithms. T o reduce computational complexity , we also formulate chan- nel estimation as a three-dimensional (3-D) harmonic re- triev al problem, which can be solv ed by a computationally efficient subspace method, namely , the impr oved multidimen- sional folding (IMDF) method [11]. On the theory side, we show that the channel and polarization parameters are identi- fiable under very mild and practical conditions – e ven when the number of paths largely exceeds the number of recei ve antennas, a practically important case that classic DP channel estimation algorithms as in [9] cannot cope with. Simulations are provided to sho wcase the effecti veness of the proposed methods. 2. SIGNAL MODEL W e consider a massiv e MIMO system, where there is one BS equipped with an M x × M y DP uniform rectangular ar- ray (URA) and one mobile station (MS) with an M r -element DP uniform linear array (ULA), which is a practical setting that is of interest to industry [2]. Throughout the paper, we consider DP array elements consisting of a pair of crossed dipoles. In the literature, this type of DP array is also known as a “cross-polarized” array [9]. The number of transmit an- tennas is M t = M x M y . The signal received by the user is giv en by x ( t ) = Hs ( t ) + n ( t ) , t = 1 , · · · , N (1) where s ( t ) ∈ C 2 M t × 1 is the transmitted signal, n ( t ) is zero- mean i.i.d.circularly symmetric complex Gaussian noise. By properly arrange elements, the downlink channel matrix can be represented as the following four -block matrix: H =  H (V r , V t ) H (V r , H t ) H (H r , V t ) H (H r , H t )  ∈ C 2 M r × 2 M t . (2) where H (V r , V t ) ∈ C M r × M t is a channel matrix between all the V -polarized transmit antennas and V -polarized receiv e an- tennas, and H (V r , H t ) ∈ C M r × M t is a channel matrix between all the H-polarized transmit antennas and V -polarized recei ve antennas; like wise for the other two blocks in (2). For notational simplicity , let p ∈ { V r , H r } and q ∈ { V t , H t } . Then, according to [1], the (p , q) subchannel ma- trix is modeled as H (p , q) = V r diag  β (p , q)  V H t (3) where ( · ) H is the conjugate transpose, V r = [ v r ( θ 1 ) · · · v r ( θ K )] , V t = [ v t ( ϑ 1 , ϕ 1 ) · · · v t ( ϑ K , ϕ K )] , and β (p , q) = [ β (p , q) 1 · · · β (p , q) K ] T stands for path-losses with ( · ) T denoting the transpose. Note that { θ k } are DOAs, { ϑ k } and { ϕ k } are azimuth and elev ation DODs, respec- tiv ely . Throughout of this paper, we assume that the receiv e and transmit antennas have the same half-wa velength inter- element spacing. Then we ha ve [ v r ( θ k )] m = e j π ( m − 1) sin( θ k ) and v t ( ϑ k , ϕ k ) = v y ,k ⊗ v x,k , where [ x ] i denotes the i th element of x , [ v x,k ] l x = e j π ( l x − 1) sin( ϕ k ) cos( ϑ k ) , l x = 0 , · · · , M x − 1 and [ v y ,k ] l y = e j π ( l y − 1) sin( ϕ k ) sin( ϑ k ) , l y = 0 , · · · , M y − 1 . Now the channel matrix in (2) can be re written as H = " V r diag  β (V r , V t )  V H t V r diag  β (V r , H t )  V H t V r diag  β (H r , V t )  V H t V r diag  β (V t , H t )  V H t # . (4) In this model, to determine the channel H , we only need to estimate K DOAs, K azimuth angles, K elev ation angles and 4 K complex path-losses. Compared to the size of the chan- nel, which is 4 M r M t , such parameterization is rather eco- nomical and is suitable for massi ve MIMO do wnlink channel estimation and limited feedback where both M t and M r (es- pecially M t ) can be very lar ge. 2.1. Challenges Although we have e xplicitly written down the channel model in (1), ho w to effecti vely estimate the parameters of interest is still unclear . Specifically , assume that H can be estimated at the receiver by matched filtering, i.e., H = XS H under a pre-selected row-orthogonal pilot sequence S , where X = [ x (1) , . . . , x ( N )] and S = [ s (1) , . . . , s ( N )] . Estimating the DO A, DOD and path-loss parameters is still very challenging. One popular type of technique to estimate parameters of the (non-DP) DD channel is described in [5–7], where the DOA and DOD domains are descretized to fine angle grids using two overcomplete angle dictionaries (codebooks), denoted by D t and D r . Then, we have H ≈ ( I 2 ⊗ D r ) G ( I 2 ⊗ D t ) H , where G is a sparse matrix that selects out the columns asso- ciated with the acti ve DODs and DO As from the dictionaries. This way , the parameter estimation problem becomes a sparse recov ery problem that can be handled by formulations such as min g k h − ( I 4 ⊗ D ∗ t ⊗ D r ) g k 2 2 + λ k g k 1 , where h = vec( H ) with vec ( · ) being the vectorization operator and g = vec( G ) ; and other sparse optimization algorithms such as orthogonal matching pursuit . The difficulty is that to ensure good spatial resolution, both D t ∈ C M t × D t and D r ∈ C M r × D r are very “fat” matrices, where D t and D r denotes the number of angle grids after quantization. Consequently , ( I 4 ⊗ D ∗ t ⊗ D r ) is 4 M t M r × 4 D t D r . If one quantizes the DO A and DOD space (ranging from − 90 ◦ to 90 ◦ ) using a resolution of one de- gree, then 4 D t D r = 131 , 044 – which poses an extremely hard sparse optimization problem. Many compromises, such as coarse quantization and hierarchical or hybrid algorithms [5, 6], hav e been employed to circumvent this issue in the lit- erature. Ho wev er , when two-dimensional antennas are de- ployed in both receiver and transmitter , the dictionary size can reach 4( D t D r ) 2 , which is hopeless. 3. PROPOSED APPR O A CH 3.1. T ensor -Based Method and Identifiability Our proposed approach starts by noticing that H is in fact a four-w ay tensor of rank (at most) K ; to see this, vectorize the four blocks in H and then stack them in a tall matrix, such that we hav e ˘ H =  V ∗ y  V ∗ x  V r  B T (5) where ( · ) ∗ denotes conjugation,  is the Khatri-Rao prod- uct, V x = [ v x, 1 · · · v x,K ] , V y = [ v y , 1 · · · v y ,K ] and B = [ β (V r , V t ) β (V r , H t ) β (H r , V t ) β (H r , H t ) ] T ∈ C 4 × K . Note that (5) is the definition of a four-way tensor of rank ≤ K in matrix form [10]. By noticing the tensor structure of H , v arious tensor de- composition algorithms such as those in [12, 13] can be di- rectly applied to estimate V x , V y , V r and B via solving the following: min V r , V x , V y , B   ˇ H −  V ∗ y  V ∗ x  V r  B T   2 F (6) where k · k F is the Frobenius norm. Note that a salient feature of tensors is that the factors are uniquely identifiable under mild conditions, as we will explain shortly . Once V x , V y , V r and B are estimated, the parameters { ˆ θ k , ˆ ϑ k , ˆ ϕ k , β (p , q) } can be computed in closed-form. Since v x,k , v y ,k and v r,k are V andermonde vectors, we may use ˆ θ k = sin − 1  1 π ∠ ( ˆ v H r,k ˆ v r,k )  (7) ˆ ϕ k = sin − 1  1 π q  ∠ ( ˆ v H x,k ˆ v x,k )  2 +  ∠ ( ˆ v H y ,k ˆ v y ,k )  2  (8) ˆ ϑ k = tan − 1  ∠ ( ˆ v H y,k ˆ v y,k ) / ∠ ( ˆ v H x,k ˆ v x,k )  (9) where ∠ ( · ) takes the phase of its argument, x and x are the vectors consisting of the first and last ( M − 1) entries of x with length M , respecti vely . An y other single-tone fre- quency estimation algorithm, e.g., [14, 15] or ML-based (pe- riodogram) methods can also be used, for better accuracy . W e should mention that by solving (6) using any of the existing tensor decomposition algorithms, we already ha ve an initial estimate of B , i.e., the path-losses. Ho we ver , since there is an intrinsic scaling ambiguity of tensor decomposi- tion, such an initial estimate may not be useful. Neverthe- less, this issue is easy to fix. Note that the array manifolds ˆ A r , ˆ A x , ˆ A y without scaling ambiguity can be constructed from { ˆ θ k , ˆ ϕ k , ˆ ϑ k } K k =1 . Then, the estimate of B without scal- ing ambiguity can be computed from the following LS prob- lem: ˆ B ← arg min B    ˘ H − ( ˆ V ∗ y  ˆ V ∗ x  ˆ V r ) B T    2 F . (10) In terms of theoretical guarantees of identifiability , we hav e the following theorem: Theorem 1 The pr oposed appr oach can uniquely identify the parameters of inter est under the DD channel model with DP arrays pr ovided that min ( M r , K ) + min( M x , K ) + min( M y , K ) + min (4 , K ) ≥ 2 K + 3 . One can easily check that { V x , V y , V r , B } meet the k - rank condition [17] provided that all the DO A, DOD and path- loss are not the same, which is a mild condition considering the random nature of multi-path. Thus, Theorem 1 essentially follows from [16]. Much better results can also be claimed, albeit in the almost surely sense – see [10]. 3.2. IMDF and Identifiability The ‘naive’ tensor-based method ignores the V andermonde structure of some of the array manifold vectors in its first step, only to impose it later . This is suboptimal. Theorem 1 in particular is a general bound that neglects the V andermonde structure in v x , v y and v r . If we take this structure into account, a better uniqueness condition can be obtained. T o this end, we rearrange the elements of ˘ H such that the re- sulting tensor is with dimension M r × M x × M y × 4 , i.e., P K k =1 v r,k ◦ v ∗ x,k ◦ v ∗ y ,k ◦ b k , where ◦ denotes the outer prod- uct and b k is the k th column of B . The above can be viewed as a multi-snapshot 3-D harmonic retriev al problem, where the number of snapshots is four, and each snapshot is written as H H H (p , q) = K X k =1 β (p , q) k v r,k ◦ v ∗ x,k ◦ v ∗ y ,k , ∀ p , q . (11) Theorem 2 The parameters  θ k , ϕ k , ϑ k , β (p , q) k  ar e all uniquely identifiable by the IMDF based pr ocedur e pr ovided that K ≤ arg max F,P r ,P x ,P y F s . t . max  ( P r − 1) P x P y , P r ( P x − 1) P y , P r P x ( P y − 1)  ≥ F 8 Q r Q x Q y ≥ F (12) wher e P r + Q r = M r + 1 , P x + Q x = M x + 1 , P y + Q y = M y + 1 . This follows by inv oking the identifiability result for the IMDF algorithm for multi-dimensional harmonic retriev al [11], which is far stronger compared to that in Theorem 1. For example, when M x = 4 , M y = 8 , and M r = 2 , the iden- tifiability of Theorem 1 is K = 7 , while the identifiability of Theorem 2 is K = 32 . Furthermore, ev en when the MS only has a single dual-polarized antenna, it can be sho wn using the IMDF based approach that the number of identifiable paths is upper bounded by K < 0 . 8187 M t . In Algorithm 1, we show the detailed procedures for esti- mating multipath parameters using IMDF . 4. NUMERICAL RESUL TS Consider a MIMO system with an 4 × 8 DP URA at the BS and a 2 -element DP ULA at the MS. This particular case is of considerable practical interest in 3GPP as a can- didate for implementation [2]. In the simulation, we assume that the multipath propagation gains are Rician distrib uted, and all the multipath parameters are randomly (uniformly) drawn. In the simulation, we assume that the multipath Algorithm 1 IMDF for DD-DP Parameter Estimation 1: Compute the least squares (LS) estimate of H , i.e., ˆ H LS , and form ˆ H H H (p , q) via (11). 2: Use Theorem 2 to pre-calculate { P x , P y , P r , Q x , Q y , Q r } , such that each ˜ H H H (p , q) can be reshaped into a P x P y P r × Q x Q y Q r matrix which is denoted as ˜ H (p , q) . 3: Perform forward-backward smoothing on the conjugate of ˜ H (p , q) to obtain ` H (p , q) , and then ∀ p ∈ { (V r , H r ) } and q ∈ { (V t , H t ) } , stack { ˜ H (p , q) , ` H (p , q) } into a P x P y P r × 8 Q x Q y Q r matrix, denoted by ˇ H . 4: Perform 3-D IMDF to ˘ H and obtain the estimates of { θ k , ϕ k , ϑ k } . 5: Use { ˆ θ k , ˆ ϕ k , ˆ ϑ k } to construct ˆ V r , ˆ V x , ˆ V y , and then esti- mate the path-loss matrix B via (10) propagation gains are Rician distributed, and all the multi- path parameters are randomly (uniformly) drawn. The BS cov ers [0 ◦ , 90 ◦ ] elev ation angular range and ( − 45 ◦ , 45 ◦ ) az- imuth angular range, while the MS only covers [ − 60 ◦ , 60 ◦ ] azimuth angular range since the elev ation angle is zero for ULA, i.e., θ k ∼ U ( − π / 3 , π / 3) , ϕ k ∼ U (0 , π/ 2) , ϑ k ∼ U ( − π / 3 , π / 3) . The non-parametric linear LS channel esti- mate is also plotted as a performance benchmark. All the results are a veraged o ver 500 Monte-Carlo trials using a com- puter with 3.2 GHz Intel Core i5-4460 and 4 GB RAM. The normalized MSE (NMSE) of channel estimates is computed from 1 500 P 500 i =1 k ˆ H i − H k 2 F / k H k 2 F where ˆ H i denotes the es- timate from the i th Monte-Carlo trial. The number of multipath randomly varies from 1 to 6. Since the channel exhibits sparse property , we include a com- pressiv e sensing (CS) based technique [7] for comparison, where each angle is quantized with 7 bits, so the resulting dictionary is with size 4 M r M t × 2 23 , which howe ver is in- feasible in a con ventional desktop. T o make this algorithm work in a fast fashion, after obtaining the LS channel esti- mate, we reshape each sub-block of the channel estimate as an M r × M x × M y tensor and average them. Then we im- plement 3-D FFT with 128 points to estimate { θ , ϑ, ϕ } , fol- lowing the so-called peak-picking technique. Finally , we up- date the path-loss matrix B via (10). W e test the performance of all the competitors under known and unknown number of multipath. For the latter, we set K = 6 to all the algorithms. Moreov er , orthogonal pilots are employed. It is observ ed from Fig. 1 that P ARAF A C outperforms the IMDF , LS and CS algorithms in both cases. Compared to Fig. 1(a), P ARAF A C, IMDF and CS suffer slight performance loss in Fig. 1(b), where the exact number of multipath is unknown. When SNR > 14 dB, we see that the NMSE of CS is even worse than the LS method. This is mainly because as SNR increases, the performance of CS is limited by the resolution ability of dictionary . (a) known K (b) unknown K Fig. 1 . NMSE of versus SNR. 5. CONCLUSION W e considered the parameter estimation problem for the DD channel model with DP arrays – which is a setup that is of particular interest to standard organizations and industry . W e proposed a tensor-based method to handle this challenging problem, which guarantees identifiability of the parameters of interest under mild and practical conditions. W e also pro- posed a reduced-complexity algorithm that is based on 3D harmonic retriev al to handle the same problem, with slight pa- rameter accuracy loss but much faster runtime performance. Numerical simulations support our analysis and sho w that the proposed procedures, esp. the IMDF-based one, are very ef- fectiv e and promising for actual implementation. 6. REFERENCES [1] A. Kammoun, H. Khanfir , Z. Altman, M. Debbah and M. Kamoun, “Preliminary results on 3-D channel mod- eling: From theory to standardization, ” IEEE Journal on Selected Ar eas in Communications , vol. 32, no. 6, pp. 1219-1229, June 2014. [2] Z. Bai, “Evolv ed univ ersal terrestrial radio access (E- UTRA); physical layer procedures, ” 3GPP , Sophia An- tipolis, T echnical Specification , 36.213 v . 11.4.0, 2013. [3] D. Zhu, J. Choi and R. W . Heath, “T wo-dimensional AoD and AoA acquisition for wideband mmW av e sys- tems with dual-polarized MIMO, ” IEEE T rans. W ir eless Comm. , , vol. 16, no. 12, pp. 7890-7905, Dec. 2017. [4] G. J. Foschini and M. J. Gans, “On limits of wireless com- munications in a fading en vironment when using multiple antennas, ” W ireless P ersonal Communications , v ol. 6, no. 3, pp. 311-335, 1998. [5] A. Alkhateeb, G. Leus and R. W . Heath,“Limited feed- back hybrid precoding for multi-user millimeter wa ve systems, ” IEEE T rans. W ir eless Comm. , vol. 14, no. 11, pp. 6481-6494, 2015. [6] P . N. Alevizos, X. Fu, N. Sidiropoulos, Y . Y ang and A. Bletsas, “Non-uniform directional dictionary-based lim- ited feedback for massiv e MIMO systems, ” Pr oc. of 15th Inter . Symp. Modeling and Optimization in Mobile, Ad Hoc, and W ireless Networks (W iOpt), Paris, pp. 1-8, 2017. [7] W .U. Bajwa, J. Haupt, A.M. Sayeed, and R. No wak, “Compressed channel sensing: A new approach to es- timating sparse multipath channels, ” Proc. of the IEEE , vol. 98, no. 6, pp. 1058-1076, 2010. [8] H. Lin, F . Gao, S. Jin and G. Y . Li, “ A ne w view of multi- user hybrid massiv e MIMO: Non-orthogonal angle divi- sion multiple access, ” IEEE J. Selected Areas in Comm., vol. 35, no. 10, pp. 2268-2280, Oct. 2017. [9] J. Li and R. T . Compton, “T wo-dimensional angle and po- larization estimation using the ESPRIT algorithm, ” IEEE T rans. Antennas Pr opagat. , vol. 40, pp. 550-555, 1992. [10] N.D. Sidiropoulos, L. De Lathauwer , X. Fu, K. Huang, E.E. Papalexakis, and C. Faloutsos, “T ensor decomposi- tion for signal processing and machine learning”, IEEE T rans. Signal Pr ocess. , vol. 65, no. 13, pp. 3551-3582, 2017 [11] J. Liu and X. Liu, “ An eigen vector-based approach for multidimensional frequency estimation with improved identifiability , ” IEEE T rans. Signal Process. , vol. 54, no. 12, pp. 4543-4556, 2006. [12] L. D. Lathauwer , B. D. Moor , and J. V andewalle, “Com- putation of the canonical decomposition by means of a simultaneous generalized schur decomposition, ” SIAM J. Matrix Anal. Appl. , vol. 26, no. 2, pp. 295-327, 2004. [13] K. Huang, N. D. Sidiropoulos, and A. P . Liav as, “ A fle x- ible and efficient algorithmic framework for constrained matrix and tensor factorization, ” IEEE T rans. Signal Pr o- cess. , vol. 64,no. 19, pp. 5052-5065, 2016. [14] C. Qian, L. Huang, H. C. So, N. D. Sidiropoulos, and J. Xie, ‘Unitary PUMA algorithm for estimating the fre- quency of a complex sinusoid, ” IEEE T rans. Signal Pr o- cess. , vol. 63, no. 20, pp. 5358-5368, 2015. [15] D. C. Rife and R. R. Boorstyn, “Single tone parameter estimation from discrete-time observ ations, ” IEEE T rans. Inf. Theory , v ol. IT -20, no. 5, pp. 591-598, 1974. [16] N. Sidiropoulos and R. Bro, “On the uniqueness of mul- tilinear decomposition of N -way arrays, ” J. Chemomet- rics , vol. 14, no. 3, pp. 229-239, 2000. [17] N. D. Sidiropoulos, and X. Liu, “Identifiability results for blind beamforming in Incoherent multipath with small delay spread, ” IEEE T rans. on Signal Pr ocess. , vol. 49, no. 1, pp. 228-236, 2001.