Linear precoding based on polynomial expansion: reducing complexity in massive MIMO

Massive multiple-input multiple-output (MIMO) techniques have the potential to bring tremendous improvements in spectral efficiency to future communication systems. Counterintuitively, the practical issues of having uncertain channel knowledge, high propagation losses, and implementing optimal non-linear precoding are solved more or less automatically by enlarging system dimensions. However, the computational precoding complexity grows with the system dimensions. For example, the close-to-optimal and relatively “antenna-efficient” regularized zero-forcing (RZF) precoding is very complicated to implement in practice, since it requires fast inversions of large matrices in every coherence period. Motivated by the high performance of RZF, we propose to replace the matrix inversion and multiplication by a truncated polynomial expansion (TPE), thereby obtaining the new TPE precoding scheme which is more suitable for real-time hardware implementation and significantly reduces the delay to the first transmitted symbol. The degree of the matrix polynomial can be adapted to the available hardware resources and enables smooth transition between simple maximum ratio transmission and more advanced RZF. By deriving new random matrix results, we obtain a deterministic expression for the asymptotic signal-to-interference-and-noise ratio (SINR) achieved by TPE precoding in massive MIMO systems. Furthermore, we provide a closed-form expression for the polynomial coefficients that maximizes this SINR. To maintain a fixed per-user rate loss as compared to RZF, the polynomial degree does not need to scale with the system, but it should be increased with the quality of the channel knowledge and the signal-to-noise ratio.


Introduction
The current wireless networks must be greatly densified to meet the exponential growth in data traffic and number of user terminals (UTs) [1]. The conventional densification approach is to decrease the inter-site distance by adding new base stations (BSs) [2]. However, the cells are subject to more interference from neighboring cells as distances shrink, which requires substantial coordination between neighboring BSs or fractional frequency reuse patterns. Furthermore, serving high-mobility UTs by small cells is Massive multiple-input multiple-output (MIMO) techniques, also known as large-scale multiuser MIMO techniques, have been shown to be viable alternatives and complements to small cells [3][4][5][6][7]. By deploying large-scale arrays with very many antennas at current macro BSs, an exceptional array gain and spatial precoding resolution can be obtained. This is exploited to achieve higher UT rates and serve more UTs simultaneously. In this paper, we consider the single-cell downlink case where one BS with M antennas serves K single-antenna UTs. As a rule of thumb, hundreds of BS antennas may be deployed in the near future to serve several tens of UTs in parallel. If the UTs are selected spatially to have a very small number of common scatterers, the user channels naturally decorrelate as M grows large [8,9] and space-division multiple access (SDMA) techniques become robust to channel uncertainty [3].
One might imagine that by taking M and K large, it becomes terribly difficult to optimize the system throughput. The beauty of massive MIMO is that this is not the case: simple linear precoding is asymptotically optimal in the regime M K 0 [3] and random matrix theory can provide simple deterministic approximations of the stochastic achievable rates [5,[10][11][12][13][14]. These so-called deterministic equivalents are tight as M grows large due to channel hardening but are usually also very accurate at small values of M and K.
Although linear precoding is computationally more efficient than its non-linear alternatives, the complexity of most linear precoding schemes is still intractable in the large (M, K) regime since the number of arithmetic operations is proportional to K 2 M. For example, both the optimal precoding parametrization in [15] and the nearoptimal regularized zero-forcing (RZF) precoding [16] require an inversion of the Gram matrix of the joint channel of all users-this matrix operation has a complexity proportional to K 2 M. A notable exception is the matched filter, also known as maximum ratio transmission (MRT) [17], whose complexity only scales as MK. Unfortunately, this precoding scheme requires roughly an order of magnitude more BS antennas to perform as well as RZF [5]. Since it makes little sense to deploy an advanced massive MIMO system and then cripple the system throughput by using interference-ignoring MRT, treating the precoding complexity problem is the main focus of this paper.
Similar complexity issues appear in multiuser detection, where the minimum mean square error (MMSE) detector involves matrix inversions [18]. This uplink problem has received considerable attention in the last two decades; see [18][19][20][21] and references therein. In particular, different reduced-rank filtering approaches have been proposed, often based on the concept of truncated polynomial expansion (TPE). Simply speaking, the idea is to approximate the matrix inverse by a matrix polynomial with J terms, where J needs not to scale with the system dimensions to maintain a certain approximation accuracy [19]. TPE-based detectors admit simple and efficient multistage/pipelined hardware implementation [18], which stands in contrast to the complicated implementation of matrix inversion. A key requirement to achieve good detection performance at small J is to find good coefficients for the polynomial. This has been a major research challenge because the optimal coefficients are expensive to compute [18]. Alternatives based on appropriate scaling [20] and asymptotic analysis [21] have been proposed. A similar TPE-based approach was used in [22] for the purpose of low-complexity channel estimation in massive MIMO systems.
In this paper, which follows our work in [23], we propose a new family of low-complexity linear precoding schemes for the single-cell multiuser downlink. We exploit TPE to enable a balancing of precoding complexity and system throughput. A main analytic contribution is the derivation of deterministic equivalents for the achievable user rates for any order J of TPE precoding. These expressions are tight when M and K grow large with a fixed ratio but also provide close approximations at small parameter values. The deterministic equivalents allow for optimization of the polynomial coefficients; we derive the coefficients that maximize the throughput. We note that this approach for precoding design is very new. The only other work is [24] by Zarei et al., of which we just became aware at the time this paper was first submitted. Unlike our work, the precoding in [24] is conceived to minimize the sum MSE of all users. Although our approach builds upon the same TPE concept as [24], the design method proposed herein is more efficient since it considers the optimization of the throughput. This metric is usually more pertinent than the sum MSE. Additionally, our work is more comprehensive in that we consider a channel model which takes into account the transmit correlation at the base station.
Our novel TPE precoding scheme enables a smooth transition in performance between MRT (J = 1) and RZF (J = min(M, K)), where the majority of the gap is bridged for small values of J. We show that J is independent of the system dimensions M and K but must increase with the signal-to-noise ratio (SNR) and channel state information (CSI) quality to maintain a fixed per-user rate gap to RZF. We stress that the polynomial structure provides a green radio approach to precoding, since it enables energy-efficient multistage hardware implementation as compared to the complicated/inefficient signal processing required to compute conventional RZF. Also, the delay to the first transmitted symbol is significantly reduced, which is of great interest in systems with very short coherence periods. Furthermore, the hardware complexity can be easily tailored to the deployment scenario or even changed dynamically by increasing and reducing J in highand low-SNR situations, respectively.

Notation
Boldface (lowercase) is used for column vectors, x, and (uppercase) for matrices, X. Let X T , X H , and X * denote the transpose, conjugate transpose, and conjugate of X, respectively, while tr(X) is the matrix trace function. The Frobenius norm is denoted as · , and the spectral norm is denoted as · 2 . A circularly symmetric complex Gaussian random vector x is denoted as x ∼ CN (x, Q), wherē x is the mean and Q is the covariance matrix. The set of all complex numbers is denoted by C, with C N×1 and C N×M being the generalizations to vectors and matrices, respectively. The M × M identity matrix is written as I M , and the zero vector of length M is denoted as 0 M×1 . For an infinitely differentiable monovariate function f (t), the th derivative at t = t 0 (i.e., d / dt f (t)| t=t 0 ) is denoted by f ( ) (t 0 ) and more concisely f ( ) , when t = 0. An analog definition is considered in the bivariate case; in particular, f (l,m) (t 0 , u 0 ) refers to the th and mth derivatives with respect to t and u at t 0 and u 0 , respectively, (i.e.,

System model
This section defines the single-cell system with flat-fading channels, linear precoding, and channel estimation errors.

Transmission model
We consider a single-cell downlink system in which a BS, equipped with M antennas, serves K single-antenna UTs. The received complex baseband signal y k ∈ C at the kth UT is given by where x ∈ C M×1 is the transmit signal and h k ∈ C M×1 represents the random channel vector between the BS and the kth UT. The additive circularly symmetric complex Gaussian noise at the kth UT is denoted by n k ∼ CN (0, σ 2 ) for k = 1, . . . , K, where σ 2 is the receiver noise variance. The small-scale channel fading is modeled as follows.

Assumption 1. The channel vector h k is modeled as
where the channel covariance matrix ∈ C M×M has bounded spectral norm 2 , as M → ∞, and z k ∼ CN (0 M×1 , I M ). The channel vector has a fixed realization for a coherence period and then takes a new independent realization. This model is known as Rayleigh block-fading.
Note that we assume that the UTs reside in a rich scattering environment described by the covariance matrix . This matrix can either be a scaled identity matrix as in [3] or describe array-specific properties (e.g., non-isotropic radiation patterns) and general propagation properties of the coverage area (e.g., for practical sectorized sites). We only consider a common covariance matrix model here, since the main focus in this publication is the precoding scheme. This simplification has been done in many recent publications. Adhikary et al. [25] have proposed to always only serve groups of UTs that share approximately equal covariance matrices, hence providing further motivation behind Assumption 1.
The application of TPE precoding to multicell systems can be found in our paper [26]. However, the models used in this paper and in [26] are incompatible and differ most prominently in the assumption whether the total transmit power increases with the number of users as in [26] or is fixed as in this paper; see (8). This seemingly negligible change has a big impact on the analysis and applicability of the models, as this assumption means that the noise term in [26] becomes asymptotically zero, while in the current work, the noise term is non-negligible. The channel estimation model in [26] and in this paper is also different, and the calculations follow very different approaches, due to the inclusion of power control later on. Another big extension in the current work is the complete complexity analysis of the TPE approach in comparison to the classical RZF approach. Only this analysis gives TPE precoding its motivation and pertinence. Finally, we want to point out that the optimization in [26] is with respect to a max-min SNR problem and the solution is not given as a closed form, while here we maximize the throughput and find a closed-form solution. Before utilizing our work, one needs to decide which model gives the most accurate asymptotic behavior for the specific type of system considered.

Assumption 2.
The BS employs Gaussian codebooks and linear precoding, where g k ∈ C M×1 denotes the precoding vector and s k ∼ CN (0, 1) is the data symbol of the kth UT.
Based on this assumption, the transmit signal in (1) is g n s n = Gs.
( 3 ) The matrix notation is obtained by letting G = [ g 1 . . . g K ] ∈ C M×K be the precoding matrix and s = [ s 1 . . . s K ] T ∼ CN (0 K×1 , I K ) be the vector containing all UT data symbols.
Consequently, the received signal (1) can be expressed as Let G k ∈ C M×(K−1) be the matrix G with column g k removed. Then, the SINR at the kth UT becomes By assuming that each UT has perfect instantaneous CSI, the achievable data rates at the UTs are r k = log 2 (1 + SINR k ), k = 1, . . . , K.

Model of imperfect channel information at transmitter
Since we typically have M ≥ K in practice, we assume that we either have a time-division duplex (TDD) protocol where the BS acquires channel knowledge from uplink pilot signaling [5] or a frequency-division duplex (FDD) protocol where temporal correlation is exploited as in [27]. In both cases, the transmitter generally has imperfect knowledge of the instantaneous channel realizations and we model this by the generic Gauss-Markov formulation; see [12,28,29]: The transmitter has an imperfect channel estimate models the independent error. The scalar parameter τ ∈ [ 0, 1] indicates the quality of the instantaneous CSI, where τ = 0 corresponds to perfect instantaneous CSI and τ = 1 corresponds to having only statistical channel knowledge.
The parameter τ depends on factors such as time/power spent on pilot-based channel estimation and user mobility. Note that we assume for simplicity that the BS has the same quality of channel knowledge for all UTs. Based on the model in (6), the matrix denotes the joint imperfect knowledge of all user channels.

Linear precoding
Many heuristic linear precoding schemes have been proposed in the literature, mainly because finding the optimal precoding (in terms of weighted sum rate or other criteria) is very computationally demanding and thus unsuitable for fading systems [30]. Among the heuristic schemes, we distinguish RZF precoding [16], which is also known as transmit Wiener filter [31], signal-toleakage-and-noise ratio maximizing beamforming [32], generalized eigenvalue-based beamformer [33], and virtual SINR maximizing beamforming [34]. The reason that RZF precoding has been proposed by different authors (under different names) is, most likely, that it provides close-to-optimal performance in many scenarios. It also outperforms classical MRT and zero-forcing beamforming (ZFBF) by combining the respective benefits of these schemes [30]. Therefore, RZF is deemed the natural starting point for this paper. Next, we provide a brief review of RZF and prior performance results in massive MIMO systems. These results serve as a starting point for Section 3.2, where we propose an alternative precoding scheme with a computational/hardware complexity more suited for large systems.

Review on RZF precoding in massive MIMO systems
Suppose we have a total transmit power constraint We stress that the total power P is fixed, while we let the number of antennas, M, and number of UTs, K, grow large.
Similar to [12], we define the RZF precoding matrix as where the power normalization parameter β is set such that G RZF satisfies the power constraint in (8) and P is a fixed diagonal matrix whose diagonal elements are power allocation weights for each user. We assume that P satisfies the following: The scalar regularization coefficient ξ can be selected in different ways, depending on the noise variance, channel uncertainty at the transmitter, and system dimensions [12,16]. In [12], the performance of each UT under RZF precoding is studied in the large (M, K) regime. This means that M and K tend to infinity at the same speed, which can be formalized as follows.

Assumption 5.
In the large (M, K) regime, M and K tend to infinity such that The user performance is characterized by SINR k in (5). Although the SINR is a random quantity that depends on the instantaneous values of the random users channels in H and the instantaneous estimate H, it can be approximated using deterministic quantities in the large (M, K) regime [10][11][12][13]. These are quantities that only depend on the statistics of the channels and are often referred to as deterministic equivalents, since they are almost surely (a.s.) tight in the asymptotic limit. This channel hardening property is essentially due to the law of large numbers. Deterministic equivalents were first proposed by Hachem et al. in [10], who have also shown their ability to capture important system performance indicators. When the deterministic equivalents are applied at finite M and K, they are referred to as large-scale approximations.
In the sequel, by deterministic equivalent of a sequence of random variables X n , we mean a deterministic sequence X n which approximates X n such that As an example, we recall the following result from [10], which provides some widely known results on deterministic equivalents. Note that we have chosen to work with a slightly different definition of the deterministic equivalents than in [10], since this better fits the analysis of our proposed precoding scheme.
and let U be any matrix with bounded spectral norm. Under Assumption 5 and for t > 0, we have The statement in (11) shows that 1 K tr(UT(t)) is a deterministic equivalent to the random quantity 1 K tr(UQ(t)). In this paper, the deterministic equivalents are essential to determine the limit to which the SINRs tend in the large (M, K) regime. For RZF precoding, as in (9), this limit is given by the following theorem.

Theorem 2.
(Adapted from Corollary 1 in [12]) Let ρ = P σ 2 and consider the notation T = T( 1 ξ ) and δ = δ( 1 ξ ). Define the deterministic scalar quantities Then, the SINRs with RZF precoding satisfies Note that all UTs obtain the same asymptotic value of the SINR since the UTs have homogeneous channel statistics. Theorem 2 holds for any regularization coefficient ξ , but the parameter can also be selected to maximize the limiting value θ of the SINRs. This is achieved by the following theorem.

Theorem 3. (Adapted from Proposition 2 in [12])
Under the assumption of a uniform power allocation, p k = P K , the large-scale approximated SINR in (12) under RZF precoding is maximized by the regularization parameter ξ , given as the positive solution to the fixed-point equation where ν(ξ ) is given by The RZF precoding matrix in (9) is a function of the instantaneous CSI at the transmitter. Although the SINRs converges to the deterministic equivalents given in Theorem 2, in the large (M, K) regime, the precoding matrix remains a random quantity that is typically recalculated on a millisecond basis (i.e., at the same pace as the channel knowledge is updated). This is a major practical issue, because the matrix inversion operation in RZF precoding is very computationally demanding in large systems [35]; the number of operation scale as O(K 2 M) and the known inversion algorithms are complicated to implement in hardware (see Section 4 for details). The matrix inversion is the key to interference suppression in RZF precoding, thus there is a need to develop less complicated precoding schemes that still can suppress interference efficiently.

Truncated polynomial expansion precoding
Motivated by the inherent complexity issues of RZF precoding, we now develop a new linear precoding class that is much easier to implement in large systems. The precoding is based on rewriting the matrix inversion by a polynomial expansion, which is then truncated. The following lemma provides a major motivation behind the use of polynomial expansions.

Lemma 1. For any positive definite Hermitian matrix X,
where the second equality holds if the parameter κ is selected such that 0 < κ < 2 max n λ n (X) . Proof 1. The inverse of an Hermitian matrix can be computed by inverting each eigenvalue, while keeping the eigenvectors fixed. This lemma follows by applying the standard Taylor series expansion x , for any |x| < 1, on each eigenvalue of the Hermitian matrix (I−κX). The condition on x corresponds to requiring that the spectral norm I − κX 2 is bounded by unity, which holds for κ < 2 max n λ n (X) . See [20] for an in-depth analysis of such properties of polynomial expansions.
This lemma shows that the inverse of any Hermitian matrix can be expressed as a matrix polynomial. More importantly, the low-order terms are the most influential ones, since the eigenvalues of (I − κX) converge geometrically to zero as grows large. This is due to each eigenvalue λ of (I − κX) having an absolute value smaller than unity, |λ| < 1, and thus λ goes geometrically to zero as → ∞. As such, it makes sense to consider a TPE of the matrix inverse using only the first J terms. This corresponds to approximating the inversion of each eigenvalue by a Taylor polynomial with J terms, hence the approximation accuracy per matrix element is independent of M and K; that is, J needs not change with the system dimensions.
TPE has been successfully applied for low-complexity multiuser detection in [18][19][20][21] and channel estimation in [22]. Next, we exploit the TPE technique to approximate RZF precoding by a matrix polynomial. Starting from G RZF in (9), we note that where (15) follows directly from Lemma 1 (for an appropriate selection of κ), (16) is achieved by truncating the polynomial (only keeping the first J terms), and (17) follows from applying the binomial theorem and gathering the terms for each exponent. Inspecting (17), we have a precoding matrix with the structure where w 0 , . . . , w J−1 are scalar coefficients. Although the bracketed term in (17) provides a potential expression for w , we stress that these are generally not the optimal coefficients when J < ∞. Also, these coefficients are not satisfying the power constraint in (8) since the coefficients are not adapted to the truncation. Hence, we treat w 0 , . . . , w J−1 as design parameters that should be selected to maximize the performance; for example, by maximizing the limiting value of the SINRs, as was done in Theorem 3 for RZF precoding. We note especially that the value of κ in (17) does not need to be explicitly known in order to choose, optimize, and implement the coefficients. We only need for κ to exist, which is always the case under Assumption 2. Besides the simplified structure, the proposed precoding matrix G TPE possesses a higher number of degrees of freedom (represented by the J scalars w ) than the RZF precoding (which has only the regularization coefficient ξ ). The precoding in (18) is coined TPE precoding and actually defines a whole class of precoding matrices for different J. For J = 1, we obtain G = w 0 √ K HP 1 2 , which equals MRT. Furthermore, RZF precoding can be obtained by choosing J = min(M, K) and coefficients based on the characteristic polynomial of ( 1 K H H H + ξ I M ) −1 (directly from Cayley-Hamilton theorem). We refer to J as the TPE order and note that the corresponding polynomial degree is J − 1. Clearly, proper selection of J enables a smooth transition between the traditional low-complexity MRT and the high-complexity RZF precoding. Based on the discussion that followed Lemma 1, we assume that the parameter J is a finite constant that does not grow with M and K.

Complexity analysis
In this section, we compare the complexities of RZF and TPE precoding in a theoretical fashion and in an implementation sense. The complexities are given as simple numbers of complex addition and multiplication operations needed for a given arithmetic operation. The number of floating point operations (flops) needed to implement these complex operations varies greatly according to the used hardware and complex number representation (i.e., polar or Cartesian). Thus, we will not attempt to give a measure in flops. Also, the ability to parallelize operations and to customize algorithm-specific circuits has a fundamental impact on the computational delays and energy consumption in practical systems.

Sum complexity per coherence period for RZF and TPE
In order to compare the number of complex operations needed for conventional RZF precoding and the proposed TPE precoding, it is important to consider how often each operation is repeated. There are two time scales: (1) operations that take place once per coherence period (i.e., once per channel realization) and (2) operations that take place every time the channel is used for downlink transmission. To differentiate between these time scales, we let T pcp data denote the number of downlink channel uses for data transmission per coherence period. Recall from (3) that the transmit signal is Gs, where the precoding matrix G ∈ C M×K changes once per coherence period and the data transmit symbols s ∈ C K×1 are different for each channel use.
The RZF precoding matrix in (9) is computed once per coherence period. There are two equivalent expressions in (9), where the difference is that the matrix inversion is either of dimension K × K or M × M. Since K ≤ M in most cases of practical interest, and especially in the massive MIMO regime, we consider the first precoding expression: 1 Assuming that 1 √ K H, ξ , β, and P 1 2 are available in advance and the Hermitian operation is "free, " we need to (1) compute the matrix-matrix multiplication (4) multiply the result with the diagonal matrix resulting from P 1 2 β. These are standard operations for matrices; thus, we obtain the numbers of complex operations as: K 2 (2M − 1), K, K 3 3 + 2K 2 M, and MK + K operations, respectively.
Step 3 is not immediately obvious, but an efficient method for this part is to compute a Cholesky factorization of 1 K H H H + ξ I K (at a cost of K 3 /3) and then solve a simple linear equation system for each row of 9). This approach is preferable to the alternative of completely inverting the matrix (again using Cholesky factorization) and then using matrix-matrix multiplication, as long as K 3 − KM > 0. Given that the alternative method has a cost of 4K 3 /3 + MK(2K − 1). It is interesting to note here that, for the case of M K, the matrix-matrix multiplication is actually more expensive than the matrix inversion (2MK 2 vs. K 3 ). 1 Once G RZF has been computed, the matrix-vector multiplication G RZF s requires M(2K − 1) operations per channel use of data transmission. In summary, RZF precoding has a total number of complex operations per coherence period of There is a second approach to looking at the RZF precoder complexity. Let the transmit signal with RZF precoding at channel use t be denoted as . Thus, one can replace the "matrix times inverse of another matrix" operation taking place each coherence period by a matrix-inverse operation per coherence period and two matrix-vector multiplications per data symbol vector. Thus, one effectively splits the previous point (3) into two parts and waits for the symbol vector to allow for the matrix-vector multiplications. This results in . Still, this complexity is dominated by the matrix-matrix multiplication inside the inverse. However, the per coherence period complexity is reduced in exchange for a slight increase in complexity per symbol. Depending on the usecase of the precoder, this change can either be advantageous or disadvantageous (see Fig. 1 and Subsection 4.2). We note that choosing to incorporate the multiplication with P 1 2 per coherence period or per symbol vector does only insignificantly change the stated outcomes. In the following, we will chose the appropriate version for each comparison. Next, we consider TPE precoding. Similar to before, we assume that 1 √ K H, w , and P 1 2 are available in advance and the Hermitian operation is "free. " Let the transmit signal vector with TPE precoding at channel use t be denoted as x (t) TPE and observe that it can be expressed as where s (t) is the vector of data symbols at channel use t andx This reveals that there is an iterative way of computing the J terms in TPE precoding. The benefit of this approach is that it can be implemented using only matrix-vector multiplications. 2 Similar to the above, we conclude that the case = 0 uses K + M(2K − 1) operations and each of the J − 1 cases of ≥ 1 needs M(2K − 1) + K(2M − 1) operations. One remarks that it is impractical and unneeded to carry out a matrix-matrix multiplication at this step. Finally, the multiplication with w and the summation require M(2J − 1) further operations. In summary, TPE precoding has a total number of arithmetic operations of When comparing RZF and TPE precoding, we note that the complexity of precomputing the RZF precoding matrix is very large, but it is only done once per coherence period. The corresponding matrix G TPE for TPE precoding is never computed separately but only indirectly as G TPE s for each data symbol vector s. Intuitively, precomputation is beneficial when the coherence period is long (compared to M and K), and the sequential computation of TPE precoding is beneficial when the system dimensions M and K are large (compared to the coherence period) or the coherence period is short. This is seen from the large dimensional complexity scaling which is O (4K 2 M) One should not forget the overhead signaling required to obtain CSI at the UTs, which makes the number of channel uses T data available for data symbols reduce with K. For example, suppose T coherence is the total coherence period and that we use a TDD protocol, where η DL is the fraction used for downlink transmission and μK channel uses (for some μ ≥ 1) are consumed by downlink pilot signals that provide the UTs with sufficient CSI. We then have T data = η DL T coherence − μK. Using this relationship, the number of arithmetic operations are illustrated numerically in Fig. 1 for η DL = 1 2 , K = 100, and μ = 2. 3 This figure shows that TPE precoding uses fewer operations than RZF precoding when the coherence period is short and the TPE order is small, while RZF is competitive for long coherence times.
We remark that all previously found results change in favor of TPE, if one uses the canonical transformation of complex to real operations by doubling all dimensions.

Remark 1. Power normalization.
In this section, we assumed that β and w (and ξ ) are known beforehand. These factors are responsible for the power normalization of the transmit signal. Depending on the chosen normalization, for example, the average per one UT in this paper requires the full precoding matrix to be known. Thus, it forbids the alternative implementation of RZF precoding detailed before. Note that this could be remedied by changing to "strict" per UT normalization. In general, we can find values for β and w , which only rely on channel statistics and are valid in the large (M, K) regime. This, and the possible fix for the alternative RZF approach, has motivated us to assume β and w as known.

Delay to the first transmission for RZF and TPE
A practically important complexity metric is the number of complex operations for the first channel use. This number can also be interpreted as the delay until the start of data transmission. This complexity can easily be found from the previous results, by choosing T data = 1. Directly looking at the massive MIMO case, we find C 1st RZF = 4MK 2 , C 1st RZF2 = 2MK 2 , and C 1st TPE = 4JMK. Hence, the first data vector is transmitted by a factor of K/(2J) earlier, 4 when TPE precoding is employed. This factor is significant and gives TPE precoding practical relevance, especially in massive MIMO systems and in very fast changing environments, i.e., when coherence periods are very short. We also remark that not wasting time during the coherence period pays off greatly, as the lost channel uses are given by the saved time multiplied by the (often large) coherence bandwidth.

Implementation complexity of RZF and TPE precoding
In practice, the number of arithmetic operations is not the main issue, but the implementation cost in terms of hardware complexity, time delays, and energy consumption. The analysis in Subsection 4.1 showed that we can only expect improvements in the sum of complex operations from TPE precoding per coherence period in certain scenarios. However, one advantage of TPE precoding is that it enables multistage hardware implementation where the computations are pipelined [20] over multiple processing cores (e.g., application-specific integrated circuits (ASICs)). This structure is illustrated in Fig. 2, where the transmitted signal x (t) is prepared in various cores (black path), while the preceding and succeeding transmit signals are computed in the "free" cores (gray paths). Each processing core performs two simple matrix-vector multiplications, each requiring approximately O(2MK) complex additions and multiplications per coherence period. This is relatively easy to implement using ASICs or FPGAs, which are know to be very energy-efficient and have low production cost. Consequently, we can select the TPE order J as large as needed to obtain a certain precoding accuracy, if we are prepared to use as many circuits of the same type as needed. Then, the delay between two consecutive transmitted symbol vectors is given only by the delay of two matrix-vector multiplications.
In comparison, the inversion of RZF precoding can only be pseudo-parallelized by using tree structures. Hence, the pipelining of the C RZF complex operations per coherence period is limited by the delay of a single processing core that implements the inverse of a matrix-matrix; this delay is most probably much larger than the two matrix-vector multiplications of TPE. The delay of a second core implementing the multiplication of the inverse Fig. 2 Simple pipelined implementation. A simple pipelined implementation of the proposed TPE precoding with J = 2, which removes the delays caused by precomputing the precoding matrix. Each block performs a simple matrix-vector multiplication, which enables highly efficient hardware implementation, and J can be increased by simply adding additional cores with the channel matrix is negligible in comparison. Like mentioned before, the precomputation of the RZF precoding matrix causes non-negligible delays that forces T pcp data to be smaller than for TPE precoding; for example, [35] describes a hardware implementation from [37] where it takes 0.15 ms to compute RZF precoding for K = 15, which translated to a loss of 0.15 ms × 200 kHz = 30 channel uses in a system with coherence bandwidth 200 kHz. Also, the number of active UTs can be much larger than this in large-scale MIMO systems [38]. TPE precoding does not cause such delays because there are no precomputations-the arithmetic operations are spread over the coherence period.
In practice, this means one can argue that only the curve pertaining to J = 1 in Fig. 1 is relevant for comparisons between TPE and RZF after implementation; if one is prepared to add (seemingly unfairly) as many computation cores as necessary to TPE.

Analysis and optimization of TPE precoding
In this section, we consider the large (M, K) regime, defined in Assumption 5. We show that SINR k , for k = 1, . . . , K, under TPE precoding converges to a limit, a deterministic equivalent, that depends only on the coefficients w , the respective attributed power p k , and the channel statistics.
Recall the SINR expression in (5) and observe that g k = where e k is the kth column of the identity matrix I K . By substituting the TPE precoding expression (18) into (5), it is easy to show that the SINR writes as where w =[ w 0 . . . w J−1 ] T and the ( , m)th elements of the for = 0, . . . , J − 1 and m = 0, . . . , J − 1. 5 Since the random matrices A k and B k are of finite dimensions, it suffices to determine a deterministic equivalent for each of their elements. To achieve this, we express them using the resolvent matrix of H. This can be done by introducing the following random functionals in t and u: By taking derivatives of X k,M (t, u) and Z k,M (t, u) at the point (t, u) = (0, 0), we obtain Substituting (24) and (25) into (20) and (21), respectively, we obtain the alternative expressions It, thus, suffices to study the asymptotic convergence of the bivariate functions X k,M (t, u) and Z k,M (t, u). This is achieved by the following new theorem and its corollary: Theorem 4. Consider a channel matrix H whose columns are distributed according to Assumption 3. Under the asymptotic regime described in Assumption 5, we have and β M (t, u) is given by

(26)
Proof 2. The proof leans heavily on lemmas presented in Appendix 1 and is detailed in Appendix 2. Corollary 1 shows that the entries of A k and B k , which depend on the derivatives of X k,M (t, u) and Z k,M (t, u), can be approximated in the asymptotic regime by T ( ) and δ ( ) , which are the derivatives of T(t) and δ(t) at t = 0. Such derivatives can be computed numerically using the iterative algorithm of [21], which is provided in Appendix 6 for the sake of completeness.
It remains to compute the aforementioned derivatives. To this end, we denote f (t) = − 1 1+tδ(t) , T (t) = −f (t)T(t), and by f ( ) , T ( ) their derivatives at t = 0. T ( ) can be calculated using the Leibniz derivation rule the respective values from Appendix 6. Rewriting (26) as and using the Leibniz rule, we obtain for any integers and m greater than 1, the expression An iterative algorithm for the computation of β ( ,m) M is given in Appendix 5.
With these derivation results on hand, we are now in the position to determine the expressions for the derivatives of the quantities of interest, namely X k,m (t, u) and b M (t, u). Using again the Leibniz derivation rule, we obtain Using these results in combination with Corollary 1, we immediately obtain the asymptotic equivalents of A k and B k : Corollary 2. Let A and B be the J × J matrices, whose entries are Then, in the asymptotic regime, for any k ∈ 1, . . . , K we have

Optimization of the polynomial coefficients
Next, we consider the optimization of the asymptotic SINRs with respect to the polynomial coefficients w = [ w 0 . . . w J−1 ] T . Using results from the previous sections, a deterministic equivalent for the SINR of the kth UT is The optimized TPE precoding should satisfy the power constraints in (8) tr G TPE G H TPE = P (27) or equivalently where the ( , m)th element of the J × J matrix C is In order to make the optimization problem independent of the channel realizations, we replace the constraint in (28) by a deterministic one, which depends only on the statistics of the channel. To find a deterministic equivalent of the matrix C, we introduce the random quantity Using the same method as for the matrices A and B, we achieve the following result: Theorem 5. Considering the setting of Theorem 4, we have the following convergence results: 3. Let C be the J × J matrix with entries given by Then, in the asymptotic regime Proof 4. The proof relies on the same techniques as before, so we provide only a sketch in Appendix 7.
Based on Theorem 5, we can consider the deterministic power constraint tr (P) w T Cw = P (30) which can be seen as an approximation of (28), in the sense that for any w satisfying (30), we have Now, the maximization of the asymptotic SINR of UT k amounts to solving the following optimization problem: The next theorem shows that the optimal solution, w opt , to (31) admits a closed-form expression. Theorem 6. Let a be a unit norm eigenvector corresponding to the maximum eigenvalue λ max of Then, the optimal value of the problem in (31) is achieved by where the scaling factor α is Moreover, for the optimal coefficients, the asymptotic SINR for the kth UT is Proof 5. The proof is given in Appendix 8.
The optimal polynomial coefficients for UT k are given in (33) of Theorem 6. Interestingly, these coefficients are independent of the user index; thus, we have indeed derived the jointly optimal coefficients. Furthermore, all users converge to the same deterministic SINR up to an UT-specific scaling factor Kp k γ tr(P) .

Remark 2.
The asymptotic SINR expressions in (35) are only functions of the statistics and the power allocation p 1 , . . . , p K . The power allocation can be optimized with respect to some system performance metric. For example, one can show that the asymptotic average achievable rate is maximized by a uniform power allocation p k = P K for all k.
Remark 3. Theorem 6 shows that the J polynomial coefficients that jointly maximize the asymptotic SINRs can be computed using only the channel statistics and the channel estimation error. The optimal coefficients are then given in closed form in (33). Numerical experiments show that the coefficients are very robust to underestimation of τ and robust to overestimation. Hence, the main feature of Theorem 6 is that the TPE precoding coefficients can be computed beforehand or at least be updated at the relatively slow rate of change of the channel statistics. Thus, the cost of the optimization step is negligible with respect to calculating the precoding itself. The performance of finite-dimensional large-scale MIMO systems is evaluated numerically in Section 6.
Remark 4. Finally, we remark that Assumption 5 prevents us from directly analyzing the scenario where K is fixed and M → ∞, but we can infer the behavior of TPE precoding based on previous works. In particular, it is known that MRT is an asymptotically optimal precoding scheme in this scenario [4]. We recall from Section 3.2 that TPE precoding reduces to MRT for J = 1. Hence, we expect the optimal coefficients to behave as w 0 = 0 and w → 0 for ≥ 1 when M → ∞. In other words, we can reduce J as M grows large and still keep a fixed performance gap to RZF precoding.

Simulation results
In this section, we compare the RZF precoding from [16] (which was restated in (9)) with the proposed TPE precoding (defined in (18)) by means of simulations. The purpose is to validate the performance of the proposed precoding scheme and illustrate some of its main properties. The performance measure is the average achievable rate of the UTs, where the expectation is taken with respect to different channel realizations and users. In the simulations, we model the channel covariance matrix as where a is chosen to be 0.1. This approach is known as the exponential correlation model [39]. More involved models could be chosen here but would make it harder to evaluate the performance and function of TPE, while not offering more insight. The sum power constraint tr G RZF/TPE G H RZF/TPE = P is applied for both precoding schemes. Unless otherwise stated, we use uniform power allocation for the UTs, since the asymptotic properties of RZF precoding are known in this case (see Theorem 3). Without loss of generality, we have set σ 2 = 1. Our default simulation model is a largescale single-cell MIMO system of dimensions M = 128 and K = 32. We first take a look at Fig. 3. It considers a TPE order of J = 3 and three different quality levels of the CSI at the BS: τ ∈ {0.1, 0.4, 0.7}. From Fig. 3, we see that RZF and TPE achieve almost the same average UT performance when a bad channel estimate is available (τ = 0.7). Furthermore, TPE and RZF perform almost identically at low SNR values, for any τ . In general, the unsurprising observation is that the rate difference becomes larger at high SNRs and when τ is small (i.e., with more accurate channel knowledge). Figure 4 shows more directly the relationship between the average achievable UT rates and the TPE order J. We consider the case τ = 0.1, M = 512, and K = 128, in order to be in a regime where TPE performs relatively bad (see Fig. 3) and the precoding complexity becomes an issue. From the figure, we see that choosing a larger value for J gives a TPE performance closer to that of RZF. However, doing so will also require more hardware; see Section 4.3. The proposed TPE precoding never surpasses the RZF performance, which is noteworthy since TPE has J degrees of freedom that can be optimized (see Section 5.1), while RZF only has one design parameter. Hence one can regard RZF precoding as an upper bound to TPE precoding in the single-cell scenario. 6 It is desirable to select the TPE order J in such a way that we achieve a certain limited rate-loss with respect to RZF precoding. Figure 5 illustrates the rate-loss (per UT) between TPE and RZF, while the number of UTs K and transmit antennas M increase with a fixed ratio (M/K = 4). The figure considers the case of τ = 0.1. We observe that the TPE order J and the system dimensions are independent in their respective effects on the rate-loss between TPE and RZF precoding. This observation is in line with previous results on polynomial expansions, for example, [19] where reduced-rank received filtering was considered. The independence between J and the system dimensions M and K (given the same ratio) is indeed a main motivation behind TPE precoding, because it implies that the order J can be kept small even when TPE precoding is applied to very large-scale MIMO systems. The intuition behind this result is that the polynomial expansion approximates the inversion of each eigenvalue with the same accuracy, irrespective of the number of eigenvalues; see Section 3.2 for details. Although the relative performance loss is unaffected by the system dimensions, we also see that J needs to be increased along with the SNR, if a constant performance gap is desired. In the simulation depicted in Fig. 6, we introduce a hypothetical case of TPE precoding (TPEopt) that optimizes the J coefficients using the estimated channel coefficients in each coherence period, instead of relying solely on the channel statistics. More precisely, the optimal coefficients in Theorem 6 are not computed using the deterministic equivalents of A, B, and C but using the original matrices from (20), (21), and (29). This plot illustrates the additional performance loss caused by precalculating the TPE coefficients based on channel statistics and asymptotic analysis, instead of carrying out the optimization step for each channel realization. The difference is virtually zero at low SNRs and high at high SNRs. Furthermore, we note that increasing the value of J has the same performance-gap-reducing effect on TPEopt, as it has on TPE (see Figs. 4 and 5). In order to preserve readability, only the curves pertaining to J = 3 are shown in Fig. 6.
Finally, to assess the validity of our results, we treat the case of non-uniform power allocation (i.e., with different values for p k ). In particular, we considered a situation where the users are divided into four classes corresponding to {c 1 , c 2 , c 4 c 4 } = {1, 2, 3, 4}, where p k = c k K in order to adhere to the scaling in Assumption 4. Figure 7 shows the theoretical large (M, K) regime (DE; based on (35)) and empirical (MC; based on (19)) average rate per UT for each class, when K = 32, M = 128, and τ = 0.1. We especially remark the very good agreement between our theoretical analysis and the empirical system performance.

Conclusions
Conventional RZF precoding provides attractive system throughput in massive MIMO systems, but its computational and implementation complexity is prohibitively high, due to the required channel matrix inversion. In this paper, we have proposed a new class of TPE precoding schemes where the inversion is approximated by truncated polynomial expansions to enable simple hardware implementation. In the single-cell downlink with M transmit antennas and K single-antenna users, this new class can approximate RZF precoding to an arbitrary accuracy by choosing the TPE order J in the interval 1 ≤ J ≤ min(M, K). In terms of implementation complexity, TPE precoding has several advantages: (1) There is no need to compute the precoding matrix beforehand (which leaves more channel uses for data transmission); (2) the delay to the first transmitted symbol is reduced significantly; (3) the multistage structure enables pipelining; and (4) the parameter J can be tailored to the available hardware.
Although the polynomial coefficients depend on the instantaneous channel realizations, we have shown that the per-user SINRs converge to deterministic values in the large (M, K) regime. This enabled us to compute asymptotically optimal coefficients using merely the statistics of the channels. The simulations revealed that the difference in performance between RZF and TPE is small at low SNRs and for large CSI errors. The TPE order J can be chosen very small in these situations, and in general, it does not need to scale with the system dimensions. However, to maintain a fixed per-user rate-loss compared to RZF, J should increase with the SNR or as the CSI quality improves.

Endnotes
1 Matrix multiplication combined with matrix inversion can be implemented using the Strassen's algorithm in [40] and the improved Coppersmith-Winograd algorithm in [41]. These are divide-and-conquer algorithms that exploit that 2 × 2 matrices can be multiplied efficiently and thereby reduce the asymptotic complexity of multipling/inverting K × K matrices to O(K 2.8074 ) and O(K 2.373 ), respectively. Unfortunately, the overhead in these algorithms is heavy and thus K needs to be at the order of several thousands to achieve a lower complexity than the Cholesky approach considered here. Hence, these alternative algorithms are unfavorable for matrices of practical sizes. 2 Intuitively, one circumvents the expensive matrix-matrix multiplication with a domino-like chain of 2J − 1 (less expensive) matrix-vector multiplications per transmitted symbol vector. This became possible by replacing the inverse of a matrix-matrix multiplication in the RZF with a sum of weighted matrix powers. 3 These parameter values correspond to symmetric downlink/uplink transmission, 2 downlink pilot symbols per UT (at different frequencies). Looking at values similar the LTE standard ( [42] Chapter 10), e.g., a coherence bandwidth of 200 kHz, and a coherence period of 5 ms one would arrive a T coherence of 1000. 4 Depending on the massive MIMO system, K can be on the order of 100 and M of the order 10K, while we will see later that J = 4 is sufficient for many cases. 5 The entries of matrices are numbered from 0, for notational convenience. 6 The optimal precoding parametrization in [15] has K − 1 parameters. To optimize some general performance metric, it is therefore necessary to let the number of design parameters scale with the system dimensions.
Using the Cauchy-Schwartz inequality, we see that which establishes the desired result.

Appendix 2: Proof of Theorem 4
Here, we present the proof of Theorem 4, which establishes the asymptotic convergence of X k,M (t, u) and Z k,M (t, u) to deterministic quantities.

Deterministic equivalent for X k,M (t, u)
We will begin the proof by looking at the random quantity X k,M (t, u). Using the notation of Lemma 2, we can write To control the quadratic form 1 K h H k Q(t) h k , we need to remove the dependency of Q(t) on vector h k . For that, we shall use the relation in (36), thereby yielding Using Lemma 3, we thus have 0, by the rank-one perturbation property in Lemma 5, we have Finally, Theorem 1 implies that The same kind of calculations can be used to deal with the quadratic form 1 K h H k Q k (t) h k , whose asymptotic limit is the same as due to the independence between the channel estimation error and the channel vector h k . Hence, Plugging the deterministic approximation of (39) and (40) into (38), we thus see that

Deterministic equivalent for Z k,M (t, u)
Finding a deterministic equivalent for Z k,M (t, u) is much more involved than for X k,M (t, u). Following the same steps as in section "Deterministic equivalent for X k,M (t, u)" Appendix 2, we decompose Z k,M (t, u) as As it will be shown next, to determine the asymptotic limit of the random variables X i (t, u), i = 1, . . . , 4, we need to find a deterministic equivalent for This is the most involved step of the proof. It will, thus, be treated separately in Appendix 3, where we establish the following lemma: Then, in the asymptotic regime described by Assumption 5, we have .

The proof of this lemma is adjourned to Appendix 3.
Let us begin by treating X 1 (t, u) The right-hand side term in the equation above can be treated using (40), thereby yielding Using Lemma 3, we can prove that Continuing, according to Lemma 7, we have Combining (42) with (43) yields Thus, in the asymptotic regime, we have Controlling the other terms X i (t, u), i = 2, 3, 4, will also include the term β(t, u). First note that X 2 (t, u) is given by Observe that Y 2 (t, u) is very similar to X 1 (t, u). The only difference is that Y 2 (t, u) is a quadratic form involving vectors h k and h k , whereas X 1 (t, u) involves only the vector h k . Following the same kind of calculations leads to we now have Similarly, X 3 (t, u) satisfies Finally, X 4 (t, u) can be treated using the same approach, thereby providing the following convergence:

Appendix 3: Proof of Lemma 7
The aim of this section is to determine a deterministic equivalent for the random quantity α M (t, u, A) = 1 K tr AQ(t)HPH H Q(u) .
The proof is technical and will make frequent use of results from Appendix 1. First, we need to control var ( α M (t, u)). This has already been treated in [10] where it was proved that var ( α M (t, u, A)) = O(K −2 ) when t = u. The same calculations hold for t = u, thus we consider in the sequel that var ( α M (t, u, A) We will only directly deal with the terms Z 1 and Z 3 , since Z 2 will be compensated by terms in Z 3 . We begin with Z 1 : .
Using Lemma 3, we can show that the first term on the right-hand side of the above equation is negligible. Therefore, Using Lemma 5, we have We now look at Z 3 , where

tr AT(t)h h H Q(t)HPH H Q(u) .
Using (37), we arrive at tr AT(t)h h H Q (t)HPH H Q(u) From (36), Z 3 can be decomposed as tr AT(t)h h H Q (t)HPH H Q (u) We sequentially deal with the terms Z 31 and Z 32 . The same arguments as those used before allow us to substitute the denominator by 1 + tδ(t), thereby yielding By Lemma 3, the quadratic forms involved in χ 2 have variance O(K −2 ), and thus can be substituted by their expected mean (see Lemma 6). We obtain tr(P) 1 K tr ( T(u)AT(t)) + o (1). (49) The term χ 1 will be compensated by Z 2 . To see that, observe that the first order of χ 1 does not change if we substitute H by H and P by P. Besides, due to Lemma 5, we can substitute Q (t) by Q(t) and Q (u) by Q(u), hence proving that

t)δ(u)tr ( T(u)AT(t))
(1 + tδ(t))(1 + uδ(u)) + o(1) Replacing A with , one finds a deterministic equivalent We notice that the objective function of (P 2 ) is independent of the norm of a. We can, therefore, select a to maximize the objective function and then adapt the norm to fit the constraint. If we discard the constraint, what remains is a classic Rayleigh quotient [45], which is maximized by the eigenvector a corresponding to the maximum eigenvalue of