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

TL;DR

This paper introduces band maximum selecting (BMS) techniques to enhance gradient pursuit algorithms for sparse channel estimation in mmWave massive MIMO systems with one-bit ADCs, addressing ill-conditioning issues.

Contribution

It proposes BMSGraSP and BMSGraHTP algorithms that improve gradient pursuit-based channel estimation under challenging conditions in mmWave MIMO systems.

Findings

BMSGraSP and BMSGraHTP outperform existing methods in accuracy.

The algorithms effectively handle ill-conditioned problems.

Enhanced channel estimation accuracy in one-bit quantized mmWave systems.

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.

Equations

$$\mathbf{Y}=\sqrt{\rho }\mathbf{HS}+\mathbf{N}$$

$$\mathbf{H}=\sum _{\ell =1}^L{\alpha }_{\ell }{\mathbf{a}}_{\mathrm{RX}}({\theta }_{\mathrm{RX},\ell }){\mathbf{a}}_{\mathrm{TX}}({\theta }_{\mathrm{TX},\ell }{)}^{\mathrm{H}}$$

$${\mathbf{a}}_{\mathrm{RX}}({\theta }_{\mathrm{RX},\ell })$$

$${\mathbf{a}}_{\mathrm{TX}}({\theta }_{\mathrm{TX},\ell })$$

$$\hat{\mathbf{Y}}$$

$$=\mathrm{Q}(\sqrt{\rho }\mathbf{HS}+\mathbf{N})$$

$$\mathrm{Q}(\mathbf{Y})=\mathrm{sign}(\mathrm{Re}(\mathbf{Y}))+j\mathrm{sign}(\mathrm{Im}(\mathbf{Y}))$$

$$\mathbf{H}={\mathbf{A}}_{\mathrm{RX}}{\mathbf{X}}^{\ast }{\mathbf{A}}_{\mathrm{TX}}^{\mathrm{H}}$$

$$\hat{\mathbf{y}}$$

$$=\mathrm{Q}(\sqrt{\rho }\mathbf{A}{\mathbf{x}}^{\ast }+\mathbf{n})$$

$${\hat{\mathbf{y}}}_{\mathrm{R}}$$

$$={\left[\begin{array}{ccc}{\hat{y}}_{\mathrm{R},1}& \cdots & {\hat{y}}_{\mathrm{R},2MT}\end{array}\right]}^{\mathrm{T}},$$

$${\mathbf{A}}_{\mathrm{R}}$$

$$={\left[\begin{array}{ccc}{\mathbf{a}}_{\mathrm{R},1}& \cdots & {\mathbf{a}}_{\mathrm{R},2MT}\end{array}\right]}^{\mathrm{T}},$$

$${\mathbf{x}}_{\mathrm{R}}^{\ast }$$

$$f(\mathbf{x})$$

$$=\sum _{i=1}^{2MT}\log \Phi (\sqrt{2\rho }{\hat{y}}_{\mathrm{R},i}{\mathbf{a}}_{\mathrm{R},i}^{\mathrm{T}}{\mathbf{x}}_{\mathrm{R}}).$$

$$\underset{\mathbf{x}\in {\mathbb{C}}^B}{\mathrm{argmax}}(f(\mathbf{x})+g(\mathbf{x}))\text{s.t.}\Vert \mathbf{x}{\Vert }_0\le L$$

$$h(\mathbf{x})$$

$$\nabla h(\mathbf{x})$$

$$={\left[\begin{array}{ccc}\nabla h(x_1)& \cdots & \nabla h(x_B)\end{array}\right]}^{\mathrm{T}},$$

$$\nabla h(\hat{\mathbf{x}}){|}_L,$$

$$\mathcal{I}=\mathrm{supp}(\nabla h(\hat{\mathbf{x}}){|}_L)$$

$$\underset{\mathbf{x}\in {\mathbb{C}}^B}{\mathrm{argmax}}h(\mathbf{x})\text{s.t.}\mathrm{supp}(\mathbf{x})\subseteq \mathcal{I},$$

$$\nabla h({\mathbf{x}}_{\mathrm{R}})$$

$$=$$

$$\mu (i,j)=\frac{|{\mathbf{a}}_i^{\mathrm{H}}{\mathbf{a}}_j|}{\Vert {\mathbf{a}}_i\Vert \Vert {\mathbf{a}}_j\Vert },$$

$$B_{\eta }(i)=\{j\mid \mu (i,j)\ge \eta \}$$

$${\mathcal{J}}_{\eta }(i)=\{j\mid j\in B_{\eta }(i),{\hat{x}}_i={\hat{x}}_j\}\setminus \{i\}$$

$$|\nabla h({\hat{x}}_j)|$$

$$\stackrel{(b)}{=}\underset{k\in \{1,\dots ,B\}}{\max }|\nabla h({\hat{x}}_k)|$$

$$\hat{\mathbf{x}}=\underset{\mathbf{x}\in {\mathbb{C}}^B}{\mathrm{argmax}}h(\mathbf{x})\text{s.t.}\mathrm{supp}(\mathbf{x})\subseteq \mathrm{supp}(\mathbf{b}{|}_L)$$

$$\mathrm{unvec}(\mathbf{Ax})$$

$$={\mathbf{A}}_{\mathrm{RX}}({\mathbf{S}}^{\mathrm{H}}({\mathbf{A}}_{\mathrm{TX}}{\mathbf{X}}^{\mathrm{H}}){)}^{\mathrm{H}},$$

$$\mathrm{unvec}({\mathbf{A}}^{\mathrm{H}}\mathbf{c})$$

$$={\mathbf{A}}_{\mathrm{RX}}^{\mathrm{H}}({\mathbf{A}}_{\mathrm{TX}}^{\mathrm{H}}(\mathbf{S}{\mathbf{C}}^{\mathrm{H}}){)}^{\mathrm{H}}$$

$$\mathrm{NMSE}=\mathbb{E}\left\{\frac{\Vert \hat{\mathbf{H}}-\mathbf{H}{\Vert }_{\mathrm{F}}^2}{\Vert \mathbf{H}{\Vert }_{\mathrm{F}}^2}\right\}$$

Full text

I Introduction

The wide bandwidth of the millimeter wave (mmWave) band provides high data rates, resulting in a significant performance gain [1, 2, 3, 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 high-resolution 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, 7, 8, 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, 11, 12], compressed sensing-based (CS-based) signal recovery techniques for mmWave massive MIMO systems with low-resolution ADCs were proposed. In [10], a sparse Bayesian learning-based (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.

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$, $\mathbf{a}$, and $\mathbf{A}$ denote a scalar, vector, and matrix. The transpose, conjugate transpose, and conjugate of $\mathbf{A}$ are denoted as $\mathbf{A}^{\mathrm{T}}$, $\mathbf{A}^{\mathrm{H}}$, and $\overline{\mathbf{A}}$. The Kronecker product of $\mathbf{A}$ and $\mathbf{B}$ is denoted as $\mathbf{A}\otimes\mathbf{B}$. The support of $\mathbf{a}$ is denoted as $\mathrm{supp}(\mathbf{a})$, which is formed by collecting all of the indices of the nonzero elements of $\mathbf{a}$. The best $s$-term approximation of $\mathbf{a}$ is denoted as $\mathbf{a}|_{s}$, which is obtained by leaving only the $s$ largest elements of $\mathbf{a}$ and hard thresholding other elements to [math]. For a set $\mathcal{A}$, the vector obtained by leaving only the elements of $\mathbf{a}$ indexed by $\mathcal{A}$ and hard thresholding the remaining elements to [math] is denoted as $\mathbf{a}|_{\mathcal{A}}$. The absolute value of a scalar $a$ and cardinality of a set $\mathcal{A}$ are denoted as $|a|$ and $|\mathcal{A}|$. The standard normal PDF and CDF are denoted as $\phi(x)=\frac{1}{\sqrt{2\pi}}e^{-\frac{x^{2}}{2}}$ and $\Phi(x)=\int_{-\infty}^{x}\frac{1}{\sqrt{2\pi}}e^{-\frac{y^{2}}{2}}dy$. The inverse Mills ratio function is defined as $\lambda(x)=\frac{\phi(x)}{\Phi(x)}$. The element-wise matrix multiplication and division are denoted as $\odot$ and $\oslash$. The element-wise standard normal PDF, CDF, and inverse Mills ratio function are $\phi(\mathbf{x})$, $\Phi(\mathbf{x})$, and $\lambda(\mathbf{x})=\phi(\mathbf{x})\oslash\Phi(\mathbf{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 $\mathbf{Y}\in\mathbb{C}^{M\times T}$ is

[TABLE]

where $\mathbf{H}\in\mathbb{C}^{M\times N}$ is the channel, $\mathbf{S}\in\mathbb{C}^{N\times T}$ is the training sequence, whose columns obey the $2$-norm constraint of $\sqrt{N}$, and $\mathbf{N}\in\mathbb{C}^{M\times T}$ is the additive white Gaussian noise (AWGN), which is distributed as $\mathrm{vec}(\mathbf{N})\sim\mathcal{CN}(\mathbf{0}_{MT},\mathbf{I}_{MT})$. The signal-to-noise ratio (SNR) is defined as $\rho$. In the mmWave band, $\mathbf{H}$ contains a small number of paths where each path is associated with its path gain, angle-of-arrival (AoA), and angle-of-departure (AoD) [15]. Then, $\mathbf{H}$ is

[TABLE]

where $L$ is the number of paths, $\alpha_{\ell}\sim\mathcal{CN}(0,1)$ is the $\ell$-th path gain, and $\theta_{\mathrm{RX},\ell}\sim\mathrm{unif}([-\pi/2,\pi/2])$ and $\theta_{\mathrm{TX},\ell}\sim\mathrm{unif}([-\pi/2,\pi/2])$ are the $\ell$-th AoA and AoD, which are independent. The steering vectors $\mathbf{a}_{\mathrm{RX}}(\theta_{\mathrm{RX},\ell})\in\mathbb{C}^{M}$ and $\mathbf{a}_{\mathrm{TX}}(\theta_{\mathrm{TX},\ell})\in\mathbb{C}^{N}$ are

[TABLE]

whose inter-element spacings are half-wavelength. The one-bit quantized received signal $\hat{\mathbf{Y}}\in\mathbb{C}^{M\times T}$ is

[TABLE]

where $\mathrm{Q}(\cdot)$ is the zero threshold one-bit quantization function defined as

[TABLE]

with $\mathrm{sign}(\cdot)$ the element-wise sign function.

The virtual channel representation [16] of $\mathbf{H}$ is

[TABLE]

where $\mathbf{A}_{\mathrm{RX}}\in\mathbb{C}^{M\times B_{\mathrm{RX}}}$ and $\mathbf{A}_{\mathrm{TX}}\in\mathbb{C}^{N\times B_{\mathrm{TX}}}$ are overcomplete discrete Fourier transform (DFT) matrices, and $\mathbf{X}^{*}\in\mathbb{C}^{B_{\mathrm{RX}}\times B_{\mathrm{TX}}}$ is the virtual channel with $B_{\mathrm{RX}}\geq M$ and $B_{\mathrm{TX}}\geq N$. For notational simplicity, we define $B=B_{\mathrm{RX}}B_{\mathrm{TX}}$. To facilitate the analysis, (II) is vectorized after plugging in (7) as

[TABLE]

where $\hat{\mathbf{y}}=\mathrm{vec}(\hat{\mathbf{Y}})$, $\mathbf{y}=\mathrm{vec}(\mathbf{Y})$, $\mathbf{A}=\mathbf{S}^{\mathrm{T}}\overline{\mathbf{A}}_{\mathrm{TX}}\otimes\mathbf{A}_{\mathrm{RX}}=\begin{bmatrix}\mathbf{a}_{1}&\cdots&\mathbf{a}_{B}\end{bmatrix}$, $\mathbf{x}^{*}=\mathrm{vec}(\mathbf{X}^{*})=\begin{bmatrix}x_{1}^{*}&\cdots&x_{B}^{*}\end{bmatrix}^{\mathrm{T}}$, and $\mathbf{n}=\mathrm{vec}(\mathbf{N})$.

To build the MAP channel estimation framework, we start with formulating the likelihood function based on the real forms of $\hat{\mathbf{y}}$, $\mathbf{A}$, and $\mathbf{x}^{*}$, which are

[TABLE]

where the complex forms and their real forms are used interchangeably in the sequel. For example, $\mathbf{x}^{*}$ and $\mathbf{x}_{\mathrm{R}}^{*}$ represent the same entity. Then, the log-likelihood function $f(\mathbf{x})$ is [17]

[TABLE]

However, $\mathbf{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 $\mathbf{x}^{*}$ is $L$-sparse with independent and identically distributed (i.i.d.) $\mathcal{CN}(0,1)$ elements, whose locations are uniformly distributed in $\{1,\dots,B\}$. The mismatch between our assumption on $\mathbf{x}^{*}$ and its true distribution reduces as we increase $B_{\mathrm{RX}}$ and $B_{\mathrm{TX}}$ due to the reduced leakage effect. In Section IV, it is shown that sufficiently large $B_{\mathrm{RX}}$ and $B_{\mathrm{TX}}$ lead to accurate channel estimation because our assumption closely approximates the true distribution of $\mathbf{x}^{*}$. With such assumption on the distribution of $\mathbf{x}^{*}$, the MAP estimate of $\mathbf{x}^{*}$ is

[TABLE]

where $g(\mathbf{x})=-\|\mathbf{x}_{\mathrm{R}}\|^{2}$ is the logarithm of the PDF of $\mathcal{CN}(\mathbf{0}_{B},\mathbf{I}_{B})$ without the constant factor. In the sequel, the objective function and its gradient in (13) are denoted as

[TABLE]

whose differentiation is with respect to $\mathbf{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 $\hat{\mathbf{x}}$ of $\mathbf{x}^{*}$ roughly as follows at each iteration. The best $L$-term approximation of $\nabla h(\hat{\mathbf{x}})$ is computed via hard thresholding as

[TABLE]

whose support

[TABLE]

is selected as the updated estimate of $\mathrm{supp}(\mathbf{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 $\hat{\mathbf{x}}$ as

[TABLE]

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

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

To illustrate how support identification fails from (17) when $\mathbf{A}$ is highly coherent, first, consider the real form $\nabla h(\mathbf{x_{\mathrm{R}}})$ of $\nabla h(\mathbf{x})$ defined as

[TABLE]

which is obtained from $\nabla\log\Phi(\mathbf{a}_{\mathrm{R}}^{\mathrm{T}}\mathbf{x}_{\mathrm{R}})=\lambda(\mathbf{a}_{\mathrm{R}}^{\mathrm{T}}\mathbf{x}_{\mathrm{R}})\mathbf{a}_{\mathrm{R}}$ and $\nabla\|\mathbf{x}_{\mathrm{R}}\|^{2}=2\mathbf{x}_{\mathrm{R}}$. Then, we establish the following observation, which can be checked from directly computing $\nabla h(x_{i})$, whose real and imaginary parts correspond to the $i$-th and $(i+B)$-th elements of $\nabla h(\mathbf{x}_{\mathrm{R}})$.

Observation 1**.**

$\nabla h(x_{i})=\nabla h(x_{j})$* if $\mathbf{a}_{i}=\mathbf{a}_{j}$ and $x_{i}=x_{j}$.*

However, Observation 1 is trivial since $\mathbf{a}_{i}\neq\mathbf{a}_{j}$ unless $i=j$. To introduce a nontrivial observation, we adopt the notion of coherence between $\mathbf{a}_{i}$ and $\mathbf{a}_{j}$, i.e.,

[TABLE]

which measures the proximity between $\mathbf{a}_{i}$ and $\mathbf{a}_{j}$ [22, 23, 24]. Then, with the $\eta$-coherence band of $i$, which is defined as [22]

[TABLE]

with $\eta\in(0,1)$, the following conjecture is established for sufficiently large $\eta$.

Conjecture 1**.**

$\nabla h(x_{i})\approx\nabla h(x_{j})$* if $j\in B_{\eta}(i)$ and $x_{i}=x_{j}$.*

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

$i=\underset{k\in\{1,\dots,B\}}{\mathrm{argmax}}\ |\nabla h(\hat{x}_{k})|$ 2. 2.

$\mathcal{J}_{\eta}(i)\cap\mathrm{supp}(\mathbf{x}^{*})=\emptyset$

where

[TABLE]

is defined as the by-product of $i$. In this example, $\mathcal{J}_{\eta}(i)$ is called the by-product of $i$ because for all $j\in\mathcal{J}_{\eta}(i)$,

[TABLE]

holds where (a) follows from Conjecture 1 and (b) from 1) even if $\mathcal{J}_{\eta}(i)\cap\mathrm{supp}(\mathbf{x}^{*})=\emptyset$ according to 2). The implication of this example is that most indices of $\mathcal{J}_{\eta}(i)$ are likely to be selected when $\nabla h(\hat{\mathbf{x}})$ is hard thresholded as in (17) due to (III), which results in an erroneous estimate of $\mathrm{supp}(\mathbf{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 $\nabla h(\hat{\mathbf{x}})$ in (17) in order to exclude the by-product indices. The proposed BMS hard thresholding function $T_{\mathrm{BMS},L}(\cdot)$ 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 $\nabla h(\hat{\mathbf{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 $\mathrm{supp}(\mathbf{x}^{*})$ in Line 6. Otherwise, the current index is excluded because it is likely to be the by-product of 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 $\nabla h(\hat{\mathbf{x}})$. To apply the BMS hard thresholding function to $\hat{\mathbf{x}}+\kappa\nabla h(\hat{\mathbf{x}})$ where $\kappa$ is the step size, simply replace $\nabla h(\hat{\mathbf{x}})$ with $\hat{\mathbf{x}}+\kappa\nabla h(\hat{\mathbf{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 non-BMS-based algorithms is that $T_{\mathrm{BMS},L}(\cdot)$ 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, $\mathcal{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 $\mathbf{b}$ since the sparsity of $\mathbf{b}$ is at most $3L$. A natural halting condition for Algorithms 2 and 3 is to halt when $\mathrm{supp}(\hat{\mathbf{x}})$ does not change from iteration to iteration [13, 14].

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

[TABLE]

to update $\hat{\mathbf{x}}$. This is called the debiasing variant of Algorithm 2, which produces a more accurate $\mathbf{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 complexity because the support of their optimization variables are constrained to $\mathcal{I}$ where $|\mathcal{I}|=O(L)$ is typically much smaller than $B\geq 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 $\mathbb{C}^{B}$. We can reduce the complexity in Line 3 of Algorithms 2 and 3 based on the fast Fourier transform (FFT) implementation of $\nabla h(\hat{\mathbf{x}})$ when $\mathbf{S}$ has an FFT-friendly structure, e.g., circularly shifted Zadoff-Chu (ZC) sequences or partial DFT matrix. From (III), note that the matrix-vector multiplications involving $\mathbf{A}$ and $\mathbf{A}^{\mathrm{H}}$ act as computational bottlenecks. These matrix-vector multiplications can be efficiently implemented using the FFT by noting that

[TABLE]

where $\mathbf{C}=\mathrm{unvec}(\mathbf{c})$. By performing matrix multiplications involving $\mathbf{A}_{\mathrm{RX}}$, $\mathbf{A}_{\mathrm{TX}}$, and $\mathbf{S}$ using the FFT, the cost of computing $\nabla h(\hat{\mathbf{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 one-bit ADCs where $M=N=64$ and $T=80$. The rows of the training signal $\mathbf{S}$ are chosen as circularly shifted ZC sequences of length $T$ [25], and the number of paths are $L=4$. We choose $B_{\mathrm{RX}}=256\gg M$ and $B_{\mathrm{TX}}=256\gg N$ for the proposed BMSGraSP and BMSGraHTP to reduce the leakage effect, but other channel estimators may use smaller $B_{\mathrm{RX}}$ and $B_{\mathrm{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

[TABLE]

where $\hat{\mathbf{H}}=\mathbf{A}_{\mathrm{RX}}\hat{\mathbf{X}}\mathbf{A}_{\mathrm{TX}}$ with $\hat{\mathbf{X}}=\mathrm{unvec}(\hat{\mathbf{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 $\mathrm{supp}(\hat{\mathbf{x}})$ does not change from iteration to iteration. We use the backtracking line search [21] to compute $\kappa$ in Line 3 of Algorithm 3. Lastly, we configure $\eta$ so as to satisfy Conjecture 1 by selecting the maximum $\eta$ satisfying

[TABLE]

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 (BG-GAMP) [12], and fast iterative shrinkage-thresholding algorithm (FISTA) [26]. BG-GAMP is a loopy BP-based iterative approximate MMSE channel estimator, which assumes that $\mathbf{x}^{*}$ is distributed as i.i.d. BG, meaning that each elements are $\mathcal{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 $\mathbf{x}^{*}$ is $g_{\mathrm{FISTA}}(\mathbf{x})=-\gamma\|\mathbf{x}\|_{1}$ without the constant factor, which is the Laplace distribution. Then, the FISTA estimate of $\mathbf{x}^{*}$ is

[TABLE]

which can be interpreted as the generalized LASSO [27] with arbitrary forms of objective functions. For a fair comparison, we configure the regularization parameter $\gamma$ so that the expected sparsity of the FISTA estimate is $3L$, a criterion suggested in [13]. Since GraSP, GraHTP, and BG-GAMP break down when the sensing matrix is ill-conditioned, we set $B_{\mathrm{RX}}=B_{\mathrm{TX}}=64$ for these channel estimators. For FISTA, we choose $B_{\mathrm{RX}}=B_{\mathrm{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_{\mathrm{RX}}$ and $B_{\mathrm{TX}}$, but increasing these parameters is forbidden since these channel estimators diverge when $\mathbf{A}$ is highly coherent. FISTA performs poorly because the Laplace distribution deviates from the true distribution of $\mathbf{x}^{*}$. In contrast, since $B_{\mathrm{RX}}$ and $B_{\mathrm{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 $\mathbf{x}^{*}$ and $c\mathbf{x}^{*}$ are indistinguishable in the high SNR regime for $c>0$ because

[TABLE]

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 $\hat{\mathbf{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 sparsity-constrained 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 & communications Technology Promotion(IITP) grant funded by the Korea government(MSIT) (No. 2016-0-00123, Development of Integer-Forcing MIMO Transceivers for 5G & Beyond Mobile Communication Systems) and by the National Research Foundation(NRF) grant funded by the MSIT of the Korea government (No. 2018R1A4A1025679).