Predictive precoding based on the Grassmannian manifold for UAV-enabled cache-assisted B5G communication systems

The unmanned aerial vehicle (UAV) can extend the network coverage and improve the system throughput for 5th generation (5G) communication systems; hence, it receives a lot of attention recently. This paper considers the problem of channel predictive precoding for UAV-enabled cache-assisted B5G multi-input multi-output (MIMO) systems. A novel channel precoder predictor is proposed, in which the prediction is conducted on a non-linear vector space—Grassmannian manifold. The predictor at the receiver utilizes the current and previous channel matrices to solve the precoder at the next time and then feeds it back to the transmitter for precoding. More specifically, two sub-matrices are extracted from the channel right singular matrices and modeled as two points on the Grassmannian manifold. Then, the geodesic between the two points is conducted. Unlike the conventional method in which the tangent vector at the previous point is parallel transported along the geodesic, we predict the next point by use of the geodesic equation directly. We analyze the computational complexity of the proposed method and demonstrate the superiority of the proposed method by comparing with the conventional one. Besides, we adopt a general Ricean channel model in the UAV MIMO system, where both the Kronecker model and Jake’s model are incorporated. The effects of various channel model parameters on the system performance in terms of the chordal error of channel predictor and the optimum step are thoroughly investigated.

the capacity of 5G systems and user equipments (UEs) as well. Therefore, it has received a lot of attention by the industry and academia.
UAV communications usually adopt the mmWave band, since the multi-input multioutput (MIMO) technology is easy to use in this band [8,9]. The MIMO system, equipped with multiple antennas at both the transmitter and receiver, makes use of spatial resources adequately and has high data throughput and reliable communication quality. Channel tracking/prediction is always a research hot spot in MIMO systems. Channel state information (CSI) is indispensable for the data decoding at the receiver, which is often obtained by channel estimation. Usually, it is completed with the aide of the pilot signal sent by the transmitter frequently. However, by doing so, quite a part of the channel resources are occupied by the pilot signal, especially for fast fading channels. To reduce the length of the pilot signal on the whole, channel tracking or prediction was proposed by some literature [10]. These methods include the Kalman filter-based method [10], extended Kalman filter-based method [11], sequential Monte Carlo filter-based method [12], and particle filter-based method [13]. Recently, there are some works on the channel tracking in UAV MIMO systems. In [14], a channel tracking method for UAV MIMO communication systems was proposed and investigated, in which the method explores the characteristics of time-varying UAV channels with the beam squint effect. In [15], to improve the quality of the UAV navigation, the authors designed a channel tracking algorithm for its flight control system, where the time-varying spatial channel is characterized by a 3D geometry-based channel model and the algorithm incorporates the outputs of multiple sensors in order to reduce the training overhead and energy consumption.
On the other front, channel predictive information can be feeded back to the transmitter for preprocessing in order to improve the system performance in terms of the capacity or the bit-error-rate (BER) [16]. For instance, in [16], the authors proposed a subspace-based channel tracking scheme for precoded MIMO orthogonal frequency division multiplexing (MIMO-OFDM) systems. The predictive CSI was used for precoding at the transmitter to improve the system throughput. With the development of UAV communications, it can be foreknown that there will be a large amount of data to be transmitted interactively between UAVs and ground base stations/terminals. For instance, with the development of forestry informatization, forestry resource management needs more data support. The UAV is suitable for collecting forestry data and sending them back to the ground control center. Besides, the UAV can be used for data offloading of possibly overloaded cellular base stations in hot spots such as sports venues and cinemas. It can also be used for periodic data distribution/collection in large-scale Internet of Things networks. Large amount of data transmission requires temporary storage space. The caching technique pre-stores the data during the off-peak traffic and improves 5G system performance in many aspects [17][18][19]. With the aid of cache, it is convenient for the UAV to store data for its own use or forward data to other UAVs or ground terminals [20][21][22]. Clearly, the feedback information of channel prediction can also be used for UAV MIMO systems to improve the quality of data transmission.
Most previous work above is conducted based on the linear vector space, i.e., Euclidean space, which can be viewed as a linear manifold. Recently, with the rise of information geometry, a few channel tracking methods based on the non-linear manifold were proposed and investigated [23][24][25][26]. For instance, in [23], the authors first modeled the (2020) 2020:128 Page 3 of 18 column spaces spanned by the right singular matrix of the MIMO channel as a point of Grassmannian manifold, and then proposed a method to track the movement of the point as well as an adaptive codebook for precoding. In [26], for the time-varying MIMO system, the authors proposed a predictive quantizer for the eigenvectors of the Gramian matrix that is created from the channel matrix, in which the quantizer operates on the compact Stiefel manifold. In this paper, we mainly predict a sub-matrix of the channel right singular matrix, which is created by selecting a few columns of the channel right singular matrix. The sub-matrix has more rows than columns and can be modeled as a point of Grassmannian manifold. When the channel is time varying, it is reasonable to assume that this sub-matrix moves along the geodesic within a short period of time. Therefore, we intend to use the theory and results of Grassmannian manifold in the channel prediction.
For the prediction methods based on Grassmannian manifold, a key parameter in the design is the choice of the search step. However, most previous work did not investigate the parameter thoroughly. Furthermore, considering the complexity of the practical UAV MIMO channel, this article incorporates both the time correlation and spatial correlation of the channel and studies the effects of various system parameters on the chordal error performance. The contributions of this article are briefly listed as follows: 1 We propose a novel channel predictive precoding method which is based on the geodesic equation on the Grassmannian manifold and analyze its computational complexity as well as the conventional one. Compared with the conventional method, although the complexity of the two methods is the same, the proposed method has better prediction performance in terms of the channel chordal error. 2 We adopt a general Ricean channel model which incorporates the time and spatial correlation by combining the Kronecker model and Jake's model. The effects of various system parameters on the system performance in terms of the chordal error of channel predictor as well as the choice of optimal step are thoroughly investigated. 3 IntheAppendix, we prove a corollary which is inferred from the geodesic equation and is helpful to verify the correctness of the geodesic equation from two points on the Grassmannian manifold.
The rest of this paper is organized as follows. The closed-loop UAV-enabled MIMO system model is described in Section 2. The process of channel precoding prediction and the corresponding problem of interest are illustrated in Section 3. The effect of various system parameters on the performance of the prediction algorithm is presented in Section 4, followed by conclusions in Section 5.
Notations: Boldface lowercase letters denote vectors, and boldface uppercase letters denote matrices. The notation E(.) represents the statistical expectation; C m×n and U m×n represent the space consisting of m×n complex matrices and the space consisting of m×n orthogonal complex matrices, respectively. For a vector x, x ∼ CN (μ, R) represents that x follows a complex Gaussian distribution with mean μ and covariance matrix R; for a matrix X, the notations X 1/2 , X H , and Tr(X) denote its square root, Hermitian transpose, and trace, respectively; besides, I m is an m × m identity matrix, I m,n with n < m is created by selecting the first n columns of I m , and U m is an m × m unitary matrix.

The UAV MIMO channel
As shown in Fig. 1, we consider a UAV-enabled MIMO system, where the base station (BS) and the UAV are equipped with N T and N R antennas, respectively [27,28]. Both the BS and UAV are equipped with cache storage units. With the aid of cache, the UAV can store data conveniently for its own use or forward them to other UAVs/ground terminals at the next time. Figure 1a is the original system diagram, from which we see that there exists the line of sight (LOS) path between BS and UAV; due to the reflections of high buildings, trees, or others, there also exist non-LOS (NLOS) paths. Hence, the channel between BS and UAV is modeled as a Ricean fading channel. For the UAV MIMO channel H(t) ∈ C N R ×N T , we not only consider its time variation, but also the spatial correlation. In [29], a spatial-temporal correlated channel model was proposed; however, only transmitter correlation was considered. Herein, we also incorporate the receiver correlation of the channel. Consequently, the UAV MIMO channel is modeled as: In Eq. (1), the NLOS path H NLOS (t) is further expressed as: where θ T and θ R denote the transmitter and receiver correlation matrix of the channel, respectively; K Rice is the Ricean factor. H ω (t) = h ij (t) ∈ C N R ×N T and it is a random matrix that has i.i.d complex Gaussian random variables with zero mean and unit variance. The autocorrelation of each entry follows Jake's model and it is given by: in which the function J 0 (.) is the well-known Bessel function of the first kind, f d is the Doppler shift, for the NLOS path, which is due to the motions of the UAV, scatters, etc. Assuming that both BS and the UAV adopt the linear antenna array, the LOS path H LOS (t) is expressed as [30]: where f d is Doppler shift for the LOS path, which is due to the motion of the UAV; d R and d T are the antenna spacing at BS and the UAV, respectively; λ is the wavelength; and α and β are the angles of arrival and departure, respectively.
Note that the channel defined in Eq. (1) is assumed to be flat fading. If the paths are resolvable, we can incorporate the OFDM technique and transform the frequency selective channel into a flat fading one.
Remark: Concerning the angles of arrival and departure, the definitions are explained as follows. Take the angle of arrival as an example. As in Fig. 2, a plane electromagnetic wave arrives at the uniform linear antenna array. We adopt the definition of [30] and define α as the direction of arrival (DOA), i.e., the angle between the incident direction of the electromagnetic wave and the x-axis. Consequently, the time delay between the two adjacent antennas is d R cos(α)/λ, which determines the form of the steering vector x R . Whereas in [15], the authors define α as the DOA. Consequently, the time delay between the two adjacent antennas is d R sin(α )/λ. Clearly, the two definitions are equivalent since sin(α ) = cos(α). Besides, it is worth noting that the latter requires the assumption that the incident direction of the wave is only in the x-y plane but the former does not require any assumptions. Therefore, we adopt the former definition.

The system operation
Before transmission, the user signal x (t) ∈ C N s ×1 , sent by the transmitter, is multiplied by a precoding matrix W (t) ∈ U N T ×N S , where t is the time index and N S is the number of data streams. After passing through fading channels, the signal at the receiver can be formulated as: in which x (t) has independent data streams with its autocorrelation function being [31,32], where the effect of noise on the communication systems can be found in [33,34]. Notation H (t) ∈ C N R ×N T denotes the MIMO flat fading channel, where the fading scenarios can be found in [35][36][37].
Assume that the channel estimation is ideal and the receiver has perfect channel state information (CSI) [38][39][40], H(t). If the channel is quasi-static fading, which means that two adjacent channel matrices can be viewed as identical, i.e., H (t + 1) ≈ H (t), the optimal precoder at time (t + 1) is given by: where U d ∈ U N S ×N S is an arbitrary unitary matrix and V d (t) ∈ U N T ×N S is conducted by selecting the first N S columns of V(t), which comes from singular value decomposition (SVD) of H(t), expressed as: The precoder expressed by (6) can not only maximize the system asymptotic capacity at high signal to noise ratio (SNR), but also minimize the mean square error of data detection. However, if the channel is fast fading, the precoder at current time t is no longer optimal for the next time t + 1. In this case, as in Fig. 1b, we need a predictor, which produces a predicted precoder V d (t + 1) for time t + 1 by the use of V d (t) and V d (t − 1), where the definition of V d (t − 1) is similar to V d (t). Consequently, the precoder at time t + 1 is given by: Clearly, the smaller the distance between V d (t + 1) and V d (t + 1), the better the prediction performance. We use the chordal distance to measure their distance, the details of which will be given in the next section. In addition, in the process of channel prediction, there is a key parameter termed as the step parameter. An aim of this article is to study the effect of various channel settings on the optimal step parameter as well as the expected chordal distance (error).

Preliminaries of the Grassmannian manifold
Grassmannian manifold, denoted as G (N T , N S ) with N S < N T , is important in mathematical differential geometry and has many engineering applications such as channel tracking [41,42], image identification [43], and emotion recognition [44] in recent years. It is the set of all N S -dimensional linear subspaces of an n-dimensional space. The corresponding definition is given by: where [X] is an equivalent class that is a set defined by [X] = XU N S : U N S is unitary . For the sake of convenience, we might as well write X ∈ G (N T , N S ), meaning that X is an equivalent class whose columns span the same p-dimensional subspace [16,45]. For two points X 1 , X 2 ∈ G (N T , N S ), the chordal distance is usually used to measure their distance and it is defined as [14,23]: where · F denotes the Frobenius norm. Besides, the concept of geodesic plays an important role in information geometry. The term geodesic comes from geodesy, the science of measuring the Earth. Initially, a geodesic was the shortest route between two points on the Earth's surface. Now, it is generalized and defined as the shortest route between two points of a manifold. The geodesic emanating from X 1 to X 2 is expressed as [41,42]: where Q 1 ∈ U N T ×N T = X 1 X ⊥ 1 , and X ⊥ 1 ∈ U N T ×(N T −N S ) is some orthogonal basis of the orthogonal complement of X 1 . The matrix B ∈ C N T ×N T is skew Hermitian, and it has the following form: where A ∈ C (N T −N S )×N S is referred to as the velocity matrix. The matrix A can be further written as: where is diagonal with its elements ϕ i , 1 ≤ i ≤ N S . The matrices U 1 ∈ U N S ×N S and U 2 ∈ U (N T −N S )×N S come from the cosine-sine decomposition of a specially constructed matrix [41,42]: in which C is diagonal with its entries cos ϕ i on the diagonal; S is also diagonal with its entries sin ϕ i , on the diagonal, 1 ≤ i ≤ N S ; and V 1 is the right singular matrix of the SVD of X H 1 X 2 , which can also be inferred from Eq. (14). With Eqs. (12)- (14), the key matrix B will be created and the geodesic Eq. 11 can be finally given afterwards. According to the original definition of X(s), we have the following corollary: We believe that the above is important since it can verify the correctness of the geodesic equation; however, it was not proved in [41,42] . Herein, we supply a proof which is shown in the Appendix.

Review of the conventional channel prediction method [41]
The matrix V d (t) can be viewed as a point of Grassmannian manifold. It is reasonable to assume that the matrix V d (t) moves or changes along the geodesic within a short period of time. The conventional prediction obtains the next point mainly by two steps: (1) parallel transport the tangent matrix at V d (t − 1) and (2) construct a new geodesic for prediction. First, given two points V d (t) and V d (t − 1), the tangent matrix at the point V d (t − 1), in the direction of the next point V d (t), is given by [41]: is the velocity matrix, and other details are easy to deduce by referring to the previous subsection.
Then, we transport the tangent matrix (t − 1) ∈ C N T ×N S parallel along the geodesic between V d (t) and V d (t − 1), and obtain a new tangent matrix [45]: where 0 ≤ s 0 ≤ 1 is the step parameter and U F F V H F is the compact SVD of the matrix (t − 1). It is clear that (t, 0) = (t − 1) and (t, 1) is the tangent matrix at the point Finally, the predictive point on the manifold at time t + 1, V Conv d (t, s), can be derived as [23,41,45]: where U E E V H E is the compact SVD of the matrix (t, 1) and s is the step parameter which can be further optimized. The optimal step s OPT can be found by solving the following equation: Once the optimal step is found, we substitute it into (18) to obtain V Conv

The proposed channel prediction method
In this subsection, we propose a novel channel prediction method, in which the next point is predicted by use of the geodesic from V d (t − 1) to V d (t) directly. First, according to the previous subsection A, the geodesic emanating from V d (t −1) to V d (t) is expressed as: 1) and B(t − 1) have similar definitions in the previous subsection A. Clearly, Second, the predictive point on the manifold at time t + 1 is given by: Similarly, the optimal step s OPT in the above formula can be found by solving Eq. (19). Note that the function V The complete prediction process is summarized by Algorithm 1.

Algorithm 1:
The proposed prediction method.

by selecting the columns of the left singular matrix that correspond to nonzero singular values. Note that it is easy to verify the number of nonzero singular values is
and obtain . − 1)). 5: Obtain F (t − 1, s + 1) according to (21); find the optimal step s OPT according to (19) and V Prop d (t + 1) afterwards.
Note that in the calculation of the matrix function exp [B (t − 1) s] , we shall use the following matrix series [46]:

Computational complexity
Usually, the number of floating point flops is used to measure of the complexity of an algorithm. We define a flop as one floating point operation and it has computational complexity O (1). Note that the matrix product X 1 X 2 requires O(mnk) flops, where X 1 ∈ C m×n and X 2 ∈ C n×k ; for an n × n matrix, both its inverse and eigen-decomposition operations require O(n 3 ) flops; for an m × n matrix with m ≥ n, the complexity of its singular value decomposition (SVD) is O(m 2 n) if the complete left singular matrix is needed; otherwise, it is O(mn 2 ) [47] .
To begin with, we analyze the computational complexity of the proposed method. Initially, the two points V d (t −1) and V d (t) are generated from the SVD of H(t −1) and H(t) .
With the above, given s, solving Assuming that the minimal search interval that stands for the accuracy of s OPT is ξ , 1/ξ times of calculating F(t − 1, s) are required to find s OPT . Hence, step 5 has complexity of O N 2 T N S /ξ . In summary, the complexity of the proposed method is O n max n 2 Then, we analyze the complexity of the conventional method briefly. Initially, we also need the two points V d (t − 1) and V d (t). In step 1, to generate V ⊥ d (t − 1), assume that the same method as in the subsection C is used. Further, we note that for an m × n matrix with m ≥ n, the compact SVD has complexity of O mn 2 because the full right singular matrix is not needed. Paying attention to the above points, it is easy to find that the complexity of the conventional method is the same as the proposed method.

Results and discussion
In this section, computer simulation is deployed to investigate the performance of the prediction algorithm. The UAV MIMO channel samples are generated according to the channel model in Section 2. Specifically, a 4 × 4 UAV MIMO system with N S = 2 is adopted. For the transmitter and receiver correlation matrices, θ T and θ R , the well-known exponential correlation model is adopted [48] , in which the (i, j)th entry of the correlation matrix is ρ |i−j| , and the constant ρ is the spatial correlation coefficient. We further denote the correlated coefficients at the transmitter and receiver as ρ T and ρ R , respectively. The angles of arrival and departure α and β are set to 30 • and 20 • , respectively; the antenna spacing d T = d R = λ. Besides, the Doppler shifts for the LOS path and NLOS path are set to be identical, i.e., f d = f d . Observe that given s > 0, the chordal error for the proposed method is always less than that for the conventional method. For case 1 in Fig. 3a, the minimum chordal error for the conventional method is about −15.42 dB at s = 0.9; the minimum chordal error for the proposed one is −18.01 dB at s = 1.0 and it is less than −15.42 dB. Therefore, the proposed method is superior to the conventional method. Similar phenomena can be observed in Fig. 3b. Figure 4 compares two methods with different values of normalized Doppler shift f d T, in which ρ T = 0.2 , ρ R = 0.3, and K Rice = 0 dB. Observe that for both methods, given the step parameter s, increasing f d T results in increasing chordal error. Increasing f d T also has the optimal step parameter s OPT decrease. For instance, for the proposed method, the optimal steps are about 1, 0.9, and 0.6, for f d T =0.01, 0.05, and 0.1, respectively. Similar to Fig. 3, the proposed method has better performance than the conventional method, no matter whether the normalized Doppler shift is large or not.  Figure 5 shows the optimal step parameter s OPT with different normalized Doppler shifts.

The effects of a few key system parameters
The system settings are consistent with those in Fig. 4. Observe that for increasing normalized Doppler shift, s OPT decreases from 1 to 0.6 gradually. The optimal step s OPT is found by solving Eq. (19), where the searching step is set to 0.1, and hence, the curve shows a stepwise downward trend.  Fig. 6 and ρ T = 0.2 for Fig. 5, respectively. From Fig. 6, we observe that given the step parameter s, increasing ρ T will decrease the chordal error and also decrease the optimal s. For instance, given s = 0.5, the chordal errors for ρ T = 0, 0.2, 0.4, 0.6, and 0.8 are −9.2 dB, −9.8 dB, −11.0 dB, −12.8 dB, and −15.8 dB, respectively; the optimal s decreases from 0.9 to 0.8 slightly when increasing ρ T from 0 to 0.8. In Fig. 7, increasing the correlation coefficient also reduces the chordal error, but conversely, the optimal s increases from 0.8 to 0.9 slightly with increasing ρ R . Therefore, these results indicate that both ρ T and ρ R are negatively related to the chordal error; ρ T is slightly negatively related to the optimal step parameter, whereas the opposite is true for ρ R . Figure 8 depicts the effect of the Ricean factor K Rice on the chordal error, where f d T = 0.05, ρ T = 0.2, and ρ R = 0.3. Observe that given s, larger K Rice results in smaller chordal error. This is possibly because that the increase of K Rice reduces the randomness of the UAV channel variation and tracking the LOS path is easier than the NLOS path. In addition, the optimal step s OPT is about 0.9 and nearly unaffected by K Rice .

Conclusions
This article studies the problem of channel predictive precoding in UAV MIMO systems. A novel predictive precoding method that is based on the geodesic in the Grassmannia manifold is proposed. Instead of transporting the tangent vector parallel along the geodesic, the proposed method predicts the next point by use of the geodesic equation directly. We compare the proposed method with the conventional one by simulation and analyze their computational complexity. Besides, the effects of various system parameters, including time and spatial correlated coefficients and the Ricean factor, on the chordal error of the predictor are also thoroughly investigated. The results reveal the following. First, the normalized Doppler shift is positively related to the chordal error, whereas the transmitter and receiver correlation coefficients and the Ricean factor are negatively related to the chordal error. Second, the normalized Doppler shift is negatively related to the optimal step parameter; the transmitter correlation coefficient is slightly negatively related to the optimal step parameter, whereas the opposite is true for the receiver correlation coefficient. Third, with the same computational complexity, the proposed predictive algorithm is superior to the conventional method in terms of the channel predictive error-the channel chordal error. The above research results will provide a reference for possible system implementation in the future. A further study is to investigate the effect of channel estimation errors on the predictor, where the two channel estimates include the previous and current estimates which are not assumed to be ideal. How to design a robust predictor is worthy of researching. Another possible study is to incorporate the codebook quantizer into the predictor design.