Beam Allocation for Millimeter-Wave MIMO Tracking Systems

TL;DR

This paper introduces a novel beam allocation strategy for millimeter-wave MIMO systems that enhances tracking performance by allowing repeated use of beam pairs and optimizing their selection through advanced algorithms.

Contribution

It proposes a new framework for beam allocation in mmWave MIMO systems, including algorithms for both orthogonal and overlapping beam pairs, with a focus on maximizing successful tracking probability.

Findings

The proposed strategy outperforms existing benchmarks in simulations.

A low-complexity solution effectively approximates the optimal beam allocation.

Derived a lower bound for the ASTP in OMP-based estimators.

Abstract

In this paper, we propose a new beam allocation strategy aiming to maximize the average successful tracking probability (ASTP) of time-varying millimeter-wave MIMO systems. In contrast to most existing works that employ one transmitting-receiving (Tx-Rx) beam pair once only in each training period, we investigate a more general framework, where the Tx-Rx beam pairs are allowed to be used repeatedly to improve the received signal powers in specific directions. In the case of orthogonal Tx-Rx beam pairs, a power-based estimator is employed to track the time-varying AoA and AoD of the channel, and the resulting training beam pair sequence design problem is formulated as an integer nonlinear programming (I-NLP) problem. By dividing the feasible region into a set of subregions, the formulated I-NLP is decomposed into a series of concave sub I-NLPs, which can be solved by recursively invoking a nonlinear branch-and-bound algorithm. To reduce the computational cost, we relax the integer constraints of each sub I-NLP and obtain a low-complexity solution via solving the Karush-Kuhn-Tucker conditions of their relaxed problems. For the case when the Tx-Rx beam pairs are overlapped in the angular space, we estimate the updated AoA and AoD via an orthogonal matching pursuit (OMP) algorithm. Moreover, since no explicit expression for the ASTP exists for the OMP-based estimator, we derive a closed-form lower bound of the ASTP, based on which a favorable beam pair allocation strategy can be obtained. Numerical results demonstrate the superiority of the proposed beam allocation strategy over existing benchmarks.

Equations

$$\mathbf{r}=\mathbf{H}\mathbf{f}x+\tilde{\mathbf{n}},$$

$$y={\mathbf{w}}^{\mathrm{H}}\mathbf{r}={\mathbf{w}}^{\mathrm{H}}\mathbf{H}\mathbf{f}x+n,$$

$$\mathbf{H}=\sqrt{N_{\mathrm{T}}{N}_{\mathrm{R}}}\sum _{\ell =1}^L{\alpha }_{\ell }{\mathbf{a}}_{\mathrm{R}}({\theta }_{\ell }){\mathbf{a}}_{\mathrm{T}}^{\mathrm{H}}({\vartheta }_{\ell }),$$

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

$${\mathbf{a}}_{\mathrm{T}}({\vartheta }_{\ell })$$

$${\theta }_{\ell }=\frac{2\pi d\sin ({\tilde{\theta }}_{\ell })}{{\lambda }_s},{\vartheta }_{\ell }=\frac{2\pi d\sin ({\tilde{\vartheta }}_{\ell })}{{\lambda }_s},\forall \ell =1,\cdots ,L,$$

$${\bar{\theta }}_m$$

$${\bar{\vartheta }}_n$$

$$\mathrm{Pr}\left\{{\theta }_{\ell }^{[\tau ]}={\bar{\theta }}_{k_1}\leftarrow {\theta }_{\ell }^{[\tau -1]}={\bar{\theta }}_{k_0}\right\}=C{\beta }^{|k_1-k_0|},$$

$$\mathrm{Pr}\left\{{\vartheta }_{\ell }^{[\tau ]}={\bar{\vartheta }}_{i_1}\leftarrow {\vartheta }_{\ell }^{[\tau -1]}={\bar{\vartheta }}_{i_0}\right\}=\tilde{C}{\tilde{\beta }}^{|i_1-i_0|},$$

$$y_{k,i}$$

$$\displaystyle\xi_{k,i}=\sum\limits_{m=1}^{\lambda_{k,i}}{\bf y}_{k,i}[m]\sim\begin{dcases}{\cal CN}\big{(}\gamma\alpha\lambda_{k,i},\sigma_{0}^{2}\lambda_{k,i}\big{)},&\text{if}~{}k=k_{1}~{}\text{and}~{}i=i_{1},\\ {\cal CN}\big{(}0,\sigma_{0}^{2}\lambda_{k,i}\big{)},&\text{otherwise}.\\ \end{dcases}$$

$$\displaystyle\Gamma_{z_{n},|\alpha|^{2}}=\Pr\Bigg{(}\bigcap\limits_{m=1,m\neq n}^{N}|\xi_{z_{n}}|^{2}>|\xi_{z_{m}}|^{2}~{}\Big{|}~{}\alpha\Bigg{)}.$$

$$\displaystyle\Gamma_{z_{n},|\alpha|^{2}}=\int\limits_{0}^{\infty}h\bigg{(}u;\lambda_{z_{n}},|\alpha|^{2}\bigg{)}\prod\limits_{m=1,m\neq n}^{N}\left(1-\exp\left(-\frac{u}{\lambda_{z_{m}}\sigma_{0}^{2}}\right)\right)\text{d}u,$$

$$\displaystyle h\Big{(}u;\lambda_{z_{n}},|\alpha|^{2}\Big{)}=\frac{1}{\lambda_{z_{n}}\sigma_{0}^{2}}\exp\left(-\frac{u+\lambda^{2}_{z_{n}}\gamma^{2}|\alpha|^{2}}{\lambda_{z_{n}}\sigma_{0}^{2}}\right)I_{0}\left(\frac{\sqrt{4\gamma^{2}|\alpha|^{2}u}}{\sigma_{0}^{2}}\right),$$

$${\Gamma }_{z_n}$$

$${\bar{\Gamma }}_1=\sum _{n=1}^N{\pi }_{z_n}\times {\Gamma }_{z_n},$$

$$\underset{\mathbf{\lambda }}{\max }$$

$$s.t.$$

$$\underset{\mathbf{\lambda }}{\max }$$

$$s.t.$$

$${\bar{\Gamma }}_1^{\mathrm{lb}}(\mathbf{\lambda })=\sum _{n=1}^N{\pi }_{s_n}\left[1-\frac{M_{\mathrm{B}}-{\lambda }_{s_n}}{{\lambda }_{s_n}^2{r}_0+{\lambda }_{s_n}}\right],$$

$${\bar{\Gamma }}_1^{\mathrm{ub}}(\mathbf{\lambda })=\sum _{n=1}^N{\pi }_{s_n}\left[1-\frac{\mathcal{F}(N-1)}{{\lambda }_{s_n}^2{r}_0+{\lambda }_{s_n}}\right].$$

$${\bar{\Gamma }}_1^{\mathrm{apx}}(\mathbf{\lambda })=\sum _{n=1}^N{\pi }_{s_n}\left[1-\frac{\mathcal{F}(N-1)}{N-1}\frac{M_{\mathrm{B}}-{\lambda }_{s_n}}{{\lambda }_{s_n}^2{r}_0+{\lambda }_{s_n}}\right],$$

$${\bar{\Gamma }}_1^{\mathrm{apx}}(\mathbf{\lambda })=\sum _{n=1}^N{\pi }_{s_n}\left[1-\frac{M_{\mathrm{B}}-{\lambda }_{s_n}}{N-1}\frac{\mathcal{F}(N-1)}{{\lambda }_{s_n}^2{r}_0+{\lambda }_{s_n}}\right],$$

$${\bar{\Gamma }}_1^{\mathrm{lb}}(\mathbf{\lambda })<{\bar{\Gamma }}_1^{\mathrm{apx}}(\mathbf{\lambda })<{\bar{\Gamma }}_1^{\mathrm{ub}}(\mathbf{\lambda }).$$

$$\underset{\mathbf{\lambda }}{\max }$$

$$s.t.$$

$${\bar{\Gamma }}_1^{\mathrm{apx}}({\mathbf{\lambda }}_N)=\sum _{n=1}^N{\pi }_{s_n}\left[1-\frac{\mathcal{F}(N-1)}{N-1}\frac{M_{\mathrm{B}}-{\lambda }_{s_n}}{{\lambda }_{s_n}^2{r}_0+{\lambda }_{s_n}}\right]<\sum _{n=1}^N{\pi }_{s_n}\le \Delta ,$$

$${\bar{\Gamma }}_1^{\mathrm{apx}}({\mathbf{\lambda }}_K)=\sum _{n=1}^K{\pi }_{s_n}\left[1-\frac{\mathcal{F}(K-1)}{K-1}\frac{M_{\mathrm{B}}-{\lambda }_{s_n}}{{\lambda }_{s_n}^2{r}_0+{\lambda }_{s_n}}\right]<\sum _{n=1}^K{\pi }_{s_n}<\Delta ,$$

$$\underset{\mathbf{\lambda }}{\min }$$

$$s.t.$$

$$\mathcal{L}=\sum _{n=1}^N\frac{{\pi }_{s_n}(M_{\mathrm{B}}-{\tilde{\lambda }}_{s_n})}{{\tilde{\lambda }}_{s_n}^2{r}_0+{\tilde{\lambda }}_{s_n}}+{\mu }_0\left(\sum _{n=1}^N{\tilde{\lambda }}_{s_n}-M_{\mathrm{B}}\right)-\sum _{n=1}^N{\mu }_n\left({\tilde{\lambda }}_{s_n}-1\right),$$

Full text

I Introduction

In recent years, with the rapid penetration of mobile broadband Internet and multimedia applications, the ever-increasing data rate demand has made current sub-6GHz bands unprecedentedly crowded. Thanks to the rich available spectrum, millimeter-wave (mmWave) communication ranging from 30GHz to 300GHz emerged as a promising solution. However, due to the high frequency, mmWave signals suffer from severe propagation loss and atmospheric absorption [2, 3]. Fortunately, the short operating wavelength enables a large antenna array to be integrated in a compact form, providing a considerable beamforming gain to compensate for the path loss in mmWave radio links [4]. Nevertheless, such beamforming requires accurate channel state information (CSI) at the transmitter and receiver, which is usually difficult to obtain in mmWave communications. To be specific, the high cost and power consumption of mixed-signal hardware with high sampling rates limit the number of radio-frequency (RF) chains in practical mmWave transceivers [5], making fully-digital processing techniques such as the traditional least-squared method impractical [6]. On the other hand, mmWave channels usually have a limited number of propagation paths, and therefore it is sufficient to only estimate channel parameters of these paths, which include angle of arrivals (AoAs), angle of departures (AoDs), and propagation gains.

There have been extensive work on channel estimation or beam training for mmWave systems over the years [7, 8, 10, 9, 11, 12, 13]. A widely used channel estimation method is to sequentially transmit highly directional training beams steering to different directions over time and pick the direction with the largest received signal-to-noise ratio (SNR) [7]. Nevertheless, this method is time-consuming and the required number of training beams is usually large in order to achieve a favorable estimate of the channel. This problem becomes even more challenging in mobile scenarios, where the channel keeps changing and the transmitter needs to frequently send training beams to update the estimation results, increasing the training overhead considerably. Therefore, an efficient and accurate beam training strategy is crucially important in mobile scenarios. In [8, 10, 9], adaptive compressed sensing (CS) algorithms have been used to estimate the mmWave channels, which essentially search potential paths using hierarchical multi-resolution codebooks. These adaptive CS algorithms can achieve a favorable estimation performance but require excessive feedback, which may significantly exacerbate the system overhead. Meanwhile, all the aforementioned channel estimation algorithms fail to capture the temporal correlation between consecutive channel realizations in mobile mmWave communications.

Recently, it has been shown in [14, 15, 16, 17, 18, 19] that the temporal correlation between channel realizations can be exploited to further improve the beam training efficiency. Specifically, the CSI in current channel realization is closely related to that of the previous one, and this relationship can be used to speed up the beam training procedure. This type of priori-aided beam training technique is referred to as channel tracking or beam tracking in the literature [14, 15, 16, 17]. To date, most existing works that investigate beam tracking techniques for mmWave channels have focused on the assumption that the values of AoA and AoD vary smoothly. In [15, 16, 17, 18, 19], the temporal variation of AoA/AoD over the considered period of time is assumed to follow a Markov process, and the AoA’s and AoD’s deviations between two consecutive channel realizations are modeled as small Gaussian random variables, based on which various Kalman filter-based beam tracking algorithms have been developed. It is also worth mentioning that the authors in [20, 21, 22, 23] have proposed to employ the mobile users’ location and trajectory information to reduce the beam training overhead. However, these strategies are limited to vehicular networks and not universal.

To incorporate the abrupt changes of mmWave channels due to blockage or other environmental obstructions, several works have proposed to employ the discrete Markov process to model the temporal variations of AoA/AoD [24, 25, 26, 27, 28]. It has been shown in [24] that the problem of tracking the time-varying AoA and AoD can be transformed into finding the support of the sparse beamspace channel, which is solved by invoking an approximate message passing algorithm. In [25, 26, 27, 28], codebook-based training beamforming vectors (beams) are adopted to reduce the design complexity. Specifically, a set of codewords (a beam codebook consists of a sequence of codewords and each codeword is a beamforming vector steering to a specific direction) that can minimize the Cramer-Rao lower bound averaged over the priori distribution of AoD are selected for beam tracking in [26] and [27], and the maximum likelihood (ML) and maximum a posteriori (MAP) criteria are respectively used to estimate the true direction of AoD. In [28], the beam tracking problem is equivalent to a partially observable Markov decision process (POMDP), where the selected training beams serving as actions of the POMDP depend on the belief vector, observation information and reward. In the high SNR regime, the observations of the beam (angular) space based on signal detection are reliable, and hence the selections of training beams using this POMDP framework are appropriate. However, the consequent high power consumption will significantly increase the dynamic range of the power amplifiers, considerably increasing the hardware cost. On the other hand, when SNR goes low, the observations of the beam space will be inaccurate, and therefore selecting the optimal training beam sequence based on this POMDP framework are unreliable and the true AoAs/AoDs can be lost.

To address the aforementioned problems, we develop a new beam pair allocation strategy aiming to maximize the average successful tracking probability (ASTP) of time-varying mmWave multiple-input multiple-output (MIMO) systems, which can work effectively for all SNR regimes. Motivated by [26, 27, 28], the temporal variations of AoA and AoD within the considered period of time are modeled as two discrete Markov processes, described by their associated transition probabilities respectively, which are assumed to be known. Highly directional transmitting (Tx) and receiving (Rx) training beams picked from two predefined codebook matrices are used to combat the severe propagation loss. To further increase the received signal powers in specific directions and consequently improve the beam tracking performance, we allow the Tx-Rx beam pairs steering to these directions to be used repeatedly in the beam training period, which is different from most existing works111It is worth mentioning that the Tx-Rx training beam pairs are also allowed to be used repeatedly in [29]. However, since an online stochastic optimization known as multi-armed bandit algorithm is considered, one feedback of the received energy is needed after each measurement, which may significantly exacerbate the system overhead.. In the following, we summarize the methodologies and main contributions of this paper:

1. In the case of orthogonal Tx-Rx beam pairs, a power-based estimator that returns the direction with the largest received signal power is employed to track the time-varying AoA and AoD, leading to a closed-form expression for the (one-step) ASTP. Since the number of repetitions of each Tx-Rx beam pair can only be an integer, selecting the optimal set of the Tx-Rx beam pairs and determining their associated repetition times to maximize the ASTP is shown to be equivalent to an integer nonlinear programming (I-NLP). Though determining the optimal Tx-Rx beam pair sequence is NP-hard, we prove that the Tx-Rx beam pair with a higher transition probability should be used more times than those with lower transition probabilities in order to achieve the maximal ASTP.

2. It is very challenging to optimize the exact repetition times of each Tx-Rx beam pair due to the complicated expression for the ASTP. Therefore, we derive a tractable approximation for the ASTP as the new objective function of the formulated I-NLP. Afterward, we divide its feasible region into a set of subregions and construct a series of concave sub I-NLPs, which can be solved via recursively invoking a nonlinear branch-and-bound (N-BB) algorithm. To avoid the computational cost of the iterative N-BB algorithm, we relax the integer constraints of each sub I-NLP and obtain a promising solution by solving the Karush-Kuhn-Tucker (KKT) conditions of these relaxed subproblems following a similar recursive manner as in the iterative N-BB algorithm.

3. For the more general scenario where the Tx-Rx training beam pairs are overlapped in the angular space, the power-based estimator performs poorly due to the non-negligible inter-beam interference. In this case, we modify the power-based estimator and propose to track the time-varying AoA and AoD via an orthogonal matching pursuit (OMP) algorithm, which essentially exploits the inter-beam interference to improve the ASTP. Moreover, since no explicit expression for the ASTP exists when the OMP-based estimator is adopted, a closed-form lower bound of the ASTP can be derived, based on which a favorable beam pair allocation strategy to maximize the ASTP is obtained. Our numerical results demonstrate the superiority of the proposed beam pair allocation strategy over the uniform and proportional allocation strategies.

The rest of the paper is organized as follows. In Section II, we describe the considered mmWave system model and the adopted beam training protocol. In Section III, orthogonal Tx-Rx training beam pairs are assumed and a power-based estimator is used to track the time-varying AoA and AoD. In Section IV, we consider a more general scenario where the Tx-Rx training beam pairs are overlapped in the angular space. Numerical results are provided in Section V, followed by the conclusions in Section VI.

Notations: Bold uppercase $\bf A$ and lowercase $\bf a$ represent matrices and column vectors respectively, and non-bold letters are scalars. ${\bf A}^{\ast}$, ${\bf A}^{\rm T}$, and ${\bf A}^{\rm H}$ represent the conjugate, transpose, conjugate transpose of $\bf A$, respectively. ${\bf A}[m,:]$, ${\bf A}[:,n]$ and ${\bf A}[m,n]$ are the $m$-th row, the $n$-th column, and the $(m,n)$-th element of $\bf A$, respectively. ${\bf A}\otimes{\bf B}$ is the Kronecker product of $\bf A$ and $\bf B$. $(b_{1}\bullet b_{2})_{N}\triangleq b_{1}+N(b_{2}-1)$. ${\cal{CN}}({\bf a},{\bf A})$ denotes a complex Gaussian distribution with mean $\bf a$ and covariance matrix $\bf A$. ${\bf I}$ is the identity matrix. ${\bf a}\triangleq\text{vec}({\bf A})$ is the vectorization operation by stacking the columns of $\bf A$ into a vector $\bf a$. \big{\{}a_{1},\cdots,a_{n}\big{\}}\big{\backslash}\big{\{}b_{1},\cdots,b_{k}\big{\}} represents the set \big{\{}a_{1},\cdots,a_{n}\big{\}} excluding \big{\{}b_{1},\cdots,b_{k}\big{\}}. $\displaystyle\binom{n}{k}$ is the number of $k$-combinations of an $n$-element set. $\exp$ is the exponential function. \Pr\big{\{}S_{2}=b_{2}\leftarrow S_{1}=b_{1}\big{\}} represents the transition probability from $S_{1}=b_{1}$ to $S_{2}=b_{2}$, while \Pr\big{(}S_{2}=b_{2}\mid S_{1}=b_{1}\big{)} is the probability of $S_{2}=b_{2}$ conditioned on $S_{1}=b_{1}$. ${\mathds{N}}^{+}$ and ${\mathds{N}}^{++}$ denote the nonnegative integer set and positive integer set, respectively.

II System Model

We consider a mmWave MIMO system, in which a base station (BS) equipped with $N_{\rm T}$ antennas communicates with a mobile station (MS) equipped with $N_{\rm R}$ antennas. Denote the $N_{\rm T}\times 1$ transmitting beamforming vector and the $N_{\rm R}\times 1$ receiving beamforming vector by $\bf f$ and $\bf w$ respectively, which are normalized to satisfy $\|{\bf f}\|^{2}=\|{\bf w}\|^{2}=1$. Moreover, the pilot symbol is denoted by $x=\sqrt{P}$, where $P$ is the power consumed per transmission in the beam training period, and the received signal is then written as

[TABLE]

where $\bf H$ is the $N_{\rm R}\times N_{\rm T}$ channel matrix between the MS and BS, and $\tilde{\bf n}$ is the $N_{\rm R}\times 1$ complex additive white Gaussian noise, i.e., $\tilde{\bf n}\sim{\cal CN}({\bm{0}},\sigma_{0}^{2}{\bf I})$. The MS adopts the receiving beamforming vector $\bf w$ to process the received signal $\bf r$, given by

[TABLE]

where $n={\bf w}^{\rm H}\tilde{\bf n}$ and $n\sim{\cal CN}(0,\sigma_{0}^{2})$ due to $\|{\bf w}\|^{2}=1$.

II-A Millimeter-Wave Channel

Since mmWave channels have a very limited number of scatters, as in [30], we use the geometric channel model to express $\bf H$ as

[TABLE]

where $L$ is the total number of propagation paths, $\alpha_{\ell},\theta_{\ell},\vartheta_{\ell}$ are the complex gain, the normalized AoA and AoD of the $\ell$-th path respectively, and $\alpha_{\ell}\sim{\cal CN}(0,\sigma_{\alpha}^{2})$, $\forall\ell=1,\cdots,L$. In addition, ${\bf a}_{\rm R}(\theta_{\ell})$ and ${\bf a}_{\rm T}(\vartheta_{\ell})$ are termed antenna array response vectors. In this paper, we assume that uniform linear arrays (ULAs) are used at the MS and BS, and therefore ${\bf a}_{\rm R}(\theta_{\ell})$ and ${\bf a}_{\rm T}(\vartheta_{\ell})$ can be respectively written as

[TABLE]

where $j$ is the imaginary unit, i.e., $j=\sqrt{-1}$. Moreover, the relationship between the normalized AoA (AoD) and the physical AoA (AoD) is expressed as [28]

[TABLE]

where $\lambda_{s}$ is the signal wavelength, $d$ is the distance between two adjacent antenna elements, $\tilde{\theta}_{\ell}$ and $\tilde{\vartheta}_{\ell}$ are the physical AoA and AoD of the $\ell$-th path respectively. By letting $d=\frac{\lambda_{s}}{2}$, we obtain that $\theta_{\ell},\vartheta_{\ell}\in[-\pi,\pi]$ when $\tilde{\theta}_{\ell},\tilde{\vartheta}_{\ell}\in[-\pi,\pi]$.

II-B Time-Varying AoA and AoD

In order to exploit the sparsity of the mmWave channels, as in [8, 11, 27, 28], we introduce two beam codebook matrices {\bf A}_{\rm R}=\big{[}{\bf a}_{\rm R}(\bar{\theta}_{1}),{\bf a}_{\rm R}(\bar{\theta}_{2}),\cdots,{\bf a}_{\rm R}(\bar{\theta}_{X_{\rm R}})\big{]} and {\bf A}_{\rm T}=\big{[}{\bf a}_{\rm T}(\bar{\vartheta}_{1}),{\bf a}_{\rm T}(\bar{\vartheta}_{2}),\cdots,{\bf a}_{\rm T}(\bar{\vartheta}_{X_{\rm T}})\big{]}, where

[TABLE]

to divide the whole angular space $[-\pi,\pi]$ into $X_{\rm R}$ and $X_{\rm T}$ directions, respectively. In addition, we follow [5, 8, 28] and assume that $\{\theta_{\ell}\}$ and $\{\vartheta_{\ell}\}$ are respectively taken from the sets $\{\bar{\theta}_{m}\}_{m=1}^{X_{\rm R}}$ and $\{\bar{\vartheta}_{n}\}_{n=1}^{X_{\rm T}}$ for simplicity222Though this on-grid assumption of the normalized AoAs/AoDs may not be rigorous, the resulting quantization errors are not significant when $X_{\rm R}$ and $X_{\rm T}$ are large enough [5, 8, 28]. The off-grid case where the values of $\{\theta_{\ell}\}$ and $\{\vartheta_{\ell}\}$ are continuous is left as one of our future works..

In mobile scenarios, as mentioned earlier, the channel realizations between two consecutive transmission blocks are correlated. Following [26, 27, 28], we model the temporal variations of each AoA (AoD) within a set of transmission blocks as a discrete Markov process (see Fig. 1 and Fig. 2), described by the following transition probability

[TABLE]

where k_{0},k_{1}\in\big{\{}1,\cdots,X_{\rm R}\big{\}}. The superscript $\tau$ denotes the $\tau$-th beam training period or transmission block, where $\tau=2,\cdots,T$. The variable $\beta\in[0,1]$ indicates the variation speed of the AoAs, and $C$ is the normalization coefficient. According to (4), when $\beta$ is small, e.g., $\beta=0.1$, the updated AoA $\theta_{\ell}^{[\tau]}$ is very likely to be in the proximity to $\theta_{\ell}^{[\tau-1]}$. On the other hand, when $\beta=1.0$, the AoAs change rapidly and $\theta_{\ell}^{[\tau]}$ will be uniformly distributed in the set $\{\bar{\theta}_{m}\}_{m=1}^{X_{\rm R}}$, which corresponds to the abrupt changes in mmWave channels. The associated transition probability of each AoD can be similarly expressed as in (4), given by

[TABLE]

where i_{0},i_{1}\in\big{\{}1,\cdots,X_{\rm T}\big{\}} and $\tilde{\beta}$ is introduced to indicate the variation speed of the AoDs. Moreover, the channel gains $\{\alpha_{\ell}\}$ are assumed to change independently from one transmission block to another [28]. Finally, it is worth highlighting that the proposed beam pair allocation strategy can be extended to other types of transition probabilities such as the ones used in [26] and [27].

II-C Beam Training Protocol

To track the time-varying AoAs and AoDs, which can result from the mobility of the MS or the reflection scatters, the BS transmits a sequence of pilot symbols via using a set of dedicated training beams to the MS periodically, which also uses a set of dedicated beams to receive them in different directions. As shown in Fig. 2, one transmission frame is assumed to consist of $T$ transmission blocks, and each transmission block is made up of $M_{\rm T}$ symbol durations, where the first $M_{\rm C}$ or $M_{\rm B}$ symbol durations are used for beam training and the rest are for data communication. For the conventional beam training protocol shown in Fig. 2(a), as no priori information is used, in each beam training period the BS and MS consume a fixed number of $M_{\rm C}$ beams in $M_{\rm C}$ symbol durations to estimate the AoAs and AoDs. While for the adopted beam training protocol depicted in Fig. 2(b), traditional channel estimation is performed in the first transmission block, since no priori information can be exploited at this block. In each of the subsequent transmission blocks, $M_{\rm B}$ training beams in $M_{\rm B}$ symbol durations are selected to execute beam tracking based on the previous estimate333When one beam training period has finished, the MS feeds back the estimated AoAs and AoDs to the BS for the subsequent data transmission and the next beam tracking procedure. and the priori transition probabilities, which will be introduced later. Moreover, as $M_{\rm B}<M_{\rm C}$, a larger fraction of time can be left for data communication, leading to a higher throughput. Finally, since channel estimation techniques have been widely investigated in the existing literature [7, 8, 10, 9, 11], we only consider the beam tracking strategy commencing from the second beam training period.

In the rest of the paper, we focus on the $\tau$-th beam training period unless otherwise specified, $\exists\tau\in\{2,\cdots,T\}$, such that $\theta_{\ell}^{[\tau-1]}=\bar{\theta}_{k_{0}}$ and $\theta_{\ell}^{[\tau]}=\bar{\theta}_{k_{1}}$ are the previous and the current AoAs of the $\ell$-th path. In accordance with (5), the previous AoD and the current AoD of the $\ell$-th path are represented by $\vartheta_{\ell}^{[\tau-1]}=\bar{\vartheta}_{i_{0}}$ and $\vartheta_{\ell}^{[\tau]}=\bar{\vartheta}_{i_{1}}$, respectively. Moreover, while below we only consider a single-path channel model and drop the subscript $\ell$ for the sake of convenience, the proposed beam tracking strategy can be readily extended to the multi-path scenario, which is discussed at the end of Section IV.

III Special Case: Unitary Codebook Matrices

In this section, to gain some insights, we follow [28] and assume that $X_{\rm T}=N_{\rm T}$ and $X_{\rm R}=N_{\rm R}$, such that the two beam codebook matrices ${\bm{A}}_{\rm T}$ and ${\bm{A}}_{\rm R}$ become two discrete Fourier transformation (DFT) matrices.

III-A Closed-Form ASTP

Similar to [26, 27, 28], we pick the Tx and Rx training beams from ${\bf A}_{\rm T}$ and ${\bf A}_{\rm R}$, respectively, and when ${\bf f}^{[m]}={\bf a}_{\rm T}(\bar{\vartheta}_{i})$ and ${\bf w}^{[m]}={\bf a}_{\rm R}(\bar{\theta}_{k})$ are chosen at the $m$-th symbol duration or measurement, $\exists m=1,\cdots,M_{\rm B}$, the received symbol can be given by

[TABLE]

where $\gamma=\sqrt{PN_{\rm T}N_{\rm R}}$ and (a) is due to the fact that ${\bf A}_{\rm T}$ and ${\bf A}_{\rm R}$ are two unitary matrices.

In order to improve the received signal power in a specific direction, its associated Tx-Rx beam pair is allowed to be used repeatedly. For convenience, we use ${\cal B}_{k,i}$ to denote the Tx-Rx beam pair ${\bf w}={\bf a}_{\rm R}(\bar{\theta}_{k})$ and ${\bf f}={\bf a}_{\rm T}(\bar{\vartheta}_{i})$, and {\cal B}\triangleq\big{\{}{\cal B}_{k,i}\mid 1\leq k\leq X_{\rm R},1\leq i\leq X_{\rm T}\big{\}} is the set consisting of all potential beam pairs. Moreover, the repetition times of ${\cal B}_{k,i}$ during one beam training period is denoted by $\lambda_{k,i}$, and the corresponding received symbols are expressed as ${\bf y}_{k,i}[1],\cdots,{\bf y}_{k,i}[\lambda_{k,i}]$, respectively. By adding up the $\lambda_{k,i}$ received symbols, we obtain that

[TABLE]

Recall that there are in total $X\triangleq X_{\rm T}\times X_{\rm R}$ distinct beam pairs in ${\cal B}$. For notational simplicity, we redefine ${\cal B}_{n}$ to represent the $n$-th beam pair, $\forall n=1,\cdots,X$, and obviously a one-to-one mapping exists between ${\cal B}_{n}$ and ${\cal B}_{k,i}$, which is denoted by $n=(k\bullet i)_{X_{\rm R}}$. Accordingly, in the rest of the paper, we use $\xi_{n}$ and $\lambda_{n}$ to replace $\xi_{k,i}$ and $\lambda_{k,i}$, respectively.

In the current beam training period, without loss of generality, we denote the selected Tx-Rx training beam pairs by ${\cal B}_{z_{1}},\cdots,{\cal B}_{z_{N}}$, where $z_{1},\cdots,z_{N}\in\{1,\cdots,X\}$ remain to be optimized with $N\leq M_{\rm B}$ since one beam pair might be used repeatedly. Following [33], a power-based estimator is introduced to estimate the updated AoA and AoD for its simplicity, and when $z_{n}\triangleq(a_{n}\bullet c_{n})_{X_{\rm R}}=(k_{1}\bullet i_{1})_{X_{\rm R}}$, the successful estimation probability is given by

[TABLE]

It can be observed from (7) that $|\xi_{z_{n}}|^{2}$ satisfies a non-central chi-squared distribution while $|\xi_{z_{m}}|^{2}$ follows an exponential distribution, $\forall m\neq n$, and therefore we can rewrite (8) as

[TABLE]

where $h(u;\lambda_{z_{n}},|\alpha|^{2})$ is given by

[TABLE]

and $I_{0}(\cdot)$ is the zero-th order modified Bessel function of the first kind. Moreover, we need to integrate $\Gamma_{z_{n},|\alpha|^{2}}$ over the exponential distribution of $|\alpha|^{2}$, which is expressed as

[TABLE]

where \lambda_{\kappa_{1},\kappa_{2},\kappa_{3},-z_{n}}\in\{\lambda_{z_{1}},\cdots,\lambda_{z_{N}}\}\big{\backslash}\{\lambda_{z_{n}}\}.

The (one-step) ASTP can be therefore expressed as

[TABLE]

where $\pi_{z_{n}}$ denotes the transition probability from \Big{\{}\bar{\theta}_{k_{0}},\bar{\vartheta}_{i_{0}}\Big{\}} to \Big{\{}\bar{\theta}_{a_{n}},\bar{\vartheta}_{c_{n}}\Big{\}}, which is directly calculated from (4) and (5).

III-B Problem Formulation

In this paper, we aim to seek the optimal Tx-Rx training beam pairs ${\cal B}_{z_{1}},\cdots,{\cal B}_{z_{N}}$ and their associated repetition times ${\bm{\lambda}}=[\lambda_{z_{1}},\cdots,\lambda_{z_{N}}]^{\rm T}$ that can maximize the ASTP, given the total number of pilot symbol durations $M_{\rm B}$. To this end, the following optimization problem is formulated:

[TABLE]

It is observed that (P1) is an integer nonlinear programming problem, which is in general NP-hard. In the subsequent theorem, we show that the domain of (P1) can be substantially reduced.

Theorem 1**.**

If the $X$ possible beam pairs in $\cal B$ are sorted in a descending order according to their associated transition probabilities, $\{{\cal B}_{1},\cdots,{\cal B}_{X}\}\rightarrow\{{\cal B}_{s_{1}},\cdots,{\cal B}_{s_{X}}\}$, in order to achieve the optimal ASTP, the numbers of used Tx-Rx beam pairs should satisfy $\lambda_{s_{1}}\geq\cdots\geq\lambda_{s_{X}}$.

Proof:

Refer to Appendix A. ∎

Thanks to Theorem 1, in the following we can use $\{{\cal B}_{s_{1}},\cdots,{\cal B}_{s_{N}}\}$ to replace $\{{\cal B}_{z_{1}},\cdots,{\cal B}_{z_{N}}\}$ and rewrite (P1) as

[TABLE]

It is still challenging to handle (P2) due to the complicated structure of its objective function in (17), and therefore we first simplify it into a more tractable form.

Lemma 1**.**

If ${\cal B}_{s_{1}},\cdots,{\cal B}_{s_{N}}$ are used in the beam training period, the ASTP with the power-based estimator is lower bounded by

[TABLE]

and upper bounded by

[TABLE]

Furthermore, $\bar{\Gamma}_{1}(\bm{\lambda})$ can be approximated by

[TABLE]

or equivalently

[TABLE]

where $\displaystyle r_{0}=\frac{PN_{\rm T}N_{\rm R}\sigma_{\alpha}^{2}}{\sigma_{0}^{2}}$ and $\displaystyle{\cal F}(N-1)=\sum\limits_{n=1}^{N-1}\frac{1}{n}$.

Proof:

Refer to Appendix B. ∎

III-C Iterative Nonlinear Branch-and-Bound Algorithm

Since ${\cal F}(N-1)<N-1$ and $M_{\rm B}-\lambda_{s_{n}}=\sum\nolimits_{m\neq n}\lambda_{s_{m}}>N-1$, it is observed that

[TABLE]

In general, while we can pick any of $\bar{\Gamma}_{1}^{\rm lb}(\bm{\lambda})$, $\displaystyle\bar{\Gamma}_{1}^{\rm ub}(\bm{\lambda})$ and $\bar{\Gamma}_{1}^{\rm apx}(\bm{\lambda})$ to replace $\bar{\Gamma}_{1}(\bm{\lambda})$ as the new objective function, in the subsequent sections, we use $\bar{\Gamma}_{1}^{\rm apx}(\bm{\lambda})$ and construct a new optimization problem, given by

[TABLE]

Though the exact value of $N$ in (25) is unknown, it has at most $M_{\rm B}$ cases, i.e., $N=1,\cdots,M_{\rm B}$, and we can thus decompose (P3) into $M_{\rm B}$ subproblems, with each one corresponding to a specific $N$. By solving these subproblems, the optimal solution to (P3) can be obtained. Moreover, since these subproblems are concave I-NLPs as demonstrated in Lemma 2, we can apply the nonlinear branch-and-bound method [39] to solve them optimally.

Lemma 2**.**

For a specific $N$, when we relax the integer variable $\lambda_{s_{n}}$ to a real variable $\tilde{\lambda}_{s_{n}}$, $\forall n=1,\cdots,N$, $\bar{\Gamma}_{1}^{\rm apx}$ becomes a concave function with respect to $\tilde{\lambda}_{s_{1}},\cdots,\tilde{\lambda}_{s_{N}}$.

Proof:

Refer to Appendix C ∎

A closer observation of $\bar{\Gamma}_{1}^{\rm apx}$ shows that we may not need to solve all the $M_{\rm B}$ subproblems of (P3). To be specific, when

[TABLE]

where $\bm{\lambda}_{N}=[\lambda_{s_{1}},\cdots,\lambda_{s_{N}}]^{\rm T}$ denotes a solution to the $N$-th subproblem, $\exists N\in\{2,\cdots,M_{\rm B}\}$, and $\Delta$ is a constant, we can see that

[TABLE]

where $K=2,\cdots,N-1$. In other words, when $\bar{\Gamma}_{1}^{\rm apx}(\bm{\lambda}_{N})<\sum\nolimits_{n=1}^{N}\pi_{s_{n}}\leq\Delta$, the maximum value of $\bar{\Gamma}_{1}^{\rm apx}(\bm{\lambda}_{K})$ will be less than $\Delta$ as well. By using this property, we can skip some subproblems to reduce the computational cost of (P3). For clarity, the proposed iterative N-BB algorithm has been summarized in Algorithm 1.

Seen from Algorithm 1, when $N=M_{\rm B}$, the associated subproblem has only one solution, given by $\{\lambda_{s_{n}}\}_{n=1}^{M_{\rm B}}=1$, and its objective value $\bar{\Gamma}_{1}^{\rm apx}(\bm{\lambda}_{M_{\rm B}})$ is taken as a temporary lower bound of $\bar{\Gamma}_{1}^{\rm apx}(\bm{\lambda})$. Next, we set $N=M_{\rm B}-1$, which corresponds to $\lambda_{s_{1}}=2$, $\{\lambda_{s_{n}}\}_{n=1}^{M_{\rm B}-1}=1$, and $\lambda_{s_{M_{\rm B}}}=0$. We evaluate $\bar{\Gamma}_{1}^{\rm apx}(\bm{\lambda}_{M_{\rm B}-1})$ and compare it with the current lower bound $\bar{\Gamma}_{1}^{\rm apx}(\bm{\lambda}_{M_{\rm B}})$. If $\bar{\Gamma}_{1}^{\rm apx}(\bm{\lambda}_{M_{\rm B}-1})$ is larger than $\bar{\Gamma}_{1}^{\rm apx}(\bm{\lambda}_{M_{\rm B}})$, we set $\bar{\Gamma}_{1}^{\rm apx}(\bm{\lambda}_{M_{\rm B}-1})$ as the new lower bound of $\bar{\Gamma}_{1}^{\rm apx}(\bm{\lambda})$. Otherwise, the current lower bound remains unchanged. We then consider $N=M_{\rm B}-2$ by letting $\lambda_{s_{M_{\rm B}-1}},\lambda_{s_{M_{\rm B}}}=0$, and $\{\lambda_{s_{n}}\}_{n=1}^{M_{\rm B}-2}\geq 1$. Optimize $\bar{\Gamma}_{1}^{\rm apx}(\bm{\lambda}_{M_{\rm B}-2})$, compare its maximal value with the current lower bound and update the lower bound if applicable. This procedure is repeated until the temporary lower bound of $\bar{\Gamma}_{1}^{\rm apx}(\bm{\lambda})$ cannot improve. The current lower bound is the global maximal value of $\bar{\Gamma}_{1}^{\rm apx}(\bm{\lambda})$, and the corresponding solution is the optimal solution to (P3).

III-D A Low-Complexity Solution

Though the iterative N-BB algorithm is able to solve (P3) optimally, its computational cost is high. In order to reduce the complexity, a suboptimal solution to (P3) is also provided by exploiting the KKT conditions. To be more specific, we first relax the $N$-th integer subproblem to a convex nonlinear optimization problem with respect to $\tilde{\lambda}_{s_{1}},\cdots,\tilde{\lambda}_{s_{N}}$, given by

[TABLE]

The associated Lagrangian of (P4) is then expressed as

[TABLE]

where $\mu_{0},\mu_{1},\cdots,\mu_{N}$ are the Lagrange multipliers, and the corresponding KKT conditions are given by

[TABLE]

It is worth mentioning that when we compute the partial derivative of $\cal L$ with respect to $\tilde{\lambda}_{s_{n}}$ in (33a), we approximate $\tilde{\lambda}_{s_{n}}^{2}r_{0}+\tilde{\lambda}_{s_{n}}$ by $\tilde{\lambda}_{\tau_{n}}^{2}r_{0}$ since the SNR at the MS $r_{0}=\frac{PN_{\rm T}N_{\rm R}\sigma_{\alpha}^{2}}{\sigma_{0}^{2}}\gg 1$. By solving the above KKT conditions, we can obtain that $\tilde{\lambda}_{s_{n}}^{\star}=\max\{1,b_{n}\}$, where

[TABLE]

and $\mu_{0}$ is chosen to guarantee that $\tilde{\lambda}_{s_{1}}^{\star}+\cdots+\tilde{\lambda}_{s_{N}}^{\star}=M_{\rm B}$. In general, these obtained solutions $\tilde{\lambda}_{s_{1}}^{\star},\cdots,\tilde{\lambda}_{s_{N}}^{\star}$ are not integers, and therefore we need to truncate them to satisfy the integer requirement. Specifically, we can round off $\tilde{\lambda}_{s_{1}}^{\star},\cdots,\tilde{\lambda}_{s_{N}}^{\star}$ to obtain an integer solution $\lambda_{s_{1}}^{\star},\cdots,\lambda_{s_{N}}^{\star}$. However, due to the rounding off operation, the constraint $\lambda_{s_{1}}^{\star}+\cdots+\lambda_{s_{N}}^{\star}=M_{\rm B}$ may be slightly violated. To tackle this problem, when $\lambda_{s_{1}}^{\star}+\cdots+\lambda_{s_{N}}^{\star}=M>M_{\rm B}$ and $K=M-M_{\rm B}$, we calculate $d_{n}=\lambda_{s_{n}}^{\star}-\tilde{\lambda}_{s_{n}}^{\star}$, $\forall n=1,\cdots,N$. If $d_{p_{1}}>d_{p_{2}}>\cdots>d_{p_{N}}$, we let $\lambda_{s_{p_{k}}}^{\star}=\lambda_{s_{p_{k}}}^{\star}-1$, $\forall k=1,\cdots,K$. On the other hand, when $\lambda_{s_{1}}^{\star}+\cdots+\lambda_{s_{N}}^{\star}=M<M_{\rm B}$ and $K=M_{\rm B}-M$, we compute $d_{n}=\tilde{\lambda}_{s_{n}}^{\star}-\lambda_{s_{n}}^{\star}$, $\forall n=1,\cdots,N$, and if $d_{p_{1}}>d_{p_{2}}>\cdots>d_{p_{N}}$, we let $\lambda_{s_{p_{k}}}^{\star}=\lambda_{s_{p_{k}}}^{\star}+1$, $\forall k=1,\cdots,K$. The other subproblems can be solved similarly and we summarize the whole procedure in Algorithm 2 for clarity.

Remark 1: It is worth mentioning that the proposed beam pair allocation strategy is deduced by assuming accurate knowledge of $\theta^{[\tau-1]}$ and $\vartheta^{[\tau-1]}$. If the previous estimates of $\theta^{[\tau-1]}$ and $\vartheta^{[\tau-1]}$ are inaccurate, the proposed strategy might worsen the $\tau$-th and the subsequent beam tracking procedures, as demonstrated in Section V. To alleviate the error propagation phenomenon incurred by the proposed beam allocation strategy, when \big{|}\xi_{z_{n}}^{[\tau-1]}\big{|}^{2}-\big{|}\xi_{z_{m}}^{[\tau-1]}\big{|}^{2}<\Omega, where \big{|}\xi_{z_{n}}^{[\tau-1]}\big{|}^{2} and \big{|}\xi_{z_{m}}^{[\tau-1]}\big{|}^{2} are the two largest received signal powers in the $(\tau-1)$-th beam training period and $\Omega$ is a pre-defined threshold, we employ uniform allocation strategy instead of invoking the proposed allocation strategy in the $\tau$-th beam training period.

IV Extension to General Non-Orthogonal Codebook Matrices

In this section, we consider a more general scenario in which $N_{\rm T}<X_{\rm T}$ and $N_{\rm R}<X_{\rm R}$, such that the two beam codebook matrices ${\bf A}_{\rm T}$ and ${\bf A}_{\rm R}$ are not DFT matrices any longer. In this case, even if the adopted Tx-Rx beam pair is not perfectly aligned with the actual AoA and AoD, the MS can still receive the pilot symbol with a considerable beamforming gain, which is different from the orthogonal case.

IV-A Power-Based Estimator

As before, we pick the columns from ${\bf A}_{\rm T}$ and ${\bf A}_{\rm R}$ as the Tx and Rx beams in each beam training period, and when ${\bf f}^{[m]}={\bf a}_{\rm T}(\bar{\vartheta}_{i})$ and ${\bf w}^{[m]}={\bf a}_{\rm R}(\bar{\theta}_{k})$ are chosen at the $m$-th symbol duration, the received symbol in (6) becomes

[TABLE]

where ${\tilde{\nu}}_{i,i_{1}}$ and $\nu_{k,k_{1}}$ are respectively expressed as

[TABLE]

It is easy to see that when $i=i_{1}$ and $k=k_{1}$, $\tilde{\nu}_{i,i_{1}}=1$ and $\nu_{k,k_{1}}=1$, and (35) reduces to (6). Following the previous description, the beam pairs ${\cal B}_{z_{1}},\cdots,{\cal B}_{z_{N}}$ with repetition times $\lambda_{z_{1}},\cdots,\lambda_{z_{N}}$ are used in the current beam training period, where $z_{p}=(a_{p}\bullet c_{p})_{X_{\rm R}}$, $p=1,\cdots,N$. Based on this assumption, the distribution of $\xi_{z_{p}}$, which is first defined in (7), is given by

[TABLE]

The successful estimation probability shown in (9) now becomes

[TABLE]

where the non-central chi-squared distribution $h(u;\lambda_{z_{n}},|\alpha|^{2})$ is already given by (10) and $Q_{1}$ is the first-order Marcum $Q$-function. Furthermore, by integrating $\Gamma_{z_{n},|\alpha|^{2}}$ over the exponential distribution of $|\alpha|^{2}$, we obtain

[TABLE]

It is challenging to derive a closed-form expression for $\Gamma_{z_{n}}$ in (40), which is therefore calculated via the numerical method. Finally, similar to (12), the ASTP can be given by

[TABLE]

IV-B OMP-Based Estimator

In general, the power-based estimator performs poorly due to the non-negligible inter-beam interference, which actually can be used to improve the beam tracking performance. To show this, we first rewrite the channel in (3) as

[TABLE]

where $\bf V$ is termed the beamspace channel representation of $\bf H$ in some literature [30]. Based on the on-grid assumption of the AoAs and AoDs, the matrix $\bf V$ is sparse with only $L$ nonzero elements, e.g., ${\bf V}[k_{1},i_{1}]$ in current transmission block.

By introducing the beamspace channel $\bf V$, we can transform the AoA and AoD estimation problem into a CS problem. To be more specific, we rewrite (6) as

[TABLE]

where ${\bf v}=\text{vec}({\bf V})$. In (IV-B), we have applied the property $\text{vec}({\bf ABC})=({\bf C}^{\rm T}\otimes{\bf A})\text{vec}({\bf B})$ in (a) and (b) is due to (36) and (37). Recall that ${\cal B}_{z_{1}},\cdots,{\cal B}_{z_{N}}$ are used, and if we collect their corresponding observations into a vector, denoted by

[TABLE]

we can obtain a CS problem, given by

[TABLE]

where {\bf n}=\Big{[}n^{[1]},\cdots,n^{[M_{\rm B}]}\Big{]}^{\rm T}. The sensing matrix $\bf A$ is written as

[TABLE]

where ${\bf a}_{z_{n}}$ is shown in (IV-B) through replacing $k$ and $i$ by $a_{n}$ and $c_{n}$, respectively, $n=1,\cdots,N$.

The MS can use the OMP algorithm [34, 35] to estimate the nonzero element in $\bf v$, which is also termed the support of $\bf v$ [24]. By adopting the OMP algorithm, the estimated support of $\bf v$ is given by

[TABLE]

where $\bm{\xi}\triangleq{\bf A}^{\rm H}{\bf y}$. This OMP-based estimator can be viewed as an improvement of the previous power-based estimator, and the successful estimation probability conditioned on $|\alpha|^{2}$ is given by

[TABLE]

where $n_{1}=(k_{1}\bullet i_{1})_{X_{\rm R}}$. Moreover, we need to integrate (47) over the exponential distribution of $|\alpha|^{2}$, given by

[TABLE]

and the ASTP of this OMP-based estimator is expressed as

[TABLE]

IV-C Problem Formulation

In this section, we replace $\bar{\Gamma}_{1}$ by $\bar{\Gamma}_{3}$ and construct a new optimization problem, given by

[TABLE]

However, different from $\bar{\Gamma}_{1}$ in (12), $\bar{\Gamma}_{3}$ has no closed-form expression, and therefore we derive a closed-form lower bound for $\Gamma_{n_{1}}$, which is shown in Lemma 3.

Lemma 3**.**

If the Tx-Rx beam pairs ${\cal B}_{z_{1}},\cdots,{\cal B}_{z_{N}}$ are used in the beam training period, where $z_{m}=(a_{m}\bullet c_{m})_{X_{\rm R}}$, $\forall m=1,\cdots,N$, and their repetition times are respectively denoted by $\lambda_{z_{1}},\cdots,\lambda_{z_{N}}$, the lower bound of $~{}\Gamma_{n_{1},|\alpha|^{2}}$ is expressed as

[TABLE]

where $n=(k\bullet i)_{X_{\rm R}}$ and $n_{1}=(k_{1}\bullet i_{1})_{X_{\rm R}}$. Accordingly, the lower bound of $~{}\Gamma_{n_{1}}$ is written as

[TABLE]

where $A_{n,n_{1}}$ and $B_{n,n_{1}}$ are respectively given by

[TABLE]

[TABLE]

Proof:

Refer to Appendix D. ∎

According to Lemma 3, the lower bound for the ASTP when using the OMP-based estimator can be written as

[TABLE]

Since we cannot directly optimize $\bar{\Gamma}_{3}$ in (50), we use $\bar{\Gamma}_{3}^{\rm lb}$ as the new objective function, and the associated optimization problem becomes

[TABLE]

Considering the complicated expression for $\bar{\Gamma}_{3}^{\rm lb}$, it is still very challenging to solve (P6) analytically. Therefore, when $M_{\rm B}$ is small, we propose to directly search for the $M_{\rm B}$ Tx-Rx training beam pairs from the two beam codebook matrices ${\bf A}_{\rm T}$ and ${\bf A}_{\rm R}$ that can maximize (55). On the other hand, when $M_{\rm B}$ is large, this exhaustive search method is prohibitive since there are $X^{M_{\rm B}}$ possible solutions. In this situation, we can exploit a heuristic algorithm, such as the differential evolution [38], to obtain a promising solution to (P6), whose details are omitted for brevity.

Remark 2: Both the power-based estimator and the OMP-based estimator can be readily extended to the multi-path scenario as in [8, 27]. To be more specific, we can estimate the $L\geq 2$ paths in an iterative fashion, and only one path is estimated via using the power-based estimator or the OMP-based estimator at each iteration. Following the idea of successive interference cancelation, the contribution of the paths that have been estimated in the previous iterations is subtracted from the received sequence before finding new paths. In addition, since jointly optimizing the $M_{\rm B}$ Tx-Rx beam pairs is very challenging when multiple paths exist, we optimize the $M_{\rm B}/L$ training beams pairs for each of the $L$ paths via solving an optimization problem similar to (P3) or (P6) separately. Though such a per-path beam allocation strategy seems to be trivial, it can still achieve a favorable beam tracking performance.

V Numerical Results

In this section, we provide numerical results to evaluate the tracking performance of the proposed beam pair allocation strategy for time-varying mmWave MIMO systems. The average tracking error probability (ATEP) is mainly used as the performance metric, which is expressed as $1-\text{ASTP}$. The geometric channel model in (3) is adopted with $L=1$ and $\sigma_{\alpha}^{2}=1$. The SNR is defined as $P/\sigma_{0}^{2}$. To guarantee a favorable angular resolution, we require that $X_{\rm T}=X_{\rm R}=64$. In addition, one transmission frame consists of $T=10$ transmission blocks [24, 28] and the temporal variations of AoA and AoD among these transmission blocks are assumed to follow two discrete Markov processes, described by (4) and (5), respectively. Moreover, the first transmission block is assumed to have exact AoA and AoD knowledge by using the traditional channel estimation algorithms such as [8]. For the remaining 9 blocks, we exploit the previously estimated AoA and AoD to conduct the current beam tracking procedure. Several benchmark methods are also introduced for comparison, which are presented as follows:

Proportional Allocation: The $M_{\rm B}$ pilot symbol durations are distributed to ${\cal B}$ in proportion to their associated transition probabilities; 2. 2.

Uniform Allocation: The $M_{\rm B}$ pilot symbol durations are uniformly distributed to ${\cal B}_{s_{1}},\cdots,{\cal B}_{s_{M_{\rm B}}}$; 3. 3.

Proposed Allocation ES: The $M_{\rm B}$ pilot symbol durations are distributed to ${\cal B}$ via solving (P2) or (P6) with an exhaustive search method; 4. 4.

Proposed Allocation BB: The $M_{\rm B}$ pilot symbol durations are distributed to ${\cal B}$ via solving (P3) with Algorithm 1; 5. 5.

Proposed Allocation KKT: The $M_{\rm B}$ pilot symbol durations are distributed to ${\cal B}$ via solving (P4) with Algorithm 2; 6. 6.

ML-Based Estimator: The time-varying AoA and AoD of a single-path channel are estimated based on the ML criterion [26]; 7. 7.

POMDP Framework: The time-varying AoA and AoD are estimated based on the belief states of the formulated POMDP in [28].

V-A Special Scenario: Orthogonal Tx-Rx Beam Pairs

The ATEP with respect to the training SNR is provided in Fig. 4, where we let $N_{\rm T}=X_{\rm T}$ and $N_{\rm R}=X_{\rm R}$. In this scenario, we see that the proportional allocation strategy outperforms the uniform allocation strategy in the low SNR region, whereas it is inferior to the latter when the training SNR is high. The proposed beam allocation strategy performs better than the two benchmarks at the whole SNR range. In addition, we note that the solution obtained by solving the $M_{\rm B}$ relaxed subproblems with KKT conditions in (P4) can achieve almost the same ATEP as that obtained via the ES method. It is also worth mentioning that the ATEP curve of the iterative N-BB algorithm is not presented in Fig.4 for clarity, since it converges to the curves of the ES and KKT-based methods.

The ATEP with respect to the AoA’s variation speed $\beta$ is shown in Fig. 4, where we can see that the proposed beam allocation strategy still performs better than the other three strategies, though the ATEPs of the four aforementioned beam allocation strategies all deteriorate when $\beta$ increases from 0.1 to 0.5. Moreover, the ATEP gap between the proposed beam allocation strategy and the uniform allocation strategy becomes narrower and narrower. The reason can be explained by observing (34), which demonstrates that the repetition times of each training beam pair is directly determined by its associated transition probability. Since the transition probability for each of the potential directions becomes more uniform when $\beta$ is large, according to (34), the proposed beam allocation strategy is asymptotically close to the uniform allocation strategy.

The ATEP with respect to the total number of training beam pairs is presented in Fig. 6, where we can see that the ATEP of the proposed allocation strategy decreases gradually when $M_{\rm B}$ increases from 30 to 50. The ATEP of the KKT-based algorithm is very close to that of the iterative N-BB algorithm, which achieves almost the same performance as the ES method, demonstrating the validity of Algorithm 1 and Algorithm 2. Moreover, we observe that the ATEP of the uniform beam allocation strategy deteriorates slightly as $M_{\rm B}$ increases. According to (11) and (12), though more directions can be measured when we increase $M_{\rm B}$, the interference in (11) also increases such that $\Gamma_{z_{n}}$ will decrease, and consequently the ATEP of the the uniform beam allocation strategy increases.

The ATEP with respect to the index of transmission block is shown in Fig. 6, in which we observe that the performance of the POMDP framework is quite robust with the process of beam tracking. The ATEPs of the ML-based estimator and the proposed power-based estimator degrade quickly compared to that of the POMDP framework. However, by exploiting the modified beam allocation strategy depicted in Remark 1, it is seen that the ATEP performance of the power-based estimator can be significantly improved, where we set the threshold $\Omega=5$ in the simulations.

V-B General Scenario

The ATEP with respect to the SNR when $N_{\rm T}<X_{\rm T}$ and $N_{\rm R}<X_{\rm R}$ is shown in Fig. 7, where we can immediately see that the proportional allocation strategy cannot work in the whole SNR region. The power-based estimator in (39) is inferior to the OMP-based estimator in (46) which exploits the inter-beam interference to estimate the AoA and AoD. In addition, we also observe that the proposed beam allocation strategy performs much better than the uniform and proportional allocation strategies. Moreover, since two adjacent beams are overlapped in the angular space when $N_{\rm T}<X_{\rm T}$ and $N_{\rm R}<X_{\rm R}$, the BS can still transmit information to the MS with a considerable beamforming gain even if the estimates of the true AoA and AoD are inaccurate. Therefore, we also provide the average beamforming gains444The average beamforming gain is expressed as {\mathbb{E}}\big{[}|{\bf a}^{\rm H}_{\rm R}\big{(}\hat{\theta}^{[\tau]}\big{)}{\bf a}_{\rm R}\big{(}\theta^{[\tau]}\big{)}{\bf a}^{\rm H}_{\rm T}\big{(}\vartheta^{[\tau]}\big{)}{\bf a}_{\rm T}(\hat{\vartheta}^{[\tau]})|^{2}\big{]}, where $\hat{\theta}^{[\tau]}$ and $\hat{\vartheta}^{[\tau]}$ represent the estimated AoA and AoD in the $\tau$-th beam training period. of the four aforementioned strategies in Fig. 7, in which we can observe that the proposed beam allocation strategy is still significantly superior to the other three strategies, especially at the low SNR regime.

The ATEP with respect to $\beta$ is depicted in Fig. 9, where we fix $\tilde{\beta}$ at 0.1. Similar to Fig. 4, while both the ATEPs of the OMP-based estimator with the proposed beam allocation strategy and uniform allocation strategy deteriorate when $\beta$ increases from 0.1 to 0.5, the former strategy still performs better than the latter. Moreover, as $\beta$ increases, we see that the ATEP gap between the two strategies also becomes narrower.

In Fig. 9, the ATEPs of the 10 transmission blocks are presented, where the AoA and AoD keep changing from one transmission block to another. It is seen that the beam tracking performance of the OMP-based estimator is significantly superior to that of the power-based estimator due to the inter-beam interference. The ML-based estimator performs better than the power-based estimator, whereas it is inferior to the OMP-based estimator. Moreover, we also observe that the ATEP of the OMP-based estimator with the proposed allocation strategy is much better than those of the benchmarks.

VI Conclusions

In this paper, we have proposed a new beam pair allocation strategy for mmWave MIMO tracking systems, which enables one Tx-Rx beam pair to be used repeatedly to improve the received signal power at that direction. We have firstly considered a special scenario in which $N_{\rm T}=X_{\rm T}$ and $N_{\rm R}=X_{\rm R}$. In this case, the Tx-Rx beam pairs are orthogonal with each other and the training beam pair sequence design problem can be approximately tackled by solving a set of concave I-NLPs. The obtained closed-form solution shows that the repetition times of each Tx-Rx beam pair is directly determined by its associated transition probability, and one beam pair with a higher transition probability should be used more times than those with lower transition probabilities. In the case of $N_{\rm T}<X_{\rm T}$ and $N_{\rm R}<X_{\rm R}$, we have derived a closed-form lower bound for the ASTP when the OMP-based estimator is used to track the time-varying AoA and AoD, based on which a favorable beam pair allocation strategy is obtained. Our numerical results have validated the superiority of the proposed allocation strategy over the existing methods.

Appendix A

In order to prove Theorem 1, we first show that $\Gamma_{z_{n}}$ is an increasing function with respect to $\lambda_{z_{n}}$, and a decreasing function with respect to $\lambda_{z_{m}}$, $\forall z_{m}\neq z_{n}$. To verify this conclusion, we relax the integer variable $\lambda_{z_{n}}$ to a real variable $\tilde{\lambda}_{z_{n}}$, and derive the partial derivative of $\Gamma_{z_{n}}$ with respect to $\tilde{\lambda}_{z_{n}}$ as

[TABLE]

where $h_{1}(\tilde{\lambda}_{z_{n}})=\tilde{\lambda}_{z_{n}}^{2}\gamma^{2}\sigma_{\alpha}^{2}+\tilde{\lambda}_{z_{n}}\sigma_{0}^{2}$ and $h_{2}(\tilde{\lambda}_{z_{n}})=h^{\prime}_{1}(\tilde{\lambda}_{z_{n}})=2\tilde{\lambda}_{z_{n}}\gamma^{2}\sigma_{\alpha}^{2}+\sigma_{0}^{2}$. It is then observed that ${\partial\Gamma_{z_{n}}}/{\partial\tilde{\lambda}_{z_{n}}}$ is larger than zero, and hence we obtain that $\Gamma_{z_{n}}(\bm{\lambda})$ is an increasing function with respect to $\tilde{\lambda}_{z_{n}}$ or $\lambda_{z_{n}}$. By relaxing $\lambda_{z_{m}}$ to a real variable $\tilde{\lambda}_{z_{m}}$ and computing the partial derivative of $\Gamma_{z_{n}}$ with respect to $\tilde{\lambda}_{z_{m}}$, we can similarly verify that $\Gamma_{z_{n}}(\bm{\lambda})$ is a decreasing function of $\lambda_{z_{m}}$.

Next we prove this theorem. Firstly, when we assume that $\pi_{z_{n}}\geq\pi_{z_{m}}$ but $\lambda_{z_{n}}<\lambda_{z_{m}}$, we can obtain that

[TABLE]

where (a) comes from the fact that $\Gamma_{z_{n}}(\bm{\lambda})$ is an increasing function with respect to $\lambda_{z_{n}}$ and (b) is due to $\exp\left(-\frac{u}{\lambda_{z_{m}}\sigma_{0}^{2}}\right)>\exp\left(-\frac{u}{\lambda_{z_{n}}\sigma_{0}^{2}}\right)$. In addition, \bm{\lambda}_{-z_{m},z_{n}}\triangleq\{\lambda_{z_{1}},\cdots,\lambda_{z_{N}}\}\big{\backslash}\{\lambda_{z_{m}},\lambda_{z_{n}}\} and $G(\bm{\lambda}_{-z_{m},z_{n}})$ is expressed as

[TABLE]

It is observed from (58) that $\Gamma_{z_{n}}<\Gamma_{z_{m}}$ when $\lambda_{z_{n}}<\lambda_{z_{m}}$. By swapping the values of $\lambda_{z_{n}}$ and $\lambda_{z_{m}}$ while keeping $\bm{\lambda}_{-z_{m},z_{n}}$ unchanged, i.e., $\lambda_{z_{n}}^{new}=\lambda_{z_{m}}$, $\lambda_{z_{m}}^{new}=\lambda_{z_{n}}$, $\lambda_{z_{p}}^{new}=\lambda_{z_{p}}$, $\forall p\neq m,n$, we can see that $\displaystyle\Gamma_{z_{n}}^{new}=\Gamma_{z_{m}}$, $\displaystyle\Gamma_{z_{m}}^{new}=\Gamma_{z_{n}}$ and $\displaystyle\Gamma_{z_{p}}^{new}=\Gamma_{z_{p}}$. The variation of the ASTP after swapping the values of $\lambda_{z_{n}}$ and $\lambda_{z_{m}}$ is given by

[TABLE]

which shows that the ASTP $\bar{\Gamma}_{1}(\bm{\lambda})$ can increase when we exchange the values of $\lambda_{z_{n}}$ and $\lambda_{z_{m}}$. By repeating this procedure, we obtain that in order to achieve the maximal $\bar{\Gamma}_{1}$, the numbers of used Tx-Rx beam pairs in $\cal B$ should satisfy $\lambda_{s_{1}}\geq\lambda_{s_{2}}\geq\cdots\geq\lambda_{s_{N}}$.

Appendix B

In order to obtain the lower bound for $\bar{\Gamma}_{1}(\bm{\lambda})$, we first prove the subsequent inequality:

[TABLE]

Induction method is used to prove (61). To be specific, when $p=1$, we can see that (61) is apparently valid. For $p=n-1$ where $n\geq 2$, we assume that (61) is true. For $p=n$, we have

[TABLE]

which proves (61). By using this inequality, we can rewrite $\Gamma_{s_{n}}$ as

[TABLE]

where $r_{0}=\frac{\gamma^{2}\sigma_{\alpha}^{2}}{\sigma_{0}^{2}}=\frac{PN_{\rm T}N_{\rm R}\sigma_{\alpha}^{2}}{\sigma_{0}^{2}}$, and the ASTP can be therefore lower bounded by

[TABLE]

In order to derive the upper bound of $\bar{\Gamma}_{1}(\bm{\lambda})$, we need to prove the subsequent relationship

[TABLE]

By exploiting the binomial theorem, we have

[TABLE]

Integrating the two sides of (66) from 0 to 1, we can obtain

[TABLE]

It is worth noting that the right-hand side (RHS) of (67) is the left-hand side (LHS) of (65). For the LHS of (67), we obtain the subsequent equality

[TABLE]

By comparing (67) and (68), it is seen that (65) has been proven. Next we consider $\Gamma_{s_{n}}$, which is

[TABLE]

where (a) follows from the fact that $1-\exp\left(-\frac{u}{\lambda_{s_{m}}\sigma_{0}^{2}}\right)$ is a decreasing function of $\lambda_{s_{m}}$, and therefore we use $\lambda_{s_{m}}=1$ to obtain an upper bound for $\Gamma_{s_{n}}$, and (b) follows from the fact that $r_{0}\gg 1$. By exploiting (69), the ASTP can be upper bounded by

[TABLE]

Moreover, we can further approximate $\Gamma_{s_{n}}$ as

[TABLE]

where $\lambda_{\rm avg}=\frac{\sum\nolimits_{m\neq n}\lambda_{s_{m}}}{N-1}=\frac{M_{\rm B}-\lambda_{s_{n}}}{N-1}$. The ASTP is therefore approximated by

[TABLE]

which completes the proof of Lemma 1.

Appendix C

To prove that $\bar{\Gamma}_{1}^{\rm apx}(\bm{\lambda})$ is a concave function with respect to $\tilde{\lambda}_{s_{1}},\cdots,\tilde{\lambda}_{s_{N}}$, we compute the first, second and third derivatives of $f(\tilde{\lambda}_{s_{n}})\triangleq(M_{\rm B}-\tilde{\lambda}_{s_{n}})/(\tilde{\lambda}_{s_{n}}^{2}r_{0}+\tilde{\lambda}_{s_{n}})$ with respect to $\tilde{\lambda}_{s_{n}}$, $\forall n=1,\cdots,N$, given by

[TABLE]

Recall that $\tilde{\lambda}_{s_{n}}\leq M_{\rm B}$, and we can observe that $df^{3}/d\tilde{\lambda}_{s_{n}}^{3}<0$, which means that $df^{2}/d\tilde{\lambda}_{s_{n}}^{2}$ is a decreasing function with respect to $\tilde{\lambda}_{s_{n}}$. Due to

[TABLE]

we can conclude that $f(\tilde{\lambda}_{s_{n}})$ is a convex function with respect to $\tilde{\lambda}_{s_{n}}$, and therefore $-f(\tilde{\lambda}_{s_{n}})$ is a concave function. By noting that $\pi_{s_{n}}>0$, $\bar{\Gamma}_{1}^{\rm apx}$ is a concave function with respect to $\tilde{\lambda}_{s_{1}},\cdots,\tilde{\lambda}_{s_{N}}$, since a nonnegative weighted sum of concave functions is concave [40]. This completes the proof of Lemma 2.

Appendix D

Recall that each row of the sensing matrix $\bf A$ can be expressed as

[TABLE]

In accordance with our previous definition $n=(k\bullet i)_{X_{\rm R}}$, the $n$-th column of $\bf A$ is written as

[TABLE]

and $\bm{\xi}[n]$ is therefore given by

[TABLE]

where $p=\lambda_{z_{1}}+\cdots+\lambda_{z_{m-1}}$. It is worth mentioning that while the “noise” terms in $\{\bm{\xi}[n]\}_{n=1}^{X}$ are correlated with each other, we ignore their correlations for tractability. The distribution of $\bm{\xi}[n]$ conditioned on $|\alpha|^{2}$ is therefore expressed as

[TABLE]

Typically, when $n=n_{1}$ or $k=k_{1}$ and $i=i_{1}$, we can write $\bm{\xi}[n_{1}]$ as

[TABLE]

By using the union bound, we can rewrite (47) as

[TABLE]

where $A_{n,n_{1}}$ and $B_{n,n_{1}}$ are respectively given by

[TABLE]

[TABLE]

By using the result of [36] and after some mathematical manipulations, we can obtain the integral expression of the Marcum $Q_{1}$ function over the exponential distribution of $|\alpha|^{2}$ as

[TABLE]

Moreover, by exploiting the integral identity [37], we obtain

[TABLE]

Finally, by integrating (82) over the exponential distribution of $|\alpha|^{2}$, we obtain

[TABLE]

which completes the proof of Lemma 3.