# Joint Newton Iteration and Neumann Series Method of Convergence-Accelerating Matrix Inversion Approximation in Linear Precoding for Massive MIMO Systems

Authors: Lin Shao; Yunxiao Zu
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Due to large numbers of antennas and users, matrix inversion is complicated in linear precoding techniques for massive MIMO systems. Several approximated matrix inversion methods, including the Neumann series, have been proposed to reduce the complexity. However, the Neumann series does not converge fast enough. In this paper, to speed up convergence, a new joint Newton iteration and Neumann series method is proposed, with the first iteration result of Newton iteration method being employed to reconstruct the Neumann series. Then, a high probability convergence condition is established, which can offer useful guidelines for practical massive MIMO systems. Finally, simulation examples are given to demonstrate that the new joint Newton iteration and Neumann series method has a faster convergence rate compared to the previous Neumann series, with almost no increase in complexity when the iteration number is greater than or equal to 2.


---

## Body

## 1. Introduction

Massive multiple-input multiple-output (MIMO) [1–3] is one of the promising technologies for the fifth-generation communication system. Compared to traditional MIMO, hundreds of antennas are equipped in base stations (BSs) in massive MIMO systems in order to achieve orders of magnitude increases in spectral and energy efficiency [4–6]. Unfortunately, the application of the massive MIMO system still faces several challenging problems in practice. For example, linear precoding techniques, such as zero-forcing (ZF) precoding, are always involved with complicated matrix inversion due to large numbers of BS antennas and users.In order to reduce the complexity, several approximated approaches have been recently proposed to avoid the exact matrix inversion. A low-complexity Richardson method has been proposed in [7], but it has an uncertain parameter which affects the convergence of the method. Conjugate gradient method was applied to reduce the complexity of data detection and precoding in the massive MIMO system with realistic antenna configurations [8]. However, many divisions are involved in the approach. Gauss-Seidel (GS) method was employed in [9] to speed up convergence. Unfortunately, the matrix inversion is not obtained directly, which will add extra complexity to the calculation of the transmitted signal after precoding. Newton iteration method converges fast and the complexity can be controlled just by the number of iterations [10]. However, a complex calculation is always required to get an initial input to ensure convergence [10]. Truncated Neumann series in [11, 12] was proposed to obtain near-optimal performance. It contains matrix multiplication and matrix addition, which are preferable in hardware implementation. Nevertheless, how to speed up the convergence rate of the Neumann series is still a problem.In this paper, we propose a new joint Newton iteration and Neumann series method, where Newton iteration method is utilized to provide an efficient searching direction for the Neumann series. Specifically, the first iteration result of Newton iteration method is employed to reconstruct the Neumann series expansion to accelerate convergence. Furthermore, a high probability convergence condition is derived to guarantee the convergence of the new approach. This condition is expected to contribute to the massive MIMO system in practice.The remainder of the paper is organized as follows. First, Section2 introduces the background including the system model, the Neumann series, and Newton iteration method. Then, the new joint Newton iteration and Neumann series method is proposed in Section 3. Simulation results are in Section 4. Finally, conclusions are drawn in Section 5.


## 2. Background

In this section, the system model, the Neumann series, and Newton iteration method are introduced.


### 2.1. System Model

Consider a downlink massive MIMO system withN antennas at the BS and K
≪
N single antenna users. The received K
×
1 vector y can be expressed as(1)
y
=
ρ
H
x
+
n
,where ρ is the signal-to-noise ratio (SNR) in the downlink, H
∈
C
K
×
N denotes the Rayleigh fading channel matrix whose entries are independent and identically distributed (i.i.d.) zero-mean unit-variance complex Gaussian variables, n
∈
C
K
×
1 is the additive white Gaussian noise vector, and x
∈
C
N
×
1 presents the transmitted signal vector after precoding and is obtained by(2)
x
=
P
s
,where s is K
×
1 vector of QAM symbols for transmission [13] and P
∈
C
N
×
K denotes the precoding matrix which specifically for ZF precoding can be represented as(3)
P
=
β
H
†
,where (4)
H
†
=
H
H
H
H
H
-
1
≜
H
H
W
-
1denotes the pseudoinverse of H and β is a scalar which is always chosen as(5)
β
=
K
tr
⁡
W
-
1
,where tr
⁡
W
-
1 denotes the trace of W
-
1. The main computation complexity for ZF precoding lies in the inversion of K
×
K matrix W, so several approximated approaches including the Neumann series have been investigated recently, which are shown in Section 1.


### 2.2. Neumann Series

Based on the Neumann series, the matrix inversion ofK
×
K matrix W is transformed to the sum of matrix polynomials as(6)
W
-
1
=
∑
n
=
0
+
∞
I
K
-
Θ
W
n
Θ
,where Θ is K
×
K matrix. Note that, according to [14], the premise condition of (6) is(7)
l
i
m
n
→
∞
I
K
-
Θ
W
n
=
0
K
.For practical use, W
-
1 is approximated as [14](8)
W
-
1
≈
∑
n
=
0
L
I
K
-
Θ
W
n
Θ
,where L is the iteration number.Remark 1.
Based on (2), (3), and (4), x is obtained as(9)
x
=
β
H
H
W
-
1
s
.Let f
≜
W
-
1
s and f can be regarded as the solution of the linear equation(10)
W
f
=
s
.Recently, Richardson method [7], conjugate gradient method [8], and GS method [9] are all applied to solve (10) in an iterative way. But they approximate W
-
1
s instead of W
-
1, while the Neumann series approximates W
-
1 directly.


### 2.3. Newton Iteration Method

Newton iteration method can be employed to estimateW
-
1 in an iterative way. Suppose Z
0 is the original estimation of W
-
1 and the following condition is satisfied:(11)
I
K
-
W
Z
0
<
1
.Then, the nth iteration estimation of W
-
1 is expressed as(12)
Z
n
=
Z
n
-
1
2
I
K
-
W
Z
n
-
1
.Newton iteration method converges fast, but formula (11) usually needs a complex calculation to get an appropriate Z
0 [10].


## 2.1. System Model

Consider a downlink massive MIMO system withN antennas at the BS and K
≪
N single antenna users. The received K
×
1 vector y can be expressed as(1)
y
=
ρ
H
x
+
n
,where ρ is the signal-to-noise ratio (SNR) in the downlink, H
∈
C
K
×
N denotes the Rayleigh fading channel matrix whose entries are independent and identically distributed (i.i.d.) zero-mean unit-variance complex Gaussian variables, n
∈
C
K
×
1 is the additive white Gaussian noise vector, and x
∈
C
N
×
1 presents the transmitted signal vector after precoding and is obtained by(2)
x
=
P
s
,where s is K
×
1 vector of QAM symbols for transmission [13] and P
∈
C
N
×
K denotes the precoding matrix which specifically for ZF precoding can be represented as(3)
P
=
β
H
†
,where (4)
H
†
=
H
H
H
H
H
-
1
≜
H
H
W
-
1denotes the pseudoinverse of H and β is a scalar which is always chosen as(5)
β
=
K
tr
⁡
W
-
1
,where tr
⁡
W
-
1 denotes the trace of W
-
1. The main computation complexity for ZF precoding lies in the inversion of K
×
K matrix W, so several approximated approaches including the Neumann series have been investigated recently, which are shown in Section 1.


## 2.2. Neumann Series

Based on the Neumann series, the matrix inversion ofK
×
K matrix W is transformed to the sum of matrix polynomials as(6)
W
-
1
=
∑
n
=
0
+
∞
I
K
-
Θ
W
n
Θ
,where Θ is K
×
K matrix. Note that, according to [14], the premise condition of (6) is(7)
l
i
m
n
→
∞
I
K
-
Θ
W
n
=
0
K
.For practical use, W
-
1 is approximated as [14](8)
W
-
1
≈
∑
n
=
0
L
I
K
-
Θ
W
n
Θ
,where L is the iteration number.Remark 1.
Based on (2), (3), and (4), x is obtained as(9)
x
=
β
H
H
W
-
1
s
.Let f
≜
W
-
1
s and f can be regarded as the solution of the linear equation(10)
W
f
=
s
.Recently, Richardson method [7], conjugate gradient method [8], and GS method [9] are all applied to solve (10) in an iterative way. But they approximate W
-
1
s instead of W
-
1, while the Neumann series approximates W
-
1 directly.


## 2.3. Newton Iteration Method

Newton iteration method can be employed to estimateW
-
1 in an iterative way. Suppose Z
0 is the original estimation of W
-
1 and the following condition is satisfied:(11)
I
K
-
W
Z
0
<
1
.Then, the nth iteration estimation of W
-
1 is expressed as(12)
Z
n
=
Z
n
-
1
2
I
K
-
W
Z
n
-
1
.Newton iteration method converges fast, but formula (11) usually needs a complex calculation to get an appropriate Z
0 [10].


## 3. Joint Newton Iteration and Neumann Series Method

It can be seen from (7) that, as the first item of the series, Θ greatly affects convergence. How to choose Θ plays a key role in the Neumann series. In [13], Θ is initially set as (
1
/
(
N
+
K
)
)
I
K and the Neumann series converges when N and K grow to infinity. In [11], Θ is then chosen as the matrix inversion of K
×
K diagonal matrix D whose entries are the main diagonal elements of W, resulting in a faster convergence rate. In this paper, Newton iteration method is employed to provide an appropriate Θ to speed up the convergence of the Neumann series. Moreover, a high probability convergence condition about N
/
K ratio is derived for the joint Newton iteration and Neumann series method. Finally, the complexity of the new approach is analyzed.


### 3.1. Joint Newton Iteration and Neumann Series Method and Its Convergence Condition

The fast convergence property of Newton iteration method inspired us to useZ
1 instead of D
-
1 as Θ to speed up the convergence rate of the Neumann series. Note that formula (11) must be satisfied before the application of Newton iteration method. According to [10], setting Z
0 to be D
-
1 makes it easy to meet (11) for massive MIMO systems. Therefore, we get the joint Newton iteration and Neumann series method as follows.First, letZ
0
=
D
-
1 and then(13)
Z
1
=
Z
0
2
I
K
-
W
Z
0
=
D
-
1
2
I
K
-
W
D
-
1
.Second, useZ
1 to reconstruct the Neumann series as follows:(14)
W
-
1
≈
∑
n
=
0
L
I
K
-
Z
1
W
n
Z
1
.It is worth pointing out that, in order to guarantee the convergence of the Neumann series withΘ
=
Z
1, the accurate lower bound of N
/
K ratio needs to be determined. For this end, a high probability convergence condition about N
/
K ratio is given in the following lemma.Lemma 2.
For downlink massive MIMO systems, the Neumann series withΘ
=
Z
1 converges with high probability when(15)
α
≜
N
K
>
1
2
-
1
2
.Proof.
∑
n
=
0
∞
I
K
-
Z
1
W
n
Z
1 converges:⇔
lim
n
→
∞
⁡
I
K
-
Z
1
W
n
=
0
K,
⇔
ρ
I
K
-
Z
1
W
<
1, where ρ
I
K
-
Z
1
W denotes the spectral radius of I
K
-
Z
1
W,
⇔
λ
A
<
1, where λ
A is any eigenvalue of A and A
=
I
K
-
Z
1
W.Note that(16)
A
=
I
K
-
Z
1
W
=
I
K
-
D
-
1
2
I
K
-
W
D
-
1
W
=
I
K
-
2
D
-
1
W
+
D
-
1
W
2
=
I
K
-
D
-
1
W
2and ∑
n
=
0
∞
I
K
-
D
-
1
W
n
D
-
1 converges ⇔
λ
B
<
1, where B
=
I
K
-
D
-
1
W.
SinceA
=
B
2 and λ
A
=
λ
B
2, ∑
n
=
0
∞
I
K
-
Z
1
W
n
Z
1 converges ⇔
∑
n
=
0
∞
I
K
-
D
-
1
W
n
D
-
1 converges.
According to [14], for ∑
n
=
0
∞
I
K
-
D
-
1
W
n
D
-
1, a high probability convergence condition in terms of α is (17)
α
>
1
2
-
1
2
.Thus, the Neumann series with Θ
=
Z
1, that is, ∑
n
=
0
∞
I
K
-
Z
1
W
n
Z
1, converges with a high probability when α
>
1
/
2
-
1
2.By (15), the maximum value of K can be calculated for a specific N to achieve a high convergence probability for the Neumann series with Θ
=
Z
1. Figure 1 shows the convergence probability with the configurations in Table 1.Table 1
MaximumK values corresponding to different BS antenna numbers under the condition of (15).
N
64
128
192
256
320
384
448
512
K
10
21
32
43
54
65
76
87Figure 1
Convergence probability of Neumann series withΘ
=
Z
1 under configurations in Table 1.As seen in Figure1, the joint Newton iteration and Neumann series method can achieve a very high probability of convergence under the condition of (15), with the convergence probability as high as 0.999. Furthermore, when N is relatively small, the performance is slightly better. Also, as Figure 1 is obtained with the smallest allowed α under (15), it can be deduced that the convergence probability will be higher as α increases.Lemma2 has potential applications in massive MIMO systems. For the two typical downlink massive MIMO configurations N
×
K
=
256
×
16 and N
×
K
=
256
×
32 with α
=
16 and α
=
8, respectively [9], by Lemma 2, it can be concluded that the new joint Newton iteration and Neumann series method is convergent in both of these scenarios.


### 3.2. Complexity Analysis

The number of complex-valued multiplications is employed as a roughly estimated complexity of an algorithm. WhenΘ
=
D
-
1 or Θ
=
(
1
/
(
N
+
K
)
)
I
K, the complexity of (8) is O
K
2 for L
=
1 and O
K
3 for L
≥
2. However, when Θ
=
Z
1, the complexity is O
K
3 for L
≥
1. As a relatively large number of iterations (e.g., L
≥
4) are usually needed to avoid too much performance loss, the complexity of the three approaches is comparable. Note that, for different Θ and L
≥
2, the complexity of the Neumann series is O
K
3 which is comparable to that of the exact matrix inversion [10]. However, the Neumann series only contains matrix addition and matrix multiplication which is strongly preferred over matrix inversion in hardware since it does not require any divisions [13].


### 3.3. Discussion

Compared with the new joint Newton iteration and Neumann series method, conjugate gradient method [8] and GS method [9] converge faster, but they approximate W
-
1
s instead of W
-
1, as stated in Remark 1. However, W
-
1 is preferred in the computation of β and other related calculations, such as the sum rate computation of ZF precoding [9] and fast matrix inversion updates [15]. In comparison, the new joint Newton iteration and Neumann series method approximates W
-
1 directly and thus can significantly reduce the total complexity in those calculations. Furthermore, it only involves matrix multiplication and matrix addition, while conjugate gradient method contains many divisions which are more difficult for hardware implementation. Richardson method [7] was applied in minimum mean square error (MMSE) signal detection, but the relaxation parameter of the method still remains unknown for ZF precoding. In other words, Richardson method [7] may not converge for ZF precoding. By contrast, the new joint Newton iteration and Neumann series method will have a very high probability of convergence if (15) is satisfied.


## 3.1. Joint Newton Iteration and Neumann Series Method and Its Convergence Condition

The fast convergence property of Newton iteration method inspired us to useZ
1 instead of D
-
1 as Θ to speed up the convergence rate of the Neumann series. Note that formula (11) must be satisfied before the application of Newton iteration method. According to [10], setting Z
0 to be D
-
1 makes it easy to meet (11) for massive MIMO systems. Therefore, we get the joint Newton iteration and Neumann series method as follows.First, letZ
0
=
D
-
1 and then(13)
Z
1
=
Z
0
2
I
K
-
W
Z
0
=
D
-
1
2
I
K
-
W
D
-
1
.Second, useZ
1 to reconstruct the Neumann series as follows:(14)
W
-
1
≈
∑
n
=
0
L
I
K
-
Z
1
W
n
Z
1
.It is worth pointing out that, in order to guarantee the convergence of the Neumann series withΘ
=
Z
1, the accurate lower bound of N
/
K ratio needs to be determined. For this end, a high probability convergence condition about N
/
K ratio is given in the following lemma.Lemma 2.
For downlink massive MIMO systems, the Neumann series withΘ
=
Z
1 converges with high probability when(15)
α
≜
N
K
>
1
2
-
1
2
.Proof.
∑
n
=
0
∞
I
K
-
Z
1
W
n
Z
1 converges:⇔
lim
n
→
∞
⁡
I
K
-
Z
1
W
n
=
0
K,
⇔
ρ
I
K
-
Z
1
W
<
1, where ρ
I
K
-
Z
1
W denotes the spectral radius of I
K
-
Z
1
W,
⇔
λ
A
<
1, where λ
A is any eigenvalue of A and A
=
I
K
-
Z
1
W.Note that(16)
A
=
I
K
-
Z
1
W
=
I
K
-
D
-
1
2
I
K
-
W
D
-
1
W
=
I
K
-
2
D
-
1
W
+
D
-
1
W
2
=
I
K
-
D
-
1
W
2and ∑
n
=
0
∞
I
K
-
D
-
1
W
n
D
-
1 converges ⇔
λ
B
<
1, where B
=
I
K
-
D
-
1
W.
SinceA
=
B
2 and λ
A
=
λ
B
2, ∑
n
=
0
∞
I
K
-
Z
1
W
n
Z
1 converges ⇔
∑
n
=
0
∞
I
K
-
D
-
1
W
n
D
-
1 converges.
According to [14], for ∑
n
=
0
∞
I
K
-
D
-
1
W
n
D
-
1, a high probability convergence condition in terms of α is (17)
α
>
1
2
-
1
2
.Thus, the Neumann series with Θ
=
Z
1, that is, ∑
n
=
0
∞
I
K
-
Z
1
W
n
Z
1, converges with a high probability when α
>
1
/
2
-
1
2.By (15), the maximum value of K can be calculated for a specific N to achieve a high convergence probability for the Neumann series with Θ
=
Z
1. Figure 1 shows the convergence probability with the configurations in Table 1.Table 1
MaximumK values corresponding to different BS antenna numbers under the condition of (15).
N
64
128
192
256
320
384
448
512
K
10
21
32
43
54
65
76
87Figure 1
Convergence probability of Neumann series withΘ
=
Z
1 under configurations in Table 1.As seen in Figure1, the joint Newton iteration and Neumann series method can achieve a very high probability of convergence under the condition of (15), with the convergence probability as high as 0.999. Furthermore, when N is relatively small, the performance is slightly better. Also, as Figure 1 is obtained with the smallest allowed α under (15), it can be deduced that the convergence probability will be higher as α increases.Lemma2 has potential applications in massive MIMO systems. For the two typical downlink massive MIMO configurations N
×
K
=
256
×
16 and N
×
K
=
256
×
32 with α
=
16 and α
=
8, respectively [9], by Lemma 2, it can be concluded that the new joint Newton iteration and Neumann series method is convergent in both of these scenarios.


## 3.2. Complexity Analysis

The number of complex-valued multiplications is employed as a roughly estimated complexity of an algorithm. WhenΘ
=
D
-
1 or Θ
=
(
1
/
(
N
+
K
)
)
I
K, the complexity of (8) is O
K
2 for L
=
1 and O
K
3 for L
≥
2. However, when Θ
=
Z
1, the complexity is O
K
3 for L
≥
1. As a relatively large number of iterations (e.g., L
≥
4) are usually needed to avoid too much performance loss, the complexity of the three approaches is comparable. Note that, for different Θ and L
≥
2, the complexity of the Neumann series is O
K
3 which is comparable to that of the exact matrix inversion [10]. However, the Neumann series only contains matrix addition and matrix multiplication which is strongly preferred over matrix inversion in hardware since it does not require any divisions [13].


## 3.3. Discussion

Compared with the new joint Newton iteration and Neumann series method, conjugate gradient method [8] and GS method [9] converge faster, but they approximate W
-
1
s instead of W
-
1, as stated in Remark 1. However, W
-
1 is preferred in the computation of β and other related calculations, such as the sum rate computation of ZF precoding [9] and fast matrix inversion updates [15]. In comparison, the new joint Newton iteration and Neumann series method approximates W
-
1 directly and thus can significantly reduce the total complexity in those calculations. Furthermore, it only involves matrix multiplication and matrix addition, while conjugate gradient method contains many divisions which are more difficult for hardware implementation. Richardson method [7] was applied in minimum mean square error (MMSE) signal detection, but the relaxation parameter of the method still remains unknown for ZF precoding. In other words, Richardson method [7] may not converge for ZF precoding. By contrast, the new joint Newton iteration and Neumann series method will have a very high probability of convergence if (15) is satisfied.


## 4. Numerical Results

The bit error rate (BER) performances of different algorithms are evaluated in order to compare their convergence rates. For simplicity, algorithm1, algorithm 2, and algorithm 3 represent the Neumann series with Θ
=
D
-
1, Θ
=
Z
1, and Θ
=
(
1
/
N
+
K
)
I
K, respectively. Moreover, ZF precoding with exact matrix inversion of W is also included as the benchmark. Newton iteration method usually needs complex calculation for initial estimation, which limits its range of application [10]. Therefore, Newton iteration method is not involved in comparison. The typical downlink massive MIMO configuration with N
×
K
=
256
×
32 is considered and 64 QAM is employed as the modulation scheme.Figures2 and 3 show the BER performance comparison between the Neumann series with different Θ. It is obvious in Figure 2 that the BER performance of both the algorithms improves with the iteration number L. However, for a given number of iterations, algorithm 2 achieves a better BER performance than algorithm 1. Therefore, algorithm 2 has a faster convergence rate than algorithm 1. Moreover, when the iteration number is relatively large (e.g., L
=
4), the BER of algorithm 2 is close to that of the exact matrix inversion, while algorithm 1 still suffers a great performance loss in comparison. In Figure 3, the advantage of algorithm 2 in convergence rate becomes more obvious. Note that the BER performance of algorithm 3 is even worse when L
=
4 than when L
=
3, which implies that algorithm 3 may not converge when α
=
8. Similarly, in Figure 4, Richardson method shows signs of nonconvergence, which indicates that the relaxation parameter configuration in [7] may not suit ZF precoding. From Figures 2, 3, and 4, it can be concluded that algorithm 2 has the fastest convergence rate.Figure 2
BER performance comparison between algorithm1 and algorithm 2 in N
×
K
=
256
×
32 massive MIMO system.Figure 3
BER performance comparison between algorithm2 and algorithm 3 in N
×
K
=
256
×
32 massive MIMO system.Figure 4
BER performance comparison between algorithm2 and Richardson method in N
×
K
=
256
×
32 massive MIMO system.


## 5. Conclusions

In this paper, a new joint Newton iteration and Neumann series method has been studied for matrix inversion computation involved in linear precoding techniques. Newton iteration method was employed to choose the initial value for the Neumann series. A high probability convergence condition was derived to ensure the convergence of the new method. Simulation results were provided to illustrate that the new matrix inversion computing method converges faster than some old ones.



---

*Source: 1745808-2016-05-31.xml*