On Matrix Completion-Based Channel Estimators for Massive MIMO Systems

Abstract

Large-scale symmetric arrays such as uniform linear arrays (ULA) have been widely used in wireless communications for improving spectrum efficiency and reliability. Channel state information (CSI) is critical for optimizing massive multiple-input multiple-output(MIMO)-based wireless communication systems. The acquisition of CSI for massive MIMO faces challenges such as training shortage and high computational complexity. For millimeter wave MIMO systems, the low-rankness of the channel can be utilized to address the challenge of training shortage. In this paper, we compared several channel estimation schemes based on matrix completion (MC) for symmetrical arrays. Performance and computational complexity are discussed and compared. By comparing the performance in different scenarios, we concluded that the generalized conditional gradient with alternating minimization (GCG-Alt) estimator provided a low-cost, robust solution, while the alternating direction method of multipliers (ADMM)-based hybrid methods achieved the best performance when the array response was perfectly known.

1. Introduction

Millimeter-wave (mmWave) wireless communications have drawn great attention in the industry and academia [1] thanks to the large bandwidth available in the 30–300 GHz band. To compensate for the significant path loss in this band and also thanks to the short wavelength, massive MIMO have been suggested for mmWave systems. In particular, large-scale symmetric antenna arrays, such as the uniform linear arrays (ULAs) and uniform planer arrays (UPAs) have been extensively considered for transmitters and receivers due to their neat structures and high gains for directional transmissions [2,3,4]. However, the coherence time in the millimeter-wave system is suggested to be short and as the number of antennas increases, the complexity of channel estimation increases. Therefore, it is challenging to acquire instantaneous channel state information (CSI) for a mmWave massive MIMO.

MmWave channels are often dominated by a small number of propagation paths, indicating that the channel is sparse in the angular domain [5]. The channel matrix can be expressed in terms of dictionary matrices, which are formed by the transmitting and receiving array response vectors, and path gains. Compressive sensing (CS) [6] can then be applied to search for the dominant paths [7,8,9,10,11,12,13]. Different measurement matrices can be used by choosing the precoders and combiners, as well as various recovery algorithms such as orthogonal matching pursuit (OMP) [7,14] and the adaptive CS [8,9,10] can be applied. In general, the above mentioned CS schemes require knowledge of the array response, which depends on the array geometry and calibration of the antenna arrays. Such knowledge can be inaccurate when there are unknown hardware impairments, e.g., due to phase and gain errors, and imperfect calibration of the antenna arrays.

In the meantime, a small number of propagation paths also indicate that the channel is low-rank [15] and can be depicted in low-dimensional subspace. Furthermore, such low-rankness is independent of the array response and calibration errors. In [16,17], both the sparsity and low-rank property of the mmWave channel are exploited to enhance CS-based channel estimators. In [16], a two-stage estimator is proposed, where the low-rankness is exploited at the first stage while sparsity in the angular domain is exploited at the second stage. In [17], the improved alternating direction method of multipliers (ADMM) method [18] is applied to exploit the low-rankness and sparsity while at the same time enhancing the performance. These estimators, however, still require knowledge of the array response vectors. They can still suffer from performance loss when there are uncertainties in the array response.

To achieve robust channel estimation, matrix completion (MC) methods exploiting only the low-rank property of the mmWave channel have recently been proposed [19,20]. The analysis in [21] shows that the rank of the channel matrix is generally very low, which is usually much smaller than the antenna dimension. In [19], the singular value projection (SVP) algorithm [22] is adopted to solve the mmWave channel estimation problem. Later, the GCG-Alt method is developed in [20] by applying the generalized conditional gradient (GCG) framework [23] together with the alternating minimum (AltMin) method [24]. There are also other widely-studied MC algorithms that can be used for mmWave channel estimation, such as the singular value thresholding algorithm (SVT) [25] and the fixed point continuation algorithm (FPC) [26]. In this paper, we discuss several mmWave channel estimators based on MC, focusing on their performance and complexity comparisons with alternative methods based on CS. We aim to examine the pros and cons for several MC estimators and the factors that influence their performances.

The rest of this paper is organized as follows. Section 2 introduces the mmWave MIMO channel model and formulates the channel estimation problem. Section 3 introduces channel estimators based on MC, including their detailed implementation. Section 4 presents simulation results, in terms of the mean squared error (MSE) and computational complexity. Section 5 concludes the paper.

Notation: ${\mathbf{A}}^{\mathrm{T}}$, ${\mathbf{A}}^{*}$, and ${\mathbf{A}}^{\mathrm{H}}$ denote transpose, conjugate, and conjugate transpose of matrix $\mathbf{A}$, respectively. ${\parallel \mathbf{A}\parallel }_1$ denotes the $l_1$ norm. $\mathbf{I}$ and $\mathbf{0}$ represent the identity matrix and zero matrix/vector, respectively. $\mathbf{A}\otimes \mathbf{B}$ and $\mathbf{A}\odot \mathbf{B}$ denote the Kronecker product and the Hadamard product, respectively. $tr\left\{\mathbf{A}\right\}$ is the trace of $\mathbf{A}$ and $\langle \mathbf{A},\mathbf{C}\rangle =tr\left\{{\mathbf{A}}^{\mathrm{H}}\mathbf{C}\right\}$ denotes the inner product of matrices $\mathbf{A}$ and $\mathbf{C}$. $E[·]$ denotes the statistical expectation and $\mathrm{abs}(·)$ represents taking element-wise absolute value. ∇ denotes the gradient of a function. For a matrix $\mathbf{A}\in {\mathbb{C}}^{M\times N},vec\left(\mathbf{A}\right)\in {\mathbb{C}}^{MN\times 1}$ denotes the vectorization of $\mathbf{A}$ and ${vec}^{-1}\left(\mathbf{A}\right)\in {\mathbb{C}}^{M\times N}$ denotes the inverse of vectorization. $\mathcal{R}(·)$ and $\mathcal{I}(·)$ denote the real and imaginary part of a number or vector, respectively. $\mathcal{CN}\left(a,b^2\right)$ denotes complex Gaussian distribution with mean a and variance $b^2$.

2. System Model

Consider a point-to-point, switch-based mmWave hybrid MIMO system, with the receiver at the mobile station (MS) shown in Figure 1. For simplicity and clarity, this paper assumes switch-based mmWave systems to investigate MC-based channel estimators. MC-based estimators can also be applied to phase shifter-based mmWave systems, when the hybrid precoders/combiners are properly designed, as shown in [20]. Therefore, the discussion in this paper can be easily extended to phase shifter-based mmWave systems. At the receiver, each of the $N_{\mathrm{MS}}$ antennas is equipped with a switch used to select one of the $N_{\mathrm{RF}}$ RF chains. The base station (BS) has the same structure with $N_{\mathrm{BS}}$ antennas and $N_{\mathrm{RF}}$ RF chains. Assume that $N_{\mathrm{s}}$ data streams are transmitted, with $N_{\mathrm{s}}\le N_{\mathrm{RF}}\le min(N_{\mathrm{BS}},N_{\mathrm{MS}})$ [14,27]. The switching operation can be represented as a precoder $\mathbf{F}$, where the nonzero entries indicate the entries of the channel matrix that are sampled. A symbol $\mathbf{s}\in {\mathbb{C}}^{\phantom{R^{{}^{\prime }}}{N}_{\mathrm{s}}\times 1}$ with $\mathbb{E}\left[{\mathbf{ss}}^{\mathrm{H}}\right]=\frac{1}{N_{\mathrm{s}}}\mathbf{I}$ is precoded, resulting in the transmitted signal $\mathbf{x}=\mathbf{Fs}$. We consider a narrow-band flat fading channel whose channel matrix $\mathbf{H}\in {\mathbb{C}}^{N_{\mathrm{MS}}\times N_{\mathrm{BS}}}$ satisfies $\mathbb{E}\left[{\parallel \mathbf{H}\parallel }_{\mathrm{F}}^2\right]=N_{\mathrm{MS}}{N}_{\mathrm{BS}}$. The received signal is expressed as:

$$\mathbf{r}=\sqrt{\rho }\mathbf{HFs}+\widehat{\mathbf{n}},$$

(1)

where $\rho$ indicates the average received power and $\widehat{\mathbf{n}}\in {\mathbb{C}}^{\phantom{R^{{}^{\prime }}}{N}_{\mathrm{MS}}\times 1}$ is a noise vector with i.i.d. entries distributed as $\mathcal{CN}(0,{\sigma }_n^2)$. Applying a combiner $\mathbf{W}$ to the received signal at the MS, the processed received signal is given by:

$$\mathbf{y}=\sqrt{\rho }{\mathbf{W}}^{\mathrm{H}}\mathbf{H}\mathbf{F}\mathbf{s}+{\mathbf{W}}^{\mathrm{H}}\widehat{\mathbf{n}}.$$

(2)

In the switch-based system, the combiner $\mathbf{W}$ has a similar structure with the precoder $\mathbf{F}$.

Following [8], the ray-clustering model of $\mathbf{H}$ is given as:

$${\displaystyle \mathbf{H}=\frac{1}{\sqrt{R}}\sum _{c=1}^C\sum _{r=1}^R{g}_{cr}{\mathbf{a}}_{\mathrm{MS}}\left({\varphi }_{cr}^{\mathrm{MS}}\right){\mathbf{a}}_{\mathrm{BS}}^{\mathrm{H}}\left({\varphi }_{cr}^{\mathrm{BS}}\right),}$$

(3)

where $C\sim \max \left\{\mathrm{Poisson}\right(\lambda ),1\}$ is the number of clusters with $\lambda$ as the mean of the Poisson distribution, and R is the number of rays in each cluster. The complex small-scale fading gain on the r-th ray of the c-th cluster is $g_{cr}$ with $g_{cr}\sim \mathcal{CN}(0,{\gamma }_c)$, where ${\gamma }_c$ is the sub-power on the c-th cluster. In Equation (3), ${\mathbf{a}}_{\mathrm{MS}}\left({\varphi }_{cr}^{\mathrm{MS}}\right)$ and ${\mathbf{a}}_{\mathrm{BS}}\left({\varphi }_{cr}^{\mathrm{BS}}\right)$ represent the array response vector for the receiver and transmitter, respectively, where ${\varphi }_{cr}^{\mathrm{MS}}$ and ${\varphi }_{cr}^{\mathrm{BS}}$ represent the corresponding azimuth AoA and AoD, which follow the Laplacian distribution [28].

Considering a uniform linear array (ULA) with distance between adjacent antennas being d, the array response is given by:

$$\begin{array}{c}\hfill {\mathbf{a}}_{\mathrm{BS}}\left({\varphi }_{cr}^{\mathrm{BS}}\right)=\frac{1}{\sqrt{N_{\mathrm{BS}}}}{\left[1,{\mathrm{e}}^{j\frac{2\pi }{{\lambda }_{\mathrm{c}}}dsin\left({\varphi }_{cr}^{\mathrm{BS}}\right)},\cdots ,{\mathrm{e}}^{j(N-1)\frac{2\pi }{{\lambda }_{\mathrm{c}}}dsin\left({\varphi }_{cr}^{\mathrm{BS}}\right)}\right]}^T\end{array}$$

(4)

where ${\lambda }_{\mathrm{c}}$ is the wavelength of the carrier wave. The array response ${\mathbf{a}}_{\mathrm{MS}}\left({\varphi }_{cr}^{\mathrm{MS}}\right)$ is constructed in the same manner as ${\mathbf{a}}_{\mathrm{BS}}\left({\varphi }_{cr}^{\mathrm{BS}}\right)$.

3. Compressive Sensing-Based Channel Estimation

It has been shown that, without considering quantization errors, the mmWave channel estimation problem can be formulated as a sparse recovery problem by modeling the channel as [29,30,31]:

$$\mathbf{H}={\mathbf{A}}_{\mathrm{MS}}{\mathbf{H}}_{\mathrm{v}}{\mathbf{A}}_{\mathrm{BS}}^{\mathrm{H}},$$

(5)

where ${\mathbf{A}}_{\mathrm{MS}}=[{\mathbf{a}}_{\mathrm{MS}}\left({\varphi }_1^{\mathrm{MS}}\right),\dots ,{\mathbf{a}}_{\mathrm{MS}}\left({\varphi }_{N_1}^{\mathrm{MS}}\right)]$ is a unitary dictionary matrix when $N_1=N_{\mathrm{MS}}$ and it is an overcomplete dictionary matrix when $N_1>N_{\mathrm{MS}}$, and ${\mathbf{A}}_{\mathrm{BS}}=[{\mathbf{a}}_{\mathrm{BS}}\left({\varphi }_1^{\mathrm{BS}}\right),\dots ,{\mathbf{a}}_{\mathrm{BS}}\left({\varphi }_{N_2}^{\mathrm{BS}}\right)]$ with $N_2\ge N_{\mathrm{BS}}$. Each column in ${\mathbf{A}}_{\mathrm{BS}}$ and ${\mathbf{A}}_{\mathrm{MS}}$ consists of a predefined array response vector. ${\mathbf{H}}_{\mathrm{v}}\in {\mathbb{C}}^{N_1\times N_2}$ is a sparse matrix with only L non-zero values, with each of its non-zero values corresponding to the complex gain of a channel path. Vectorization of the channel matrix (5) produces:

$$\mathrm{vec}\left(\mathbf{H}\right)=({\mathbf{A}}_{\mathrm{BS}}^{*}\otimes {\mathbf{A}}_{\mathrm{MS}})\mathbf{x},$$

(6)

where $\mathbf{x}=\mathrm{vec}\left({\mathbf{H}}_{\mathrm{v}}\right)$ is a $N_1{N}_2\times 1$ sparse vector with L non-zero values. We define $\mathsf{\Psi }=({\mathbf{A}}_{\mathrm{BS}}^{*}\otimes {\mathbf{A}}_{\mathrm{MS}})$ as a $N_{\mathrm{BS}}{N}_{\mathrm{MS}}\times N_1{N}_2$ dictionary matrix. Sparse recovery schemes can then be used to estimate the channel, which transforms the task of estimating $\mathbf{H}$ to estimating the non-zero coefficients in $\mathbf{x}$.

A widely used method to estimate $\mathbf{x}$ from the the received signal is orthogonal matching pursuit (OMP) [7,14]. Using OMP, L path directions from the $N_1{N}_2$ candidates in the dictionary are determined. The computational complexity of the OMP method is approximately $O\left(NL{N}_1{N}_2\right)$, where N is the length of the received signal. In general, a larger dictionary leads to better performance but also higher computational complexity.

The above mentioned CS approach uses a discretized approximation of the channel. It may suffer from the off-grid issue if the physical propagation paths are off the assumed grid of the angles. In this case, the number of non-zero entries in the beamspace channel ${\mathbf{H}}_{\mathrm{v}}$ may not be exactly equal to L, leading to a power leakage. Another challenge is that the knowledge of the array response is required, which may be imperfect in practice due to unknown hardware impairments and imperfect calibrations.

4. Matrix Completion-Based Channel Estimation

In this section, we introduce MC-based estimation methods for the mmWave channel by exploiting the low-rankness of the channel matrix. By appropriately choosing the training scheme with proper precoders and combiners, the received signal provides noisy observations of a subset of the entries of $\mathbf{H}$:

$${\left[\tilde{\mathbf{Y}}\right]}_{i,j}=\left\{\begin{array}{cc}{\left[\tilde{\mathbf{H}}\right]}_{i,j},\hfill & (i,j)\in \mathsf{\Omega }\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(7)

where $\tilde{\mathbf{H}}=\mathbf{H}+\mathbf{N}$ is the perturbed channel matrix, $\mathbf{N}$ is a noise matrix, $\mathsf{\Omega }$ denotes a sample domain, and ${\left[\tilde{\mathbf{H}}\right]}_{i,j}$ is the $\left(i,j\right)$th entry of $\tilde{\mathbf{H}}$. Define $p=N/\left(N_{\mathrm{BS}}{N}_{\mathrm{MS}}\right)$ as the sampling density, where N is the total number of samples observed.

It is discussed in [19] that when the mmWave channel matrix has strong non-coherent characteristics, it can be recovered from a subset samples of the channel matrix. We can thus formulate the channel estimation problem as a low-rank matrix completion problem as:

$$\underset{\widehat{\mathbf{H}}}{min}\mathrm{rank}\left(\widehat{\mathbf{H}}\right),\mathrm{s}.\mathrm{t}.\left|\right|P_{\mathsf{\Omega }}\left(\widehat{\mathbf{H}}\right)-P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right){\left|\right|}_F^2\le {\delta }_n^2,$$

(8)

where ${\delta }_n^2$ is the error tolerance parameter and the sampling operator $P_{\mathsf{\Omega }}(·)$ is defined as:

$${\left[P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)\right]}_{i,j}=\left\{\begin{array}{cc}{\left[\tilde{\mathbf{H}}\right]}_{i,j}\hfill & (i,j)\in \mathsf{\Omega }\hfill \\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}\right.$$

(9)

where ${\left[\tilde{\mathbf{H}}\right]}_{i,j}$ denotes the $(i,j)$-th entry of $\tilde{\mathbf{H}}$. The sampling operator $P_{\mathsf{\Omega }}$ significantly influences the performance of the algorithm [32]. Bernoulli and uniform sampling models are proposed in [32] and a uniform spatial sampling model (USS), which improves the performance, is proposed in [33]. With USS, $N/N_{\mathrm{BS}}$ samples are taken for each column of the target matrix.

4.1. MC Estimators

In the following, we discuss several MC estimators that can be used to solve the problem in Equation (8).

4.1.1. SVT Estimator

Before presenting the SVT algorithm, let us define the matrix shrinkage operator:

$${\mathcal{S}}_{\tau }\left(\mathbf{X}\right)={\mathbf{U}}_{\mathrm{X}}{\mathcal{S}}_{\tau }\left({\mathsf{\Sigma }}_{\mathrm{X}}\right){\mathbf{V}}_{\mathrm{X}}^{\mathrm{H}},$$

(10)

where ${\mathsf{\Sigma }}_{\mathrm{X}}$ denotes the singular value matrix of $\mathbf{X}$ and ${\mathcal{S}}_{\tau }\left({\mathsf{\Sigma }}_{\mathrm{X}}\right)$ is the element-wise shrinkage operator:

$${\mathcal{S}}_{\tau }\left(x\right)=\left\{\begin{array}{cc}0,\hfill & x\le \tau \hfill \\ x-\tau ,\hfill & x>\tau \hfill \end{array}\right.$$

(11)

where $\tau$ is the threshold.

The SVT algorithm [25] can be applied to provide a heuristic solution to Equation (8), which consists of two major steps,

$$\left\{\begin{array}{c}{\widehat{\mathbf{H}}}^k={\mathcal{S}}_{\tau }\left({\mathbf{X}}^{k-1}\right)\hfill \\ {\mathbf{X}}^k={\mathbf{X}}^{k-1}+\delta P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}-{\widehat{\mathbf{H}}}^k\right)\hfill \end{array}\right.$$

(12)

where $\tau >0$, $\delta$ is a step size, and $k=1,2,\dots$. The iteration is stopped when a stopping criterion is met or a maximum number of iterations $J_{\mathrm{SVT}}$ is reached. Singular value decomposition (SVD) is required at each iteration. Some comments regarding the implementation of the SVT algorithm to the mmWave channel estimation problem are as follows:

• The threshold is set as $\tau =5\sqrt{N_{\mathrm{MS}}{N}_{\mathrm{BS}}}$ following [25];
• The stepsize is set as $\delta =1.2/p$;
• Assuming the initialization ${\mathbf{X}}^0=\mathbf{0}$, ${\widehat{\mathbf{H}}}^k=\mathbf{0}$ for a small $k<k_0$. As such, ${\mathbf{X}}^k=k\delta {\mathcal{P}}_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right),k=1,\dots ,k_0$. The algorithm can begin with computing ${\widehat{\mathbf{H}}}^{k_0}$ to save work. From [25], the integer $k_0$ is determined by

  $$\frac{\tau }{\delta {∥P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)∥}_2}\in \left(k_0-1,k_0\right].$$

SVT Estimator is shown in below Algorithm 1:

Algorithm 1 SVT Estimator

Require:  $P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)$, $\delta$, $ϵ$, $\tau$, $J_{\mathrm{SVT}}$, $k_0$.

1: Set ${\mathbf{X}}^0=k_0\delta P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)$; 2: for$k=1$ to $J_{\mathrm{SVT}}$ do 3:     Set ${\widehat{\mathbf{H}}}^k={\mathcal{S}}_{\tau }\left({\mathbf{X}}^{k-1}\right)$ 4:     Set ${\mathbf{X}}^k={\mathbf{X}}^{k-1}+\delta P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}-{\widehat{\mathbf{H}}}^k\right)$ 5:     if $\frac{\parallel P_{\mathsf{\Omega }}({\widehat{\mathbf{H}}}^k-\tilde{\mathbf{H}}){\parallel }_{\mathrm{F}}}{\parallel P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right){\parallel }_{\mathrm{F}}}\le ϵ$ then break; 6:     end if 7: end for 8: return$\widehat{\mathbf{H}}={\widehat{\mathbf{H}}}^k$

From [25], SVT is effective for completing large matrices with low ranks. Its performance degrades as the rank increases. The computational complexity of the SVT algorithm mainly arises from Step 3.

4.1.2. FPC Estimator

The FPC algorithm [26] reformulates the MC problem using the nuclear norm, which is the summation of the singular values,

$$\underset{\widehat{\mathbf{H}}}{min}\mu {∥\widehat{\mathbf{H}}∥}_{*}+\frac{1}{2}{\parallel P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)-P_{\mathsf{\Omega }}\left(\widehat{\mathbf{H}}\right)\parallel }_{\mathrm{F}}^2,$$

(13)

where $\mu >0$ is the regularization parameter. The algorithm consists of two steps similar to SVT:

$$\left\{\begin{array}{c}{\mathbf{Y}}^k={\widehat{\mathbf{H}}}^k-\delta (P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)-P_{\mathsf{\Omega }}\left(\widehat{\mathbf{H}}\right))\hfill \\ {\widehat{\mathbf{H}}}^{k+1}={\mathcal{S}}_{\delta {\mu }_m}\left({\mathbf{Y}}^k\right)\hfill \end{array}\right.$$

(14)

where the threshold of the singular value thresholding operator is set as a variable $\delta {\mu }_m$ rather than a fixed value as that in SVT. A continuous strategy [34] is used to accelerate the convergence by adapting ${\mu }_m$. The details are presented in Algorithm 2. Some comments are as below:

• We can set the step size $\delta \in \left(0,2/{\lambda }_{max}\left(P_{\mathsf{\Omega }}{\left(\tilde{\mathbf{H}}\right)}^{\mathrm{H}}{P}_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)\right)\right)$ according to [26], where ${\lambda }_{max}$ is the maximum eigenvalue;
• ${\mu }_{m+1}$ decreases as

  $${\mu }_{m+1}=max\left\{{\mu }_m{\eta }_{\mu },{\mu }_{\mathrm{final}}\right\},m=1,2,\cdots ,M,$$

  where M, determined by $({\mu }_{\mathrm{final}},{\eta }_{\mu })$, controls the step size and the estimation accuracy, ${\mu }_{\mathrm{final}}$ is small (e.g., ${\mu }_{\mathrm{final}}=0.01$), and the parameter $0<{\eta }_{\mu }<1$ determines the decreasing rate for consecutive ${\mu }_m$.

Algorithm 2 FPC Estimator

Require: $P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)$, $ϵ$, $J_{\mathrm{FPC}}$, $\delta$, ${\mu }_{\mathrm{final}}$, and ${\eta }_{\mu }$

1: Initialization: ${\widehat{\mathbf{H}}}^0=\mathbf{0}$, $m=0$, ${\mu }_m={∥P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)∥}_2$ 2: while${\mu }_m>{\mu }_{\mathrm{final}}$do 3:     ${\mu }_m=max({\mu }_m{\eta }_{\mu },{\mu }_{\mathrm{final}})$ 4:     for $k=1:J_{\mathrm{FPC}}$ do 5:         ${\mathbf{Y}}^k={\widehat{\mathbf{H}}}^k-\delta \left(P_{\mathsf{\Omega }}(\tilde{\mathbf{H}}-{\widehat{\mathbf{H}}}^k)\right)$ 6:         ${\widehat{\mathbf{H}}}^k={\mathcal{S}}_{\delta {\mu }_m}\left({\mathbf{Y}}^k\right)$ 7:         if $\frac{{∥{\widehat{\mathbf{H}}}^{k+1}-{\widehat{\mathbf{H}}}^k∥}_{\mathrm{F}}}{{∥{\widehat{\mathbf{H}}}^k∥}_{\mathrm{F}}}<ϵ$ then break; 8:         end if 9:     end for 10: end while 11: return$\widehat{\mathbf{H}}={\widehat{\mathbf{H}}}^k$

The main computational cost of the FPC algorithm is in Step 6 of Algorithm 2 due to the SVD. In addition, the FPC algorithm needs to choose the step size $\delta$ by calculating the maximum eigenvalue of $P_{\mathsf{\Omega }}{\left(\tilde{\mathbf{H}}\right)}^{\mathrm{H}}{P}_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)$ and has a higher computational complexity per iteration than that of the SVT algorithm.

4.1.3. SVP Estimator

The SVP algorithm [22] is based on the projected gradients and is detailed in Algorithm 3. This algorithm requires that the rank L of the channel matrix to be known. The step size can be chosen empirically as $\eta =1/(1+{\delta }_0)p$ with $0<{\delta }_0<1/3$. The stopping criterion is based on the norm of the difference in the sampled channel matrix, where the small threshold $0<ϵ<1/2$ can be set such as $ϵ={10}^{-3}$. The SVP algorithm also needs to calculate the SVD in Step 4, which is the most computationally expensive step of the algorithm.

Algorithm 3 SVP Estimator

Require: $P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)$, L, $\eta$, $ϵ$

1: $\mathbf{Initialization}$: ${\widehat{\mathbf{H}}}^0=\mathbf{0},t=0$ 2: while$\parallel P_{\mathsf{\Omega }}({\widehat{\mathbf{H}}}^t-\tilde{\mathbf{H}}){\parallel }_{\mathrm{F}}\le ϵ$do 3:     ${\mathbf{X}}^{t+1}\leftarrow {\widehat{\mathbf{H}}}^t-\eta \left(P_{\mathsf{\Omega }}({\widehat{\mathbf{H}}}^t-\tilde{\mathbf{H}})\right)$ 4:     Compute the L principal singular vectors of ${\mathbf{X}}^{t+1}$ : ${\mathbf{U}}_{\mathrm{L}},{\mathsf{\Sigma }}_{\mathrm{L}},{\mathbf{V}}_{\mathrm{L}}$. 5:     ${\widehat{\mathbf{H}}}^{t+1}={\mathbf{U}}_{\mathrm{L}}{\mathsf{\Sigma }}_{\mathrm{L}}{\mathbf{V}}_{\mathrm{L}}^H$ 6:     $t=t+1$; 7: end while 8: return$\widehat{\mathbf{H}}={\widehat{\mathbf{H}}}^k$

4.1.4. GCG-Alt Estimator

In [20], a generalized conditional gradient framework with alternating minimization (GCG-Alt) is developed for the MC problem. The nuclear norm is used to promote low-rankness as:

$$\underset{\widehat{\mathbf{H}}}{min}\frac{1}{2}\parallel P_{\mathsf{\Omega }}(\widehat{\mathbf{H}}-\tilde{\mathbf{H}}){\parallel }_{\mathrm{F}}^2+\mu {\parallel \widehat{\mathbf{H}}\parallel }_{*}$$

(15)

where $\mu >0$ is a regularization factor. Let:

$$f\left(\widehat{\mathbf{H}}\right)\triangleq \frac{1}{2}{∥P_{\mathsf{\Omega }}\left(\widehat{\mathbf{H}}\right)-P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)∥}_{\mathrm{F}}^2.$$

(16)

The gradient direction of the kth iteration of $f\left(\widehat{\mathbf{H}}\right)$ [23]:

$${\mathbf{X}}_k={\mathbf{u}}_{k-1}{\mathbf{v}}_{k-1}^{\mathrm{H}},$$

(17)

where $\left({\mathbf{u}}_{k-1},{\mathbf{v}}_{k-1}\right)$ is the top singular vector pair of $\nabla f\left({\widehat{\mathbf{H}}}_{k-1}\right)=P_{\mathsf{\Omega }}\left({\widehat{\mathbf{H}}}_{k-1}-\tilde{\mathbf{H}}\right)$ which can be found iteratively. The channel matrix is updated as:

$${\widehat{\mathbf{H}}}_k=(1-{\eta }_k){\widehat{\mathbf{H}}}_k+{\theta }_k{\mathbf{X}}_k,$$

(18)

where ${\eta }_k\in [0,1]$ is the step size and ${\theta }_k$ is adaptively chosen.

By using a property of nuclear norm [20], the optimization problem can be reformulated as:

$$\tilde{\varphi }(\mathbf{U},\mathbf{V})\triangleq f\left(\mathbf{U}{\mathbf{V}}^{\mathrm{H}}\right)+\frac{1}{2}\mu \left({\parallel \mathbf{U}\parallel }_{\mathrm{F}}^2+{\parallel \mathbf{V}\parallel }_{\mathrm{F}}^2\right),$$

(19)

where $\mathbf{U}\in {\mathbb{C}}^{N_{\mathrm{MS}}\times \widehat{r}}$ and $\mathbf{V}\in {\mathbb{C}}^{N_{\mathrm{BS}}\times \widehat{r}}$ with $\widehat{r}$ being the rank of $\widehat{\mathbf{H}}$. Alternating minimization can then be used to optimize Equation (19). The details of solving the alternate minimization problem can be found in [20]. The overall algorithm is summarized in Algorithm 4.

Algorithm 4 GCG-Alt Estimator

Require: $P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right),\mu ,ϵ,{ϵ}_{\mathrm{a}}$

1: Initialization:${\mathbf{U}}_0=⌀,{\mathbf{V}}_0=⌀,k=0,{\widehat{\mathbf{H}}}_k=\mathbf{0},{ϵ}_0=\infty$, ${\delta }_k^2={∥P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)∥}_F^2$ 2: while${\delta }_k^2>(N+\sqrt{8N}){\sigma }^2$do 3:     Compute the top singular vector pair $({\mathbf{u}}_k,{\mathbf{v}}_k)$ of $P_{\mathsf{\Omega }}\left({\widehat{\mathbf{H}}}_k-\tilde{\mathbf{H}}\right)$ 4:     $k=k+1$ 5:     ${\eta }_k\leftarrow 2/(k+1)$ 6:     ${\mathbf{x}}_{k\mathsf{\Omega }}=vec\left(P_{\mathsf{\Omega }}\left({\mathbf{u}}_k{\mathbf{v}}_k^{\mathrm{H}}\right)\right)$ 7:     ${\tilde{\mathbf{h}}}_{\mathsf{\Omega }}=vec\left(P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)\right)$ 8:     ${\widehat{\mathbf{h}}}_{k\mathsf{\Omega }}=vec\left(P_{\mathsf{\Omega }}\left({\widehat{\mathbf{H}}}_{k-1}\right)\right)$ 9:     ${\theta }_k=\frac{2\mathcal{R}\left({\mathbf{x}}_{k\mathsf{\Omega }}^{\mathrm{H}}{\tilde{\mathbf{h}}}_{\mathsf{\Omega }}\right)-\left(1-{\eta }_{\mathrm{k}}\right){\mathbf{x}}_{k\mathsf{\Omega }}^{\mathrm{H}}{\widehat{\mathbf{h}}}_{k\mathsf{\Omega }}-2\mu }{2{\mathbf{x}}_{k\mathsf{\Omega }}^{\mathrm{H}}{\mathbf{x}}_{k\mathsf{\Omega }}}$ 10:     ${\mathbf{U}}_k=\left[\sqrt{1-{\eta }_k}{\mathbf{U}}_{k-1},\sqrt{{\theta }_k}{\mathbf{u}}_{k-1}\right]$ 11:     ${\mathbf{V}}_k=\left[\sqrt{1-{\eta }_k}{\mathbf{V}}_{k-1},\sqrt{{\theta }_k}{\mathbf{v}}_{k-1}\right]$ 12:     $\mathbf{Initialization}$: $i=0,{ϵ}_k^0=\infty$, $\left({\mathbf{U}}_k^0,{\mathbf{V}}_k^0\right)\leftarrow \left({\mathbf{U}}_k,{\mathbf{V}}_k\right)$. 13:     while ${ϵ}_k^i>{ϵ}_{\mathrm{a}}$ do 14:         $i=i+1$ 15:         Find ${\mathbf{V}}_k^i$ that minimizes Equation (19) given $\mathbf{U}={\mathbf{U}}_k^i$; 16:         Find ${\mathbf{U}}_k^i$ that minimizes Equation (19) given $\mathbf{V}={\mathbf{V}}_k^i$; 17:         Calculate $${ϵ}_k^i=\frac{\tilde{\varphi }\left({\mathbf{U}}_k^{i-1},{\mathbf{V}}_k^{i-1}\right)-\tilde{\varphi }\left({\mathbf{U}}_k^i,{\mathbf{V}}_k^i\right)}{\tilde{\varphi }\left({\mathbf{U}}_k^{i-1},{\mathbf{V}}_k^{i-1}\right)}$$ 18:     end while 19:     $({\mathbf{U}}_k,{\mathbf{V}}_k)\leftarrow ({\mathbf{U}}_k^i,{\mathbf{V}}_k^i)$ 20:     Calculate ${ϵ}_k=\frac{{∥{\widehat{\mathbf{H}}}_k∥}_{\mathrm{F}}^2-{∥{\widehat{\mathbf{H}}}_{k-1}∥}_{\mathrm{F}}^2}{{∥{\widehat{\mathbf{H}}}_{k-1}∥}_{\mathrm{F}}^2}$ 21:     ${\delta }_k^2={∥P_{\mathsf{\Omega }}({\widehat{\mathbf{H}}}_k-\tilde{\mathbf{H}})∥}_{\mathrm{F}}^2$ 22: end while 23: return$\widehat{\mathbf{H}}={\mathbf{U}}_k{\mathbf{V}}_k^{\mathrm{H}}$

The above MC methods have different computational complexities. SVT, SVP, and FPC all need SVD, which can be implemented using the PROPACK [35] based on the iterative Lanczos bidiagonalization algorithm with partial reorthogonalization. The FPC has a higher complexity as SVD is repeated for different values of ${\mu }_m$. The GCG-Alt has the least complexity as the full SVD is not required. SVP is effective for large matrix completion problems with very low ranks, while the FPC, SVT, and GCG-Alt estimators allow higher ranks. The SVP estimator requires rank knowledge, while the FPC and GCG-Alt estimators implicitly determine the rank through choosing regularization parameters or thresholds.

4.2. MC-Based Hybrid Estimators

Next we discuss two MC-based hybrid methods that jointly exploit the sparsity and low-rankness of the channel.

4.2.1. ADMM Estimator

In [17], the low-rankness of $\tilde{\mathbf{H}}$ and the sparsity of the beamspace channel ${\tilde{\mathbf{H}}}_{\mathrm{v}}$ are jointly exploited and an ADMM method is proposed. Leveraging the side information that $\mathbf{H}$ has a sparse virtual representation given by Equation (5), the channel estimation problem is formulated following [36] as:

$$\begin{array}{cc}& {\displaystyle \underset{\widehat{\mathbf{H}},{\widehat{\mathbf{H}}}_{\mathrm{v}}}{min}}{\tau }_{\mathrm{L}}\parallel \widehat{\mathbf{H}}{\parallel }_{*}+{\tau }_{\mathrm{S}}{\parallel {\widehat{\mathbf{H}}}_{\mathrm{v}}\parallel }_1\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& P_{\mathsf{\Omega }}\left(\widehat{\mathbf{H}}\right)=P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)\mathrm{and}\widehat{\mathbf{H}}={\mathbf{A}}_{\mathrm{MS}}{\widehat{\mathbf{H}}}_{\mathrm{v}}{\mathbf{A}}_{\mathrm{BS}}^{\mathrm{H}}\hfill \end{array}$$

(20)

where the nuclear norm and $l_1$-norm together with the regularization parameters ${\tau }_{\mathrm{L}}$ and ${\tau }_{\mathrm{S}}$ are used to promote low-rankness and sparsity, respectively. The above problem is then reformulated by incorporating the constraints as penalty terms as:

$$\begin{array}{ll}& \underset{\widehat{\mathbf{H}},\mathbf{E},{\widehat{\mathbf{H}}}_{\mathrm{v}},\mathbf{C}}{min}{\tau }_{\mathrm{L}}\parallel \widehat{\mathbf{H}}{\parallel }_{*}+{\tau }_{\mathrm{S}}{∥{\widehat{\mathbf{H}}}_{\mathrm{v}}∥}_1+\frac{1}{2}{\parallel \mathbf{C}\parallel }_{\mathrm{F}}^2+\frac{1}{2}{∥P_{\mathsf{\Omega }}(\mathbf{E}-\widehat{\mathbf{H}})∥}_{\mathrm{F}}^2\\ \mathrm{s}.\mathrm{t}.& \widehat{\mathbf{H}}=\mathbf{E}\mathrm{and}\mathbf{C}=\mathbf{E}-{\mathbf{A}}_{\mathrm{MS}}{\widehat{\mathbf{H}}}_{\mathrm{v}}{\mathbf{A}}_{\mathrm{BS}}^{\mathrm{H}}.\end{array}$$

(21)

where $\mathbf{E}\in {\mathbb{C}}^{N_{\mathrm{MS}}\times N_{\mathrm{BS}}}$ and $\mathbf{C}$ are two auxiliary matrix variables. This problem is then solved by using ADMM which involves the iterative updates of the variables and Lagrangian multipliers. The augmented Lagrangian function of Equation (21) is given by:

$$\begin{array}{lll}{\mathcal{L}}_1\left(\widehat{\mathbf{H}},\mathbf{E},{\widehat{\mathbf{H}}}_{\mathrm{v}},\mathbf{C},{\mathbf{Z}}_1,{\mathbf{Z}}_2\right)& \triangleq & {\tau }_{\mathrm{L}}\parallel \widehat{\mathbf{H}}{\parallel }_{*}+{\tau }_{\mathrm{S}}\parallel {\widehat{\mathbf{H}}}_{\mathrm{v}}{\parallel }_1+\frac{1}{2}{\parallel \mathbf{C}\parallel }_F^2\hfill \\ & & +\frac{1}{2}{∥P_{\mathsf{\Omega }}(\mathbf{E}-\widehat{\mathbf{H}})∥}_F^2+tr\left({\mathbf{Z}}_1^{\mathrm{H}}(\widehat{\mathbf{H}}-\mathbf{E})\right)+\frac{t}{2}{\parallel \widehat{\mathbf{H}}-\mathbf{E}\parallel }_{\mathrm{F}}^2\\ & & +tr\left({\mathbf{Z}}_2^H\left(\mathbf{E}-{\mathbf{A}}_{\mathrm{MS}}{\widehat{\mathbf{H}}}_{\mathrm{v}}{\mathbf{A}}_{\mathrm{BS}}^{\mathrm{H}}-\mathbf{C}\right)\right)\\ & & +\frac{t}{2}{∥\mathbf{E}-{\mathbf{A}}_{\mathrm{MS}}{\widehat{\mathbf{H}}}_{\mathrm{v}}{\mathbf{A}}_{\mathrm{BS}}^{\mathrm{H}}-\mathbf{C}∥}_{\mathrm{F}}^2,\end{array}$$

(22)

where ${\mathbf{Z}}_1$ and ${\mathbf{Z}}_2\in {\mathbb{C}}^{N_{\mathrm{MS}}\times N_{\mathrm{BS}}}$ are dual variables (the Lagrange multipliers) and $t>0$ is the step size. The estimator is summarized in Algorithm 5, where:

• $\tau ={\rho }_1∥P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)∥$ with ${\rho }_1=\frac{3N}{N_{\mathrm{BS}}{N}_{\mathrm{MS}}}$ in Step 3;
• ${\mathbf{z}}_i$ denotes the vectorization of ${\mathbf{Z}}_i$, and similarly for other variables;
• $\mathbf{A}\triangleq {\sum }_{i=1}^{N_{\mathrm{MS}}}diag{\left({\left[{\mathsf{\Omega }}_{*}\right]}_i\right)}^T\otimes {\mathbf{I}}_{ii}$ where ${\mathsf{\Omega }}_{*}\in {\{0,1\}}^{N_{\mathrm{MS}}\times N_{\mathrm{BS}}}$ is composed of N ones and $N_{\mathrm{BS}}{N}_{\mathrm{MS}}-N$ zeros, the value 1 indicates the position of a sample from the channel matrix, and ${\left[{\mathsf{\Omega }}_{*}\right]}_i$ denotes the i-th row of ${\mathsf{\Omega }}_{*}$, and $\mathbf{I}$ is the $N_{\mathrm{MS}}\times N_{\mathrm{MS}}$ matrix that the value at its $(i,i)$-th position is 1 and the remaining position is 0 [17];
• The parameters in Equation (20) are chosen empirically as ${\tau }_{\mathrm{L}}=t{∥P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)∥}_2$ and ${\tau }_{\mathrm{S}}=\frac{0.1}{1-10log\left({\sigma }_n^2\right)}$, where ${\sigma }_n^2$ is the noise power.

The computational cost of the ADMM algorithm is mainly due to the SVD in Step 3 and the matrix operations in Steps 4, 5, 6, and 7.

Algorithm 5 ADMM Estimator

Require:  $P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)$, t, ${\tau }_{\mathrm{L}}$, ${\tau }_{\mathrm{S}}$, $J_{\mathrm{ADMM}}$, $\mathbf{A}$, $\mathsf{\Psi }$, $\tau$.

1: $\mathbf{Initialization}$: ${\widehat{\mathbf{H}}}^0={{\widehat{\mathbf{H}}}_{\mathrm{v}}}^0={\mathbf{C}}^0={\mathbf{E}}^0={\mathbf{Z}}_2^0={\mathbf{Z}}_1^0=\mathbf{0}$; their vectorizations ${\widehat{\mathbf{h}}}^0={{\widehat{\mathbf{h}}}_{\mathrm{v}}}^0={\mathbf{c}}^0={\mathbf{e}}^0={\mathbf{z}}_2^0={\mathbf{z}}_1^0=\mathbf{0}$; 2: for$\ell =0,1,\dots ,J_{\mathrm{ADMM}}-1$do 3:     Update ${\widehat{\mathbf{H}}}^{\ell +1}={\mathbf{U}}_{\mathrm{A}}{\mathcal{S}}_{\tau }\left({\mathbf{S}}_{\mathrm{A}}\right){\mathbf{V}}_{\mathrm{A}}^{\mathrm{H}}$, where $\left[{\mathbf{U}}_{\mathrm{A}},{\mathbf{S}}_{\mathrm{A}},{\mathbf{V}}_{\mathrm{A}}\right]={\mathbf{E}}^{\ell }-\frac{1}{t}{\mathbf{Z}}_1^{\ell }$ 4:     Update $$\begin{array}{cc}\hfill {\mathbf{e}}^{\ell +1}=& {\left({\mathbf{A}}^{\mathrm{H}}\mathbf{A}+2t\mathbf{I}\right)}^{-1}\left({\mathbf{z}}_1^{\ell }+t{\widehat{\mathbf{h}}}^{\ell +1}+{\mathbf{A}}^{\mathrm{H}}vec\left(P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)\right)+{\mathbf{z}}_2^{\ell }+t{\mathbf{c}}^{\ell }+t\mathsf{\Psi }{\widehat{\mathbf{h}}}_{\mathrm{v}}^{\ell }\right).\hfill \end{array}$$ 5:     Update $$\begin{array}{ccc}\hfill {\widehat{\mathbf{h}}}_{\mathrm{v}}^{\ell +1}=& sign\left(\mathcal{R}\left({\mathbf{v}}^{\ell +1}\right)\right)\circ max\left(\left|\mathcal{R}\left({\mathbf{v}}^{\ell +1}\right)\right|-{\tau }_{\mathrm{S}}^{\prime },0\right)\hfill & \hfill +sign\left(\mathcal{I}\left({\mathbf{v}}^{\ell +1}\right)\right)\circ max\left(\left|\mathcal{I}\left({\mathbf{v}}^{\ell +1}\right)\right|-{\tau }_{\mathrm{S}}^{\prime },0\right)\end{array}$$ with ${\mathbf{V}}^{\ell +1}\triangleq {\mathbf{A}}_{\mathrm{MS}}^H\left(\frac{1}{t}{\mathbf{Z}}_2^{\ell }-{\mathbf{C}}^{\ell }+{\mathbf{E}}^{\ell +1}\right){\mathbf{A}}_{\mathrm{BS}}$ and ${\tau }_{\mathrm{S}}^{\prime }\triangleq {\tau }_{\mathrm{S}}/t$. 6:     Update $${\mathbf{C}}^{\ell +1}=\frac{t}{t+1}\left({\mathbf{E}}^{\ell +1}-{\mathbf{A}}_{\mathrm{MS}}{\widehat{\mathbf{H}}}_v^{\ell +1}{\mathbf{A}}_{\mathrm{BS}}^{\mathrm{H}}+\frac{1}{t}{\mathbf{Z}}_2^{\ell }\right)$$ 7:     Update $${\mathbf{Z}}_1^{\ell +1}={\mathbf{Z}}_1^{\ell }+t\left(P_{\mathsf{\Omega }}({\mathbf{E}}^{\ell +1}-\tilde{\mathbf{H}})\right)$$ $${\mathbf{Z}}_2^{\ell +1}={\mathbf{Z}}_2^{\ell }+t\left({\mathbf{E}}^{\ell +1}-{\mathbf{A}}_{\mathrm{MS}}{\widehat{\mathbf{H}}}_v^{\ell +1}{\mathbf{A}}_{\mathrm{BS}}^{\mathrm{H}}-{\mathbf{C}}^{\ell +1}\right)$$ 8: end for; 9: return$\widehat{\mathbf{H}}={\widehat{\mathbf{H}}}^{J_{\mathrm{ADMM}}-1}$

4.2.2. Two-Stage Estimator

The two-stage estimator [16] also exploits the sparsity and low-rankness of the channel matrix. In the first stage, MC is applied to provide a denoised channel estimation based on the low rankness of the channel matrix and then the second stage employs CS to refine the estimation based on the array response and the virtual representation of the channel matrix.

The SVT algorithm [25] introduced above solves the low-rank matrix completion problem:

$$\underset{\widehat{\mathbf{H}}}{min}\mathrm{rank}\left(\widehat{\mathbf{H}}\right),\mathrm{s}.\mathrm{t}.\left|\right|P_{\mathsf{\Omega }}\left(\widehat{\mathbf{H}}\right)-P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right){\left|\right|}_F^2\le {\delta }_n^2.$$

(23)

The fast iterative shrinkage threshold algorithm (FISTA) proposed by [37] is used to solve the sparse vector recovery problem:

$$\underset{\mathrm{vec}\left({\widehat{\mathbf{H}}}_{\mathrm{v}}\right)}{min}\parallel \mathsf{\Psi }\mathrm{vec}\left({\widehat{\mathbf{H}}}_{\mathrm{v}}\right)-\mathrm{vec}\left(\widehat{\mathbf{H}}\right){\parallel }_{\mathrm{F}}^2+\lambda {\parallel \mathrm{vec}\left({\widehat{\mathbf{H}}}_{\mathrm{v}}\right)\parallel }_1.$$

(24)

The estimator is summarized in Algorithm 6, where:

• The parameters of the first stage, SVT, is the same with Algorithm 1;
• $\mathsf{\Psi }={\mathbf{A}}_{\mathrm{BS}}^{*}\otimes {\mathbf{A}}_{\mathrm{MS}}$;
• $\lambda$ is a constant stepsize, e.g., $\lambda =0.001$;
• ${\lambda }_{max}$ is the top eigenvalue of ${\mathsf{\Psi }}^{\mathrm{H}}\mathsf{\Psi }$.

The complexity of SVT in the first stage has been anlyzed. The cost of the FISTA algorithm mainly consists of applying the sensing matrix in Step 4, which has a complexity of $O\left(N_{\mathrm{BS}}^2{N}_{\mathrm{MS}}^2\right)$.

Algorithm 6 Two-Stage Estimator

Require: $P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)$

1: Use SVT to recover $\widehat{\mathbf{H}}$ as $$\begin{array}{cc}\hfill & \underset{\widehat{\mathbf{H}}}{min}\parallel \widehat{\mathbf{H}}{\parallel }_{*}\mathrm{s}.\mathrm{t}.\left|\right|P_{\mathsf{\Omega }}\left(\widehat{\mathbf{H}}\right)-P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right){\left|\right|}_F^2\le {\delta }_n^2\hfill \end{array}$$ Require: $\mathsf{\Psi },\mathrm{vec}\left(\widehat{\mathbf{H}}\right),\lambda ,J_{\mathrm{FISTA}}$, the top eigenvalue ${\lambda }_{max}$ of ${\mathsf{\Psi }}^{\mathrm{H}}\mathsf{\Psi }$ 2: Initialize ${\mathbf{y}}_1={\mathbf{x}}_1=\mathbf{0}$, $t_1=1$ 3: for$i=1:J_{\mathrm{FISTA}}$do 4:     ${\mathbf{c}}_i={\mathbf{y}}_i-(1/{\lambda }_{max}){\mathsf{\Psi }}^{\mathrm{H}}(\mathsf{\Psi }{\mathbf{y}}_i-\mathrm{vec}\left(\widehat{\mathbf{H}}\right))$ 5:     ${\mathbf{x}}_{i+1}=(max(abs\left({\mathbf{c}}_i\right)-2\lambda {\lambda }_{max},0))sign\left({\mathbf{c}}_i\right)$ 6:     $t_{i+1}=\frac{1+\sqrt{1+4t_i^2}}{2}$ 7:     ${\mathbf{y}}_{i+1}={\mathbf{x}}_{i+1}+\left(\frac{t_i-1}{t_{i+1}}\right)\left({\mathbf{x}}_{i+1}-{\mathbf{x}}_i\right)$ 8:     ${\mathbf{x}}_i\leftarrow {\mathbf{x}}_{i+1}$ 9:     ${ϵ}_{i+1}={\parallel \mathsf{\Psi }{\mathbf{x}}_i-\mathrm{vec}\left(\widehat{\mathbf{H}}\right)\parallel }_2$ 10:     if $abs({ϵ}_{i+1}-{ϵ}_i)\le {10}^{-6}$ then 11:         Break 12:     end if 13: end for 14: return$\widehat{\mathbf{H}}={\mathbf{A}}_{\mathrm{MS}}{vec}^{-1}\left({\mathbf{x}}_i\right){\mathbf{A}}_{\mathrm{BS}}^{\mathrm{H}}$

5. Numerical Results

Consider switch-based MIMO systems over a mmWave channel at 90 GHz. When not otherwise specified, the number of clusters $C\sim max\left(\mathrm{Poisson}\right(1.8),1)$, and the number of rays $R\sim \mathcal{U}\left[1,20\right]$ (the total number of paths is $L=CR$); the AoDs and AoAs follow Laplace distributions with a standard deviation of ${15}^{\circ }$ [38]; the sub-power of the clusters ${\gamma }_{\mathrm{c}}=1$; and the ULA at the MS and BS has $N_{\mathrm{MS}}=32$ antennas and $N_{\mathrm{BS}}=128$ antennas, respectively, and $N_{\mathrm{RF}}=4$ RF chains. Here, the CS method based on OMP is compared with MC-based estimators. The parameter settings are as follows:

• OMP: The unitary dictionary is set with $N_1=N_{\mathrm{MS}}=32$ and $N_2=N_{\mathrm{BS}}=128$. The stopping threshold is set as ${ϵ}_{\mathrm{O}\mathrm{M}\mathrm{P}}=0.025{\sigma }^2,0.05{\sigma }^2,0.075{\sigma }^2,0.1{\sigma }^2,0.125{\sigma }^2,0.15{\sigma }^2$ with $\mathrm{PNR}=0,5,10,15,20,25\mathrm{dB}$ [20];
• SVT: $\delta =3.2$, $ϵ={10}^{-4}$, $\tau =5\sqrt{N_{\mathrm{MS}}{N}_{\mathrm{BS}}}$, $J_{\mathrm{SVT}}=100$, and $k_0=5$;
• SVP: $\eta =0.5$ and $ϵ={10}^{-4}$;
• FPC: ${\mu }_{\mathrm{final}}=0.01$, $J_{\mathrm{FPC}}=100$, $ϵ={10}^{-4}$, ${\eta }_{\mu }=0.25$, and $\delta =1.99$;
• GCG-Alt: $\mu ={\sigma }_{\mathrm{n}}^2$, $ϵ=0.01$, and ${ϵ}_{\mathrm{a}}=0.1$;
• ADMM: ${\tau }_{\mathrm{L}}=t{\parallel P_{\mathsf{\Omega }}\left(\tilde{\mathbf{H}}\right)\parallel }_2$, $t=0.005$, ${\tau }_{\mathrm{S}}=\frac{0.1}{1-10log\left({\sigma }_{\mathrm{n}}^2\right)}$, $\delta =3.2$;
• Two-Stage: $J_{\mathrm{FISTA}}=100$, $\lambda =0.001$.

The normalized MSE (NMSE) is defined as:

$$\mathrm{NMSE}=\frac{\parallel \widehat{\mathbf{H}}{-\mathbf{H}\parallel }_F^2}{{\parallel \mathbf{H}\parallel }_{\mathrm{F}}^2}$$

where $\widehat{\mathbf{H}}$ is an estimate of the channel matrix $\mathbf{H}$.

5.1. Comparison of NMSE When There Are No Hardware Impairments

Figure 2 compares different estimators in terms of NMSE with different PNRs, which is defined as $\mathrm{PNR}=10log10\left(\frac{\rho }{{\sigma }_{\mathrm{n}}^2}\right)$, when $N_{\mathrm{MS}}=32$, $N_{\mathrm{BS}}=128$, and $N_{\mathrm{RF}}=4$. The SVP algorithm performed worse than others. The ADMM algorithm performed the best, and the NMSE of the ADMM estimator was at least 5 dB better than other estimators. The NMSE with $N_{\mathrm{BS}}=64$, $N_{\mathrm{MS}}=64$ and $N_{\mathrm{RF}}=8$ is shown in Figure 3. From Figure 2 and Figure 3, the relative performance of the different estimators is similar with different numbers of antennas.

We next show the performance with different number of channel paths. Here, we assume the number of clusters $C=1$, the numbers of rays R is changed from 1 to 22, and so the total number of channel paths $L=R$. From Figure 4, as the numbers of paths L increases, the NMSE of different algorithms degrades. The SVP algorithm is the most sensitive to the number of paths while the OMP estimator is the most robust. The sampling ratio is also critical for performance. From Figure 5, the NMSE improves substantially with increasing sampling ratio. It is noted that CS-based schemes such as the OMP and the two-stage algorithms are more advantageous when the sampling ratio is low.

5.2. NMSE Comparison When There Are Hardware Impairments

In practice, it is inevitable to have impairments of the antenna elements, which are typically time-varying, e.g., due to temperature changes or hardware aging [39]. Therefore, the array response may be severely impacted. Due to mechanical reasons and uncertainty regarding the precise position of the antenna phase center, the actual antenna position may deviate from the assumed ideal array shape. Following [20], we define the gain and phase error vector at the BS as:

$${\mathbf{e}}_{\mathrm{BS}}={\left[{\beta }_1{\mathrm{e}}^{j{\omega }_1},{\beta }_2{\mathrm{e}}^{j{\omega }_2},\cdots ,{\beta }_{N_{\mathrm{BS}}}{\mathrm{e}}^{j{\omega }_{N_{\mathrm{BS}}}}\right]}^{\mathrm{T}}$$

(25)

where ${\omega }_i$ represents the phase errors and ${\beta }_i$ denotes the amplitude gain of each antenna element. The gain and phase error vector ${\mathbf{e}}_{\mathrm{MS}}$ at the MS is similar to ${\mathbf{e}}_{\mathrm{BS}}$. Considering such errors, the received signal in Equation (2) can be expressed as:

$$\tilde{\mathbf{y}}={\mathbf{W}}^{\mathrm{H}}{\mathbf{E}}_{\mathrm{MS}}\mathbf{H}{\mathbf{E}}_{\mathrm{BS}}^{\mathrm{H}}\mathbf{F}s+{\mathbf{W}}^{\mathrm{H}}{\mathbf{E}}_{\mathrm{MS}}\widehat{\mathbf{n}}$$

(26)

where ${\mathbf{E}}_{\mathrm{MS}}$ is a diagonal matrix with ${\mathbf{e}}_{\mathrm{MS}}$ as the diagonal elements, and ${\mathbf{E}}_{\mathrm{BS}}$ is defined similarly.

We carry out simulations to examine the robustness of the different approaches when there are phase and gain errors in the array response. Those errors are assumed to be uniformly distributed within certain range and characterized by the level of phase and gain errors, respectively, following [20]. It is found that MC estimators were not affected by the phase or gain errors since the estimators were independent of the array response vectors. The sparsity-based methods use array response which depends on the phase and gain information, and can result in poor channel estimation when unknown gain and phase errors are present. These were validated by simulation experiments. From Figure 6 and Figure 7 where NMSE achieved with different levels of phase and gain errors are compared, the MC estimators based on SVP, FPC, GCG-Alt, and SVT, were insensitive to phase errors or gain errors. The estimators exploiting CS, such as the one based on OMP and hybrid methods (the ADMM and two-stage estimators), were more sensitive to phase errors or gain errors.

5.3. Computational Complexity

Finally, we compare the computational complexity of different algorithms in Figure 8 and Figure 9. It can be seen that in general, the hybrid estimators (the ADMM and two-stage estimators) exhibited higher complexity because they both involved SVD and the application of the sensing matrices were of large sizes. In particular, the ADMM algorithm had the highest complexity. The MC estimators had moderate complexity that did not vary significantly with the PNRs and the FPC exhibited higher complexity than SVT and SVP as more SVD operations were required. The GCG-Alt algorithm had the lowest complexity. This was because the GCG-Alt algorithm had a fast convergence rate and also avoided the SVD operations used in the other MC algorithms such as SVP, SVT, and FPC. The computational complexity of the OMP estimator increased with the PNR because the numbers of paths recovered by the OMP increased as the PNR increased.

5.4. Performance with Line of Sight (LoS) Propagation

In the above, we have focused on the comparison for the channels where the different paths had the same average power. In practical applications, there may exist line of sight (LoS) propagation, where a path contributes a significant portion of the power gain [40,41]. In this case, the channel model can be modified from Equation (3) to:

$${\displaystyle \mathbf{H}=\frac{1}{\sqrt{R}}\sum _{c=1}^C\sum _{r=1}^R{g}_{cr}{\mathbf{a}}_{\mathrm{MS}}\left({\varphi }_{cr}^{\mathrm{MS}}\right){\mathbf{a}}_{\mathrm{BS}}^{\mathrm{H}}\left({\varphi }_{cr}^{\mathrm{BS}}\right)+\beta {\mathbf{a}}_{\mathrm{MS}}\left({\varphi }_{\mathrm{LoS}}^{\mathrm{MS}}\right){\mathbf{a}}_{\mathrm{BS}}^{\mathrm{H}}\left({\varphi }_{\mathrm{LoS}}^{\mathrm{BS}}\right).}$$

(27)

In our simulations we set $g_{cr}\sim \mathcal{CN}(0,0.5)$ and $\beta =\sqrt{C/2}$. The NMSE results are shown in Figure 10. With a LoS path which dominates the power gain, the effective rank of the channel matrix may be reduced. The various channel estimators considered are more effective in this case, leading to lower NMSE as compared with the case where the LoS path is absent. The results also demonstrate the effectiveness of the MC-based methods in realistic scenarios.

6. Conclusions

In this paper, we compared the performance of several MC-based channel estimators for mmWave massive MIMO systems. It was observed that the hybrid ADMM algorithm exhibited the best performance in general, which jointly exploited the low rank property of the channel matrix and the sparsity in the angular domain. However, it also exhibited the highest complexity among the estimators compared. The MC-based estimators (using GCG-Alt, SVT, SVP, or FPC) were robust against array impairments as they did not rely on array response vectors. Among them the GCG-Alt estimator exhibited the lowest complexity, better performance, and provided a competitive solution when the arrays were not perfectly calibrated.

In this work, we considered a point-to-point mmWave system. The estimators could also be applied to multiuser systems when orthogonal training schemes are deployed. The comparison of these methods when nonorthogonal training is used would be an interesting study for future work.