# Improved Approximate Expectation Propagation Massive MIMO Detector with Second-Order Richardson Iteration

Authors: Qian Deng; Xuehui Chen; Wanting Fu; Xiaopeng Liang; Shilong Xie; Feng Shu; Yuan Yuan Wu
Journal: Wireless Communications and Mobile Computing
 (2022)

---

## Abstract

The expectation propagation (EP) detector achieves significantly better performance than the linear detectors (such as minimum mean squared error detector) in massive MIMO systems, which has drawn great attention recently. EP’s approximation (EPA) algorithm simplifies the update formula of the EP algorithm by reexpressing the moment matching condition so that the number of matrix inversions in the EP algorithm is reduced to one. However, the expense is that the EPA algorithm requires higher accuracy for this inversion; otherwise, the bit-error-rate (BER) performance will suffer serious losses. To tackle this issue, the SORI iterative algorithm is introduced to obtain the high-precision result of this inversion to ensure the good BER performance of the EPA algorithm. First, the new expression of the SORI iterative algorithm is derived under the equivalent real-valued system. Second, the improved EPA-SORI algorithm is then introduced by the SORI algorithm, which is used to approximate the initial value of the EPA algorithm under the real-valued system. Finally, by designing the initial solution and the relaxation factor of the EPA-SORI iterative algorithm, the convergence rate can be quickly increased without increasing the complexity. Simulation and complexity results exhibit that in various massive MIMO system configurations, the proposed EPA-SORI algorithm can achieve the same BER performance as the Exact EP algorithm with significantly lower complexity. At the same time, compared with MMSE and the existing EPA algorithms, the proposed EPA-SORI algorithm has a better performance-complexity trade-off advantage, which is more obvious in scenarios with high modulation order and a large number of users.


---

## Body

## 1. Introduction

Future wireless communication requires a high data rate and a tremendous amount of connection for emerging applications such as the Internet of things (IoT) [1, 2]. The access of massive IoT devices to the network will lead to tremendous growth in the data volume of mobile communication services, and the wireless network capacity will face unprecedented challenges [3–5]. Aiming at this challenge, main solutions include the usage of larger bandwidth, higher-order MIMO [6, 7], higher-order modulation, more effective coding, and so on. Among these solutions, through deep utilization of spatial dimensions, massive MIMO technology attains enhanced wireless communication capacity and spectral efficiency and has become a key technology for 5 G/B5G wireless communication [8]. With enormous system dimensions and the use of higher-order modulation, signal detection faces a challenge in terms of computational burden and hardware implementation [9].The traditional optimal signal detector, maximum likelihood (ML) detector, faces the problem of exponential increase in computational complexity for massive MIMO systems [10]. In contrast, linear detection algorithms (such as MMSE, ZF) have reduced computational complexity. Especially when the loading factor ζ≪1 (ζ≜Nt/Nr, where Nt and Nr represent the number of single antenna users and the number of base station antennas, respectively), the minimum mean squared error (MMSE), and Zero Forcing (ZF) linear detection algorithms can achieve near-optimal system performance. However, as a large number of IoT devices are connected to the cellular network, more users need to be served in a cell (e.g., mobile phones, unmanned aerial vehicles (UAVs) [11, 12], sensors [13], vehicle to vehicle (V2V) [14]). Unfortunately, as the number of users increases, the performance of the linear detection algorithm suffers severe degradation. Compared with ML algorithm, other complex detectors (e.g., belief propagation (BP) [15, 16], approximate message passing (AMP) [17]) can achieve excellent performance. However, the convergence speed of the iterative update in the BP algorithm will decrease, and the update formula during the iteration process becomes more complicated, because massive MIMO contains a large number of ring structures when the number of users increases. To address the above issues, an EP algorithm is proposed [18, 19]. For a random value of ζ and high-order modulation, the EP algorithm not only shows significantly better performance than other algorithms such as BP, AMP and MMSE, but also has great flexibility and strong robustness, which has attracted wide attention. However, during each iteration of EP, it is necessary to perform a full matrix inversion with a complexity as high as O (N3). In addition, the huge computational cost makes it difficult to implement on hardware. To address this problem, EP’s approximation algorithm (EPA) simplifies the update formula of the EP algorithm by reexpressing the moment matching condition, so that the number of matrix inversions in the EP algorithm is reduced to one [20]. However, the expense is that the EPA algorithm requires higher accuracy for this inversion; otherwise, the bit-error-rate (BER) performance will suffer serious performance losses. Recently, some methods based on matrix polynomial decomposition are proposed to apply to the EPA algorithm (such as EPA-NSA and EPA-wNSA) [9, 20]. When loading factor ζ≪1, the EPA-NSA algorithm shows good performance. But with the increase of ζ, the EPA-NSA algorithm shows slow convergence or even nonconvergence, resulting in serious degradation of system performance. Although EPA-wNSA algorithm has an improved performance compared with the EPA-NSA algorithm, the degree of improvement is very limited. In addition, the above-mentioned algorithms like EPA-NSA, EPA-wNSA, etc., in spite of avoiding the direct inversion of the matrix and reducing some complexity, calculation of the Gram matrix with a complexity up to O (Nt2Nr) still needs to be involved to obtain high-precision signal detection. Some famous linear iterative method, such as Gauss-Seidel [21–23], successive over relaxation [24, 25], suffer from calculating the Gram matrix and low parallelism. Therefore, they can not directly be applied to alleviate the computational burden of EP iteration. In order to solve these complex issues, an improved algorithm EPA-SORI is proposed in this paper. The EPA-SORI algorithm introduces the SORI iterative algorithm in EPA to obtain a high-precision result of this one inversion to ensure that the EPA algorithm has good error bit rate performance. The contributions of this paper are as follows:(1)
We first deduce a new expression of the SORI iterative algorithm in the real value system. Then, the SORI algorithm under the real value system is further applied to the EPA algorithm. Furthermore, LLR approximation is provided to enhance the accuracy of the EPA-SORI detector.(2)
According to the random matrix theory, under the real-valued system, promising initial solution and relaxation factors are used to further enhance the convergence rate and accuracy and then reduce the computational complexity. Furthermore, a theoretical analysis of the convergence speed of the proposed EPA-SORI algorithm is presented. Theoretical analysis proves that the convergence speed of the proposed EPA-SORI algorithm is significantly higher than the recently reported EPA-wNSA algorithm.(3)
In the iterative process, the proposed EPA-SORI algorithm not only requires no matrix inversion operations but also avoids direct calculation of the Gram matrix, effectively reducing the complexity of the entire algorithm. At the same time, the proposed EPA-SORI algorithm is high-parallel and hardware-friendly.(4)
Simulation and complexity results show that the bit error rate performance of the EPA-SORI algorithm is much better than MMSE, and the complexity is much lower than MMSE. Compared with existing EPA algorithms (such as EPA-INSA, EPA-wNSA, etc.), the proposed EPA-SORI can achieve performance close to Exact EP with lower complexity and higher convergence rate. Furthermore, with high modulation order and a large number of users, EPA-SORI will show a more obvious performance-complexity trade-off advantage.


### 1.1. Notation

Matrices and column vectors are represented by uppercase and lowercase boldface letters, respectively. The element ini-th row and j-th column of matrix A is denoted by Ai,j. ⋅T, ⋅H, ⋅−1, ⋅2, and ⋅ denote transpose, conjugate transpose, inversion, 2−norm and determinant, respectively. Also, the probability distribution of s is denoted by ps. The real part and imaginary part are denoted by Re⋅ and Im⋅ , respectively. Ny:μ,Σ represents the Gaussian probability distribution with mean μ and variance Σ.


## 1.1. Notation

Matrices and column vectors are represented by uppercase and lowercase boldface letters, respectively. The element ini-th row and j-th column of matrix A is denoted by Ai,j. ⋅T, ⋅H, ⋅−1, ⋅2, and ⋅ denote transpose, conjugate transpose, inversion, 2−norm and determinant, respectively. Also, the probability distribution of s is denoted by ps. The real part and imaginary part are denoted by Re⋅ and Im⋅ , respectively. Ny:μ,Σ represents the Gaussian probability distribution with mean μ and variance Σ.


## 2. System Model

We consider a massive MU-MIMO system in whichNr antennas are deployed at the base station to serve Nt users. s˜=s˜1,s˜2,…,s˜NtT is a Nt×1-dimensional transmission signal vector, where s˜i∈Ω˜ and Ω˜ are the symbol set of the M− QAM constellation (e.g., M=16/64/256, etc.). We assume that H˜∈ℂNr×Nt is a flat Rayleigh fading channel. The received signal y˜∈ℂNr×1 at the receiver can be modeled by the following:(1)y˜=H˜s˜+n˜,where n˜ is the additive white Gaussian noise (AWGN), and its elements satisfy a complex Gaussian distribution with mean 0 and variance σ˜n2. Here, the received signal can be reexpressed by the equivalent real-valued augmented vectors as follows:(2)Rey˜Imy˜=ReH˜−ImH˜ImH˜ReH˜Res˜Ims˜+Ren˜Imn˜,where y=Rey˜Imy˜T∈ℝ2Nr×1, s=Res˜Ims˜T∈Ω2Nt×1 and Ω represents the set of real and imaginary parts of the point set on the constellation in M− QAM modulation. n=Ren˜Imn˜T, where n∈N0,σn2I2Nr satisfy a real Gaussian distribution with mean 0 and variance σn2=σ˜n2/2. Thus, (1) can be modeled by y=Hs+n, where H=ReH˜−ImH˜ImH˜ReH˜.


## 3. EP and EPA Algorithms

EP algorithm is a reasoning method based on Bayesian inference, which is used to estimate the value ofx under the condition that y is known when a joint distribution pxy is given [26–28]. If we use pxy to estimate x directly, the complexity will increase exponentially with the dimensions of x and y, which consume huge hardware resources in massive MIMO scenario. Therefore, an approximate distribution is used to estimate pxy. In model y=Hs+n, the joint posterior distribution of the y and s is ps|y∝Ny:Hs,σn2I2Nr⋅ps. To facilitate the analysis, according to [7], when the transmitted symbols are independent of each other, the EP detector uses a nonstandard Gaussian distribution to replace the prior distribution ps of each transmitting antenna, that is, p^si=exp−1/2Λisi2+γisi. And then, a posterior distribution p^s|y whose distribution satisfies the exponential family approximate. p^s|y can be expressed as follows:(3)p^s|y=Ny:Hs,σn2I2Nr∏i=12Ntp^si=Ny:Hs,σn2I2Nr∏i=12Ntexp−12Λisi2+γisi.The cumulative multiplication in (3) can be transformed into the following form:(4)p^s|y=Ny:Hs,σn2I2Nrexp−12sTΛs+γTs,where γ=γ1,γ2…,γ2NtT∈ℝ2Nt×1 is a real-valued column vector, and Λ=diagΛ1,Λ2…,Λ2Nt is a diagonal matrix. According to (4), the expression of mean μ and variance Σ of p^s|y is as follows:(5)Σ=1σn2HTH+Λ−1,(6)μ=Σ1σn2HTy+γ.Having constructedp^s|y, it is necessary to use the moment matching technique so that p^s|y and ps|y are as similar as possible. However, the complexity of p^s|y is very high, and it is difficult to directly obtain Λ and γ. Aiming at this problem, a sequential EP algorithm is proposed in [9, 20]. It is assumed that p^si|y is the i− th edge of p^s|y:(7)p^\isi|y=p^si|yexp−0.5Λisi2+γisi.γi0=0 and Λi0=Es are adopted as the initial solutions of (7), where Es denotes the mean symbol energy. Then, the alternative distribution expression is expressed as follows:(8)prlsi|y=p^l/isi|y⋅psi,where p^l/isi is expressed as follows:(9)p^l/isi∝exp−0.5Visi2+ρisi.Then, the first-order momentsηil and the second-order moments χil of prlsi|y can be obtained by the following:(10)ηil=Eprlsi|ysi=∑si∈Ωsiexp−0.5Vil−1si2+ρil−1si∑u∈Ωexp−0.5Vil−1u2+ρil−1u,(11)χil=Eprlsi|ysi2−Eprlsi|y2si=∑si∈Ωsiexp−0.5Vil−1si2+ρil−1si∑u∈Ωexp−0.5Vil−1u2+ρil−1u−∑si∈Ωsiexp−0.5Vil−1si2+ρil−1si∑u∈Ωexp−0.5Vil−1u2+ρil−1u2.Assuming thattil, hi2l are the mean and variance of ql/isi, the values of Vil and ρil for l iterations can be obtained as follows:(12)Vil=1hi2l=1−∑iilEs∑iil,(13)ρil=tilhi2l=μil∑iil.Next, the updatedΛ and γ can be obtained by matching the first-order and second-order moments of prlsi|y and the cavity marginal distribution p^\isi|y. Since the calculation process of the EP algorithm is too cumbersome, and each iteration involves the matrix inversion calculation in equation (5), which brings difficulties to the signal detection in massive MIMO scenarios. To simplify the calculation process and reduce computational complexity, EP’s approximate algorithm EPA is proposed in [20]. The EPA algorithm reexpresses the moment matching conditions and no longer uses explicit matrix matching to calculate the results of each iteration. It simplifies the update formula in the iteration process of the EP algorithm and uses a fixed matrix to estimate the matrix V before the iteration. Then, the approximate value is used to estimate the value of Σ, and finally, we will obtain an approximate algorithm EPA which is different from the Exact EP iteration process. The EPA algorithm reduces the EP algorithm that requires multiple matrix inversion calculations to one time, which can effectively avoid the impact of multiple inaccurate approximate matrix inversion results. However, the expense is that the EPA algorithm requires a complete and high-precision matrix inversion to obtain an initial value of iteration. The implementation of EPA is detailed in Algorithm 1.Algorithm 1: The EPA algorithm.
(1) Input:y, H, L, σn2, Es, ϑ.(2) Output:s^=tL.(3)yMF=HTy.(4)W=HTH+σn2/EsI2Nt.(5)s=W−1y.(6)Di=hiThi,i=1,2,…,2Nt.(7)D=diagD1,D2,…,D2Nt.(8) fori=1,2,...,2Nt do(9)ti0=si/1−DiEs;(10) end for(11) repeat(12) fori=1,2,…,2Nt do(13)ηi=argminui∈Ωui−til−12;(14) end for(15)hl=σn−2Hη.(16)m=b−HThl.(17)tl=D−1Tm+ηl.(18)tl=1−ϑtl−1+ϑtl, the damping factor ϑ∈0,1.(19)l=l+1;(20) until convergence orl>L.


## 4. The Proposed EPA-SORI Detector

MMSE solution is usually used as the iterative initial solution of EP detection [20], then we also use MMSE solution as the iterative initial solution t0 of the proposed EPA-SORI algorithm. Define W≜HTH+σn2/EsI2Nt, and the initial iteration t0 can be obtained by the following:(14)t0=Σ0μ0=W−1yMF.From (14), the initial iteration solution problem can be treated as a solution of the equation: Wt0=yMF. Note that the result of this equation is required to be highly accurate; otherwise, the bit error rate performance will suffer serious degradation, so we propose to apply SORI iteration to obtain a high-precision iterative initial solution. In the complex value system, the corresponding SORI algorithm can be constructed as follows:(15)s¯l+1=s¯l−1+α1β1y˜MF−W˜s¯l+α1s¯l−s¯l−1,where l≥1, s¯ denotes the estimated value of W˜−1y˜MF, and α1 is the relaxation factor, β1=2/λ˜max+λ˜min, λ˜max and λ˜min are the maximum and minimum eigenvalues of W˜, respectively, where W˜=H˜HH˜+σ˜n2/EsINt. Note that the EPA algorithm is implemented in a real-valued system; thus, we need to apply the SORI algorithm to a real-valued system. Here, the SORI method can be reexpressed as follows:(16)s^l+1=s^l−1+α2β2yMF−Ws^l+α2s^l−s^l−1,where l≥1, s^∈ℝ2Nr×1 denotes the estimated value of t0, and α2 denotes the relaxation factor, β2=2/λ˜max+λ˜min, and λmax and λmin are the maximum eigenvalue and minimum eigenvalue of W, respectively. The relationship between yMF and y˜MF satisfies: yMF=Rey˜MFImy˜MFT. In recently reported literature [29], the SORI algorithm is mainly applied to signal detection under complex-valued systems and cannot be directly applied to real-valued systems. Thus, it is necessary to further deduce the maximum and minimum eigenvalues of W. As long as the eigenvalues of HTH can be obtained, the eigenvalues of W can be obtained.Lemma 1.
In the real-valued system, the maximum and minimum eigenvalues ofW are as follows:(17)λmax⟶Nr1+SizeCHSizeRH2λmin⟶Nr1−SizeCHSizeRH2,where SizeC⋅ and SizeR⋅ denote the number of columns and rows.Proof.
Assuming that one of the eigenvalues ofHTH is λ1, and one of the eigenvalues of H˜HH˜ is λ. Based on the properties of eigenvalues, we have the following:(18)λINt−H˜HH˜=0,(19)λ1I2Nt−HTH=0.
Here,H˜HH˜ and HTH can be expressed as follows:(20)H˜HH˜=Re2H˜+Im2H˜,(21)HTH=Re2H˜+Im2H˜00Re2H˜+Im2H˜.
Substitute (20) and (21) into (18) and (19), and we have the following:(22)λINt−Re2H˜+Im2H˜=0,(23)λ1INt−H˜HH˜00λ1INt−H˜HH˜=0.Simplify (23), and we have the following:(24)λ1INt−Re2H˜+Im2H˜2=0.The comparison between (24) and (22) shows that the eigenvalues of the two are the same, i.e., λmax⟶Nr1+Nt/Nr2,λmin⟶Nr1−Nt/Nr2. In addition, we can get the same conclusion by using random matrix theory [30, 31]. According to random matrix theory [30], in a real-valued system, Hi,j∼N0,1/2. Thus, 1/NrHi,j∼N0,1/2Nr. Hence, we have λmaxHTH=Nr1+2Nt/2Nr2 and λminHTH=Nr1−2Nt/2Nr2. Thus, Lemma 1 is proved.
The next step is to derive the relaxation factorα2 in the real-valued system. According to [29], α2≜2/1+1−ρ2Gri, where ρGri denotes the spectral radius of Gri, and Gri=I2Nt−β2W. Given that b and a are the maximum and minimum eigenvalues of Gri, ρGri=b, thus α2 can be obtained by the following formula:(25)α2=21+1−b2.
Since(26)a=−β2λmax+1=λmin−λmaxλmin+λmax,(27)b=−β2λmin+1=λmax−λminλmin+λmax.
Thus, we have the following:(28)α2=21+1−λmax−λmin/λmin+λmax2.
It can be seen thatα2 is only related to λmin and λmin.
According to the typical properties of the massive MIMO channel,H˜ satisfies asymptotic orthogonality, and the Gram matrix W˜ is dominant diagonally. From (20) and (21), we have W=W˜00W˜, thus, the matrix W is also dominant diagonally. Therefore, W could satisfy the approximation as follows:(29)Wi,j≈λmax+λmin2,i=j,0,i≠j.
It can be seen from (29), when β2=2/λmax+λmin, β2yMF can be used to approximate W−1yMF, the SORI algorithm can choose the initial iteration solution of s^0=0, s^1=β2yMF. The SORI algorithm is summarized in Algorithm 2.
Algorithm 2: The SORI algorithm.
(1) Input:y, H, K, σn2,Es.(2) Output:s^K+1.(3)yMF=HTy(4)λmax⟶Nr1+Nt/Nr2,λmin⟶Nr1−Nt/Nr2.(5)β2=2/λmax+λmin.(6)α2=2/1+1−λmax−λmin/λmin+λmax2.(7)s^0=0, s^1=β2yMF.(8) forl=1,2,…,K do(9)hk=Hs^l;(10)s^l+1=s^l−1+α2β2yMF−HThk−σn2Es−1s^l+α2s^l−s^l−1.(11) end for
In the proposed EPA-SORI algorithm,s^0=0 and s^1=D−1yMF are used as the initial iteration solutions to further improve the convergence speed. Since D is the diagonal matrix of the W matrix, the complexity will not increase. The EPA-SORI algorithm is summarized in Algorithm 3, and an intuitive diagram of the processing of the EPA-SORI algorithm is presented in Figure 1.
Algorithm 3: Proposed EPA-SORI algorithm.
(1) Input:y, H, L, K, σn2, Es, ϑ.(2) Output:tL.(3)yMF=HTy(4)b=σn−2yMF.(5)Di=hiThi,i=1,2,…,2Nt.(6)D=diagD1,D2,…,D2Nt.(7)λmax⟶Nr1+Nt/Nr2,λmin⟶Nr1−Nt/Nr2.(8)β2=2/λmax+λmin(9)α2=2/1+1−λmax−λmin/λmin+λmax2(10)s^0=0, s^1=D−1yMF.(11) forl=1,2,...,K do(12)hk=Hs^l;(13)s^l+1=s^l−1+α2β2yMF−HThk−σn2Es−1s^l+α2s^l−s^l−1.(14) end for(15)s=s^K+1,Γ2=diagD1−1,D2−1,…,D2Nt−1(16) fori=1,2,…,2Nt do(17)ti0=si/1−Γi2Es(18) end for(19) repeat(20) fori=1,2,…,2Nt do(21)ηi=argminui∈Ωui−til−12(22) end for(23)hl=σn−2Hη(24)m=b−HThl(25)tl=Γ2Tm+ηl(26)tl=1−ϑtl−1+ϑtl, the damping factor ϑ∈0,1.(27)l=l+1.(28) until convergence orl>L.Figure 1
Functional diagram of EPA-SORI processing and detector.


## 5. Convergence Performance Analysis

This section mainly analyzes the convergence performance of the EPA-SORI algorithm and compares it with that of the EPA-wNSA algorithm. Note that the convergence performance of the EPA-SORI and EPA-wNSA algorithms mainly depends on the convergence performance of the initialization part (i.e., SORI and wNSA).Lemma 2.
The iterative spectral radius of the proposed EPA-SORI and the EPA-wNSA satisfy the following relationship:(30)ρGsori−ρGw=1+δζ−1ζ−1<0.Proof.
First,εl is denoted as the SORI estimation error after l iterations:(31)εl=s^−s^l.
From (16), the relationship between εlεl+1 and εl−1εl can be obtained by the following:(32)εlεl+1=0I2Nt1−α2I2Ntα2Griεl−1εl.
Thus,Gsori denotes the iterative matrix of the SORI iterative algorithm which can be expressed as follows [29]:(33)Gsori=0I2Nt1−α2I2Ntα2Gri.
Assuming thatλri and λsori are the eigenvalues of Gri and Gsori, respectively. b1 and a1 are the maximum and minimum eigenvalues of Gsori, respectively. At the same time, b and a are the maximum and minimum eigenvalues of Gri, respectively. Thus, we have λriI2Nt−Gri=0 and λsoriI4Nt−Gsori=0. Then, after substituting (33) into λsoriI4Nt−Gsori=0, we have λsori2+α2−1/λsoriα2I2Nt−Gri=0. Thus, compared with λriI2Nt−Gri=0, it can be derived as follows:(34)λsori2−λriλsori−1α2−1=0.
Next, by substitutingb into (34), the maximum eigenvalue b1 of Gsori can be obtained by b1=α2b±m/2, where m=α2b2−4α2−1. And since α2=2/1+1−b2, m can be calculated as follows:(35)α2b2−4α2−1=4b21+1−b22−41−1−b21+1−b2=4b2−41+1−b21−1−b21+1−b22=0.
Based on the above derivation, we haveb1=α2b/2. And after simplification, the spectrum radius of Gsori can be calculated by the following:(36)ρGsori=b1=α2−1≈NtNr=ζ.
Assuming thatGw denotes the iterative matrix of the wNSA iterative algorithm in a real-valued system, which is given by the following:(37)ρGw=1−δρD−1W,where δ denotes the weight factor of EPA-wNSA, and D satisfies the following approximation:(38)Di,j≈Nr+σn2Es,i=j,0,i≠j.
Thus, we have the following:(39)ρGw=1−δρWNr+σn2/Es.
(39) can be further simplified by substituting the minimum eigenvalue λmin of W. Furthermore, σn2/Es can be neglected since it is much smaller than Nr. Thus, we have the following:(40)ρGw=1−δ1−ζ2.
Comparing the spectrum radiusGsori and Gw, we have the following:(41)ρGsori−ρGw=1+δζ−1ζ−1,where 0<δ<1,0<ζ<1; hence 1+δζ−1ζ−1<0. Thus, Lemma 2 is proved. The convergence rate R of the iterative algorithm is closely related to the spectral radius of the iterative matrix G, i.e., R=−lgρG. Hence, the smaller the spectral radius of the iteration matrix, the greater the convergence speed. From (30), Gsori has a smaller spectrum radius than Gw; it means that the proposed EPA-SORI algorithm exhibits favorable convergence performance.


## 6. Computational Complexity Analysis

In this section, the computational complexity is given by the number of real-valued multiplications (RMULs). As shown in Figure1, the entire calculation process is divided into two parts: the initialization process and the iterative processes of the EPA-SORI algorithm.Firstly, in the initialization process, the number of RMULs required byyMF and b are 4NtNr and 2Nt, respectively. Since D is a diagonal matrix, thus the complexity is 4NtNr. Please note here that compared with matrix multiplication, the complexity of λmax, λmin and α2 can be ignored. Thus, the initial solution s^0 does not need to be calculated, and the number of RMULs involved in calculating s^1 is 4Nt. The next step is the SORI iteration part of the initialization process. The RMULs of hk and s^l+1 for each iteration is approximately 4NtNr and 4NtNr+6Nt, respectively. Therefore, the RMULs complexity involved in calculating the SORI iteration is K8NtNr+6Nt, where K denotes the number of iterations of SORI. The final step of this process is the initialization of t, whose computational complexity is 4Nt.Secondly, the RMULs complexity involved in calculating the EPA iteration isL8NtNr+8Nt, where L denotes the number of EPA iterations.Finally, the overall complexity of the proposed EPA-SORI algorithm is approximately as follows:(42)K+18NtNr+6Nt+L8NtNr+8Nt+4Nt.Table1 shows the computational complexity of EPA-SORI algorithm, MMSE [10], EPA-INSA [9], EPA-wNSA [32], and the Exact EP [18] algorithm. From Table 1, note that Gram matrix calculations with a complexity of up to 4Nt2Nr cannot be avoided in many of the reported methods, such as MMSE, EPA-wNSA, and the Exact EP. Additionally, MMSE and the Exact EP also involve the inversion of a matrix with a complexity of up to 4Nt3. In contrast, the proposed EPA-SORI algorithm only involves operations with complexity of about 8NtNrK+1+8NtNrL. This is because in the SORI algorithm, direct calculation of the Gram matrix is avoided by splitting calculation. For example, to calculate HTHs^1, we first calculate Hs^1, and then multiply the result with HT. Here, we have a much lower calculation complexity, which is 4NtNr+4Nt. Compared with these methods of calculating Gram matrix first and then multiplying with vector, the complexity, in this case, is greatly reduced. Since K and L are much smaller than Nt and Nr for massive MIMO systems, the complexity of the EPA-SORI algorithm is much lower than that of other algorithms.Table 1
Computational complexities comparison.
Detector schemeComputational complexityMMSE [10]4NtNt2+NtNr+NrExact EP [18]4NtNt2+NtNr+Nr+L4Nt2+8Nt+4NtEPA-INSA [9]4Nt2Nr+4KNt2+4NtNr+L4Nt2+8Nt+4NtEPA-wNSA [32]4Nt2Nr+4K+4Nt2+4NtNr+L4Nt2+8Nt+6NtProposed EPA-SORIK+18NtNr+6Nt+L8NtNr+8Nt+4Nt


## 7. Simulation Result

In this section, the BER performance results of the EPA-SORI algorithm are presented and compared with EPA-wNSA [32], MMSE [10], EP-INSA [9], and the Exact EP [18] algorithms. To fully demonstrate and verify the BER performance of the proposed EPA-SORI, simulations are performed under different modulation methods (i.e., 16/64/256QAM) and different loading factors (i.e., ζ=0.5/0.25). For some damping factor ϑ∈0,1, we set ϑ=0.5 in EPA-wNSA and EPA-SPRI algorithms according to [20]. Assuming that the base station is able to obtain perfect channel state information (CSI) and complete signal detection based on the obtained CSI.To further exhibit the convergence performance of EPA-SORI, Error-vector magnitude (EVM), which is defined asEVM=s^−s2/s2×100% [33], is considered in Figure 2. As presented in Figure 2, EPA-INSA diverges for three different modulation methods when Nr=128, Nt=64. In contrast, the convergence performance of EPA-wNSA can be improved by the optimal choice of the weighted factor, but the degree of improvement is very limited [32]. At the same time, the proposed EPA-SORI algorithm could fast converge to the accurate Exact EP algorithm with obviously fewer iterations. As the modulation order increases, the benefit brought by EPA-SORI is more obvious. Therefore, it is verified that the advantage of the proposed algorithm is in fast convergence.Figure 2
The error-vector magnitude comparison for different detectors under different modulations with fixedNt=64 and Nr=128.As shown in Figures3-5, when the Massive MIMO system is configured as 32×64, 32×128, or 64×128, the loading factor ζ is 0.5 and 0.25. In each system configuration, we consider three modulation methods: 16QAM, 64QAM, and 256QAM, and different algorithms are compared and analyzed. In Figure 4, at BER 10−3 with Nr=128, Nt=32 for 256-QAM, EPA-SORI K=3,L=7 has a better performance than that of MMSE 0.9dB. And in Figure 3, at BER 10−3 with Nr=64, Nt=32 for 256-QAM, EPA-SORI K=9,L=7 outperforms MMSE 2.2dB. It can be seen that the BER of EPA-SORI algorithm is obviously superior to that of the MMSE, and as the loading factor ζ increases, the advantage will grow further. In Figures 4 and 5, EPA-INSA K=6,L=7 has a good performance when ζ=0.25, but does not converge when ζ=0.5. Meanwhile, choosing the weight factor of EPA-wNSA helps improve the convergence performance of EPA-wNSA but makes it difficult for further improvement, especially at a low ζ. In contrast, EPA-SORI uses SORI iteration to solve the only one-time matrix inversion in EPA, which enables it to fast converge to the accurate Exact EP algorithm with low complexity. For example, in Figure 4, at BER 10−3 with Nr=128,Nt=32 for 64-QAM, EPA-SORI K=3,L=7 outperforms EPA-wNSA K=5,L=7 0.3 dB, and for 256-QAM, EPA-SORI K=3,L=7 outperforms EPA-wNSA K=5,L=73.1dB. In Figure 5, at BER 10−3 with Nr=128, Nt=64 for 64-QAM, with only 8 iterations used by EPA-SORI K=8,L=7, its performance is outperforming that of EPA-wNSA algorithm K=21,L=7 2 dB which requires 21 iterations 0.2 dB. And for 256-QAM, the advantage will grow further. In other words, the performance of the EPA-SORI algorithm is always superior to that of EPA-wNSA under the same Massive MIMO system configuration.Figure 3
The BER performance comparison of the proposed EPA-SORI, MMSE, EPA-INSA, EPA-wNSA, and the Exact EP for different modulations. At the same time, the system configurations withNt=32 and Nr=64 are considered.Figure 4
The BER performance comparison of the proposed EPA-SORI, MMSE, EPA-INSA, EPA-wNSA, and the Exact EP for different modulations, whenNt=32 and Nr=128.Figure 5
The BER performance comparison of the proposed EPA-SORI, MMSE, EPA-INSA, EPA-wNSA, and the Exact EP for different modulations whenNt=64 and Nr=128.Furthermore, a clear overview of the performance-complexity trade-off under different modulation methods and system configurations is provided in Figure6. From Figure 6, the EPA-SORI algorithm can achieve not only significantly better performance than MMSE with a significantly lower computational complexity, but also almost the same performance as the Exact EP. For example, in Figures 6(b) and 6(c), at BER 10−3 with Nr=64, Nt=32 for 64-QAM and 256-QAM, EPA-SORI K=9,L=7 outperforms MMSE 2.4dB and 3.8dB, respectively, which is very close to that of Exact EP. However, EPA-SORI K=9,L=7 only consumes 67% and 63% of computational cost of MMSE and Exact EP, respectively. It is also observed from Figure 6, EPA-wNSA can improve the BER performance by increasing the number of iterations K at the cost of a substantial increase in computational complexity. In contrast, the EPA-SORI algorithm can achieve BER performance close to Exact EP with fewer iterations, and its complexity is much lower than the Exact EP and EPA-wNSA. For example, in Figure 6(b), at BER 10−3 with Nr=64, Nt=32 for 64-QAM, by increasing the number of iterations K, the performance of EPA-wNSA K=29,L=7 is increased by 0.9 dB compared with EPA-wNSA K=16,L=7, but the complexity is increased by 11% compared with EPA-wNSA K=16,L=7. However, EPA-SORI K=8,L=7 can achieve a performance close to the Exact EP when K=8, but it only consumes 66% of complexity of EPA-wNSA K=29,L=7. Next, we compare the EPA-SORI, EP-wNSA and Exact EP algorithms in Figures 6(a) and 6(d). In Figure 6(a), at BER 10−3 with Nr=128, Nt=32 for 256-QAM, EPA-SORI K=3,L=7 outperforms EPA-wNSA K=5,L=7 2.3 dB, and the complexity is 65% of EP-wNSA K=10,L=7 and 52% of that of the Exact EP. In addition, in Figure 6(d), at BER 10−3 with Nr=128, Nt=64 for 256-QAM, EPA-SORI K=8,L=7 outperforms EPA-wNSA K=21,L=7 2 dB. At the same time, it only consumes 42% of the computational cost of EP-wNSA K=27,L=7 and 29% of that of Exact EP. In summary, compared with MMSE, Exact EP and the recently reported EPA-wNSA algorithms, the proposed EPA-SORI algorithm has a better performance-complexity trade-off advantage, which is more obvious in scenarios with high modulation order and a large number of users.Figure 6
The performance-complexity trade-off comparison for different detectors under different system configurations. Exact EP is set as the performance benchmark to compare the SNR loss of MMSE, EPA-wNSA detectors. EPA-INSA detector is not compared in this figure because it fails to converge results in very poor performance, which cannot achieve the illustrated BER level.
(a)(b)(c)(d)


## 8. Conclusion

In this paper, we propose a novel data-detection scheme, EPA-SORI detector, which can achieve the same BER performance as the Exact EP algorithm with significantly lower complexity in various massive MIMO system configurations. At the same time, compared with MMSE and the existing EPA algorithms, the proposed EPA-SORI algorithm has a better performance-complexity trade-off advantage, which is more obvious in scenarios with high modulation order and a large number of users. The proposed algorithm avoids the direct calculation of the Gram matrix. At the same time, several effective techniques (i.e., the iteration initial solution and the optimal relaxation factor) are adopted to further enhance the convergence rate and accuracy.In future work, there will be many potential applications. The proposed design can be extended to other more complex scenarios, such as the extension of EPA-SORI to decentralized architectures [34–36]. Also, it can be further combined with deep learning methods to improve performance [37, 38]. Finally, we will investigate the proposed design to more realistic channel scenarios in our future work.



---

*Source: 5336187-2022-02-16.xml*