© Springer Nature Singapore Pte Ltd. and Science Press, Beijing, China  2019

Leibo Liu, Guiqiang Peng and Shaojun WeiMassive MIMO Detection Algorithm and VLSI Architecturehttps://doi.org/10.1007/978-981-13-6362-7_2

2. Linear Massive MIMO Detection Algorithm

Leibo Liu¹  , Guiqiang Peng²   and Shaojun Wei³  

(1)

Institute of Microelectronics, Tsinghua University, Beijing, China

(2)

Institute of Microelectronics, Tsinghua University, Beijing, China

(3)

Institute of Microelectronics, Tsinghua University, Beijing, China

 

 

Leibo Liu (Corresponding author)

Email: liulb@tsinghua.edu.cn

 

Guiqiang Peng

Email: pgq13@mails.tsinghua.edu.cn

 

Shaojun Wei

Email: wsj@tsinghua.edu.cn

Massive MIMO signal detection is the key technology of next generation wireless communication (such as 5G) [1], and how to detect the transmitted signal from the mass MIMO system efficiently and accurately is of vital importance. As for massive MIMO signal detection, there are many algorithms to implement the signal detection. Generally, these algorithms can be divided into the linear detection algorithm and the nonlinear detection algorithm according to different calculation methods [2]. Although the linear detection algorithm is less accurate than the nonlinear detection algorithm, it is still a practical signal detection method for massive MIMO system in some cases due to its low complexity. In the linear detection algorithm, the difficulty people often encounter is the calculation to find the inverse matrix of a large-scale matrix, especially when the scale of a massive MIMO system is very large, the algorithm complexity is very high and the corresponding hardware is difficult to implement. Therefore, this chapter introduces several typical linear iterative algorithms for massive MIMO signal detection. Using these algorithms, the iterations between vectors or matrices can be effectively used to avoid direct inversion of large-scale matrices and reduce complexity of the linear detection algorithm. In the following sections, we will introduce Neumann series approximation (NSA) algorithm, the Chebyshev iteration algorithm, the Jacobi iteration algorithm and the Conjugate gradient (CG) algorithm respectively. And the optimization methods of the Chebyshev iteration algorithm, the Jacobi iteration algorithm and the CG algorithm are also introduced to get better linear detection algorithms. In addition, this chapter also compares the complexity and accuracy of these algorithms with those of other massive MIMO signal detection algorithms.

2.1 Analysis of Linear Detection Algorithm

Massive MIMO signal detection algorithms are usually divided into the linear detection algorithm and the nonlinear detection algorithm. The nonlinear detection algorithm, as the name implies, refers to the algorithm that adopts the nonlinear algorithm to recover the transmitted signal s from the received signal y. Such algorithms usually have high accuracy, but the computation complexity is also high. For example, the maximum likelihood (ML) detection is a typical nonlinear detection algorithm [3]. In theory, the ML detection is very ideal for massive MIMO signal detection with high accuracy. However, in the ML detection, the number of cycles required for computation depends largely on the modulation order q and the number of the user antennas (N_t). The total number of cycles is denoted as . Undoubtedly, this result is disastrous undoubtedly because the total number of cycles will still increase considerably even if the modulation order or the number of user antennas increases very little. Therefore, ML detection is not applicable in practice although it is a very ideal detection method in theory, especially in massive MIMO signal detection. In general, the nonlinear detection algorithm is more accurate than the linear detection algorithm in massive MIMO signal detection, but with higher complexity. We will introduce the nonlinear detection algorithm in detail Chap. 4.

Corresponding to the nonlinear detection algorithm, the linear detection algorithm usually estimates the signal s by operating a matrix. The common linear detection algorithms include the zero-forcing (ZF) detection algorithm and the minimum mean square error (MMSE) detection method [4], both of which transform massive MIMO signal detection into linear matrix equation solution by deforming the channel matrix H, that is . According to the massive MIMO channel model in Sect. 1.​2.​1, Eq. (2.1) for both the received signal and the transmitted signal is

(2.1)

where y is the received signal, H is the channel matrix, s is the transmitted signal, and n is the additive noise. The linear detection algorithm of massive MIMO is to simultaneously left multiply the both sides of Eq. (2.1) by conjugate transpose H^H of the channel matrix, neglecting the additive noise n, so as to obtain Eq. (2.2)

(2.2)

If Eq. (2.3) is true, then we obtain Eq. (2.4).

(2.3)

(2.4)

The detection for the transmitted signal s is implemented. However, there is an error in Eq. (2.4) due to the existence of the additive noise n. Based on the above idea, the transmitted signal s can be estimated by a matrix W which makes Eq. (2.5) hold

(2.5)

where the denotes the estimated transmitted signal. In this way, the massive MIMO linear detection is transformed into the estimation of matrix W.

The ZF is a common linear detection algorithm. Its main idea is to neglect the additive noise n in the massive MIMO channel model in the analysis, which will make the massive MIMO detection algorithm much simpler and easier to implement. However, considering that the noise is usually not negligible in actual situation, the result obtained by using this algorithm may not be the optimal solution. For a massive MIMO system with a scale of , the signal received by a receiving antenna at the base station can be expressed as

(2.6)

where , is the ith element of the transmitted signal s, . Now we define a vector that follows Eq. (2.7)

(2.7)

where , will be acted as rows to form a matrix , There is obviously from Eq. (2.7). And combining with Eq. (2.4), we have

(2.8)

In this way, the transmitted signal s can be estimated as

(2.9)

Obviously, when the additive noise is , is strictly satisfied. Because the ZF detection algorithm meets the conditions in Eq. (2.7), it can eliminate the interference between the data sent by different transmitting antennas, and can get relatively accurate detection results when the signal-to-noise ratio is relatively high.

Although there are some deficiencies in the accuracy of the ZF detection algorithm, its derivation process also provides some ideas whether the influence of noise n can be added to the W matrix to simplify the solution of the transmitted signal s using the same method of solving linear matrix. On this basis, a MMSE detection algorithm is proposed.

The MMSE detection algorithm is another typical linear detection algorithm. Its basic idea is to make the estimated signal as close to the real value as possible. In the MMSE detection algorithm, the adopted objective function is

(2.10)

By solving the matrix W based on Eq. (2.10), the following equation can be obtained:

(2.11)

Find the partial derivatives for Eq. (2.11) and set it equal to zero, and get

(2.12)

By solving Eq. (2.12), we can get

(2.13)

where N₀ is the spectral density of noise, E_s is the spectral density of signal, and the Gram matrix is defined. Similar to the ZF detection algorithm, we can deduce the MMSE detection algorithm. It also makes the estimated signal , which is different from the ZF detection algorithm. in the ZF detection algorithm and in the MMSE detection algorithm. Both algorithms can estimate the transmitted signal s in massive MIMO detection.

Whether in the ZF detection algorithm or in the MMSE detection algorithm, matrix inverse is involved. The scale of the channel matrix H is usually very large in massive MIMO system. But the matrix inversion is more complex, which is usually difficult to achieve in the actual signal detection circuit. In order to avoid the huge complexity of matrix inversion, many algorithms have been put forward to reduce the algorithm complexity so as to reduce the perplexities caused by large-scale matrix inversion. Among these algorithms, the typical ones are the NSA algorithm, the Chebyshev iteration algorithm, the Jacobi iteration algorithm, the CG algorithm and so on. The above four common algorithms will be introduced below in detail.

2.2 Neumann Series Approximation Algorithm

In the massive MIMO system, the number of receiving antennas is usually much greater than that of the user antennas [5], that is, N_r is much larger than N_t. Since the elements in channel matrix H are independent identically distributed, and the real and imaginary parts are subject to the Gaussian distribution with parameter N(0, 1), both Gram matrix G and MMSE matrix are diagonally dominant matrices. Gram matrix G tends to be a scalar matrix when Nr tends to infinity, i.e., [5]. Using this property, we can simplify the inverse process of matrix in the linear massive MIMO detection algorithm.

Because of the diagonal dominance of the MMSE matrix A, if taking the diagonal elements of matrix A out and denoting it as matrix D, it is very easy to seek the inverse of matrix D. As mentioned above, the condition that A tends to D is that the number of user antennas N_t tends to infinity. However that is not actually possible. Even it is often far from that condition. Thus, we use to approximate the matrix A and finding the inverse of the matrix A will lead to relatively large errors.

2.2.1 Algorithm Design

Using Neumann series [6] can get the exact matrix inverse. The Neumann series expands the inverse matrix A⁻¹ of the MMSE matrix A into

(2.14)

In the above Eq. (2.14) the matrix X need to satisfy

(2.15)

The matrix A is decomposed into the sum of diagonal elements and non-diagonal elements . Bring the sum into Eq. (2.14), we can get

(2.16)

According to Eq. (2.15), if the condition is satisfied, the expansion of Eq. (2.16) will definitely converge.

Expanding the matrix A⁻¹ into a sum of numerous terms is not practical. In order to apply Neumann series to the linear massive MIMO detection algorithm, it is necessary to use NSA to estimate Eq. (2.16) [7]. The main idea of solving inverse matrix with NSA is to take out the first K items of Neumann series in Eq. (2.16), then the expression of calculating the first K terms Neumann series is

(2.17)

By using Eq. (2.17), Neumann series with a certain number of terms can be calculated to change the infinite number of terms into the finite number of terms, so as to reduce the computational complexity and estimate the inverse of MMSE matrix A. By approximating A⁻¹, an approximate MMSE equilibrium matrix can be obtained. A⁻¹ can be expressed in different expressions according to the number of selected items K. When K = 1, , now . When K = 2, , the operation complexity is . When K = 3, is

(2.18)

The computational complexity of Eq. (2.18) is , which is comparable to the actual computational complexity of the inverse matrix, but the approximate operation of Eq. (2.18) is less. When K > 4, the complexity of the inverse computation of the actual matrix may be lower than that of the approximate algorithm.

2.2.2 Error Analysis

Obviously, the sum of the preceding K items in Eq. 2.17 can cause errors, which can be expressed as

(2.19)

By using Eq. (2.17) to estimate the transmitted signal and substituting Eq. (2.17) into Eq. (2.19), the expression about the error of the transmitted signal of Eq. (2.20) can be obtained

(2.20)

By taking the second norm of the second term of Eq. (2.20), we have

(2.21)

We can see from Eq. (2.21) that if the condition of inequality (2.22) is satisfied in Eq. (2.21), the approximate error exponent approaches zero as the number of terms K increases, and it can be proved that inequality (2.22) is a sufficient condition for the convergence of formula (2.16).

(2.22)

Now it is necessary to prove that Eq. (2.16) converges when the scale of the massive MIMO system satisfies the condition that N_r is much greater than N_t, the elements in the channel matrix H are independent and identically distributed, and all of them obey the complex Gaussian distribution with parameter N (0, 1). More specifically, it is necessary to prove that the condition of convergence of Neumann series and the condition of minimal error in Eq. (2.21) are only related to N_t and N_r. Here is a theorem and corresponding proof below.

Theorem 2.2.1

When , the elements in the channel matrix H are independent of each other and satisfy the complex Gaussian distribution with its variance being 1, the following expression is obtained

(2.23)

Proof

To prove Theorem 2.2.1, we need to give three other lemmas and their proofs.

Lemma 2.2.1

Let , independent and identically distributed, and make them satisfy the complex Gauss distribution with its variance being 1.

(2.24)

Proof

First we can get

(2.25)

The operation process in Eq. (2.25) can be described as the following steps. The nonzero term can be expressed as and after the Quadratic term is expanded, in which there are a total of items and items . According to , , we can get the conclusion of Lemma 2.2.1.

Lemma 2.2.2

Let , and is independent and identically distributed, obey the complex Gauss distribution with the variance being 1, and , then

(2.26)

Proof

First, g is rewritten as

(2.27)

Among them, is independent and identically distributed, also follows the real Gauss distribution with the mean being 0 and the variance being 1. Therefore, obeys the inverse distribution with free degrees. The inverse distribution with free degrees corresponds to the inverse Gaussian distribution with free degrees. The fourth inverse distribution can be obtained by Eq. (2.28):

(2.28)

Then we have the conclusion from Eq. (2.26).

Lemma 2.2.3

Let , and the elements in the channel matrix H satisfy the complex Gaussian distribution with the independent identically distributed zero mean and unit variance.

(2.29)

Proof

The normalized Gram matrix G corresponds to matrix A, . Therefore, the element of ith row and jth column in matrix A can be expressed as

(2.30)

where is the element of ith row and jth column of the matrix G. Therefore, we can find the inequality in Eq. (2.31):

(2.31)

Then use the Cauchy–Schwartz inequality for Eq. (2.31), there is

(2.32)

Use Lemmas 2.2.2 and 2.2.3 in the calculation of first and second moments respectively, get the following expression:

(2.33)

Now let uss prove Theorem 2.2.1. Using Markov’s inequality, we obtain

(2.34)

Combining the upper bound of and in Lemma 2.2.3, we can get the conclusion of Theorem 2.2.1.

As we can see from Eq. (2.23), when , the probability of the condition (2.22) satisfied is greater. This theorem shows that Neumann series converges with a certain probability, the greater the ratio of , the greater the probability of convergence. In addition, the theorem also provides α condition for minimizing the error of residual estimation and when alpha is less than 1, the greater the number of terms K selected, the greater the probability of convergence.

2.2.3 Complexity and Block Error Rate

The NSA reduces the computational complexity of inverse matrix. Now we discuss the advantages and limitations of the NSA from two aspects of computational complexity and block error rates. Here, we only consider the case that is 64, 128, and 256 respectively.

In the exact algorithm to solve the inverse matrix, the Cholesky decomposition (CHD) algorithm [8] has lower complexity than other exact solutions, such as direct matrix inversion, QR decomposition, LU decomposition [3, 9], etc. Therefore, the CHD algorithm can be selected as the object to be compared with the NSA. The complexity of the CHD algorithm for inverse matrix solution is in the in the range, while the complexity of the NSA for inverse matrix is in the range and the range respectively when K = 1 and K = 2. When K > 3, the complexity of the NSA mainly comes from the multiplication between large matrices, and the complexity increases linearly with the value of K. For example, when K = 3, there is one multiplication between large matrices once in the algorithm, and when K = 4, there are two multiplications between large matrices, that is, when K > 3, you need to calculate K − 2 multiplications between large matrices in the NSA. So when K > 3, the complexity of the NSA algorithm is . It can be seen that when , the NSA has no advantage in complexity compared with the CHD algorithm.

The complexity of an algorithm mainly depends on the number of real number multiplications in the algorithm. Figure 2.1 depicts the variation curve of the number of real number multiplications with the number of antennas in the CHD algorithm and the NSA algorithm with different K values. Figure 2.1 shows that the complexity of the NSA algorithm is lower than that of the CHD algorithm at , while the complexity of the Newman series approximation algorithm is higher when K > 3.

Fig. 2.1

The Relationship between the number of the user antennas and the number of the real number multiplications in the algorithm.

© [2018] IEEE. Reprinted, with permission, from Ref. [7]

NSA is the solution of approximating matrix inversion by taking the first K terms of the Neumann series. Obviously, the more the number of terms is taken, the closer it is to the exact result, but the cost is the increase in complexity. Thus, the precision and complexity are a pair of contradictions. In order to compare the block error rate between Neuman series approximation and Holesky decomposition algorithm, the uplink of massive MIMO system is selected here. At the base station, MMSE detection using the above NSA and CHD algorithm is adopted, and is defined. Figure 2.2 shows at a different , the block error rate of the NSA and the CHD algorithm when is equal to 4, 8, and 12 respectively.

Fig. 2.2

Block error Rate Curve at a , b , c .

© [2018] IEEE. Reprinted, with permission, from Ref. [7]

As we can see from Fig. 2.2, when K = 1 or K = 2, the NSA algorithm has a large block error rate, when the number of antennas at the base station is large, it can make up for part of block error rate. Considering the requirement of 10% block error rate [10] in LTE, it is not suitable for practical applications when K = 1 and K = 2 with modulation order of 64-QAM. The simulation result shows that the block error rate is less than 10⁻² when K = 1, , , and the number of terms required by the NSA algorithm is fewer when the modulation order is 16-QAM. When the modulation order is 64-QAM and K = 3, the result of the NSA algorithm is close to that of the CHD algorithm. For example, when K = 3, and K = 3, , , and the block error rate is 10⁻², the SNR loss of NSA algorithm is less than 0.25 dB. Therefore, when the ratio of the massive MIMO system is large, the NSA algorithm is used to find the inverse of the matrix and the term K is 3, so that the lower block error rate can be reduced under the condition of low computational complexity.

In summary, in the massive MIMO system, when the ratio is small, the CHD algorithm and other exact arithmetic are required to seek the inverse of the matrix. When ratio is large, the NSA algorithm can be used to approximate the inversion. By using the NSA algorithm, we can find the relatively accurate inverse matrix results with low computational complexity. This makes the NSA an efficient and accurate method for massive MIMO detection in some specific cases.

2.3 Chebyshev Iteration Algorithm

2.3.1 Algorithm Design

The Chebyshev iteration algorithm [11] is an algorithm for solving matrix equation by using iteration computation to avoid large matrix inversion. Its basic iteration form is

(2.35)

where σ is the correction matrix and K is the number of iterations. The σ can be expressed as

(2.36)

(2.37)

where denotes the residual vector, ρ^(K) and φ^(K) are the Chebyshev polynomial parameters for two iterations, and β is an iterative parameter related to the eigenvalue of matrix A. Therefore, the Chebyshev iteration can be used to solve the linear equation in the MMSI of massive MIMO detection, so as to avoid the operation complexity caused by large matrix inversion. In this section, for convenience, set , so there is .

Although the Chebyshev iteration can be used in the MMSE of massive MIMO detection, it still faces some challenges. First of all, since the parameters such as β, ρ^(K) and φ^(K) are related to the eigenvalues of matrix A, it is difficult to calculate these parameters. Second, at the beginning of iteration, it is necessary to solve matrix A. Matrix A involves multiplication of large-scale matrix, which will consume a lot of hardware resources. Third, different initial values of iteration affect the convergence rate of the algorithm, and how to determine a good initial value is also the challenge of the algorithm. To solve the above problems, the Chebyshev iteration is optimized in this section to make it more suitable for massive MIMO signal detection.

According to Eq. (2.35), the iteration form in MMSE can be written as

(2.38)

where is the initial value of iteration, and the parameter satisfies:

(2.39)

In order to reduce computational complexity, matrix A and in Eq. (2.39) can be split into

(2.40)

In Eqs. (2.36) and (2.37), the parameters ρ^(K) and φ^(K) can be expressed as

(2.41)

(2.42)

where T is Chebyshev polynomial, and α is the parameter related to the eigenvalues of matrix A. According to Chebyshev polynomial [11], the expression of T is

(2.43)

(2.44)

(2.45)

Combining formula (2.41) and (2.42), we can get Eqs. (2.46)–(2.49):

(2.46)

(2.47)

(2.48)

(2.49)

where and satisfy:

(2.50)

(2.51)

where λ_max and λ_min are the maximum and minimum eigenvalues of matrix A respectively. Since computing the eigenvalues of matrix A is complicated, an approximation is adopted here. As N_r and N_t increase, λ_max and λ_min can be approximately expressed

(2.52)

(2.53)

So far, all the parameters used in the Chebyshev iteration can be expressed in the scale parameters of the channel matrix H in the massive MIMO system. According to Eq. (2.38), we still need an iterative initial value if we want to use the Chebyshev iterative algorithm to estimate the signal s. Theoretically, although using any initial values we can get the final estimation, the convergence rates of the algorithm corresponding to different initial values are not the same. A good initial value can make the algorithm converge faster and generate the desired results, achieving twice the result with half the effort.

As described in Sect. 2.2, in the massive MIMO system, the number of receiving antennas is often much larger than the number of the user antennas, i.e., , the elements in the channel matrix H are subject to Gaussian distribution that is independent and identically distributed and with parameters N(0, 1), so that matrix A is a diagonally dominant matrix and satisfies

(2.54)

So the initial value can be approximated as

(2.55)

This initial value enables the Chebyshev iteration to achieve a faster convergence rate, and the computational complexity of the initial value is very low. In addition, the computation of the initial value can be executed in parallel.

In massive MIMO signal detection, it is often necessary to output the log-likelihood ratio for the use of the next stage circuit, so it is necessary to discuss how to use Chebyshev iteration to find the approximate log-likelihood ratio. The estimated transmitted signal can be expressed as

(2.56)

Set and , they can be used to solve equivalent channel gain and NPI respectively. In combination with Eqs. (2.37) and (2.38), the estimated received can be expressed as.

(2.57)

Set , the iteration in Eq. (2.57) can be approximated as the inverse of matrix A, X, and Y. For example

(2.58)

where , and all are diagonal matrices. Similarly, matrices X and Y are shown as

(2.59)

(2.60)

The equivalent channel gain μ_i and NPI can be approximately expressed as

(2.61)

(2.62)

The signal-to-interference-plus-noise ratio (SINR) can be calculated by combining Eqs. (2.58)–(2.62). However, this algorithm has a high computational complexity. An algorithm based on initial eigenvalue solution is presented here to solve LLR, which can reduce the computational complexity of LLR [12]. The LLR is shown as below.

(2.63)

where the parameter meets

(2.64)

where indicates the SINR of the ith user, and are the set of modulation points in the constellation diagram, and the bth bits of and are 0 and 1 respectively. For the sake of convenience, write as expressed as a form of linear equation. Here, the approximate SINR no longer depends on the result of and .

Based on the above analysis, we can obtain the optimized Chebyshev iterative algorithm that approximates an MMSE detector in a massive MIMO system and name it as parallelizable Chebyshev iteration (PCI). The algorithm [13] is shown in Algorithm 2.1.

Algorithm 2.1

The parallelizable Chebyshev iteration (PCI) algorithm for soft-output MMSE detection

2.3.2 Convergence

Now let us discuss the problem of the convergence rate of the Chebyshev iteration. After K steps, the approximation error can be expressed as [14]

(2.65)

In the above equation, satisfies

(2.66)

So the error can be expressed as

(2.67)

From the previous description, we know that matrix A is a diagonally dominant matrix, so there is

(2.68)

where S is a complex matrix of and meets the condition . J is an upper triangular matrix. Now two lemmas are proposed and their corresponding proofs are given, followed by corresponding conclusions.

Lemma 2.3.1

In the massive MIMO system, there is

(2.69)

where is the Kth normalized Chebyshev polynomial.

Proof

The Chebyshev polynomials in Eq. (2.45) can be rewritten as [15]

(2.70)

Combining Eqs. (2.66) and (2.70), the Chebyshev polynomial is converted to

(2.71)

On this basis, using the identity to transform Eq. (2.71) into

(2.72)

where and meet:

(2.73)

(2.74)

From Eqs. (2.50)–(2.53), the parameters satisfy and . As a result. Now satisfies

(2.75)

When the iteration number K increases, the value is limited, so there is

(2.76)

Considering that is an operator similar to that in Eq. (2.72), set . When K is very large, we can get

(2.77)

If the eigenvalues of A satisfy , then the value of keeps a constant, so the eigenvalues of V and satisfy

(2.78)

From Eqs. (2.50)–(2.53), Eqs. (2.76)–(2.78). The following equation is obtained:

(2.79)

Therefore, the conclusion of Lemma 2.3.1 is obtained.

According to Eqs. (2.52) and (2.53), when remains constant but increases, the maximum eigenvalue and the minimum eigenvalue of matrix A will approach to . It is proved in Lemma 2.3.1 that the value is less than 1 and will decrease as the ratio of to increases. According to Eqs. (2.67) and (2.68), the estimation error is very small when the number of iterations K is limited. And from Lemma 2.3.1, we know that the estimation error will decrease with the increase of the number of iterations K. Therefore, we can get

(2.80)

(2.81)

In other words, in massive MIMO signal detection, using Chebyshev iterative algorithm to estimate the transmitted signal s, the calculation error is very small, even close to zero.

Lemma 2.3.2

In the massive MIMO system, and are satisfied, where , and are respectively the normalized Chebyshev polynomials corresponding to the Chebyshev iteration algorithm, the CG algorithm and the NSA algorithm.

Proof

The convergence rate of the Chebyshev iteration algorithm is

(2.82)

In order to make the convergence faster, should be as great as possible in Eq. (2.82). The problem is transformed to find the minimization of the maximum value by Eqs. (2.76) and (2.78), i.e.,

(2.83)

Set , and combining with Eq. (2.51), we can find the minimization of the maximum value . Thus can be summarized as

(2.84)

, the of the CG algorithm can be summarized as

(2.85)

We can solve from Eq. (2.85). This inequality indicates that compared with CG algorithm and Chebyshev iteration algorithm, has a smaller maximum value using Chebyshev iteration algorithm.

In combination with Eqs. (2.50)–(2.53) and (2.73), the of Chebyshev iteration algorithm can be approximated to

(2.86)

According to literature [7], The of NSA algorithm is

(2.87)

Equation (2.87) indicates that compared with the Chebyshev iteration algorithm and the NSA algorithm, has a smaller maximum value.

Combined with Eqs. (2.68), (2.82) and (2.83), a smaller maximum value of results in a faster convergence rate. Therefore, as we know from Lemma 2.3.2, the convergence rate of Chebyshev iteration algorithm is faster than CG algorithm and NSA algorithm.

2.3.3 Complexity and Parallelism

After discussing the convergence and convergence rate of the Chebyshev iteration algorithm, we will analyze the computational complexity of the Chebyshev iteration algorithm. As computational complexity mainly refers to the number of multiplications in the algorithm, we can count the number of multiplications required by the algorithm to evaluate the computational complexity. In the Chebyshev iteration algorithm, the first part of the calculation is the parameter calculation based on Eqs. (2.46)–(2.53). Because of the fixed scale of the massive MIMO system, these parameters only need to be computed once and stored in memory as constants when calculating multiple sets of data so the actual number of multiplications is negligible. The second part is to match the filter vector and find the approximate initial solution by using the eigenvalues. The matrix needs to be multiplied by vector y to get the vector , then calculate of its value. There are actual multiplications in this process. There are three steps in the third part of the computation. First, multiply the matrix H by the transmitted vector . Then solve the residual vector and finally find the correction vector . These three steps require , and real number multiplications respectively. The last part of the computation is the calculation of approximate log-likelihood ratio. This step requires multiplications when the modulation order is 64-QAM. Therefore, the total number of multiplications required by the Chebyshev iteration algorithm is .

As mentioned before, the method to directly find the inverse matrix will result in the computational complexity increase. Such computation is time-consuming, which makes hardware implementation difficult. A large-scale matrix multiplication takes 36 clock cycles when many pieces of hardware are parallel. These defects affect the energy and area efficiency of the detector hardware. In addition, the PCI contains the initial value, whose computation process can be considered as one iteration. Therefore, in order to balance, the iteration number K of PCI contains a calculation of the initial value and K-1 iterations. The comparison of computational complexity of different methods is listed in Table 2.1. The computational complexity of accurate MMSE detection methods such as the CHD algorithm is . The computational complexity of NSA increases with the increase of the number of iterations K, and the specific manifestation is that when K < 3 the complexity is , when K = 3 the complexity increases to and the computational complexity of NSA is higher than that of accurate MMSE detection when K > 3. Usually, to ensure the detection accuracy, the number of the iterations of NSA K should be greater than 3. In addition, the NSA in literature [7] adopts the explicit method of large-scale matrix multiplication as matrix multiplied by matrix H, which consumes a lot of computing resources, i.e. . Since the results of direct inversion can be used repeatedly in downlink, if the parallelizable Chebyshev algorithm is used for calculation, the complexity will be doubled, but still lower than the results in literature [7]. The following methods directly compute matrix A and indirectly compute A⁻¹, including implicit version of the Neumann series (INS) approximation and implicit version of CG [16].The proposed PCI can also be modified to compute the explicit version of the matrix A. Table 2.1 lists the associated computational complexity of PCI, which achieves lower or equal complexity compared with literatures [12], [16] and [17]. Comparing with PCI under the same conditions, the least square conjugate gradient least square (CGLS) algorithm [18] and the optimized coordinate descent (OCD) algorithm [19] implement large-scale matrix multiplication and inversion in a completely implicit manner. Compared with the CGLS algorithm and the OCD algorithm, the PCI has advantages in computational complexity.

Table 2.1

Computational complexity analysis

Algorithm	K = 2	K = 3	K = 4	K = 5
NSA [7]	2NrNt2 + 6NtNr + 4Nt2 + 2Nt	2NrNt2 + 6NtNr + 2Nt3 + 4Nt2	2NrNt2 + 6NtNr + 6Nt3	2NrNt2 + 6NtNr + 10Nt3-6Nt2
INS [17]	2NrNt2 + 6NtNr + 10Nt2 + 2Nt	2NrNt2 + 6NtNr +14Nt2 + 2Nt	2NrNt2 + 6NtNr +18Nt2 + 2Nt	2NrNt2 + 6NtNr +22Nt2 + 2Nt
GAS [12]	2 NrNt2 + 6NtNr +10Nt2-2Nt	2NrNt2 + 6NtNr +14Nt2-6Nt	2NrNt2 + 6NtNr +18Nt2-10Nt	2NrNt2 + 6NtNr +22Nt2-14Nt
CG [16]	2 NrNt2 + 6NtNr +8Nt2 + 33Nt	2NrNt 2 + 6NtNr + 12Nt2 + 49Nt	2NrNt 2 + 6NtNr + 16Nt2 + 65Nt	2NrNt 2 + 6NtNr + 20Nt2 + 81Nt
CGLS [18]	16NtNr + 20Nt2 + 32Nt	20NtNr + 28Nt2 + 48Nt	24NtNr + 36Nt2 + 64Nt	28NtNr + 44Nt2 + 80Nt
OCD [19]	16NtNr +4Nt	20NtNr + 6Nt	32NtNr + 8Nt	40NtNr + 10Nt
PCI	2NrNt 2 + 6NtNr + 8Nt2 + 8Nt	2NrNt 2 + 6NtNr + 12Nt2 + 12Nt	2NrNt 2 + 6NtNr + 16Nt2 + 16Nt	2NrNt 2 + 6NtNr + 20Nt2 + 20Nt
	12NtNr + 7Nt	20NtNr + 13Nt	28NtNr + 19Nt	36NtNr + 25Nt

Now consider the parallelism of PCI. In the Gauss–Seidel (GAS) algorithm [12], there is a strong correlation between elements, which means that the method has low parallelism. In the calculation of GAS algorithm, when calculating at the Kth iteration, it requires , and , of the previous K iterations. Moreover, in the CHD algorithm as well as NSA, GAS, and CG methods, there is a strong correlation between large-scale matrix inverse and multiplication, which requires the calculation of large-scale matrix multiplication first. This results in their reduced parallelism. PCI parallelism is an important problem in algorithm design and hardware implementation. According to the algorithm, when calculating and correcting vector and residual vector , each element in them can be done in parallel. Besides this, we can find that the implicit method reduces the correlation between large-scale matrix multiplication and inverse calculation, and improves the parallelism of the algorithm.

2.3.4 Bit Error Rate

To evaluate the performance of PCI, BER’s simulation results are compared with NSA, CG, GAS, OCD and the Richardson (RI) algorithm below. Furthermore, the BER performance comparison also includes the MMSE algorithm based on the CHD. SNR is defined at the receiving antenna [7]. The number of iterations K of PCI is an initial value calculation of the K-1th iterations.

Figure 2.3 shows the BER performance comparisons (PCI and other methods) and the simulation results when , with SNR set to 14 dB. As we can see from the figure, the PCI only uses a small number of iterations as increases, and achieves the near-optimal performance (compared with accurate MMSE). This result demonstrates the reasonability of the reduced computational complexity. Figure 2.3 also exhibits the performance of the NSA. When the number of iterations is relatively small, the detection accuracy loss cannot be ignored. However when the number of iterations is large, the system consumes a lot of hardware resources. Hence, this figure validates that PCI in massive MIMO system is superior to that in NSA.

Fig. 2.3

BER simulation results of various algorithms at N = 0, SNR = 14 dB.

© [2018] IEEE. Reprinted, with permission, from Ref. [13]

Figure 2.4 shows the BER comparisons between the PCI using the initial value for the iteration and the traditional zero vector initial solution. In this simulation, the numbers of antennas are 16 and 128 respectively. The initial solution based on the eigenvalue approximation achieves low detection accuracy loss with the same number of iterations. The simulation result in the PCI when K = 3 is very close to the BER when K = 4, which means that the initial value for the iteration reduces the amount of computation while maintaining the similar detection accuracy.

Fig. 2.4

BER performance comparisons between PCI after updating initial value and PCI for conventional zero vector value.

© [2018] IEEE. Reprinted, with permission, from Ref. [13]

Figure 2.5 exhibits the performance of PCI, CG, OCD, GAS, RI, and NSA. It also provides the MMSE similar to the CHD algorithm for reference. The comparison of the three simulation results shows that PCI achieves the near-optimal performance under different antenna configurations. To implement the same BER in PCI, SNR is required to be almost identical to GAS and OCD methods, but smaller than RI and NSA. When the ratio of to is small, CG is slightly better than PCI. When the ratio of to is larger, PCI has better detection performance. For example, in Fig. 2.5 (c), when K = 3, the SNR required to implement 10⁻⁶ BER is 17.25 dB, which is close to the precise MMSE (16.93 dB), GAS (17.57 dB), OCD (17.67 dB), and CG (17.71 dB). By comparison, the SNRs required for RI and NSA are 18.45 and 19.26 dB, respectively.

Fig. 2.5

BER comparisons between PCI and other algorithms. a N_r = 64, N_t = 16, b N_r = 90, N_t = 16, c N_r = 128, N_t = 16, d N_r = 162, N_t = 16

© [2018] IEEE. Reprinted, with permission, from Ref. [13]

2.3.5 Analysis on Channel Model Impact

The channel of the massive MIMO system also affects the algorithm. The Kronecker channel model [20] is often used to evaluate the performance because it is more practical than the Rayleigh fading channel. In the Kronecker channel model, the elements of the channel matrix satisfy , where d(z) indicates the channel fading (such as path fading and shadow fading). Another significant feature of the Kronecker channel model is the consideration of channel correlation. Specifically, R_r and R_t denote the channel correlation of the receiving antenna and the transmitted antenna respectively. This part is also based on the Kronecker channel model, and the channel H can be expressed as

(2.88)

where H_i.i.d. is a random matrix, the elements of which are independent and identically distributed and are subject to the complex Gaussian distribution with a zero mean and a unit variance.

In the simulation, the radius of each hexagonal region is r = 500 m, and the users’ locations are independent and random. The independent shadow fading C satisfies

(2.89)

Considering the channel fading variance , where , κ and ||·|| are the base station location, path loss index and Euclidean norm respectively. The simulation adopts the following assumption: κ = 3.7, and the transmitted power is .

Now we will discuss the influence of the Kronecker channel model on eigenvalue approximation. When the channel is independent and identically distributed Rayleigh fading, the value of the diagonal element of matrix A approximates to β, ξ is the channel correlation coefficient. When the correlation coefficient increases (the Kronecker channel model), the approximate error also increases slightly. Therefore, the eigenvalue approximation is still applicable to more practical channel models, such as the Kronecker channel model. The influence of different channel models on computational complexity is also considered. In order to satisfy the requirement of low approximation error, the number of iterations of the proposed PCI should be slightly increased for enhanced channel correlation and increased eigenvalue approximation error. Therefore, there is a slight increase in the computational complexity of PCI, which is a limitation of this method. As discussed earlier, there are three ways to compute large-scale matrix multiplication and inverse. When channel frequency is flat and changes slowly [21], as channel hardening becomes obvious, the result of matrix multiplication and inverse of explicit and partially implicit methods can be partially reused, such as NSA, INS, GAS, and CG. However, these methods also have some limitations. For example, the high computational complexity at the beginning of the detector can affect subsequent operations (such as matrix decomposition and reciprocal of diagonal elements), which need to start only after matrix multiplication is completed. In addition, the hardware utilization of large-scale multiplication is not high, which reduces the energy and area efficiency of detector. These methods (i.e., CGLS, OCD, and PCI) use implicit methods to solve large-scale matrix multiplication and inversion. The computation of these methods is much lower and the parallelism is higher than the explicit method. However, the results cannot be reused in the case of low frequency selectivity due to small-scale fading [22] average and channel hardening, which is another limitation of PCI. Finally, Fig. 2.6 shows the impact of large-scale fading and spatial correlation on the MIMO channel (Kronecker channel model), which is an important problem in the actual MIMO system. Simulation results show that compared with MMSE, the accuracy loss of PCI is less. As the channel correlation increases (the channel correlation coefficient increases), the number of iterations of all methods increases to reduce the approximation error. With the same number of iterations, PCI achieves lower error compared to the NSA and RI methods. Although PCI has similar detection performance compared with CG, OCD, and GAS, its main advantages are higher parallelism and lower computational complexity than the other three methods. In a word, PCI is superior to other methods under more practical channel model conditions.

Fig. 2.6

BER comparisons of various algorithms under the Kronecker channel model.

© [2018] IEEE. Reprinted, with permission, from Ref. [13]

2.4 Jacobi Iteration Algorithm

2.4.1 Weighted Jacobi Iteration and Convergence

This section introduces an optimized Jacobi iteration algorithm named the weighted Jacobi iteration (WeJi) algorithm. As discussed earlier, matrix G and matrix W are diagonally dominant matrices. Here, we can decompose matrix W into , where P is a diagonal matrix, and Q is a matrix with the diagonal elements as 0. Using the WeJi to solve the linear equation, the transmitted signal can be estimated as

(2.90)

where , is an iteration matrix, K is the number of iterations, and is the initial solution. In addition, , it plays a crucial role in the convergence and convergence rate of WeJi. In WeJi, the range of the parameter ω is set as [9]. Because P is a diagonal matrix, its inverse matrix is very easy to find, the computational complexity of WeJi is greatly reduced.

The initial value of the iteration will affect the convergence and convergence rate of the iteration, so we need to find a good initial value when we use the WeJi to solve the massive MIMO signal detection problem. Here, the initial value with low computational complexity can be obtained by using NSA. Therefore, we can set the initial value of iteration as Eq. (2.91), which can make WeJi converge at a faster rate.

(2.91)

In addition, it is because of the increase of algorithm convergence rate that this initial value also reduces hardware resource consumption and increases data throughput rate.

According to Eq. (2.90), when the iteration number K tends to be infinite, the error of the transmitted signal estimated by using the WeJi is [7]

(2.92)

Here, obviously . So the convergence rate of WeJi is

(2.93)

where is the spectral radius of the matrix B. We can see that when is very small the convergence rate of the algorithm is higher. As for the convergence rate of the WeJi, two lemmas and their corresponding proofs are given here.

Lemma 2.4.1

In massive MIMO systems, , and are the iterative matrices of the WeJi and NSA respectively.

Proof

The spectral radius of is defined as

(2.94)

In the WeJi, and ω are close to 1. can be converted to , so Eq. 2.94 can also be written as

(2.95)

In the Newman series approximation, there are

(2.96)

Combined with Eq. (2.95), we can get

(2.97)

Lemma 2.4.1 shows that the WeJi converges faster than that of the Newman series approximation. Without loss of generality, l₂ norm is used to estimate the error of iteration, for example

(2.98)

According to the above expression (2.98), if is satisfied, the approximate error of WeJi will exponentially approach 0 with the increase of the number of iterations K.

Lemma 2.4.2

In the massive MIMO system, the probability of satisfies:

(2.99)

Proof

We can get the following formula according to the Markov inequality:

(2.100)

Note that if the parameter ω satisfies and approaches 1 in WeJi, then is close to 0, and satisfies , so the effect of can be ignored. The formula (2.99) satisfies

(2.101)

So the probability of is related to . Consider , the element of the ith row and jth column in the matrix A satisfies

(2.102)

As a result, can be expressed as

(2.103)

Then using the Cauchy–Schwartz inequality for Eq. (2.103), and we have

(2.104)

Obviously, there are two key items and . In massive MIMO system, the diagonal elements of the matrix A are close to , so term can be approximated as

(2.105)

Now consider . We can express it according to formula (2.102) as follows:

(2.106)

Let X and μ satisfy Eqs. (2.107) and (2.108), where is the average of .

(2.107)

(2.108)

Therefore, there is

(2.109)

The matrix C is the covariance matrix. Note that the elements in the channel matrix H are independent and identically distributed and follow , so when , and are also independent and identically distributed and obey . Therefore, expression (2.109) can be written as follows:

(2.110)

The expression above shows that when , obeys and obeys . Therefore, the probability density function of the random variable is

(2.111)

where Γ is the gamma function [23]. We can get , . When , since and are independent, and , , . After neglecting the zero items, can be expressed as

(2.112)

By substituting Eqs. (2.105) and (2.112) into Eqs. (2.101) and (2.104), we can obtain the conclusion in Lemma 2.4.2.

Lemma 2.4.2 demonstrates that when is fixed, the probability of increases with the increase of . Because of in the massive MIMO system, the probability of is close to 1.

2.4.2 Complexity and Frame Error Rate

Some analysis needs to be done for the WeJi. The WeJi first performs the calculation of the matrices and . In order to facilitate the calculation for the WeJi, the initial solution should be gotten as soon as possible. So the vector T and the matrix R should be ready in the allocated time. The initial value of the WeJi needs to be calculated based on Eq. (2.91). It is worth noting that the architecture of the initial solution can be reused in the next iteration block when considering the hardware design. Finally, the iteration of the final value is executed in Eq. (2.90) of the algorithm. In the iteration section, the matrix multiplication is implemented by vector multiplication. All elements of the vector can be executed in parallel. The computational complexity of the complex matrix inversion can reduce the number of iterations. In addition, the weighted parameters reduce the number of iterations and achieve similar performance, which also reduces the computational complexity of the detector. Now we are going to compare the WeJi with the recently developed algorithms in terms of computational complexity, parallelism, and hardware design realizability. Because the MMSE and the WeJi need to compute matrices G and , the work mainly focuses on the computational complexity [7, 12, 24, 25] of the matrix inversion and the LLR calculation. The computational complexity is estimated based on the number of real number multiplications required, and each complex multiplication requires four real number multiplications. The computational complexity of the first calculation comes from the product of the diagonal matrix P⁻¹ and the matrix Q and the vector , the results are and respectively. The computational complexity of the second calculation comes from the multiplication of the iterative matrices B and F, which involves real number multiplications. The computation complexity of the third part of comes from the calculation of initial solution. The computation complexity of the last part comes from the computation of channel gain, NPI variance and LLR. Therefore, the total number of multiplications required by the WeJi is . Figure 2.7 shows the comparison of the numbers of the real number multiplications between the WeJi and other methods. The WeJi has lower computational complexity compared with that of GAS, SOR, and SSOR methods. When K = 2, the computational complexity of NSA is relatively low. In general, K should not be less than 3 in NSA to ensure the accuracy of detection. When K = 3, NSA shows higher computation. As a result, the reduction in computation of NSA is negligible.

Fig. 2.7

Comparison of the numbers of the actual multiplications between the WeJi and other algorithms.

© [2018] IEEE. Reprinted, with permission, from Ref. [29]

On the other hand, we need to consider the hardware implementation of the WeJi so that it can be executed as parallel as possible. The solution of the WeJi in Eq. (2.90) can be written as

(2.113)

The computation of only requires the elements of the previous iterations, so all the elements of can be computed in parallel. However, in the GAS, the successive over-relaxation (SOR) method [26] and the symmetric successive over-relaxation (SSOR) method [27], each transmitted signal has a strong correlation in the iterative steps. When calculating , we need of the Kth iteration and of the (K-1)th iteration. This means that the computation for each element cannot be executed in parallel. Therefore, neither the GAS [25] nor the SOR [26] can achieve a high data throughput rate and their throughput rates are far lower than that of the WeJi.

It was noted that the detection method is also proposed based on the Jacobi iteration in the literature [28]. Compared with this method, the WeJi described in this section achieves better performance in the following three aspects. First, the WeJi is a method based on hardware architecture design consideration, that is, the hardware implementation is fully considered in the process of algorithm optimization and improvement. In the process of algorithm design for the WeJi, the detection accuracy, computational complexity, parallelism, and hardware reusability are considered. On the contrary, hardware implementation problems are not considered in the literature [28], such as parallelism and hardware reuse. Second, the initial iterative solutions in WeJi in this section are different from that in literature [28]. The initial value in literature [28] is fixed. By contrast, the method described in this section takes into account the characteristics of the massive MIMO system, including a computational method for the initial solution. According to Eq. (2.91), the initial iterative solution is close to the final result, so the number of iterations can be reduced and the hardware consumption can be minimized. Furthermore, the method of the initial solution for the WeJi is similar to the algorithm in the later iterative steps, and the hardware resources can be reused. Because the Gram matrix G computation will occupy a large number of clock cycles before the iteration, the reuse of hardware resources will not affect the system throughput rate. Third, compared with the literature [28], the WeJi introduces a weighed factor, as shown in Eq. (2.90), so that improves the accuracy of the solution and consequently reduces hardware resource consumption. In addition, the same unit can be reused to increase unit utilization during the pre-iteration and the iteration, and this reuse does not affect the data throughput of hardware.

Next, we will discuss the frame error rate (FER) of the WeJi and the latest other signal detection algorithms. The FER performance of the exact matrix inversion (the CHD) algorithm is also used as a comparison object. In comparison, we consider modulation scheme of 64-QAM. The channel is assumed to be an independent and identically distributed Rayleigh fading matrix. The output (LLR) is adopted by Viterbi decoding. In the receiver, the LLR is the soft-input of the viterbi decoding. As for 4G and 5G, we have discussed a kind of parallel cascade convolution code, the Turbo code [1]. The Turbo scheme currently used in 4G is also an important encoding scheme for 5G and is widely used. Furthermore, these emulation settings are often used in many massive MIMO detection algorithms and architectures in the 5G communications.

Figure 2.8 shows the FER performance curves for the WeJi [29], NSA [7], RI [14], intra-iterative interference cancelation (IIC) [30], CG [18], GAS [12, 25], OCD [31] and MMSE [8, 24]. In Fig. 2.8a, the algorithm in the 64 * 8 massive MIMO system with 1/2 bit rate is simulated. Figure 2.8b shows the FER performance of the 128 * 8 massive MIMO system with 1/2 bit rate. To demonstrate the advantages of the proposed method at higher bit rates, Fig. 2.8c shows the performance of a 128 *8 massive MIMO system with a 3/4 bit rate. These simulation results show that the WeJi can achieve near-optimal performance at different MIMO scales and bit rates. To achieve the same FER, the SNR required by the WeJi is almost the same as that required by the MMSE, but lower than that required by the OCD, CG, GAS, IIC, RI, and NSA. By Fig. 2.8, the proposed WeJi can achieve better FER performance than the existing technical methods in different MIMO scales.

Fig. 2.8

Performance diagram of various algorithms. a N_r = 64, N_t = 8, 1/2 bit rate, b N_r = 128, N_t = 8, 1/2 bit rate, c N_r = 128, N_t = 8, 3/4 bit rate.

© [2018] IEEE. Reprinted, with permission, from Ref. [29]

2.4.3 Analyses on Channel Model Effects

The previous simulation results are obtained based on the Rayleigh fading channel model. In order to prove the superiority of the proposed algorithm in a more real channel model, Fig. 2.9 shows the effects of large-scale fading and spatial correlation of MIMO channels on the FER performance of different algorithms. The Kronecker channel model [20] was used to evaluate the FER performance of algorithm, because it is more practical than the independent and equally distributed Rayleigh fading channel model. The Kronecker channel model assumes that transmission and reception are separable, and the measurements show that the Kronecker model is a good approximation to the nonline-of-sight scenario. Therefore, this model is widely used in the literature. In the channel model, the elements of the channel matrix satisfy , where is an arbitrary function that interprets channel attenuation such as shadow and path loss. Consider the channel attenuation variance , where , , κ and denote the user’s location, base station’s location, path loss index and Euclidean norm respectively. The independent shadow fading C satisfies . Combined with the correlation matrix (R_r), the Kronecker channel matrix H can be written as

Fig. 2.9

FER performance of Kronecker channel model.

© [2018] IEEE. Reprinted, with permission, from Ref. [29]

(2.114)

is a random matrix whose elements are independent and identically distributed. It is a complex Gaussian distribution with the zero mean and unit variance. Exponential correlation is a model used to generate correlation matrices. The elements in correlation matrix R_r can be written as

(2.115)

where ξ is the correlation factor between adjacent branches. The users in the same cell are evenly distributed in the hexagon with radius r = 500 m. Now we assume κ = 3.7, , the transmitting power for the simulation. The correlation factors are 0.2, 0.5, and 0.7 respectively. Figure 2.9 shows that in order to achieve the same FER, the SNR required by the proposed algorithm is also smaller than CG, RI and NSA, which proves that the algorithm can maintain its advantages in a real model.

2.5 Conjugate Gradient Algorithm

2.5.1 Algorithm Design

CG is an iterative algorithm for solving linear matrix equations [9]. Its approximate solution is presented as follows:

(2.116)

where p is an auxiliary vector, K is the number of iterations, and can be calculated in Eq. (2.117)

(2.117)

Here r is the residual vector, which is expressed by Eq. (2.118)

(2.118)

By using the CG algorithm, we can solve the linear matrix equation. Hence this method can be applied in the massive MIMO signal detection. However, the conventional CG algorithm still has some deficiencies. First of all, in the original algorithm, although every element in the vector can be multiplied and accumulated in parallel with each row element of the matrix when matrix times vector, there is a strong data dependence between the steps, so the computation must be carried out step by step according to the order of the algorithm, and the degree of parallelism is not high. Secondly, the conventional CG algorithm does not provide information about the initial value of iteration, but the zero vector is usually used in the iteration. Obviously, zero vector can only satisfy the feasible requirements, but it is not optimal. Therefore, finding a better initial iterative value can make the algorithm converge faster and reduce the number of iterations, thus indirectly reducing the computational complexity of the algorithm. The CG algorithm has been improved to satisfy better parallelism and convergence for the above two points. The improved CG algorithm is named as three-term-recursion conjugate gradient (TCG) algorithm [32].

The conventional CG algorithm is equivalent to Lanczos orthogonalization algorithm [9], so it is presented as

(2.119)

Here, polynomial q_j can be used to make . When , , , we can get and Eq. (2.120)

(2.120)

Since , and are orthogonal, i.e., , We can derive from Eq. (2.120)

(2.121)

(2.122)

And since , and , we can get

(2.123)

We can also derive from Eq. (2.120)

(2.124)

By processing the massive MIMO system and utilizing the above derivation process, we can get the TCG algorithm that is applied in the massive MIMO signal detection.

Algorithm 2.2

The parallelizable Chebyshev iteration algorithm used in the massive MIMO to detect the minimum mean square error method

We can see from the Algorithm 2.2 that no data dependence exists between and and between and . So they can perform computation at the same time, increasing the parallelism of the algorithm. Besides since there is also operation of matrix multiplied by the vector, each element in the vector can still be multiplied and accumulated with each row element of the matrix in the hardware design.

2.5.2 Convergence

After discussing the CG and the optimization of CG, it is necessary to study the convergence of CG. First, the polynomial is defined, shown as the following (2.125):

(2.125)

(2.126)

When the above constraint conditions are extended to the case of , the polynomial (2.125) becomes

(2.127)

The polynomial can also be expressed as

(2.128)

If is defined, it can be derived from Eq. (2.127)

(2.129)

We can find that satisfies from the above expression (2.129).

(2.130)

where . If is for an exact solution vector, the expression of CG convergence is obtained from (2.131):

(2.131)

The observation (2.131) shows that denotes the exact solution vector, which is a definite invariant vector. So when the initial value of iteration is determined, on the right side of (2.131) can be regarded as a constant. When the number of iterations K increases, exponentially approaches 0. This indicates that CG can converge faster when K is large, and the error of CG algorithm decreases with the increase of number of iterations K.

2.5.3 Initial Iteration Value and Search

Iteration requires corresponding initial value vectors in the TCG. As mentioned earlier, how to find a good initial value vector is crucial to the convergence rate of the algorithm. When introducing PCI in Sect. 2.3, we mentioned an eigenvalue-based initial value algorithm, which can be used in the TCG to make algorithm converge faster. Besides this, another quadrant-based initial value algorithm will be discussed in this section.

As discussed above, if the influence of noise n is ignored in the massive MIMO system, the real part and the imaginary part of the received signal y, the transmitted signal s and the channel matrix H are separated and written in the real number form, and the relationship between them can be expressed as

(2.132)

(2.133)

where the subscript R indicates the expression in the form of real numbers. By using real number form in ZF, the transmitted signal s can be estimated as

(2.134)

In the channel matrix H_R in the real number form, each element is independent and identically distributed, and obeys the standard Gaussian distribution, so the matrix can be approximated as a diagonal matrix. And the diagonal elements are all nonnegative since each diagonal element is the sum of squares. If is approximated as a diagonal matrix, then is also a diagonal matrix, and its diagonal elements are also nonnegative. Now consider the relationship between the ith element of vector and the jth element of vector . According to Eq. (2.134), we can get

(2.135)

where is nonnegative, so both and are positive or negative at the same time. By transforming the solution from real number form to complex form, we can deduce that and are in the same quadrant. Based on this conclusion, a new iteration initial value can be proposed.

Now we are considering the ith element of the transmitted vector in the massive MIMO system with modulation order of 64-QAM. Since and are in the same quadrant, assuming that is in the first quadrant, then we make . The coordinate of the hollow circle is (4,4) as in Fig. 2.10. Since will eventually be located in the first quadrant, and the average distance between the point (4, 4) and all constellation points in the first quadrant is less than the average distance between the point (0, 0) and all constellation points in the first quadrant, the point (4, 4) is closer to the final solution so that algorithm 2.3 can converge as soon as possible. The quadrant-based iterative initial value algorithm can be obtained based on this principle.

Fig. 2.10

Schematic diagram of quadrant—based initial value algorithm

Algorithm 2.3

Quadrant-based initial value algorithm

Taking the quadrant-based initial iteration value algorithm as the initial value result of TCG, the TCG can converge faster. According to the error requirements of the massive MIMO signal detection, the number of iterations can be reduced, and the computational complexity of the algorithm can be reduced from another angle, which makes the TCG even better.

In the massive MIMO signal detection algorithm, no matter whether it is a linear detection algorithm or a nonlinear detection algorithm, each algorithm will eventually obtain the transmitted vector through its computation and search for the constellation points with the smallest Euclidean distance from them according to each element in . These constellation points will serve as the final estimated transmitted vector . It is worth noting that at present, all modulation modes of the MIMO system are two-dimensional modulation, that is, the modulated constellation points are all located in the plane rectangular coordinate system. Therefore, when calculating the minimum Euclidean distance constellation point of each element in vector , it is actually calculating the constellation point with the minimum distance from the element plane in . Based on this, we can simplify the traditional method, in which the Euclidean distance between the middle element and each constellation point in the constellation diagram is first obtained, and then the size is compared. As shown in Fig. 2.11, it is assumed that an element is in the position of the hollow circle in the diagram. By dividing the constellation diagram into several parts with the dotted line in the diagram, the nearest constellation point to the hollow circle can be determined according to the region in which the hollow circle is located, that is, the constellation points in the region are just the constellation points with the smallest distance from the element. By using this analysis, the operation of finding the constellation points of the minimum Euclidean distance of the elements in the transmitted vector after the iteration can be simplified correspondingly, and a point finding algorithm based on rounding is obtained. Using algorithm 2.4 to locate constellation points, we can quickly find the final result at very low computational complexity.

Fig. 2.11

Schematic diagram of rounding off-based point seeking method

Algorithm 2.4

Rounding off-based point seeking method

2.5.4 Complexity and Parallelism

In the TCG algorithm, the first part is a series of initial value calculations, including , W and z₀. These calculations involve the matrix multiplied by the vector y, the matrix multiplied by the matrix H to get the matrix W, the matrix W multiplied by the vector . The computational complexities of them are , and respectively. The complexity of matrix W is only half (that is ) because it is a symmetric matrix. Since these parameters are computed only once in the pre-iteration process, the complexity is counted only once. The calculation of the second part is the multiplication of matrices and vectors in the iteration, that is, the vector z is multiplied by the matrix W. The computational complexity is . The calculation of the third part is the calculation of two inner products of and , and their complexity is and respectively. The last part is the update for and z whose complexities are and respectively. Note that when the number of iterations K is 1, the value of ρ is 1, so the step 13 and the step 14 in the algorithm do not need to be calculated. Therefore the total complexity of the algorithm is . Table 2.2 lists the complexity of various algorithms at different iterations. We can see from the table that the complexity of the TCG algorithm is low. For example, when , , K = 2, the OCD algorithm has the lowest complexity of 16,416 real number multiplications, and the TCG algorithm has 21,632 real number multiplications, which is higher than that of the OCD algorithm, while the complexity of other algorithms is higher than that of the TCG algorithm.

Table 2.2

Complexity comparison

	K = 2	K = 3	K = 4	K = 5
NSA [7]	2NrNt2 + 6NtNr + 4Nt2 + 2Nt	2NrNt2 + 6NtNr + 2Nt3 + 4Nt2	2NrNt2 + 6NtNr + 6Nt3	2NrNt2 + 6NtNr + 10Nt3-6Nt2
INS [17]	2NrNt2 + 6NtNr + 10Nt2 + 2Nt	2NrNt2 + 6NtNr + 14Nt2 + 2Nt	2NrNt2 + 6NtNr + 18Nt2 + 2Nt	2NrNt2 + 6NtNr + 22Nt2 + 2Nt
GAS [12]	2NrNt2 + 6NtNr + 10Nt2 − 2Nt	2NrNt2 + 6NtNr + 14Nt2 − 6Nt	2NrNt2 + 6NtNr + 18Nt2 − 10Nt	2NrNt2 + 6NtNr + 22Nt2 − 14Nt
CG [16]	2NrNt2 + 6NtNr + 8Nt2 + 33Nt	2NrNt2 + 6NtNr + 12Nt2 + 49Nt	2NrNt2 + 6NtNr + 12Nt2 + 49Nt	2NrNt2 + 6NtNr + 20Nt 2 + 81Nt
CGLS [18]	24NtNr + 20Nr2 + 8Nr + 44Nt	32NtNr + 28Nr2 + 12Nr + 66Nt	40NtNr + 36Nr2 +16Nr + 88Nt	48NtNr + 44Nr2 + 20Nr + 110Nt
OCD [19]	16NtNr + 4Nt	24NtNr + 6Nt	32NtNr + 8Nt	40NtNr + 10Nt
TCG	2NrNt2 + 4NtNr + 12Nt2 + 48Nt	2NrNt2 + 4NtNr + 16Nt2 + 80Nt	2NrNt2 + 4NtNr +20Nt2 + 112Nt	2NrNt2 + 4NtNr + 24Nt2 + 144Nt

The parallelism of the TCG algorithm is also very important. According to Step 2, Step 3, Step 6, and Step 7 of the algorithm, each row element of the matrix is multiplied by the vector when the matrix is multiplied by the vector, and each multiply accumulate operation can be carried out simultaneously. There is no data dependence between the eighth and ninth steps of the algorithm and between the fourteenth and fifteenth steps of the algorithm except the parallel computation of the matrix multiplied by vectors, so the computation can be performed simultaneously. This algorithm has parallelism in two cases. Compared with other algorithms, the parallelism between steps has great advantages.

As previously analyzed, the TCG has significant advantages over other algorithms in complexity and parallelism. In addition, the TCG algorithm optimizes the initial value of iteration and the final point finding method, so that this algorithm can bring the maximum benefit. In general, complexity and accuracy are in contradiction, and lower complexity often means lower precision. Therefore, although the above performance shows the excellent performance of the TCG algorithm in this respect, it does not show that it is a relatively comprehensive algorithm, and it is also required to consider the accuracy of the algorithm.

2.5.5 Symbol Error Rate

Figure 2.12 shows the effect of the initial value on the performance of the CG algorithm in the massive MIMO system with the size of 128 × 8. The initial value is calculated by the quadrant division initial value algorithm. Obviously, for CG algorithm with the same SNR, the initial value algorithm can achieve lower symbol error rate (SER) for the same number of iterations, and when the number of iterations K = 2, the initial value algorithm has a greater impact on the performance of the algorithm.

Fig. 2.12

Effect of the initial value algorithm on the SER of the CG algorithm

Figure 2.13 shows the SNR and SER curves of various algorithms for different number of iterations in the massive MIMO system with a size of 128*8 and with modulation order of 64-QAM. It can be clearly seen from the graph that, the SNR of NSA [7] is 10.17 dB, and the SNR of MMSE [8, 24], CG [12] and OCD [31] is around 9.66, 10.10, and 10.11 dB when the number of iterations K = 2 and SER is 10⁻². However the SNR of WeJi [29] is about 10.42 db, and the SNR of NSA [7] is larger than 20 dB. Figure 2.14 shows the influence of various algorithms on SER in MIMO systems of different sizes. We can see that the SER of TCG algorithm is obviously lower than that of the NSA and WeJi [29] in MIMO systems with different sizes.

Fig. 2.13

SER curves of various algorithms for different number of iterations in massive MIMO system with 128 * 8

Fig. 2.14

Influence of different algorithms on SER in MIMO systems with different scales

Therefore, the TCG has obvious advantages compared with the three algorithms of NSA, CG and CGLS. The TCG is better than the other three methods in complexity and parallelism, while the performance of SER is still worse than the three methods under the same SNR. Compared with the PCI algorithm, the TCG algorithm has lower complexity and better parallelism than the PCI. Besides having the parallelism of the matrix multiplied by vectors in the PCI, there is also parallelism between the steps in the algorithm. Compared with the OCD algorithm, the complexity of the TCG algorithm is higher than that of the OCD algorithm, and the OCD algorithm also has better performance of SER. However, the OCD method also has its inherent shortcomings: Because the OCD algorithm has too strong data dependence, the obtained data through calculation need to be stored in the register and each calculation needs to read data from the register, the parallelism is so poor that the operation can only be executed step by step sequentially.

Compared with the NSA, CG, CGLS, PCI, OCD, and other algorithms, the TCG algorithm is not the best in all performance parameters, but it overall achieves a better compromise. In the actual massive MIMO signal detection, the appropriate algorithm can be selected according to the actual application requirements.

References

1. 1.
   Andrews JG, Buzzi S, Wan C et al (2014) What Will 5G Be? IEEE J Sel Areas Commun 32(6):1065–1082
2. 2.
   Kim SP, Sanchez JC, Rao YN et al (2006) A comparison of optimal MIMO linear and nonlinear models for brain-machine interfaces. J Neural Eng 3(2):145–161
3. 3.
   Burg A, Borgmann M, Wenk M et al (2005) VLSI implementation of MIMO detection using the sphere decoding algorithm. IEEE J Solid-State Circuits 40(7):1566–1577
4. 4.
   Trimeche A, Boukid N, Sakly A et al (2012) Performance analysis of ZF and MMSE equalizers for MIMO systems. In: International conference on design & technology of integrated systems in nanoscale Era, pp 1–6
5. 5.
   Rusek F, Persson D, Lau BK et al (2012) Scaling up MIMO: opportunities and challenges with very large arrays. Sig Process Mag IEEE 30(1):40–60
6. 6.
   Teaching and Research Office of Computational Mathematics (2000) Fundamental of Numeric Analysis. Tongji University, Tongji University Press
7. 7.
   Wu M, Yin B, Wang G et al (2014) Large-scale MIMO Detection for 3GPP LTE: algorithms and FPGA implementations. IEEE J Selected Topics Signal Process 8(5):916–929
8. 8.
   Auras D, Leupers R, Ascheid GH (2014) A novel reduced-complexity soft-input soft-output MMSE MIMO detector: algorithm and efficient VLSI architecture. In: IEEE international conference on communications, pp 4722–4728
9. 9.
   Golub GH, Van Loan CF (1996) Matrix computations. Mathe Gazette 47(5 Series II):392–396
10. 10.
    Shahab MB, Wahla MA, Mushtaq MT (2015) Downlink resource scheduling technique for maximized throughput with improved fairness and reduced BLER in LTE. In: International conference on telecommunications and signal processing, pp 163–167
11. 11.
    Zhang C, Li Z, Shen L et al (2017) A low-complexity massive mimo precoding algorithm based on chebyshev iteration. IEEE Access 5(99):22545–22551
12. 12.
    Dai L, Gao X, Su X et al (2015) Low-complexity soft-output signal detection based on Gauss-Seidel method for uplink multiuser large-scale MIMO System. IEEE Trans Veh Technol 64(10):4839–4845
13. 13.
    Peng G, Liu L, Zhang P et al (2017) Low-computing-load, high-parallelism detection method based on Chebyshev Iteration for massive MIMO Systems with VLSI architecture. IEEE Trans Signal Process 65(14):3775–3788MathSciNet
14. 14.
    Gao X, Dai L, Ma Y et al (2015) Low-complexity near-optimal signal detection for uplink large-scale MIMO systems. Electron Lett 50(18):1326–1328
15. 15.
    Gutknecht MH, Röllin S (2000) The Chebyshev iteration revisited. Parallel Comput 28(2):263–283MathSciNetzbMATH
16. 16.
    Yin B, Wu M, Cavallaro JR et al (2014) Conjugate gradient-based soft-output detection and precoding in massive MIMO systems. In: Global communications conference, pp 3696–3701
17. 17.
    Čirkić M, Larsson EG (2014) On the complexity of very large multi-user MIMO Detection
18. 18.
    Yin B, Wu M, Cavallaro JR et al (2015) VLSI design of large-scale soft-output MIMO detection using conjugate gradients. In: IEEE international symposium on circuits and systems, pp 1498–1501
19. 19.
    Wu M, Dick C, Cavallaro JR et al (2016) FPGA design of a coordinate descent data detector for large-scale MU-MIMO. In: IEEE international symposium on circuits and systems, pp 1894–1897
20. 20.
    Werner K, Jansson M (2009) Estimating MIMO channel covariances from training data under the Kronecker model. Sig Process 89(1):1–13zbMATH
21. 21.
    Sun Q, Cox DC, Huang HC et al (2002) Estimation of continuous flat fading MIMO channels. IEEE Trans Wireless Commun 1(4):549–553
22. 22.
    Rappaport TS (2002) Wireless Communications—Principles and Practice. Second Edition. (The Book End). 8(1):33–38
23. 23.
    Li X, Chen CP (2013) Inequalities for the gamma function. J Inequal Pure Appl Mathematics 8(1):554–563MathSciNet
24. 24.
    Prabhu H, Rodrigues J, Liu L et al (2017) A 60 pJ/b 300 Mb/s 128 × 8 Massive MIMO Precoder-Detector in 28 nm FD-SOI
25. 25.
    Wu Z, Zhang C, Xue Y et al (2016) Efficient architecture for soft-output massive MIMO detection with Gauss-Seidel method. In: IEEE international symposium on circuits and systems, pp 1886–1889
26. 26.
    Zhang P, Liu L, Peng G et al (2016) Large-scale MIMO detection design and FPGA implementations using SOR method. In: IEEE international conference on communication software and networks, pp 206–210
27. 27.
    Quan H, Ciocan S, Qian W et al (2015) Low-complexity MMSE signal detection based on WSSOR method for massive MIMO Systems. In: IEEE international symposium on broadband multimedia systems and broadcasting, pp 193–202
28. 28.
    Kong BY, Park IC (2016) Low-complexity symbol detection for massive MIMO uplink based on Jacobi method. In: IEEE international symposium on personal, indoor, and mobile radio communications, pp 1–5
29. 29.
    Peng G, Liu L, Zhou S, et al (2017) A 1.58 Gbps/W 0.40 Gbps/mm² ASIC implementation of MMSE detection for $128x8$ 64-QAM Massive MIMO in 65 nm CMOS. IEEE Trans Circuits Syst I Regular Papers PP(99):1–14
30. 30.
    Chen J, Zhang Z, Lu H et al (2016) An intra-iterative interference cancellation detector for large-scale MIMO communications based on convex optimization. IEEE Trans Circuits Syst I Regul Pap 63(11):2062–2072
31. 31.
    Wu M, Dick C, Cavallaro JR et al (2016) High throughput data detection for massive MU-MIMO-OFDM using coordinate descent. IEEE Trans Circuits Syst I Regul Pap 63(12):2357–2367
32. 32.
    Tongxiang Gu (2015) Iterative methods and pretreatment technology. Science Press, Beijing