---
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 \gg 1$, estimating the channel covariance matrices of the users is also challenging in terms of both 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 our proposed technique is computationally very competitive and yields quite a comparable performance while it does not require the knowledge of the statistics of the users channel vectors.' 
bibliography:
- 'references.bib'
title: 'Low-Complexity Statistically Robust Precoder/Detector Computation for Massive MIMO Systems[^1]'
---

Massive MIMO, Kaczmarz Algorithm, Detection, Precoding. 

Introduction
============

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 a high sum spectral efficiency, while the large antenna diversity achieved by using $M \gg K$ antennas provides [*channel hardening*]{}, such that also the per-user rates are very stable and robust against statistical variations of the propagation environment [@Marzetta-TWC10; @Intro_MM_Ashikhmin; @Intro_MM_Rusek; @Intro_MM_Larsson].  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 \to \infty$, the intra cellular interference vanishes so does the effect of noise and fast (small-scale) fading [@Marzetta-TWC10].  Moreover, due to channel hardening and to the fact that for $M$ much larger than $K$ the users channel vectors tend to become almost 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 [@Marzetta-TWC10; @downlink_MM_linear; @caire2003achievable; @weingarten2006capacity]. 

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 [@uplink_decorrelator; @uplink_decorrelator_sync; @uplink_decorrelator_async]), regularized ZF (RZF), where the regularization factor is tuned to the level of the noise at the users receivers, and *linear minimum mean squared error* (MMSE) detector, which achieves the optimal Signal-to-Interference plus Noise Ratio (SINR) among all linear receivers [@uplink_MMSE] (we refer to [@MIMO_detection_survey] for a more detailed review of the available techniques for MIMO precoding and detection).  Computing the corresponding linear precoding/detection matrices from the estimated users channel vectors requires matrix inversion and matrix-matrix multiplication, which requires $O(K^3)$ and $O(MK^2)$ operations, respectively.  This imposes a very high computational complexity on massive MIMO implementations since the number of BS antennas $M$ and the number of served users $K$ can be very large.  Several works have recently reported hardware chips that are designed to specifically compute beamforming matrices in massive MIMO, see e.g. [@MM-implementation] for a chip design for $M=128$ antennas and $K=8$ users.  Many implementations of massive MIMO, however, are still using general purpose signal processing hardwares (such as DSP and FPGA).  As a result, a significant research effort has been recently devoted to developing efficient and low-complexity 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 [@Moshavi1996; @muller2001design; @Sessler2005; @Zarei; @Bjornson; @2017Benzin] 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 users channels.  In particular, as seen in [@Zarei; @Bjornson; @2017Benzin], they require the knowledge of each users channel covariance matrix.  This may not be practical, especially in massive MIMO regime where the covariance matrix of each user channel vector has very large dimension $M \times M$ and needs to be estimated and tracked (i.e., updated over time) from the user channel measurements.  Estimation and tracking of the users channel covariance matrix in the massive MIMO regime has also been intensely investigated in recent works (e.g., see [@haghighatshoar2017massive; @haghighatshoar2016low]) 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 [@uplink_MMSE_reduced_rank].  In a recent line of work [@ghods2017optimally], the *approximate message passing* (AMP) algorithm, initially developed for sparse signal recovery in [@bayati2011dynamics], 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) [@Kaczmarz].  A randomized version of KA was recently proposed and analyzed for solving consistent OD SLE [@Strohmer2006].  With the advent of stochastic gradient techniques in machine learning, KA has been revitalized [@needell2014stochastic] and applied to other problems such as solving quadratic equations [@chi2016kaczmarz].  Our work in this paper is mainly based on the randomized KA proposed in [@Strohmer2006].  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 a computational viewpoint, while it does not require the knowledge of the statistics of the users channel vectors. 

Notation
--------

We represent scalar constants by non-boldface letters [(e.g., $x$ or $X$)]{}, sets by calligraphic letters (e.g., $\Xc$), vectors by boldface small letters (e.g., $\xv$), and matrices by boldface capital letters (e.g., $\Xm$).  We denote the $i$-th row and the $j$-th column of a matrix $\bfX$ with the row-vector $\Xm_{i,.}$ and the column-vector $\Xm_{.,j}$ respectively.  We represent the components of a vector $\bfx=(x_1, \dots, x_n)^\transp$ by $x_i$, $i=1,2, \dots, n$. We denote a sequence of vectors by upper/lower indices, e.g., $\bfx^1, \bfx^2, \cdots$ or $\bfx_1, \bfx_2, \cdots$. We indicate the Hermitian conjugate and the transpose of a matrix $\bfX$ by $\bfX^\herm$ and $\Xm^\transp$ respectively, with the same notation being used for vectors and scalars.  We denote the inner product between two matrices (and similarly two vectors) $\bfX$ and $\bfY$ by $\inp{\bfX}{\bfY}=\trace(\bfX^\herm \bfY)$. We use $\|\bfX\|_\sfF=\sqrt{\inp{\bfX}{\bfX}}$ for the Frobenius norm of a matrix $\bfX$ and $\|\bfx\|$ for the $l_2$-norm of a vector $\bfx$. For two $m\times n$ and $m'\times n$ matrices $\bfX$ and $\bfY$, we denote by $[\bfX;\bfY]$ the $(m+m')\times n$ matrix obtained by stacking the rows of $\bfX$ on top of the rows of $\bfY$. The identity matrix of order $p$ is represented by $\bfI_p$. For an integer $k>0$, we use the shorthand notation $[k]$ for $\{1,2,\dots, k\}$.

Problem Statement {#sec:definition}
=================

Consider a multiuser massive MIMO system with a BS with $M\gg 1$ antennas serving $K$ single-antenna users.  We denote the *true* channel vector of user $i \in [K]$ to the $M$ BS antennas with $\bfh^i \in \bC^M$, where $\bfh^i=(h^i_1, \dots, h^i_M)^\transp$ with $h^i_k$ denoting the complex channel gain to the $k$-th antenna.  We denote the *true* $M\times K$ channel matrix consisting of the channel vectors of all the users by $\bfH=[\bfh^1, \dots, \bfh^K]$. We also denote a possibly quantized and noisy *estimate* of this matrix available at the base station (BS) by the $M\times K$ matrix $\bfQ$. It is worthwhile to mention that, in this paper, we always assume that all the computations done at the BS are based on the estimate matrix $\bfQ$. In [*time-division duplexing*]{} (TDD) systems, the matrix $\bfQ$ is obtained during the training period of each channel coherence block[^2] by transmitting orthogonal pilots from the users to the BS in the UL (exploiting UL-DL channel reciprocity [@Marzetta-TWC10; @shepard2012argos]).  In [*frequency-division duplexing*]{} (FDD) systems, in contrast, the BS obtains the estimated channel matrix $\bfQ$ by pilot transmission in the DL and some form of analog or quantized feedback from the users to the BS in the UL (see, e.g., [@caire2010multiuser]). 

Detection and precoding vectors in the UL and the DL {#sec:up}
----------------------------------------------------

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: $$\label{eq:ul1}
\bfy=\bfH\bfs+\bfn,$$ where $\bfs=(s_1, \dots, s_K)^\transp \in \bC^K$ is the vector of symbols sent by the $K$ users, where $\bfH$ is the true channel matrix of the users of size $M \times K$ as before, and where $\bfn \sim \cg(0, \sigma^2 \bfI_M)$ is the additive zero-mean white Gaussian noise (AWGN) at the BS receiver.  Without loss of generality, we assume that the users symbols have mean zero and variance $\bE[|s_i|^2] = P$, and incorporate the pathloss effects due to possibly different distance of the users form the BS as part of the channel vectors $\{\bfh^k: k \in [K]\}$. As stated before, we will consider only the case where the BS detects symbols $\bfs$ from the observations $\bfy$ 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: $$\begin{aligned}
\hat{\bfs}_\text{ZFD}&=(\bfQ^\herm\bfQ)^{-1}\bfQ^\herm\bfy\label{eq:ulZF},\\
\hat{\bfs}_\text{MMSED}&=(\bfQ^\herm\bfQ+\frac{1}{\SNR}\bfI_K)^{-1}\bfQ^\herm\bfy,\label{eq:ulMMSE}\end{aligned}$$ where $\snr=\frac{P}{\sigma^2}$ 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\left( MK^2\right)$ operations, and matrix inversion which requires approximately $O\left(K^3\right)$ operations.  Note that as an alternative to ZFD and MMSED in and respectively, one can apply Maximum Ratio Combining (MRC) $\hat{\bfs}=\bfQ^\herm\bfy$ for detecting the users symbols in the UL.  MRC requires only matrix-vector multiplication and, hence, needs roughly $O(MK)$ operations, which is much less than that needed for ZFD or MMSED.  However, in cases where high spectral efficiency is required (high-SNR conditions) or where inter-cell interference is not the dominating effect, MRC performs much worse than ZFD and MMSED in terms of achievable per-user rates, unless the number of antennas is very large [@RZF_Hoydis; @Huh11; @ZF_MRC_MM] [^3].  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 [@downlink_MM_linear].  The precoding vectors in the DL can be computed according to [@RZF_Hoydis; @RZF_Peel] as: $$\begin{aligned}
\mathbf{x}_{\text{ZFBF}}&=\beta_\text{ZFBF}\bfQ(\bfQ^\herm \bfQ)^{-1}\bfs\label{eq:dlZF}\\
\bfx_{\text{RZFBF}}&=\beta_\text{RZFBF} \bfQ(\bfQ^\herm \bfQ+\xi \bfI)^{-1}\bfs\label{eq:dlMMBF},\end{aligned}$$ where $\beta_\text{ZFBF}$ and $\beta_\text{RZFBF}$ are scaling parameters to impose a given BS average transmit power, and where $\xi$ is the regularization parameter [@RZF_Peel].  Each user $k \in [K]$ will then receive the following signal $$y_k={\bfh^k}^\herm\bfx+n_k,$$ where $\bfx$ can be either $\bfx_{\text{ZF}}$ or $\bfx_{\text{RZFBF}}$, and where $n_k\sim\cg(0, \sigma^2)$ is the AWGN at the user.  Computing precoding vectors of the ZFBF and RZFBF in and respectively has the same order of complexity as in the corresponding UL scenario and . 

Low-Complexity Beamforming via Kaczmarz Algorithm {#sec:KA_intro}
-------------------------------------------------

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., $\hat{\bfs}_\text{ZFD}, \hat{\bfs}_\text{MMSED}$ in and , or finding the suitable $M$-dim precoding signal to be the transmitted to the users in the DL, i.e., $\mathbf{x}_{\text{ZFBF}}, \bfx_{\text{RZFBF}}$ in and , can be posed as finding the solution $\bfw^*$ of an SLE of the form $\bfA \bfw=\bfb$ where $\bfA$ depends on the estimated channel state $\bfQ$ and $\bfb$ depends on the users signal $\bfs$ in a DL and on the noisy received signal from the users $\bfy$ as in in a UL scenario.  We find the optimal solution $\bfw^*$ of $\bfA \bfw=\bfb$ via an iterative procedure known as Kaczmarz Algorithm (KA).  This technique was initially proposed by Kaczmarz for efficiently solving a consistent overdetermined (OD) SLE [@Kaczmarz].  KA works as follows.  At the start of each iteration $t$, the algorithm has an estimate $\bfw^t$ of the optimal solution $\bfw^*$ and selects one of the rows, say row number $r(t)$, of $\bfA$. If the estimate $\bfw^t$ satisfies the $r(t)$-th equation, namely, if $\inp{\bfa_{r(t)}}{\bfw^t}= b_{r(t)}$, where $\bfa_{r(t)}=(\bfA_{r(t),.})^\herm$ denotes the conjugate of the $r(t)$-th row of $\bfA$, then KA keeps $\bfw^t$ as it is.  Otherwise, it updates $\bfw^t$ along $\bfa_{r(t)}$ to make the $r(t)$-th equation consistent.  This is summarized in Algorithm \[tab:Kac\].  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 $\bfw$ (in $l_2$-norm) to the current solution $\bfw^t$ that satisfies the selected equation or mathematically speaking $\bfw^{t+1}=\argmin_{\bfw} \|\bfw - \bfw^t\|^2 \text{ s.t. }  \inp{\bfa_{r(t)}}{\bfw}= b_{r(t)}$.

Initialize $\bfw^0$. Update $\bfw^t$ as $\bfw^{t+1}=\bfw^t+\frac{ b_{r(t)}- \inp{\bfa_{r(t)}}{ \bfw^t} }{\| \bfa_{r(t)}\|^2}\bfa_{r(t)}$.

The convergence speed of KA is known to depend highly on the order of selecting the rows of $\bfA$, 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 $\bfA$ are selected randomly according to some probability distribution that might depend on $\bfA$. Such a randomized version of KA was recently proposed and analyzed by Strohmer and Vershynin in [@Strohmer2006], where they showed that the performance of the resulting randomized KA depends on the condition number of the matrix $\bfA$. It was also shown in [@Strohmer2006] 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 [@Strohmer2006] 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 [@Needell2010].  We investigate the effect of this residual term in the massive MIMO setup in Section \[sec:results:residual\] via numerical simulations, where we illustrate that it generally incurs a considerable loss in per-user rate.  Also, to avoid this effect, we develop a new variant of KA that always solves consistent equations after a suitable transformation, where we also show that the resulting updates of the new algorithm converge to the desired least-squares solution.  Our method resembles [@extendedKaczmarz], which uses a step prior to KA to remove the inconsistency of the equations but has much less computational complexity. 

Unified Mathematical Framework {#sec:math frame}
==============================

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 $\bfA \bfx=\bfb$ where $\bfA\in \bC^{m\times n}$ is the matrix of linear equations, where $\bfx \in \bC^n$ is the vector of unknown values, and where $\bfb\in \bC^m$ is the vector of known coefficients.  We say that the SLE $\bfA \bfx=\bfb$ is *consistent* if it has a solution.  We always assume that $\bfA$ has full row-rank in the UD case ($m<n$) and full column-rank in the OD case ($m\geq n$).  As we will see, these conditions are satisfied (with probability $1$) in massive MIMO applications studied in this paper. 

We will consider the randomized version of KA explained in Algorithm\[tab:Kac\], where at each iteration the algorithm selects a row of $\bfA$ randomly and independently from the previous iterations and according to a given probability distribution $\bfp=(p_1, \dots, p_m)^\transp$, where $p_i\in [0,1]$ is the probability of selecting row $i\in [m]$ and where $\sum_{i\in [m]}p_i=1$. We denote the random rank-1 projection operator produced by the random row selection by $\scrP_R=\frac{\bfa_R \bfa_R^\herm}{\|\bfa_R\|^2} \in \bC^{n\times n}$, where $R \in [m]$ is a random variable denoting the index of the selected row, where $\bP[R=i]=p_i$ and where we denote by $\bfa_i=(\bfA_{i,.})^\herm \in \bC^m$ the conjugate of the $i$-th row of $\bfA$. We define the average of the projection operator $\scrP_R$ as $$\begin{aligned}
\label{avg_proj}
\overline{\scrP}=\bE_R[\scrP_R]=\sum_{i=1}^m p_i \frac{\bfa_i \bfa_i^\herm}{\|\bfa_i\|^2}.\end{aligned}$$ Note that $\overline{\scrP}$ is an $n \times n$ *positive semi-definite* (PSD) matrix, with $\tr(\overline{\scrP})=\sum_{i \in [m]}p_i=1$. As we will see in the next sections:

\(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 $\clX \subset \bC^n$ generated by conjugate of the rows of $\bfA$, i.e., the column span of $\bfA^\herm$.

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

\[avg\_gain\] Let $\bfA$ be an $m\times n$ matrix.  Let $\bfp \in \bR_+^m$ be a probability distribution over the rows of $\bfA$ and define the average projection operator $\overline{\scrP}$ produced by the rows of $\bfA$ as in .  Let $\clX$ be the subspace produced by the span of columns of $\bfA^\herm$. The average gain of the matrix $\bfA$ over the subspace $\clX$ is defined as: $$\begin{aligned}
    \kappa_{\clX}(\bfA, \bfp)=\min_{\bfx\in \clX, \bfx \neq 0} \frac{\bfx^\herm \overline{\scrP} \bfx}{\|\bfx\|^2}. 
    \label{eq:opt_gain}
    \end{aligned}$$

Note that $\kappa_{\clX}(\bfA, \bfp)$ depends both on the matrix $\bfA$ and the probability distribution $\bfp$. Moreover, as $\tr(\overline{\scrP})=1$ and as the subspace $\clX$ has dimension $\min\{m,n\}$ and $\overline{\scrP}$ acts on this subspace, we have that $\kappa_{\clX}(\bfA, \bfp) \in [0,\frac{1}{\min\{m,n\}}]$. We will show that the closer $\kappa_\clX(\bfA, \bfp)$ is to $\frac{1}{\min\{m,n\}}$, the faster KA converges across the iterations.  Note that from , for a given matrix $\bfA$, $\kappa_{\clX}(\bfA, \bfp)=\lambda_{\min}(\overline{\scrP}| \clX)$ is given by the minimum eigen-value of $\overline{\scrP}$ over the subspace $\clX$. As $\overline{\scrP}$ is a linear function of $\bfp$ and as $\lambda_{\min}(\overline{\scrP}| \clX)$ is a concave function of $\overline{\scrP}$ [@tao2012topics], $\kappa_{\clX}(\bfA, \bfp)$ is a concave function of $\bfp$ over the convex set of probability distributions over the row set of $\bfA$. Therefore, finding the best distribution $\bfp$ maximizing $\kappa_{\clX}(\bfA, \bfp)$ 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 addressed in this paper.  This might be, however, very useful for situations where one solves many equations with the same matrix $\bfA$, whereby the complexity of computing the optimal $\bfp$ is amortized over time.  This is not the case in massive MIMO applications since the matrix $\bfA$ depends on the estimated channel matrix $\bfQ$, which changes from one coherence time to another.  Instead, we will use the suboptimal probability distribution $\bfp$ proposed by [@Strohmer2006], where $p_i=\frac{\|\bfa_i\|^2}{\|\bfA\|_\sfF^2}$ and the probability of selecting each row $i\in [m]$ scales proportionally to its squared $l_2$-norm $\|\bfa_i\|^2$. 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 $\bfp$. Replacing this probability distribution in results in the following average operator: $$\begin{aligned}
\label{avg_proj2}
\overline{\scrP}&=\bE_R[\scrP_R]=\sum_{i=1}^m p_i \frac{\bfa_i \bfa_i^\herm}{\|\bfa_i\|^2}=\sum_{i=1}^m \frac{\|\bfa_i\|^2}{\|\bfA\|_\sfF^2} \frac{\bfa_i \bfa_i^\herm}{\|\bfa_i\|^2}=\sum_{i=1}^m \frac{\bfa_i \bfa_i^\herm}{\|\bfA\|_\sfF^2}=\frac{\bfA^\herm \bfA}{\|\bfA\|_\sfF^2}.\end{aligned}$$ From Definition \[avg\_gain\], this yields the following gain parameter for the matrix $\bfA$, which will be used in the rest of the paper. 

\[min\_gain\] Let $\bfA$ be the $m\times n$ matrix as defined before and let $\clX$ be a linear subspace of $\bC^n$ produced by the column span of $\bfA^\herm$. We define the normalized minimum gain of the matrix $\bfA$ along the subspace $\clX$ as: $$\begin{aligned}
    \kappa_\clX(\bfA)=\min_{\bfx\in \clX, \bfx \neq 0} \frac{\|\bfA \bfx\|^2}{\|\bfA\|_\sfF^2 \|\bfx\|^2},
    \end{aligned}$$ $\lozenge$

Note that $\kappa_\clX(\bfA) \in [0,\frac{1}{\min\{m,n\}}]$ as before, and the closer $\kappa_\clX(\bfA)$ becomes to $\frac{1}{\min\{m,n\}}$, the faster KA converges across the iterations.  Fig.\[fig:opt\_subopt\] compares the sub-optimal gain $\kappa_\clX(\bfA)$ with the best average gain $\kappa_{\clX}(\bfA, \bfp^*)$ obtained by the optimal probability distribution $\bfp^*$. We consider a massive MIMO scenario where the matrix $\bfA$ has i.i.d.  $\cg(0,1)$ components with $M=256$ rows corresponding to BS antennas but different number of columns $K\in \{10,15,25\}$ corresponding to the number of served users.  KA for such a matrix $\bfA$ arises in an UL scenario where the users have spatially white channel vectors.  From well-known results in random matrix theory [@random_mtx_verdu], it is not difficult to show that with a very high probability (on the realization of the matrix $\bfA$) $$\begin{aligned}
\kappa_\clX(\bfA)\approx\frac{(\sqrt{M}-\sqrt{K})^2}{MK}\label{dumm_rev}
=\frac{(1-\sqrt{\rho})^2}{K},\end{aligned}$$ where $\rho:=\frac{K}{M}$ denotes the ratio between the number of served users $K$ and the number of BS antennas $M$, typically called the *loading factor* of the massive MIMO system.  Considering the fact that $\kappa_{\clX}(\bfA, \bfp^*)\leq \frac{1}{\min\{M,K\}}=\frac{1}{K}$, it is seen from that, for a fixed loading factor $\rho$, the matrix gain $\kappa_\clX(\bfA)$ gives an approximation of the best gain $\kappa_{\clX}(\bfA, \bfp^*)$ up to a multiplicative factor (see also Fig.\[fig:opt\_subopt\]). 

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. 

\[main\_thm\] Let $\bfA$ be an $m \times n$ matrix and let $\bfb \in \bC^m$. Let $\clX$ be the subspace of $\bC^n$ generated by the column span of $\bfA^\herm$ and suppose that $\bfA \bfx = \bfb$ has a solution $\bfx^* \in \clX$. Let $\bfx^0 \in \clX$ be an arbitrary initialization and let $\bfx^t$ be the estimate at time $t$ of the randomized KA with the suboptimal distribution as explained before starting from $\bfx^0$. Then, $$\begin{aligned}
        \bE\big [\|\bfx^t - \bfx^*\|^2\big ] \leq (1-\kappa_\clX(\bfA))^t \|\bfx^0 - \bfx^*\|^2,
        \end{aligned}$$ where the expectation is taken over the randomized row selection in KA, and where $\kappa_\clX(\bfA)$ denotes the minimum gain of the linear operator $\bfA$ over the subspace $\clX$ as in Definition \[min\_gain\].  $\square$

From the update equation of KA, we have: $$\begin{aligned}
            \bfx^{t+1}- \bfx^*&=\bfx^t - \frac{\inp{\bfa_{r(t)}}{\bfx^t} - b_{r(t)}}{\|\bfa_{r(t)}\|^2} \bfa_{r(t)} - \bfx^*=\bfx^t - \bfx^* - \frac{\inp{\bfa_{r(t)}}{\bfx^t - \bfx^*}}{\|\bfa_{r(t)}\|^2} \bfa_{r(t)}\nonumber\\
            &=\big(\bfI_n - \frac{\bfa_{r(t)} \bfa_{r(t)}^\herm }{\|\bfa_{r(t)}\|^2}\big) (\bfx^t - \bfx^*)=(\bfI_n-\scrP_{r(t)})(\bfx^t - \bfx^*),
            \end{aligned}$$ where $r(t) \in [m]$ denotes the index of the random row selected by KA at iteration $t$, where $b_{r(t)}$ is the component of $\bfb$ at position $r(t)$, where we used the fact that $\bfx^*$ is a solution of $\bfA \bfx=\bfb$, thus, $b_{r(t)}=\inp{\bfa_{r(t)}}{\bfx^*}$, and where $\scrP_{r(t)}:=\frac{\bfa_{r(t)} \bfa_{r(t)}^\herm}{\|\bfa_{r(t)}\|^2}$ is the random rank-1 projection operator onto the subspace spanned by $\bfa_{r(t)}$ as defined before.  It is not difficult to check that $\bfI_n-\scrP_{r(t)}$ is also a projection operator, thus, $(\bfI_n-\scrP_{r(t)})^\herm (\bfI_n-\scrP_{r(t)})=(\bfI_n-\scrP_{r(t)})$. This yields: $$\begin{aligned}
            \|\bfx^{t+1} - \bfx^*\|^2=(\bfx^{t}-\bfx^*)^\herm (\bfI_n-\scrP_{r(t)}) (\bfx^{t}-\bfx^*). 
            \end{aligned}$$ Note that $\bfx^{t}-\bfx^*$ is independent of $\bfI_n-\scrP_{r(t)}$ since it depends on the random selection of the rows at iterations before $t$. Using this independence, conditioning on all row selections before $t$, and taking the expected value over $r(t)$ yields: $$\begin{aligned}
            \bE\Big [&\|\bfx^{t+1}-\bfx^*\|^2 \big |r(1), \dots, r(t-1)\Big ]
            =(\bfx^t - \bfx^*)^\herm \bE[\bfI_n - \scrP_{r(t)}] (\bfx^t - \bfx^*)\\
            &=(\bfx^t - \bfx^*)^\herm (\bfI_n - \overline{\scrP}) (\bfx^t - \bfx^*)
            = \|\bfx^t - \bfx^*\|^2 - (\bfx^t - \bfx^*)^\herm \overline{\scrP} (\bfx^t - \bfx^*)\\
            &\stackrel{(a)}{\leq} (1-\kappa_\clX(\bfA))\|\bfx^t - \bfx^*\|^2,
            \end{aligned}$$ where in $(a)$ we used the fact that the minimum gain of $\overline{\scrP}$ for the mentioned suboptimal distribution is given by $\kappa_\clX(\bfA)$ as in Definition \[min\_gain\].  Finally, taking expected value with respect to $r(1), \dots, r(t-1)$ and applying the induction we obtain: $$\begin{aligned}
            \bE\big[\|\bfx^{t+1}-\bfx^*\|^2\big]& \leq \big(1-\kappa_\clX(\bfA) \big)\bE\big[\|\bfx^{t}-\bfx^*\|^2\big] \leq \big(1-\kappa_\clX(\bfA) \big)^{t+1} \|\bfx^0 - \bfx^*\|^2,
            \end{aligned}$$ where we used the fact that the initial point $\bfx^0$ is selected deterministically and dropped the expectation.  This completes the proof. 

From Theorem \[main\_thm\], it is seen that the convergence speed of KA depends on the parameter $\kappa_\clX(\bfA)$, where a larger $\kappa_\clX(\bfA)$ guarantees a faster convergence.  Moreover, since $\kappa_\clX(\bfA) \leq \frac{1}{\min\{m,n\}}$, at least $O(\min\{m,n\})$ iterations are needed for the convergence.  $\lozenge$

Since the random variable $\|\bfx^t - \bfx^*\|^2$ is non-negative, from the well-known Markov’s inequality [@grimmett2001probability], its expected value also provides a probabilistic bound on its tail behavior, i.e., for any $\zeta>0$, we have: $$\begin{aligned}
                    \bP\big [\|\bfx^t - \bfx^*\|^2 \geq \zeta\big] &\leq \frac{\bE\big[\|\bfx^{t}-\bfx^*\|^2\big]}{\zeta} \leq \big(1-\kappa_\clX(\bfA) \big)^t \frac{\|\bfx^0 - \bfx^*\|^2}{\zeta},
                    \end{aligned}$$ where $\kappa_\clX(\bfA)$ is as defined before.  $\lozenge$

In the following, we will design several variants of KA depending on the specific beamforming scenario and apply Theorem \[main\_thm\] to analyze the convergence speed of the resulting algorithm.  We will also explain how the convergence speed depends on the system 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. 

Kaczmarz Algorithm for different Beamforming Scenarios {#sec:scenarios}
======================================================

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 \[sec:definition\] that we denote the perfect channel state of the users by an $M\times K$ matrix $\bfH$ and the noisy and possibly quantized estimate of $\bfH$ available at the BS by the $M\times K$ matrix $\bfQ$.

MMSE/ZF detection in the UL {#sec:mmse_ul}
---------------------------

We first consider a UL scenario.  From what we have previously discussed in Section \[sec:up\], the ZFD in and MMSED in can be derived, by setting $\xi=0$ and $\xi=\frac{1}{\snr}$ respectively, from the following general expression: $$\begin{aligned}
\label{MMSED_eq}
\widehat{\bfs}=(\bfQ^\herm \bfQ + \xi \bfI_K)^{-1} \bfQ^\herm \bfy.\end{aligned}$$ Let us now focus on the general case where $\xi\ne0$. We make the following observation, that we prove here for the sake of completeness. 

\[prop:mmse\_bf\_ul\] Let $\bfQ$, $\bfy$ and $\widehat{\bfs}$ be as in .  Then, the signal estimate $\widehat{\bfs}$ is the optimal solution of the optimization problem $\argmin_{\bfw \in \bC^K} \|\bfQ \bfw - \bfy\|^2 + \xi \|\bfw\|^2$. $\square$

Taking the derivative of the cost function gives $2\bfQ^\herm (\bfQ \bfw - \bfy) + 2\xi \bfw$. Setting this equal to zero yields the optimal solution $\bfw^*=(\bfQ^\herm \bfQ + \xi \bfI_K)^{-1} \bfQ^\herm \bfy$ which coincides with $\widehat{\bfs}$.

Note that we can write the $l_2$ regularized least-squares cost function $\|\bfQ \bfw - \bfy\|^2 + \xi \|\bfw\|^2$ in Proposition \[prop:mmse\_bf\_ul\] as $\|\bfB \bfw - \bfy_0\|^2$, where $\bfB=[\bfQ; \sqrt{\xi} \bfI_K]$ is an $(M+K)\times K$ matrix and where $\bfy_0=[\bfy; 0]$ is the $(M+K)$-dim vector obtained by appending a $K$-dim zero vector to the noisy received signal $\bfy$. The signal estimate $\widehat{\bfs}$ in can also be written as the least-squares solution of the newly defined cost function and is given by $$\begin{aligned}
\label{s_eq_mmse}
\widehat{\bfs}=(\bfB^\herm \bfB)^{-1} \bfB^\herm \bfy_0= (\bfQ^\herm \bfQ + \xi \bfI_K)^{-1} \bfQ^\herm \bfy.\end{aligned}$$

Let us consider the SLE $\bfB \bfw =\bfy_0$. Note that this SLE is OD and is always inconsistent when $\xi>0$. Furthermore, it is generally inconsistent even when $\xi=0$ due to the presence of noise in the received measurements $\bfy$. Thus, a direct application of the randomized KA will result in a residual error as we discussed in Section \[sec:KA\_intro\] and as we will illustrate via numerical simulations in Section \[sec:results:residual\].  To avoid this, we solve this SLE in two steps.  We first define a new vector $$\begin{aligned}
\label{y0_MMSED}
\widehat{\bfy}_0= \bfB \widehat{\bfs}\stackrel{(i)}{=}\bfB (\bfB^\herm \bfB)^{-1}\bfB^\herm \bfy_0=\bfB (\bfB^\herm \bfB)^{-1}\bfQ^\herm \bfy,\end{aligned}$$ where in $(i)$ we used .  Note that $\widehat{\bfy}_0$ lies in the subspace $\clX$ of $\bC^{M+K}$ spanned the columns of $\bfB$. Moreover, from $$\begin{aligned}
\label{mmse_eq_ul_dumm}
\bfB^\herm \widehat{\bfy}_0=(\bfB^\herm \bfB) (\bfB^\herm \bfB)^{-1}\bfQ^\herm \bfy=\bfQ^\herm \bfy,\end{aligned}$$ it is seen that $\widehat{\bfy}_0$ satisfies the UD SLE $\bfB^\herm \bfz=\bfb$ for the unknown $\bfz \in \bC^{M+K}$ and the set of known coefficients $\bfb=\bfQ^\herm \bfy$. We first apply the randomized KA to this SLE to recover $\widehat{\bfy}_0$. We then have the following result. 

\[prop:MMSED\_ul\_2\] Let $\bfz^t$ be the estimate obtained at iteration $t$ of KA applied to the UD SLE $\bfB^\herm \bfz=\bfb$ where $\bfb=\bfQ^\herm \bfy$. Assume that KA starts from the zero initialization $\bfz^0=0$ and let $\widehat{\bfy}_0$ be as in .  Then $$\begin{aligned}
    \bE\big [\|\bfz^t - \widehat{\bfy}_0\|^2 \big ] \leq (1-\kappa_\clX(\bfB^\herm))^t \|\widehat{\bfy}_0\|^2,
    \end{aligned}$$ where $\kappa_\clX(\bfB^\herm)$ is the minimum gain of the matrix $\bfB^\herm$ over the subspace $\clX$ generated by the columns of $\bfB$. $\square$

From the update equation for KA and applying a simple induction, it is seen that starting from $\bfz^0=0$, all the estimates $\bfz^t$ produced by KA lie in the subspace $\clX$ produced by the columns of $\bfB$. Thus, the result simply follows by applying Theorem \[main\_thm\]. 

Once $\widehat{\bfy}_0$ is recovered, we solve the OD SLE $\bfB \widehat{\bfs}=\widehat{\bfy}_0$ in to find the estimate $\widehat{\bfs}$. In contrast with the initial SLE $\bfB \widehat{\bfs}=\bfy$, which was inconsistent, the OD SLE in is always consistent since $\widehat{\bfy}_0$ always lies in the subspace spanned by the columns of $\bfB$. Note that since $\bfB \widehat{\bfs}=\widehat{\bfy}_0$ is consistent and $\bfB=[\bfQ; \sqrt{\xi}\bfI_K]$ has $\sqrt{\xi}\bfI_K$ as its submatrix, $\widehat{\bfs}$ is simply given by the last $K$ components of $\widehat{\bfy}_0$ divided by $\sqrt{\xi}$. As a result, there is no need to solve the second SLE $\bfB \widehat{\bfs}=\widehat{\bfy}_0$ per se.  Moreover, denoting by $\bfz^t$ the estimate of KA for the first SLE $\bfB^\herm \bfz=\bfb$ with $\bfb=\bfQ^\herm \bfy$, we can obtain an estimate $\widehat{\bfs}^t$ at each iteration $t$ of KA by dividing the last $K$ components of $\bfz^t$ by $\sqrt{\xi}$. By introducing $\bfu^t\in \bC^M$ and $\bfv^t\in \bC^K$ such that $\bfz^t=[\bfu^t; \sqrt{\xi} \bfv^t]$, we can summarize KA for $\bfB^\herm \widehat{\bfy}_0=\bfQ^\herm \bfy$ as in Algorithm \[tab:Kac\_MMSED\].  A direct inspection shows that, due to the submatrix $\sqrt{\xi} \bfI_K$ in $\bfB=[\bfQ; \sqrt{\xi} \bfI_K]$, at each iteration of KA only a single randomly selected component of the estimate $\bfv^t$ is updated, where $\bfv^t$ converges to the least-squares estimate of the symbols $\widehat{\bfs}$ corresponding to the MMSED.  In particular, due to the continuity of $\widehat{\bfs}$ with respect to $\xi$ as $\xi \to 0$, one can obtain the ZFD estimate $\widehat{\bfs}$ by setting $\xi=0$ in Algorithm \[tab:Kac\_MMSED\]. 

$\bfQ\in \bC^{M \times K}$, received noisy UL signal $\bfy \in \bC^M$, and MMSED parameter $\xi\geq 0$. .  Compute $\bfb=\bfQ^\herm \bfy \in \bC^K$. Define $\bfu^t \in \bC^M$ and $\bfv^t \in \bC^K$ with $\bfu^0=0$, $\bfv^0=0$. Pick a row $r(t)$ of $\bfQ^\herm$ denoted by $\bfq_{r(t)}^\herm$ with a probability proportional to $\|\bfq_{r(t)} \|^2+ \xi$ given by $\frac{\|\bfq_{r(t)} \|^2+ \xi}{\|\bfQ\|_\sfF^2 + K \xi}$. Compute the residual $\gamma^t:=\frac{ b_{r(t)}- \inp{\bfq_{r(t)}}{\bfu^t} - \xi\, v^{t}_{r(t)}}{\|  \bfq_{r(t)}\|^2 + \xi}$. Update $\bfu^{t+1}=\bfu^t + \gamma^t \bfq_{r(t)}$. Update $v_{r(t)}^{t+1}=v^t_{r(t)}+ \gamma^t$, and $v_{j}^{t+1}=v^t_{j}$ for $j \neq r(t)$. Set $\bfw=\bfv^{T-1}$. Set the signal estimate to $\widehat{\bfs}=\bfw$.

To analyze the convergence speed of KA, from Proposition \[prop:MMSED\_ul\_2\], we need to compute the gain of the matrix $\bfB^\herm$ along the subspace $\clX$ spanned by the conjugate of the rows of $\bfB^\herm$, that is, the column span of $\bfB$. Since every vector in $\clX$ can be written as $\bfB \bfw$ for some $\bfw\in \bC^K$, a direct calculation shows that $\kappa_\clX(\bfB^\herm)$ is given by: $$\begin{aligned}
\kappa_\clX(\bfB^\herm)&=\min_{\bfx \in \clX, \bfx\neq 0} \frac{\|\bfB^\herm \bfx\|^2}{\|\bfB^\herm\|_\sfF^2 \|\bfx\|^2}=\min_{\bfw \in \bC^K, \bfw\neq 0} \frac{\|\bfB^\herm \bfB \bfw\|^2}{\|\bfB^\herm\|_\sfF^2 \|\bfB \bfw\|^2}\nonumber\\
&=\frac{\lambda_{\min} (\bfB^\herm \bfB)}{\|\bfB\|_\sfF^2}=\frac{\lambda_{\min} (\bfQ^\herm \bfQ) + \xi}{ \|\bfQ\|_\sfF^2 + K\xi}.\label{MMSED_DL}\end{aligned}$$ Note that $\frac{\lambda_{\min} (\bfQ^\herm \bfQ) }{ \|\bfQ\|_\sfF^2} \leq \frac{1}{K}$ and it is not difficult to check that for any $\xi\ge0$, we have $$\begin{aligned}
 \frac{\lambda_{\min} (\bfQ^\herm \bfQ) }{ \|\bfQ\|_\sfF^2} \leq \frac{\lambda_{\min} (\bfQ^\herm \bfQ) + \xi}{ \|\bfQ\|_\sfF^2 + K\xi}  \leq \frac{1}{K}.\end{aligned}$$ This implies that the algorithm converges faster for larger values of $\xi\ge0$ (i.e., smaller UL SNR) and the slowest convergence corresponds to the ZFD with $\xi=0$. When the channel vectors of the users have white complex Gaussian components (no spatial correlation), from well-known results in random matrix theory [@random_mtx_verdu], it results that $$\begin{aligned}
\kappa_\clX(\bfB^\herm)&\geq  \frac{\lambda_{\min} (\bfQ^\herm \bfQ)}{\|\bfQ\|_\sfF^2}= \frac{(\sqrt{M}-\sqrt{K})^2}{MK}
=\frac{(1-\sqrt{\rho})^2}{K}.\end{aligned}$$ Therefore, from Proposition \[prop:MMSED\_ul\_2\] it results that $O(\frac{K}{(1-\sqrt{\rho})^2})$ iterations are sufficient to reasonably estimate the symbol sequence $\widehat{\bfs}$ for the MMSED for all ranges of $\xi> 0$ and for the ZFD for $\xi=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 complexity scales as $O(\frac{KM}{(1-\sqrt{\rho})^2})$. To this we should also add $O(KM)$ operations needed for computing the sub-optimal row-sampling probability distribution and $O(KM)$ operations needed for calculating $\bfQ^\herm \bfy$ in Algorithm \[tab:Kac\_MMSED\]. 

RZF/ZF Precoding in the DL {#sec:MMBF_dl}
--------------------------

In this section, we consider RZF precoding in the DL.  Suppose $\bfs=(s_1, \dots, s_K)^\transp$ 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 $\bfy:=\bfQ(\bfQ^\herm \bfQ + \xi \bfI_K)^{-1} \bfs$ via its $M$ antennas.  We observe that $\bfy$ lies in the subspace spanned by the columns of $\bfQ$. We design a new KA that finds the $K$-dim signal $\bfw=(\bfQ^\herm \bfQ + \xi \bfI_K)^{-1}\bfs$, from which we can obtain $\bfy=\bfQ \bfw$ via a matrix-vector multiplication.  We first prove the following result. 

\[prop:mmbf\_dl\] Let $\bfQ$, $\bfs$ and $\bfw=(\bfQ^\herm \bfQ + \xi \bfI_K)^{-1}\bfs$ be as introduced before.  Then, $\bfw$ is the optimal solution of the following optimization problem $\bfw=\argmin_{\bfx \in \bC^K}  \|\bfQ \bfx\|^2+ \xi \|\bfx - \bfs_\xi\|^2 $, where $\bfs_\xi=\frac{\bfs}{\xi}$. $\square$

Taking the derivative with respect to $\bfx$ and setting it equal to $0$ yields $\bfQ^\herm \bfQ \bfx + \xi \bfx - \bfs=0$. This gives the optimal solution as $\bfx^*=(\bfQ^\herm \bfQ + \xi \bfI_K)^{-1} \bfs$, which coincides with $\bfw=(\bfQ^\herm \bfQ + \xi \bfI_K)^{-1}\bfs$.

Channel state estimate $\bfQ\in \bC^{M \times K}$, users symbols $\bfs \in \bC^K$, and RZFBF parameter $\xi\geq 0$. .  Define $\bfu^t \in \bC^M$ and $\bfv^t \in \bC^K$ with $\bfu^0=0$, $\bfv^0=0$. Pick a row $r(t)$ of $\bfQ^\herm$ denoted by $\bfq_{r(t)}^\herm$ with a probability proportional to $\|\bfq_{r(t)} \|^2+ \xi$ given by $\frac{\|\bfq_{r(t)} \|^2+ \xi}{\|\bfQ\|_\sfF^2 + K \xi}$. Compute the residual $\gamma^t:=\frac{ s_{r(t)}- \inp{\bfq_{r(t)}}{\bfu^t} - \xi\, v^{t}_{r(t)}}{\|  \bfq_{r(t)}\|^2 + \xi}$. Update $\bfu^{t+1}=\bfu^t + \gamma^t \bfq_{r(t)}$. Update $v_{r(t)}^{t+1}=v^t_{r(t)}+ \gamma^t$, and $v_{j}^{t+1}=v^t_{j}$ for $j \neq r(t)$. Set $\bfw=\bfv^{T-1}$. Set the desired DL signal to $\bfy=\bfQ \bfw$.

Proposition \[prop:mmbf\_dl\] implies that we can obtain $\bfw=(\bfQ^\herm \bfQ + \xi \bfI_K)^{-1}\bfs$ as the optimal solution of $\|\bfQ \bfx\|^2 + \xi \|\bfx - \bfs_\xi\|^2$. We first write this as the least-squares minimization $\|\bfB \bfx-\bfb\|^2$ where $\bfB=[\bfQ; \sqrt{\xi}\bfI_K]$ is an $(M+K)\times K$ matrix and where $\bfb=[0;\sqrt{\xi}\bfs_\xi]$ is an $(M+K)$-dim vector with $\sqrt{\xi}\bfs_\xi=\frac{\bfs}{\sqrt{\xi}}$ in its last positions.  The SLE $\bfB \bfw=\bfb$ is OD and should be solved for the vector $\bfw \in \bC^K$, from which we obtain the desired signal $\bfy=\bfQ \bfw$ to be transmitted in the DL.  This SLE is inconsistent unless $\bfs=0$. Thus, we follow the same procedure as in the design of MMSED for the UL.  We first define $\bfz=\bfB \bfw$ and solve the UD SLE $$\begin{aligned}
\label{mmbf_z_dumm}
\bfB^\herm \bfz= \bfB^\herm \bfb=\sqrt{\xi} \frac{\bfs}{\sqrt{\xi}}=\bfs,\end{aligned}$$ for $\bfz$ and then solve the OD but consistent SLE $\bfB \bfw=\widehat{\bfz}$, where $\widehat{\bfz}$ is the estimate of $\bfz$ obtained from KA applied to .  Note that, as in the MMSED, due to the presence of $\sqrt{\xi} \bfI_K$ submatrix in $\bfB=[\bfQ; \sqrt{\xi} \bfI_K]$, the solution of the latter SLE is simply obtained by dividing the last $K$ components of $\widehat{\bfz}$ by $\sqrt{\xi}$. Denoting by $\bfz^t=[\bfu^t; \sqrt{\xi} \bfv^t]$ the estimate at iteration $t$ of KA, we can summarize this in Algorithm \[tab:Kac\_MMBF\]. 

In terms of speed of convergence, a direct calculation shows that the minimum gain of the matrix $\bfB^\herm$ over the subspace $\clX$ spanned by the columns of $\bfB$ is given by: $$\begin{aligned}
\kappa_\clX (\bfB^\herm)&=\frac{ \lambda_{\min} (\bfB^\herm \bfB)}{\|\bfB\|_\sfF^2}=\frac{ \lambda_{\min} (\bfQ^\herm \bfQ) + \xi}{\|\bfQ\|_\sfF^2 + K \xi}.\end{aligned}$$ 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 computational complexity of KA scales as $O(\frac{KM}{(1-\sqrt{\rho})^2})$. To this we should also add $O(KM)$ operations needed for computing the sampling probability distribution and $O(KM)$ operation needed for calculating $\bfQ\bfw$ in Algorithm \[tab:Kac\_MMBF\]. 

A direct inspection in Algorithm \[tab:Kac\_MMSED\] and \[tab:Kac\_MMBF\] also reveals that there is a symmetry between the MMSED for the UL and RZFBF for the DL: In the former, first the noisy received signal $\bfy$ is multiplied by $\bfQ^\herm$ and then the resulting vector $\bfb=\bfQ^\herm \bfy$ enters KA, which yields the estimate of the users symbols $\widehat{\bfs}$, whereas in the latter, first the users symbols $\bfs$ to be transmitted in the DL enters KA and then the output of KA is multiplied by $\bfQ$ 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/multiple-access channel [@vishwanath2003duality]. 

Computation of the Detector/Precoder Matrix {#sec:matrix_design}
===========================================

Up to now, we have used our proposed KA to compute an estimate of the users symbols from the received signal in the UL or to compute a suitable precoding vector in the DL.  More specifically, we did not compute any UL/DL detection/precoding matrix from the available estimate of the channel state directly.  In this section, our goal is to show that our proposed KA can be generalized to compute the detection/precoding matrix as well.  This results in a saving in computational complexity when the channel coherence block is quite large such that the detection/precoding matrix corresponding to each coherence block is computed only once and is then used for all the time-frequency signal dimensions inside the same coherence block.  In the following, we focus on a specific coherence block and denote the estimate of the channel state of the users over this coherence block by the $M\times K$ matrix $\bfQ$ as before.  We consider the following procedure to compute the corresponding detection/precoding matrix. 

In the DL, we use Algorithm \[tab:Kac\_MMBF\] with the input $\bfQ$ and run $K$ KAs in parallel where the input to the $i$-th KA is $\bfs_i=\bfe_i \in \bC^K$, where $\bfe_i$ denotes the $i$-th canonical basis, which has $1$ as its $i$-th component and is $0$ elsewhere.  Note that the $K$ parallel KAs may share their randomness, i.e., they may use the same random row index at each iteration, or they may use their own independent randomness.  Let $\bfw_i \in \bC^K$ denote the output of the $i$-th KA after a suitable convergence.  Then, from the last step of Algorithm \[tab:Kac\_MMBF\], namely, post-multiplication by $\bfQ$, the precoding vector $\bfy_i \in \bC^M$ in the DL produced by the $i$-th KA (corresponding to the symbol vector $\bfs_i=\bfe_i$) would be $\bfy_i=\bfQ \bfw_i$. However, due to the convergence of KA, $\bfy_i$ should be (approximately) equal to the output of the RZF in the DL with the input $\bfs_i=\bfe_i$, i.e., $$\begin{aligned}
\label{dumm_bf_mat1}
\bfy_i=\bfQ(\bfQ^\herm \bfQ + \xi \bfI_K)^{-1} \bfs_i=\bfQ(\bfQ^\herm \bfQ + \xi \bfI_K)^{-1}\bfe_i,\end{aligned}$$ which corresponds to the $i$-th column of the precoding matrix $\bfQ(\bfQ^\herm \bfQ + \xi \bfI_K)^{-1}$. Since $\bfQ$ is assumed to have a full column rank, from and the fact that $\bfy_i=\bfQ \bfw_i$, it immediately results that $\bfw_i$, after a suitable convergence, should correspond to the $i$-th column of $(\bfQ^\herm \bfQ + \xi \bfI_K)^{-1}$. Therefore, $(\bfQ^\herm \bfQ + \xi \bfI_K)^{-1}$ can be well approximated by the $K \times K$ matrix $\bfW=[\bfw_1, \dots, \bfw_K] \in \bC^{K \times K}$ consisting of the output produced by $K$ parallel KAs with corresponding inputs $\{\bfs_i=\bfe_i: i\in [K]\}$. Also, the RZF matrix for the DL $\bfQ(\bfQ^\herm \bfQ + \xi \bfI_K)^{-1}$ can be approximated as $\bfG=\bfQ \bfW$. Note that if the $K$ KAs can be run in parallel, computing the precoding matrix $\bfG$ will require $O(MK)$ operations including the $O(MK)$ computational complexity of finding the row selection probability (which is the same for all KAs and is shared among them).  Also, in terms of computational complexity, it is beneficial to keep the precoding matrix $\bfG$ in the factorized form $\bfG=\bfQ \bfW$ and compute the precoding vector as follows $$\begin{aligned}
\bfy=\bfG \bfs = \bfQ \bfW \bfs =(\bfQ (\bfW \bfs)),\end{aligned}$$ where the parentheses in the last expression show the order of multiplication.  Otherwise, we will need an additional $O(M K^2)$ operations to compute $\bfG$ directly, which increases the complexity by a factor of $O(K)$. Now consider the UL scenario.  Since the detection matrix for the UL, $(\bfQ^\herm \bfQ + \xi \bfI_K)^{-1} \bfQ^\herm$ is the conjugate of the precoding matrix for the DL, we can obtain an approximation of the detection matrix as $\bfG^\herm=\bfW^\herm \bfQ^\herm$, where $\bfW$ denotes the $K \times K$ matrix produced at the output of the $K$ parallel KA in the DL scenario explained before[^4].  Analogous to the DL, it is beneficial in terms of complexity to keep $\bfG^\herm$ in the factorized from and compute the estimate of users symbols from the received signal $\bfy \in \bC^M$ in the UL by $$\begin{aligned}
\widehat{\bfs}=\bfG^\herm \bfy= \bfW^\herm \bfQ^\herm \bfy=(\bfW^\herm(\bfQ^\herm \bfy)),\end{aligned}$$ where again the parentheses in the last expression show the order of multiplication.  Recall that, assuming that $K$ KAs can be run in parallel, computing the detection matrix $\bfG^\herm$ (of course in its factorized form $\bfG^\herm=\bfW^\herm \bfQ^\herm$) requires $O(MK)$ operations. 

Achievable Ergodic Rate {#sec:erg rate}
=======================

In this section, we derive lower and upper bounds on the achievable ergodic rate of our proposed KA.  We first prove an important (although an easy-to-see) property of KA. 

\[lin\_prop\] Consider an arbitrary SLE $\bfA \bfx =\bfb$ where $\bfA$ is the $m \times n$ matrix of equations, where $\bfx \in \bC^n$ is the set of unknowns, and where $\bfb \in \bC^m$ is the set of known coefficients.  Let $\bfx^t$ be the estimate of KA at iteration $t$ starting from the zero initialization $\bfx^0=0$. Then, $\bfx^t=\scrG^t(\bfA) \bfb$ where $\scrG^t(\bfA) \in \bC^{n\times m}$ is a linear operator that depends on $\bfA$ and the internal randomization of KA until the iteration $t$ but not on $\bfb$. $\square$

We use induction on $t$. For $t=0$, we have $\bfx^0=0$ which is trivially a linear function (zero function) of $\bfb$. Let us assume that the induction hypothesis is true for the $t$-th iteration, i.e., $\bfx^t=\scrG^t(\bfA) \bfb$, and prove the claim for $t+1$. Suppose $r(t) \in [m]$ is the index of the row selected at iteration $t$ and let us denote the conjugate of this row by $\bfa_{r(t)}$. Then, we have $$\begin{aligned}
    \bfx^{t+1}&=\bfx^t+\frac{ b_{r(t)}- \inp{\bfa_{r(t)}}{ \bfx^t} }{\| \bfa_{r(t)}\|^2}\bfa_{r(t)}
    \stackrel{(i)}{=} (\bfI_n - \frac{ \bfa_{r(t)} \bfa_{r(t)}^\herm }{\| \bfa_{r(t)}\|^2} ) \bfx^t+\frac{ \bfa_{r(t)} \bfe_{r(t)}^\herm}{\| \bfa_{r(t)}\|^2}\bfb\\
    &\stackrel{(ii)}{=}\left (  (\bfI_n- \frac{ \bfa_{r(t)} \bfa_{r(t)}^\herm }{\| \bfa_{r(t)}\|^2} ) \scrG^t(\bfA) +  \frac{ \bfa_{r(t)} \bfe_{r(t)}^\herm}{\| \bfa_{r(t)}\|^2} \right ) \bfb
    =:\scrG^{t+1}(\bfA) \bfb,\label{q_proof_2}
    \end{aligned}$$ where in $(i)$ we defined $\bfe_{r(t)} \in \bC^m$ as the canonical vector that has $1$ at the index $r(t)$ and is zero elsewhere and replaced $b_{r(t)}=\bfe_{r(t)}^\herm \bfb$, and where in $(ii)$ we used the induction hypothesis that $\bfx^t=\scrG^t(\bfA) \bfb$ is a linear function of $\bfb$ with $\scrG^t(\bfA)$ not depending on $\bfb$. This implies that $\bfx^{t+1}=\scrG^{t+1}(\bfA) \bfb$, where it is seen that the linear operator $\scrG^{t+1}(\bfA) $ depends on $\bfA$ and on the internal randomization of KA until iteration $t+1$ but does not depend on $\bfb$. This proves the induction step and completes the proof. 

Using Proposition \[lin\_prop\], we will derive lower and upper bounds on the achievable ergodic rate when the KA is run for a specific number of iterations $t$, which can be generally much lower than $O(K)$ iterations typically needed for its convergence (see, e.g., Section \[sec:scenarios\]).  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 $\bfy=\bfH \bfs + \bfn$, where $\bfs \in \bC^K$ is the symbol sequence transmitted by the users, where $\bfH$ is the $M\times K$ matrix of the true channel state, and where $\bfn$ is the AWGN noise.  From Proposition \[lin\_prop\], it immediately results that, assuming a zero initialization, the signal estimate $\widehat{\bfs}_t$ produced by KA at a specific iteration $t$ will be a random linear function of $\bfy$. More precisely, we have $\widehat{\bfs}_t =\scrG^t(\bfQ) \bfy$, where the $K\times M$ matrix $\scrG^t(\bfQ)$ depends on the noisy channel state $\bfQ$ and the internal randomization of KA but does not depend on the received signal $\bfy$. Letting $\scrG^t(\bfQ)=[\bfg_1^t, \dots, \bfg_k^t]^\herm \in \bC^{K \times M}$ and considering the fact that the symbols estimate produced by KA is given by $\widehat{\bfs}_t =\scrG^t(\bfQ) \bfy$, we can interpret $\bfg_k^t \in \bC^M$, $k \in [K]$, as the detection vector corresponding to the $k$-th user at $t$-th iteration.  Note that the resulting detection matrix $\scrG^t(\bfQ)$ is generally far from MMSED matrix but converges to it when the KA is run for sufficiently many iterations $t$. Overall, one 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 $\bfT^{\text{up}}= \scrG^t(\bfQ) \bfH\in \bC^{K \times K}$, where for simplicity we dropped the explicit dependence of $\bfT^{\text{up}}$ on the iteration $t$.

We use the upper and lower bounds on the achievable ergodic rate developed in [@fundamental_MM].  In particular, by treating the interference as noise and by coding across several coherence blocks, an upper bound on the ergodic capacity is given by[^5] $$\begin{aligned}
\overline{R}_k=\bE\left [ \log\left(1+ \frac{|\bfT^{\text{up}}_{k,k}|^2}{\sum_{k'\neq k} |\bfT^{\text{up}}_{k,k'}|^2 +\sigma^2_{k}}\right) \right ],\label{upper_bound_dumm}\end{aligned}$$ where $\sigma^2_{k}=\frac{1}{\snr}\|\bfg_k^t\|^2$, with $\bfg_k^t$ denoting the conjugate of the $k$-th row of $\scrG^t(\bfQ)$, and where the expectation is taken over all the randomness of the channel state $\bfH$ and the detection matrix $\scrG^t(\bfQ)$ including the randomness due to KA.  This yields the upper bound $\overline{R}=\sum_{k \in [K]} \overline{R}_k$ on the achievable sum rate.  We will also use the following lower bound on the ergodic rate from [@fundamental_MM]: $$\begin{aligned}
\underline{R}_k=\log \left ( 1+ \frac{|\bE[\bfT^{\text{up}}_{k,k}]|^2}{\var(\bfT^{\text{up}}_{k,k})+ \sum_{k'\neq k} \bE[|\bfT^{\text{up}}_{k,k'}|^2]  +\bE[\sigma^2_{k}]}  \right),\label{lower_bound_dumm}\end{aligned}$$ where $\var(.)$ denotes the variance of a random variable.  This yields the lower bound $\underline{R}=\sum_{k \in [K]} \underline{R}_k$ on the ergodic rate.  Following a similar argument, we obtain the corresponding upper/lower bounds on the ergodic capacity for the DL. 

In the next section, we use Monte-Carlo (MC) simulations to empirically evaluate these lower and upper bounds on the performance of our algorithm.  Therefore, we need to be able to compute, at each instance of MC simulation, the underlying matrix $\scrG^t(\bfQ)$ for a specific iteration $t$. Let us consider a specific instance of MC simulation and let us denote by $\clR_t:=\{r(l): l=1,2, \dots, t\}$ the random indices of the rows selected by KA up to iteration $t$ while decoding users symbols in a UL scenario.  We can apply the iterative formula in the proof of Proposition \[lin\_prop\] to compute $\scrG^t(\bfQ)$ for each iteration $t$ explicitly.  As an alternative, we can apply the technique proposed in Section \[sec:matrix\_design\] for computing the detector/precoder matrix.  For example, in the UL scenario considered here, while running KA to decode users symbols from the noisy input in a specific instance of MC simulation, we also run in parallel $K$ KAs with corresponding inputs $\{\bfs_i=\bfe_i: i \in [K]\}$ as explained in Section \[sec:matrix\_design\].  However, we impose the additional restriction that all the $K$ KAs use the same row indices as in the original KA used for decoding users symbols.  In this way, we can compute, at each instance of MC simulation, $\scrG^t(\bfQ)$ along with the estimate of users symbols.  Repeating the same procedure for each MC simulation, we produce independent realizations of the random variables needed for computing upper/lower bounds in and , respectively.  A similar technique can be applied to a DL scenario. 

Simulation results {#sec:results}
==================

In this section, we evaluate the performance of our proposed technique empirically via numerical simulations. 

Channel Model
-------------

For all the simulations, we assume a correlated channel model for every user $k \in [K]$ where the true channel vector of the user $k$ is given by $\bfh^k= \Phim^{1/2}\bfz^k$ where $\bfz^k \sim \cg\left( 0,\bfI_M\right)$ and where $\Phim$ is the covariance matrix of $\bfh^k$ with entries given by $\left[ \bf \Phi\right] _{i,j}=a^{\left| i-j\right| }$ for some correlation parameter $a \in [0,1]$. We assume that the entries of $\bfz^k$, thus, those of $\bfh^k$, vary independently from one coherence period to another.  We also assume that the estimated channel matrix $\bfQ$, available at the BS in each coherence period, is a noisy version of the true channel matrix $\bfH$ given by $\bfQ=\sqrt{1-\tau^2}\bfH+\tau\bfN$ where $\bfN$ is a matrix with i.i.d.  entries $\left[ \bfN\right] _{m,k}\sim \cg\left( 0,1\right) $, and where $\tau \in [0,1]$ is the estimation quality parameter with $\tau=0$ yielding a perfect channel estimation. 

Kaczmarz Residual Effect {#sec:results:residual}
------------------------

In this section, we investigate via numerical simulations the effect of the residual error in a KA directly applied to an OD SLE.  We focus on a UL scenario where the noisy received signal at the BS is given by $\bfy=\bfH \bfs^\circ + \bfn^\circ$, where $\bfH$ is the true $M\times K$ channel matrix, where $\bfs^\circ \in \bC^K$ is the users symbols, and where $\bfn^\circ$ is the AWGN at BS antennas.  As explained in Sections \[sec:definition\] and \[sec:scenarios\], neglecting the effect of the noise $\bfn^\circ$ and applying the KA to the OD SLE $\bfH {\bfs}=\bfy$ with the unknown ${\bfs}$ and the noisy input $\bfy$ results in a residual error such that the decoded users symbols do not converge to the estimate produced by the MMSE detector.  This has been illustrated via numerical simulations in Fig.\[fig:compare\_error\_norm\], where it is seen that KA (Algorithm \[tab:Kac\]) directly applied to OD SLE $\bfH \bfs= \bfy$ yields a considerably large estimation error, whereas our proposed KA yields an estimation error that approaches that of the MMSE detector after sufficiently many iterations.  Note that for the simulations, we assume that the true channel state $\bfH$ is available at the BS to avoid the additional error caused by the noisy channel state and to pinpoint only the residual error due to KA.  In Fig.\[fig:compare\_error\_norm\], we also compare the performance of our algorithm with that of another method proposed in [@Herman], which applies the KA to the UD SLE $[\bfH, \sqrt{\xi} \bfI_K] \bfz=\bfy$ with the noisy input $\bfy$, where $\xi>0$ is the noise parameter and where $\bfz=[\bfs; \bfn] \in \bC^{M+K}$ contains both the symbol vector $\bfs\in \bC^K$ and the noise $\bfn\in \bC^M$. The underlying idea in [@Herman] is that the SLE $[\bfH, \sqrt{\xi} \bfI_K] \bfz=\bfy$ is always consistent, thus, the proposed KA avoids the effect of the residual (similarly to our proposed KA).  However, a direct calculation based on our convergence analysis of KA in Section \[sec:math frame\] reveals that the average gain of the resulting matrix $\bfB:=[\bfH, \sqrt{\xi} \bfI_K]$ in [@Herman] along the subspace $\clX$ produced by its rows (see, e.g., Definition \[avg\_gain\]) is given by $\kappa_\clX(\bfB)=\frac{\xi}{\|\bfH\|_\sfF^2 + K \xi}$, which is much lower than the gain $\frac{\lambda_{\min} (\bfH^\herm \bfH) + \xi}{\|\bfH\|_\sfF^2 + K \xi}$ corresponding to the matrix $[\bfH^\herm, \sqrt{\xi} \bfI_K]$ in our proposed KA (this is obtained from assuming the perfect channel state $\bfQ=\bfH$).  As a result, we expect that the KA in [@Herman] be extremely slower than our proposed KA, especially for high SNR ($\xi \to 0$). 

This is easily verified via numerical simulations, as illustrated in Fig.\[fig:compare\_error\_norm\].  It is seen that at a relatively low SNR ($\text{SNR}=0$dB), KA in [@Herman] ultimately convergess to the estimation error of the original MMSED but it converges much slower than our proposed KA.  At high SNR ($\text{SNR}=20$dB) where $\xi \to 0$, in contrast, we can identify two different convergence phases.  At first, KA in [@Herman] performs almost identically to the KA applied to the OD SLE $\bfH \bfs= \bfy$ (as expected due to $\xi \to 0$) and reduces the estimation error rapidly to the saturation limit of OD SLE $\bfH \bfs= \bfy$ caused by the noise residual effect.  Then, the algorithm spends the rest of the iterations on reducing/compensating the residual error until it approaches the estimation error of the original MMSED, but this happens very slowly due to the small value of $\xi$.

We also perform numerical simulations to investigate the effect of the residual error directly on the achievable ergodic rate in Massive MIMO.  Fig.\[fig:simple\_Kac\] illustrates the numerical results.  As expected, the residual error reduces the achievable rate considerably.  It is also seen that although the MMSED in [@Herman] generally performs better than the simple KA (Algorithm \[tab:Kac\]) at low SNRs, it completely loses its efficiency at high SNRs because of its extremely slow convergence.  Our proposed method, however, is able to converge to the performance of the original MMSED in all ranges of SNRs. 

Downlink {#sec:results:DL}
--------

In the DL, we compare the performance of our method with that of the *Truncated Polynomial Expansion* (TPE) technique proposed in [@Bjornson] in terms of spectral efficiency and total number of operations as a criterion for affordable computational complexity.  Fig.\[fig:sim:down1\] 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.  \[fig:sim:down1\] 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 [@Bjornson].  However, it has also the following crucial advantages compared with [@Bjornson]:

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.  \[fig:sim:down1\] 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.  \[fig:sim:down1\]. 

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.   When the channel response available at the BS changes quite fast with time (channel aging), or in OFDM-based systems where the channel response is available in only a subset of subcarriers (signal dimensions), the channel response needs to be, somehow, interpolated in time or in frequency respectively, to obtain the (approximate) channel response in other signal dimensions.  A scheme that computes directly the precoding vectors one-by-one, such as our proposed method, is much better suited to any channel tracking/prediction in time, interpolation in the OFDM subcarriers, etc.  than a scheme such as TPE that computes the precoding/detection matrix from a subset of signal dimensions and then keeps it fixed for the remaining signal dimensions inside a coherence block.  As a result, our proposed method can be combined with suitable channel prediction/interpolation techniques to obtain much better performances in situations where block fading channel model does not approximate the true channel very well. 

4.   The proposed algorithm does not require the knowledge of channel estimation quality parameter $\tau$.

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. 

Uplink {#sec:results:UL}
------

In the UL, we compare the performance of our algorithm with a very recent *approximate message passing* (AMP) technique proposed in [@ghods2017optimally; @Maleki_Studer_ISIT2015], especially with [@Maleki_Studer_ISIT2015] 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.  \[fig:sim:up1\] 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 ZFD and MMSED up to $20$dB SNR respectively.[^6] It is also seen that the AMP MMSED [@Maleki_Studer_ISIT2015] is able to converge to the original MMSED when the channel is uncorrelated ($a \approx 0$), with an order of $24MK$ operations.[^7] However, for a slightly correlated channel ($a\geq0.6$), neither the standard AMP [@Maleki_Studer_ISIT2015] nor the nonparametric AMP [@ghods2017optimally] are able to suitably recover the users 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. 

Gap-to-Capacity versus Computational Complexity
-----------------------------------------------

In this part, we compare our proposed method with the TPE [@Bjornson] and the AMP [@Maleki_Studer_ISIT2015; @ghods2017optimally] in terms of *normalized gap-to-capacity* versus the number of operations $t$ given by $\frac{S_{\max}-S_t}{S_{\max}}$, where $S_t$ denotes the upper/lower bound of per-user spectral efficiency of a specific beamforming/detection scheme after $t$ operations and where $S_{\max}$ denotes the maximum possible upper/lower bound of per-user spectral efficiency achieved by original linear schemes. 

This is illustrated in Fig.  \[fig:sim:gap\] 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 $\text{Gap}=10^{-2}$ with nearly $40MK$ operations.  It is also worthwhile to mention that the AMP algorithm [@Maleki_Studer_ISIT2015; @ghods2017optimally] is quite efficient for uncorrelated channels but fail in correlated channels. 

Conclusion {#sec:conclusion}
==========

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 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 very competetive with the schemes proposed in the literature based on random matrix theory and approximate message passing in terms of computations and has a comparable performance. 

[^1]: A short version of this paper was presented in the Workshop on Smart Antennas (WSA), Bochum, Germany, 2018 [@WSA_paper_published]. 

[^2]: The channel coherence block is the region in the time-frequency plane where the small-scale fading process generating the coefficients of $\bfH$ can be considered constant.  For a channel with physical coherence time $\Delta T_c$ and coherence bandwidth $\Delta W_c$, the channel coherence block spans approximately $\lceil \Delta T_c \times \Delta W_c \rceil$ complex signal dimensions. 

[^3]: In cases where these conditions are not met, i.e., when the system is dominated by inter-cell interference and/or when the CSI quality is poor due to high mobility, the advantage of ZFD/MMSED over MRC may disappear, as shown for example in [@fundamental_MM Section 6.3].  In those cases, it is obvious that MRC becomes a simple and attractive alternative to any more sophisticated detection scheme. 

[^4]: Note that, here, we decided to use Algorithm \[tab:Kac\_MMBF\] (initially designed for the DL) to compute beamforming matrices for both the UL and the DL.  Instead, we could apply Algorithm \[tab:Kac\_MMSED\] to $K$ parallel KAs with corresponding inputs $\{\bfb_i=\bfe_i: i \in [K]\}$ (see the definition of the variable $\bfb$ in Algorithm \[tab:Kac\_MMSED\]) to compute the beamforming matrix for the UL. 

[^5]: This upper bound is obtained in [@caire2017ergodic] under the condition that the true channel state $\bfH$ is available at the BS such that the channel interference matrix $\bfT^{\text{up}}$ is fully known to the BS.  This yields an upper bound on the achievable ergodic rate in our scenario, where only a noisy version of $\bfH$ is assumed to be known at the BS, thus, the channel interference matrix $\bfT^{\text{up}}= \scrG^t(\bfQ)^\herm \bfH$ is not perfectly known at the BS (since $\bfH$ is not known). 

[^6]: for the upper bounds, $32MK$ iterations suffice. 

[^7]: 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 $4MK$ additions and multiplications per iteration.