Matrix Decomposition for Massive MIMO Detection

Massive multiple-input multiple-output (MIMO) is a key technology for fifth generation (5G) communication system. MIMO symbol detection is one of the most computationally intensive tasks for a massive MIMO baseband receiver. In this paper, we analyze matrix decomposition algorithms for massive MIMO systems, which were traditionally used for small-scale MIMO detection due to their numerical stability and modular design. We present the computational complexity of linear detection mechanisms based on QR, Cholesky and LDL-decomposition algorithms for different massive MIMO configurations. We compare them with the state-of-art approximate inversion-based massive MIMO detection methods. The results provide important insights for system and very large-scale integration (VLSI) designers to select appropriate massive MIMO detection algorithms according to their requirement.


I. INTRODUCTION
Massive MIMO is a key technology for fifth generation (5G) communication systems to achieve very high performance within the available radio spectrum. It is an extension of conventional small-scale MIMO, where a large number of antennas are used at the base station (BS) which serves a large number of users to achieve high data throughput and spectral efficiency [1], [2]. Massive MIMO can provide uniform good service to the user terminals with a high mobility environment. However, the benefit of massive MIMO systems come with the disadvantage of computational complexity. The complexity of symbol detection algorithms grows exponentially as the number of antennas increase in a massive MIMO system [3]. Due to the high number of antennas, conventional linear detectors also require large matrix inversion. Hence, in the past decade, a new class of detectors based on approximate inversion has become a popular choice in very large-scale integration (VLSI) implementation.
Approximate inversion-based detectors (AID) utilize channel hardening properties of massive MIMO systems to calculate the Gramian matrix inversion. This class of detectors work well for certain configurations of massive MIMO systems, for example, when the ratio between numbers of BS antennas and users is large [4]. However, their performance starts to become unstable when this ratio gets smaller. Therefore, the focus has been again shifting towards exact inversion-based linear detectors for their guaranteed performance and robustness [5]- [8]. As the semiconductor technology has also matured greatly over the past decade, the focus has been shifting towards applying better system design than saving every possible logic gate. Telecommunication industries are not always interested in saving logic area by adopting an unstable and risky solutions for their products. Therefore, we envision that exact inversionbased solutions will be more popular in the future for massive MIMO VLSI community.
Matrix decomposition algorithms have been extensively utilized for the matrix inversion procedure of small-scale MIMO detection. They provide better numerical stability over straightforward inversion methods. In addition, they help to achieve a modular design, where the whole procedure can be divided and distributed between different developers. A few popular matrix decomposition algorithms for linear detection are QR, Cholesky and LDL decomposition. Despite the importance of the matrix decomposition algorithms, their complexity analysis for massive MIMO and comparison with existing AID is lacking in the literature. An analysis of matrix decomposition algorithms for small-scale MIMO detection can be found in [9]. A comparison of explicit vs. implicit approximate inversion has been provided in [10]. However, [9] does not address the complexity from massive MIMO perspective and [10] does not focus on the complexity of matrix decomposition algorithms.
In this paper, we address this timely topic and analyze the computation complexity of matrix decomposition algorithms for massive MIMO. We also compare them with the state-ofart AID for the massive MIMO. Several hard output simulations to support our premise regarding different detection mechanisms are also presented in this paper. We believe these results will provide an important guideline to the VLSI designers to select the appropriate algorithms that is suitable to the product requirements.
The rest of the paper is organized in the following way: In Section II, we discuss the system model of a massive MIMO system. We also discuss linear detection, non-linear detection in this section. In Section III, we introduce AID and provide simulation results to demonstrate their instability for different antenna configurations. In Section IV, we present different matrix decomposition algorithms, and their application in different types of detection methods with simulation results. In Section V, we present the complexity analysis and compare with state-of-the-art AID methods. The conclusion is drawn in Section VI.

II. SYSTEM MODEL AND DETECTION METHODS
A massive multi-user MIMO base-station (BS) with N antennas is considered. We assume U single antenna users are transmitting towards the BS, where U ≤ N . Assuming that the channel between the users and BS antennas is frequency flat, the relationship between the transmit and receive vector can be characterized as where y ∈ C B is a received signal vector, x ∈ C U is a transmit symbol vector, H ∈ C B×U is a channel matrix, and n ∈ C B is a circularly symmetric complex white Gaussian noise vector with zero mean and σ 2 noise variance. The system is presented in Fig. 7.
Communication Channel

A. Linear Detection
A MIMO symbol detector determines the transmitted symbol vector x from the received signal vector y. The two most popular linear detection are zero-forcing (ZF) and linear minimum mean-square error (MMSE) equalization. The ZF algorithm does not consider the effect of noise vector n, inverts the channel matrix H to determine the transmitted vector. Therefore, the ZF detection can be expressed as where H † is a pseudo-inverse of H. The ZF detector requires an inversion of the Gramian matrix, G, where G ZF = H H H. The MMSE detector is an improvement over the ZF, which takes the noise into account. MMSE detection can be expressed asx where I U is the U × U identity matrix. The Gramian matrix is modified with a regularization by noise variance for MMSE, i.e.,

B. Non-linear Detection
The earliest massive MIMO detectors are iterative nonlinear detectors, which are initialized with a linear detector, and in subsequent iterations tries to update the results. Example of such detectors include likelihood ascent search (LAS), reactive tabu search (RTS), CHEMP etc. We present a nonlinear detector called ADMM-based Infinity-Norm Detection (ADMIN), proposed in [11], to show how these non-linear detectors for massive MIMO works in general.
ADMIN solves a box-constrained detection problem with a convex optimization method called alternating direction method of multipliers (ADMM). The ADMM method can be used to solve a convex problem by breaking it into smaller problems in an iterative fashion. The first iteration of ADMIN has been presented as, where β is an scaled form of σ 2 I. The values of z and λ λ λ vectors are calculated with two other iterations of ADMIN. When z and λ λ λ are initialized with zeros, (4) resembles the MMSE equation. Thus, the ADMIN detection method is also based on an inversion of a regularized Gramian matrix,

III. INSTABILITY OF APPROXIMATE INVERSION-BASED DETECTION
The exact inversion of the Gramian can be complex when the number of users increases. For example, Gramian of a 16user system (U = 16) will be a 16 × 16 matrix. Several AID became popular in the VLSI community over the past decade due to their less complexity. We present three such detectors in this section and analyze their error-rate performance with Matlab simulations.

A. Neumann Series Approximation
Gramian matrix G can be decomposed into a diagonal matrix (X) and off-diagonal matrix (E) as G = X + E. The Neuman series approximation (NSA) [12] of such a system can be expressed as NSA only use inversions of the diagonal matrix, X, which can be achieved with reciprocals of the diagonal elements. Both precision and complexity of the inversion increases with more iterations.

B. Gauss-Seidel Method
In Gauss-Seidel (GS) method, G can be decomposed as G = D+L+R, where D, L and R are the diagonal component, the strictly lower triangular component, and strictly upper triangular component, respectively. The GS can be used to estimate the transmitted signal vectorx aŝ wherex MF = H H y is a matched filter [13].

C. Conjugate Gradient Method
Conjugate Gradient (CG) is another approximate method used in the MIMO detection which can be formulated aŝ where p (n) is the conjugate direction with respect to the Gramian matrix and α (n) is a scalar parameter which is commonly known as the step size [14].

D. Error-rate Performance
We present error-rate simulation results for exact inversionbased MMSE, and approximate inversion-based NSA, GS, and CG detection methods. We have considered t = 3 iterations of NSA, GS and CG detection methods. The bit error rate (BER) of the detectors with respect to signal-to-noise (SNR) ratios is shown in Figs 2-4. Here, we used 10,000 Monte-Carlo trials for all simulations. The modulation scheme for these simulations is 64QAM. We consider an i.i.d. Rayleigh fading channel between the BS and users. In Fig. 2, the BER for 256-antenna BS and 16 users is presented. Here, the NSA, GS, and CG present very similar performance and perform as well as the exact inversion-based MMSE. In Fig. 3, simulations have been taken for 32 antenna BS and for 16 users. Here, all AID performs poorly and can not detect received symbol vectors properly. In Fig. 4, simulations have been taken for 64 BS antennas and 16 users. Here, NSA and CG detectors can not detect received symbol. However, GS detection method provides good performance in this scenario where the MMSE provides only 2 dB gain over the GS in this scenario. It is evident from the simulations above that all AID works very well and can provide a competitive solution compared to MMSE when the ratio between the BS antennas and number of users is high. Therefore, for such a system, it is possible to select low complexity AID methods. However, RF chains and front-end logic associated with the antennas are expensive. Therefore, having large number of antennas for a very small number of users are not always feasible from design and cost perspective. When the ratio of numbers between BS antennas and users becomes smaller, performance of AID starts to deteriorate as seen in Figs. 3 and 4. In extreme cases, when the ratio becomes close or equal to 1, these detectors do not function at all and fail to detect incoming symbol vectors. In a nutshell, the performance of these detectors is not stable and as the industry has to support various configurations for their product, such detectors might not be the ideal solution for them. Therefore, exact inversion-based detectors are still going to be an attractive solution for industries for years to come.

IV. MATRIX DECOMPOSITION ALGORITHMS
Instead of doing a straightforward matrix inversion, several matrix decomposition methods can be used which are more numerically stable. They also bring a modular design, where an inversion process can be split into parts. Here, we mainly focus on inverting the regularized Gramian matrix, which is used in the MMSE equalization. In this section, we analyze how three popular matrix decomposition algorithms, i.e. QR, Cholesky and LDLT can be applied to find the MMSE solution.

A. MMSE with QR
QR decomposition algorithm decomposes a U × U matrix A as A = QR. Here, Q is an U × U unitary matrix with orthogonal columns and R is an U ×U upper triangular matrix with non-zero diagonal elements. QR decomposition can be applied to obtain the MMSE solution as There are several algorithms to compute the QR decomposition, such as, Gram-Schmidt process, Householder transformations, Givens rotation etc. We focus on Gram-Schmidt process in this paper.

B. MMSE with Cholesky
Cholesky is another popular decomposition algorithm which utilizes a U × U lower triangular matrix L to decompose a U × U matrix A as A = LL H . The MMSE can be calculated asx Here,z can be calculated by forward substitution andx can be solved by backward substitution.

C. MMSE with LDL
LDL is another popular decomposition algorithm which decomposes a U × U matrix A as A = LDL H , where D is a diagonal matrix and L is a lower triangular matrix with zeros as diagonal. The MMSE can be calculated with the LDL asx Similar to Cholesky,z can be calculated by forward substitution andx can be solved by backward substitution.

D. Error-rate Performance
We present error-rate performance of MMSE, ADMIN and single-input multiple-output (SIMO) in Fig. 5 and 6. Similar to earlier simulations, 10,000 Monte-Carlo trials are used. The communication channel between BS and users is assumed to be i.i.d. Rayleigh Fading channel. The MMSE utilizes QR and ADMIN utilizes LDL decomposition in these simulations. We take t = 5 iterations of ADMIN in this simulations. In Fig. 5, the detectors are simulated for 32 BS antennas supporting 32 users transmitting with 64-QAM. It can be noticed that even the exact matrix inversion-based MMSE can not provide a satisfactory performance with high power, i.e., can not bring the error-rate below 10 −2 with approximately 35 dB of SNR. ADMIN provides about 5 dB gain over MMSE, which is a significant margin for communication systems using higher order modulation. Fig. 4 shows, when the ratio of numbers between BS antennas and users is small, non-linear detectors are required for necessary performance gain. Most of these non-linear detectors are based on exact inversion and matrix decomposition algorithms will be crucial for their implementation. The non-linear detectors can also provide very high gain compared to MMSE for other scenarios. In Fig. 5, ADMIN with t = 5 iterations can provide around 10 dB gain over MMSE algorithm. In this case, the difference between the SIMO bound and ADMIN is less than 10 dB.

V. COMPLEXITY ANALYSIS
We analyze the complexity of matrix decomposition algorithms in this section. First, we explain in detail how we computed the complexity of QR decomposition based on classical Gram-Schmidt process. We calculate the complexity of Gram-Schmidt applied on a U × U matrix because the Gramian is always square in the MMSE. Real multiplications dictate the overall complexity of an algorithm as they take significantly more logic than additions or subtractions. Therefore, we focus on real multiplications for comparing the algorithms. Classical Gram-Schmidt algorithm is presented in Algorithm 1.
In line 3 of Algorithm 1, norm of a vector is computed as r i,i = q i 2 . For norm calculation, we need a total of U 2 number of complex multiplications. Here, U number of multiplications are required for a single norm and we need a total of U norms. For q i = q i /r i,i , a total of U reciprocals Algorithm 1 Classical Gram-Schmidt input: A outputs: Q, R 1: Q = A 2: for i = 1, . . . , U 3: for j = i + 1, . . . , U 6: q j = q j − r i,j q i 8: end 9: end needed to compute 1 ri,i and these real-valued reciprocals will be multiplied U 2 times to the complex valued q i . Thus, a total 2U 2 real multiplications are needed for line 4. Next, r i,j = q H i q j implies a vectorvector multiplication where each vector has U complex elements. Therefore, we need U number of complex multiplications for a single q j . Here, we have U number of complex multiplications for r i,j q i . Similarly, we get a total of U (U −1) 2 U multiplication for r i,j q i . The operations are presented in Table I. 3 Therefore, we get the total number of complex multiplication (CM) for Gram-Schmidt as We assume, a complex-complex scalar multiplication hardware unit utilizes four real multiplication units and a real-complex scalar multiplication unit utilize two real multiplication units. Therefore, we get 4U 3 real multiplication from the U 3 complex multiplication. We add this with number of real multiplications (RM) for line 4, which leads to Similarly, we calculate the number of multiplications required for Cholesky decomposition from Algorithm 2. A total of U (U −1) 2 complex multiplication is required to compute the diagonal elements of L (line 3). The off-diagonal elements require a total of U 3 −3U 2 +2U 6 complex multiplication (line 5). Therefore, the total number of complex multiplication is Algorithm 2 Cholesky Decomposition input: A outputs: L 1: for i = 1, . . . , U 2: j = 1, . . . , i − 1 3: Computation of Cholesky also requires a U number of recipro- , which are then multiplied in line 5. Therefore, we require an additional U (U −1) real multiplications. Therefore, the total number of real-multiplication is The complexity of LDL decomposition is very similar to that of Cholesky with an additional 4U (U −1) real multiplications. We provide Table II with more detail of the operations required for the matrix decomposition algorithms. The number of real multiplication, addition and subtractions for QR is significantly higher than Cholesky and LDL decomposition. The complexity of Cholesky and LDL are almost similar except the additional multiplication required for LDL and the square root operations required for Cholesky.

A. Comparison against approximate inversion-based detectors
A comparison of the matrix decomposition algorithms and the AID methods is presented in Table III. All the matrix decomposition algorithms and NSA scale with U 3 while GS and CG methods scale with U 2 in terms of operations. The cubic complexity of NSA originates from the matrix multiplications associated with every iteration t in (X −1 E) t . For example, for t = 3, we have to compute a P * P * P where P = X −1 E. In Fig. 7, matrix decomposition algorithms and the AID methods are compared in terms of numbers of real multiplication required for each algorithm. Here, the AIDs applied t = 3 iterations each. From this figure, the QR decomposition and the NSA are significantly more complex than the other algorithms. Cholesky and LDL provide similar complexity, which is also supported by the operation numbers of Table II. GS and CG methods require the least number of multiplications for their execution.
As Cholesky and LDL provide significantly less complexity over the QR decomposition, they are ideal choice for implementation platforms. These algorithms even provide lower complexity than NSA. The error-rate performance of the CG method is not satisfactory according to Figs. 3 and 4. On the other hand, GS method is very promising as it outperforms other AIDs in terms of error rate while using low number of operations. From Fig. 7, the number of operations required for GS is less than half of Cholesky and LDL. It should be noted that the channel model utilized in this work is i.i.d.  8  8  2432  604  576  Cholesky  8  8  392  42  36  LDL  -8  560  42  36   16   Gram-Schmidt  16  16  17920  4532  4352  Cholesky  16  16  2960  210  136  LDL  -16  3680  210  136   32   Gram-Schmidt  32  32  137216  34660  33792  Cholesky  32  32  22816  930  528  LDL  -32  25792  930   and therefore, the performance of the GS method might suffer in a more realistic channel model. We plan to compare GS, Cholesky, LDL in a more realistic channel model for our future works.

VI. CONCLUSION
We presented an analysis of matrix decomposition algorithms from massive MIMO context and their suitability for VLSI design of massive MIMO detectors. We compared massive MIMO detection algorithms based on approximate inversion and conventional matrix decomposition algorithms. The results will provide a guideline to select the appropriate algorithm for the research community. We concluded that Cholesky and LDL are viable solutions for massive MIMO detection. We also found that the GS can be an attractive approximate inversion-based detector which might provide a good balance between error-rate performance and complexity.