Uplink Sparse Channel Estimation for Hybrid Millimeter Wave Massive MIMO Systems by UTAMP-SBL

The compressive sensing (CS)-based sparse channel estimator is recognized as the most effective solution to the excessive pilot overhead in massive MIMO systems. However, due to the complex signal processing in the wireless communication systems, the measurement matrix in the CS-based channel estimation is sometimes “unfriendly” to the channel recovery. To overcome this problem, in this paper, the state-of-the-art sparse Bayesian learning using approximate message passing with unitary transformation (UTAMP-SBL), which is robust to various measurement matrices, is leveraged to address the multi-user uplink channel estimation for hybrid architecture millimeter wave massive MIMO systems. Specifically, the sparsity of channels in the angular domain is exploited to reduce the pilot overhead. Simulation results demonstrate that the UTAMP-SBL is able to achieve effective performance improvement than other competitors with low pilot overhead.


Introduction
Massive MIMO is a key technology for the fifth-generation (5G) communication systems [1]. Thanks to the shorter wavelength of the millimeter wave (mmWave) signal, the numerous antennas are packed into a compact-size array, which facilitates the commercial deployment of massive MIMO systems [2], but under the consideration of the high power consumption of Analog-to-Digital Converters (ADCs), the hybrid beamforming architecture, that divides the precoder and the combiner into the analog and digital domains, is developed to solve it and regarded as an effective alternative [3,4]. In this hybrid architecture system, the precoding is crucial and dependent on the acquisition of accurate channel state information (CSI). However, the number of pilots used for channel estimation will increase linearly with that of antennas and users in the wireless systems. As a result, the channel estimation with fewer pilots is a significant challenge in the hybrid mmWave massive MIMO systems [3].
The compressive sensing (CS)-based channel estimator is capable of exploiting the sparsity of channels to reduce the pilot overhead greatly and a lot of studies on this topic have been done. Some studies focus on the sparsity in the delay domain. In [5], the orthogonal matching pursuit (OMP) algorithm is utilized to estimate the wideband mmWave delay-domain sparse channels. In [6], an adaptive structured subspace pursuit algorithm at the user is proposed to estimate the delay-domain MIMO channels with the spatio-temporal common sparsity. The authors of [7] propose a sparse channel recovery algorithm based on the vector approximate message passing (VAMP) for the massive MIMO delay-domain channels. Some works concentrate on the sparsity in the angular domain. For example, a distributed sparsity adaptive matching pursuit (DSAMP) algorithm is proposed in [8], whereby the spatially common sparsity is exploited. In [9], a novel sparse Bayesian learning (SBL) approach is utilized to recover the angular-domain block-sparse channels. In [10], a structured turbo-CS algorithm is proposed to estimate the angulardomain massive MIMO channels which are modeled by a Markov chain prior. Furthermore, the work [11] even jointly exploits the 3-D clustered structure of channels in the angulardelay domain and proposes an approximate message passing (AMP)-based estimation algorithm. However, all of above angular-domain sparse channel estimation schemes are based on the assumption that the angles exactly lie on the grids, or ignore the power leakage caused by grid mismatch. Recently, there are many off-grid channel estimation techniques. In [12], the authors utilize the low-rank structure along with the sparsity in angular domain to improve the channel estimation performance, and the off-grid angles can be recovered with their algorithm successfully. In [13], the angles are treated as random parameters and the grid-less quantized variational Bayesian channel estimation algorithm is proposed for antenna array systems with low resolution ADCs. Besides, the multi-dimensional variational line spectral estimation algorithm proposed in [14] can be effectively applied for multi-dimensional off-grid angle estimation. In addition, a novel super-resolution downlink channel estimation approach developed from the SBL is provided in [15], where the sampled angular grid points are treated as the underlying parameters.
The performance of these CS-based algorithms is affected by the measurement matrix which is usually associated with the signal processing operations in the system. If these operations are not carefully designed, the measurement matrix perhaps becomes detrimental to channel recovery. Recently, a CS algorithm, termed the SBL using approximate message passing with unitary transformation (UTAMP-SBL), is proposed by Luo et al. and it outperforms the state-of-the-art AMP-based SBL algorithms, the Gaussian generalized AMP-based SBL (GGAMP-SBL), in terms of robustness, speed and recovery accuracy for difficult measurement matrices [16]. In this paper, we apply UTAMP-SBL to the sparse channel estimation to improve the performance. Specifically, the multi-user uplink channel estimation for hybrid architecture mmWave massive MIMO systems is studied as an example. And the angular-domain sparsity of the mmWave channels is fully exploited in our estimation. Additionally, this algorithm can be extended to other sparse channel estimation for the performance improvement. Simulation results verify that the UTAMP-SBL outperforms the other competitors.
The remainder of this paper is organized as follows. In Section 2, the hybrid millimeter wave massive MIMO system model is introduced. In Section 3, the uplink sparse channel estimation with UTAMP-SBL is described. In Section 4, the Cramér-Rao bound (CRB) is provided as a performance benchmark. Simulation results are provided and the performance is discussed in Section 5. Conclusions are given in Section 6.
Notations: H and ⊗ denote the conjugate transpose operation and the Kronecker product. and represent the componentwise vector multiplication and the componentwise vector division. vec(A) denotes vectorizing the matrix A as a vector. I, 1 and 0 denote the identity matrix, the all-one column vector and the all-zero column vector, respectively. diag(a) returns a diagonal matrix with the elements of vector a on its main diagonal. N C (x; µ x , Σ x ) represents the Gaussian distribution of the complex vector x with mean µ x and covariance matrix Σ x . Ga(γ; , η) denotes a Gamma distribution with shape parameter and rate parameter η.
] denotes the expectation of the function f (x) with respect to probability density b(x). Finally, ∝ denotes equality up to a constant scale factor.

Millimeter Wave MIMO Channel Model
We consider a mmWave MIMO-OFDM system where the base station (BS) is equipped with N BS antennas and serves K multi-antenna user equipments (UEs), and each UE has N UE antennas. The frequency-domain channel between the BS and the kth UE at the pth (p = 1, 2, · · · , P) subcarrier can be modeled as [17] where L is the number of physical paths, α l,k is the complex path gain, f s is the sampling rate, τ l,k is the path delay, a BS (θ l,k ) and a UE (φ l,k ) are the steering vectors at the BS and the UE, respectively, θ l,k and φ l,k are the azimuth angles of departure and arrival (AoD/AoA) uniformly distributed in [0, 2π). For each UE in the system, we generate one line of sight (LoS) path and L − 1 non-LoS (NLoS) paths. The path gains follow the complex Gaussian distribution and the ratio of the power of LOS path to that of NLoS path is 10 dB [18]. In the OFDM system, the cyclic prefix (CP) is introduced to mitigate the inter-symbol interference and we have τ max f s ≤ N cp , where τ max is the maximum path delay and N cp is the length of CP. So we generate τ 1,k f s = 1 and τ l,k f s ∈ {2, 3, · · · , N cp }, l = 1. For the steering vectors, assuming the uniform linear arrays are adopted at the BS and the UE, they can be shown as where λ is the wavelength and d = λ/2 is the antenna spacing.

Uplink Pilot Transmission
We consider the typical hybrid analog digital precoding and combining architecture in the mmWave massive MIMO systems, where the BS and each UE are assumed to have N RF BS RF chains and N RF UE RF chains, respectively. We focus on the uplink pilot training. The kth UE transmits the pilot s p,t,k ∈ C N s ×1 at the pth pilot subcarrier and the tth OFDM symbol to the BS, where N s is the number of data streams. At the BS, the receiver receives pilots from all UEs in the same time-frequency resources and the received signal can be written as where F BB p,t,k ∈ C N RF UE ×N s and F RF t,k ∈ C N UE ×N RF UE are the digital precoder and the analog precoder of the kth user, Z BB p,t ∈ C N RF BS ×N s and Z RF t ∈ C N BS ×N RF BS are the digital combiner and the analog combiner of the BS and n p,t ∈ C N BS ×1 is the additive white Gaussian noise (AWGN). Note that the analog precoder and the analog combiner are frequency-flat, while the digital precoder and the digital combiner are different for each subcarrier. This is because the RF phase shifters can provide constant response over the wide frequency band when the fully connected network is used.

Sparse Channel Estimation Formulation
In the multi-user uplink channel estimation, the channels of K users, H p,k ∈ C N BS ×N UE , k = 1, 2, . . . , K, need to be estimated from the received signal y p,t . It is clear that there are KN BS N UE non-zero elements in these K channel matrices. Because there are only N s elements in y p,t , we have to arrange KN BS N UE /N s pilots for channel estimation and the pilot overhead is unacceptable. Fortunately, we benefit from the inherent sparse characteristic of millimeter wave channels and apply the CS algorithms in the channel estimation to reduce the pilot overhead.
In order to formulate the sparse channel estimation problem, we first rewrite the frequency-domain channel matrix H p,k ∈ C N BS ×N UE as where H ω p,k ∈ C G BS ×G UE is the virtual angle-domain channel matrix, A BS ∈ C N BS ×G BS and A UE ∈ C N UE ×G UE are the partial discrete Fourier transform (DFT) matrices, G BS and G UE are the numbers of grids at the BS and the UE. Because of the domination of the LOS path in the multipath channel, the virtual angle-domain mmWave channel appears the sparsity, i.e., there are only a few non-zero elements in the channel matrix H ω p,k and other components are zero. Actually, these zeros are not strictly equal to zero but are only close to zero owing to the power leakage problem [19], which is ignored in this paper for simplification and it will be disscussed in the future research. The frequency-domain channel H p,k and the transformed virtual angle-domain channel H ω p,k are shown in Figure 1, where the frequency-domain channel is generated according to (1) and the virtual angledomain channel is derived from (4). Obviously, from Figure 1b, the power of the channel is concentrated in a few virtual angles and these channel components can be estimated from the low-overhead pilot. Substituting (4) in (3), we have where Z p,t = Z RF t Z BB p,t ∈ C N BS ×N s is the combiner at the BS, x p,t,k = F RF t,k F BB p,t,k s p,t,k ∈ C N UE ×1 is the pilot after precoding at the UE andn p,t = Z H p,t n p,t ∈ C N s ×1 is the noise after combining at the BS. Note that usually the channel H ω p,k is invariant in some successive OFDM symbols. Rewrite (5) in matrix form as Using the result vec(ABC) (6) can be equivalently expressed as where . We stack the received signals that are transmited on G successive symbols, and we have where y p = y T p,1 , y T p,2 , · · · , y T p,G G UE andn p = n T p,1 , n T p,2 , · · · , n T p,G T ∈ C GN s ×1 . As shown above, in the channel estimation, the measurement matrix Φ p derives from pilots, precoders, combiners and the partial DFT matrices, so it is a difficult matrix in all probability. With the difficult measurement matrix, the existing sparse channel estimation algorithms can fail to recover channels. As a result, the UTAMP-SBL which is robust to difficult measurement matrix is a wise choice.

Uplink Channel Estimation with UTAMP-SBL
Focus on the channel estimation of the single subcarrier and the subscript p is omitted for conciseness, y = Φh ω +n.
Let M = GN s , N = KG BS G UE , so the sizes of y, Φ and h ω are M × 1, M × N and N × 1, respectively, where M N. In order to characterize the sparsity of h ω , the Gaussian prior model is adopted, i.e., where p h ω (h ω |γ) = N C (h ω ; 0, diag(γ)) and p γ (γ) = Ga(γ; , η). Define the singular value decomposition (SVD) of Φ is Φ = UΣV, where U is an orthogonal matrix and U H U = I. Perform unitary transformation to (9) and we have where r = U H y ∈ C M×1 , w = U Hn ∈ C M×1 and w follows N C (w; 0, σI). Define z = ΣVh ω and the likelihood function of z is p r (r|z; σ) = N C (r; z, σI).
According to the Bayesian rule, the joint probability density function is where the Dirac function factor δ(·) is adopted first in [20] to facilitate the derivation of the algorithm. The factor graph of the factorization (13) is shown in Figure 2. The UTAMP-SBL algorithm will be derived from message-passing on this graph and applied to the channel estimation. Firstly, we consider the massages between node p r and node z. According to the sumproduct algorithm, the massage from function node p r to variable node z is m p r →z (z) = p r (r|z; σ) = N C (z; r,σI). The massage from variable node z to function node p r is defined as m z→p r (z) = N C z; µ, diag(τ µ ) . So the belief b(z) = m z→p r (z)m p r →z (z) can be calculated as b(z) = N C (z;ẑ, diag(v z )), where the mean and variance are v z = 1 1 τ µ + 1 Then, we focus on the massages between node p h ω and node h ω . Similarly, the massage from function node p h ω to variable node h ω is defined as m p h ω →h ω (h ω ) = N C (h ω ; 0, diag(γ)). And the massage from variable node h ω to function node p h ω is defined as where the mean and variance are For the sparse channel estimation, the mean of the belief b(h ω ) is the estimated channel. Next, the variable µ, τ µ , π and τ π mentioned above can be updated with the UTAMP which is derived from the AMP. In the UTAMP, the means and variances are calculated by [16] where the intermediate variables s and τ s are given by [16] Besides, there are several model parameters need to be learned at each iteration. For the noise precisionσ, it can be updated with the expectation-maximization (EM) algorithm [16] For the variable γ, according to the sum-product algorithm, the componentwise massage from function node p h ω n to variable node γ n is The message from variable node γ to function node So,γ n is the mean of this distribution, i.e., when we set η = 0. For the parameterˆ , it is difficult to renewed, but there is an effective way found in [16] and we use it directly. As a result,ˆ is calculated bŷ The process of the channel estimation with UTAMP-SBL is given in Algorithm 1. In this algorithm, the iteration will be halted if the number of it reaches the maximum value T max or ĥ ω(t) −ĥ ω(t−1) / ĥ ω(t) < 10 −12 .

: until halt
The computational complexity of the UTAMP-SBL is dominated by the complex multiplications required for matrix-vector operations. At the stage of initialization, the complexity of SVD is O(M 2 N) [21]. The matrix-vector product of U H y is O(M 2 ). At the stage of iteration, the complexity of calculating τ µ and µ are respectively O(M 2 + 2M) and O(2MN + M) in the step 3 of Algorithm 1. The calculations of remaining steps only involve the component-wise vector multiplication or scalar operations which bring a small amount of complex multiplication compared with matrix-vector product. Therefore, the computational complexity of them is omitted. When the iteration reaches T, the total complexity is O(M 2 N + M 2 + T(M 2 + 2M + 2MN + M)). Discarding the low-order terms, the total complexity is O(M 2 N).

Cramér-Rao Bound Analysis
In this section, the CRB is provided as a performance benchmark of the channel estimation. From Section 3, the sparse channel estimation is modeled as an SBL problem and the CRBs for SBL derived in [22] can be adopted directly.

Simulation Results
The simulation parameters are chosen as follows, N BS = G BS = 64, N UE = G UE = 16, N RF BS = N RF UE = 4, N s = 4, K = 4, P = 32 and L = 4. The pilots, precoders and combiners are generated according to [17]. In order to evaluate the performance, the normalized MSE between the estimated channelĥ ω and the real channel h ω is introduced. We compare the NMSE performance for the OMP [5], SBL [9], EM-BBG-VAMP [7], SuRe-CSBL [15], TurboCS [10], UTAMP-SBL, and the CRB derived in Section 4. Figure 3 shows the NMSE versus the pilot overhead at SNR = 20 dB. As defined in Section 3, N is the dimension of the received signal y and M is the dimension of the sparse channel h ω . In this simulation, the pilot overhead is denoted as M/N. We can find that the off-grid algorithm, SuRe-CSBL, could not achieve the desired performance when the pilot is insufficient in spite of it is the winners in [15]. The SBL, EM-BBG-VAMP and UTAMP-SBL work well with 9.38% pilot overhead compared with other algorithms. Thus, when the pilot overhead is extremely low, the SuRe-CSBL, TurboCS and OMP are not sensible choices. It is also shown that the UTAMP-SBL performs distinctly better than the other estimators with the same pilot overhead. From another side, it is also concluded that the UTAMP-SBL can accurately estimate the channels with fewer pilots. For example, when the NMSE performance meets −10 dB, the OMP, SBL and EM-BBG-VAMP respectively need 22.66%, 21.88%, 16.41% pilot overhead. However, the UTAMP-SBL only needs 11.72% overhead, which means the pilot overhead is reduced by 28.58% compared with the EM-BBG-VAMP.  Figure 4 shows the NMSE versus the SNR with 25.00% pilot overhead. Due to insufficient of pilots, the SuRe-CSBL and TurboCS could not estimate the channels well, even if their NMSEs are lower than CRB at low SNRs. It is observed that the UTAMP-SBL has the best performance compared with the other competitors. This is because the UTAMP-SBL is robust to different types of difficult measurement matrices, such as non-zero mean, rankdeficient, correlated, or ill-conditioned matrix [16]. This feature is crucial for the practical application of this algorithm. In the channel estimation, the measurement matrix Φ p is usually related to the complex signal processing operations as described in Section 2.3, so it can be a difficult matrix. Therefore, the UTAMP-SBL with strong robustness performs better in the channel estimation. Moreover, the matrix inversion step of SBL causing the high complexity is avoided by replacing the E-step in the EM with UTAMP, and the algorithm complexity is reduced [16].   Figure 5 provides the NMSE versus the number of iterations at different SNRs with 25.00% pilot overhead, which illustrates the convergences of the algorithms. Since the algorithms, SuRe-CSBL, TurboCS, are not suitable for our channel estimation when the pilot overhead is 25.00%, only the convergences of the OMP, SBL, EM-BBG-VAMP and UTAMP-SBL are compared. When SNR = −10 dB, the OMP fails to converge even if the number of iterations reaches 300, and the NMSE performance decreases with the increasing iterations. This shows that the OMP is not suitable for the low SNR cases. The SBL, EM-BBG-VAMP and UTAMP-SBL require about 100 iterations to converge. When SNR = 10 dB, the OMP is still with the worst convergence and the UTAMP-SBL is with the best convergence. However, when SNR = 30 dB, the EM-BBG-VAMP converges slightly faster than the UTAMP-SBL. From the results, the EM-BBG-VAMP converges within 20 iterations while the UTAMP-SBL requires about 40 iterations. And these two algorithms are both faster than the SBL and the OMP. Additionally, in this figure, only the performance of the UTAMP-SBL becomes better with the iteration at all SNRs, while that of other algorithms sometimes deteriorates with the iteration. This further demonstrates the advantage of the UTAMP-SBL in convergence. From above observation and analysis, the UTAMP-SBL is superior to the comparison algorithms in the convergence speed and the NMSE performance. OMP 10dB SBL 10dB EM-BBG-VAMP 10dB UTAMP-SBL 10dB OMP 30dB SBL 30dB EM-BBG-VAMP 30dB UTAMP-SBL 30dB 10dB 30dB Figure 5. The NMSE versus the number of iterations at different SNRs with 25.00% pilot overhead.

Conclusions
In this paper, the UTAMP-SBL algorithm is applied to estimate the sparse channels for improving the performance. In order to evaluate the effect of this algorithm in practical channel estimation, we study the multi-user uplink channel estimation for hybrid architecture mmWave massive MIMO systems. Simulation results demonstrate the superiority of the UTAMP-SBL in terms of NMSE performance and pilot overhead reduction.