Sparse Subspace Decomposition for Millimeter Wave MIMO Channel Estimation

TL;DR

This paper introduces a sparse subspace decomposition method for millimeter wave MIMO channel estimation, effectively addressing the challenges of large array sizes and hybrid architectures with low overhead.

Contribution

The paper proposes a novel SSD technique formulated as an optimization problem, improving channel estimation accuracy for large-scale millimeter wave MIMO systems.

Findings

Outperforms existing channel estimation techniques

Achieves low overhead in large array scenarios

Effective for hybrid precoding architectures

Abstract

Millimeter wave multiple-input multiple-output (MIMO) communication systems must operate over sparse wireless links and will require large antenna arrays to provide high throughput. To achieve sufficient array gains, these systems must learn and adapt to the channel state conditions. However, conventional MIMO channel estimation can not be directly extended to millimeter wave due to the constraints on cost-effective millimeter wave operation imposed on the number of available RF chains. Sparse subspace scanning techniques that search for the best subspace sample from the sounded subspace samples have been investigated for channel estimation.However, the performance of these techniques starts to deteriorate as the array size grows, especially for the hybrid precoding architecture. The millimeter wave channel estimation challenge still remains and should be properly addressed before the system can be deployed and used to its full potential. In this work, we propose a sparse subspace decomposition (SSD) technique for sparse millimeter wave MIMO channel estimation. We formulate the channel estimation as an optimization problem that minimizes the subspace distance from the received subspace samples. Alternating optimization techniques are devised to tractably handle the non-convex problem. Numerical simulations demonstrate that the proposed method outperforms other existing techniques with remarkably low overhead.

Tables

Table 1. TABLE I: Execution Time (in seconds) for Algorithms.

$N_r\times N_t$	Channel
Uses ($K$)	MF	NNM	SSD	SSD-T
$16\times 64$	100	5.23	228.83	15.09	2.41
$16\times 64$	150	10.60	483.09	25.05	3.78
$8\times 16$	20	0.32	3.47	1.19	0.17
$8\times 16$	60	1.20	12.60	3.91	0.57

Equations

$$\mathbf{r}={\mathbf{W}}_{\mathrm{D}}^{\ast }{\mathbf{W}}_{\mathrm{A}}^{\ast }\mathbf{H}{\mathbf{F}}_{\mathrm{A}}{\mathbf{F}}_{\mathrm{D}}\mathbf{s}+{\mathbf{W}}_{\mathrm{D}}^{\ast }{\mathbf{W}}_{\mathrm{A}}^{\ast }\mathbf{n},$$

$${\mathbf{y}}_k={\mathbf{W}}_k^{\ast }\mathbf{H}{\mathbf{F}}_k{\mathbf{s}}_k+{\mathbf{W}}_k^{\ast }{\mathbf{n}}_k,k=1,\dots ,K,$$

$$k_{opt}={\mathrm{argmax}}_{k=1,\dots ,K}{\Vert {\mathbf{y}}_k\Vert }_2^2.$$

$${\mathbf{y}}_k={\mathbf{W}}_k^{\ast }\mathbf{H}{\mathbf{f}}_k+{\mathbf{W}}_k^{\ast }{\mathbf{n}}_k,k=1,2,\dots ,K,$$

$$K>d(N_{\mathrm{r}}+N_{\mathrm{t}}-d)/N_{\mathrm{RF}}.$$

$$\displaystyle\!\!\big{(}\widehat{{\mathbf{W}}},\widehat{{\mathbf{F}}}\big{)}\!\!=\!\mathop{\mathrm{argmax}}\limits_{{\mathbf{W}},{\mathbf{F}}}\!\left\|{\mathbf{W}}^{*}{\mathbf{H}}{\mathbf{F}}\right\|_{F}^{2}\ \!\text{s.t.}\ {\mathbf{W}}^{*}{\mathbf{W}}\!\!=\!{\mathbf{I}}_{N_{\mathrm{d}}},{\mathbf{F}}^{*}{\mathbf{F}}\!\!=\!{\mathbf{I}}_{N_{\mathrm{d}}},\!\!\!\!$$

$$\displaystyle\begin{aligned} \big{(}\widehat{{\mathbf{U}}},\widehat{{\mathbf{V}}}\big{)}=&\mathop{\mathrm{argmin}}_{{\mathbf{U}},{\mathbf{V}}}\min_{\mathbf{\Lambda}}\left\|{\mathbf{H}}-{\mathbf{U}}\mathbf{\Lambda}{\mathbf{V}}^{*}\right\|_{F}^{2}\\ &\text{~{}subject to}~{}~{}~{}{\mathbf{U}}^{*}{\mathbf{U}}={\mathbf{I}}_{d},~{}{\mathbf{V}}^{*}{\mathbf{V}}={\mathbf{I}}_{d},\end{aligned}$$

$$\displaystyle\begin{aligned} \big{(}\widehat{{\mathbf{U}}},\widehat{{\mathbf{V}}},\widehat{\mathbf{\Lambda}}\big{)}=&\mathop{\mathrm{argmin}}_{{\mathbf{U}},{\mathbf{V}},\mathbf{\Lambda}}\sum_{k=1}^{K}\left\|{\mathbf{y}}_{k}-{\mathbf{W}}_{k}^{*}{\mathbf{U}}\mathbf{\Lambda}{\mathbf{V}}^{*}{\mathbf{f}}_{k}\right\|^{2}_{2}\\ &\text{~{}subject to}~{}~{}~{}{\mathbf{U}}^{*}{\mathbf{U}}={\mathbf{I}}_{d},~{}{\mathbf{V}}^{*}{\mathbf{V}}={\mathbf{I}}_{d}.\end{aligned}$$

$$\mathcal{A}(\mathcal{F},\mathcal{W},\mathbf{H})={\left[({\mathbf{W}}_1^{\ast }\mathbf{H}{\mathbf{f}}_1{)}^T,\dots ,({\mathbf{W}}_K^{\ast }\mathbf{H}{\mathbf{f}}_K{)}^T\right]}^T\in {\mathbb{C}}^{K{N}_{\mathrm{RF}}\times 1}.$$

$$\displaystyle\begin{aligned} \big{(}\widehat{{\mathbf{U}}},\widehat{{\mathbf{V}}},\widehat{\mathbf{\Lambda}}\big{)}=&\mathop{\mathrm{argmin}}_{{\mathbf{U}},{\mathbf{V}},\mathbf{\Lambda}}\left\|{\mathbf{y}}-{\mathcal{A}}({\mathcal{F}},{\mathcal{W}},{\mathbf{U}}\mathbf{\Lambda}{\mathbf{V}}^{*})\right\|^{2}_{2}\\ &\text{~{}subject to}~{}~{}~{}{\mathbf{U}}^{*}{\mathbf{U}}={\mathbf{I}}_{d},~{}{\mathbf{V}}^{*}{\mathbf{V}}={\mathbf{I}}_{d}.\end{aligned}$$

$$\begin{array}{r}\hat{\mathbf{H}}=\underset{\mathbf{H}}{\mathrm{argmin}}{\Vert \mathbf{H}\Vert }_{\ast }\text{subject to }\mathbf{y}-\mathcal{A}(\mathcal{F},\mathcal{W},\mathbf{H})\in \mathcal{C}.\end{array}$$

$$\displaystyle\begin{aligned} \big{(}\widehat{{\mathbf{U}}},\widehat{{\mathbf{V}}}\big{)}=\mathop{\mathrm{argmin}}\limits_{{\mathbf{U}},{\mathbf{V}}}\left\|{\mathbf{y}}-{\mathcal{A}}({\mathcal{F}},{\mathcal{W}},{\mathbf{U}}{\mathbf{V}})\right\|_{2}^{2},\end{aligned}$$

$$\displaystyle\begin{aligned} \big{(}\widehat{{\mathbf{U}}},\widehat{{\mathbf{V}}},\widehat{\mathbf{\Lambda}}\big{)}=&\mathop{\mathrm{argmin}}\limits_{{\mathbf{U}},\mathbf{\Lambda},{\mathbf{V}}}\left\|{\mathbf{y}}-{\mathcal{A}}({\mathcal{F}},{\mathcal{W}},{\mathbf{U}}\mathbf{\Lambda}{\mathbf{V}}^{*})\right\|_{2}^{2}\\ &\text{subject to~{}~{}}\mathop{\mathrm{tr}}({\mathbf{U}}^{*}{\mathbf{U}})\leq d,~{}\mathop{\mathrm{tr}}({\mathbf{V}}^{*}{\mathbf{V}})\leq d.\end{aligned}$$

$$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\text{(S1)}\begin{cases}{\mathbf{U}}_{(\ell+1)}\!=\!&\!\!\!\!\mathop{\mathrm{argmin}}\limits_{{\mathbf{U}}}\big{\|}{\mathbf{y}}\!-\!{\mathcal{A}}({\mathcal{F}},{\mathcal{W}},{\mathbf{U}}\mathbf{\Lambda}_{(\ell)}{\mathbf{V}}_{(\ell)}^{*})\big{\|}_{2}^{2}\\ &\!\!\!\!\!\!\text{subject to~{}}\mathop{\mathrm{tr}}({\mathbf{U}}^{*}{\mathbf{U}})\leq d,\end{cases}$$

$$\text{(S2)}\{\begin{array}{ll}{\mathbf{V}}_{(\ell +1)}=& \underset{\mathbf{V}}{\mathrm{argmin}}{\Vert \mathbf{y}-\mathcal{A}(\mathcal{F},\mathcal{W},{\mathbf{U}}_{(\ell +1)}{\mathbf{\Lambda }}_{(\ell )}{\mathbf{V}}^{\ast })\Vert }_2^2\\ & \text{subject to }\mathrm{tr}({\mathbf{V}}^{\ast }\mathbf{V})\le d,\end{array}$$

$$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\text{(S3)~{}}\!\mathbf{\Lambda}_{(\ell+1)}\!=\!\mathop{\mathrm{argmin}}\limits_{\mathbf{\Lambda}}\!\big{\|}{\mathbf{y}}\!-\!{\mathcal{A}}({\mathcal{F}},{\mathcal{W}},{\mathbf{U}}_{(\ell+1)}\mathbf{\Lambda}{\mathbf{V}}_{(\ell+1)}^{*})\big{\|}_{2}^{2},$$

$$\begin{array}{r}\hat{\mathbf{U}}=\underset{\mathbf{U}}{\mathrm{argmin}}{\Vert \mathbf{y}-\mathcal{A}(\mathcal{F},\mathcal{W},\mathbf{U}\mathbf{\Lambda }{\mathbf{V}}^{\ast })\Vert }_2^2\text{s.t. }\mathrm{tr}({\mathbf{U}}^{\ast }\mathbf{U})\le d,\end{array}$$

$$\mathbf{B}=[((\mathbf{\Lambda }{\mathbf{V}}^{\ast }{\mathbf{f}}_1{)}^T\otimes {\mathbf{W}}_1^{\ast }{)}^T,\dots ,((\mathbf{\Lambda }{\mathbf{V}}^{\ast }{\mathbf{f}}_K{)}^T\otimes {\mathbf{W}}_K^{\ast }{)}^T{]}^T$$

$$\hat{\mathbf{U}}=\underset{\mathbf{U}}{\mathrm{argmin}}{\Vert \mathbf{y}-\mathbf{B}\mathrm{vec}(\mathbf{U})\Vert }_2^2\text{subject to }{\Vert \mathrm{vec}(\mathbf{U})\Vert }_2^2\le d,$$

$$\{\begin{array}{l}{\mathbf{B}}^{\ast }\mathbf{B}\mathrm{vec}(\mathbf{U})-{\mathbf{B}}^{\ast }\mathbf{y}+\lambda \mathrm{vec}(\mathbf{U})=0,\text{(first-order condition)}\\ \lambda \ge 0,\text{(dual constraint)}\\ \lambda [{\Vert \mathrm{vec}(\mathbf{U})\Vert }_2^2-d]=0,\text{(complementary slackness)}\text{par}\end{array}$$

$$\hat{\mathbf{\Lambda }}=\underset{\mathbf{\Lambda }}{\mathrm{argmin}}{\Vert \mathbf{y}-\mathcal{A}(\mathcal{F},\mathcal{W},\mathbf{U}\mathbf{\Lambda }{\mathbf{V}}^{\ast })\Vert }_2^2,$$

$$\mathrm{diag}(\hat{\mathbf{\Lambda }})=({\mathbf{P}}^{\ast }\mathbf{P}{)}^{-1}{\mathbf{P}}^{\ast }\mathbf{y},$$

$$\mathcal{A}(\mathcal{F},\mathcal{W},\mathbf{U}\mathbf{\Lambda }{\mathbf{V}}^{\ast })$$

$$\displaystyle\textstyle\widehat{\mathbf{\Lambda}}=\mathop{\mathrm{argmin}}\limits_{\mathbf{\Lambda}}\big{\|}{\mathbf{y}}-\sum\limits_{j=1}^{d}[\mathbf{\Lambda}]_{j,j}{\mathcal{A}}({\mathcal{F}},{\mathcal{W}},[{\mathbf{U}}]_{j}[{\mathbf{V}}]_{j}^{*})\big{\|}_{2}^{2}.$$

$$\mathrm{diag}(\hat{\mathbf{\Lambda }})=({\mathbf{P}}^{\ast }\mathbf{P}{)}^{-1}{\mathbf{P}}^{\ast }\mathbf{y},$$

$$\displaystyle\textstyle\text{NMSE}=\mathbb{E}\Big{[}{\big{\|}{\mathbf{H}}-\widehat{{\mathbf{H}}}\big{\|}^{2}_{F}}/{\left\|{\mathbf{H}}\right\|^{2}_{F}}\Big{]}.$$

$$\displaystyle\textstyle{\big{\|}{\mathbf{y}}-{\mathcal{A}}({\mathcal{F}},{\mathcal{W}},\widehat{{\mathbf{H}}}_{(\ell)})\big{\|}_{2}^{2}}\big{/}(K\!N_{\mathrm{RF}})<\sigma^{2}.$$

$$\mathbf{H}=\sqrt{\frac{N_{\mathrm{r}}{N}_{\mathrm{t}}}{L}}\sum _{l=1}^L{\alpha }_l{\mathbf{a}}_{\mathrm{r}}({\varphi }_{\mathrm{r},l}){\mathbf{a}}_{\mathrm{t}}^{\ast }({\varphi }_{\mathrm{t},l}).$$

Full text

I Introduction

The operating frequency of modern wireless systems is steadily shifting to the millimeter wave (e.g.,​ 28-100 GHz) bands which can provide a much wider bandwidth. These higher frequencies will force systems to use large arrays to generate narrow beams in order to overcome the pathloss and atmospheric impairments encountered at these frequencies. There is particular interest in using advanced multiple antenna concepts at millimeter wave frequencies, but these techniques cannot be directly applied due to the prohibitively complex baseband signal processing overhead​ [1].

To alleviate the increased overhead, hybrid precoding architectures [2, 3] have been explored for use at millimeter wave frequencies [4, 5]. In a hybrid architecture, a large antenna-array is driven by a limited number of RF chains mapped to the antennas using analog processing, i.e., a phase shifter network that linearly processes the RF signal. Digital processing of the small number of RF chains is then possible. This architecture, however, does not allow the digital baseband to directly access the entire channel, since it only processes analog pre-processed samples. This causes a subspace sampling limitation problem [6], which limits the hybrid system’s ability to track and adapt to the channel.

To address the subspace sampling limitation, earlier work mainly focused on hierarchical beam scanning techniques [6, 7, 8]. The transmitter and receiver iteratively search for the best subspace pair by scanning all the subspace pairs chosen in the hierarchical subspace sounding codebooks [6, 7]. However, the performances of these techniques will be degraded as the antenna size grows, especially when combined with hybrid precoding, due to the increased search overhead.

There have been multiple channel measurement campaigns demonstrating that practical millimeter wave channels present substantially sparse scattering [9, 10]. In this kind of environment, the size of the useful subspace is considerably smaller than the channel dimension. Hence, channel estimation should be adapted to the sparse channel subspace in order to address the subspace sampling limitation problem. Recently, subspace estimation using the Arnoldi method [11] and the support recovery leveraging the virtual channel representation [12] have been proposed. With the exception of the work in [11, 12], little is known on how to efficiently adapt channel estimation to the sparse but large-dimensional millimeter wave multiple-input multiple-output (MIMO) channel.

In this paper, we formulate millimeter wave MIMO channel estimation as a subspace distance minimization problem. This problem is related to the general rank minimization problem [13], which is NP-hard. The rank minimization problem [13] can be relaxed by using nuclear norm minimization (NNM) [14, 15]. Though NNM provides robust performance, it does not guarantee a low rank solution and does not scale to large-dimensional applications due to computational issues. A matrix factorization (MF) algorithm is proposed in [16], which is faster than NNM, to ensure a low rank solution. The general challenge of the techniques in [14, 15, 16] is that when noisy observations are made, the reconstruction is often matched to the noise, causing the critical over-matching problem.

In this paper, we propose a sparse subspace decomposition technique for the large-dimensional millimeter wave MIMO channel estimation. The devised technique is general enough to allow channel estimation in different classes of sparse channels. We design an alternating optimization technique that tractably handles the formulated sparse decomposition problem. The algorithm has low complexity and converges to a stationary point. Moreover, based on intuition from the strong law of large numbers, we develop a simple method to avoid the over-matching issue arising in the low signal-to-noise ratio (SNR) regime. The proposed algorithm outperforms the conventional methods [15, 16] with substantially reduced computational overhead.

Notations

A bold lower case letter ${\mathbf{a}}$ is a vector, a bold capital letter ${\mathbf{A}}$ is a matrix. $[{\mathbf{A}}]_{i,j}$ is the $i$th row and $j$th column entry of ${\mathbf{A}}$, and $[{\mathbf{A}}]_{i}$ is the $i$th column of ${\mathbf{A}}$. ${\mathbf{A}}^{*},{\mathbf{A}}^{T},{\mathbf{A}}^{\!-1}$, tr(${\mathbf{A}}$), $\|{\mathbf{A}}\|_{F}$, $\|{\mathbf{a}}\|_{2}$, $\|{\mathbf{A}}\|_{*}$, respectively, are the conjuagate transpose, transpose, inverse, trace of ${\mathbf{A}}$, Frobenius norm of ${\mathbf{A}}$, $\l_{2}$-norm of ${\mathbf{a}}$, and the nuclear norm of ${\mathbf{A}}$, which is the sum of its singular values. ${\mathbf{I}}_{M}\!\in\!{\mathbb{R}}^{M\times M}$ is the identity matrix. ${\mathbf{A}}\otimes{\mathbf{B}}$ is the Kronecker product of ${\mathbf{A}}$ and ${\mathbf{B}}$. Let $\mathop{\mathrm{vec}}({\mathbf{A}})$ be the operator that stacks columns of ${\mathbf{A}}$ into a column vector. Let $\mathop{\mathrm{diag}}({\mathbf{A}})$ be the operator that collects the diagonal elements of a square matrix ${\mathbf{A}}$ and forms a column vector.

II System Model

In this section, we provide signal model, motivations, and general statements of the channel estimation problem in the millimeter wave MIMO systems.

II-A Signal Model

Consider a point-to-point MIMO hybrid precoding system with $N_{\mathrm{t}}$ transmit and $N_{\mathrm{r}}$ receive antennas, where each side is equipped with $N_{\mathrm{RF}}$ RF chains. The number of data streams is $N_{\mathrm{d}}$, and we assume $N_{\mathrm{d}}\!\leq\!N_{\mathrm{RF}}\!\leq\!\min(N_{\mathrm{r}},N_{\mathrm{t}})$. The transmitter employs the analog precoder ${\mathbf{F}}_{\mathrm{A}}\!\in\!{\mathbb{C}}^{N_{\mathrm{t}}\times N_{\mathrm{RF}}}$ and digital precoder ${\mathbf{F}}_{\mathrm{D}}\!\in\!{\mathbb{C}}^{N_{\mathrm{RF}}\times N_{\mathrm{d}}}$. At the receiver, ${\mathbf{W}}_{\mathrm{D}}\!\in\!{\mathbb{C}}^{N_{\text{RF}}\times N_{\mathrm{d}}}$ and ${\mathbf{W}}_{\mathrm{A}}\!\in\!{\mathbb{C}}^{N_{\mathrm{r}}\times N_{\mathrm{RF}}}$, respectively, denote the digital and analog combiners. The entries of ${\mathbf{F}}_{\mathrm{A}}$ and ${\mathbf{W}}_{\mathrm{A}}$ are constrained such that $|[{\mathbf{F}}_{\mathrm{A}}]_{i,j}|\!\!=\!\!\frac{1}{\sqrt{N_{\mathrm{t}}}}$ and $|[{\mathbf{W}}_{\mathrm{A}}]_{i,j}|\!\!=\!\!\frac{1}{\sqrt{N_{\mathrm{r}}}}$, $\forall i,j$, which are imposed due to analog processing. The channel input and output relation is therefore given by

[TABLE]

where ${\mathbf{H}}\in{\mathbb{C}}^{N_{\mathrm{r}}\times N_{\mathrm{t}}}$ is the sparse millimeter wave MIMO channel with $\mathop{\mathrm{rank}}({\mathbf{H}})=L$. We assume $\frac{L}{\min(N_{\mathrm{r}},N_{\mathrm{t}})}\approx 0$. Notice that, in general, $L$ is not known a priori. Instead, we assume that a rank upper bound of ${\mathbf{H}}$, which we define as $d$ (i.e., $\mathop{\mathrm{rank}}({{\mathbf{H}}})\!\leq\!d$), is known a priori. The vector ${\mathbf{n}}\in{\mathbb{C}}^{N_{\mathrm{r}}\times 1}$ in (1) is the additive noise with each entry independent and identically distributed (i.i.d.) as ${\mathbf{n}}\!\sim\!{\mathcal{C}}{\mathcal{N}}(\mathbf{0},\sigma^{2}{\mathbf{I}}_{N_{\mathrm{r}}})$. The vector ${\mathbf{s}}\in{\mathbb{C}}^{N_{\mathrm{d}}\times 1}$ is the transmit signal satisfying $\mathbb{E}[\left\|{\mathbf{F}}_{\mathrm{A}}{\mathbf{F}}_{\mathrm{D}}{\mathbf{s}}\right\|_{2}^{2}]\!=\!1$. The SNR is, thus, $1/\sigma^{2}$.

Extracting a high quality channel estimate is the key to facilitating advanced MIMO precoding techniques. However, conventional MIMO channel estimation does not directly extend to a millimeter wave hybrid precoding architecture since the digital baseband only accesses the analog compressed channel ${\mathbf{W}}_{\mathrm{A}}^{*}{\mathbf{H}}{\mathbf{F}}_{\mathrm{A}}\!\in\!{\mathbb{C}}^{N_{\mathrm{RF}}\times N_{\mathrm{RF}}}$, rather than the entire ${\mathbf{H}}$ as shown in (1). This limitation is the major hurdle when estimating the large-dimensional millimeter wave MIMO channel. Without direct access to ${\mathbf{H}}$, the channel sounding problem is converted to a subspace sampling problem.

II-B Subspace Scanning

CSI acquisition techniques using subspace sampling were investigated in [6, 7, 8]. Provided a subspace pair $({\mathbf{W}}_{k},{\mathbf{F}}_{k})$, where ${\mathbf{W}}_{k}\!\in\!{\mathbb{C}}^{N_{\mathrm{r}}\times N_{\mathrm{d}}}$ and ${\mathbf{F}}_{k}\!\in\!{\mathbb{C}}^{N_{\mathrm{t}}\times N_{\mathrm{d}}}$, the subspace sample ${\mathbf{y}}_{k}\!\!\in\!{\mathbb{C}}^{N_{\mathrm{d}}\times 1}$​ at the $k$th channel use is

[TABLE]

where ${\mathbf{s}}_{k}\in{\mathbb{C}}^{N_{\mathrm{d}}\times 1}$ is typically an all one vector with proper normalization. Each subspace pair $({\mathbf{W}}_{k},{\mathbf{F}}_{k})$ is chosen from a pre-designed subspace sampling codebook [6, 7, 8]. In particular, if both the digital and analog parts are utilized, we have ${\mathbf{F}}_{k}\!=\!{\mathbf{F}}_{\mathrm{A},k}{\mathbf{F}}_{\mathrm{D},k}$ and ${\mathbf{W}}_{k}\!=\!{\mathbf{W}}_{\mathrm{A},k}{\mathbf{W}}_{\mathrm{D},k}$, $k=1,\ldots,K$.

After collecting $K$ subspace samples, the best pair $({\mathbf{W}}_{k_{opt}},{\mathbf{F}}_{k_{opt}})$ is selected by maximizing the received power,

[TABLE]

The best pair $({\mathbf{W}}_{k_{opt}},{\mathbf{F}}_{k_{opt}})$ can be directly used or further processed to generate the precoder and combiner for data transmission. This framework is investigated in [6, 8] for beamforming and [7] for precoding.

II-C Motivations and General Statement of Technique

Obviously, the decesion in (3) is highly susceptible to noise. When the sampling SNR is low, the noise realizations can lead to incorrect subspace decisions. This is because the subspace scanning approach only looks at the largest power of $\{{\mathbf{y}}_{k}\}_{k=1}^{K}$ rather than correlating the observation with the aligned subspace pairs. Since the millimeter wave channel is sparse [9, 10], advanced approaches must also leverage this fact to enhance the channel estimation. Our proposed approach will therefore take sparse subspace information into account by taking the subspace correlation with the noisy observations.

To this end, we slightly modify the subspace sampling framework in (2) to make full use of $N_{\mathrm{RF}}$ subspace dimension at baseband. At the $k$th channel use, we have

[TABLE]

where ${\mathbf{y}}_{k}\!\in\!{\mathbb{C}}^{N_{\mathrm{RF}}\times 1}$, ${\mathbf{W}}_{k}\!\in\!{\mathbb{C}}^{N_{\mathrm{r}}\times N_{\mathrm{RF}}}$, and ${\mathbf{f}}_{k}\!\in\!{\mathbb{C}}^{N_{\mathrm{t}}\times 1}$. Unlike (2), the receiver in (4) collects $N_{\mathrm{RF}}$ subspace samples per channel use using the subspace pair $({\mathbf{W}}_{k},{\mathbf{f}}_{k})$ where ${\mathbf{W}}_{k}\!=\!{\mathbf{W}}_{\mathrm{A},k}{\mathbf{W}}_{\mathrm{D},k}$ and ${\mathbf{f}}_{k}\!=\!{\mathbf{F}}_{\mathrm{A},k}{\mathbf{F}}_{\mathrm{D},k}{\mathbf{s}}_{k}$ with ${\mathbf{W}}_{\mathrm{D},k}\!\in\!{\mathbb{C}}^{N_{\mathrm{RF}}\!\times\!N_{\mathrm{RF}}}$ and ${\mathbf{F}}_{\mathrm{D},k}\!\in\!{\mathbb{C}}^{N_{\mathrm{RF}}\!\times\!N_{\mathrm{RF}}}$. After collecting $K$ subspace samples, we form ${\mathbf{y}}=[{\mathbf{y}}_{1}^{T},{\mathbf{y}}_{2}^{T},\ldots,{\mathbf{y}}_{K}^{T}]^{T}\!\in{\mathbb{C}}^{(K\!N_{\mathrm{RF}})\times 1}$.

Channel Use Overhead $K$

The subspace sampling in (4) must be done within a channel coherence block $T$. Typically, we assume the number of channel uses $K$ such that $K\!\ll\!T$. Let the minimum required channel uses be $K_{\min}$. Assuming the noiseless scenario, $K_{\min}$ that guarantees the perfect reconstruction of ${\mathbf{H}}\!\in\!{\mathbb{C}}^{N_{\mathrm{r}}\times N_{\mathrm{t}}}$ with $\mathop{\mathrm{rank}}({\mathbf{H}})\!\leq\!d$ is given by $K_{\min}\!=\!d(N_{\mathrm{r}}+N_{\mathrm{t}}-d)$ [15]. Since the sampling in (4) can collect in total $K\!N_{\mathrm{RF}}$ samples, we have the minimum channel uses $K_{\min}\!=\!\frac{d(N_{\mathrm{r}}+N_{\mathrm{t}}-d)}{N_{\mathrm{RF}}}$. Now, considering the noise, the channel use overhead should satisfy

[TABLE]

III Sparse Reconstruction Objective Functions

In this section, we formulate the sparse subspace estimation problem as an optimization problem under the sparse channel constraints. This problem will be solved in the next section.

First off, suppose the precoding and combining gain maximization problem for ${\mathbf{H}}$ with $\mathop{\mathrm{rank}}({\mathbf{H}})\leq d$, i.e.,

[TABLE]

where ${\mathbf{W}}\in{\mathbb{C}}^{N_{\mathrm{r}}\times N_{\mathrm{d}}}$ and ${\mathbf{F}}\in{\mathbb{C}}^{N_{\mathrm{t}}\times N_{\mathrm{d}}}$. The (6) can be written in an alternative form given by

[TABLE]

where ${\mathbf{U}}\!\in\!{\mathbb{C}}^{N_{\mathrm{r}}\times d}$ and ${\mathbf{V}}\!\in\!{\mathbb{C}}^{N_{\mathrm{t}}\times d}$. The $\mathbf{\Lambda}\!\in\!{\mathbb{C}}^{d\times d}$ is a diagonal matrix with the diagonal elements being sorted in descending order. The relation between (6) and (7) is that $[\widehat{{\mathbf{W}}}]_{i}\!=\![\widehat{{\mathbf{U}}}]_{i}$ and $[\widehat{{\mathbf{F}}}]_{i}\!=\![\widehat{{\mathbf{V}}}]_{i}$, $i\!=\!1,2,\cdots,N_{\mathrm{d}}$. The formulation in (7) is important because it shows that maximizing the precoding gain in (6) is equivalent to approximating the low rank matrix ${\mathbf{H}}$ with ${\mathbf{U}}\mathbf{\Lambda}{\mathbf{V}}^{*}$. One challenge is how to efficiently use the $K$ subspace samples in (4) to make a correct approximation.

Inspired by (7), we now formulate our objective function under the noisy observation that minimizes the subspace distance from the received $K$ subspace samples in (4). Mathematically,

[TABLE]

The channel estimate will be $\widehat{{\mathbf{H}}}=\widehat{{\mathbf{U}}}\widehat{\mathbf{\Lambda}}\widehat{{\mathbf{V}}}^{*}$. However, dealing with the summation in (8) is not convenient. To make the problem intuitive, we introduce an affine map ${\mathcal{A}}:{\mathcal{F}}\times{\mathcal{W}}\times{\mathbb{C}}^{N_{\mathrm{r}}\times N_{\mathrm{t}}}\mapsto{\mathbb{C}}^{K\!N_{\mathrm{RF}}\times 1}$, where ${\mathcal{F}}=\{{\mathbf{f}}_{k}\}_{k=1}^{K}$ and ${\mathcal{W}}=\{{\mathbf{W}}_{k}\}_{k=1}^{K}$. Then, ${\mathcal{A}}({\mathcal{F}},\!{\mathcal{W}},\!{\mathbf{H}})$ yields

[TABLE]

The problem in (8) is now rewritten by

[TABLE]

The problem in (10) has much overlap with the rank minimzation problems [13, 14, 15, 16]. The original rank minimization problem [13] is, however, non-convex and is infeasible to be directly solved (i.e., NP-hard). One alternative approach is to use the nuclear norm heuristic [14]. The nuclear norm minimization (NNM) problem leveraging the millimeter wave channel subspace samples in (4) can be written as

[TABLE]

The set ${\mathcal{C}}$ is a convex set which can be adjusted based on the noise power [15]. This problem is convex and solvable.

Though NNM provides robust performance, it does not guarantee the low rank solution, leading to inaccurate channel estimate under the noisy observation. The most significant drawback is that the computational complexity becomes prohibitively high as $N_{\mathrm{r}}$ and $N_{\mathrm{t}}$ grow.

Much like our setting, the matrix factorization (MF) technique [16] is proposed to guarantee the low rank solution and can be formulated as

[TABLE]

where ${\mathbf{U}}\!\in\!{\mathbb{C}}^{N_{\mathrm{r}}\times d}$ and ${\mathbf{V}}\!\in\!{\mathbb{C}}^{N_{\mathrm{t}}\times d}$. An alternating optimization can be used to find a local minimum of (11). The computation complexity is lower than NNM. However, since ${\mathbf{U}}$ and ${\mathbf{V}}$ are not constrained to lie on the unitary subspace, when it is compared with (10) and NNM, it exhibits the worst performance (as shown in Section V).

IV Sparse Millimeter Wave Channel Estimation

Our method will address two major challenges: low-rank guarantee and reduced complexity with robust performance.

IV-A Sparse Subspace Decomposition (SSD)

The objective function in (10) finds the column and row subspaces of ${\mathbf{H}}$, i.e., ${\mathbf{U}}$ and ${\mathbf{V}}$, respectively, along with the power allocation matrix $\mathbf{\Lambda}$ to meet ${\mathbf{U}}^{*}{\mathbf{U}}={\mathbf{I}}_{d}$ and ${\mathbf{V}}^{*}{\mathbf{V}}={\mathbf{I}}_{d}$. Unfortunately, these semi-unitary constraints are not tractable to handle since they are non-convex. Though Cayley Transform [17] can be employed to deal with the orthogonality constraints, the computational complexity is prohibitively high. Instead, we will consider the convex relaxation of (10) as

[TABLE]

Notice that the problem (12) is still non-convex since those optimization parameters are coupled each other. Nevertheless, the problem (12) can be suboptimally but tractably solved by using alternating minimization techniques. In particular, optimizing one parameter by fixing the other two parameters in (12) is convex. We can iteratively optimize ${\mathbf{U}}$, ${\mathbf{V}}$, and $\mathbf{\Lambda}$ by solving the following three subproblems:

• (S1)

Fix the row subspace ${\mathbf{V}}$ and power allocation $\mathbf{\Lambda}$, optimize the column subspace ${\mathbf{U}}$,

• (S2)

Fix the column subspace ${\mathbf{U}}$ and power allocation $\mathbf{\Lambda}$, optimize the row subspace ${\mathbf{V}}$,

• (S3)

Fix the row subspace ${\mathbf{V}}$ and column subspace ${\mathbf{U}}$, optimize the power allocation $\mathbf{\Lambda}$.

A formal description of the alternating minimization is provided in Algorithm 1. Since (S1), (S2), and (S3) are all convex, we can obtain the global optima of each subproblem. Moreover, due to the covexity, the objective function $\|{\mathbf{y}}-{\mathcal{A}}({\mathcal{F}},{\mathcal{W}},{\mathbf{U}}\mathbf{\Lambda}{\mathbf{V}}^{*})\|_{2}^{2}$ converges over the iterations.

Now, the optimal solutions of the subproblems are of interest. We limit our discussion to the column subspace optimization problem (S1), bearing in mind that the same applies to the row subspace optimization (S2). For simplicity, we omit the iteration index $\ell$ attached to the variables throughout this subsection. The following lemma provides the solution to (S1).

Lemma 1

Suppose the following quadratic programming,

[TABLE]

where ${\mathbf{y}}\in{\mathbb{C}}^{K\!N_{\mathrm{RF}}\times 1}$, ${\mathbf{U}}\in{\mathbb{C}}^{N_{\mathrm{r}}\times d}$, ${\mathbf{V}}\in{\mathbb{C}}^{N_{\mathrm{t}}\times d}$, and $\mathbf{\Lambda}\in{\mathbb{C}}^{d\times d}$. Let ${\mathbf{B}}$ be

[TABLE]

where ${\mathbf{B}}\in{\mathbb{C}}^{(K\!N_{\mathrm{RF}})\times(dN_{\mathrm{r}})}$. Then, the optimal solution $\widehat{{\mathbf{U}}}$ is given by either $\mathop{\mathrm{vec}}(\widehat{{\mathbf{U}}})=({\mathbf{B}}^{*}{\mathbf{B}})^{-1}{\mathbf{B}}^{*}{\mathbf{y}}$ such that $\|\mathop{\mathrm{vec}}(\widehat{{\mathbf{U}}})\|_{2}^{2}\leq d$ or $\mathop{\mathrm{vec}}(\widehat{{\mathbf{U}}})\!=\!({\mathbf{B}}^{*}{\mathbf{B}}\!+\!\lambda{\mathbf{I}}_{dN_{\mathrm{r}}})^{-1}{\mathbf{B}}^{*}{\mathbf{y}}$ such that $\ g(\lambda)\!\triangleq\!\|\mathop{\mathrm{vec}}(\widehat{{\mathbf{U}}})\|_{2}^{2}\!=\!d$, where $\lambda>0$ is the unique solution to the fixed point equation $g(\lambda)\!=\!d$, in which $g(\lambda)$ is monotonically decreasing in $\lambda>0$111The $\lambda>0$ that satisfies $\|\mathop{\mathrm{vec}}({\widehat{{\mathbf{U}}}})\|_{2}^{2}=d$ can be tractably found by using a bisection method..

Proof:

Given the affine map ${\mathcal{A}}$ in (9), using the Kronecker product equality, $\mathop{\mathrm{vec}}({\mathbf{W}}_{k}^{*}{\mathbf{U}}\mathbf{\Lambda}{\mathbf{V}}^{*}{\mathbf{f}}_{k})=((\mathbf{\Lambda}{\mathbf{V}}^{*}{\mathbf{f}}_{k})^{T}\otimes{\mathbf{W}}_{k}^{*})\mathop{\mathrm{vec}}({\mathbf{U}})$, for $k=1,\ldots,K$ and collecting them as a column vector yield an equivalent problem to (16) as

[TABLE]

where ${\mathbf{B}}$ is defined in (17). This problem is convex and applying KKT condition [18] for (18) gives

[TABLE]

where $\lambda$ is the Lagrange multiplier for the constraint in (18). By the complementary slackness, the optimal solution is found by considering two cases, i.e., $\lambda=0$ and $\left\|\mathop{\mathrm{vec}}({{\mathbf{U}}})\right\|_{2}^{2}-d=0$.

When $\lambda\!=\!0$, we have $\mathop{\mathrm{vec}}(\widehat{{\mathbf{U}}})\!=\!({\mathbf{B}}^{*}{\mathbf{B}})^{-1}{\mathbf{B}}^{*}{\mathbf{y}}$ with $\|\mathop{\mathrm{vec}}({\widehat{{\mathbf{U}}}})\|_{2}^{2}\!\leq\!d$. According to (5), $K\!N_{\mathrm{RF}}\!\geq\!{d(N_{\mathrm{r}}+N_{\mathrm{t}}-d)}\!\geq\!dN_{\mathrm{r}}$, hence, ${\mathbf{B}}^{*}{\mathbf{B}}$ is invertible.

When $\lambda\!>\!0$, $\left\|\mathop{\mathrm{vec}}({{\mathbf{U}}})\right\|_{2}^{2}\!=\!d$ holds and the optimal $\widehat{{\mathbf{U}}}$ should meet the two conditions $\mathop{\mathrm{vec}}({\mathbf{U}})=({\mathbf{B}}^{*}{\mathbf{B}}+\lambda{\mathbf{I}})^{-1}{\mathbf{B}}^{*}{\mathbf{y}},~{}\lambda>0$ and $\left\|\mathop{\mathrm{vec}}({{\mathbf{U}}})\right\|_{2}^{2}=d$. This concludes the proof.∎∎

We turn our attention to the subspace power allocation problem (S3).

Lemma 2

Suppose the subspace power allocation problem

[TABLE]

where $\mathbf{\Lambda}\!\in\!{\mathbb{C}}^{d\times d}$ is the diagonal matrix, ${\mathbf{y}}\!\in\!{\mathbb{C}}^{K\!N_{\mathrm{RF}}\times 1}$, ${\mathbf{U}}\!\in\!{\mathbb{C}}^{N_{\mathrm{r}}\times d}$, and ${\mathbf{V}}\!\in\!{\mathbb{C}}^{N_{\mathrm{t}}\times d}$. Then the optimal solution to (19) is

[TABLE]

where ${\mathbf{P}}\!\in\!{\mathbb{C}}^{K\!N_{\mathrm{RF}}\times d}$, and $[{\mathbf{P}}]_{j}\!=\!{\mathcal{A}}({\mathcal{F}},{\mathcal{W}},[{\mathbf{U}}]_{j}[{\mathbf{V}}]_{j}^{*}),~{}\forall j$.

Proof:

Decompose the affine map ${\mathcal{A}}$ in (9) as

[TABLE]

where $[\mathbf{\Lambda}]_{j,j}$ in (20) is the power allocation to the corresponding subspace $[{\mathbf{U}}]_{j}[{\mathbf{V}}]_{j}^{*}$. Using (20), the problem (19) can be rewritten as

[TABLE]

Solving the least square problem in (21) results in

[TABLE]

where ${\mathbf{P}}\!\in\!{\mathbb{C}}^{K\!N_{\mathrm{RF}}\times d}$, and $[{\mathbf{P}}]_{j}\!=\!{\mathcal{A}}({\mathcal{F}},{\mathcal{W}},[{\mathbf{U}}]_{j}[{\mathbf{V}}]_{j}^{*}),~{}j=1,\cdots,d$. This concludes the proof. ∎∎

IV-B ​Sparse Subspace​ Decomposition with​ Thresholding ​(SSD-T)

In our channel estimation problem, a normalized mean square error (NMSE), i.e.,

[TABLE]

is evaluated to confirm the accuracy of an estimator. Obviously, there is a gap between the objective function in (12) and NMSE in (22). Though our objective function monotonically decreases per iteration to a local minimum point, the NMSE in (22) would exhibit a non-monotonic behavior. This is because the objective in (12) contains the disturbance, i.e., the noise ${\mathbf{n}}$. The decoupling between the objective and NMSE becomes critical in a low SNR, especially when $\|{\mathbf{y}}\|_{2}^{2}\approx\|\tilde{{\mathbf{n}}}\|_{2}^{2}$, where $\tilde{{\mathbf{n}}}=[({\mathbf{W}}_{1}^{*}{\mathbf{n}}_{1})^{T},\cdots,({\mathbf{W}}_{K}^{*}{\mathbf{n}}_{K})^{T}]^{T}$. In this case, the objective function is approximated as $\|\tilde{{\mathbf{n}}}\!-\!{\mathcal{A}}({\mathcal{F}},{\mathcal{W}},\widehat{{\mathbf{H}}})\|_{2}^{2}$. Minimizing $\|\tilde{{\mathbf{n}}}\!-\!{\mathcal{A}}({\mathcal{F}},{\mathcal{W}},\widehat{{\mathbf{H}}})\|_{2}^{2}$ results in large deviation from the true ${\mathbf{H}}$.

To illustrate this in a low SNR scenario, we display the curves of the objective function of (12) and NMSE of the proposed SSD (i.e., Algorithm 1) across iteration $\ell$ in Fig. 1 with [math] dB SNR. As expected, the objective function of (12) decreases monotonically. However, NMSE rather increases after a few iterations because the algorithm finds $\widehat{{\mathbf{H}}}$ close to the noise, resulting in the over-matching.

One remedy is to terminate the iteration rather than making $\|{\mathbf{y}}-{\mathcal{A}}({\mathcal{F}},{\mathcal{W}},\widehat{{\mathbf{H}}})\|_{2}^{2}$ as small as possible. Specifically, provided the output of Algorithm 1 at iteration $\ell$, i.e., $\widehat{{\mathbf{H}}}_{(\ell)}$, the stopping condition can be given by

[TABLE]

The underlying intuition of (23) is that by the strong law of large number the left hand side (l.h.s.) of (23) will converge to the variance of noise, i.e., $\sigma^{2}$, as $K$ grows, given that the output of Algorithm 1 is close to the true ${\mathbf{H}}$. The inequality in (23) signifies that we treat the case when the l.h.s. of (23) is smaller than $\sigma^{2}$ as the over-matching and stop the iteration. The formal description of the SSD combined with the thresholding in (23), namely SSD with thresholding (SSD-T), is provided in Algorithm 2. By reducing the number of iterations, Algorithm 2 will decrease the computational complexity while enhancing the performance in a low SNR.

V Numerical Simulations

In this section, we numerically evaluate the NMSE of the proposed SSD and SSD-T algorithms and compare them with the matrix factorization (MF) [16] and nuclear norm minimization (NNM) [15] techniques. Also, we benchmark the computational overheads of SSD, SSD-T, MF, and NNM. The numerical simulation setting is first discussed.

(1) Channel model: We assume the prevalent physical channel representation that models sparse millimeter wave MIMO channels [1, 10]. Assume there are $L$ propagation paths between the transmitter and receiver, where the channel is modeled via

[TABLE]

The $\alpha_{l}$ is the complex gain of the $l$th path, i.i.d. as $\alpha_{l}\sim{\mathcal{C}}{\mathcal{N}}(0,1)$, $\phi_{\mathrm{t},l}$ and $\phi_{\mathrm{r},l}$ are angles of departure and arrival at the transmitter and receiver, respectively. The ${\mathbf{a}}_{\mathrm{t}}(\phi_{\mathrm{t},l})$ and ${\mathbf{a}}_{\mathrm{r}}(\phi_{\mathrm{r},l})$ are the uniform linear array (ULA) response vectors at the transmitter and receiver, respectively, where the inter-element spacing of ULA is set to half of the wavelength. The maximum number of iteration of the proposed SSD (i.e., Algorithm 1) is 30. We further assume $L\!=\!2$ paths, $N_{\mathrm{RF}}\!=\!4$ RF chains, and the number of data stream $N_{\mathrm{d}}\!=\!L$. We assume the rank upper bound $d\!=\!3$.

(2) Performance Evaluation: The NMSE statistics across different SNRs and channel uses ($K$) are evaluated. Each curve is obtained after averaging over $100$ channel realizations.

Table I shows the execution time statistics (in seconds on a $3.4$ GHz CPU) for $8\times 16$ and $16\times 64$ channel dimensions when SNR = 10 dB. It can be seen from Table I that both proposed algorithms exhibit lower complexity. The SSD-T shows the fastest execution while the NNM becomes substantially slower, as the antenna dimension grows.

Fig. 2 illustrates the NMSE curves of the proposed SSD and SSD-T compared with MF and NNM for $16\times 64$ channel and $K\!=\!100$ channel uses. It can be seen that both SSD and SSD-T outperform the MF across a broad range of SNRs. NNM achieves robust performance at the expense of the huge computational overhead. On the other hand, Fig. 2 reveals that SSD-T achieves similar performance as NNM in the low SNR regime, while it starts to outperform NNM in the high SNR, verifying that SSD-T has the ability to resolve the complexity and accuracy tradeoffs.

Fig. 3 presents NMSE curves at $10$ dB and $30$ dB SNRs across different channel uses. Fig. 3 clearly shows that both SSD and SSD-T achieve similar performances with NNM at $10$ dB SNR, while keeping the complexity remarkably low. At $30$ dB SNR, both SSD and SSD-T outperform NNM. The gap of the proposed approaches, compared with MF, is clear in Fig. 3. Evidently, increasing $K$ in the high SNR has the effect of enhancing the accuracy of SSD and SSD-T, while only marginal improvement is observed for NNM.

VI Conclusions

In this paper, we investigated the estimation of sparse millimeter wave MIMO channels using hybrid precoding systems. Based on the equivalency of optimizing the unitary precoder/combiner of the sparse channel matrix and approximating the channel via a sparse subspace decomposition, we proposed the SSD algorithm. The SSD consists of three constituent estimators: one for the row subspace, one for the column subspace, and the other for the subspace power allocation. To avoid unexpected over-matching issues, we devised the SSD-T algorithm. Our simulations showed that SDD and SDD-T achieve similar or better accuracy than the best performing existing sparse matrix estimation benchmarks (e.g., NNM algorithm) with remarkably low overhead.