## Gradient Pursuit-Based Channel Estimation for MmWave Massive MIMO Systems with One-Bit ADCs

In-soo Kim and Junil Choi Department of Electrical Engineering POSTECH Pohang, Korea Email: { insookim, junil } @postech.ac.kr

Abstract -In this paper, channel estimation for millimeter wave (mmWave) massive multiple-input multiple-output (MIMO) systems with one-bit analog-to-digital converters (ADCs) is considered. In the mmWave band, the number of propagation paths is small, which results in sparse virtual channels. To estimate sparse virtual channels based on the maximum a posteriori (MAP) criterion, sparsity-constrained optimization comes into play. In general, optimizing objective functions with sparsity constraints is NP-hard because of their combinatorial complexity. Furthermore, the coarse quantization of one-bit ADCs makes channel estimation a challenging task. In the field of compressed sensing (CS), the gradient support pursuit (GraSP) and gradient hard thresholding pursuit (GraHTP) algorithms were proposed to approximately solve sparsity-constrained optimization problems iteratively by pursuing the gradient of the objective function via hard thresholding. The accuracy guarantee of these algorithms, however, breaks down when the objective function is ill-conditioned, which frequently occurs in the mmWave band. To prevent the breakdown of gradient pursuit-based algorithms, the band maximum selecting (BMS) technique, which is a hard thresholder selecting only the 'band maxima,' is applied to GraSP and GraHTP to propose the BMSGraSP and BMSGraHTP algorithms in this paper.

ADCs were proposed. In [10], a sparse Bayesian learningbased (SBL-based) approximate maximum a posteriori (MAP) channel estimator was proposed, which sought the channel estimate using the variational Bayesian (VB) method. The generalized expectation consistent signal recovery (GEC-SR) [11] and generalized approximate message passing (GAMP) [12] algorithms are iterative approximate minimum mean squared error (MMSE) estimators, which are based on the turbo principle and loopy belief propagation (BP). However, the accuracy guarantee of these algorithms breaks down when the sensing matrix is ill-conditioned, which occurs when the angular domain grid resolution of the virtual channel representation is too high.

## I. INTRODUCTION

The wide bandwidth of the millimeter wave (mmWave) band provides high data rates, resulting in a significant performance gain [1]-[4]. Also, the small wavelength enables the use of large arrays at the transmitter and receiver, which is widely known as massive multiple-input multiple-output (MIMO) systems. The power consumption of an analog-to-digital converter (ADC), however, increases linearly with the sampling rate and exponentially with the ADC resolution, which makes highresolution ADCs impractical for mmWave massive MIMO systems [5]. To reduce the impractically high power consumption, one possible solution is to use low-resolution ADCs, which recently gained popularity [6]-[9]. In this paper, we consider the extreme scenario of mmWave massive MIMO systems with one-bit ADCs.

The number of propagation paths in the mmWave band is small, which naturally leads to sparse virtual channels. Therefore, channel estimation for mmWave massive MIMO systems with one-bit ADCs can be formulated as sparse signal recovery with quantized measurements. In [10]-[12], compressed sensing-based (CS-based) signal recovery techniques for mmWave massive MIMO systems with low-resolution

In this paper, we propose gradient pursuit-based iterative approximate MAP channel estimators for mmWave massive MIMO systems with one-bit ADCs. In the mmWave band, the MAP channel estimation framework can be cast into a sparsity-constrained optimization problem. To approximately solve such problem iteratively, we adopt the gradient support pursuit (GraSP) [13] and gradient hard thresholding pursuit (GraHTP) [14] algorithms, which are gradient pursuit-based CS techniques. Similar to the aforementioned CS-based algorithms, however, these algorithms break down when the grid resolution is too high because the resulting objective function is ill-conditioned. To prevent such breakdown, we propose the band maximum selecting (BMS) technique, which is a hard thresholder selecting only the 'band maxima.' The BMS technique is then applied to GraSP and GraHTP, which results in the proposed BMSGraSP and BMSGraHTP algorithms. According to the simulation results, BMSGraSP and BMSGraHTP outperform other channel estimators including GAMP, which shows the superiority of our proposed techniques.

Notation: a , a , and A denote a scalar, vector, and matrix. The transpose, conjugate transpose, and conjugate of A are denoted as A T , A H , and A . The Kronecker product of A and B is denoted as A B ⊗ . The support of a is denoted as supp( a ) , which is formed by collecting all of the indices of the nonzero elements of a . The best s -term approximation of a is denoted as a | s , which is obtained by leaving only the s largest elements of a and hard thresholding other elements to 0 . For a set A , the vector obtained by leaving only the elements of a indexed by A and hard thresholding the remaining elements to 0 is denoted as a | A . The absolute value of a scalar a and cardinality of a set A are denoted as | a | and |A| . The standard normal PDF and CDF are denoted as φ x ( ) = 1 √ 2 π e -x 2 2 and Φ( x ) = ∫ x -∞ √ 1 2 π e -y 2 2 dy . The inverse Mills ratio function is defined as λ x ( ) = φ x ( ) Φ( x ) . The element-wise matrix multiplication and division are denoted as /circledot and /circledivide . The element-wise standard normal PDF, CDF, and inverse Mills ratio function are φ ( x ) , Φ( x ) , and λ ( x ) = φ ( x ) /circledivide Φ( x ) .

## II. SYSTEM MODEL

Consider a mmWave massive MIMO system with one-bit ADCs, whose N -antenna transmitter and M -antenna receiver are equipped with uniform linear arrays (ULAs). In the channel estimation phase, a training sequence of length T is transmitted to the receiver. By collecting the signals over the T time slots, the received signal Y ∈ C M T × is

$$\tilde { \nu } = \sqrt { \rho } \text{HS} + N$$

where H ∈ C M N × is the channel, S ∈ C N × T is the training sequence, whose columns obey the 2 -norm constraint of √ N , and N ∈ C M T × is the additive white Gaussian noise (AWGN), which is distributed as vec( N ) ∼ CN ( 0 MT , I MT ) . The signalto-noise ratio (SNR) is defined as ρ . In the mmWave band, H contains a small number of paths where each path is associated with its path gain, angle-of-arrival (AoA), and angleof-departure (AoD) [15]. Then, H is

$$\mathbf H = \sum _ { \ell = 1 } ^ { L } \alpha _ { \ell } \mathbf a _ { \mathrm R X } ( \theta _ { \mathrm R X, \ell } ) \mathbf a _ { \mathrm T X } ( \theta _ { \mathrm T X, \ell } ) ^ { H } \quad \quad ( 2 ) \\ \text{re $L$ is the number of maths. $\alpha_{\varepsilon}$ $\sim \mathcal{CN}(0.1)$ is the $\ell$}$$

$$\dots, \dots, \nu, \\ \mathbf a _ { \text{RX} } ( \theta _ { \text{RX}, \ell } ) & = \frac { 1 } { \sqrt { M } } \left [ 1 \ \cdots \ \ e ^ { - j \pi ( M - 1 ) \sin ( \theta _ { \text{RX}, \ell } ) } \right ] ^ { \text{T} }, \ \ ( 3 ) \ \text{tr} \ \text{thed} \\ \mathbf a _ { \text{TX} } ( \theta _ { \text{TX}, \ell } ) & = \frac { 1 } { \sqrt { N } } \left [ 1 \ \cdots \ \ e ^ { - j \pi ( N - 1 ) \sin ( \theta _ { \text{TX}, \ell } ) } \right ] ^ { \text{T} }, \ \ ( 4 ) \ \text{and} \ \ B _ { \text{T} } \\ \quad \ \text{it is a b}$$

where L is the number of paths, α /lscript ∼ CN (0 , 1) is the /lscript -th path gain, and θ RX ,/lscript ∼ unif([ -π/ , π/ 2 2]) and θ TX ,/lscript ∼ unif([ -π/ , π/ 2 2]) are the /lscript -th AoA and AoD, which are independent. The steering vectors a RX ( θ RX ,/lscript ) ∈ C M and a TX ( θ TX ,/lscript ) ∈ C N are

$$\iota _ {, \ell } ) = \frac { 1 } { \sqrt { N } } \left [ 1 \ \dots \ \ e ^ { - j \pi ( N - 1 ) \sin ( \theta _ { T X, \ell } ) } \right ] ^ { T }, \ \ ( \\ \dashv \text{-element spacings are half-wavelength. The one-l}$$

$$\hat { Y } & = Q ( \mathbf Y ) & \quad \stackrel { \cdot \cdot } { \text{on tl} } \\ & = Q ( \sqrt { \rho } \mathbf H \mathbf S + \mathbf N ) & \quad \stackrel { \cdot \cdot } { \text{on tl} }$$

whose inter-element spacings are half-wavelength. The one-bit quantized received signal ˆ Y ∈ C M T × is where Q( ) · is the zero threshold one-bit quantization function defined as

$$Q ( Y ) = \sin ( R e ( Y ) ) + j \sin ( I m ( Y ) ) \quad \ \ ( 6 ) \, \stackrel { \chi \chi \chi } { \text{object} }$$

with sign( ) · the element-wise sign function.

The virtual channel representation [16] of H is

$$H = A _ { R X } X ^ { * } A _ { T X } ^ { H }$$

where A RX ∈ C M B × RX and A TX ∈ C N × B TX are overcomplete discrete Fourier transform (DFT) matrices, and X ∗ ∈

C B RX × B TX is the virtual channel with B RX ≥ M and B TX ≥ N . For notational simplicity, we define B = B RX B TX . To facilitate the analysis, (5) is vectorized after plugging in (7) as

$$\hat { y } & = Q ( y ) \\ & = Q ( \sqrt { \rho } \mathbf A \mathbf x ^ { * } + \mathbf n )$$

where ˆ = vec( ˆ ) y Y , y = vec( Y ) , A = S A T TX ⊗ A RX = [ a 1 · · · a B ] , x ∗ = vec( X ∗ ) = [ x ∗ 1 · · · x ∗ B ] T , and n = vec( N ) .

To build the MAP channel estimation framework, we start with formulating the likelihood function based on the real forms of ˆ y , A , and x ∗ , which are

$$\begin{array} { l l l } \text{of} \, \hat { y }, A, \text{ and} \, x ^ { * }, \text{ which are} \\ \hat { y } _ { R } = [ \text{Re} ( \hat { y } ) ^ { T } & \text{Im} ( \hat { y } ) ^ { T } ] ^ { T } \\ = \left [ \hat { y } _ { R, 1 } & \cdots & \hat { y } _ { R, 2 M T } \right ] ^ { T }, & ( 9 ) \\ A _ { R } = \left [ \text{Re} ( A ) & - \text{Im} ( A ) \right ] \\ \text{Im} ( A ) & \text{Re} ( A ) \right ] \\ = \left [ a _ { R, 1 } & \cdots & a _ { R, 2 M T } \right ] ^ { T }, & ( 10 ) \\ x ^ { * } _ { R } = \left [ \text{Re} ( x ^ { * } ) ^ { T } & \text{Im} ( x ^ { * } ) ^ { T } \right ] ^ { T } & ( 11 ) \\ \text{the complex forms and their real forms are used inter-} \\ \text{ably in the sequel. For example,} \, x ^ { * } \text{ and} \, x ^ { * } _ { R } \text{ represent} \\ \end{array}$$

where the complex forms and their real forms are used interchangeably in the sequel. For example, x ∗ and x ∗ R represent the same entity. Then, the log-likelihood function f ( x ) is [17]

$$e \text{ same empty. } \text{ then, use log-timevolume function } J \left ( x \right ) \text{ is } [ 1 / j ] \\ f ( x ) & = \log \Pr \left [ \hat { y } = Q ( \sqrt { \rho } \hat { x } + \mathbf n ) \, | \, x \right ] \\ & = \sum _ { i = 1 } ^ { 2 M T } \log \Phi ( \sqrt { 2 \rho } \hat { j } _ { R, i } \mathbf a _ { R, i } ^ { T } x _ { R } ). \quad \text{ (12)} \\ \text{over, } x ^ { * } \text{ is approximately sparse due to the leakage effect}$$

However, x ∗ is approximately sparse due to the leakage effect [18], which complicates the channel estimation problem. To formulate the MAP channel estimation framework, we assume that x ∗ is L -sparse with independent and identically distributed (i.i.d.) CN (0 , 1) elements, whose locations are uniformly distributed in { 1 , . . . , B } . The mismatch between our assumption on x ∗ and its true distribution reduces as we increase B RX and B TX due to the reduced leakage effect. In Section IV, it is shown that sufficiently large B RX and B TX lead to accurate channel estimation because our assumption closely approximates the true distribution of x ∗ . With such assumption on the distribution of x ∗ , the MAP estimate of x ∗ is

$$\underset { x \in C ^ { B } } { \arg \max } \, \left ( f ( x ) + g ( x ) \right ) \, \text{s.t.} \, \| \mathbf x \| _ { 0 } \leq L \quad \ \ ( 1 3 ) \\.$$

where g ( x ) = -‖ x R ‖ 2 is the logarithm of the PDF of CN ( 0 B , I B ) without the constant factor. In the sequel, the objective function and its gradient in (13) are denoted as

$$h ( x ) & = f ( x ) + g ( x ), \quad \quad ( 1 4 ) \\ \nabla h ( x ) & = \nabla f ( x ) + \nabla g ( x ) \\ & = \left [ \nabla h ( x _ { 1 } ) \ \cdots \ \ \nabla h ( x _ { B } ) \right ] ^ { \text{T} }, \quad \ \ ( 1 5 ) \\ \text{ose differentiation is with respect to $x$.}$$

whose differentiation is with respect to x .

## III. PROPOSED BMSGRASP AND BMSGRAHTP ALGORITHMS

In general, solving (13) is NP-hard due to its sparsity constraint. To approximately optimize sparsity-constrained objective functions iteratively by pursuing the gradient of the objective function, GraSP [13] and GraHTP [14] algorithms were proposed in the field of CS, which generalize the well-known compressive sampling matching pursuit (CoSaMP) [19] and hard thresholding pursuit (HTP) [20] algorithms to objective functions with arbitrary forms.

To solve (13), GraSP and GraHTP update the L -sparse current estimate ˆ x of x ∗ roughly as follows at each iteration. The best L -term approximation of ∇ h (ˆ) x is computed via hard thresholding as

$$\nabla h ( \hat { \mathbf x } ) | _ { L }, \quad \quad \ \ ( 1 6 ) \ \text{Co}$$

whose support

$$\mathcal { I } = \text{supp} ( \nabla h ( \hat { \mathbf x } ) | _ { L } ) \, \quad \, \quad ( 1 7 ) \, \stackrel { \cdots } { \mathbf r } \, \text{ Graf} \, \text{ 1 7 } \, \text{ }$$

is selected as the updated estimate of supp( x ∗ ) . This step can be interpreted as support identification, or joint AoA and AoD estimation. Then, the elements corresponding to the identified support are updated by selecting the updated ˆ x as

$$\underset { x \in \mathbb { C } ^ { B } } { \arg \max } \, h ( x ) \ s. t. \ s u p ( x ) \subseteq \mathcal { I }, \quad \ \ ( 1 8 ) \quad 2 ) \underset { \quad \ \text{where} } { \sum } \,. \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \colon \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \ddots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots$$

which can be interpreted as path gain estimation. Note that (18) is a convex optimization problem because h ( x ) is concave with its support constraint being convex. The concavity of h ( x ) follows from the fact that f ( x ) is concave because Φ( ) · is log-concave; similarly, g ( x ) is concave due to the convexity of 2 -norm [21].

The accuracy guarantee of GraSP and GraHTP, however, breaks down when h ( x ) does not have a stable restricted Hessian [13] or is not strongly convex and smooth [14]. The problem is that for large B RX and B TX , the resulting A is ill-conditioned because its columns become highly coherent, which leads to unfavorable h ( x ) . Specifically, with highly coherent A , GraSP and GraHTP are likely to fail to accurately identify the support of x ∗ from (17).

To illustrate how support identification fails from (17) when A is highly coherent, first, consider the real form ∇ h ( x R ) of ∇ h ( x ) defined as

$$& \vee h ( x ) \text{ defined as } & \text{ hard } \text{th} \\ & \quad \nabla h ( x _ { R } ) & \text{ product} \\ & \quad = \left [ R e ( \nabla h ( x ) ) ^ { T } \ \text{ Im} ( \nabla h ( x ) ) ^ { T } \right ] ^ { T } & \text{ T _ { BMS, L } } \\ & \quad = \sum _ { i = 1 } ^ { 2 M T } \lambda ( \sqrt { 2 \rho } \hat { y } _ { R, i } a _ { R, i } ^ { T } x _ { R } ) \sqrt { 2 \rho } \hat { y } _ { R, i } a _ { R, i } - 2 x _ { R } } & \text{ present} \\ & \quad = \text{$A_{R}^{\prime}(\lambda(\sqrt{2\rho} \hat { y } _ { R } \odot A_{R} x _ { R } ) \odot \sqrt{2\rho} \hat { y } _ { R } ) - 2 x _ { R }, \quad ( 1 9 ) & \text{ as the } \text{ index } i \\ & \text{which is obtained from } \nabla \log \Phi ( a _ { R } ^ { T } x _ { R } ) \, = \, \lambda ( a _ { R } ^ { T } x _ { R } ) a _ { R } \text{ and } \text{ check} \\ & \nabla \| x _ { R } \| ^ { 2 } = 2 x _ { R }. \text{ Then, we establish the following observation, } \text{ this pair} \\ & \quad \dots \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cd.$$

which is obtained from ∇ log Φ( a T R x R ) = λ ( a T R x R ) a R and ∇‖ x R ‖ 2 = 2 x R . Then, we establish the following observation, which can be checked from directly computing ∇ h x ( i ) , whose real and imaginary parts correspond to the i -th and ( i + B ) -th elements of ∇ h ( x R ) .

$$\Gamma \quad & \text{Observation $1$. $\nabla h(x_{i})=\nabla h(x_{j})$ if $a_{i}=\mathfrak{a}_{j}$ and $x_{i}=x_{j}$.} \\ \quad \ \. \quad \. \quad \. \quad \. \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \ddots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots$$

/negationslash

However, Observation 1 is trivial since a i = a j unless i = j . To introduce a nontrivial observation, we adopt the notion of coherence between a i and a j , i.e.,

$$\mu ( i, j ) = \frac { | a _ { i } ^ { H } a _ { j } | } { \| a _ { i } \| \| a _ { j } \| }, \quad \text{(2)} \\ \text{rec the non-unimit} \, \text{that was an} \, \circ \, \text{and} \, \circ \, \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Gamma \Delta$$

which measures the proximity between a i and a j [22]-[24]. Then, with the η -coherence band of i , which is defined as [22]

$$B _ { \eta } ( i ) = \{ j \ | \ \mu ( i, j ) \geq \eta \} \quad \quad \quad ( 2 1 )$$

with η ∈ (0 , 1) , the following conjecture is established for sufficiently large η .

$$\text{Conjecture } 1. \ \nabla h ( x _ { i } ) \approx \nabla h ( x _ { j } ) \ i f \ j \in B _ { \eta } ( i ) \ a n d \ x _ { i } = x _ { j }. \\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots.$$

Now, based on Conjecture 1, we illustrate how GraSP and GraHTP fail to accurately identify the support of x ∗ from (17) when A is highly coherent. To proceed, we consider the following example where x ∗ and ˆ y are realized with the current estimate ˆ x so as to satisfy

$$1 ) \ i = \underset { \substack { \arg \max \\ k \in \{ 1, \dots, B \} } } { \arg \max } | \nabla h ( \hat { x } _ { k } ) | \\ 2 ) \ \mathcal { J } _ { \eta } ( i ) \cap \text{supp} ( \mathbf x ^ { * } ) = \emptyset \\ \text{where}$$

where

$$\mathcal { J } _ { \eta } ( i ) = \{ j \ | \ j \in B _ { \eta } ( i ), \hat { x } _ { i } = \hat { x } _ { j } \} \, \langle \{ i \} \quad \ \ ( 2 2 )$$

is defined as the by-product of i . In this example, J η ( ) i is called the by-product of i because for all j ∈ J η ( ) i ,

$$| \nabla h ( \hat { x } _ { j } ) | & \overset { ( a ) } { \approx } | \nabla h ( \hat { x } _ { i } ) | \\ & \overset { ( b ) } { = } \max _ { k \in \{ 1, \dots, B \} } | \nabla h ( \hat { x } _ { k } ) | \quad \text{(3)} \\ \text{where } \omega \text{ follow from Confucture } 1 \text{ and } 1 \text{$\mathfrak{w}$ from } 1 \text{}$$

holds where (a) follows from Conjecture 1 and (b) from 1) even if J η ( ) i ∩ supp( x ∗ ) = ∅ according to 2). The implication of this example is that most indices of J η ( ) i are likely to be selected when ∇ h (ˆ) x is hard thresholded as in (17) due to (23), which results in an erroneous estimate of supp( x ∗ ) because of 2). Based on the observation established from this example, we propose the BMS technique to remedy such problem.

The BMS technique provides a guideline of how to hard threshold ∇ h (ˆ) x in (17) in order to exclude the byproduct indices. The proposed BMS hard thresholding function T BMS ,L ( ) · is an L -term hard thresholding function developed based on Conjecture 1. The details of the BMS technique are presented in Algorithm 1. In Line 3, the index of the maximum element of ∇ h (ˆ) x is selected among the unchecked index set as the current index. The by-product testing set of the current index is formed in Line 4. In Line 5, the current index is checked whether it is greater than the by-product testing set. In this paper, we refer to Line 5 as the band maximum criterion. If the current index is indeed the 'band maximum,' which satisfies the band maximum criterion, this index is selected as the estimate of supp( x ∗ ) in Line 6. Otherwise, the current index is excluded because it is likely to be the by-product of

## Algorithm 1 BMS hard thresholding technique

$$\text{Algorithm} \, 1 \, BMS hard thresholding technique \, \text{$\frac{1}{2}$} \\ \text{Input} \colon \, \dot { x }, \nabla h ( \dot { x } ), L \\ \text{Output} \colon \, T _ { \text{BMS}, L } ( \nabla h ( \dot { x } ) ) \, \text{$1$} \, \mathcal { S } = \emptyset, \mathcal { I } = \{ 1, \dots, B \} \\ \text{2$} \, \text{while} \, | \mathcal { S } | < L \, \text{do} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \\ \text{another index rather than the ground truth index. In Line} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$}$$

another index rather than the ground truth index. In Line 8, the unchecked index set is updated.

At this point, we emphasize that Algorithm 1 is applied to ∇ h (ˆ) x . To apply the BMS hard thresholding function to ˆ + x κ ∇ h (ˆ) x where κ is the step size, simply replace ∇ h (ˆ) x with ˆ + x κ ∇ h (ˆ) x in the input, output, and Lines 3, 5, and 10 of Algorithm 1. This variant can be derived based on the same logic using Conjecture 1. Now, the BMSGraSP and BMSGraHTP algorithms are proposed to solve (13).

BMSGraSP and BMSGraHTP are the variants of GraSP and GraHTP. The difference between our BMS-based and nonBMS-based algorithms is that T BMS ,L ( ) · is used as a hard thresholder instead of the naive best L -term approximation as in (17). The details of the proposed BMSGraSP and BMSGraHTP are given in Algorithms 2 and 3. Lines 3, 4, and 5 of Algorithms 2 and 3 proceed based on the same logic. In Line 3, the gradient of the objective function is computed. Then, I is selected from the support of the hard thresholded gradient of the objective function in Line 4, which corresponds to joint AoA and AoD estimation. In Line 5, the objective function is maximized subject to the support constraint, which can be interpreted as path gain estimation. This step can be solved via convex optimization since the objective function is concave with the support constraint being convex. Additionally, Line 6 of Algorithm 2 hard thresholds b since the sparsity of b is at most 3 L . A natural halting condition for Algorithms 2 and 3 is to halt when supp(ˆ) x does not change from iteration to iteration [13], [14].

Remark 1: Instead of merely hard thresholding b in Line 6 of Algorithm 2, we can solve the following convex optimization problem

$$\hat { x } = & \underset { x \in \mathbb { C } ^ { B } } { \arg \max } \, h ( x ) \ s. t. \ \text{supp} ( x ) \subseteq \underset { \real } { \longrightarrow } \, \text{supp} ( b | _ { L } ) \quad ( 2 4 ) \\ \dots \dots \wedge \, \wedge \, \pi \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \,.$$

to update ˆ x . This is called the debiasing variant of Algorithm 2, which produces a more accurate x ∗ [13]. The complexity, however, increases.

Remark 2: The convex optimization problems in BMSGraSP and BMSGraHTP, which consist of Line 5 of Algorithms 2 and 3, can be solved with relatively low com-

## Algorithm 2 BMSGraSP algorithm

Input:

Output: ˆ x h ( ) · , L

1: ˆ = x 0 B

- 2: while halting condition do
- 3: z = ∇ h (ˆ) x
- 4: I = supp( T BMS 2 , L ( z )) ∪ supp(ˆ) x

x

- 5: b = argmax B h ( x ) s.t. supp( x ) ⊆ I

C

6:

∈

ˆ =

x

b

L

- 7: end while

|

## Algorithm 3 BMSGraHTP algorithm

Input: · Output: ˆ x

h

( )

,

L

1: ˆ = x 0 B

2: while halting condition do

3: z = ˆ + x κ ∇ h (ˆ) x

4: I = supp( T BMS ,L ( z ))

5: ˆ = argmax x h ( x ) s.t. supp( x ) ⊆ I

x ∈ C B

6: end while

plexity because the support of their optimization variables are constrained to I where |I| = O L ( ) is typically much smaller than B ≥ MN in mmWave massive MIMO systems; L is small due to the small number of paths, whereas M and N are large due to the large arrays. The same logic holds for (24). Therefore, the complexity of Algorithms 2 and 3 is dominated by Line 3, which requires the gradient of the objective function defined on C B . We can reduce the complexity in Line 3 of Algorithms 2 and 3 based on the fast Fourier transform (FFT) implementation of ∇ h (ˆ) x when S has an FFT-friendly structure, e.g., circularly shifted Zadoff-Chu (ZC) sequences or partial DFT matrix. From (19), note that the matrix-vector multiplications involving A and A H act as computational bottlenecks. These matrix-vector multiplications can be efficiently implemented using the FFT by noting that

$$\cdots \, \dots \, \dots \, \cdots \, \cd$$

where C = unvec( ) c . By performing matrix multiplications involving A RX , A TX , and S using the FFT, the cost of computing ∇ h (ˆ) x can be reduced significantly.

## IV. SIMULATION RESULTS

In this section, the performance of BMSGraSP and BMSGraHTP are evaluated in terms of the accuracy and achievable rate. We consider a mmWave massive MIMO system with onebit ADCs where M = N = 64 and T = 80 . The rows of the training signal S are chosen as circularly shifted ZC sequences of length T [25], and the number of paths are L = 4 . We

Fig. 1. NMSE vs. SNR with M = N = 64 , T = 80 , and L = 4 , while B RX and B TX vary from algorithm to algorithm.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00508v1-figure-9425d08c-6d5b-4b4d-8c61-8497908ef952.jpg)

choose B RX = 256 /greatermuch M and B TX = 256 /greatermuch N for the proposed BMSGraSP and BMSGraHTP to reduce the leakage effect, but other channel estimators may use smaller B RX and B TX to prevent the breakdown caused by the resulting ill-conditioned sensing matrix. The accuracy of a channel estimator is measured based on its normalized MSE (NMSE), which is defined as

$$\begin{array} {$$

where ˆ H = A RX ˆ XA TX with ˆ X = unvec(ˆ) x .

For BMSGraSP, we consider its debiasing variant, which replaces Line 6 of Algorithm 2 with (24). We set Algorithms 2 and 3 to halt when supp(ˆ) x does not change from iteration to iteration. We use the backtracking line search [21] to compute κ in Line 3 of Algorithm 3. Lastly, we configure η so as to satisfy Conjecture 1 by selecting the maximum η satisfying

$$\min _ { i \in \{ 1, \dots, B \} } | B _$$

The other channel estimators to be compared are GraSP [13] and GraHTP [14], which are non-BMS-based gradient pursuit-based algorithms, the Bernoulli-Gaussian-GAMP (BGGAMP) [12], and fast iterative shrinkage-thresholding algorithm (FISTA) [26]. BG-GAMP is a loopy BP-based iterative approximate MMSE channel estimator, which assumes that x ∗ is distributed as i.i.d. BG, meaning that each elements are CN (0 , 1) with probability L/B but zero otherwise. FISTA is an accelerated proximal gradient descent method-based iterative MAP channel estimator based on the assumption that the logarithm of the PDF of x ∗ is g FISTA ( x ) = - ‖ ‖ γ x 1 without the constant factor, which is the Laplace distribution. Then, the FISTA estimate of x ∗ is

$$\underset { x \in \mathbb { C } ^ { B } } { \arg \$$

which can be interpreted as the generalized LASSO [27] with arbitrary forms of objective functions. For a fair compari-

Fig. 2. Achievable rate lower bound [12] vs. SNR with M = N = 64 , T = 80 , and L = 4 , while B RX and B TX vary from algorithm to algorithm.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00508v1-figure-a02c692f-f0b1-46b3-9374-b9e655494ad7.jpg)

son, we configure the regularization parameter γ so that the expected sparsity of the FISTA estimate is 3 L , a criterion suggested in [13]. Since GraSP, GraHTP, and BG-GAMP break down when the sensing matrix is ill-conditioned, we set B RX = B TX = 64 for these channel estimators. For FISTA, we choose B RX = B TX = 256 as in our BMSGraSP and BMSGraHTP.

The NMSE of various channel estimators is shown in Fig. 1 at different SNRs. In Fig. 1, our channel estimators outperform other channel estimators in the medium and high SNR regimes. The poor performance of GraSP, GraHTP, and BG-GAMP is caused by the leakage effect due to the small B RX and B TX , but increasing these parameters is forbidden since these channel estimators diverge when A is highly coherent. FISTA performs poorly because the Laplace distribution deviates from the true distribution of x ∗ . In contrast, since B RX and B TX are large, BMSGraSP and BMSGraHTP do not suffer from the leakage effect. As a side note, we mention that all channel estimators experience performance degradation as the SNR enters the high SNR regime. This phenomenon is due to the coarse quantization of one-bit ADCs, which results in magnitude information loss. To illustrate this phenomenon, note that x ∗ and c x ∗ are indistinguishable in the high SNR regime for c &gt; 0 because

$$Q ( \sqrt { \rho } \mathbf A x ^ { * } + \math$$

where (a) and (b) follow from the high SNR regime assumption and (6), which implies that the information in c is lost. To combat such performance degradation, the concept of dithering was suggested in [28], but this is beyond the scope of this paper.

We also compare the achievable rate lower bound of various channel estimators in Fig. 2 at different SNRs. This achievable rate lower bound, which was derived in [12], is obtained by selecting the precoders and decoders based on ˆ H and applying the Bussgang decomposition [29] in conjunction with the fact that the Gaussian noise is the worst-case noise. According to Fig. 2, our channel estimators outperform other channel estimators, which agrees with Fig. 1.

## V. CONCLUSION

In this paper, we proposed gradient pursuit-based iterative approximate MAP channel estimators for mmWave massive MIMO systems with one-bit ADCs. In the mmWave band, the MAP channel estimation framework can be cast into a sparsityconstrained optimization problem, which is NP-hard to solve. To approximately solve such problem iteratively by pursuing the gradient of the objective function, the GraSP and GraHTP algorithms were proposed in the field of CS, which generalize the well-known CoSaMP and HTP algorithms. GraSP and GraHTP, however, break down when the objective function is ill-conditioned, which is likely to occur in the mmWave band. As a solution to such breakdown, the BMS technique, which hard thresholds the gradient of the objective function based on the band maximum criterion, was proposed in this paper. The BMS technique was applied to GraSP and GraHTP to produce the proposed BMSGraSP and BMSGraHTP algorithms. The simulation results showed that the BMSGraSP and BMSGraHTP algorithms outperform other channel estimators.

## ACKNOWLEDGMENT

This work was partly supported by Institute for Information &amp; communications Technology Promotion(IITP) grant funded by the Korea government(MSIT) (No. 2016-0-00123, Development of Integer-Forcing MIMO Transceivers for 5G &amp; Beyond Mobile Communication Systems) and by the National Research Foundation(NRF) grant funded by the MSIT of the Korea government (No. 2018R1A4A1025679).

## REFERENCES

- [1] A. L. Swindlehurst, E. Ayanoglu, P. Heydari, and F. Capolino, 'Millimeter-wave massive MIMO: the next wireless revolution?' IEEE Communications Magazine , vol. 52, no. 9, pp. 56-62, Sept. 2014.
- [2] T. S. Rappaport, S. Sun, R. Mayzus, H. Zhao, Y. Azar, K. Wang, G. N. Wong, J. K. Schulz, M. Samimi, and F. Gutierrez, 'Millimeter wave mobile communications for 5G cellular: it will work!' IEEE Access , vol. 1, pp. 335-349, 2013.
- [3] F. Boccardi, R. W. Heath, A. Lozano, T. L. Marzetta, and P. Popovski, 'Five disruptive technology directions for 5G,' IEEE Communications Magazine , vol. 52, no. 2, pp. 74-80, Feb. 2014.
- [4] S. Rangan, T. S. Rappaport, and E. Erkip, 'Millimeter-wave cellular wireless networks: potentials and challenges,' Proceedings of the IEEE , vol. 102, no. 3, pp. 366-385, Mar. 2014.
- [5] B. Le, T. W. Rondeau, J. H. Reed, and C. W. Bostian, 'Analog-to-digital converters,' IEEE Signal Processing Magazine , vol. 22, no. 6, pp. 69-77, Nov. 2005.
- [6] C. Moll´ en, J. Choi, E. G. Larsson, and R. W. Heath, 'Uplink performance of wideband massive MIMO with one-bit ADCs,' IEEE Transactions on Wireless Communications , vol. 16, no. 1, pp. 87-100, Jan. 2017.
- [7] L. Fan, S. Jin, C. Wen, and H. Zhang, 'Uplink achievable rate for massive MIMO systems with low-resolution ADC,' IEEE Communications Letters , vol. 19, no. 12, pp. 2186-2189, Dec. 2015.
- [8] J. Zhang, L. Dai, S. Sun, and Z. Wang, 'On the spectral efficiency of massive MIMO systems with low-resolution ADCs,' IEEE Communications Letters , vol. 20, no. 5, pp. 842-845, May 2016.
- [9] S. Jacobsson, G. Durisi, M. Coldrey, U. Gustavsson, and C. Studer, 'Throughput analysis of massive MIMO uplink with low-resolution ADCs,' IEEE Transactions on Wireless Communications , vol. 16, no. 6, pp. 4038-4051, Jun. 2017.
- [10] Y. Ding, S. Chiu, and B. D. Rao, 'Bayesian channel estimation algorithms for massive MIMO systems with hybrid analog-digital processing and low-resolution ADCs,' IEEE Journal of Selected Topics in Signal Processing , vol. 12, no. 3, pp. 499-513, Jun. 2018.
- [11] H. He, C. Wen, and S. Jin, 'Bayesian optimal data detector for hybrid mmWave MIMO-OFDM systems with low-resolution ADCs,' IEEE Journal of Selected Topics in Signal Processing , vol. 12, no. 3, pp. 469483, Jun. 2018.
- [12] J. Mo, P. Schniter, and R. W. Heath, 'Channel estimation in broadband millimeter wave MIMO systems with few-bit ADCs,' IEEE Transactions on Signal Processing , vol. 66, no. 5, pp. 1141-1154, Mar. 2018.
- [13] S. Bahmani, B. Raj, and P. T. Boufounos, 'Greedy sparsity-constrained optimization,' Journal of Machine Learning Research , vol. 14, no. Mar, pp. 807-841, 2013.
- [14] X.-T. Yuan, P. Li, and T. Zhang, 'Gradient hard thresholding pursuit.' Journal of Machine Learning Research , vol. 18, pp. 166-1, 2017.
- [15] T. A. Thomas, H. C. Nguyen, G. R. MacCartney, and T. S. Rappaport, '3D mmWave channel model proposal,' in 2014 IEEE 80th Vehicular Technology Conference (VTC2014-Fall) , Sept. 2014, pp. 1-6.
- [16] A. M. Sayeed, 'Deconstructing multiantenna fading channels,' IEEE Transactions on Signal Processing , vol. 50, no. 10, pp. 2563-2579, Oct. 2002.
- [17] J. Choi, J. Mo, and R. W. Heath, 'Near maximum-likelihood detector and channel estimator for uplink multiuser massive MIMO systems with one-bit ADCs,' IEEE Transactions on Communications , vol. 64, no. 5, pp. 2005-2018, May 2016.
- [18] J. Mo, P. Schniter, N. G. Prelcic, and R. W. Heath, 'Channel estimation in millimeter wave MIMO systems with one-bit quantization,' in 2014 48th Asilomar Conference on Signals, Systems and Computers , Nov. 2014, pp. 957-961.
- [19] D. Needell and J. A. Tropp, 'CoSaMP: iterative signal recovery from incomplete and inaccurate samples,' Applied and computational harmonic analysis , vol. 26, no. 3, pp. 301-321, 2009.
- [20] S. Foucart, 'Hard thresholding pursuit: an algorithm for compressive sensing,' SIAM Journal on Numerical Analysis , vol. 49, no. 6, pp. 25432563, 2011.
- [21] S. Boyd and L. Vandenberghe, Convex optimization . Cambridge university press, 2004.
- [22] A. Fannjiang and W. Liao, 'Coherence pattern-guided compressive sensing with unresolved grids,' SIAM Journal on Imaging Sciences , vol. 5, no. 1, pp. 179-202, 2012.
- [23] N. Jindal, 'MIMO broadcast channels with finite rate feedback,' in GLOBECOM '05. IEEE Global Telecommunications Conference, 2005. , vol. 3, Nov. 2005, pp. 5 pp.-.
- [24] A. Shrivastava and P. Li, 'Asymmetric LSH (ALSH) for sublinear time maximum inner product search (MIPS),' in Advances in Neural Information Processing Systems , 2014, pp. 2321-2329.
- [25] D. Chu, 'Polyphase codes with good periodic correlation properties (Corresp.),' IEEE Transactions on Information Theory , vol. 18, no. 4, pp. 531-532, Jul. 1972.
- [26] A. Beck and M. Teboulle, 'A fast iterative shrinkage-thresholding algorithm for linear inverse problems,' SIAM journal on imaging sciences , vol. 2, no. 1, pp. 183-202, 2009.
- [27] R. Tibshirani, 'Regression shrinkage and selection via the lasso,' Journal of the Royal Statistical Society. Series B (Methodological) , pp. 267-288, 1996.
- [28] N. Liang and W. Zhang, 'Mixed-ADC massive MIMO,' IEEE Journal on Selected Areas in Communications , vol. 34, no. 4, pp. 983-997, Apr. 2016.
- [29] J. J. Bussgang, 'Crosscorrelation functions of amplitude-distorted Gaussian signals,' 1952.
In-soo Kim Junil Choi
Gradient Pursuit-Based Channel Estimation for MmWave Massive MIMO Systems with One-Bit ADCs
In this paper, channel estimation for millimeter wave (mmWave) massive
multiple-input multiple-output (MIMO) systems with one-bit analog-to-digital
converters (ADCs) is considered. In the mmWave band, the number of propagation
paths is small, which results in sparse virtual channels. To estimate sparse
virtual channels based on the maximum a posteriori (MAP) criterion,
sparsity-constrained optimization comes into play. In general, optimizing
objective functions with sparsity constraints is NP-hard because of their
combinatorial complexity. Furthermore, the coarse quantization of one-bit ADCs
makes channel estimation a challenging task. In the field of compressed sensing
(CS), the gradient support pursuit (GraSP) and gradient hard thresholding
pursuit (GraHTP) algorithms were proposed to approximately solve
sparsity-constrained optimization problems iteratively by pursuing the gradient
of the objective function via hard thresholding. The accuracy guarantee of
these algorithms, however, breaks down when the objective function is
ill-conditioned, which frequently occurs in the mmWave band. To prevent the
breakdown of gradient pursuit-based algorithms, the band maximum selecting
(BMS) technique, which is a hard thresholder selecting only the "band maxima,"
is applied to GraSP and GraHTP to propose the BMSGraSP and BMSGraHTP algorithms
in this paper.