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

In millimeter wave (mmWave) massive multiple-input multiple-output (MIMO) systems, one-bit analog-to-digital converters (ADCs) are employed to reduce the impractically high power consumption, which is incurred by the wide bandwidth and large arrays. In practice, the mmWave band consists of a small number of paths, thereby rendering sparse virtual channels. Then, the resulting maximum a posteriori (MAP) channel estimation problem is a sparsity-constrained optimization problem, which is NP-hard to solve. In this paper, iterative approximate MAP channel estimators for mmWave massive MIMO systems with one-bit ADCs are proposed, which are based on the gradient support pursuit (GraSP) and gradient hard thresholding pursuit (GraHTP) algorithms. The GraSP and GraHTP algorithms iteratively pursue the gradient of the objective function to approximately optimize convex objective functions with sparsity constraints, which are the generalizations of the compressive sampling matching pursuit (CoSaMP) and hard thresholding pursuit (HTP) algorithms, respectively, in compressive sensing (CS). However, the performance of the GraSP and GraHTP algorithms is not guaranteed when the objective function is ill-conditioned, which may be incurred by the highly coherent sensing matrix. In this paper, the band maximum selecting (BMS) hard thresholding technique is proposed to modify the GraSP and GraHTP algorithms, namely the BMSGraSP and BMSGraHTP algorithms, respectively. The BMSGraSP and BMSGraHTP algorithms pursue the gradient of the objective function based on the band maximum criterion instead of the naive hard thresholding. In addition, a fast Fourier transform-based (FFT-based) fast implementation is developed to reduce the complexity. The BMSGraSP and BMSGraHTP algorithms are shown to be both accurate and efficient, whose performance is verified through extensive simulations.


In-soo Kim and Junil Choi
Abstract-In millimeter wave (mmWave) massive multipleinput multiple-output (MIMO) systems, one-bit analog-to-digital converters (ADCs) are employed to reduce the impractically high power consumption, which is incurred by the wide bandwidth and large arrays. In practice, the mmWave band consists of a small number of paths, thereby rendering sparse virtual channels. Then, the resulting maximum a posteriori (MAP) channel estimation problem is a sparsity-constrained optimization problem, which is NP-hard to solve. In this paper, iterative approximate MAP channel estimators for mmWave massive MIMO systems with one-bit ADCs are proposed, which are based on the gradient support pursuit (GraSP) and gradient hard thresholding pursuit (GraHTP) algorithms. The GraSP and GraHTP algorithms iteratively pursue the gradient of the objective function to approximately optimize convex objective functions with sparsity constraints, which are the generalizations of the compressive sampling matching pursuit (CoSaMP) and hard thresholding pursuit (HTP) algorithms, respectively, in compressive sensing (CS). However, the performance of the GraSP and GraHTP algorithms is not guaranteed when the objective function is ill-conditioned, which may be incurred by the highly coherent sensing matrix. In this paper, the band maximum selecting (BMS) hard thresholding technique is proposed to modify the GraSP and GraHTP algorithms, namely the BMSGraSP and BMSGraHTP algorithms, respectively. The BMSGraSP and BMSGraHTP algorithms pursue the gradient of the objective function based on the band maximum criterion instead of the naive hard thresholding. In addition, a fast Fourier transform-based (FFTbased) fast implementation is developed to reduce the complexity. The BMSGraSP and BMSGraHTP algorithms are shown to be both accurate and efficient, whose performance is verified through extensive simulations.

I. INTRODUCTION
In millimeter wave (mmWave) massive multiple-input multiple-output (MIMO) systems, the wide bandwidth of the mmWave band in the range of 30-300 GHz offers a high data rate, which guarantees a significant performance gain [1]- [4]. However, the power consumption of analogto-digital converters (ADCs) is scaled quadratically with the sampling rate and exponentially with the ADC resolution, which renders high-resolution ADCs impractical for mmWave systems [5]. To reduce the power consumption, low-resolution ADCs were suggested as a possible solution, which recently gained popularity [6]- [9]. Coarsely quantizing the received signal using low-resolution ADCs results in an irreversible loss of information, which might cause a significant performance degradation. In this paper, we consider the extreme scenario of using one-bit ADCs for mmWave systems.
In practice, the mmWave band consists of a small number of propagation paths, which results in sparse virtual channels. In the channel estimation point of view, sparse channels are favorable because the required complexity and measurements can be reduced. Sparsity-constrained channel distributions, however, cannot be described in closed forms, which makes it difficult to exploit Bayesian channel estimation. In [10], [11], channel estimators for massive MIMO systems with onebit ADCs were proposed, which account for the effect of the coarse quantization. The near maximum likelihood (nML) channel estimator [10] selects the maximizer of the likelihood function subject to the L 2 -norm constraint as the estimate of the channel, which is solved using the projected gradient descent method [12]. However, the channel sparsity was not considered in [10]. In [11], the Bussgang linear minimum mean squared error (BLMMSE) channel estimator was derived by linearizing one-bit ADCs based on the Bussgang decomposition [13]. The BLMMSE channel estimator is an LMMSE channel estimator for massive MIMO systems with one-bit ADCs, whose assumption is that the channel is Gaussian. Therefore, the sparsity of the channel is not taken into account in [11] either.
To take the channel sparsity into account, iterative approximate MMSE estimators for mmWave massive MIMO systems with one-bit ADCs were proposed in [14], [15]. The generalized expectation consistent signal recovery (GEC-SR) algorithm in [14] is an iterative approximate MMSE estimator based on the turbo principle [16], which can be applied to any nonlinear function of the linearly-mapped signal to be estimated. Furthermore, the only constraint on the distribution of the signal to be estimated is that its elements must be independent and identically distributed (i.i.d.) random variables. Therefore, the GEC-SR algorithm can be used as an approximate MMSE channel estimator for any ADC resolutions ranging from one-bit to high-resolution ADCs. However, the inverse of the sensing matrix is required at each iteration, which is impractical in massive MIMO systems in the complexity point of view.
The generalized approximate message passing-based (GAMP-based) channel estimator for mmWave massive MIMO systems with low-resolution ADCs was proposed in [15], which is another iterative approximate MMSE channel estimator. In contrast to the GEC-SR algorithm, only matrix-vector multiplications are required at each iteration, which is favorable in the complexity point of view. As in the GEC-SR algorithm, the GAMP-based algorithm can be applied to any ADC resolutions and any channel distributions as long as the elements of channel are i.i.d. random variable. The performance of the GEC-SR and GAMP algorithms, however, cannot be guaranteed when the sensing matrix is ill-conditioned, which frequently occurs in the mmWave band. To prevent the sensing matrix from becoming ill-conditioned, the GAMP-based channel estimator constructs the virtual channel representation using discrete Fourier transform (DFT) matrices, whose columns are orthogonal. However, such virtual channel representation results in a large gridding error, which leads to performance degradation.
Our goal is to propose an iterative approximate maximum a posteriori (MAP) channel estimator for mmWave massive MIMO systems with one-bit ADCs. Due to the sparse nature, the MAP channel estimation problem is converted to a sparsity-constrained optimization problem, which is NP-hard to solve [17]. To approximately solve such problems iteratively, the gradient support pursuit (GraSP) and gradient hard thresholding pursuit (GraHTP) algorithms were proposed in [17] and [18], respectively. The GraSP and GraHTP algorithms pursue the gradient of the objective function at each iteration by hard thresholding. These algorithms are the generalizations of the compressive sampling matching pursuit (CoSaMP) [19] and hard thresholding pursuit (HTP) [20] algorithms, respectively, in compressive sensing (CS).
With highly coherent sensing matrix, however, the GraSP and GraHTP algorithms do not perform appropriately since the objective function becomes ill-conditioned. To remedy such break down, we exploit the band maximum selecting (BMS) hard thresholding technique, which is then applied to the GraSP and GraHTP algorithms to propose the BMSGraSP and BMSGraHTP algorithms, respectively. The proposed BMSbased algorithms perform hard thresholding for the gradient of the objective function based on the proposed band maximum criterion, which tests whether an index is the ground truth index or the by-product of another index. To reduce the complexity of the BMS-based algorithms, a fast Fourier transform-based (FFT-based) fast implementation of the objective function and gradient is proposed. The BMSbased algorithms are shown to be both accurate and efficient, which is verified through extensive simulations.
The rest of this paper is organized as follows. In Section II, mmWave massive MIMO systems with one-bit ADCs are described. In Section III, the MAP channel estimation framework is formulated. In Section IV, the BMS hard thresholding technique is proposed, which is applied to the GraSP and GraHTP algorithms. In addition, an FFT-based fast implementation is proposed. In Section V, the results and discussion are presented, and the conclusions are followed in Section VI.
Notation: a, a, and A denote a scalar, vector, and matrix, respectively. a 0 , a 1 , and a represent the L 0 -, L 1 -, and L 2 -norm of a, respectively. A F is the Frobenius norm of A. The transpose, conjugate transpose, and conjugate of A are denoted as A T , A H , and A, respectively. The element-wise vector multiplication and division of a and b are denoted as a ⊙ b and a ⊘ b, respectively. The sum of all of the elements of a is denoted as sum(a). The vectorization of A is denoted as vec(A), which is formed by stacking all of the columns of A. The unvectorization of a is denoted as unvec(a), which is the inverse of vec(A). The Kronecker product of A and B is denoted as A ⊗ B. The support of a is denoted as supp(a), which is the set of indices 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 formed by leaving only the s largest (in absolute value) elements of a and replacing the other elements with 0. Similarly, the vector obtained by leaving only the elements of a indexed by the set A and replacing the other elements with 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|, respectively. The set difference between the sets A and B, namely A ∩ B c , is denoted as A \ B. φ(a) and Φ(a) are element-wise standard normal PDF and CDF functions of a, whose i-th elements are 1 2 dx, respectively. The m×1 zero vector and m×m identity matrix are denoted as 0 m and I m , respectively.

A. System Model
As shown in Fig. 1, consider a mmWave massive MIMO system with uniform linear arrays (ULAs) at the transmitter and receiver. The N -antenna transmitter transmits a training sequence of length T to the M -antenna receiver. Therefore, the received signal Y = y [1] which is the collection of the t-th received signal y[t] ∈ C M over t ∈ {1, . . . , T }. In the mmWave band, the channel H ∈ C M×N consists of a small number of paths, whose parameters are the path gains, angle-of-arrivals (AoAs), and angle-of-departures (AoDs) [21]. Therefore, H is where L is the number of paths, α ℓ ∈ C is the path gain of the ℓ-th path, and θ RX,ℓ ∈ [−π/2, π/2] and θ TX,ℓ ∈ [−π/2, π/2] are the AoA and AoD of the ℓ-th path, respectively. The steering vectors a RX (θ RX,ℓ ) ∈ C M and a TX (θ TX,ℓ ) ∈ C N are a TX (θ TX,ℓ ) = 1 √ N 1 · · · e −jπ(N −1) sin(θ TX,ℓ ) T where the inter-element spacing is half-wavelength. The train-  At the receiver, the real and imaginary parts of the received signal are quantized by one-bit ADCs. The quantized received signalŶ isŶ where Q(·) is the one-bit quantization function, whose threshold is zero. Therefore, Q(Y) is where sign(·) is the element-wise sign function. The goal is to estimate H by estimating {α ℓ } L ℓ=1 , {θ RX,ℓ } L ℓ=1 , and {θ TX,ℓ } L ℓ=1 fromŶ.

B. Virtual Channel Representation
In the mmWave channel model in (2), to Y renders a nonlinear channel estimation problem. To convert the nonlinear channel estimation problem to a linear channel estimation problem, we adopt the virtual channel representation [22].
The virtual channel representation of H is where the dictionary pair A RX ∈ C M×BRX and A TX ∈ C N ×BTX is the collection of B RX ≥ M steering vectors and B TX ≥ N steering vectors, respectively. Therefore, A RX and A TX are whose gridding AoAs {ω RX,i } BRX i=1 and AoDs {ω TX,j } BTX j=1 are selected so as to form overcomplete DFT matrices. The gridding AoAs and AoDs are the B RX and B TX points from [−π/2, π/2], respectively, to discretize the AoAs and AoDs because the ground truth AoAs and AoDs are unknown.
To make a dictionary pair of the overcomplete DFT matrix form, the gridding AoAs and AoDs are given as ω RX,i = sin −1 (−1 + 2/B RX · (i − 1)) and ω RX,j = sin −1 (−1 + 2/B TX · (j − 1)), respectively. We prefer overcomplete DFT matrices because they are relatively well-conditioned, and DFT matrices are friendly to the FFT-based implementation, which will be discussed in Section IV. The virtual channel is the nearest to (θ RX,ℓ , θ TX,ℓ ) but zero otherwise. In general, the error between H and A RX X * A H TX is inversely proportional to B RX and B TX . To approximate H using (7) with negligible error, the dictionary pair must be dense, namely B RX ≫ M and B TX ≫ N .
Before we proceed, we provide a supplementary explanation on the approximation in (7). In this work, we attempt to estimate the L-sparse X * in (7) because the sparse assumption on X * is favorable when the goal is to formulate the channel estimation problem as a sparsity-constrained problem. The cost of assuming that X * is L-sparse is, as evident, the approximation error shown in (7). Alternatively, the approximation error can be perfectly removed by considering X * satisfying H = A RX X * A H TX , i.e., equality instead of approximation. One well-known X * satisfying the equality is the minimum Frobenius norm solution, i.e., Such X * , however, has no evident structure to exploit in channel estimation, which is the reason why we assume that X * is L-sparse at the cost of the approximation error in (7).
In practice, the arrays at the transmitter and receiver are typically large to compensate the path loss in the mmWave band, whereas the number of line-of-sight (LOS) and near LOS paths is small [23]. Therefore, X * is sparse when the dictionary pair is dense because only L elements among B RX B TX elements are nonzero where L ≪ M N ≪ B RX B TX . In the sequel, we use the shorthand notation B = B RX B TX .
To facilitate the channel estimation framework, we vectorize (1) and (5) in conjunction with (7). First, note that where the gridding error E ∈ C M×T represents the mismatch in (7). 1 Then, the vectorized received signal 1 In practice, X * may be either approximately sparse or exactly sparse to formulate (10). If X * is approximately sparse, the leakage effect is taken into account so the mismatch in (7) becomes zero, namely vec(E) = 0 M T . In contrast, the mismatch in (7) must be taken into account with a nonzero E when X * is exactly sparse. Fortunately, E is inversely proportional to B RX and B TX . Therefore, we adopt the latter definition of X * and propose our algorithm ignoring E assuming that B RX ≫ M and B TX ≫ N . The performance degradation from E will become less as B RX and B TX become sufficiently large. where The vectorized quantized received signalŷ = vec(Ŷ) The goal is to estimate L-sparse x * fromŷ.

III. PROBLEM FORMULATION
In this section, we formulate the channel estimation problem using the MAP criterion. To facilitate the MAP channel estimation framework, the real counterparts of the complex forms in (16) are introduced. Then, the likelihood function of x * is derived.
The real counterparts are the collections of the real and imaginary parts of the complex forms. Therefore, the real which are the collections of the real and imaginary parts of y, A, and x * , respectively. In the sequel, we use the complex forms and the real counterparts interchangeably. For example, x * and x * R refer to the same entity. Before we formulate the likelihood function of x * , note that e is hard to analyze. However, e is negligible when the dictionary pair is dense. Therefore, we formulate the likelihood function of x * without e. The price of such oversimplification is negligible when B RX ≫ M and B TX ≫ N , which is to be shown in Section V where e = 0 MT . To derive the likelihood function of x * , note that given x * . Then, from (20) in conjunction with (16), the loglikelihood function f (x) is [10] f (x) = log Pr ŷ = Q( If the distribution of x * is known, the MAP estimate of x * is argmax where g MAP (x) is the logarithm of the PDF of x * . In practice, however, g MAP (x) is unknown. Therefore, we formulate the MAP channel estimation framework based on {α ℓ } L ℓ=1 , {θ RX,ℓ } L ℓ=1 , and {θ TX,ℓ } L ℓ=1 where we assume the followings: , and {θ TX,ℓ } L ℓ=1 are independent. Then, the MAP estimate of x * considering the channel sparsity is argmax ignoring the constant factor. However, note that only the optimization problems (22) and (23) are equivalent in the sense that their solutions are the same, not g MAP (x) and g(x). In the ML channel estimation framework, the ML estimate of x * is In the sequel, we focus on solving (23) because (23) reduces to (24) when g(x) = 0. In addition, we denote the objective function and the gradient in (23) as h(x) and ∇h(x), respectively. Therefore, where the differentiation is with respect to x.

IV. CHANNEL ESTIMATION VIA GRADIENT PURSUIT
In this section, we propose the BMSGraSP and BMS-GraHTP algorithms to solve (23), which are the variants of the GraSP [17] and GraHTP [18] algorithms, respectively. Then, an FFT-based fast implementation is proposed. In addition, we investigate the limit of the BMSGraSP and BMSGraHTP algorithms in the high SNR regime in one-bit ADCs.

A. Proposed BMSGraSP and BMSGraHTP Algorithms
Note that h(x) in (23) is concave because f (x) and g(x) are the sums of the logarithms of Φ(·) and φ(·), respectively, which are log-concave [24]. However, (23) is not a convex optimization problem because the sparsity constraint is not convex. Furthermore, solving (23) is NP-hard because of its combinatorial complexity. To approximately optimize convex objective functions with sparsity constraints iteratively by pursuing the gradient of the objective function, the GraSP and GraHTP algorithms were proposed in [17] and [18], respectively.
To solve (23), the GraSP and GraHTP algorithms roughly proceed as follows at each iteration when x is the current estimate of x * where the iteration index is omitted for simplicity.
First, the best L-term approximation of ∇h(x) is computed, which is where T L (·) is the L-term hard thresholding function. Here, T L (·) leaves only the L largest elements (in absolute value) of ∇h(x), and sets all the other remaining elements to 0. Then, after the estimate of supp(x * ) is updated by selecting i.e., I is the set of indices formed by collecting the L indices of ∇h(x) corresponding to its L largest elements (in absolute value), the estimate of x * is updated by solving the following optimization problem which can be solved using convex optimization because the support constraint is convex [24]. The GraSP and GraHTP algorithms are the generalizations of the CoSaMP [19] and HTP [20] algorithms, respectively. This follows because the gradient of the squared error is the scaled proxy of the residual.
To solve (23) using the GraSP and GraHTP algorithms, h(x) is required either to have a stable restricted Hessian [17] or to be strongly convex and smooth [18]. These conditions are the generalizations of the restricted isometry property (RIP) in CS [25], which means that h(x) is likely to satisfy these conditions when A is either a restricted isometry, wellconditioned, or incoherent. In practice, however, A is highly coherent because the dictionary pair is typically dense to reduce the mismatch in (7).
To illustrate how the GraSP and GraHTP algorithms fail to solve (23) when A is highly coherent, consider the real counterpart of ∇h(x). The real counterpart ∇h( is the inverse Mills ratio function. 2 Then, the following observation holds from directly computing ∇h(x i ), whose real and imaginary parts are the i-th and (i + B)-th elements of ∇h(x R ), respectively.
However, Observation 1 is meaningless because a i = a j unless i = j. To establish a meaningful observation, consider the coherence between a i and a j , which reflects the proximity between a i and a j according to [26], [27] 2 The element-wise vector division in the inverse Mills ratio function is meaningless because the arguments of the inverse Mills ratio function are scalars in (30). The reason we use the element-wise vector division in the inverse Mills ratio function will become clear in (37), whose arguments are vectors.  Then, using the η-coherence band, which is [26] where η ∈ (0, 1), we establish the following conjecture when η is sufficiently large.
At this point, we use Conjecture 1 to illustrate how the GraSP and GraHTP algorithms fail to estimate supp(x * ) from (28) by naive hard thresholding when A is highly coherent. To proceed, consider the following example, which assumes that x * andŷ are realized with x representing the current estimate of x * so as to satisfy 1) i = argmax is the by-product of i. Here, J is called the by-product of i because which follows from Conjecture 1, holds despite x * j = 0 for all j ∈ J . In other words, the by-product of i refers to the fact that ∇h(x i ) and ∇h(x j ) are indistinguishable for all j ∈ J according to (34), but the elements of x * indexed by J are 0 according to 3). Therefore, when we attempt to estimate supp(x * ) by hard thresholding ∇h(x), the indices in J will likely be erroneously selected as the estimate of supp(x * ) because ∇h(x j ) and the maximum element of ∇h(x), which is ∇h(x i ) according to 2), are indistinguishable for all j ∈ J .
To illustrate how (28) cannot estimate supp(x * ) when A is highly coherent, consider another example where ∇h(x) and T L (∇h(x)) are shown in Figs. 2 and 3, respectively. In this example, supp(x * ) is widely spread, whereas most of supp(T L (∇h(x))) are in the coherence band of the index of the maximum element of ∇h(x). This shows that hard thresholding ∇h(x) is not sufficient to distinguish whether an index is the ground truth index or the by-product of another index. To solve this problem, we propose the BMS hard thresholding technique.
The BMS hard thresholding function T BMS,L (·) is an Lterm hard thresholding function, which is proposed based on Conjecture 1. The BMS hard thresholding technique is presented in Algorithm 1. Line 3 selects the index of the maximum element of ∇h(x) from the unchecked index set as the current index. Line 4 constructs the by-product testing set. Line 5 checks whether the current index is greater than the by-product testing set. In this paper, Line 5 is referred to as 3 We use the term "ground truth" to emphasize that the ground truth x * is the true virtual channel which actually gives the quantized received signalŷ from (16), whereas x merely represents the point where ∇h(x) is computed to estimate supp(x * ) via hard thresholding. the band maximum criterion. If the band maximum criterion is satisfied, the current index is selected as the estimate of supp(x * ) in Line 6. Otherwise, the current index is not selected as the estimate of supp(x * ) because the current index is likely to be the by-product of another index rather than the ground truth index. Line 8 updates the unchecked index set.
Note that Algorithm 1 is a hard thresholding technique applied to ∇h(x). If the BMS hard thresholding technique is applied to x + κ∇h(x) where κ is the step size, ∇h(x) is replaced by x + κ∇h(x) in the input, output, and Lines 3, 5, and 10 of Algorithm 1. This can be derived using the same logic based on Conjecture 1. Now, we propose the BMSGraSP and BMSGraHTP algorithms to solve (23).
The BMSGraSP and BMSGraHTP algorithms are the variants of the GraSP and GraHTP algorithms, respectively. The difference between the BMS-based and non-BMS-based algorithms is that the hard thresholding function is T BMS,L (·) instead of T L (·). The BMSGraSP and BMSGraHTP algorithms are presented in Algorithms 2 and 3, respectively. Lines 3, 4, and 5 of Algorithms 2 and 3 roughly proceed based on the same logic. Line 3 computes the gradient of the objective function. Line 4 selects I from the support of the hard thresholded gradient of the objective function. Line 5 maximizes the objective function subject to the support constraint. This can be solved using convex optimization because the objective function and support constraint are concave and convex, respectively. In addition, b is hard thresholded in Line 6 of Algorithm 2 because b is at most 3L-sparse. A natural halting condition of Algorithms 2 and 3 is to halt when the current and previous supp(x) are the same. The readers who are interested in a more in-depth analyses of the GraSP and GraHTP algorithms are referred to [17] and [18], respectively.
Remark 1: Instead of hard thresholding b, we can solvẽ which is a convex optimization problem, to obtainx in Line 6 of Algorithm 2. This is the debiasing variant of Algorithm 2 [17]. The advantage of the debiasing variant of Algorithm 2 is a more accurate estimate of x * . However, the complexity is increased, which is incurred by solving (35).

Remark 2:
In this paper, we assume that only h(x) and ∇h(x) are required at each iteration to solve (23) using Algorithms 2 and 3, which can be accomplished when the first order method is used to solve convex optimization problems in Line 5 of Algorithms 2 and 3. An example of such first order method is the gradient descent method with the backtracking line search [24].

B. Fast Implementation via FFT
In practice, the complexity of Algorithms 2 and 3 is demanding because h(x) and ∇h(x) are required at each iteration, which are high-dimensional functions defined on C B where B ≫ M N . In recent works on channel estimation and data detection in the mmWave band [14], [15], [28], the FFT-based implementation is widely used because H can be approximated by (7) using overcomplete DFT matrices. In this paper, an FFT-based fast implementation of h(x) and ∇h(x) is proposed, which is motivated by [14], [15], [28].
To facilitate the analysis, we convert the summations in h(x) and ∇h(x R ) to matrix-vector multiplications by algebraically manipulating (21) and (30). Then, we obtain where we see that the bottlenecks of h(x) and ∇h(x) come from the matrix-vector multiplications involving A R and A T R resulting from the large size of A. For example, the size of A is 5120 × 65536 in Section V where M = N = 64, B RX = B TX = 256, and T = 80. (37) with c ∈ C MT being the complex form of c R . From the fact that

To develop an FFT-based fast implementation of the matrixvector multiplications involving
we now attempt to compute Ax and A H c via the FFT. Then, Ax and A H c are unvectorized according to where X = unvec(x) ∈ C BRX×BTX and C = unvec(c) ∈ C M×T . If the matrix multiplication involving S can be implemented using the FFT, e.g., Zadoff-Chu (ZC) [29] or DFT [11] training sequence, (40) and (41) can be implemented using the FFT because A RX and A TX are overcomplete DFT matrices. For example, each column of A TX X H in (40) can be computed using the B TX -point FFT with pruned outputs, i.e., retaining only N outputs, because we constructed A TX as an overcomplete DFT matrix.
In particular, the matrix multiplications involving A TX , S H , and A RX in (40)  Therefore, the complexity of Algorithms 2 and 3 is reduced when h(x) and ∇h(x) are implemented using the FFT operations.
Remark 3: Line 5 of Algorithms 2 and 3 is equivalent to solving where f I (x I ) = log Pr ŷ = Q( and A I ∈ C MT ×|I| is the collection of a i with i ∈ I. Therefore, only h I (x I ) and ∇h I (x I ) are required in Line 5 of Algorithms 2 and 3, which are low-dimensional functions defined on C |I| where |I| = O(L). If h I (x I ) and ∇h I (x I ) are computed based on the same logic in (40) and (41) but A replaced by A I , the complexity of Algorithms 2 and 3 is reduced further because the size of the FFT is reduced in Line 5.

V. RESULTS AND DISCUSSION
In this section, we evaluate the performance of Algorithms 2 and 3 from different aspects in terms of the accuracy, achievable rate, and complexity. Throughout this section, we consider a mmWave massive MIMO system with one-bit ADCs, whose parameters are M = N = 64 and T = 80. The rest vary from simulation to simulation, which consist of B RX , B TX , and L. In addition, we consider S, whose rows are the circular shifts of the ZC training sequence of length T as in [15], [33]. Furthermore, H is either random or deterministic. If H is random, α ℓ ∼ CN (0, 1), θ RX,ℓ ∼ unif([−π/2, π/2]), and θ TX,ℓ ∼ unif([−π/2, π/2]) are independent. Otherwise, we consider different H from simulation to simulation.
The MSEs of {α ℓ } L ℓ=1 , {θ RX,ℓ } L ℓ=1 , and {θ TX,ℓ } L ℓ=1 are where (α ℓ ,θ RX,ℓ ,θ TX,ℓ ) corresponds to some nonzero element ofX = unvec(x) ∈ C BRX×BTX . The normalized MSE (NMSE) of H is whereH = A RXX A TX . In (45)-(48), the symbol˜is used to emphasize that the quantity is an estimate. Throughout this section, we consider the debiasing variant of Algorithm 2. The halting condition of Algorithms 2 and 3 is to halt when the current and previous supp(x) are the same. The gradient descent method is used to solve convex optimization problems, which consist of (35) and Line 5 of Algorithms 2 and 3. The backtracking line search is used to compute the step size in the gradient descent method and κ in Line 3 of Algorithm 3. In addition, η is selected so that Conjecture 1 is satisfied. In this paper, we select the maximum η satisfying min i∈{1,2,...,B} The BE hard thresholding technique was proposed in [26], which was applied to the orthogonal matching pursuit (OMP) algorithm [34]. In this paper, we apply the BE hard thresholding technique to the GraSP algorithm, which results in the BEGraSP algorithm. In this example, B RX = B TX = 256 for the BMS-based and BE-based algorithms. We assume that L = 8 and H is deterministic where α ℓ = (0.8 + 0.1(ℓ − 1))e j π 4 (ℓ−1) . However, {θ RX,ℓ } L ℓ=1 and {θ TX,ℓ } L ℓ=1 vary from simulation to simulation, which are either widely spread (Fig. 4) or closely spread (Fig. 5). In Figs. 4 and 5, the notion of widely and closely spread paths refer to the fact that the minimum 2-norm distance between the paths are either relatively far or close, i.e., min i =j (θ RX,i −θ RX,j , θ TX,i −θ TX,j ) 2 of Fig. 4, which is (π/18) 2 + (π/18) 2 , is greater than that of Fig. 5, which is (π/36) 2 + (π/36) 2 . The path gains, AoAs, and AoDs are assumed to be deterministic because the CRB is defined for deterministic parameters only [35]. A variant of the CRB for random parameters is known as the Bayesian CRB, but adding the Bayesian CRB to our work is left as a future work because applying the Bayesian CRB to nonlinear measurements, e.g., one-bit ADCs, is not as straightforward.
According to Figs. 4 and 5, the BMS-based algorithms succeed to estimate both widely spread and closely spread paths, whereas the BE-based algorithms fail to estimate closely spread paths. This follows because the BE hard thresholding technique was derived based on the assumption that supp(x * ) is widely spread. In contrast, the BMS hard thresholding technique is proposed based on Conjecture 1 without any assumption on supp(x * ). This means that when supp(x * ) is closely spread, the BE hard thresholding technique cannot properly estimate supp(x * ) because the BE hard thresholding technique, by its nature, excludes the elements near the maximum element of x * from its potential candidate. The BMS hard thresholding technique, in contrast, uses the elements near the maximum element of x * to construct the by-product testing set only, i.e., Line 4 of Algorithm 1. Therefore, the BMS-based algorithms are superior to the BE-based algorithms when the paths are closely spread. The Cramér-Rao bounds (CRBs) of MSE({α ℓ } L ℓ=1 ), MSE({θ RX,ℓ } L ℓ=1 ), and MSE({θ TX,ℓ } L ℓ=1 ) are provided, which were derived in [36]. The gaps between the MSEs and their corresponding CRBs can be interpreted as a performance limit incurred by the discretized AoAs and AoDs. To overcome such limit, the AoAs and AoDs must be estimated based on the off-grid method, which is beyond the scope of this paper. . The CRB is provided as a reference, which was derived in [36].
In addition, note that MSE({α ℓ } L ℓ=1 ) and NMSE(H) worsen as the SNR enters the high SNR regime. To illustrate why x * cannot be estimated in the high SNR regime in one-bit ADCs, note that in the high SNR regime with c > 0, which means that x * and cx * are indistinguishable because the magnitude information is lost by one-bit ADCs. The degradation of the recovery accuracy in the high SNR regime with one-bit ADCs is an inevitable phenomenon, as observed from other previous works on low-resolution ADCs [11], [14], [15], [33], [37]. In Figs. 6 and 7, we compare the performance of Algorithms 2, 3, and other estimators when H is random. The Bernoulli Gaussian-GAMP (BG-GAMP) algorithm [15] is an iterative approximate MMSE estimator of x * , which was derived based on the assumption that x * i is distributed as CN (0, 1) with probability L/B but zero otherwise, namely the BG distribution. The fast iterative shrinkage-thresholding algorithm (FISTA) [38] is an iterative MAP estimator of x * , which was derived based on the assumption that the logarithm of the PDF of x * is g FISTA (x) = −γ x 1 ignoring the constant factor, namely the Laplace distribution. Therefore, the estimate of x * is which is solved using the accelerated proximal gradient descent method [38]. The regularization parameter γ is selected so that the expected sparsity of (51) is 3L for a fair comparison, which was suggested in [17]. In this example, L = 4, whereas B RX and B TX vary from algorithm to algorithm. In particular, we select B RX = B TX = 256 for Algorithms 2, 3, and the FISTA, whereas B RX = M and B TX = N for the BG-GAMP algorithm. In Fig. 6, we compare the accuracy of Algorithms 2, 3, and other estimators at different SNRs using NMSE(H). . The CRB is provided as a reference, which was derived in [36].
According to Fig. 6, Algorithms 2 and 3 outperform the BG-GAMP algorithm and FISTA as the SNR enters the medium SNR regime. The accuracy of the BG-GAMP algorithm is disappointing because the mismatch in (7) is inversely proportional to B RX and B TX . However, increasing B RX and B TX is forbidden because the BG-GAMP algorithm diverges when A is highly coherent. The accuracy of the FISTA is disappointing because the Laplace distribution does not match the distribution of x * . Note that (23), which is the basis of Algorithms 2 and 3, is indeed the MAP estimate of x * , which is in contrast to the FISTA. According to Fig. 6, NMSE(H) worsens as the SNR enters the high SNR regime, which follows from the same reason as in Figs. 4 and 5.
In Fig. 7, we compare the achievable rate lower bound of Algorithms 2, 3, and other estimators at different SNRs when the precoders and combiners are selected based onH. The achievable rate lower bound shown in Fig. 7 is presented in [15], which was derived based on the Bussgang decomposition [13] in conjunction with the fact that the worst-case noise is Gaussian. According to Fig. 7, Algorithms 2 and 3 outperform the BG-GAMP algorithm and FISTA, which is consistent with the result in Fig. 6.
In Fig. 8, we compare the complexity of Algorithms 2, 3, and other estimators at different B RX and B TX when H is random. To analyze the complexity, note that Algorithms 2, 3, and the FISTA require h(x) and ∇h(x) at each iteration, whose bottlenecks are Ax and A H c, respectively, while the BG-GAMP algorithm requires Ax and A H c at each iteration. Therefore, the complexity is measured based on the number of complex multiplications performed to compute Ax and A H c, which are implemented based on the FFT. In this example, L = 4, whereas SNR is either 0 dB or 10 dB.
In this paper, the complexity of the BG-GAMP algorithm is used as a baseline because the BG-GAMP algorithm is widely used. The normalized complexity is defined as the number of complex multiplications performed divided by the per-iteration  complexity of the BG-GAMP. For example, the normalized complexity of the FISTA with B RX = B TX = 256 is 160 when the complexity of the FISTA with B RX = B TX = 256 is equivalent to the complexity of the 160-iteration BG-GAMP algorithm with B RX = B TX = 256. In practice, the BG-GAMP algorithm converges in 15 iterations when A is incoherent [39]. In this paper, an algorithm is said to be as efficient as the BG-GAMP algorithm when the normalized complexity is below the target threshold, which is 15. As a sidenote, our algorithms, namely the BMSGraSP and BMSGraHTP, requires 2.1710 and 2.0043 iterations in average, respectively, across the entire SNR range.
According to Fig. 8, the complexity of the FISTA is impractical because the objective function of (51) is a highdimensional function defined on C B where B ≫ M N . In contrast, the complexity of Algorithms 2 and 3 is dominated by (42), whose objective function is a low-dimensional function defined on C |I| where |I| = O(L). The normalized com-

VI. CONCLUSIONS
In the mmWave band, the channel estimation problem is converted to a sparsity-constrained optimization problem, which is NP-hard to solve. To approximately solve sparsityconstrained optimization problems, the GraSP and GraHTP algorithms were proposed in CS, which pursue the gradient of the objective function. The GraSP and GraHTP algorithms, however, break down when the objective function is ill-conditioned, which is incurred by the highly coherent sensing matrix. To remedy such break down, we proposed the BMS hard thresholding technique, which is applied to the GraSP and GraHTP algorithms, namely the BMSGraSP and BMSGraHTP algorithms, respectively. Instead of directly hard thresholding the gradient of the objective function, the BMSbased algorithms test whether an index is the ground truth index or the by-product of another index. We also proposed an FFT-based fast implementation of the BMS-based algorithms, whose complexity is reduced from O(M 4 ) to O(M 2 log M ). In the simulation, we compared the performance of the BMSbased, BE-based, BG-GAMP, and FISTA algorithms from different aspects in terms of the accuracy, achievable rate, and complexity. The BMS-based algorithms were shown to outperform other estimators, which proved to be both accurate and efficient. Our algorithms, however, addressed only the flat fading scenario, so an interesting future work would be to propose a low-complexity channel estimator capable of dealing with the wideband scenario.