Fast and Robust Expectation Propagation MIMO Detection via Preconditioned Conjugated Gradient
Luca Schmid Dominik Sulz Laurent Schmalen

 



Submitted version, April 29, 2024
Schmid et al.: XXX

Fast and Robust Expectation Propagation MIMO Detection via Preconditioned Conjugated Gradient
    Luca Schmid1, Dominik Sulz2, and Laurent Schmalen1
1Communications Engineering Lab, Karlsruhe Institute of Technology, 76187 Karlsruhe, Germany

2Mathematical Institute, University of Tübingen, 72076 Tübingen, Germany

Email: 1, 2

    April 29, 2024
===========================================================================================================================================================================================================================================


We study the expectation propagation (EP) algorithm for symbol detection in massive multiple-input multiple-output (MIMO) systems. The EP detector shows excellent performance but suffers from a high computational complexity due to the matrix inversion, required in each EP iteration to perform marginal inference on a Gaussian system. We propose an inversion-free variant of the EP algorithm by treating inference on the mean and variance as two separate and simpler subtasks: 
We study the preconditioned conjugate gradient algorithm for obtaining the mean, which can significantly reduce the complexity and increase stability by relying on the Jacobi preconditioner that proves to fit the EP characteristics very well.
For the variance, we use a simple approximation based on linear regression of the Gram channel matrix.
Numerical studies on the Rayleigh-fading channel and on a realistic 3GPP channel model reveal the efficiency of the proposed scheme, which offers an attractive performance-complexity tradeoff and even outperforms the original EP detector in high multi-user inference cases where the matrix inversion becomes numerically unstable.
This work has received funding in part from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 101001899) and in part from the German Federal Ministry of Education and Research (BMBF) within the project Open6GHub (grant agreement 16KISK010).


    Expectation propagation, massive MIMO detector, conjugate gradient, low complexity, 6G.




§ INTRODUCTION

Massive MIMO is discussed as a key technology for future cellular communication systems due to its ability to serve multiple users on the same time and frequency resource with a high spectral efficiency by leveraging diversity and multiplexing gains whose advantages scale with the number of antennas <cit.>.
Due to strong MUI, efficient symbol detection becomes a crucial bottleneck in achieving the promising theoretical guarantees in practice. Optimal detection based on the ML criterion quickly becomes intractable for growing system sizes. On the other hand, linear schemes with reduced complexity, like the LMMSE detector, show poor performance.

Many proposed MIMO detectors with a promising performance-complexity tradeoff are based on iterative algorithms that can be described as MP on probabilistic graphical models <cit.>. 
A high-accuracy MP detector is based on the EP algorithm <cit.>. To obtain global consistency between local approximations, the EP detector iteratively constructs a Gaussian approximation to the posterior distribution of the transmitted MIMO symbols to match the moments in each MIMO dimension.
It thereby achieves excellent performance for various antenna configurations <cit.>.

A major drawback of the EP detector is the cubic complexity of the required matrix inversion in each EP iteration, which becomes a prohibitive computational burden for future massive MIMO systems with 32-256 antennas.
To avoid the expensive computation of the matrix inverse, the authors in <cit.> reformulate the EP update equations using a fixed point approximation. 
The proposed approximation is quite coarse and the resulting matrix inversion-free detector suffers from slow convergence rate and severe performance degradation <cit.>. Therefore, the authors suggest falling back to the matrix inverse only for the initialization of the algorithm. The improved detector EPANet is proposed in <cit.>, where deep learning techniques are considered to optimize the detector in an end-to-end manner. The EPANet algorithm works very well for low-MUI scenarios but the performance gap to the EP detector increases in the presence of high MUI.
On the other hand, various works in literature study methods for approximating the matrix inverse, e.g., based on the Neumann-series <cit.>. However, the complexity reduction is limited if a good approximation is intended. 

In contrast to the works mentioned above, we argue that the EP detector does not require the explicit computation of the full matrix inverse. Instead, the problem can be divided into the simpler subtasks of finding the diagonal elements of the inverse as well as solving a linear system of equations. 
The latter is a rigorously studied problem and there exists a variety of efficient approximative methods such as the GS or the CG algorithm <cit.>. These methods were already considered for the MIMO detection problem to approximate the LMMSE detector <cit.> and in <cit.>, the error of the CG algorithm was studied theoretically in the context of EP for the asymptotic case.

In this work, we compare the suitability of different approximative linear solvers in relieving the computational burden in the EP detector.
Leveraging an intrinsic property of EP, we propose a well-suited preconditioner to substantially reduce the required number of CG iterations.
For the diagonal elements of the inverse, we rely on a simple, yet effective approximation based on linear regression of the Gram channel matrix.
Numerical studies with varying MUI scenarios on the Rayleigh fading channel as well as on a more realistic 3GPP 3D MIMO UMa channel model demonstrate the effectiveness of the novel method and provide valuable insights into both applied approximations.
The proposed algorithm proves to be numerically more stable while achieving competitive results to the original EP detector with lower complexity.



§ SYSTEM MODEL AND PRELIMINARIES



 §.§ System Model
  
We consider the multi-user MIMO channel, where N_t single-antenna users simultaneously transmit information to a base station that is equipped with N_r antennas. The observation y∈ℂ^N_r at the receiver is modeled by

    y = Hx + n,

where H∈ℂ^N_r× N_t is the channel matrix, x∈ℳ^N_t is the vector of transmitted symbols and n∼𝒞𝒩(n:0,2σ^2 I_N_r) is an additive circular-symmetric complex Gaussian noise vector. The transmit symbols are independently and uniformly sampled from a QAM constellation ℳ⊂ℂ of size |ℳ|=M. We define the SNR in dB at the receiver as

    𝗌𝗇𝗋 := 10 log_10( 𝔼[ ‖Hx‖_2^2 ]/𝔼[ ‖n‖_2^2 ]).

The complex-valued system can be decomposed into an equivalent real-valued representation

    y = Hx + n,   H = 
         [   Re{H} - Im{H};   Im{H}   Re{H} ]∈ℝ^2N_r× 2N_t

with y = (Re{y}, Im{y})^ T, x = (Re{x}, Im{x})^ T and n = (Re{n}, Im{n})^ T, which will become convenient later. The complex-valued QAM constellation ℳ is accordingly transformed to a real-valued PAM constellation ℳ.

The goal of MIMO detection is to infer the transmit symbols x from the observation y at the receiver.
The likelihood distribution can be expressed as

    p(y|x) ∏_n=1^2 N_t1_{x_n ∈ℳ},     x∈ℝ^2 N_t

where the indicator functions 1_{x_n ∈ℳ} constrain the real-valued variables x_n to takes discrete values of the constellation ℳ.
Optimal detection with respect to the SER is based on the symbol-wise ML estimate, i.e., on the marginal distributions of eq. (<ref>):

    x̂_n,ML  = _x_n ∈ℝ p(y | x_n)  1_{x_n ∈ℳ}.

We assume the availability of CSI, i.e., the receiver has perfect knowledge of H and σ^2.
For the modeling of H, we follow two approaches:



  §.§.§ Rayleigh-fading Channel Model

Each element of H is independently sampled from a circular complex standard Gaussian distribution 𝒞𝒩(0,1). Albeit the rather unrealistic assumption of uncorrelated channel coefficients in practice, this channel model is well-studied, enabling a good comparison of this work to other methods in the literature, like <cit.>.



  §.§.§ 3GPP 3D MIMO UMa NLOS Channel Model

For a more realistic model, we consider an UMa scenario from the 3GPP 3D MIMO channel model <cit.> with a NLOS propagation.
We define a configuration similar to <cit.>: N_t single-antenna users are randomly dropped within a radius of 10-500m from a base station in a 120°-cell sector and follow a linear trajectory with a constant speed of 30km/h (8.3m/s). The base station is equipped with N_r/2 dual-polarized antennas that are installed in a rectangular planar array at 25m height. The center frequency is set to 2.53GHz, and the bandwidth of 20MHz is divided among 1024 sub-carriers.
As opposed to <cit.>, this work does not focus on the exploitation of correlations in the frequency and time domain. To sample more independent channel realizations, we subsample the sub-carriers with a rate of 32 and only sample 1 OFDM symbol per channel realization.
Furthermore, we assume perfect power control, i.e., 

    1/N_r∑_i=1^N_r‖H_ij‖_2^2 = 1,    j=1,…,N_t.




 §.§ Expectation Propagation for MIMO Detection

Exact ML detection in (<ref>) is typically intractable because marginalization on p(x|y) becomes prohibitive for N_t≫ 1.
EP is an approximate inference method that simplifies exact Bayesian inference by approximating posterior distributions with members of an exponential family <cit.>. 
In the remainder of this section, we discuss the application of the EP method for MIMO detection as it was proposed by Céspedes et al. in <cit.>. 
For an excellent in-depth treatment of the general EP method, we refer the reader to <cit.> and <cit.>.

To simplify the distribution in (<ref>), we approximate 1_{x_n ∈ℳ} by a 1-dimensional Gaussian in canonical form

    1_{x_n ∈ℳ}≈exp( γ_n x_n - 1/2Λ_n x_n^2 ) =: t_n(x_n),

with free parameters γ_n and Λ_n>0.
This significantly facilitates inference, since the resulting distribution

    q(x) := 𝒩(y: Hx, σ^2I_2N_r) ∏_n=1^2N_t t_n(x_n)

belongs to the Gaussian exponential family ℱ_G for which inference is tractable <cit.>.
In this context, the meaning of tractable inference is twofold: First, we can compute the moment parameters associated with the sufficient statistics of q(x), i.e., the covariance matrix and mean

    Σ   = ( σ^-2H^ TH + diag(Λ) )^-1, 
    μ   = Σ( σ^-2H^ Ty + γ),

where we have introduced the vector notations Λ=(Λ_1,…,Λ_2N_t)^ T and γ=(γ,…,γ_2N_t)^ T. 
Second, marginal inference for ℱ_G is tractable and is directly given by
q(x_n) = 𝒩(x_n: μ_n, Σ_nn).
EP is an iterative method to determine the parameters γ_n and Λ_n > 0, n=1,…,2N_t such that the KL divergence D_KL(p(x|y) ‖ q(x) ) is minimized. As shown in <cit.>, this is equivalent to matching the expected sufficient statistics of p(x|y) and q(x), also known as moment matching. 
The EP algorithm approaches the moment matching condition by sequentially matching the moments of marginal distributions in the factorization in (<ref>).
In each iteration ℓ, the parameter updates γ_n^(ℓ+1) and Λ_n^(ℓ+1) are chosen such that the moments of the marginal Gaussian distribution

    q^(ℓ)(x_n)  exp( γ_n^(ℓ+1) x_n - 1/2Λ_n^(ℓ+1) x_n^2 )/exp( γ_n^(ℓ) x_n - 1/2Λ_n^(ℓ) x_n^2 )

match with the mean μ_p_n^(ℓ) and variance σ_p_n^2(ℓ) of the marginal tiled distribution

    q^(ℓ)(x_n)  1_{x_n ∈ℳ}/exp( γ_n^(ℓ) x_n - 1/2Λ_n^(ℓ) x_n^2 ),

for all n=1,…,2N_t, respectively.
We denote by q^(ℓ)(x_n)=𝒩(x_n:μ_n^(ℓ),σ_n^2(ℓ)) the marginals of q^(ℓ)(x) in iteration ℓ.
Solving the moment matching conditions for γ_n^(ℓ+1) and Λ^(ℓ+1) leads to the EP update equations

    Λ_n^(ℓ+1)   = β( 1/σ_p_n^2(ℓ) - 1-σ_n^2(ℓ)Λ_n^(ℓ)/σ_n^2(ℓ)) + (1-β) Λ_n^(ℓ), 
    γ_n^(ℓ+1)   = β( μ_p_n^(ℓ)/σ_p_n^2(ℓ) - μ_n^(ℓ)/σ_n^2(ℓ) + γ_n^(ℓ)) + (1-β)γ_n^(ℓ),

where β∈ [0,1] is a smoothing parameter.
After convergence or a fixed number of EP iterations L, the obtained approximation q(x) ≈ p(x|y) significantly simplifies the MIMO detection problem in (<ref>) to a nearest neighbor decision

    x̂_n,EP = _x_n ∈ℳ‖ x_n - μ_n^(L)‖_2^2.
 



§ THE EPICG ALGORITHM
 
In this section we give a detailed description of the proposed Expectation Propagation integrated Conjugate Gradient (EPiCG) algorithm.
To perform marginal inference on the multivariate Gaussian distribution q(x), i.e., to obtain the mean μ_n^(ℓ) and variance σ_n^2(ℓ) of the marginals q^(ℓ)(x_n), the EP detector computes the covariance matrix Σ in (<ref>) for each iteration ℓ.
The matrix inversion operation is unfavorable from a numerical point of view since it can lead to numerical instabilities and further imposes a complexity of the order 𝒪((2N_t)^3). This is typically unaffordable for modern massive MIMO system sizes and thus inhibits widespread employment of the EP detector in practice. In addition, the solution can become less accurate for ill-conditioned channel matrices, limiting the detection performance <cit.>.

We argue that marginal inference on q(x) does not require the computation of the full matrix inverse. 
We only require the diagonal elements of the inverse to obtain the variance of the marginals q^(ℓ)(x_n); their mean is provided as the solution of a linear system of equations.[Note that the GaBP algorithm can provide the mean and variance of the marginals in one MP algorithm and thus seems predestined for this problem. However, due to the underlying fully connected graph, we observed non-convergent behavior in most cases.]
Based on this observation, we propose to divide the  marginal inference problem into two separate and simpler sub-tasks.
In the following, we discuss both subtasks and propose approximative solutions with reduced complexity. Note that although we are treating both problems separately, the intrinsic moment matching condition of the EP algorithm combines both tasks in each EP iteration ℓ.



 §.§ Marginal Inference: Mean

As an alternative to (<ref>), the mean values μ_n^(ℓ) can be computed without the matrix inversion by solving the linear system of equations

    A^(ℓ)μ^(ℓ) = b^(ℓ),   ℓ=0,…,L

with 

    A^(ℓ)   := σ^-2H^ TH + diag(Λ^(ℓ)), 
    b^(ℓ)   := σ^-2H^ Ty + γ^(ℓ).

Solving a linear system is generally numerically more stable than computing the inverse, especially if the inverse is only part of a larger computation <cit.>, as it is the case in (<ref>).
Several stable and rigorously studied numerical methods exist to solve linear systems. We focus on iterative methods since they provide low-complexity approximative solutions upon early termination.
Besides the well-known GS method <cit.>, we consider Krylov subspace methods in the following. 



  §.§.§ Krylov Subspace Methods

In numerical linear algebra, Krylov subspace methods are a well-understood tool to iteratively find approximations to solutions of linear systems. The key idea is to find the best approximation of the solution in a subspace, the so-called Krylov space. In each iteration, this subspace is enlarged and a better approximation is constructed in the larger Krylov subspace. Standard algorithms for the Krylov subspace method are the CG algorithm and Lanczos/Arnoldi-based methods such as the FOM and BiCG <cit.>. The latter methods involve a solution of a k × k linear system in the k^th step, which becomes inefficient for large k. The CG algorithm, which was originally published in <cit.>, is specifically tailored to symmetric positive definite problems like the matrices in (<ref>), and each iteration has the same computational complexity. However, the CG algorithm is known to suffer from slow convergence if the problem is badly conditioned <cit.>. We can substantially accelerate convergence by applying a preconditioner to the CG that is specifically tailored to the characteristics of our considered problem. For more details on preconditioning, we refer the reader to <cit.>.


  §.§.§ pCG Method

Suppose we have a linear system Ax = b for a symmetric, positive definite matrix A∈ℝ^N × N and a matrix B, such that B≈A^-1. The matrix B is called the preconditioner and has to be chosen adequately. For an arbitrary initialization x^(0), g^(0) = Ax^(0) - b, d^(0) = Bg^(0) and ρ^(0) = ⟨Bg^(0),g^(0)⟩ the pCG algorithm reads for k=0,1,… as 

    x^(k+1)   = x^(k) - α^(k)d^(k),       g^(k+1) = g^(k) - α^(k)Ad^(k), 
    d^(k+1)   = Bg^(k+1) + β^(k)d^(k),       ρ^(k+1) = ⟨Bg^(k+1),g^(k+1)⟩,

where α^(k) = ρ^(k)/⟨Ad^(k),d^(k)⟩ and β^(k) = ρ^(k+1)/ρ^(k). The computational cost per pCG step is dominated by the matrix-vector product Ad^(k), which leads to a global complexity in the order of 𝒪(kN^2). The pCG algorithm reaches (except for round-off errors) the exact solution of the linear system after k=N iterations. For approximative solutions with k<N, a reasonable stopping criterion is √(ρ^(k))≤tol·b_2.

We apply the pCG algorithm to solve the system in (<ref>).
To obtain fast convergence of the pCG method, the choice of a suitable preconditioner is crucial. In this special setting, the convergence of the EP algorithm leads to decreasing marginal variances which implies large entries in Λ^(ℓ) (cf. (<ref>)) and hence strong diagonal dominance of A^(ℓ). This behavior is also observed in our experiments. Hence, it is natural to choose a diagonal, so-called, Jacobi preconditioner <cit.>. For A^(ℓ)=(A_nm)_n,m=1^2N_t, a good choice is B=diag(1/A_nn)_n=1^2N_t. This choice for B requires no matrix inversion and allows for very efficient matrix-vector products in the pCG algorithm.



 §.§ Marginal Inference: Variance

The marginal variances σ_n^2(ℓ) are obtained in each iteration ℓ of the EP algorithm by extracting the diagonal elements of the covariance matrix Σ^(ℓ) = A^-1(ℓ) in (<ref>).
To avoid the full matrix inversions for reduced complexity, we rely on an approximation similar to <cit.>, which is based on the Neumann-series expansion terminated after the first term:

    σ^2(ℓ)_n = Σ^(ℓ)_nn = 1/A^(ℓ)_nn + ∑_m ≠ nA^(ℓ)_nm A^(ℓ)_mn/A^(ℓ)_mm1/A^2(ℓ)_nn≈1/A^(ℓ)_nn.

Inspired by <cit.>, we introduce correction parameters to improve the approximation in (<ref>):

    σ̃^2(ℓ)_n := max{α_1^(ℓ)1/A^(ℓ)_nn + α_2^(ℓ), 1/A^(ℓ)_nn}.

Compared to <cit.>, the parameters apply a linear transformation on 1/A_nn (instead of on A_nn), and we lower bound σ̃^2_n by 1/A_nn instead of by a constant (ϵ = 10^-12).
Optimizing α_1^(ℓ) and α_2^(ℓ) in each iteration ℓ based on a labeled dataset 𝒟^(ℓ)_train = { (A^(ℓ), A^-1(ℓ))_d, d=1,…,D} of representative matrix-inverse pairs can be interpreted as linear regression[In numerical studies, we found that a piecewise exponential regression with two intervals can fit the data much better. However, the overall gains in terms of EP detection performance were insignificant.] between the dependent variables (A^-1(ℓ))_nn=σ^2(ℓ)_n and the regressors 1/A^(ℓ)_nn.
Following the method of least squares, we obtain

    ( α_1^(ℓ), α_2^(ℓ)) = _(α_1, α_2) ∈ℝ^2∑_𝒟^(ℓ)_train∑_n=1^2N_t( α_1 1/A^(ℓ)_nn + α_2 - σ^2(ℓ)_n )^2,

which can be minimized in closed form <cit.>.
We observe that this approximation works remarkably well for EP iterations ℓ>0 and for ℓ=0 if the Gram channel matrix H^ TH is strongly diagonal dominant. In the remaining cases, the initialization with the inverse Σ^(0) has shown to significantly increase accuracy.



 §.§ The Proposed EPiCG Algorithm

Combining the discussions in Sec. <ref> and <ref>, we propose the EPiCG algorithm  with two different initializations as our main contribution of this work.
We apply the original EP detector (cf. Sec. <ref>) with the following modifications:

    
  * EPiCG: Compute the means μ_n^(ℓ) by solving the linear system in (<ref>) with the pCG method and Jacobi-preconditioning. Approximate the variances σ_n^2(ℓ) by the strategy proposed in Sec. <ref>.
    
  * EPiCG-Σ^(0): For ℓ = 0, initialize μ_n^(0) and σ_n^2(0) via the inverse A^-1(0)=Σ^(0) as in the original EP detector. Otherwise for ℓ > 0, apply the EPiCG algorithm.





§ EXPERIMENTAL STUDY

We evaluate the EPiCG algorithm for MIMO detection and focus in particular on the analysis and distinction between the two proposed approximations. In the following numerical experiments, we use a 16-QAM constellation and fix the number of receive antennas to N_r=128. To cover different MUI scenarios, we study both N_t=64 and N_t=128 single-antenna users. The results are based on 5000 channel matrices H that were randomly sampled from the channel models introduced in Sec. <ref>, respectively.



 §.§ Approximation of the Marginal Mean
 
We start by evaluating the suitability of different approximative methods for solving the linear system (<ref>) in each EP iteration ℓ. Figure <ref> plots the RRE ‖A^(ℓ)x̂^(k) - b^(ℓ)‖_2/‖b^(ℓ)‖_2 over the steps k of the GS, CG and pCG method, respectively.
The results are shown for the Rayleigh fading channel with N_t=128 and 𝗌𝗇𝗋=19dB. 
We compare the methods with two non-iterative approaches from literature: solving the linear system using the inverse x̂ = A^-1(ℓ)b^(ℓ) as proposed in the original EP detector <cit.>, and the approximation used in the EPA algorithm, i.e., avoiding to solve the linear system completely by assuming x̂_n = μ_p_n^(ℓ). However, this assumption is not correct, but only applies when the fixed points of EP are reached <cit.>.

In EP iteration ℓ=0, the GS solver is outperformed by the CG and pCG method for k>16. The latter behave almost identically. The benefit of the Jacobi preconditioner becomes apparent for the subsequent EP iterations ℓ>0, where the matrix A^(ℓ) becomes more and more diagonal-dominant. 
While the conventional CG solver converges similarly or slower compared to ℓ=0, the pCG method significantly increases the convergence rate in each subsequent EP iteration ℓ. 
Since the computational complexity scales quadratically for each step of the considered solvers, the plots indicate a precision-complexity tradeoff.
For instance in ℓ=1, the pCG solver only requires k=16 steps to reach the same RRE=10^-4 as the inverse-based solution of the EP detector. 
Note that the RRE is inherently available within the pCG steps which offers great flexibility for controlling this tradeoff.

The approximation of the EPA algorithm is very coarse in ℓ=0 but quickly improves with the convergence of the EP algorithm, i.e., with increasing ℓ. It approximatively requires k=4 iterations of the pCG algorithm to reach the RRE of the EPA algorithm in each EP iteration ℓ.



 §.§ Detection Performance

We compare the performance of the proposed EPiCG algorithm to the LMMSE detector, the EPANet algorithm <cit.> as well as to the original EP detector <cit.>. We fix β=0.1 for the EP and EPiCG algorithm and we perform L=10 iterations for all iterative detection schemes. We also consider the initialization variant of the EPiCG-Σ^(0) algorithm for the EPANet algorithm as suggested in <cit.>, accordingly denoted by EPANet-Σ^(0).
We fix the stopping criterion of the pCG algorithm to tol(ℓ)=(10^-3,10^-4,10^-5,10^-5,…,10^-5). For the Rayleigh-fading channel, the precision is thereby similar to the matrix inverse (cf. Fig. <ref>) which enables a good comparison since the main difference between the EPiCG algorithm and the EP detector stems from the variance approximation.
To create 𝒟_train, we randomly generate D=100 transmission scenarios with 𝗌𝗇𝗋 values sampled from a uniform distribution with bounds chosen such that the SER of the resulting EPiCG detector performance is in the target SER interval [10^-2;10^-1].


Figure <ref> evaluates the SER over the SNR for the Rayleigh-fading channel. In the low-MUI scenario with N_t=64, both the EPANet and the EPiCG algorithm approach the EP detector performance. 
In the high-MUI scenario with N_t=N_r, the EPANet detector behaves similarly to the LMMSE detector and only reaches competitive performance when being initialized with Σ^(0). The EPiCG detector follows the waterfall behavior of the EP detector with an 𝗌𝗇𝗋 offset of 0.75dB for SER=10^-2. We observed that the remaining gap between the EPiCG algorithm compared to the EP detector comes from the variance approximation in iteration ℓ=0 as the dominant source of error. The EPiCG-Σ^(0) algorithm consequently closes this gap.



We further evaluate the performance of the considered detectors on the 3GPP 3D MIMO UMa NLOS channel model in Fig. <ref>. The significantly increased correlation in this channel model and the consequently decreased diagonal dominance of the Gram channel matrices generally degrades the performance of both the EPANet detector and the EPiCG algorithm. For N_t=64, the EPiCG-Σ^(0) algorithm reduces the 𝗌𝗇𝗋 gap to EP detection from 4dB to 2dB at an SER=10^-2, compared to the EPiCG and EPANet-Σ^(0) detectors.
In the high-MUI scenario with N_t=128, the EPiCG-Σ^(0) algorithm outperforms the EP detector. We conjecture that the performance of the EP detector degrades due to an unstable solution of eq. (<ref>) via the matrix inverse. To verify this, we replace the computation of the mean in the EP detector by the pCG algorithm with tol=10^-5 and only use the diagonal elements of the matrix inverse to obtain the variances. The result is denoted by EP-pCG in Fig. <ref> and confirms the conjecture.



 §.§ Complexity Analysis

Finally, we evaluate the proposed EPiCG detection algorithms from a computational complexity perspective. Table <ref> provides the total number of pCG steps that are required within the EPiCG algorithm over all EP iterations, i.e., to solve L linear systems of dimension 2N_t× 2N_t. Since we use a precision-based stopping criterion instead of a fixed number of pCG steps, the numbers represent rounded mean values over all 5000 evaluation samples.
The EPiCG-Σ^(0) variant uses the matrix inverse for the computation of the marginal means in EP iteration ℓ=0 and requires approximately half of the pCG steps compared with the EPiCG algorithm.
Note that optimizing the distribution of the pCG steps over the EP iterations might have the potential to improve further the performance-complexity tradeoff, which we leave for future work.



§ CONCLUSION

We propose the EpiCG algorithm as a faster and more robust method for massive MIMO detection.
Replacing the matrix inversion for marginal inference in the EP algorithm by the application of the pCG method reduces the computational complexity from 𝒪(N^3) to 𝒪(kN^2) operations in each EP iteration.
This offers a flexible performance-complexity tradeoff, while the accuracy is only mildly reduced due to an approximation of the marginal variance. 
Simulations on realistic MIMO channels illustrate that the pCG algorithm furthermore stabilizes the computation compared to the matrix inversion for highly-correlated channels with strong MUI.


10

url@samestyle








yang_fifty_2015
S. Yang and L. Hanzo, “Fifty years of MIMO detection: The road to
  large-scale MIMOs,” IEEE Commun. Surveys Tuts., vol. 17,
  no. 4, pp. 1941–1988, 2015.

wu_low-complexity_2014
S. Wu, L. Kuang, Z. Ni, J. Lu, D. Huang, and Q. Guo, “Low-complexity iterative
  detection for large-scale multiuser MIMO-OFDM systems using approximate
  message passing,” IEEE J. Sel. Topics Signal Process., vol. 8,
  no. 5, pp. 902–915, Oct. 2014.

minka_expectation_2013
T. P. Minka, “Expectation propagation for approximate Bayesian
  inference,” in Proc. UAI, Seattle, WA, USA, 2001, pp. 362–369.

cespedes_expectation_2014
J. Céspedes, P. M. Olmos, M. Sánchez-Fernández, and F. Perez-Cruz,
  “Expectation propagation detection for high-order high-dimensional
  MIMO systems,” IEEE Trans. Commun., vol. 62, no. 8, pp.
  2840–2849, Aug. 2014.

tan_low-complexity_2019
X. Tan, Y.-L. Ueng, Z. Zhang, X. You, and C. Zhang, “A low-complexity
  massive MIMO detection based on approximate expectation
  propagation,” IEEE Trans. Veh. Technol., vol. 68, no. 8, pp.
  7260–7272, Aug. 2019.

ge_improving_2021
Y. Ge, X. Tan, Z. Ji, Z. Zhang, X. You, and C. Zhang,
  “enImproving approximate expectation propagation
  massive MIMO detector with deep learning,”
  enIEEE Wireless Commun. Lett., vol. 10,
  no. 10, pp. 2145–2149, 2021.

zhang_expectation_2018
Y. Zhang, Z. Wu, C. Li, Z. Zhang, X. You, and C. Zhang, “Expectation
  propagation detection with Neumann-series approximation for
  massive MIMO,” in Proc. IEEE SiPS, Cape Town, South Africa,
  2018.

GvL13
G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed.1em
  plus 0.5em minus 0.4emPhiladelphia, PA: Johns Hopkins University
  Press, 2013.

dai_low-complexity_2015
L. Dai, X. Gao, X. Su, S. Han, C.-L. I, and Z. Wang, “Low-complexity
  soft-output signal detection based on Gauss-Seidel method for
  uplink multiuser large-scale MIMO systems,” IEEE Trans.
  Veh. Technol., vol. 64, no. 10, 2015.

wei_learned_cg_2020
Y. Wei, M.-M. Zhao, M. Hong, M.-j. Zhao, and M. Lei, “Learned conjugate
  gradient descent network for massive MIMO detection,”
  IEEE Trans. Signal Process., vol. 68, pp. 6336–6349, 2020.

takeuchi_rigorous_2017
K. Takeuchi and C.-K. Wen, “Rigorous dynamics of expectation-propagation
  signal detection via the conjugate gradient method,” in Proc. IEEE
  SPAWC, Sapporo, Japan, Jul. 2017.

khani_adaptive_2020
M. Khani, M. Alizadeh, J. Hoydis, and P. Fleming, “Adaptive neural signal
  detection for massive MIMO,” IEEE Trans. Wireless Commun.,
  vol. 19, no. 8, pp. 5635–5648, Aug. 2020.

3gpp:channel_model
Study on 3D Channel Model for LTE, 3GPP TR, 2015.

bishop_pattern_2006
C. M. Bishop, Pattern Recognition and Machine Learning.1em plus
  0.5em minus 0.4emNew York: Springer, 2006.

Seeger2005ExpectationPF
M. W. Seeger, “Expectation propagation for exponential families,” University
  of California, Berkeley, CA, USA, Technical Report, 2006.

kosasih_graph_2022
A. Kosasih, V. Onasis, V. Miloslavskaya, W. Hardjawana, V. Andrean, and
  B. Vucetic, “enGraph neural network aided
  MU-MIMO detectors,” enIEEE J. Sel.
  Areas Commun., vol. 40, no. 9, pp. 2540–2555, Sep. 2022.

croz1992stability
J. J. D. Croz and N. J. Higham, “Stability of methods for matrix inversion,”
  IMA J. of Numerical Analysis, vol. 12, no. 1, pp. 1–19, 1992.

SiSz2007
V. Simoncini and D. B. Szyld, “Recent computational developments in Krylov
  subspace methods for linear systems,” Numerical Linear Algebra with
  Applications, vol. 14, no. 1, pp. 1–59, 2007.

hestenes1952cg
M. R. Hestenes and E. Stiefel, “Methods of conjugate gradients for solving
  linear systems,” Journal of Research of the National Bureau of
  Standards, vol. 49, no. 1, 1952.




