1

Low-Complexity Statistically Robust
Precoder/Detector Computation for Massive

arXiv:1711.11405v2 [cs.IT] 21 Dec 2017

MIMO Systems
Mahdi Nouri Boroujerdi∗ , Saeid Haghighatshoar† , Member, IEEE, Giuseppe
Caire†, Fellow, IEEE
∗

School of Electrical and Computer Engineering, University of Tehran, Tehran, Iran, Email: {mhd.nouri}@ut.ac.ir
†

Communications and Information Theory Group, Technische Universität Berlin
Email: {saeid.haghighatshoar, caire}@tu-berlin.de

Abstract
Massive MIMO is a variant of multiuser MIMO in which the number of antennas at the base station
(BS) M is very large and typically much larger than the number of served users (data streams) K. Recent
research has widely investigated the system-level advantages of massive MIMO and, in particular, the
beneficial effect of increasing the number of antennas M . These benefits, however, come at the cost of
a dramatic increase in hardware and computational complexity. This is partly due to the fact that the BS
needs to compute precoding/receiving vectors in order to coherently transmit/detect data to/from each
user, where the resulting complexity grows proportionally to the number of antennas M and the number
of served users K. Recently, different algorithms based on tools from asymptotic random matrix theory
and/or approximated message passing have been proposed to reduce such complexity. The underlying
assumption in all these techniques, however, is that the exact statistics (covariance matrix) of the channel
vectors of the users is a priori known. This is far from being realistic, especially taking into account
that, in the high-dim regime of M ≫ 1, the estimation of the user channel covariance matrices is also
challenging both in terms of computation and storage requirements. In this paper, we propose a novel
technique for computing the precoder/detector in a massive MIMO system. Our method is based on the
randomized Kaczmarz algorithm and does not require a priori knowledge of the statistics of the users’
channel vectors. We analyze the performance of our proposed algorithm theoretically and compare its
performance with that of other techniques based on random matrix theory and approximate message
passing via numerical simulations. Our results indicate that the proposed technique has a comparable

2

performance while being highly robust to variations in the channel statistics and being computationally
very competitive.

Index Terms
Massive MIMO, Kaczmarz Algorithm, Detection, Precoding.

I. I NTRODUCTION
We consider a massive MIMO system where a base station (BS) is equipped with an array
consisting of M antennas and serves K users simultaneously, using spatial multiplexing. In
such a system, the spatial multiplexing gain of K achieves high sum spectral efficiency, while
the large antenna diversity achieved by using M ≫ K antennas provides channel hardening,
such that also the per-user rates are very stable and robust against statistical variations of the
propagation environment [2, 2–5]. Using large number of antennas at the BS has many other
interesting advantages, for instance, as the number of antennas tends to infinity, i.e. M → ∞,
the intra cellular interference vanishes so does the effect of noise and fast (small-scale) fading
[3]. Moreover, due to channel hardening and to the fact that for M much larger than K the user
channel vectors tend to become mutually orthogonal, linear precoding/detection in the Uplink
(UL) and the Downlink (DL) is sufficient to achieve a performance comparable to that of the
optimal nonlinear transmitter/receiver [3, 6–8].
In this paper, we focus on linear schemes such as zero forcing (ZF) precoder, also known as
decorrelator (when used as BS receiver in the UL [9–11]), regularized ZF (RZF), where the
regularization factor is tuned to the level of the noise at the user receivers, and linear minimum
mean squared error (MMSE) detector, which achieves the optimal Signal-to-Interference plus
Noise Ratio (SINR) among all linear receivers [12] (we refer to [13] for a more detailed review
of the available techniques for MIMO precoding and detection). Computing the corresponding
linear precoding/detection matrices from the estimated user channel vectors requires matrix
inversion and matrix-matrix multiplication, which consists of O(K 3 ) and O(MK 2 ) operations,
respectively. This imposes a very high computational complexity to massive MIMO implementations
since the number of BS antennas M and the number of served users K can be very large.
Recently a significant research effort has been devoted to develop efficient and low-complexity

3

algorithms for designing precoding/detection matrices in massive MIMO. Most of the works
exploit underlying statistical properties of the channel states of the users and asymptotic results
in random matrix theory to reduce the complexity. An important class of these algorithm use the
truncated polynomial expansion (TPE) technique [14–19] in order to approximate the inverse
matrix via a Taylor series expansion with a small number of terms. This reduces the complexity
by a factor growing proportionally to the number of antennas M (or users K). TPE-based
techniques, however, require pre-calculation of the polynomial coefficients based on the statistics
of the channel vectors, thus, require statistical information of the user channels. In particular,
as seen in [17–19], they require the knowledge of each user channel covariance matrix. This
may not be practical, especially in massive MIMO regime, where the covariance matrix of each
user channel has very large dimension M × M and needs to be estimated and tracked (i.e.,
updated over time) from the user channel measurements. Estimation and tracking of the user
channel covariance matrix in the massive MIMO regime has also been intensely investigated in
recent works (e.g., see [20, 21]) and it is known to be quite a challenging and computationally
demanding problem. A similar approach for linear MMSE detection consists of using a reduced
rank filter, as in [22]. In a recent line of work [23], the approximate message passing (AMP)
algorithm, initially developed for sparse signal recovery in [24], has been used to reduce the
computational complexity of massive MIMO precoding and detection.
In this paper, we propose a novel technique for designing precoding/detection matrices in a
massive MIMO system. Our method is based on the randomized Kaczmarz Algorithm (KA)
and does not require any knowledge of the statistics of the users’ channels. KA was initially
proposed by Kaczmarz as an iterative technique for solving over-determined (OD) set of linear
equations (SLE) [25]. A randomized version of KA was recently proposed and analyzed for
solving consistent OD SLE [26]. With the advent of stochastic gradient techniques in machine
learning, KA has been revitalized [27] and applied to other problems such as solving quadratic
equations [28]. Our work in this paper is mainly based on the randomized KA proposed in
[26]. We extend this algorithm to the under-determined (UD) SLE and use it to design new
variants of KA for massive MIMO applications. We also analyze the performance of our proposed
algorithm theoretically and compare its performance with that of other relevant techniques based
on random matrix theory and approximate message passing via numerical simulations. Our results
indicate that the proposed scheme has a comparable performance and is very competitive from

4

a computational viewpoint, while being highly robust to variations in channel statistics.
A. Notation
We represent scalar constants by non-boldface letters (e.g., x or X), sets by calligraphic letters
(e.g., X ), vectors by boldface small letters (e.g., x), and matrices by boldface capital letters (e.g.,
X). We denote the i-th row and the j-th column of a matrix X with the row-vector Xi,. and the
column-vector X.,j respectively. We represent the components of a vector x = (x1 , . . . , xn )T by
xi , i = 1, 2, . . . , n. We denote a sequence of vectors by upper/lower indices, e.g., x1 , x2 , · · · or

x1 , x2 , · · · . We indicate the Hermitian conjugate and the transpose of a matrix X by XH and

XT respectively, with the same notation being used for vectors and scalars. We denote the inner

product between two matrices (and similarly two vectors) X and Y by hX, Yi = tr(XH Y). We
use kXkF = hX, Xi for the Frobenius norm of a matrix X and kxk for the l2 -norm of a vector

x. For two m × n and m′ × n matrices X and Y, we denote by [X; Y] the (m + m′ ) × n matrix
obtained by stacking the rows of X on top of the rows of Y. The identity matrix of order p is
represented by Ip . For an integer k > 0, we use the shorthand notation [k] for {1, 2, . . . , k}.
II. P ROBLEM S TATEMENT
Consider a multiuser massive MIMO system with a BS with M ≫ 1 antennas serving K
single-antenna users. We denote the true channel vector of user i ∈ [K] to the M BS antennas

with hi ∈ CM , where hi = (hi1 , . . . , hiM )T with hik denoting the complex channel gain to the

k-th antenna. We denote the whole true M × K channel matrix by H = [h1 , . . . , hK ], consisting
of the channel vectors of all the users arranged by columns. We also denote an estimate of this
matrix available at the base station (BS) by the M × K matrix Q. It is worthwile to clarify that
all the computation done in the BS are only based on the estimate matrix Q. In Time-Division
Duplexing (TDD) systems, matrix Q is obtained during the training period of each channel
coherence block1 by transmitting orthogonal pilots from the users in the UL (exploiting UL-DL
channel reciprocity [3, 29]). In Frequency-Division Duplexing (FDD) systems, the BS obtains
1

The channel coherence block is the region in the time-frequency plane where the small-scale fading process generating the

coefficients of H can be considered constant. For a channel with physical coherence time ∆Tc and coherence bandwidth ∆Wc ,
the channel coherence block spans approximately ⌈∆Tc × ∆Wc ⌉ complex signal dimensions.

5

the estimated channel matrix by DL pilot transmission and some form of analog or quantized
feedback from the users to the BS in the UL (e.g., see [30]).
A. Detection and precoding vectors in the UL and the DL
In the UL, all the K users transmit their data symbols to the BS. The M-dim received vector
at the BS is given by:
y = Hs + n,

(1)

where s = (s1 , . . . , sK )T ∈ CK is the vector of symbols sent by the K users and n ∼

CN (0, σ 2 IM ) is the additive zero-mean white Gaussian noise (AWGN) at the BS receiver.
Without loss of generality we assume that the user symbols have mean zero and variance
E[|si |2 ] = P , and incorporate the pathloss effects due to possibly different distance of the users

form the BS as part of the channel vectors {hk }. As said before, we consider the case where
the BS detects symbols s from the observations y via a linear receiver. Two widely used linear
detectors are zero forcing detector (ZFD) and minimum mean squared error detector (MMSED)
which yield the following symbol estimates:
ŝZFD = (QH Q)−1 QH y,
ŝMMSED = (QH Q +
where SNR =

P
σ2

1
IK )−1 QH y,
SNR

(2)
(3)

denotes the per-user transmit Signal-to-Noise Ratio (SNR). In the case of

massive MIMO, where K is very large (e.g., some tens) and M is even much larger (e.g., some
hundreds), even implementing these simple linear detectors is computationally challenging since
they require large-dim matrix-matrix multiplication, which needs roughly O (MK 2 ) operations,
and matrix inversion which requires approximately O (K 3 ) operations.
A similar computational issue arises in the DL, where the BS transmits data to the users by
coherently beamforming data to the users via suitable linear precoding matrices. We consider
linear precoding schemes such as zero forcing beamforming (ZFBF) and regularized zero forcing
beamforming (RZFBF), which have been shown to achieve near-optimal performance in massive
MIMO [6]. The precoding vectors in the DL can be computed according to [31, 32] as:
xZFBF = βZFBF Q(QH Q)−1 s
xRZFBF = βRZFBF Q(QH Q + ξI)−1s,

(4)
(5)

6

where βZFBF and βRZFBF are scaling parameters to impose a given BS average transmit power
and ξ is the regularization parameter [32]. Each user k ∈ [K] will then receive the following
signal
H

yk = hk x + nk ,

(6)

where x can either be xZF or xRZFBF and nk ∼ CN (0, σ 2 ) is the AWGN at the user. Computing
precoding vectors of the ZFBF and RZFBF in (4) and (5) has the same order of complexity as
in the corresponding UL scenario (2) and (3).
B. Low-Complexity Beamforming via Kaczmarz Algorithm
Our goal in this paper is to reduce the complexity of the detector/precoder computation in
the massive MIMO. We first show that detecting the users’ signals in the UL, i.e., ŝZFD , ŝMMSED
in (2) and (3), or finding the suitable M-dim precoded signal to be the transmitted to the users
in the DL, i.e., xZFBF , xRZFBF in (4) and (5), can be posed as finding the solution w∗ of an SLE
of the form Aw = b where A depends on the estimated channel state Q and b depends on
the users’ signal s in a DL and on the noisy received signal from the users in a UL scenario.
We find the optimal solution w∗ of Aw = b via an iterative procedure known as Kaczmarz
Algorithm (KA). This technique was initially proposed by Kaczmarz for efficiently solving a
consistent overdetermined (OD) SLE [25].
KA works as follows. At the start of each iteration t, the algorithm has an estimate wt of the
optimal solution w∗ and selects one of the rows, say row number r(t), of A. If the estimate
wt satisfies the r(t)-th equation, namely, har(t) , wti = br(t) , where ar(t) = (Ar(t),. )H denotes
the conjugate of the r(t)-th row of A, then KA keeps wt . Otherwise, it updates wt along ar(t)

to make the r(t)-th equation consistent. This is summarized in Algorithm 1. Further insight
into the performance of KA is gained by noting that each Kaczmarz update indeed obeys the
minimum energy perturbation principle: At each iteration t, KA finds the closest vector w (in
l2 -norm) to the current solution wt that satisfies the selected equation or mathematically speaking
wt+1 = arg minw kw − wt k2 s.t. har(t) , wi = br(t) .
The convergence speed of KA is known to depend highly on the order of selecting the
rows of A, referred to as the update schedule. In particular, the algorithm can be quite slow
under a round-robin schedule that selects the rows periodically. This has motivated the study
of randomized variants of KA in which the rows of A are selected randomly according to

7

Algorithm 1 Kaczmarz Algorithm
1:

Initialize w0 .

2:

for t = 0, 1, . . . , T − 1 do

3:

Pick a row r(t) of A denoted by the row vector aHr(t) at each iteration t. In [26], row i

is selected randomly with a probability proportional to its l2 -norm squared kai k2 .
4:

Update wt as wt+1 = wt +

br(t) −har(t) ,wt i
ar(t)
kar(t) k2

such that after the update, har(t) , wt+1 i =

br(t) and r(t)-th equation is fulfilled.
5:

end for

some probability distribution that might depend on A. Such a randomized version of KA was
recently proposed and analyzed by Strohmer and Vershynin in [26], where they showed that the
performance of the resulting randomized KA depends on the condition number of the matrix
A. It was also shown in [26] that, for a wide class of matrices, the proposed algorithm can
potentially outperform other well-known techniques such as conjugate gradient.
The traditional KA and also the randomized variant proposed in [26] are applicable to OD
but consistent cases, where there is a solution satisfying all the equations. For massive MIMO
applications considered in this paper, in contrast, we will need variants of KA for both over-determined
and under-determined scenarios. Moreover, in the over-determined cases the equations are almost
always inconsistent. For example, as we will see this inconsistency arises in a UL scenario due
to the presence of noise in the received signal at the BS. In fact, massive MIMO linear detection
in such cases requires finding the least-squares solution of an appropriate SLE. However, an
ad-hoc application of KA results in a residual term and does not yield the desired solution [33].
We avoid this residual effect in our case by developing a new variant of KA that always solves
consistent equations after suitable transformation, where we also show that the resulting updates
converge to the desired least-squares solution. Our method resembles [34], which uses a step
prior to KA to remove the inconsistency, but has much less computational complexity.
III. U NIFIED M ATHEMATICAL F RAMEWORK
In this section, we propose a unified approach that we will use to analyze the performance of
all variants of KA studied in the sequel. We will focus on the SLE given by Ax = b where
A ∈ Cm×n is the matrix of linear equations, where x ∈ Cn is the vector of unknown values, and

8

where b ∈ Cm is the vector of known coefficients. We say that the SLE Ax = b is consistent
if it has a solution. We always assume that A has full row-rank in the UD case (m < n) and
full column-rank in the OD case (m ≥ n). As we will see, these conditions are satisfied (with
probability 1) in the wireless communication applications considered here.
We will consider the randomized version of KA explained in Algorithm 1, where at each
iteration the algorithm selects a row of A randomly and independently from the previous
iterations and according to a given probability distribution p = (p1 , . . . , pm )T , where pi ∈ [0, 1]
P
is the probability of selecting row i ∈ [m] and where i∈[m] pi = 1. We denote the random

rank-1 projection operator produced by the random row selection by PR =

aR aH
R
kaR k2

∈ Cn×n , where

R ∈ [m] is a random variable denoting the index of the selected row, where P[R = i] = pi and

where we denote by ai = (Ai,. )H ∈ Cm the conjugate of the i-th row of A. We define the
average of the projection operator PR as
P = ER [PR ] =

m
X
i=1

pi

ai aHi
.
kai k2

Note that P is an n × n positive semi-definite (PSD) matrix, with tr(P) =
we will see in the next sections:

(7)
P

i∈[m]

pi = 1. As

(1) We will design KA such that, depending on the specific scenario, all the estimates generated
across all the iterations of KA belong to the subspace X ⊂ Cn generated by conjugate of the

rows of A, i.e., the column span of AH .

(2) The crucial parameter that will control the convergence speed of KA will be the average
gain of the matrix A over the subspace X defined as follows.

Definition 1 (average gain): Let A be the m × n matrix. Let p ∈ Rm
+ be a probability

distribution over the rows of A and define the average projection operator P produced by
the rows of A as in (7). Let X be the subspace produced by the span of columns of AH . The
average gain of the matrix A over the subspace X is defined as:
xH Px
.
x∈X ,x6=0 kxk2

κX (A, p) = min

(8)

Note that κX (A, p) depends both on the matrix A and the probability distribution p. Moreover,
as tr(P) = 1 and as the subspace X has dimension min{m, n} and P acts on this subspace,
1
we have that κX (A, p) ∈ [0, min{m,n}
]. We will show that the closer κX (A, p) is to

1
min{m,n}

the

9
M = 256 antennas, K = 10 users
eigen values
0.1
minimum eigen value
1
upper bound K
0.08
suboptimal κX (A)

0.15

Eigen Values of P

M = 256 antennas, K = 15 users

0.1

eigen values
minimum eigen value
1
upper bound K
suboptimal κX (A)

0.06

·10−2

8

M = 256 antennas, K = 25 users
eigen values
minimum eigen value
1
upper bound K
suboptimal κX (A)

6

4

0.04

0.05

2
0.02

0

0

2

4

6

8

10

0

0

2

4

6

8

10

12

14

0

0

5

10

15

20

25

Fig. 1: Comparison of the value κX (A, p) for the optimal probability distribution p∗ (obtained by solving the
SDP) with that of the suboptimal one for a matrix A with i.i.d. CN (0, 1) components and with M = 256 rows
(BS antennas) and K ∈ {10, 15, 25} columns (users). The plots illustrate the eigen-value distribution of the average
operator P for the optimal distribution p∗ , where the minimum eigen value corresponds to κX (A, p∗ ). The plots

also illustrate the upper bound
p with pi =

2

kai k
kAk2F

1
K

=

1
min{M,K}

and the suboptimal κX (A) achieved with the suboptimal distribution

. It is seen that for a matrix A with i.i.d. CN (0, 1) components, the suboptimal distribution has

also a good κX (A), which is quite close to the best κX (A, p∗ ).

faster KA converges across the iterations. Note that from (7) for a given matrix A, κX (A, p) =
λmin (P|X ) is given by the minimum eigen-value of P over the subspace X . As P is a linear
function of p and as λmin(P|X ) is a concave function of P [35], κX (A, p) is a concave
function of p over the convex set of probability distributions over the row set of A. Therefore,
finding the best distribution p maximizing κX (A, p) can be posed as maximizing a concave
function over a convex set. This can be formulated as a semi-definite programming (SDP) and
solved with an affordable complexity via convex optimization techniques. However, this typically
incurs the same order of complexity as directly computing the detection/precoding vectors and
is not suitable for the massive MIMO applications which are addressed in this paper. This might
be useful for situations where one solves many equations with the same matrix A, whereby
the complexity of computing p is amortized over time. This is not the case in massive MIMO
applications since the matrix A depends on the estimated channel matrix Q, which changes
from one coherence time to another. Instead, we will use the suboptimal probability distribution
p proposed by [26], where pi =

kai k2
kAk2F

and the probability of selecting each row i ∈ [m] scales

proportionally to its squared l2 -norm kai k2 . Note that computing this probability distribution
incurs a complexity of O(mn), which is definitely much lower than directly solving an SDP to
find the optimal p. Replacing this probability distribution in (7), results in the following average

10

operator:
P = ER [PR ] =

m
X
i=1

=

ai aHi
pi
kai k2

(9)

m
m
X
X
AH A
ai aHi
kai k2 ai aHi
=
=
.
2
2
2
2
kAk
ka
k
kAk
kAk
i
F
F
F
i=1
i=1

(10)

From Definition 1, this yields the following gain parameter for the matrix A, which will be used
in the following.
Definition 2: Let A be the m × n matrix as defined before and let X be a linear subspace of

Cn produced by the column span of AH . We define the normalized minimum gain of the matrix
A along the subspace X as:

kAxk2
,
2
x∈X ,x6=0 kAk2
F kxk

κX (A) = min

(11)
♦

1
Note that κX (A) ∈ [0, min{m,n}
] as before, and the closer κX (A) becomes to

1
min{m,n}

the

faster KA converges across the iterations. Fig. 1 compares κX (A) with the best average gain
κX ,p (A) obtained via optimizing over the probability distribution p. We consider a massive
MIMO scenario where the matrix A has i.i.d. CN (0, 1) components with M = 256 rows
corresponding to BS antennas but different number of columns K ∈ {10, 15, 25} corresponding
to the number of served users. KA for such a matrix A arises in an UL scenario where the users
have spatially white channel vectors. It is seen from Fig. 1 that κX (A) gives an approximation
of the best gain up to a multiplicative factor. This, as we will see, results in a KA that converges
slower than the optimally-tuned KA up to the same multiplicative factor.
For the rest of the paper, we will always use the suboptimal row distribution, which can
be computed in O(mn) iterations. We now prove the following result that gives a unified
performance guarantee for all the cases that will be considered in this paper.
Theorem 1: Let A be an m × n matrix and let b ∈ Cm . Let X be the subspace of Cn

generated by the column span of AH and suppose that Ax = b has a solution x∗ ∈ X . Let

x0 ∈ X be an arbitrary initialization in X and let xt be the estimate at time t of the randomized

KA with the suboptimal distribution as explained before starting from x0 . Then,


E kxt − x∗ k2 ≤ (1 − κX (A))t kx0 − x∗ k2 ,

(12)

11

where the expectation is taken over the randomized row selection in KA, and where κX (A)
denotes the minimum gain of the linear operator A over the subspace X as in Definition 2. 
Proof: From the update equation of KA, we have:
t+1

x

har(t) , xt i − br(t)
−x =x −
ar(t) − x∗
2
kar(t) k
∗

t

= xt − x∗ −
= In −

(13)

har(t) , xt − x∗ i
ar(t)
kar(t) k2

(14)

(xt − x∗ )

(15)

ar(t) aHr(t) 
kar(t)

k2

= (In − Pr(t) )(xt − x∗ ),

(16)

where r(t) ∈ [m] denotes the index of the random row selected by KA at iteration t, where br(t)

is the component of b at position r(t), where we used the fact that x∗ is a solution of Ax = b,
thus, br(t) = har(t) , x∗ i, and where Pr(t) :=

ar(t) aH
r(t)
kar(t) k2

is the random rank-1 projection operator

onto the subspace spanned by ar(t) as defined before. It is not difficult to check that In − Pr(t)

is also a projection operator, thus, (In − Pr(t) )H (In − Pr(t) ) = (In − Pr(t) ). This yields:
kxt+1 − x∗ k2 = (xt − x∗ )H (In − Pr(t) )(xt − x∗ ).

(17)

Note that xt − x∗ is independent of In − Pr(t) since it depends on the random selection of
the rows at iterations before t. Using this independence, conditioning on all the row selections
before t, and taking the expected value over r(t) yields:


E kxt+1 − x∗ k2 r(1), . . . , r(t − 1)

(18)

= (xt − x∗ )H E[In − Pr(t) ](xt − x∗ )

(19)

= (xt − x∗ )H (In − P)(xt − x∗ )

(20)

= kxt − x∗ k2 − (xt − x∗ )H P(xt − x∗ )

(21)

(a)

≤ (1 − κX (A))kxt − x∗ k2 ,

(22)

where in (a) we used the fact that the minimum gain of P for the mentioned suboptimal
distribution is given by κX (A) as in Definition 2. Finally taking expected value with respect to
r(1), . . . , r(t − 1) and applying the induction we obtain:



 

E kxt+1 − x∗ k2 ≤ 1 − κX (A) E kxt − x∗ k2

12

≤ 1 − κX (A)

t+1

kx0 − x∗ k2 ,

(23)

where we used the fact that the initial point x0 is selected deterministically and dropped the
expectation.
Remark 1: From Theorem 1, it is seen that the convergence speed of KA depends on the
parameter κX (A), where a larger κX (A) guarantees a faster convergence. Moreover, since
κX (A) ≤

1
,
min{m,n}

at least an order of O(min{m, n}) iterations are needed for the convergence.

♦
Remark 2: Since the random variable kxt −x∗ k2 is non-negative, from the Markov’s inequality,
its expected value also provides a probabilistic bound on its tail behavior, i.e., for any ζ > 0,
we have:
E[kxt − x∗ k2 ]
ζ
t kx0 − x∗ k2
,
≤ 1 − κX (A)
ζ

P[kxt − x∗ k2 ≥ ζ] ≤

where κX (A) is as defined before.

(24)
♦

In the following, we will design several variants of KA depending on the specific scenario and
apply Theorem 1 to analyze the performance of the resulting algorithm. We will also explain
how the convergence speed depends on the systems parameters such as the number of antennas
M, the number of served users K, and the spatial correlation of the channel vectors of the users.
IV. K ACZMARZ A LGORITHM

FOR DIFFERENT

B EAMFORMING S CENARIOS

In this section, we consider a multiuser massive MIMO system with a BS having M antennas
and K single-antenna users. We remind from Section II that the perfect channel state of the
users is dnoted by an M × K matrix H and the noisy and possibly quantized estimate of H
available at the BS is denoted by the M × K matrix Q.
A. MMSE/ZF detection in the UL
We first consider a UL scenario. From what we have previously discussed in Section II-A,
the ZFD in (2) and MMSED in (3) can be derived, by setting ξ = 0 and ξ =

1
SNR

respectively,

from the following general expression:
b
s = (QH Q + ξIK )−1 QH y.

(25)

13

Let us now focus on the general case where ξ 6= 0.
We make the following observation, that we prove here for the sake of completeness.
Proposition 1: Let Q, y and b
s be as in (25). Then, the signal estimate bs is the optimal

solution of the optimization problem arg minw∈CK kQw − yk2 + ξkwk2.



Proof: Taking the derivative of the cost function gives 2QH (Qw − y) + 2ξw. Setting this

equal to zero yields the optimal solution w∗ = (QH Q + ξIK )−1 QH y which coincides with bs.

Note that we can write the l2 regularized least-squares cost function kQw − yk2 + ξkwk2 in
√
Proposition 1 as kBw − y0 k2 where B = [Q; ξIK ] is an (M + K) × K matrix and where
y0 = [y; 0] is the (M + K)-dim vector obtained by appending a K-dim zero vector to the noisy
received signal y. The signal estimate bs in MMSED can also be written as the least-squares

solution of the newly defined cost function and is given by

bs = (BH B)−1 BH y0 = (QH Q + ξIK )−1 QH y.

(26)

Let us consider the SLE Bw = y0 . Note that this SLE is OD and is always inconsistent when
ξ > 0. Furthermore, it is generally inconsistent even when ξ = 0 due to the presence of noise
in the received measurements y. Thus, a direct application of the randomized KA will result in
a residual error as we discussed in Section II-B. To avoid this, we solve this SLE in two steps.
We first define a new vector
(i)

b0 = Bbs = B(BH B)−1 BH y0 = B(BH B)−1 QH y,
y

(27)

b0 = (BH B)(BH B)−1 QH y = QH y,
BH y

(28)

b0 lies in the subspace X of CM +K spanned the column
where in (i) we used (26). Note that y
of B. Moreover, from

b0 satisfies the UD SLE BH z = b for the unknown z ∈ CM +K and the set of
it is seen that y

b0 . We
known coefficients b = QH y. We first apply the randomized KA to this SLE to recover y

then have the following result.

Proposition 2: Let zt be the estimate obtained at iteration t of KA applied to the UD SLE

BH z = b where b = QH y. Assume that KA starts from the zero initialization z0 = 0 and let
b0 be as in (27). Then
y



b0 k2 ≤ (1 − κX (BH ))t kb
E kzt − y
y0 k2 ,

(29)

14

where κX (BH ) is the minimum gain of the matrix BH over the subspace X generated by the


columns of B.

Proof: From the update equation for KA and applying a simple induction, it is seen that
starting from z0 = 0, all the estimates zt produced by KA lie in the subspace X produced by
the columns of B. Thus, the result simply follows by applying Theorem 1.
b0 is recovered, we solve the OD SLE Bbs = y
b0 in (27) to find the estimate bs. In contrast
Once y

with the initial SLE Bbs = y, which was inconsistent, the OD SLE in (27) is always consistent
b0 always lies in the subspace spanned by the columns of B. Note that since Bbs = y
b0
since y
√
√
is consistent and B = [Q; ξIK ] has ξIK as its submatrix, b
s is simply given by the last K
√
b0 divided by ξ. As a result, there is no need to solve the second SLE Bbs = y
b0
components of y

per se. Moreover, denoting by zt the estimate of KA for the first SLE BH z = b with b = QH y,
we can obtain an estimate b
st at each iteration t of KA by dividing the last K components of
√
√
zt by ξ. By introducing ut ∈ CM and vt ∈ CK such that zt = [ut ; ξvt ], we can summarize

b0 = QH y as in Algorithm 2. A direct inspection shows that, due to the submatrix
KA for BH y
√
√
ξIK in B = [Q; ξIK ], at each iteration of KA only a single randomly selected component

of the estimate vt is updated, where vt converges to the least-squares estimate of the symbols bs

corresponding to the MMSED. In particular, one can obtain the ZFD estimate b
s by setting ξ = 0

in Algorithm 2 .

To analyze the convergence speed of KA, from Proposition 2, we need to compute the gain

of the matrix B along the subspace X spanned by the rows of B. Since every vector in X can

be written as Bw for some w ∈ CK , a direct calculation shows that κX (BH ) is given by:
kBH xk2
x∈X ,x6=0 kBH k2 kxk2
F

κX (BH ) = min
=
=
Note that

λmin (QH Q)
kQk2F

≤

1
K

kBH Bwk2
2
w∈CK ,w6=0 kBH kF kBwk2
min

λmin (QH Q) + ξ
λmin (BH B)
.
=
kBk2F
kQk2F + Kξ

(30)

and it is not difficult to check that for any ξ > 0, we have
1
λmin(QH Q) + ξ
λmin (QH Q)
≤ .
≤
2
2
kQkF
kQkF + Kξ
K

(31)

This implies that the algorithm converges faster for larger values of ξ > 0 (lower UL SNR)
and the slowest convergence corresponds to the ZFD with ξ = 0. When the channel vectors of

15

Algorithm 2 Kaczmarz Algorithm for the UL
1:

Input: estimate of channel state Q ∈ CM ×K , received noisy UL signal y ∈ RM , and
MMSED parameter ξ > 0.

2:

Compute b = QH y ∈ CK .

3:

Define ut ∈ CM and vt ∈ CK with u0 = 0, v0 = 0.

4:

for k = 0, 1, . . . , T − 1 do

5:

Pick a row r(t) of QH denoted by qHr(t) with a probability proportional to kqr(t) k2 + ξ

given by

kqr(t) k2 +ξ
.
kQk2F +Kξ

t
br(t) −hqr(t) ,ut i−ξ vr(t)

6:

Compute the residual γ t :=

7:

Update ut+1 = ut + γ t qr(t) .

8:

t+1
t
Update vr(t)
= vr(t)
+ γ t , and vjt+1 = vjt for j 6= r(t).

9:

kqr(t) k2 +ξ

.

end for

10:

Output: Set w = vT .

11:

Output: Set the signal estimate to b
s = w.

the users have white complex Gaussian components (no spatial correlation), from well-known
results in random matrix theory, it results that

√
√
√
(1 − ρ)2
( M − K)2
λmin(QH Q)
=
,
=
κX (B ) ≥
kQk2F
MK
K
H

where we define ρ :=

K
M

(32)

as the loading factor of the massive MIMO system. Therefore, from

Proposition 2 it results that O( (1−K√ρ)2 ) iterations are sufficient to reasonably estimate the symbol
sequence b
s for the MMSED receiver for all ranges of ξ > 0. In particular, it is seen that the
required estimation time grows linearly with the number of users K. Since each step of KA

requires multiplying two M-dim vectors, thus, requiring M operations, the total computational
KM
√ 2 ). To this we should also add O(KM) operations needed for
complexity scales as O( (1−
ρ)

computing the sub-optimal row-sampling probability distribution and O(KM) operation needed
for calculating QH y in Algorithm 2.
B. RZF/ZF Precoding in the DL
In this section, we consider RZF precoding in the DL. Suppose s = (s1 , . . . , sK )T is the
symbol sequence BS intends to send to the users in the DL. In the RZFBF in the DL, the BS
transmits the M-dim vector y := Q(QH Q + ξIK )−1 s via its M antennas. We observe that y lies

16

in the subspace spanned by the columns of Q. We design a new KA that finds the K-dim signal
w = (QH Q + ξIK )−1 s, from which we can obtain y = Qw via a matrix-vector multiplication.
We prove the following result.
Proposition 3: Let Q, s and w = (QH Q + ξIK )−1 s be as introduced before. Then, w is the
optimal solution of the following optimization problem w = arg minx∈CK kQxk2 + ξkx − sξ k2 ,
where sξ = ξs .



Proof: Taking the derivative with respect to x and setting it equal to 0 yields QH Qx +
ξx − s = 0. This gives the optimal solution as x∗ = (QH Q + ξIK )−1 s, which coincides with
w = (QH Q + ξIK )−1 s.

Proposition 3 implies that we can obtain w = (QH Q + ξIK )−1 s as the optimal solution of
kQxk2 + ξkx − sξ k2 . We first write this as the least-squares minimization kBx − bk2 where
√
√
B = [Q; ξIK ] is an (M + K) × K matrix and where b = [0; ξsξ ] is an (M + K)-dim vector
√
with ξsξ = √sξ in its last positions. The SLE Bw = b is OD and should be solved for the
vector w ∈ CK , from which we obtain the desired signal y = Qw to be transmitted in the DL.
This SLE is inconsistent unless s = 0. Thus, we follow the same procedure as in the design of
MMSED for the UL. We first define z = Bw and solve the UD SLE
p s
(33)
BH z = BH b = ξ √ = s,
ξ
for the z and then solve the OD but consistent SLE Bw = b
z, where b
z is the estimate of z
√
obtained from KA applied to (33). Note that, as in the MMSED, due to the presence of ξIK
√
submatrix in B = [Q; ξIK ], the solution of the latter SLE is simply obtained by dividing the
√
√
last K components of b
z by ξ. Denoting zt = [ut ; ξvt ] the estimate at iteration t of KA, we
can summarize this in Algorithm 3.

In terms of speed of convergence, a direct calculation shows that the minimum gain of the

matrix BH over the subspace X spanned by the columns of B is given by:
κX (B) =

λmin (QH Q) + ξ
λmin (BH B)
.
=
kBk2F
kQk2F + Kξ

(34)

In particular, when the channel vectors of the users have white complex Gaussian components
(no spatial correlation), a direct calculation as in the case of MMSED shows that the total
KM
√ 2 ). To this we should also add O(KM)
computational complexity of KA scales as O( (1−
ρ)

operations needed for computing the sampling probability distribution and O(KM) operation
needed for calculating Qw in Algorithm 3.

17

Algorithm 3 Kaczmarz Algorithm for the DL
1:

Input: estimate of channel state Q ∈ CM ×K , user symbols s ∈ RK , and RZFBF parameter
ξ > 0.

2:

Define ut ∈ CM and vt ∈ CK with u0 = 0, v0 = 0.

3:

for k = 0, 1, . . . , T − 1 do

Pick a row r(t) of QH denoted by qHr(t) with a probability proportional to kqr(t) k2 + ξ

4:

given by

kqr(t) k2 +ξ
.
kQk2F +Kξ

t
sr(t) −hqr(t) ,ut i−ξ vr(t)

5:

Compute the residual γ t :=

6:

Update ut+1 = ut + γ t qr(t) .

7:

t+1
t
Update vr(t)
= vr(t)
+ γ t , and vjt+1 = vjt for j 6= r(t).

8:

end for

9:

Output: Set w = vT .

10:

kqr(t) k2 +ξ

.

Output: Set the desired DL signal to y = Qw.

A direct inspection in Algorithm 2 and 3 also reveals that there is a symmetry between the
MMSED for the UL and RZFBF for the DL: in the former the noisy received signal y is
multiplied by QH and the resulting vector b = QH y enters KA, which yields the estimate of
the user symbols b
s, whereas in the latter the user symbols s to be transmitted in the DL enters

KA and the output of KA is multiplied by Q to yield the DL signal. Moreover, the speed of
convergence is the same for both cases. This highly resembles the well-known UL-DL duality

in the underlying vector Gaussian broadcast/multiaccess channel [36].
V. C OMPUTATION

OF THE

D ETECTOR /P RECODER M ATRIX

Up to now, we always used our proposed KA to compute an estimate of the user symbols from the
received signal in the UL or to compute suitable precoding vectors in the DL. More specifically,
we did not compute any UL/DL detection/precoding matrix from the available channel state
directly. In this section, we show that our proposed KA can be generalized to compute also
the matrices. This results in a saving in computational complexity when the channel coherence
block is quite large such that the detection/precoding matrices can be computed once and used for
all time-frequency signal dimensions in the same coherence block. We will need the following
result.

18

Proposition 4: Consider an arbitrary SLE Ax = b where A is the m×n matrix of equations,
where x ∈ Cn is the set of unknowns, and where b ∈ Cm is the set of known coefficients.

Let xt be the estimate of KA at iteration t starting from the zero initialization x0 = 0. Then,

xt = G t (A)b where G t (A) is a linear operator that depends on A and the internal randomization


of KA but not on b.

Proof: We use induction on t. For t = 0, we have x0 = 0 which is trivially a linear function
(zero function) of b. Let us assume that the induction hypothesis is true for the t-th iteration,
i.e., xt = G t (A)b, and prove the claim for t + 1. Suppose r(t) ∈ [m] is the index of the row
selected at iteration t and let us denote the conjugate of this row by ar(t) . Then, we have
t+1

x

br(t) − har(t) , xt i
=x +
ar(t)
kar(t) k2
t

(i)

= (I −

(ii)

=

ar(t) aHr(t)
kar(t) k2

(I −

t

)x +

ar(t) aHr(t)
kar(t)

=: G t+1 (A)b,

k2

(35)

ar(t) eHr(t)
kar(t) k2

)G t (A) +

b

ar(t) eHr(t)
kar(t) k2

(36)
!

b

(37)
(38)

where in (i) we defined er(t) as the canonical vector that has 1 at the index r(t) and is zero
elsewhere and replaced br(t) = eHr(t) b, and where in (ii) we used the induction hypothesis that
xt = G t (A)b is a linear function of b where G t (A) does not depend on b. This implies that
xt+1 = G t+1 (A)b, where it is seen that the linear operator G t+1 (A) depends on A and on
the index of the random row r(t) selected by KA but does not depend on b. This proves the
induction step and completes the proof.
We use Proposition 4 to design detection/precoding matrices for the UL and the DL as follows.
In the DL, we use Algorithm 3 and run K KAs in parallel where the input to the i-th KA is
s = ei where ei denotes the i element of the canonical basis, which has 1 in its i-th component
and is 0 elsewhere. Let wi ∈ CK denote the output of the i-th KA after suitable convergence
and let W = [w1 , . . . , wK ] ∈ CK×K . Then, from Algorithm 3 and the underlying linearity of

KA proved in Proposition 4, we immediately obtain an approximation of the DL beamfoming
matrix as G = QW, where Q is the noisy estimate of the channel state available at the BS.
Note that multiplication by Q is due to post-multiplication with Q in Algorithm 3. It is seen that
the output of the matrix G for an input ei is given by Gei = Qwi , which as expected coincides

19

with the DL precoding vector computed by the i-th KA. Note that if all the K KAs can be run
in parallel, computing the precoding matrix G will require O(MK) iterations. Also, in terms
of computational complexity, it is beneficial to keep the precoding matrix G in the factorized
form G = QW, and compute the precoding vector as follows
y = Gs = QWs = (Q(Ws)),

(39)

where the parentheses in the last expression show the order of multiplication. Otherwise, we
will need an additional O(MK 2 ) operations to compute G directly, which incurs an order of
O(K) extra complexity.
A similar procedure yields the detection matrix for the UL. For the UL, we apply Algorithm
2 and consider K parallel KAs, where the input to the i-KA is b = ei ∈ CK (see the vector
b in Algorithm 2), where ei denotes the i-th canonical basis. Let wi be the output of the i-th
KA after convergence and let W = [w1 , . . . , wK ] ∈ CK×K . Then, from the underlying linearity
of KA proved in Proposition 4, it results that the desired detection matrix for the UL is given
by G = WQH (see the pre-multiplication by QH in Algorithm 2). Analogous to the DL, it is
beneficial in terms of complexity to keep G in the factorized from and compute the estimate of
user symbols in the UL by
b
s = Gy = WQH y = (W(QH y)),

(40)

where again the parentheses in the last expression show the order of multiplication. As in the
DL, assuming that K KAs can be run in parallel, computing the detection matrix G (of course
in its factorized form G = WQH ) requires O(MK) iterations.
VI. ACHIEVABLE E RGODIC R ATE
In this section, we derive lower and upper bounds on the achievable ergodic rate of our proposed
KA algorithm in section IV. For simplicity, we focus on a UL scenario, with a similar argument
holding for a DL scenario. In the UL, the received signal at the BS is given by y = Hs+n, where
s ∈ CK is the symbol sequence transmitted from the users, where H is the M × K matrix of the
true channel state, and where n is the AWGN noise. From (27) and (28) and our explanation in
Section IV-A and also Proposition 4, it can easily be seen that the signal estimate bs produced by
KA will be a linear function of y assuming zero initilization. More precisely, we have bs = G (Q)y

20

where G (Q) depends on the noisy channel state and the internal randomization of KA but is
independent of the received signal y. Denoting Gup = G (Q) = [g1 , . . . , gk ]H ∈ CK×M with gk
being the beamforming vector of the k-th user, we can write the time-varying (varying across
coherence times due to channel randomness and varying across slots inside a coherence block
due to internal randomness of KA) interference matrix as Tup = Gup H.
We use the upper and lower bounds on the achievable ergodic rate developed in [37]. In
particular, by treating the interference as noise and by coding across several coherence blocks,
an upper bound on the ergodic capacity is given by2
!#
"
2
|Tup
k,k |
,
Rk = E log 1 + P
up 2
2
k ′ 6=k |Tk,k ′ | + σk
where σk2 =

1
gH g
SNR k k

(41)

with gk denoting the k-th row of Gup , and where the expectation is taken

over all the randomness of the channel state H and the beamforming vector Gup including the
P
randomness due to KA. This yields the upper bound R = k∈[K] Rk on the achievable sum
rate. We will also use the following lower bound on the ergodic rate from [37]:
!
2
]|
|E[Tup
P k,k up 2
Rk = log 1 +
,
2
)
+
Var(Tup
k,k
k ′ 6=k E[|Tk,k ′ | ] + E[σk ]

(42)

where Var(.) denotes the variance of a random variable. This yields the lower bound R =
P
k∈[K] Rk on the ergodic rate. Following a similar argument, we obtain the corresponding

upper/lower bounds on the ergodic capacity for the DL. We use Monte Carlo simulations in the
next section to empirically evaluate these lower and upper bounds on the performance of our
algorithm.
VII. S IMULATION

RESULTS

In this section, we evaluate the performance of our proposed technique empirically via numerical
simulations.
A. Channel Model
For all the simulations, we assume a correlated channel model for every user k ∈ [K] where

the true channel vector of the user k is given by hk = Φ1/2 zk where zk ∼ CN (0, IM ) and
2

This upper bound is obtained in [38] under the condition that the true channel sate H is known at the BS. This is an upper

bound for our scenario, where only a noisy version of H is assumed to be known at the BS.

21

where Φ is the covariance matrix of hk with entries given by [Φ]i,j = a|i−j| for some correlation
parameter a ∈ [0, 1]. We assume that the entries of zk , thus, those of hk , vary independently
from one coherence period to another. We also assume that the estimated channel matrix Q,
available at the BS in each coherence period, is a noisy version of the true channel matrix H
√
given by Q = 1 − τ 2 H + τ N where N is a matrix with i.i.d. entries [N]m,k ∼ CN (0, 1),
and where τ ∈ [0, 1] is the estimation quality parameter with τ = 0 yielding a perfect channel
estimation.
B. Downlink
In the DL, we compare the performance of our method with that of the Truncated Polynomial
Expansion (TPE) technique proposed in [18] in terms of spectral efficiency and total number
of operations as a criterion for affordable computational complexity. Fig. 2 illustrates the upper
and the lower bound on the achievable per-user spectral efficiency versus SNR of our proposed
algorithm. It is seen from Fig. 2 that for both ZFBF and RZFBF, 32MK operations guarantee
a perfect match with the original ZFBF and RZFBF up to 20 dB SNR respectively, where for
lower SNRs, a lower number of iterations is sufficient. It is seen that our proposed algorithm has
a slightly more complexity than the TPE proposed in [18]. However, it has also the following
crucial advantages compared with [18]:
1) In TPE, matrix inversion is approximated via Taylor series expansion with a finite number
of terms J, known as the TPE order. So for a certain value of J, one needs to compute J
coefficients of the expansion in advance. For higher signal to noise ratios, typically higher
value of J is needed. This can be seen in Fig. 2 where for lower SNR values, lower J
can also make a good estimate. Moreover, when J changes, one needs to compute all
the coefficients from scratch, which incurs additional complexity. This is in contrast to
our proposed method where only the number of iterations is needed to be increased or
decreased for higher or lower values of SNR respectively. This can also be seen in Fig. 2.
2) By changing the number of scheduled users, all the computations of the polynomial
coefficients in TPE must be done from scratch. In our proposed method, however, everything
is calculated each time a signal is going to be sent in the DL, or received in the UL, and
hence no extra calculations are required.
3) The proposed algorithm does not require the knowledge of channel estimation quality

Average per-User Spectral Efficiency [bit/s/Hz]

22

Downlink, τ = 0.1, M = 256, K = 32
10MK operations
original ZFBF-ub
original RZFBF-ub
original ZFBF-lb
original RZFBF-lb
P ZFBF-ub-10
P RZFBF-ub-10
TPE RZFBF-ub-10
P ZFBF-lb-10
P RZFBF-lb-10
TPE RZFBF-lb-10

8
6
4
2
0
−20

−10

0
SNR [dB]

18MK operations
8
6
4
2

10

0
20 −20

original ZFBF-ub
original RZFBF-ub
original ZFBF-lb
original RZFBF-lb
P ZFBF-ub-18
P RZFBF-ub-18
TPE RZFBF-ub-18
P ZFBF-lb-18
P RZFBF-lb-18
TPE RZFBF-lb-18

−10

0
SNR [dB]

32MK operations
8
6
4

original ZFBF-ub
original RZFBF-ub
original ZFBF-lb
original RZFBF-lb
P ZFBF-ub-32
P RZFBF-ub-32
P ZFBF-lb-32
P RZFBF-lb-32

2

10

0
20 −20

−10

0
SNR [dB]

10

20

Fig. 2: Upper and lower bound of the per user ergodic spectral efficiency in the DL versus SNR for a loading
factor of ρ =

K
M

= 18 , channel estimation quality factor τ = 0.1 and uncorrelated channels (a = 0), is compared for

the proposed method against original linear precoders (ZFBF and RZFBF) and TPE scheme [18].The letters ’P’, ’lb’
and ’ub’ used in the figure denote proposed method, lower bound and upper bound respectively. The numbers at the
end show computational complexity as a factor of M K which finally shows the number of arithmetic operations.
The computational complexity budget is assumed to be 10M K according to J = 3 in TPE, 18M K according to
J = 5 in TPE and 32M K.

parameter τ .
It is also worthwhile to mention that, the number of operations in both the proposed method
and the TPE grows with the coherence period. For very large coherence periods traditional
matrix-matrix multiplication will outperform both TPE and the proposed method, in terms of
number of operations, although this happens only when the channel varies really slowly.
C. Uplink
In the UL, we compare the performance of our algorithm with a very recent approximate
message passing (AMP) technique proposed in [23, 39], especially with [39] where the authors
used AMP to build a nonparametric MMSE detector for massive MIMO that is able to obtain an
estimate of the signal and noise power. Fig. 3 illustrates our simulation results for the UL scenario
in terms of upper and lower bounds on the achievable per-user spectral efficiency. It is seen that
in uncorrelated channels (a = 0), for both ZFD and MMSED, 40MK operations are sufficient
for our proposed method to achieve a perfect match to the upper and lower bounds of the original

23

Uplink, τ = 0, M = 256, K = 32
a=0, 24MK operations
10

Average per-User Spectral Efficiency [bit/s/Hz]

8
6
4
2
0
−20

a=0, 40MK operations
10

original ZFD-ub
original MMSED-ub
original ZFD-lb
original MMSED-lb
P ZFD-ub-24
P MMSED-ub-24
S-AMP-MMSED-ub-24
P ZFD-lb-24
P MMSED-lb-24
S-AMP MMSED-lb-24

−10

0
SNR [dB]

8
6
4
2
10

20

0
−20

a=0.6, 24MK operations
10
8
6
4
2
0
−20

0
SNR [dB]

10
8
6
4
2
10

−10

0
SNR [dB]

10

20

a=0.6, 40MK operations

original ZFD-ub
original MMSED-ub
original ZFD-lb
original MMSED-lb
P ZFD-ub-24
P MMSED-ub-24
S-AMP-MMSED-ub-24
P ZFD-lb-24
P MMSED-lb-24
S-AMP-MMSED-lb-24

−10

original ZFD-ub
original MMSED-ub
original ZFD-lb
original MMSED-lb
P ZFD-ub-40
P MMSED-ub-40
P ZFD-lb-40
P MMSED-lb-40

20

0
−20

original ZFD-ub
original MMSED-ub
original ZFD-lb
original MMSED-lb
P ZFD-ub-40
P MMSED-ub-40
S-AMP-MMSED-ub-40
P ZFD-lb-40
P MMSED-lb-40
N-AMP MMSED-lb-40

−10

0
SNR [dB]

10

20

Fig. 3: Upper and lower bound of the per user ergodic spectral efficiency in the UL versus SNR for a loading factor
of ρ =

K
M

= 18 , channel estimation quality factor τ = 0 and correlation parameter a = {0, 0.6} is compared for the

proposed method against original linear detectors (ZFD and MMSED), standard AMP [39] and nonparametric AMP
[40]. The letters ’P’, ’S’, ’N’, ’lb’ and ’ub’ used in the figure denote proposed method, standard, nonparametric,
lower bound and upper bound respectively. The numbers at the end show computational complexity as a factor of
M K which finally shows the number of arithmetic operations. The computational complexity budget is assumed
to be {24M K, 40M K}.

ZFD and MMSED up to 20 dB SNR respectively.3 It is also seen that the AMP MMSED [39]
is able to converge to the original MMSED when the channel is uncorrelated (a ≈ 0), with
3

for the upper bounds, 32M K iterations suffice.

24

an order of 24MK operations.4 However, for a slightly correlated channel (a ≥ 0.6), neither
the standard AMP [39] nor the nonparametric AMP [23] are able to suitably recover the user
symbols even for a large number of iterations whereas our proposed algorithm becomes slower
due to channel correlation but still converges by slightly increasing the number of iterations.
D. Gap-to-Capacity versus Computational Complexity
In this part, we compare our proposed method with the TPE [18] and the AMP [23, 39] in terms
of normalized gap-to-capacity versus the number of operations t given by

Smax −St
,
Smax

where St

denotes the upper/lower bound of per-user spectral efficiency of a specific beamforming/detection
scheme after t operations and where Smax denotes the maximum possible upper/lower bound of
per-user spectral efficiency achieved by original linear schemes. This is illustrated in Fig. 4 for
different beamforming schemes in the DL and detection schemes in the UL for different channel
correlations a = 0 and a = 0.6. It is seen that for uncorrelated channels, i.e a = 0, the proposed
method can achieve Gap = 10−2 with nearly 40MK operations. It is also worthwhile to mention
that the AMP [39, 40] are quite efficient for uncorrelated channels but fail in correlated channels.
VIII. C ONCLUSION
In this paper, we proposed a novel statistically robust approach for reducing the computational
complexity of beamforming in massive MIMO systems. We showed that designing all variants of
linear precoders/ detectors for the DL and the UL can be posed as the solution of an appropriate
set of linear equations and applied a randomized variant of Kaczmarz algorithm to solve the
resulting equations. We also proposed a unified framework to analyze the performance of our
proposed algorithm. Our algorithm does not require knowledge of the statistics of users’ channel
vectors and its convergence speed depends can be fully specified in terms of the condition
number of the channel matrix, which depends on the number of BS antennas M, the number
of served users K, as well as the statistical correlation of the channel vectors of the users.
Our numerical simulations show that our algorithm is more robust than other schemes proposed
in the literature based on random matrix theory and approximate message passing and has a
4

This translates into 6 iterations in AMP, since AMP carries out two matrix-vector operations in each iteration which will

results in an order of 4M K additions and multiplications per iteration.

25

Gap, τ = 0, M = 256, K = 32, SNR = 20 dB
Uplink, a = 0
P ZFD-ub
P MMSED-ub
S-AMP-MMSED-ub
P ZFD-lb
P MMSED-lb
S-AMP-MMSED-lb
N-AMP-MMSED-lb

101

Gap to Expected Spectral Efficiency (normalized)

Downlink, a = 0

100

10−1

100

10−1

10−2
2M K 8M K

101

P ZFBF-ub
P RZFBF-ub
TPE RZFBF-ub
P ZFBF-lb
P RZFBF-lb
TPE RZFBF-lb

101

16M K 24M K 32M K 40M K

10−2
2M K 8M K

16M K 24M K 32M K 40M K

Complexity Budget

Complexity Budget

Uplink, a = 0.6

Downlink, a = 0.6

P ZFD-ub
S-AMP-MMSED-ub
P MMSED-lb
N-AMP-MMSED-lb

P MMSED-ub
P ZFD-lb
S-AMP-MMSED-lb

100

100

10−1

10−1

10−2
2M K 8M K

16M K 24M K 32M K 40M K

P ZFBF-ub
P RZFBF-ub
TPE RZFBF-ub
P ZFBF-lb
P RZFBF-lb
TPE RZFBF-lb

101

10−2
2M K 8M K

Complexity Budget

16M K 24M K 32M K 40M K

Complexity Budget

Fig. 4: Gap to the maximum spectral efficiency for the proposed scheme is compared against that of the TPE
[18], standard AMP [39] and nonparametric AMP [40]. Loading factor, SNR, channel estimation quality parameter
and channel correlation parameter are ρ =

K
M

= 81 , 20 dB, τ = 0 and a = {0, .06} respectively. The legends ‘P’,

‘S’, ‘N’ in the figure refer to our proposed method, standard and nonparametric respectively.

competitive performance.
R EFERENCES
[1] M. N. Boroujerdi, S. Haghighatshoar, and G. Caire, “Low-complexity and statistically robust beamformer
design for massive mimo systems,” submitted to Workshop on Smart Antennas (WSA) 2018.
[2] L. Lu, G. Y. Li, A. L. Swindlehurst, A. Ashikhmin, and R. Zhang, “An overview of massive mimo: Benefits
and challenges,” IEEE Journal of Selected Topics in Signal Processing, vol. 8, no. 5, pp. 742–758, Oct 2014.
[3] T. L. Marzetta, “Noncooperative cellular wireless with unlimited numbers of base station antennas,” IEEE
Trans. on Wireless Commun., vol. 9, no. 11, pp. 3590–3600, Nov. 2010.

26

[4] F. Rusek, D. Persson, B. K. Lau, E. G. Larsson, T. L. Marzetta, O. Edfors, and F. Tufvesson, “Scaling up
mimo: Opportunities and challenges with very large arrays,” IEEE Signal Processing Magazine, vol. 30, no. 1,
pp. 40–60, Jan 2013.
[5] E. G. Larsson, O. Edfors, F. Tufvesson, and T. L. Marzetta, “Massive mimo for next generation wireless
systems,” IEEE Communications Magazine, vol. 52, no. 2, pp. 186–195, February 2014.
[6] H. Yang and T. L. Marzetta, “Performance of conjugate and zero-forcing beamforming in large-scale antenna
systems,” IEEE Journal on Selected Areas in Communications, vol. 31, no. 2, pp. 172–179, February 2013.
[7] G. Caire and S. Shamai, “On the achievable throughput of a multiantenna Gaussian broadcast channel,” IEEE
Transactions on Information Theory, vol. 49, no. 7, pp. 1691–1706, 2003.
[8] H. Weingarten, Y. Steinberg, and S. S. Shamai, “The capacity region of the Gaussian multiple-input
multiple-output broadcast channel,” IEEE Transactions on Information Theory, vol. 52, no. 9, pp. 3936–3964,
2006.
[9] K. S. Schneider, “Optimum detection of code division multiplexed signals,” IEEE Transactions on Aerospace
and Electronic Systems, vol. AES-15, no. 1, pp. 181–185, Jan 1979.
[10] R. Lupas and S. Verdu, “Linear multiuser detectors for synchronous code-division multiple-access channels,”
IEEE Transactions on Information Theory, vol. 35, no. 1, pp. 123–136, Jan 1989.
[11] M. K. Varanasi, “Noncoherent detection in asynchronous multiuser channels,” IEEE Transactions on
Information Theory, vol. 39, no. 1, pp. 157–176, Jan 1993.
[12] Z. Xie, R. T. Short, and C. K. Rushforth, “A family of suboptimum detectors for coherent multiuser
communications,” IEEE Journal on Selected Areas in Communications, vol. 8, no. 4, pp. 683–690, May
1990.
[13] S. Yang and L. Hanzo, “Fifty years of mimo detection: The road to large-scale mimos,” IEEE Communications
Surveys Tutorials, vol. 17, no. 4, pp. 1941–1988, Fourthquarter 2015.
[14] S. D. L. Moshavi Shimon, Kanterakis Emmanuel G., “Multistage linear receivers for ds-cdma systems,”
International Journal of Wireless Information Networks, vol. 3, no. 1, pp. 1–17, Jan 1996. [Online].
Available: https://doi.org/10.1007/BF02106658
[15] R. R. Muller and S. Verdú, “Design and analysis of low-complexity interference mitigation on vector channels,”
IEEE Journal on Selected Areas in Communications, vol. 19, no. 8, pp. 1429–1441, 2001.
[16] G. M. A. Sessler and F. K. Jondral, “Low complexity polynomial expansion multiuser detector for cdma
systems,” IEEE Transactions on Vehicular Technology, vol. 54, no. 4, pp. 1379–1391, July 2005.
[17] S. Zarei, W. Gerstacker, R. R. Mller, and R. Schober, “Low-complexity linear precoding for downlink
large-scale mimo systems,” in 2013 IEEE 24th Annual International Symposium on Personal, Indoor, and
Mobile Radio Communications (PIMRC), Sept 2013, pp. 1119–1124.
[18] A. Kammoun, A. Müller, E. Björnson, and M. Debbah, “Linear precoding based on polynomial expansion:
Large-scale multi-cell mimo systems,” IEEE Journal of Selected Topics in Signal Processing, vol. 8, no. 5,

27

pp. 861–875, 2014.
[19] A. Benzin, G. Caire, Y. Shadmi, and A. Tulino, “Truncated Polynomial Expansion Downlink Precoders and
Uplink Detectors for Massive MIMO,” ArXiv e-prints, Nov. 2017.
[20] S. Haghighatshoar and G. Caire, “Massive mimo channel subspace estimation from low-dimensional
projections,” IEEE Transactions on Signal Processing, vol. 65, no. 2, pp. 303–318, 2017.
[21] ——, “Low-complexity massive mimo subspace estimation and tracking from low-dimensional projections,”
arXiv preprint arXiv:1608.02477, 2016.
[22] M. L. Honig and W. Xiao, “Performance of reduced-rank linear interference suppression,” IEEE Transactions
on Information Theory, vol. 47, no. 5, pp. 1928–1946, Jul 2001.
[23] R. Ghods, C. Jeon, G. Mirza, A. Maleki, and C. Studer, “Optimally-tuned nonparametric linear equalization
for massive mu-mimo systems,” arXiv preprint arXiv:1705.02985, 2017.
[24] M. Bayati and A. Montanari, “The dynamics of message passing on dense graphs, with applications to
compressed sensing,” IEEE Transactions on Information Theory, vol. 57, no. 2, pp. 764–785, 2011.
[25] S. Kaczmarz, “Angenäherte auflsung von systemen linearer gleichungen,” Bulletin International de l’Acadmie
Polonaise des Sciences et des Lettres. Classe des Sciences Mathmatiques et Naturelles. Srie A, Sciences
Mathmatiques, pp. 335–357, 1937.
[26] T. Strohmer and R. Vershynin, A Randomized Solver for Linear Systems with Exponential Convergence. Berlin,
Heidelberg: Springer Berlin Heidelberg, 2006, pp. 499–507, https://doi.org/10.1007/11830924 45.
[27] D. Needell, R. Ward, and N. Srebro, “Stochastic gradient descent, weighted sampling, and the randomized
kaczmarz algorithm,” in Advances in Neural Information Processing Systems, 2014, pp. 1017–1025.
[28] Y. Chi and Y. M. Lu, “Kaczmarz method for solving quadratic equations,” IEEE Signal Processing Letters,
vol. 23, no. 9, pp. 1183–1187, 2016.
[29] C. Shepard, H. Yu, N. Anand, E. Li, T. Marzetta, R. Yang, and L. Zhong, “Argos: Practical many-antenna base
stations,” in Proceedings of the 18th annual international conference on Mobile computing and networking.
ACM, 2012, pp. 53–64.
[30] G. Caire, N. Jindal, M. Kobayashi, and N. Ravindran, “Multiuser MIMO achievable rates with downlink
training and channel state feedback,” IEEE Transactions on Information Theory, vol. 56, no. 6, pp. 2845–2866,
2010.
[31] J. Hoydis, S. ten Brink, and M. Debbah, “Massive mimo in the ul/dl of cellular networks: How many antennas
do we need?” IEEE Journal on Selected Areas in Communications, vol. 31, no. 2, pp. 160–171, February
2013.
[32] C. B. Peel, B. M. Hochwald, and A. L. Swindlehurst, “A vector-perturbation technique for near-capacity
multiantenna multiuser communication-part i: channel inversion and regularization,” IEEE Transactions on
Communications, vol. 53, no. 1, pp. 195–202, Jan 2005.
[33] D. Needell, “Randomized kaczmarz solver for noisy linear systems,” BIT Numerical Mathematics, vol. 50,

28

no. 2, pp. 395–403, Jun 2010. [Online]. Available: https://doi.org/10.1007/s10543-010-0265-5
[34] A. Zouzias and N. Freris, “Randomized extended kaczmarz for solving least-squares,” arXiv preprint
arXiv:1205.5770v3, 2013.
[35] T. Tao, Topics in random matrix theory. American Mathematical Society Providence, RI, 2012, vol. 132.
[36] S. Vishwanath, N. Jindal, and A. Goldsmith, “Duality, achievable rates, and sum-rate capacity of gaussian
mimo broadcast channels,” IEEE Transactions on Information Theory, vol. 49, no. 10, pp. 2658–2668, 2003.
[37] H. Y. H. Q. N. Thomas L. Marzetta, Erik G. Larsson, Fundamentals of Massive MIMO. Cambridge: cambridge
university press, 2016.
[38] G. Caire, “On the ergodic rate lower bounds with applications to massive mimo,” arXiv preprint
arXiv:1705.03577, 2017.
[39] C. Jeon, R. Ghods, A. Maleki, and C. Studer, “Optimality of large mimo detection via approximate message
passing,” in 2015 IEEE International Symposium on Information Theory (ISIT), June 2015, pp. 1227–1231.
[40] R. Ghods, C. Jeon, G. Mirza, A. Maleki, and C. Studer, “Optimally-tuned nonparametric linear
equalization for massive MU-MIMO systems,” CoRR, vol. abs/1705.02985, 2017. [Online]. Available:
http://arxiv.org/abs/1705.02985

