An Improved Jacobi-Based Detector for Massive MIMO Systems

Abstract

Massive multiple-input-multiple-output (MIMO) is one of the key technologies in the fifth generation (5G) cellular communication systems. For uplink massive MIMO systems, the typical linear detection such as minimum mean square error (MMSE) presents a near-optimal performance. Due to the required direct matrix inverse, however, the MMSE detection algorithm becomes computationally very expensive, especially when the number of users is large. For achieving the high detection accuracy as well as reducing the computational complexity in massive MIMO systems, we propose an improved Jacobi iterative algorithm by accelerating the convergence rate in the signal detection process.Specifically, the steepest descent (SD) method is utilized to achieve an efficient searching direction. Then, the whole-correction method is applied to update the iterative process. As the result, the fast convergence and the low computationally complexity of the proposed Jacobi-based algorithm are obtained and proved. Simulation results also demonstrate that the proposed algorithm performs better than the conventional algorithms in terms of the bit error rate (BER) and achieves a near-optimal detection accuracy as the typical MMSE detector, but utilizing a small number of iterations.

1. Introduction

Massive multiple-input-multiple-output (MIMO) is an emerging technology for communication application which contributes a promising technology for the wireless sensor networks (WSNs) [1,2,3] and the fifth generation (5G) wireless communications [4]. In such systems, the base station (BS) was equipped with hundreds of antennas serving tens of single-antenna users in the same frequency band [5,6,7]. Benefit from the massive MIMO system is capable of achieving higher multiplexing and diversity gains compared with the conventional small-scale MIMO system [8,9].

However, due to the increasing dimension, the signal detection in the uplink may become a challenge of the massive MIMO system [10]. According to [11], the maximum likelihood (ML) detector is an optimal algorithm. However, its substantially high computational complexity contributes the bottleneck of massive MIMO systems. Therefore, assorted detection algorithms have been proposed. Some nonlinear detection algorithms such as sphere decoder (SD) and its variants have been proposed to achieve near-optimal performance with low complexity [12,13]. However, the complexity is still unfavorable because of the large dimension of the systems or the high modulation order. Thus, one can resort to the linear detection algorithm, such as the minimum mean square error (MMSE) with the near-optimal detection accuracy [14]. However, the inverse of the large-dimensional covariance matrices required by MMSE would result in a prohibitively high complexity.

Given this challenge, low-complexity approximate matrix inversion has drawn substantial attention. The typical algorithms can be divided into two categories. One typical category is based on the series expansion, such as the Neumann series expansion (NSE) detection algorithm [15]. The NSE algorithm utilizes the first few terms of the series expansion to approximate matrix inverse. However, its performance suffers from a significant loss with the scaling up massive MIMO and only marginal reduction in complexity can be achieved.

Another category is based on the iterative algorithms derived from linear equations. To expand, it converts the matrix inversion process into solving linear equations, where the Gauss Seidel (GS) based iterative algorithm is one of the typical algorithms and utilizes the most up-to-date values at each iteration, leading to the faster convergence rate and lower complexity than NSE detection algorithm [16]. Furthermore, the symmetric successive over-relaxation (SSOR) algorithm [17] can also be applied, where different weights are put in the iteration structure. To enhance the robustness of the SSOR algorithm, the optimal parameter of weights is determined by exploiting channel characteristics, yielding a fast convergence and a high accuracy. Unfortunately, each iteration of the GS and SSOR calculations is serial in nature, where each component of the new iteration depends on all of the previously calculated components. Fortunately, the Jacobi algorithm is a simple iterative algorithm with high parallelism [18]. Some improved Jacobi-based algorithms have been proposed, one of which is the damped Jacobi (DJ) algorithm [19]. The algorithm has great convergence performance when the number antennas at the BS is far greater than the number of antennas of the users. However, in the uplink communications, the increase of the users will make the requirements of DJ algorithm no longer satisfied, resulting in a great performance degradation.

Motivated by the above concerns, we focus our attention on further improving the convergence rate and calculation accuracy of the Jacobi algorithm with low-complexity calculation. Based on the Jacobi algorithm, we invoke the steepest descent (SD) method [20] and whole-correction method [21] to modify the iteration process. The major contributions of our work can be listed as below: (1) The SD method [20] is invoked in Jacobi algorithm to improve the convergence rate. As the SD method provides an efficient searching direction, the combined SD and Jacobi algorithm is capable of satisfying the convergence condition. Furthermore, compared to the conventional Jacobi algorithm, the combined algorithm is more efficient to accommodate a different number of users. (2) Note that, when the iterative approximate solution is close to the exact solution, the SD method will no longer provide an efficient searching direction, and the whole-correction method [21] is employed to improve the convergence. By exploiting the obtained, the more accurate solution can be attained by the whole-correction method.

The remainder of the paper is structured as follows. We briefly introduce the system model and present a review of previously proposed approaches in Section 2. In Section 3, a hybrid Jacobi-based iterative method is proposed, where the analysis of convergence and computational complexity of the proposed Jacobi iterative method are presented. In Section 4, numerous simulation results are demonstrated. Finally, the entire paper is concluded in Section 5.

Notation: In this paper, bold-face, lower-case letters represent the column vectors (e.g., $\mathbf{a}$), while the bold-face, upper-case letters refer to the matrices (e.g., $\mathbf{A}$). The superscript of matrix $\mathbf{A}$, ${\mathbf{A}}^T$, ${\mathbf{A}}^H$, ${\mathbf{A}}^{-1}$ and ${\mathbf{A}}^{+}$ indicate the transpose, the Hermitian transpose, the inverse and the Moore–Penrose inverse of $\mathbf{A}$, respectively. Furthermore, ${\mathbf{I}}_K$ denotes the $K\times K$ identify matrix, $\left[\mathbf{a}\mathbf{b}\right]$ denotes the inner product of the vectors, and $\parallel ·\parallel$ stands for the Euclidean norm of a vector.

2. System Model

2.1. System Model

Consider an uplink large-scale MIMO system where the BS is equipped with $N_r$ receiving antennas serving $N_t$ single-antenna users and $N_r\gg N_t$ [22]. The received signal vector, which is denoted by $\mathbf{y}={[y_1,y_2,\cdots ,y_{N_r}]}^T$, at the BS is expressed as

$$\mathbf{y}=\mathbf{Hx}+\mathbf{n},$$

(1)

where $\mathbf{x}={[x_1,x_2,\cdots ,x_{N_t}]}^T$ denotes the transmitting signal set and $x_i\in \Omega$ is the transmitted signal from the ith user with $\Omega$ being a complex signal set obtained by modulating and mapping the source message $\mathbf{s}$. $\mathbf{H}\in {\mathbf{C}}^{N_r\times N_t}$ represents a flat-fading channel matrix which has a complex-valued independent and identically distributed (i.i.d.) Gaussian entries following $CN\left(0,1\right)$, and $\mathbf{n}={[n_1,n_2,\cdots ,n_{N_r}]}^T$ is the complex-valued additive white Gaussian noise vector with zero mean and ${\sigma }^2$ variance, respectively.

The purpose of MIMO signal detection is to estimate the transmitted vector $\mathbf{x}$. The channel matrix $\mathbf{H}$ and the received signal vector $\mathbf{y}$ are utilized for the detection. The structure of the signal detection is illustrated in Figure 1.

2.2. MMSE Detector

The MMSE approach achieves a near-optimal detector in massive MIMO systems, where the resulting estimated symbol vector $\overline{\mathbf{x}}$ by utilizing MMSE detector is given by

$$\overline{\mathbf{x}}={\left({\mathbf{H}}^H\mathbf{H}+{\sigma }^2{\mathbf{I}}_{N_t}\right)}^{-1}{\mathbf{H}}^H\mathbf{y}.$$

(2)

In Equation (2), $\mathbf{A}={\mathbf{H}}^H\mathbf{H}+{\sigma }^2{\mathbf{I}}_{N_t}$ is the symmetric positive definite matrix and $\mathbf{b}={\mathbf{H}}^H\mathbf{y}$ is regarded as the output of the matched filter, respectively. The computational complexity of direct matrix inversion ${\mathbf{A}}^{-1}$ is on the order of $O\left(N_t^3\right)$, which incurs a heavy computational burden in massive MIMO systems when the large-dimensional matrix inversion needs to be performed for large number of users.

2.3. Conventional Jacobi Algorithm

To avoid the direct matrix inverse of $\mathbf{A}$, we can convert the MMSE detection into solving the linear equation as

$$\mathbf{A}\overline{\mathbf{x}}=\mathbf{b}.$$

(3)

By utilizing Jacobi algorithm [18], the matrix $\mathbf{A}$ can be decomposed as $\mathbf{A}=\mathbf{D}+\mathbf{E}$, where $\mathbf{D}$ and $\mathbf{E}$ are respectively considered as the diagonal and off-diagonal components parts of $\mathbf{A}$. Then, the iterative solution of the Jacobi algorithm is given by

$${\overline{\mathbf{x}}}^{\left(k\right)}={\mathbf{D}}^{-1}\left[\left(\mathbf{D}-\mathbf{A}\right){\overline{\mathbf{x}}}^{\left(k-1\right)}+\mathbf{b}\right],$$

(4)

where k denotes the iteration index. Since the convergence rate of the iteration is highly based on the initial setting, we consider the initial solution ${\overline{\mathbf{x}}}^{\left(0\right)}$ as ${\overline{\mathbf{x}}}^{\left(0\right)}=\left({\mathbf{D}}^{-1}-{\mathbf{D}}^{-1}\mathbf{E}{\mathbf{D}}^{-1}\right)\mathbf{b}$ in this work [23].

3. Proposed Algorithm for MIMO Detection

In this section, an improved Jacobi algorithm for MIMO detection is proposed in uplink massive MIMO systems. Specifically, the proposed algorithm is divided into two stages. In the first stage, the SD method is applied towards an efficient searching direction in the Jacobi algorithm. In the second stage, the whole-correction method is utilized to further improve the convergence rate of the algorithm. Then, a fast convergence of the proposed algorithm is proved following an analysis of the complexity of the algorithm.

3.1. Improved Jacobi Algorithm

For achieving a fast convergence of the estimation of $\mathbf{x}$, the SD method is applied to provide an efficient searching direction for Jacobi algorithm [20], which is called the SD-Jacobi iteration in this work. The combined iteration consists of two steps:

(1) During the first half of the iteration, utilize the SD method to compute symbol vector as

$${\overline{\mathbf{x}}}^{\left(k+1/2\right)}={\overline{\mathbf{x}}}^{\left(k\right)}+{\alpha }_k{\mathbf{r}}^{\left(k\right)},$$

(5)

where ${\mathbf{r}}^{\left(k\right)}=\mathbf{b}-\mathbf{A}{\overline{\mathbf{x}}}^{\left(k\right)}$ is the residual vector, and ${\alpha }_k$ represents the variable step size, respectively. The step size ${\alpha }_k$, which depends on ${\mathbf{r}}^{\left(k\right)}$, is determined by

$${\alpha }_k=\frac{\left[{\mathbf{r}}^{\left(k\right)},{\mathbf{r}}^{\left(k\right)}\right]}{\left[\mathbf{A}{\mathbf{r}}^{\left(k\right)},{\mathbf{r}}^{\left(k\right)}\right]}.$$

(6)

(2) During the second half of the iteration, conduct Jacobi iteration, which is given by

$${\overline{\mathbf{x}}}^{\left(k+1\right)}={\mathbf{D}}^{-1}\left[\left(\mathbf{D}-\mathbf{A}\right){\overline{\mathbf{x}}}^{\left(k+1/2\right)}+\mathbf{b}\right]={\mathbf{D}}^{-1}\left(\mathbf{D}-\mathbf{A}\right){\overline{\mathbf{x}}}^{\left(k+1/2\right)}+{\mathbf{D}}^{-1}\mathbf{b}=\mathbf{G}{\overline{\mathbf{x}}}^{\left(k+1/2\right)}+{\mathbf{D}}^{-1}\mathbf{b},$$

(7)

where $\mathbf{G}={\mathbf{D}}^{-1}\left(\mathbf{D}-\mathbf{A}\right)=\mathbf{I}-{\mathbf{D}}^{-1}\mathbf{A}$.

Substituting Equation (5) into Equation (7) yields

$$\begin{array}{cc}\hfill {\overline{\mathbf{x}}}^{\left(k+1\right)}& =\mathbf{G}\left({\overline{\mathbf{x}}}^{\left(k\right)}+{\alpha }_k{\mathbf{r}}^{\left(k\right)}\right)+{\mathbf{D}}^{-1}\mathbf{b}={\overline{\mathbf{x}}}^{\left(k\right)}+{\alpha }_k{\mathbf{r}}^{\left(k\right)}-{\mathbf{D}}^{-1}\mathbf{A}\left({\overline{\mathbf{x}}}^{\left(k\right)}+{\alpha }_k{\mathbf{r}}^{\left(k\right)}\right)+{\mathbf{D}}^{-1}\mathbf{b}\hfill \\ & ={\overline{\mathbf{x}}}^{\left(k\right)}+{\alpha }_k{\mathbf{r}}^{\left(k\right)}+{\mathbf{D}}^{-1}\left(\mathbf{b}-\mathbf{A}{\overline{\mathbf{x}}}^{\left(k\right)}-{\alpha }_k\mathbf{A}{\mathbf{r}}^{\left(k\right)}\right)={\overline{\mathbf{x}}}^{\left(k\right)}+{\alpha }_k{\mathbf{r}}^{\left(k\right)}+{\mathbf{D}}^{-1}\left({\mathbf{r}}^{\left(k\right)}-{\alpha }_k\mathbf{A}{\mathbf{r}}^{\left(k\right)}\right)\hfill \\ & ={\overline{\mathbf{x}}}^{\left(k\right)}+{\alpha }_k{\mathbf{r}}^{\left(k\right)}+{\mathbf{D}}^{-1}\left({\mathbf{r}}^{\left(k\right)}-{\alpha }_k{\mathbf{p}}^{\left(k\right)}\right),\hfill \end{array}$$

(8)

where $\left\{\begin{array}{c}{\mathbf{r}}^{\left(k\right)}=\mathbf{b}-\mathbf{A}{\overline{\mathbf{x}}}^{\left(k\right)}\\ {\alpha }_k=\frac{\left[{\mathbf{r}}^{\left(k\right)},{\mathbf{r}}^{\left(k\right)}\right]}{\left[{\mathbf{p}}^{\left(k\right)},{\mathbf{r}}^{\left(k\right)}\right]}\\ {\mathbf{p}}^{\left(k\right)}=\mathbf{A}{\mathbf{r}}^{\left(k\right)}.\end{array}\right.$

Equation (8) shows that, at each iteration, the involved SD method improves the convergence rate of the Jacobi iteration. However, when the iterative estimation is close to the solution, the SD method no longer provides an efficient iteration direction. Hence, a whole-correction method is applied to continue improving the convergence of the SD-Jacobi iteration afterwards, which is introduced below.

The whole-correction method is designed to solve the optimal convergence problem in a linear system [21]. It assigns the different weights to a series of approximate solutions of the iterations to obtain the optimal solution, which is presented below.

Suppose ${\overline{\mathbf{x}}}_1,{\overline{\mathbf{x}}}_2,\cdots ,{\overline{\mathbf{x}}}_m\left(m>1\right)$ are the approximate iterative solutions of the symbol vector in Equation (2), which are assumed mutually exclusive. In this work, let $\underset{m}{\underbrace{{\overline{\mathbf{x}}}_1={\overline{\mathbf{x}}}^{\left(k-m+1\right)},{\overline{\mathbf{x}}}_2={\overline{\mathbf{x}}}^{\left(k-m+2\right)},\cdots ,{\overline{\mathbf{x}}}_m={\overline{\mathbf{x}}}^{\left(k\right)}}}$, where m is the number of selected solutions, which is discussed later. Considering the different weights to the iterative solutions, the vector $\overline{\mathbf{x}}$, which represents the updated ${\overline{\mathbf{x}}}^{\left(k\right)}$, can be expressed as

$$\begin{array}{c}\overline{\mathbf{x}}=a_1{\overline{\mathbf{x}}}_1+a_2{\overline{\mathbf{x}}}_2+\cdots +a_m{\overline{\mathbf{x}}}_m,\hfill \\ a_1+a_2+\cdots +a_m=1,\hfill \end{array}$$

(9)

where the weight $a_i$ denotes the contribution ratio of ${\overline{\mathbf{x}}}_i$ to $\overline{\mathbf{x}}$.

The whole-correction method leads to a more accurate solution of the vector $\overline{\mathbf{x}}$ from the approximate solution set of ${\overline{\mathbf{x}}}_i$. Therefore, the updated solution $\overline{\mathbf{x}}$ should satisfy ${∥\mathbf{b}-\mathbf{A}\overline{\mathbf{x}}∥}_2<\underset{1\le i\le m}{min}{∥\mathbf{b}-\mathbf{A}{\overline{\mathbf{x}}}_i∥}_2$.

We are now in position to determine the weight $a_i$. Let’s select $q\in \left\{1,2,\cdots ,m\right\}$ randomly. Then, we can obtain

$$\overline{\mathbf{x}}={\overline{\mathbf{x}}}_q+\sum _{i\ne q}{a}_i\left({\overline{\mathbf{x}}}_i-{\overline{\mathbf{x}}}_q\right).$$

(10)

The residual vector regarding $a_q$ is ${\mathbf{r}}_q=\mathbf{b}-\mathbf{A}{\overline{\mathbf{x}}}_q$. Thus, the residual vector of the updated $\overline{\mathbf{x}}$ is rewritten as

$$\mathbf{r}\left(\overline{\mathbf{x}}\right)=\mathbf{b}-\mathbf{A}\overline{\mathbf{x}}=(\mathbf{b}-\mathbf{A}{\overline{\mathbf{x}}}_q)-\sum _{i\ne q}{a}_q\left(\mathbf{A}{\overline{\mathbf{x}}}_i-\mathbf{A}{\overline{\mathbf{x}}}_q\right)={\mathbf{r}}_q+\sum _{i\ne q}{a}_q\left({\mathbf{r}}_i-{\mathbf{r}}_q\right).$$

(11)

In this case, the optimization problem can be rewritten as

$${∥\mathbf{b}-\mathbf{A}\overline{\mathbf{x}}∥}_2=∥{\mathbf{r}}_q+\sum _{i\ne q}{a}_q\left({\mathbf{r}}_i-{\mathbf{r}}_q\right)∥={∥{\mathbf{r}}_q+{\mathbf{Q}}_q{\mathbf{y}}_q∥}_2=min,$$

(12)

where

$$\begin{array}{c}{\mathbf{Q}}_q=\left({\mathbf{r}}_1-{\mathbf{r}}_q,\cdots ,{\mathbf{r}}_{q-1}-{\mathbf{r}}_q,{\mathbf{r}}_{q+1}-{\mathbf{r}}_q,\cdots ,{\mathbf{r}}_m-{\mathbf{r}}_q\right),\hfill \\ {\mathbf{y}}_q={\left(a_1,\cdots ,a_{q-1},a_{q+1},\cdots a_m\right)}^T.\hfill \end{array}$$

The minimum norm solution of the optimization problem in Equation (12) is given by

$$\begin{array}{c}{\mathbf{y}}_q=-{\mathbf{Q}}_q^{+}{\mathbf{r}}_q,\hfill \\ a_q=1-{\displaystyle \sum _{i\ne q}}{a}_i.\hfill \end{array}$$

(13)

Therefore, an improved solution of $\overline{\mathbf{x}}$ can be achieved when the contribution ratio is obtained. It is verified later that this improved solution is more accurate than any arbitrary approximate solutions of the iterations. Based on the above discussion, the pseudo code of the proposed algorithm is presented in Algorithm 1, where the algorithm terminates at the K-th iteration. Note that an important issue of the proposed algorithm is related to a suitable choice of the number of approximate solutions m, which highly influences the convergence rate and the computation complexity of the algorithm. In the next subsection, we will consider the value m from two aspects.

Algorithm 1: Proposed Improved Jacobi Iterative Algorithm

Step 1) Compute the initial approximate solution of the transmitted signal $\mathbf{x}$, mark it as ${\overline{\mathbf{x}}}^{\left(0\right)}$.

Step 2) (if $m>2$) Compute $m-2$ approximate solutions of $\mathbf{x}$ using Equation (7), mark as

$\underset{m-2}{\underbrace{{\overline{\mathbf{x}}}^{\left(1\right)},\cdots ,{\overline{\mathbf{x}}}^{\left(m-2\right)}}}$, to obtain enough mutual solution which are required by the

whole-correction method.

For $k=1:m-2$

Use SD-Jacobi iteration to obtain ${\overline{\mathbf{x}}}^{\left(k\right)}$;

End for

Step 3) Employe the whole-correction method to update the approximate solution.

For $k=m-1:K$

Use SD-Jacobi iteration to obtain ${\overline{\mathbf{x}}}^{\left(k\right)}$;

Obtain the updated ${\overline{\mathbf{x}}}^{\left(k\right)}$ by $\underset{m}{\underbrace{{\overline{\mathbf{x}}}^{\left(k-m+1\right)},\cdots ,{\overline{\mathbf{x}}}^{\left(k\right)}}}$;

End for

3.2. Parameter Selection

In general from the existing iterative MIMO detection algorithms aim at the small number of iterations. In this case, number of approximate solutions obtained by the iterations is also limited. Thus, in this part, we only consider the value of m as 2, 3 and 4. Firstly, we focus on the complexity of calculating the correction coefficient. In each iteration process, we need to compute the correction coefficient $a_1,\cdots ,a_m$, which can be written as:

$$\begin{array}{c}{\mathbf{y}}_q=-{\mathbf{Q}}_q^{+}{\mathbf{r}}_q,\hfill \\ a_q=1-{\displaystyle \sum _{i\ne q}}{a}_i,\hfill \end{array}$$

(14)

where $q\in \left\{1,\cdots ,m\right\}$. For illustrative purposes, we consider the q as $q=1$, then Equation (14) can be rewritten as:

$$\left[\begin{array}{ccc}{a}_2& \cdots & a_m\end{array}\right]=-{\left[\begin{array}{ccc}{\mathbf{r}}_2-{\mathbf{r}}_1& \cdots & {\mathbf{r}}_m-{\mathbf{r}}_1\end{array}\right]}^{+}{\mathbf{r}}_1.$$

(15)

Case 1: for $m>2$, it should be noted that we need to calculate a Moore–Penrose inverse to obtain a correction coefficient, the required number of complex-valued multiplication is $O\left(m{N}_t^2\right)$ [24].

Case 2: when $m=2$, the correction coefficient $a_1$ and $a_2$ can be rewritten as

$$a_2=-{\left({\mathbf{r}}_2-{\mathbf{r}}_1\right)}^{+}{\mathbf{r}}_1=-\frac{{\left({\mathbf{r}}_2-{\mathbf{r}}_1\right)}^H{\mathbf{r}}_1}{{∥{\mathbf{r}}_1-{\mathbf{r}}_2∥}_2^2},a_1=-\frac{{\left({\mathbf{r}}_1-{\mathbf{r}}_2\right)}^H{\mathbf{r}}_2}{{∥{\mathbf{r}}_1-{\mathbf{r}}_2∥}_2^2}.$$

(16)

The calculation of correction coefficient requires performing two vector multiplication and one norm calculation with the complexity is $O\left(N_t\right)$. The second case requires far fewer computations than the first case to realize coefficient calculation. Therefore, in terms of computational complexity, m = 2 is more advantageous for the iteration process than m = 3 and 4. Next, the detector performance is shown in terms of the bit error rate (BER) performance over a Rayleigh-fading channel to illustrate the effect of m value on performance. We consider an uplink massive MIMO scenario with 128 × 32 antenna configuration and 64-quadrature amplitude modulation (QAM) modulation scheme. In addition, the number of iterations is set as K = 3.

Figure 2 and

Table 1. Performance comparison for different m.

m Value	Signal-to-Noise Ratio (Target ${10}^{-5}$)	Computational Complexity of Moore–Penrose Inverse
m = 2	$11.56$ decibel (dB)	$O\left(N_t\right)$
m = 3	$11.44$ dB	$O\left(3N_t^2\right)$
m = 4	$11.55$ dB	$O\left(4N_t^2\right)$

show that the proposed algorithm with m = 3 and m = 4 does not result in an appreciable gain in performance with respect to the case m = 2. We can observe that the proposed algorithm with m = 3 only outperforms the proposed algorithm with m = 2 in 0.12 dB for BER = ${10}^{-5}$ and its complexity is $3N_t$ times higher. The proposed algorithm m = 4 even has the same BER of the proposed algorithm with m = 2. As a consequence, we choose the value of m is 2. In Algorithm 2, the pseudo code of the proposed algorithm with m = 2 is presented, where the algorithm terminates at the K-th iteration.

3.3. Convergence Proof

Taking the iterative process into account, the overall convergence performance is split into two parts for theoretical analysis. For the first part, we research the convergence performance of the combined SD and Jacobi iterative process and compare it with the conventional Jacobi algorithm.

Theorem 1.

A necessary and sufficient condition for the iterative algorithm ${\mathbf{x}}^{\left(k\right)}=\mathbf{B}{\mathbf{x}}^{\left(k-1\right)}+\mathbf{c}$ to be convergent for all initial vectors ${\mathbf{x}}^{\left(0\right)}$ is that $\rho \left(\mathbf{B}\right)<1$ [25].

As the theorem shows, the convergence velocity of the iterative algorithm is closely related to the spectral radius of the iteration matrix. Therefore, we derive the iteration matrix of the combined algorithm firstly. The iterative form of iterative process of the combined algorithm is concluded as below:

$$\begin{array}{cc}\hfill {\overline{\mathbf{x}}}_2^{\left(k\right)}& ={\overline{\mathbf{x}}}^{\left(k-1\right)}+{\alpha }_{k-1}{\mathbf{r}}^{\left(k-1\right)}+{\mathbf{D}}^{-1}\left({\mathbf{r}}^{\left(k-1\right)}-{\alpha }_{k-1}{\mathbf{p}}^{\left(k-1\right)}\right)\hfill \\ & =\left(\mathbf{I}-{\alpha }_{k-1}\mathbf{A}-{\mathbf{D}}^{-1}\mathbf{A}+{\alpha }_{k-1}{\mathbf{D}}^{-1}{\mathbf{A}}^2\right){\overline{\mathbf{x}}}^{\left(k-1\right)}+{\alpha }_{k-1}\mathbf{b}+{\mathbf{D}}^{-1}\mathbf{b}-{\alpha }_{k-1}{\mathbf{D}}^{-1}\mathbf{Ab},\hfill \end{array}$$

(17)

where $\mathbf{B}=\mathbf{I}-{\alpha }_{k-1}\mathbf{A}-{\mathbf{D}}^{-1}\mathbf{A}+{\alpha }_{k-1}{\mathbf{D}}^{-1}{\mathbf{A}}^2$ is considered as the iteration matrix. Afterward, analyze the spectral radius of $\mathbf{B}$, hereinafter $\rho \left(\mathbf{B}\right)$:

$$\begin{array}{cc}\hfill \rho \left(\mathbf{B}\right)& =\rho \left(\mathbf{I}-{\alpha }_{k-1}\mathbf{A}-{\mathbf{D}}^{-1}\mathbf{A}+{\alpha }_{k-1}{\mathbf{D}}^{-1}{\mathbf{A}}^2\right)\hfill \\ & =\rho \left(\left(\mathbf{I}-{\alpha }_{k-1}\mathbf{A}\right)\left(\mathbf{I}-{\mathbf{D}}^{-1}\mathbf{A}\right)\right).\hfill \end{array}$$

(18)

Since the MMSE filtering matrix $\mathbf{A}$ is a positive definite matrix, the diagonal matrix $\mathbf{D}$ of $\mathbf{A}$ can be considered as the positive definite. According to that, we can infer that the norm of $\rho \left(\mathbf{I}-{\mathbf{D}}^{-1}\mathbf{A}\right)$ and $\rho \left(\mathbf{I}-{\alpha }_{k-1}\mathbf{A}\right)$, $k=1:K$ are both less than one [26]. Furthermore, we can get $\rho \left(\mathbf{I}-{\mathbf{D}}^{-1}\mathbf{A}\right)<1$ and $\rho \left(\mathbf{I}-{\alpha }_{k-1}\mathbf{A}\right)<1$, respectively. Thus, it is easy to infer $\rho \left(\mathbf{B}\right)<1$. According to Theorem 1, we can assert that the iteration of the combined algorithm is convergent.

Algorithm 2: Proposed improved Jacobi-based algorithm with m = 2

Input: (1) $\mathbf{H}$, channel matrix

(2) $\mathbf{y}$, received signal matrix

(3) K, number of iterations

Output: detected signal ${\overline{\mathbf{x}}}^{\left(K\right)}$

// Initialization

1: $\mathbf{A}={\mathbf{H}}^{\mathrm{H}}\mathbf{H}+{\sigma }^2{\mathbf{I}}_{N_t},\mathbf{b}={\mathbf{H}}^H\mathbf{y}$;

2: ${\overline{\mathbf{x}}}^{\left(0\right)}=\left({\mathbf{D}}^{-1}-{\mathbf{D}}^{-1}\mathbf{E}{\mathbf{D}}^{-1}\right)\mathbf{b}$;

3: ${\mathbf{r}}^{\left(0\right)}=\mathbf{b}-\mathbf{A}{\overline{\mathbf{x}}}^{\left(0\right)}$;

// Iterative Process(repeat K times)

4: for $k=1$ to K (iteration times) do

5:  if k = 1; the value of ${\mathbf{r}}^{\left(0\right)}$ is obtained in the initialization

6:  else;${\mathbf{r}}^{\left(k-1\right)}=a_1{\mathbf{r}}_1+a_2{\mathbf{r}}_2$, the value of ${\mathbf{r}}_1,{\mathbf{r}}_2$ is obtained in the previous iteration

7:  end if

8:  ${\mathbf{p}}^{\left(k-1\right)}=\mathbf{A}{\mathbf{r}}^{\left(k-1\right)}$;

9:  ${\alpha }_{k-1}=\frac{\left[{\mathbf{r}}^{\left(k-1\right)},{\mathbf{r}}^{\left(k-1\right)}\right]}{\left[{\mathbf{p}}^{\left(k-1\right)},{\mathbf{r}}^{\left(k-1\right)}\right]}$;

10:  ${\overline{\mathbf{x}}}^{\left(k\right)}={\overline{\mathbf{x}}}^{\left(k-1\right)}+{\alpha }_{k-1}{\mathbf{r}}^{\left(k-1\right)}+{\mathbf{D}}^{-1}\left({\mathbf{r}}^{\left(k-1\right)}-{\alpha }_{k-1}{\mathbf{p}}^{\left(k-1\right)}\right)$;

11:  ${\overline{\mathbf{x}}}_1={\overline{\mathbf{x}}}^{\left(k-1\right)},{\overline{\mathbf{x}}}_2={\overline{\mathbf{x}}}^{\left(k\right)}$;

12:  ${\mathbf{r}}_1={\mathbf{r}}^{\left(k-1\right)},{\mathbf{r}}_2=\mathbf{b}-\mathbf{A}{\overline{\mathbf{x}}}_2$;

13:  $a_1=-\frac{{\left({\mathbf{r}}_1-{\mathbf{r}}_2\right)}^H{\mathbf{r}}_2}{{∥{\mathbf{r}}_1-{\mathbf{r}}_2∥}_2^2},a_2=1-a_1$;

14:  ${\overline{\mathbf{x}}}^{\left(k\right)}=a_1{\overline{\mathbf{x}}}_1+a_2{\overline{\mathbf{x}}}_2$;

15: end for

In addition, the spectral radius of the Jacobi iterative matrix is $\rho \left(\mathbf{I}-{\mathbf{D}}^{-1}\mathbf{A}\right)$. Comparing Label (18) with the spectral radius of the conventional Jacobi algorithm, we obtain the following inequality:

$$\begin{array}{c}\rho \left(\left(\mathbf{I}-{\alpha }_{k-1}\mathbf{A}\right)\left(\mathbf{I}-{\mathbf{D}}^{-1}\mathbf{A}\right)\right)<\rho \left(\mathbf{I}-{\mathbf{D}}^{-1}\mathbf{A}\right).\end{array}$$

(19)

It is worth pointing out that the smaller spectral radius of the iterative matrix, the faster the convergence of the scheme. Therefore, we can infer that the convergence velocity of the proposed algorithm is faster than the conventional Jacobi algorithm.

For the second part, we focus on the convergence of the whole correction method, that is, prove the convergence of the whole correction method [21]. Supposing that the residuals value correspond to ${\overline{\mathbf{x}}}_1,{\overline{\mathbf{x}}}_2$ is ${\mathbf{r}}_1=\mathbf{b}-\mathbf{A}{\overline{\mathbf{x}}}_1,{\mathbf{r}}_2=\mathbf{b}-\mathbf{A}{\overline{\mathbf{x}}}_2$; then, the residual of the correction solution is $\mathbf{r}=\mathbf{b}-\mathbf{A}\overline{\mathbf{x}}$. Assume the residual value ${\mathbf{r}}_1,{\mathbf{r}}_2,\mathbf{r}$ as the space vector which takes the origin of the coordinate as the starting point, and the ${\mathbf{r}}_1,{\mathbf{r}}_2,\mathbf{r}$ as the end point. Select $i,k\in \left\{1,2\right\}$ randomly, and $i\ne k$, for ease of presentation and description, let i = 1, k = 2, then we can get

$$\begin{array}{cc}\hfill {\left({\mathbf{r}}_1-{\mathbf{r}}_2\right)}^H\mathbf{r}& ={\left({\mathbf{r}}_1-{\mathbf{r}}_2\right)}^H\left({\mathbf{r}}_2+a_1\left({\mathbf{r}}_1-{\mathbf{r}}_2\right)\right)\hfill \\ & ={\left({\mathbf{r}}_1-{\mathbf{r}}_2\right)}^H\left({\mathbf{r}}_2+{\mathbf{Q}}_2{\mathbf{y}}_2\right)\hfill \\ & ={\left({\mathbf{r}}_1-{\mathbf{r}}_2\right)}^H\left({\mathbf{r}}_2-{\mathbf{Q}}_2{\mathbf{Q}}_2^{+}{\mathbf{r}}_2\right)\hfill \\ & =\left({\left({\mathbf{r}}_1-{\mathbf{r}}_2\right)}^H-{\left({\mathbf{r}}_1-{\mathbf{r}}_2\right)}^H{\mathbf{Q}}_2{\mathbf{Q}}_2^{+}\right){\mathbf{r}}_2\hfill \\ & =\left({\left({\mathbf{r}}_1-{\mathbf{r}}_2\right)}^H-{\left({\mathbf{Q}}_2{\mathbf{Q}}_2^{+}\left({\mathbf{r}}_1-{\mathbf{r}}_2\right)\right)}^H\right){\mathbf{r}}_2\hfill \\ & =0.\hfill \end{array}$$

(20)

Thus, we can conclude that $\mathbf{r}$ is perpendicular to the hyperplane determined by points ${\mathbf{r}}_1$, ${\mathbf{r}}_2$, which means

$${∥\mathbf{r}∥}_2\le {∥{\mathbf{r}}_1∥}_2,{∥\mathbf{r}∥}_2\le {∥{\mathbf{r}}_2∥}_2.$$

(21)

As a consequence, the proposed algorithm exhibits good convergence performance.

3.4. Computational Complexity Analysis

The computational complexity is one of the necessary aspects in measuring the performance of detection algorithms for massive MIMO systems. In this part, we compute the proposed algorithm complexity in each step and analyze the complexity of the proposed algorithm and the conventional algorithms. Since all the algorithms we mentioned need to calculate the matrix $\mathbf{A}={\mathbf{H}}^H\mathbf{H}+{\sigma }^2{\mathbf{I}}_{N_t}$ and the output matched filter $\mathbf{b}$, we only concentrate on the complexity of the remaining parts. In addition, the complexity is defined as the required complex multiplications.

(1) compute ${\overline{\mathbf{x}}}^{\left(0\right)}$: As the computation of ${\overline{\mathbf{x}}}^{\left(0\right)}$ involves a multiplication of a $N_t\times N_t$ diagonal matrix with a $N_t\times N_t$ off-diagonal matrix and a multiplication of a $N_t\times N_t$ off-diagonal matrix with a $N_t\times 1$ vector, the complexity of initial solution is $2N_t^2$.

(2) compute ${\mathbf{r}}^{\left(k\right)}$: The computation complexity of ${\mathbf{r}}^{\left(k\right)}$ should be evaluated by two scalar multiplication with a $N_t\times 1$ vector, yielding the total complexity of 2$N_t$.

(3) compute ${\overline{\mathbf{x}}}_2$: In Algorithm 2, the computation of ${\overline{\mathbf{x}}}_2$ includes a multiplication of a $N_t\times N_t$ matrix $\mathbf{A}$ with a $N_t\times 1$ vector $\mathbf{r}$ : $\mathbf{p}$, two inner products of two $N_t\times 1$ vectors: ${\alpha }_{k-1}$, two scalar multiplication with $N_t\times 1$ vectors and a $N_t\times N_t$ diagonal matrix with a $N_t\times 1$ vector. In conclusion, the complexity in this step is $N_t^2+5N_t$.

(4) compute $\overline{\mathbf{x}}$: having obtained the value of ${\overline{\mathbf{x}}}_1$, ${\overline{\mathbf{x}}}_2$, the computation of $\overline{\mathbf{x}}$ includes a multiplication of a $N_t\times N_t$ matrix $\mathbf{A}$ with a $N_t\times 1$ vector ${\overline{\mathbf{x}}}_2$: ${\mathbf{r}}_2$, one multiplication of two $N_t\times 1$ vectors: $a_1$ and two scalar multiplication with $N_t\times 1$ vectors. In this step, the total complexity is $N_t^2+5N_t$.

In summary, the total calculated quantity required by the proposed algorithm is $(2N_t^2+10N_t)K+2N_t^2$.

From

Table 2. Computational complexity comparison.

Detectors	Complexity (Iteration Times K)
MMSE	$\frac{5}{6}{N}_t^3+\frac{3}{4}{N}_t^2+\frac{4}{3}{N}_t$
GS [16]	$\left(K+5\right)N_t^2$
SSOR [17]	$\left(2N_t^2+2N_t\right)K+2N_t^2$
Jacobi [19]	$N_t^2(K+2)$
Proposed	$\left(2N_t^2+10N_t\right)K+2N_t^2$

and Figure 3, due to the number of iterations usually being relatively small, the proposed algorithm has approximately the same low computational complexity as the conventional GS and SSOR algorithms. Although the DJ algorithm requires the fewest calculation among these algorithms, it requires much larger iterations to achieve the same accuracy with the proposed algorithm. Additionally, for any number of iterations, the proposed algorithm can reduce the computational complexity of the MMSE detection algorithm from $O\left(N_t^3\right)$ to $O\left(N_t^2\right)$. This advantage of computational complexity is reflected distinctly, especially when the number of users is large.

4. Simulation Results

Numerical simulations are designed to assess the performance of the proposed improved Jacobi iteration algorithm. We provide the simulation results on the BER performance to compare the proposed algorithm with the conventional algorithms. The BER curve of the typical MMSE is used as the benchmark for comparison. The simulation environment is assumed to be the uplink massive MIMO system with the 64-QAM modulation scheme.

First, we show simulation results for a massive MIMO system in an uncorrelated Rayleigh fading channel. Figure 4 illustrates the BER performance over a range of signal-to-noise (SNR) between the DJ algorithm and the improved Jacobi-based algorithm, when the complex channels of size $N_r\times N_t=128\times 32$. In addition, the number of iteration is denoted by K. As we can see, with the number of iterations is increased, both algorithms can improve the BER performance. However, the DJ algorithm requires massive iterations to meet the BER performance requirement. In comparison, the proposed algorithm can achieve near-optimal performance with a small number of iterations. In addition, compared with the exact MMSE algorithm, the proposed approximated algorithm can achieve a satisfying performance with the relatively small iterations K (i.e., $K=3$). In addition, when K = 3, the BER performance of the MMSE algorithm and that of the proposed algorithm has a difference within 0.1 dB.

Additionally, the BER performance of the proposed algorithm is compared with the conventional detection algorithms GS and SSOR in Figure 5, when $N_r\times N_t=128\times 32$. It is clearly noticed that the BER performances of all algorithms are closer to that of MMSE algorithm when the number of iterations increases. Obviously, the proposed algorithm can achieve the exact performance of the MMSE algorithm when the number of iteration is 3, which requires the computation complexity is $O\left(8N_t^2\right).$ The conventional algorithm GS and SSOR algorithm require at least four iterations to achieve the near MMSE performance, which means that the required complexity is $O\left(9N_t^2\right)$ and $O\left(10N_t^2\right)$, respectively. As a consequence, the proposed Jacobi-based algorithm outperforms the conventional algorithms in massive MIMO systems.

Figure 6 reveals that the advantages of the proposed Jacobi-based algorithm come to an effect when the number of users increases. Note that we can set the number of the receiving antennas at BS as 128 and SNR = 10 dB is applied to the channel. As seen in Figure 6, when the number of single-antenna users increases, the performance of the all algorithms suffers a non-negligible degradation. The proposed algorithm still can achieve almost the same BER performance as the MMSE algorithm with small iterations (i.e., $K=3$), regardless of the number of users. In addition, the GS and SSOR algorithms require more iterations (i.e., $K=4$) to achieve the near BER performance of MMSE algorithm.

Finally, for realistic MIMO systems, the spatial correlation plays an important role in the BER performance. Thus, we consider the influence of channel correlation on the performance of BER. The exponential correlation model is formulated in [27], and r represents the correlation coefficient. Figure 7 shows that the convergence performance of all algorithms degrades with the serious channel correlation. Furthermore, it is worth noting that the proposed algorithm can still converge to the MMSE algorithm with a small number of iterations (i.e., $K=3$ when $r=0.1$, $K=4$, when $r=0.4$). When $r=0.1$, to ensure the performance of the final approximation, the required complexity of the proposed Jacobi-based algorithm, GS algorithm and SSOR algorithm is $O\left(8N_t^2\right)$ (i.e., $K=3$), $O\left(9N_t^2\right)$ (i.e., $K=4$) and $O\left(10N_t^2\right)$ (i.e., $K=4$), respectively. Similarly, when $r=0.4$, the required complexity of the proposed Jacobi-based algorithm, GS algorithm and SSOR algorithm is $O\left(10N_t^2\right)$ (i.e., $K=4$), $O\left(11N_t^2\right)$ (i.e., $K=6$) and $O\left(14N_t^2\right)$ (i.e., $K=6$), respectively. Consequently, the proposed algorithm can still enjoy lower complexity than the GS and SSOR algorithm to achieve the near exact performance of the MMSE algorithm with serious channel correlation.

5. Conclusions

In this paper, we propose an improved Jacobi iterative algorithm in signal detection in the massive MIMO system. The performance improvement has been achieved with the fast convergence of the iteration process, which is conducted by integrating an SD method and the whole-correction method into the traditional Jacobi iterative algorithm. We evaluate the performance of the proposed algorithm by a theoretical analysis and Monte Carlo simulations. Based on the theoretical analysis, the proposed algorithm has been proved to have good convergence and low computational complexity. Simulation results show that the proposed Jacobi-based algorithm performs better than the conventional algorithms. Especially, with the small number of iterations, the improved Jacobi-based algorithm performs almost the same as the MMSE algorithm in terms of BER, but achieving a much lower computationally complexity compared to one of the latter. Future work includes the joint channel estimation and signal detection based on the proposed algorithm in the massive MIMO system and the related hardware design for a real application system.