#  Reflective block Kaczmarz algorithms for least squares

Changpeng Shao<EMAIL_ADDRESS>Academy of Mathematics and Systems
Science
Chinese Academy of Sciences
Beijing, 100190 China

###### Abstract

In [Steinerberger, Q. Appl. Math., 79:3, 419-429, 2021] and [Shao, SIAM J.
Matrix Anal. Appl. 44(1), 212-239, 2023], two new types of Kaczmarz
algorithms, which share some similarities, for consistent linear systems were
proposed. These two algorithms not only compete with many previous Kaczmarz
algorithms but, more importantly, reveal some interesting new geometric
properties of solutions to linear systems that are not obvious from the
standard viewpoint of the Kaczmarz algorithm. In this paper, we
comprehensively study these two algorithms. First, we theoretically analyse
the algorithms given in [Steinerberger, Q. Appl. Math., 79:3, 419-429, 2021]
for solving least squares. Second, we extend the two algorithms to block
versions and provide their theoretical convergence rates. Our numerical
experiments also verify the efficiency of these algorithms. Third, as a
theoretical complement, we address some key questions left unanswered in
[Shao, SIAM J. Matrix Anal. Appl. 44(1), 212-239, 2023].

Key words. least squares; block Kaczmarz algorithm; block Householder matrix;
reflection.

##  1 Introduction

The Kaczmarz method is a useful iterative algorithm for solving linear systems
of equations [11]. Due to its simple iterative scheme, various variants of the
Kaczmarz method have been studied in the past, e.g., [7, 2, 4, 14, 16, 17, 25,
27, 15, 9, 18, 10, 22, 21, 19, 13, 12]. In this paper, we focus on two new
types of Kaczmarz algorithms proposed recently in [22, 21].

Let $A$ be an $m\times n$ real matrix and ${\bf b}$ be an $m\times 1$ real
vector. We use $A_{i}$ to denote the $i$-th row of $A$. The Kaczmarz algorithm
for solving the linear system $A{\bf x}={\bf b}$, as studied in [22, 21],
reads as follows: arbitrarily choose ${\bf x}_{0}\in\mathbb{R}^{n}$, update

${\bf x}_{k+1}={\bf x}_{k}+2\frac{b_{i_{k}}-A_{i_{k}}^{{\tt T}}{\bf
x}_{k}}{\|A_{i_{k}}\|^{2}}A_{i_{k}},\quad k=0,1,2,\ldots,$ (1.1)

where in [22], $i_{k}\in[m]$ is chosen with probability
$\|A_{i_{k}}\|^{2}/\|A\|_{\tt F}^{2}$, while in [21], $i_{k}=(k\mod m)+1$.
Here $\|A\|_{\tt F}$ is the Frobenius norm of $A$. So, the former is
randomised, and the latter is deterministic. As shown in [22, 21], these two
new types of Kaczmarz algorithms can outperform previous ones, such as the
randomised Kaczmarz algorithm [25] and the block Kaczmarz algorithm [17].
Additionally, these two Kaczmarz algorithms, especially the one presented in
[21], demonstrate some new and interesting geometric properties that are not
obvious from the standard Kaczmarz algorithm viewpoint. These properties
connect solutions of linear systems to centers of certain spheres. For
example, a main result in [21] states that the iterative process (1.1)
produces a series of points on $2m$ high-dimensional spheres, which can
essentially be decomposed into $m$ pairs. For each pair, the midpoint of the
centers of the spheres is exactly the solution of $A{\bf x}={\bf b}$ when it
is consistent. Thus, the problem of solving a consistent linear system of
equations is closely relevant to the problem of finding the center of a sphere
from knowing many points on the sphere, and the latter problem can be solved
various ways.

The two papers [22, 21] provided a solid starting point for these algorithms;
however, many questions remain unanswered. For example, the analysis of the
algorithm given in [22] is mainly for invertible linear systems. However,
inconsistent linear systems are more common in practice. In [21], the author
thoroughly studied the iterative procedure (1.1) for both consistent and
inconsistent linear systems. The big difference between [21] and [22] is the
convergence rate. The convergence rate of the algorithm studied in [21]
depends on a new quantity $\eta(A)$ that is closely related to the condition
number $\kappa(A)$ of $A$, while the convergence rate of [22] (including many
previous ones) depends on $\kappa_{\tt F}(A):=\|A\|_{\tt F}\|A^{-1}\|$. As
$\kappa(A)\leq\kappa_{\tt F}(A)$, the algorithm given in [21] could converges
faster. However, in [21], few properties of $\eta(A)$ were found.

In this paper, we comprehensively study these two types of new Kaczmarz
algorithms. The contributions of this paper are threefold: First, we
theoretically analyse the algorithms given in [22] for solving least squares,
rather than just invertible linear systems. Second, we extend the algorithms
from [22, 21] to block versions and provide their theoretical convergence
rates. Third, as a theoretical complement of our previous work, we address
some key questions mentioned in [21], particularly regarding the properties of
$\eta(A)$.

We describe these results in more detail below. To avoid confusion with the
standard Kaczmarz method and some previous ones (where the constant 2 is
replaced by 1 in (1.1), making each step an orthogonal projection), we refer
to the iterative process in (1.1) as the reflective Kaczmarz algorithm.

###  1.1 Randomised reflective Kaczmarz algorithm for least squares

In [22], the author proved that if $A$ is invertible and ${\bf
x}_{*}=A^{-1}{\bf b}$, then

${\mathbb{E}}\left[\left\|{\bf x}_{*}-\frac{1}{N}\sum_{k=0}^{N-1}{\bf
x}_{k}\right\|^{2}\right]\leq\frac{1+\|A\|_{\tt F}^{2}\|A^{-1}\|^{2}}{N}\|{\bf
x}_{*}-{\bf x}_{0}\|^{2}.$ (1.2)

In this paper, we show a similar result for general $A$. Now let ${\bf
x}_{*}=A^{+}{\bf b}$ be a solution of the least square problem
$\arg\min\|A{\bf x}-{\bf b}\|$ and assume that ${\rm Rank}(A)=n$, then (See
Theorem 3.3)

${\mathbb{E}}\left[\left\|{\bf x}_{*}-\frac{1}{N}\sum_{k=0}^{N-1}{\bf
x}_{k}\right\|^{2}\right]\leq\frac{1+\|A\|_{\tt F}^{2}\|A^{-1}\|^{2}}{N}\|{\bf
x}_{*}-{\bf x}_{0}\|^{2}+O(\|{\bf c}\|^{2}/\sigma_{n}^{2}),$ (1.3)

where ${\bf c}={\bf b}-A{\bf x}_{*}$ is orthogonal to the column space of $A$
and $\sigma_{n}$ is the minimal singular value. The convergence rates (1.2)
and (1.3) imply that, to ensure a good approximation, $N$ should be as large
as $\|A\|_{\tt F}^{2}\|A^{-1}\|^{2}$. This coincides with many previous
Kaczmarz algorithms, such as [25]. The result (1.3) is similar to the one
obtained in [27, Theorem 3.7] for standard randomised Kaczmarz algorithm for
least squares.

###  1.2 Block reflective Kaczmarz algorithm

In block Kaczmarz algorithm, a subset of rows of $A$ is used at each step of
iteration. Block Kaczmarz algorithms can converge faster than the standard
ones, as has been verified in many previous works [7, 14, 17, 15, 19]. To
establish a block version of (1.1), we consider block Householder matrix,
defined as follows: Let $M$ be a $r\times n$ matrix with full row rank, the
block Householder matrix in terms of $M$ is defined by
${\mathcal{H}}_{M}:=I_{n}-2M^{{\tt T}}(MM^{{\tt T}})^{-1}M.$ Let
$Z\subseteq[m]$, we set $A_{Z}=\sum_{i\in Z}{\bf e}_{i}A_{i}^{\tt
T}\in\mathbb{R}^{|Z|\times n}$ as the submatrix of $A$ whose rows are indexed
by $Z$. The corresponding block Householder matrix is ${\mathcal{H}}_{A_{Z}}$.
The iterative procedure now becomes

${\bf x}_{k+1}={\mathcal{H}}_{A_{Z}}{\bf x}_{k}+2A_{Z}^{+}{\bf b}_{Z}={\bf
x}_{k}-2A_{Z}^{+}(A_{Z}{\bf x}_{k}-{\bf b}_{Z}),\quad{\bf b}_{Z}=\sum_{i\in
Z}b_{i}{\bf e}_{i}\in\mathbb{R}^{|Z|}.$ (1.4)

In the above, the second formula is a natural generalisation of (1.1).
Regarding the first one, our goal is to have the iterative procedure satisfy
the following condition: ${\bf x}_{k+1}-{\bf x}_{*}={\mathcal{H}}_{A_{Z}}({\bf
x}_{k}-{\bf x}_{*})$, where ${\bf x}_{*}$ is a solution of $A{\bf x}={\bf b}$.
Simplifying this leads to the first formula.

Regarding (1.4), there are two types of methods for choosing $Z$, each leading
to a different type of block Kaczmarz algorithm.

  * •

Algorithm 1: Let $\\{1,2,\ldots,m\\}=Z_{1}\cup\cdots\cup Z_{p}$ be a fixed
partition of the row indices. Arbitrarily choose ${\bf x}_{0}$, then update

${\bf x}_{k+1}={\bf x}_{k}-2A_{Z_{i}}^{+}A_{Z_{i}}{\bf
x}_{k}+2A_{Z_{i}}^{+}{\bf b}_{Z_{i}},\quad k=0,1,2,\ldots,$ (1.5)

where $i\in[p]$ is chosen with probability $\|A_{Z_{i}}\|_{\tt
F}^{2}/\|A\|_{\tt F}^{2}$.

  * •

Algorithm 2: Fix an integer $q\leq\min(m,n)$. Arbitrarily choose ${\bf
x}_{0}$, then update

${\bf x}_{k+1}={\bf x}_{k}-2A_{Z}^{+}A_{Z}{\bf x}_{k}+2A_{Z}^{+}{\bf
b}_{Z},\quad k=0,1,2,\ldots,$ (1.6)

where $Z=\\{i_{1},\ldots,i_{q}\\}\subseteq[m]$ is obtained in a way such that
$i_{j}$ is collected with probability $\|A_{i_{j}}\|^{2}/\|A\|_{\tt F}^{2}$.

In Section 4, we analyse the convergence rates of the above two algorithms for
both consistent and inconsistent linear systems, see Theorems 4.4, 4.6, 4.9
and 4.11 for more details. The results indicate that these two algorithms have
similar convergence rates. We also numerically test their performance in
Section 6. The numerical results affirm that the efficiency of randomised
block reflective Kaczmarz algorithm than the standard one.

###  1.3 More theoretical results on deterministic reflective Kaczmarz
algorithm

In the third part of this paper (see Section 5), our main focus in on the
deterministic version of the iterative process (1.1), which was previously
studied in [21]. Here, we explore more theoretical results about this
algorithm.

  * •

First, we extend this algorithm to a block version and analyse the
corresponding convergence rate. Our main results (i.e., Theorems 5.7 and 5.9)
indicate that all key results proved in [21] also hold in the block version.

  * •

Second, we analyse a quantity, denoted as $\eta(A)$, see Definition 5.6, that
is closely related to the convergence rate. In [21], a upper bound of
$\eta(A)$ is given, i.e., $\eta(A)=O(\kappa(A)^{2}\log m)$. Here we provide
more properties about this quantity. For example, in Theorem 5.19, we will
prove that

$\frac{|\eta^{2}-\sigma^{-2}_{\min}(A)|^{2}}{\eta^{2}\sigma^{-2}_{\min}(A)\Big{(}1+\eta^{2}\Big{)}\Big{(}1+\sigma^{-2}_{\min}(A)\Big{)}}\leq
4\|L\|^{2},$ (1.7)

where $L$ is the lower triangular matrix such that $L+L^{\tt T}=AA^{\tt T}-I$
and $\sigma_{\min}(A)$ is the minimal nonzero singular value of $A$.
Particularly, if $A$ has restricted isometry property with restricted isometry
constant $\delta$, then the above upper bound if $O(\delta^{2}\log^{2}m)$. We
also study the behavior of $\eta(A)$ under left or right orthogonal actions,
see Proposition 5.16 and Theorem 5.23. We will show that $\eta(A)$ is
invariant under right orthogonal actions, but not under left orthogonal
actions. In the latter case, we give a quantitative description of the change.
The result is similar to (1.7). This property is in sharp contrast to singular
values, which are invariant under both left and right orthogonal actions. We
also provide an equation to compute $\eta(A)$, which is similar to the
characteristic equation for singular values, see Theorem 5.18. Interestingly,
this equation was indeed previously found by Coxeter [5] in 1939.

  * •

Third, we provide further results on connections to solving least square
problems. As analysed in [21], the iterative process (1.1) indeed solves
$A^{\tt T}L^{-1}A{\bf x}=A^{\tt T}L^{-1}{\bf b}$ rather than $A^{\tt T}A{\bf
x}=A^{\tt T}{\bf b}$. Some reasons were provided. In this paper, we give more
results on this and present a method to fix this, see Corollary 5.27 and
Proposition 5.28.

##  2 Preliminaries

###  2.1 Notation

For a matrix $A=(a_{ij})$, its Frobenius norm is defined as $\|A\|_{\tt
F}=(\sum_{i,j}|a_{ij}|^{2})^{1/2}$, its operator norm is the maximum singular
value, denoted as $\|A\|$ or $\sigma_{\max}(A)$. We use $\sigma_{\min}(A)$ to
denote its minimum nonzero singular value, and $A^{+}$ to denote its
pseudoinverse. The condition number is $\kappa=\|A\|\|A^{+}\|$, and the scaled
condition number is $\kappa_{\tt F}=\|A\|_{\tt F}\|A^{+}\|$. For a matrix $A$,
we use $A^{\dagger},A^{\tt T}$ to denote its complex conjugate and transpose
respectively. For any two square matrices $A,B$ of the same dimension, with
$A\preceq B$, we mean for any ${\bf v}$ we have ${\bf v}^{\tt T}A{\bf
v}\leq{\bf v}^{\tt T}B{\bf v}$. We use $\\{{\bf e}_{1},\ldots,{\bf e}_{n}\\}$
to represent the standard basis of $\mathbb{R}^{n}$. For any two vectors ${\bf
a},{\bf b}$, we use $\langle{\bf a},{\bf b}\rangle$ to represent their inner
product. Let ${\bf v}$ be a vector; sometimes we will use Dirac notation
$|{\bf v}\rangle$ to represent the unit vector ${\bf v}/\|{\bf v}\|$ when it
makes the notation more concise and neat. Let $m$ be an integer, we use $[m]$
to denote the set $\\{1,2,\ldots,m\\}$.

###  2.2 Perturbation of generalised eigenvalues

Let $A,B\in\mathbb{C}^{m\times m}$ be two matrices, in the generalised
eigenvalue problem [24], we aim to find all $(\alpha,\beta)\neq(0,0)$ with
$\alpha,\beta\in\mathbb{C}$ and ${\bf x}\in\mathbb{C}^{m}$ such that

$\beta A{\bf x}=\alpha B{\bf x}.$

The matrix pair $(A,B)$ is called regular if $\det(\beta A-\alpha B)$ is not
identically zero. The distance between the eigenvalues
$(\alpha_{1},\beta_{1})$ and $(\alpha_{2},\beta_{2})$ is measured by the
chordal metric

$\rho((\alpha_{1},\beta_{1}),(\alpha_{2},\beta_{2}))=\frac{|\alpha_{1}\beta_{2}-\alpha_{2}\beta_{1}|}{\sqrt{|\alpha_{1}|^{2}+|\beta_{1}|^{2}}\sqrt{|\alpha_{2}|^{2}+|\beta_{2}|^{2}}}.$

######  Proposition 2.1 (see [23]).

Let $(A,B)$ and $(\widetilde{A},\widetilde{B})=(A+E,B+F)$ be regular matrix
pairs, and let $(\alpha,\beta)$ be an eigenvalue of $(A,B)$, then there is an
eigenvalue $(\tilde{\alpha},\tilde{\beta})$ of $(\widetilde{A},\widetilde{B})$
satisfying

$\rho((\alpha,\beta),(\tilde{\alpha},\tilde{\beta}))=\frac{\sqrt{\|A\|^{2}+\|B\|^{2}}}{\gamma(A,B)}\left(\frac{\sqrt{\|E\|^{2}+\|F\|^{2}}}{\sqrt{\|A\|^{2}+\|B\|^{2}}}\right)^{1/m},$

where

$\gamma(A,B)=\max_{|\alpha|^{2}+|\beta|^{2}=1}\sigma_{\min}(\beta A-\alpha
B).$

######  Remark 2.2.

In the above, if $\beta=0,\alpha=-1$, then $\sigma_{\min}(\beta A+\alpha
B)=\sigma_{\min}(B)$. So $\gamma(A,B)\geq\sigma_{\min}(B)$. Similarly,
$\gamma(A,B)\geq\sigma_{\min}(A)$.

Let $Z=(A,B),W=(C,D)$ be regular matrix pairs with eigenvalues
$(\alpha_{i},\beta_{i})$ and $(\gamma_{i},\delta_{i})$, respectively. Define
the generalized spectral variation of $W$ with respect to $Z$ by

$S_{Z}(W):=\max_{i}\min_{j}\rho((\alpha_{j},\beta_{j}),(\gamma_{i},\delta_{i})).$

View $Z=(A,B)$ as a matrix, set $P_{Z}=Z^{+}Z$. We define a metric as follows:

$d_{2}(Z,W):=\|P_{Z}-P_{W}\|.$

By [8, Equation (1.13)],

$d_{2}(Z,W)\leq\|(ZZ^{\dagger})^{-1/2}\|\|Z-W\|.$ (2.1)

######  Proposition 2.3 (Theorem 2.1 in [8]).

Let $Z=(A,B)$ be a diagonalisable pair, i.e., there exists matrices $S,T$ such
that both $SAT$ and $SBT$ are diagonals. Let $W$ be a regular pair, then

$S_{Z}(W)\leq\|T\|\|T^{-1}\|d_{2}(Z,W).$

##  3 Randomised reflective Kaczmarz algorithm for least squares

Let $A\in\mathbb{R}^{m\times n}$ and ${\bf b}\in\mathbb{R}^{m}$, for least
squares, we aim to find $\tilde{{\bf x}}_{*}$ such that $\|\tilde{{\bf
x}}_{*}-{\bf x}_{*}\|\leq\varepsilon\|{\bf x}^{*}\|$, where ${\bf
x}_{*}:=\arg\min\|A{\bf x}-{\bf b}\|$. In this section, we study the following
randomised iterative scheme, which we call the reflective Kaczmarz algorithm,
for solving least squares: Arbitrarily choose ${\bf x}_{0}$, then update via

${\bf x}_{k+1}={\bf x}_{k}+2\frac{b_{i_{k}}-A_{i_{k}}^{{\tt T}}{\bf
x}_{k}}{\|A_{i_{k}}\|^{2}}A_{i_{k}},\quad k=0,1,2,\ldots,$ (3.1)

where $i_{k}\in[m]$ is chosen with probability $\|A_{i_{k}}\|^{2}/\|A\|_{\tt
F}^{2}$. Our main goal of this section is to prove the convergence rate of the
above iterative for least squares, see Theorem 3.3.

For convenience, we set

$\mathcal{R}_{i}=I_{n}-2\frac{A_{i}^{\tt T}A_{i}}{\|A_{i}\|^{2}}$

as the reflection generated by the $i$-th row of $A$, then

${\bf x}_{k+1}=\mathcal{R}_{i_{k}}{\bf
x}_{k}+2\frac{b_{i_{k}}}{\|A_{i_{k}}\|^{2}}A_{i_{k}}.$

Note that when $A$ is invertible, we have ${\bf x}_{k+1}-{\bf
x}_{*}=\mathcal{R}_{i_{k}}({\bf x}_{k}-{\bf x}_{*})$, and so the vectors ${\bf
x}_{0},{\bf x}_{1},\ldots$ lie on the sphere centered at the solution ${\bf
x}_{*}$ with radius $\|{\bf x}_{0}-{\bf x}_{*}\|$. One simple way to
approximate the centre is using the average of $\\{{\bf x}_{0},{\bf
x}_{1},\ldots,{\bf x}_{N-1}\\}$ for some $N$. In [22], Steinerberger computed
that it suffices to choose $N=O(\|A\|_{\tt
F}^{2}\|A^{-1}\|^{2}\log(1/\varepsilon))$ when $A$ is invertible. Below, we
study the theoretical results of this algorithm for solving least squares.

We decompose ${\bf b}=A{\bf x}_{*}+{\bf c}$, where ${\bf x}_{*}=A^{+}{\bf b}$
is an optimal solution of $\arg\min\|A{\bf x}-{\bf b}\|$, and ${\bf c}$ is
orthogonal to the column space of $A$ (i.e., $A^{\tt T}{\bf c}=0$). We denote
the error in the $k$-th step of iteration as

${\bf e}_{k}={\bf x}_{k}-{\bf x}_{*},$ (3.2)

then

${\bf e}_{k+1}={\bf e}_{k}+2\frac{A_{i_{k}}^{{\tt T}}{\bf
x}_{*}+c_{i_{k}}-A_{i_{k}}^{{\tt T}}{\bf
x}_{k}}{\|A_{i_{k}}\|^{2}}A_{i_{k}}=\mathcal{R}_{i_{k}}{\bf
e}_{k}+\frac{2c_{i_{k}}}{\|A_{i_{k}}\|^{2}}A_{i_{k}}.$ (3.3)

It is easy to check that

${\mathbb{E}}_{i}\left[\frac{c_{i}}{\|A_{i}\|^{2}}A_{i}\right]=\sum_{i=1}^{m}\frac{c_{i}}{\|A\|_{\tt
F}^{2}}A_{i}=\frac{{\bf c}^{\tt T}A}{\|A\|_{\tt F}^{2}}=0.$ (3.4)

Here we used the fact that $A^{\tt T}{\bf c}=0$.

######  Lemma 3.1.

The expectation value of $\|{\bf e}_{k}\|^{2}$ satisfies:

${\mathbb{E}}\left[\|{\bf e}_{k}\|^{2}\right]=\|{\bf
e}_{0}\|^{2}+\frac{4k\|{\bf c}\|^{2}}{\|A\|_{\tt F}^{2}}.$ (3.5)

###### Proof.

We first compute the expectation value with respect to $i_{k}$ using (3.3).

$\displaystyle{\mathbb{E}}_{i_{k}}\left[\|{\bf e}_{k+1}\|^{2}\right]$
$\displaystyle=$
$\displaystyle{\mathbb{E}}_{i_{k}}\left[\|\mathcal{R}_{i_{k}}{\bf
e}_{k}\|^{2}\right]+{\mathbb{E}}_{i_{k}}\left[\frac{4c_{i_{k}}^{2}}{\|A_{i_{k}}\|^{2}}\right]+{\mathbb{E}}_{i_{k}}\left[\frac{4c_{i_{k}}}{\|A_{i_{k}}\|^{2}}\langle
A_{i_{k}},\mathcal{R}_{i_{k}}{\bf e}_{k}\rangle\right]$ $\displaystyle=$
$\displaystyle\|{\bf e}_{k}\|^{2}+\frac{4\|{\bf c}\|^{2}}{\|A\|_{\tt F}^{2}}.$

Here in the last step we used (3.4) and the fact that $\langle
A_{i_{k}},\mathcal{R}_{i_{k}}{\bf e}_{k}\rangle=-\langle A_{i_{k}},{\bf
e}_{k}\rangle$. Inductively, we obtain the claimed result. ∎

######  Lemma 3.2.

Assume that $A$ has full column rank, then for any $\ell\geq 0$

$\Big{|}{\mathbb{E}}[\langle{\bf e}_{k},{\bf
e}_{k+\ell}\rangle]\Big{|}\leq\Big{(}1-2\kappa_{\tt
F}^{-2}\Big{)}^{\ell}\Big{(}\|{\bf e}_{0}\|^{2}+\frac{4k\|{\bf
c}\|^{2}}{\|A\|_{\tt F}^{2}}\Big{)},$

where $\kappa_{\tt F}=\|A\|_{\tt F}\|A^{+}\|$.

###### Proof.

We first prove the case when $\ell=1$. By (3.3) and (3.4),

$\displaystyle{\mathbb{E}}_{i_{k}}[\langle{\bf e}_{k},{\bf e}_{k+1}\rangle]$
$\displaystyle=$ $\displaystyle{\mathbb{E}}_{i_{k}}[\langle{\bf
e}_{k},\mathcal{R}_{i_{k}}{\bf e}_{k}\rangle]+\langle{\bf
e}_{k},{\mathbb{E}}_{i_{k}}\left[\frac{2c_{i_{k}}}{\|A_{i_{k}}\|^{2}}A_{i_{k}}\right]\rangle$
$\displaystyle=$ $\displaystyle\langle{\bf
e}_{k},\sum_{i=1}^{m}\frac{\|A_{i}\|^{2}}{\|A\|_{\tt
F}^{2}}(I-2|A_{i}\rangle\langle A_{i}|){\bf e}_{k}\rangle$ $\displaystyle=$
$\displaystyle\langle{\bf e}_{k},(I-\frac{2A^{\tt T}A}{\|A\|_{\tt F}^{2}}){\bf
e}_{k}\rangle.$

Since $A$ has full column rank, we have111If $A$ does not have full column
rank, then we can only obtain $|{\mathbb{E}}_{i_{k}}[\langle{\bf e}_{k},{\bf
e}_{k+1}\rangle]|\leq\|{\bf e}_{k}\|^{2}$ here. This is not enough for
convergence.

$|{\mathbb{E}}_{i_{k}}[\langle{\bf e}_{k},{\bf
e}_{k+1}\rangle]|\leq(1-\frac{2}{\|A^{+}\|^{2}\|A\|_{\tt F}^{2}})\|{\bf
e}_{k}\|^{2}=(1-2\kappa_{\tt F}^{-2})\|{\bf e}_{k}\|^{2}.$

By Lemma 3.1, this further means

$|{\mathbb{E}}[\langle{\bf e}_{k},{\bf e}_{k+1}\rangle]|\leq(1-2\kappa_{\tt
F}^{-2}){\mathbb{E}}[\|{\bf e}_{k}\|^{2}]=(1-2\kappa_{\tt F}^{-2})(\|{\bf
e}_{0}\|^{2}+\frac{4k\|{\bf c}\|^{2}}{\|A\|_{{\tt F}}^{2}}).$

More generally, for any $\ell\geq 1$

$\displaystyle{\mathbb{E}}[\langle{\bf e}_{k},{\bf e}_{k+\ell}\rangle]$
$\displaystyle=$
$\displaystyle{\mathbb{E}}\Big{[}{\mathbb{E}}_{i_{k}}[\langle{\bf
e}_{k},\mathcal{R}_{i_{k}}{\bf e}_{k+\ell-1}\rangle]\Big{]}$ $\displaystyle=$
$\displaystyle{\mathbb{E}}[\langle{\bf
e}_{k},\sum_{i=1}^{m}\frac{\|A_{i}\|^{2}}{\|A\|_{\tt
F}^{2}}(I-2|A_{i}\rangle\langle A_{i}|){\bf e}_{k+\ell-1}\rangle]$
$\displaystyle=$ $\displaystyle{\mathbb{E}}[\langle{\bf
e}_{k},(I-\frac{2A^{\tt T}A}{\|A\|_{\tt F}^{2}}){\bf e}_{k+\ell-1}\rangle].$

By induction,

$\displaystyle{\mathbb{E}}[\langle{\bf e}_{k},{\bf
e}_{k+\ell}\rangle]={\mathbb{E}}[\langle{\bf e}_{k},(I-\frac{2A^{\tt
T}A}{\|A\|_{\tt F}^{2}})^{\ell}{\bf e}_{k}\rangle].$

Thus,

$|{\mathbb{E}}[\langle{\bf e}_{k},{\bf e}_{k+\ell}\rangle]|\leq(1-2\kappa_{\tt
F}^{-2})^{\ell}(\|{\bf e}_{0}\|^{2}+\frac{4k\|{\bf c}\|^{2}}{\|A\|_{\tt
F}^{2}}).$

This completes the proof. ∎

######  Theorem 3.3 (Convergence rate of reflective Kaczmarz algorithm for
least squares).

Assume that $A\in\mathbb{R}^{m\times n}$ has full column rank and $\sigma_{n}$
is the minimal nonzero singular value. Let ${\bf b}=A{\bf x}_{*}+{\bf c}$ be
an orthogonal decomposition, where ${\bf x}_{*}=A^{+}{\bf b}$. Let $\\{{\bf
x}_{0},{\bf x}_{1},{\bf x}_{2},\ldots\\}$ be the points generated by the
process (3.1). Then

$\displaystyle{\mathbb{E}}\left[\,\left\|{\bf
x}_{*}-\frac{1}{N}\sum_{k=0}^{N-1}{\bf x}_{k}\right\|^{2}\,\right]$
$\displaystyle\leq$ $\displaystyle\frac{1+\kappa_{\tt F}^{2}}{N}\|{\bf
x}_{*}-{\bf x}_{0}\|^{2}+\left(2+2\kappa_{\tt F}^{2}+10\frac{\kappa_{\tt
F}^{6}}{N^{2}}\right)\frac{\|{\bf c}\|^{2}}{\|A\|_{\tt F}^{2}}$ (3.6)
$\displaystyle=$ $\displaystyle\frac{1+\kappa_{\tt F}^{2}}{N}\|{\bf
x}_{*}-{\bf x}_{0}\|^{2}+O\left(\frac{\|{\bf
c}\|^{2}}{\sigma_{n}^{2}}\right).$

###### Proof.

Using the notation (3.2), it suffices to estimate

$\displaystyle{\mathbb{E}}\left[\,\left\|\sum_{k=0}^{N-1}{\bf
e}_{k}\right\|^{2}\,\right]$ $\displaystyle=$
$\displaystyle{\mathbb{E}}\left[\sum_{k=0}^{N-1}\|{\bf
e}_{k}\|^{2}+\sum_{k\neq\ell}\langle{\bf e}_{k},{\bf e}_{\ell}\rangle\right]$
$\displaystyle=$ $\displaystyle N\|{\bf e}_{0}\|^{2}+\frac{2N(N-1)\|{\bf
c}\|^{2}}{\|A\|_{\tt
F}^{2}}+2\sum_{k=0}^{N-2}\sum_{\ell=k+1}^{N-1}{\mathbb{E}}[\langle{\bf
e}_{k},{\bf e}_{\ell}\rangle]$ $\displaystyle\leq$ $\displaystyle N\|{\bf
e}_{0}\|^{2}+\frac{2N(N-1)\|{\bf c}\|^{2}}{\|A\|_{{\tt
F}}^{2}}+2\sum_{k=0}^{N-2}\sum_{\ell=k+1}^{N-1}(1-2\kappa_{\tt
F}^{-2})^{\ell-k}\left(\|{\bf e}_{0}\|^{2}+\frac{4k\|{\bf c}\|^{2}}{\|A\|_{\tt
F}^{2}}\right).$

The second step is based on Lemma 3.1 and the third step is caused by Lemma
3.2. It is easy to compute that for any $0\leq q\leq 1$,

$\displaystyle\sum_{k=0}^{N-2}\sum_{l=k+1}^{N-1}q^{l-k}$ $\displaystyle=$
$\displaystyle\frac{q^{N+1}-q+Nq(1-q)}{(1-q)^{2}}~{}~{}\leq~{}~{}\frac{N}{1-q},$
$\displaystyle\sum_{k=0}^{N-2}\sum_{l=k+1}^{N-1}kq^{l-k}$ $\displaystyle=$
$\displaystyle\frac{qN^{2}(1-q)^{2}+q(3-q)+2q(1-q^{N})}{2(1-q)^{3}}~{}~{}\leq~{}~{}\frac{N^{2}}{2(1-q)}+\frac{5}{2(1-q)^{3}}.$

Thus, replacing $q=1-2\kappa_{\tt F}^{-2}$ gives

$\displaystyle{\mathbb{E}}\left[\,\left\|\sum_{k=0}^{N-1}{\bf
e}_{k}\right\|^{2}\,\right]\leq N(1+\kappa_{\tt F}^{2})\|{\bf
e}_{0}\|^{2}+(2N(N-1)+2N^{2}\kappa_{\tt F}^{2}+10\kappa_{\tt
F}^{6})\frac{\|{\bf c}\|^{2}}{\|A\|_{\tt F}^{2}}.$

Dividing both sides by $N$ yields the claimed bound. ∎

The result above coincides with the main theorem of [22] for invertible linear
systems. The error bound obtained in Theorem 3.3 is also very close to that of
the standard Kaczmarz method for inconsistent linear systems, e.g., see [27,
Theorem 3.7], [16, Theorem 2.1]. All are influenced by a convergence horizon
$\|{\bf c}\|^{2}/\sigma_{n}^{2}$. From Theorem 3.3, to ensure that the error
is bounded by $\varepsilon^{2}$, we need to choose $N=O(\kappa_{\tt
F}^{2}/\varepsilon^{2})$. As analyzed in [22, 21], the dependence on
$\varepsilon$ can be reduced to $\log(1/\varepsilon)$ using an idea similar to
the binary method. Essentially, we run the algorithm for a while until the
error is smaller than half of $\|{\bf x}_{*}-{\bf x}_{0}\|^{2}$, compute the
average of ${\bf x}_{0},\ldots,{\bf x}_{N-1}$, and then start the algorithm
again with this average as a new initial point.

##  4 Block reflective Kaczmarz algorithms

Block Kaczmarz methods have been widely studied in the past, e.g., see [7, 14,
17, 15, 19]. Their performance is generally better than that of the standard
Kaczmarz method. We below extend the algorithm (3.1) to block versions. Below,
we will theoretically prove that block reflective Kaczmarz algorithms can be
more efficient than the standard reflective Kaczmarz algorithm. Additionally,
in Section 6, the numerical tests will also verify this.

First, we recall the concept of block Householder matrix.

######  Definition 4.1 (Block Householder matrix).

Let $M$ be a $r\times n$ matrix with full row rank, the block Householder
matrix in terms of $M$ is defined by

${\mathcal{H}}_{M}:=I_{n}-2M^{{\tt T}}(MM^{{\tt T}})^{-1}M=I_{n}-2M^{{\tt
T}}(M^{{\tt T}})^{+}=I_{n}-2M^{+}M.$ (4.1)

So ${\mathcal{H}}_{M}$ is symmetric and
$M{\mathcal{H}}_{M}=-M,{\mathcal{H}}_{M}M^{+}=-M^{+}$.

Let $Z\subseteq[m]$, we set $A_{Z}=\sum_{i\in Z}{\bf e}_{i}A_{i}^{\tt
T}\in\mathbb{R}^{|Z|\times n}$ as the submatrix of $A$ whose rows are indexed
by $Z$. The corresponding block Householder matrix is ${\mathcal{H}}_{A_{Z}}$
or simply ${\mathcal{H}}_{Z}$ when it does not cause confusion. The iterative
procedure in block form we are interested in is

${\bf x}_{k+1}={\mathcal{H}}_{Z}{\bf x}_{k}+2A_{Z}^{+}{\bf b}_{Z}={\bf
x}_{k}-2A_{Z}^{+}(A_{Z}{\bf x}_{k}-{\bf b}_{Z}),\quad{\bf b}_{Z}=\sum_{i\in
Z}b_{i}{\bf e}_{i}\in\mathbb{R}^{|Z|}.$ (4.2)

The two algorithms, depending on the choice of $Z$, that we will consider have
been described in the introduction. For convenience, we restate them below.

  * •

Algorithm 1: Let $[m]=Z_{1}\cup\cdots\cup Z_{p}$ be a fixed partition of the
row indices. Arbitrarily choose ${\bf x}_{0}$, then update

${\bf x}_{k+1}={\bf x}_{k}-2A_{Z_{i}}^{+}A_{Z_{i}}{\bf
x}_{k}+2A_{Z_{i}}^{+}{\bf b}_{Z_{i}},\quad k=0,1,2,\ldots,$ (4.3)

where $i\in[p]$ is chosen with probability $\|A_{Z_{i}}\|_{\tt
F}^{2}/\|A\|_{\tt F}^{2}$.

  * •

Algorithm 2: Fix an integer $q\leq\min(m,n)$. Arbitrarily choose ${\bf
x}_{0}$, then update

${\bf x}_{k+1}={\bf x}_{k}-2A_{Z}^{+}A_{Z}{\bf x}_{k}+2A_{Z}^{+}{\bf
b}_{Z},\quad k=0,1,2,\ldots,$ (4.4)

where $Z=\\{i_{1},\ldots,i_{q}\\}\subseteq[m]$ is obtained in a way such that
$i_{j}$ is collected with probability $\|A_{i_{j}}\|^{2}/\|A\|_{\tt F}^{2}$.

In the following, we study the convergence rate of the above two algorithms
for consistent and inconsistent linear systems separately.

###  4.1 For consistent linear systems

In this subsection, we analyse the algorithm for consistent linear systems
$A{\bf x}={\bf b}$. Let ${\bf x}_{*}=A^{+}{\bf b}$ be a solution with minimum
norm and ${\bf e}_{k}={\bf x}_{k}-{\bf x}_{*}$ be the error at the $k$-th step
of iteration. Then it is easy to check that ${\bf
e}_{k+1}={\mathcal{H}}_{Z_{i}}{\bf e}_{k}$ when $A{\bf x}={\bf b}$ is
consistent.

We first analyse Algorithm 1.

######  Lemma 4.2.

Assume that $A$ has full column rank. Denote $\gamma=\min_{1\leq i\leq
p}\frac{\|A_{Z_{i}}\|_{\tt F}}{\|A_{Z_{i}}\|}$, then the operator norm of
${\mathbb{E}}_{i}[{\mathcal{H}}_{Z_{i}}]$ satisfies

$\|{\mathbb{E}}_{i}[{\mathcal{H}}_{Z_{i}}]\|\leq
1-\frac{2\gamma^{2}}{\kappa_{\tt F}^{2}}.$

###### Proof.

The result follows from a standard calculation:

$\displaystyle{\mathbb{E}}[{\mathcal{H}}_{Z_{i}}]$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{p}\frac{\|A_{Z_{i}}\|_{\tt F}^{2}}{\|A\|_{\tt
F}^{2}}(I_{n}-2A_{Z_{i}}^{+}A_{Z_{i}})$ $\displaystyle=$ $\displaystyle
I_{n}-\frac{2}{\|A\|_{\tt F}^{2}}\sum_{i=1}^{p}\|A_{Z_{i}}\|_{\tt
F}^{2}A_{Z_{i}}^{\tt T}(A_{Z_{i}}A_{Z_{i}}^{\tt T})^{-1}A_{Z_{i}}$
$\displaystyle\preceq$ $\displaystyle I_{n}-\frac{\min_{i}\|A_{Z_{i}}\|_{\tt
F}^{2}\sigma_{\min}((A_{Z_{i}}A_{Z_{i}}^{\tt T})^{-1})}{\|A\|_{\tt
F}^{2}}A^{\tt T}A$ $\displaystyle\preceq$
$\displaystyle\left(1-\frac{2}{\kappa_{\tt F}^{2}}\min_{1\leq i\leq
p}\frac{\|A_{Z_{i}}\|_{\tt F}^{2}}{\|A_{Z_{i}}\|^{2}}\right)I_{n}.$

In the first inequality, $\sigma_{\min}$ refers to the minimal singular value.
In the second inequality, we used the facts that
$\sigma_{\min}((A_{Z_{i}}A_{Z_{i}}^{\tt
T})^{-1})=\sigma_{\max}(A_{Z_{i}})^{-2}$ and $A^{\tt
T}A\succeq\sigma_{\min}(A^{\tt T}A)I_{n}$ since $A$ has full column rank. ∎

######  Lemma 4.3.

Assume that $A$ has full column rank, then for any $\ell\geq 0$, we have

$\big{|}{\mathbb{E}}[\langle{\bf e}_{k},{\bf
e}_{k+\ell}\rangle]\big{|}=\left(1-\frac{2\gamma^{2}}{\kappa_{\tt
F}^{2}}\right)^{\ell}\|{\bf e}_{0}\|^{2}.$

###### Proof.

The claim is trivial when $\ell=0$. We next consider the case that $\ell=1$.
It is easy to see from ${\bf e}_{k+1}={\mathcal{H}}_{Z_{i}}{\bf e}_{k}$ that
$\|{\bf e}_{k}\|=\|{\bf e}_{0}\|$ for all $k$. By Lemma 4.2, we have

$\displaystyle\big{|}{\mathbb{E}}[\langle{\bf e}_{k},{\bf
e}_{k+1}\rangle]\big{|}=\big{|}{\mathbb{E}}[\langle{\bf
e}_{k},{\mathbb{E}}_{i}[{\mathcal{H}}_{Z_{i}}]{\bf
e}_{k}\rangle]\big{|}\leq\|{\mathbb{E}}_{i}[{\mathcal{H}}_{Z_{i}}]\|\|{\bf
e}_{k}\|^{2}\leq(1-\frac{2\gamma^{2}}{\kappa_{\tt F}^{2}})\|{\bf
e}_{0}\|^{2}.$

Generally, when $\ell\geq 2$ we have

$\displaystyle{\mathbb{E}}[\langle{\bf e}_{k},{\bf e}_{k+\ell}\rangle]$
$\displaystyle=$ $\displaystyle{\mathbb{E}}[\langle{\bf
e}_{k},{\mathbb{E}}_{i}[{\mathcal{H}}_{Z_{i}}]{\bf e}_{k+\ell-1}\rangle]$
$\displaystyle=$ $\displaystyle{\mathbb{E}}[\langle{\bf
e}_{k},{\mathbb{E}}_{i}[{\mathcal{H}}_{Z_{i}}]^{\ell}{\bf e}_{k}\rangle]$
$\displaystyle\leq$
$\displaystyle\|{\mathbb{E}}_{i}[{\mathcal{H}}_{Z_{i}}]\|^{\ell}\|{\bf
e}_{k}\|^{2}$ $\displaystyle\leq$
$\displaystyle(1-\frac{2\gamma^{2}}{\kappa_{\tt F}^{2}})^{\ell}\|{\bf
e}_{0}\|^{2}.$

This finishes the proof. ∎

######  Theorem 4.4 (Convergence rate of Algorithm 1 for consistent linear
systems).

Assume that the linear system $A{\bf x}={\bf b}$ is consistent and $A$ has
full column rank. Let ${\bf x}_{*}=A^{+}{\bf b}$ be a solution and $\\{{\bf
x}_{0},{\bf x}_{1},{\bf x}_{2},\ldots\\}$ be a series of points generated by
the procedure (4.3). Then

${\mathbb{E}}\left[\,\left\|{\bf x}_{*}-\frac{1}{N}\sum_{k=0}^{N-1}{\bf
x}_{k}\right\|^{2}\,\right]\leq\frac{1}{N}\left(1+\frac{\kappa_{\tt
F}^{2}}{2\gamma^{2}}\right)\|{\bf x}_{*}-{\bf x}_{0}\|^{2}.$

###### Proof.

Recall that ${\bf e}_{k}={\bf x}_{k}-{\bf x}_{*}$. So

$\left\|\sum_{k=0}^{N-1}{\bf e}_{k}\right\|^{2}=\sum_{k=0}^{N-1}\|{\bf
e}_{k}\|^{2}+2\sum_{k=0}^{N-2}\sum_{\ell=k+1}^{N-1}\langle{\bf e}_{k},{\bf
e}_{l}\rangle=N\|{\bf
e}_{0}\|^{2}+2\sum_{k=0}^{N-2}\sum_{\ell=k+1}^{N-1}\langle{\bf e}_{k},{\bf
e}_{\ell}\rangle.$

Denote $q=1-\frac{2\gamma^{2}}{\kappa_{\tt F}^{2}}$. By Lemma 4.3, we have

${\mathbb{E}}\left[\,\sum_{k=0}^{N-2}\sum_{\ell=k+1}^{N-1}\langle{\bf
e}_{k},{\bf
e}_{\ell}\rangle\,\right]\leq\sum_{k=0}^{N-2}\sum_{\ell=k+1}^{N-1}q^{\ell-k}\|{\bf
e}_{0}\|^{2}\leq\frac{N}{1-q}\|{\bf e}_{0}\|^{2}.$

In the above estimation, we used an inequality in the proof of Theorem 3.3.
Putting it all together, we obtain the claimed bound. ∎

From the above theorem, to ensure a good approximation of ${\bf x}_{*}$, we
have $N\approx\kappa_{\tt F}^{2}/\gamma^{2}$. Generally, $\gamma\geq 1$. In
the standard version without using blocks, $\gamma=1$, so the above theorem
matches Theorem 3.3 when ${\bf c}=0$. The above theorem also explains why the
block version converges faster, especially when $\gamma\gg 1$.

Next, we analyse Algorithm 2.

######  Lemma 4.5.

Let $Z$ be a collection of $q$ row indices, then for any ${\bf u}$ that lies
in the column space of $A$, we have

${\mathbb{E}}_{Z}[\langle{\bf u}|{\mathcal{H}}_{Z}|{\bf
u}\rangle]\leq\left(1-\frac{2q}{\|A\|^{2}\|A^{+}\|^{4}\|A\|_{\tt
F}^{2}}\right)\|{\bf u}\|^{2}=\left(1-\frac{2q}{\kappa^{2}\kappa_{\tt
F}^{2}}\right)\|{\bf u}\|^{2}.$

Consequently, when $A$ has a full column rank, we have

$\|{\mathbb{E}}_{Z}[{\mathcal{H}}_{Z}]\|\leq 1-\frac{2q}{\kappa^{2}\kappa_{\tt
F}^{2}}.$

###### Proof.

We use $A_{Z}$ to denote the submatrix of $A$ whose rows are indexed by $Z$.
In our notation, ${\mathcal{H}}_{Z}$ refers to the block Householder matrix
${\mathcal{H}}_{A_{Z}}$. We use $\sigma_{\max}(A_{Z})$ to denote the maximal
singular value of $A_{Z}$. Note that

$\displaystyle\langle{\bf u}|{\mathcal{H}}_{Z}|{\bf u}\rangle=\|{\bf
u}\|^{2}-2\langle{\bf u}|A_{Z}^{{\tt T}}(A_{Z}A_{Z}^{{\tt T}})^{+}A_{Z}|{\bf
u}\rangle\leq\|{\bf u}\|^{2}-\frac{2}{\sigma_{\max}^{2}(A_{Z})}\|A_{Z}{\bf
u}\|^{2}.$

Regarding $\sigma_{\max}(A_{Z})$, notice that $A_{Z}A_{Z}^{{\tt
T}}=\sum_{i,j\in Z}\langle A_{i},A_{j}\rangle{\bf e}_{i}{\bf e}_{j}^{\tt T}$
is a submatrix of $AA^{{\tt T}}$, so we have
$\sigma_{\max}^{2}(A_{Z})=\|A_{Z}A_{Z}^{{\tt T}}\|\leq\|A\|^{2}.$ It is also
easy to see that

${\mathbb{E}}[\|A_{Z}{\bf u}\|^{2}]=\langle{\bf
u}|{\mathbb{E}}\left(\sum_{i\in Z}A_{i}^{\tt T}A_{i}\right)|{\bf
u}\rangle\geq\frac{|Z|\min_{i\in Z}\|A_{i}\|^{2}}{\|A\|_{\tt
F}^{2}}\langle{\bf u}|A^{\tt T}A|{\bf u}\rangle\geq\frac{|Z|\min_{i\in
Z}\|A_{i}\|^{2}}{\|A^{+}\|^{2}\|A\|_{\tt F}^{2}}\|{\bf u}\|^{2},$

where $|Z|=q$ is the size of $Z$. In the last inequality, we used the fact
that ${\bf u}$ lies in the column space of $A$. Putting it together, we obtain
the claimed result. ∎

In Lemma 4.2, we introduced a parameter $\gamma=\min_{1\leq i\leq
p}\frac{\|A_{Z_{i}}\|_{\tt F}}{\|A_{Z_{i}}\|}$. It is easy to see that
$\gamma^{2}\geq\min_{1\leq i\leq p}\frac{|Z_{i}|}{\kappa(Z_{i})^{2}}$. Here
$\kappa(Z_{i})$ is the condition number of $A_{Z_{i}}$, a submatrix of $A$ by
choosing rows of $A$ indexed by $Z_{i}$. In the proof of Lemma 4.5, we bounded
the operator norm of $A_{Z}$ via that of $A$, and so obtained
$|Z|/\kappa(A)^{2}$ in the end. These two correspond to each other as we
always have $\min_{1\leq i\leq
p}\frac{|Z_{i}|}{\kappa(Z_{i})^{2}}\geq\frac{\min_{1\leq i\leq
p}|Z_{i}|}{\kappa(A)^{2}}$. But still, the bound in Lemma 4.2 is slightly
better.

With the above lemma, we can prove a similar result to Lemma 4.3. With these
two results, we can similarly prove the following theorem, which is comparable
to Theorem 4.4.

######  Theorem 4.6 (Convergence rate of Algorithm 2 for consistent linear
systems).

Assume the linear system $A{\bf x}={\bf b}$ is consistent and $A$ has full
column rank. Let ${\bf x}_{*}=A^{+}{\bf b}$ be a solution and $\\{{\bf
x}_{0},{\bf x}_{1},{\bf x}_{2},\ldots\\}$ be a series of points generated by
the procedure (4.4). Then

${\mathbb{E}}\left[\,\left\|{\bf x}_{*}-\frac{1}{N}\sum_{k=0}^{N-1}{\bf
x}_{k}\right\|^{2}\,\right]\leq\frac{1}{N}\left(1+\frac{\kappa^{2}\kappa_{\tt
F}^{2}}{2q}\right)\|{\bf x}_{*}-{\bf x}_{0}\|^{2},$ (4.5)

where $q$ is the integer specified in (4.4).

###### Proof.

The proof is similar to that of Theorem 4.4. Note that

$\displaystyle{\mathbb{E}}\left[\,\left\|\sum_{k=0}^{N-1}{\bf
e}_{k}\right\|^{2}\,\right]$ $\displaystyle=$
$\displaystyle{\mathbb{E}}\left[\,\sum_{k=0}^{N-1}\|{\bf
e}_{k}\|^{2}\,\right]+2\sum_{k=0}^{N-2}\sum_{\ell=k+1}^{N-1}{\mathbb{E}}[\langle{\bf
e}_{k},{\bf e}_{\ell}\rangle]$ $\displaystyle=$ $\displaystyle N\|{\bf
e}_{0}\|^{2}+2\sum_{k=0}^{N-2}\sum_{\ell=k+1}^{N-1}{\mathbb{E}}[\langle{\bf
e}_{k},{\bf e}_{\ell}\rangle].$

By Lemma 4.5, for any $c\geq 1$, we have

$\displaystyle{\mathbb{E}}[\langle{\bf e}_{k}|{\bf e}_{k+c}\rangle]$
$\displaystyle=$ $\displaystyle{\mathbb{E}}\Big{[}{\mathbb{E}}_{Z}[\langle{\bf
e}_{k}|{\mathcal{P}}_{Z}|{\bf e}_{k+c-1}\rangle]\Big{]}$ $\displaystyle=$
$\displaystyle\langle{\bf
e}_{k}|({\mathbb{E}}_{Z}[{\mathcal{H}}_{Z}])^{c}|{\bf e}_{k}\rangle$
$\displaystyle\leq$
$\displaystyle\|{\mathbb{E}}_{Z}[{\mathcal{H}}_{Z}]\|^{c}\|{\bf e}_{k}\|^{2}$
$\displaystyle\leq$ $\displaystyle\left(1-\frac{2q}{\kappa^{2}\kappa_{\tt
F}^{2}}\right)\|{\bf e}_{0}\|^{2}.$

Therefore,

$\displaystyle{\mathbb{E}}\left[\,\left\|\frac{1}{N}\sum_{k=0}^{N-1}{\bf
e}_{k}\right\|^{2}\,\right]$ $\displaystyle\leq$
$\displaystyle\frac{1}{N}\|{\bf
e}_{0}\|^{2}+\frac{2}{N^{2}}\sum_{k=0}^{N-2}\sum_{\ell=k+1}^{N-1}\left(1-\frac{2q}{\kappa^{2}\kappa_{\tt
F}^{2}}\right)^{\ell-k}\|{\bf e}_{0}\|^{2}$ $\displaystyle\leq$
$\displaystyle\frac{1}{N}\left(1+\frac{\kappa^{2}\kappa_{\tt
F}^{2}}{2q}\right)\|{\bf e}_{0}\|^{2},$

which is as claimed. ∎

###  4.2 For inconsistent linear systems

In this subsection, we consider Algorithm 1 and Algorithm 2 for solving
inconsistent linear systems, i.e., for solving least squares. The error
analysis is similar but slightly harder than that of Theorems 3.3, 4.4 and
4.6. So we will make an extra but reasonable assumption in the analysis to
simplify the argument. We will explain the reason after the proof of Lemma
4.7. Our main results (i.e., Theorems 4.9 and 4.11) coincide with Theorems 4.4
and 4.6 when the linear system is consistent.

Let ${\bf x}_{*}=A^{+}{\bf b}$ be an optimal solution of $\arg\min\|A{\bf
x}-{\bf b}\|$, then there is a ${\bf c}$ such that ${\bf b}$ admits an
orthogonal decomposition ${\bf b}=A{\bf x}_{*}+{\bf c}$. For any
$Z\subseteq[m]$, we also have ${\bf b}_{Z}=A_{Z}{\bf x}_{*}+{\bf c}_{Z}$. Here
${\bf b}_{Z},{\bf c}_{Z}$ are the sub-vectors of ${\bf b},{\bf c}$
respectively restricted to indices in $Z$, see (4.2). Let the error at the
$k$-th step of iteration be ${\bf e}_{k}={\bf x}_{k}-{\bf x}_{*}$.

As usual, we start from Algorithm 1.

######  Lemma 4.7.

Assume that ${\mathbb{E}}[\|{\bf e}_{k}\|]\leq E$ for all $k$. Let
$R=\frac{\|{\bf c}\|}{\|A\|_{\tt F}}\max_{i}\|A_{Z_{i}}\|_{\tt
F}\|A_{Z_{i}}^{+}\|$, then

$\displaystyle{\mathbb{E}}[\|{\bf e}_{k}\|^{2}]$ $\displaystyle\leq$
$\displaystyle\|{\bf e}_{0}\|^{2}+4kR(R+E),$ (4.6)
$\displaystyle{\mathbb{E}}[\langle{\bf e}_{k},{\bf e}_{k+\ell}\rangle]$
$\displaystyle\leq$ $\displaystyle(1-\frac{2\gamma^{2}}{\kappa_{\tt
F}^{2}})^{\ell}\|{\bf e}_{0}\|^{2}+4kR(R+E)(1-\frac{2\gamma^{2}}{\kappa_{\tt
F}^{2}})^{\ell}+\frac{RE\kappa_{\tt
F}^{2}}{\gamma^{2}}\left(1-(1-\frac{2\gamma^{2}}{\kappa_{\tt
F}^{2}})^{\ell}\right),$ (4.7)

where $\gamma$ is defined in Lemma 4.2.

###### Proof.

From (4.3), we have

${\bf e}_{k+1}={\bf x}_{k}-{\bf x}_{*}-2A_{Z_{i}}^{+}A_{Z_{i}}{\bf
x}_{k}+2A_{Z_{i}}^{+}(A_{Z_{i}}{\bf x}_{*}+{\bf
c}_{Z_{i}})={\mathcal{H}}_{Z_{i}}{\bf e}_{k}+2A_{Z_{i}}^{+}{\bf c}_{Z_{i}}.$

Note that

$\displaystyle{\mathbb{E}}_{i}\left[\|A_{Z_{i}}^{+}{\bf
c}_{Z_{i}}\|^{2}\right]$ $\displaystyle=$
$\displaystyle{\mathbb{E}}_{i}\left[{\bf c}_{Z_{i}}^{\tt
T}(A_{Z_{i}}A_{Z_{i}}^{\tt T})^{-1}{\bf c}_{Z_{i}}\right]$ $\displaystyle=$
$\displaystyle\frac{\sum_{i}\|A_{Z_{i}}\|_{\tt F}^{2}}{\|A\|_{\tt F}^{2}}{\bf
c}_{Z_{i}}^{\tt T}(A_{Z_{i}}A_{Z_{i}}^{\tt T})^{-1}{\bf c}_{Z_{i}}$
$\displaystyle\leq$ $\displaystyle\frac{\|{\bf c}\|^{2}}{\|A\|_{\tt
F}^{2}}\max_{i}\|A_{Z_{i}}\|_{\tt F}^{2}\|A_{Z_{i}}^{+}\|^{2}.$

Namely, ${\mathbb{E}}_{i}\left[\|A_{Z_{i}}^{+}{\bf c}_{Z_{i}}\|^{2}\right]\leq
R^{2}$. As ${\mathcal{H}}_{Z_{i}}A_{Z_{i}}^{+}=-A_{Z_{i}}^{+}$, we have

$\displaystyle{\mathbb{E}}[\|{\bf e}_{k+1}\|^{2}]$ $\displaystyle=$
$\displaystyle{\mathbb{E}}\Big{[}{\mathbb{E}}_{i}[\|{\bf
e}_{k}-2A_{Z_{i}}^{+}{\bf c}_{Z_{i}}\|^{2}]\Big{]}$ $\displaystyle=$
$\displaystyle{\mathbb{E}}\Big{[}\|{\bf
e}_{k}\|^{2}+4{\mathbb{E}}_{i}[\|A_{Z_{i}}^{+}{\bf
c}_{Z_{i}}\|^{2}]-4\langle{\bf e}_{k},{\mathbb{E}}_{i}[A_{Z_{i}}^{+}{\bf
c}_{Z_{i}}]\rangle\Big{]}$ $\displaystyle\leq$
$\displaystyle{\mathbb{E}}[\|{\bf e}_{k}\|^{2}]+(4R^{2}+4RE).$

Inductively, we obtain the first claim.

We below prove the second claim. Denote
${\mathcal{H}}={\mathbb{E}}_{i}[{\mathcal{H}}_{Z_{i}}]$ for convenience. Then
by Lemma 4.2, $\rho:=\|{\mathcal{H}}\|\leq 1-\frac{2\gamma^{2}}{\kappa_{\tt
F}^{2}}.$ We can compute that

$\displaystyle{\mathbb{E}}[\langle{\bf e}_{k},{\bf e}_{k+\ell}\rangle]$
$\displaystyle=$ $\displaystyle{\mathbb{E}}\Big{[}\Big{\langle}{\bf
e}_{k},{\mathbb{E}}_{i}[{\mathcal{H}}_{Z_{i}}{\bf
e}_{k+\ell-1}+2A_{Z_{i}}^{+}{\bf c}_{Z_{i}}]\Big{\rangle}\Big{]}$
$\displaystyle\leq$ $\displaystyle{\mathbb{E}}[\langle{\bf
e}_{k}{\mathcal{H}},{\bf e}_{k+\ell-1}\rangle]+2RE$ $\displaystyle=$
$\displaystyle{\mathbb{E}}\Big{[}\Big{\langle}{\bf
e}_{k}{\mathcal{H}},{\mathbb{E}}_{i}[{\mathcal{H}}_{Z_{i}}{\bf
e}_{k+\ell-2}+2A_{Z_{i}}^{+}{\bf c}_{Z_{i}}]\Big{\rangle}\Big{]}+2RE$
$\displaystyle\leq$ $\displaystyle{\mathbb{E}}[\langle{\bf
e}_{k}{\mathcal{H}}^{2},{\bf e}_{k+\ell-1}\rangle]+2RE(1+\rho).$

Continuous this process, we will end up with

$\displaystyle{\mathbb{E}}[\langle{\bf e}_{k},{\bf e}_{k+\ell}\rangle]$
$\displaystyle\leq$ $\displaystyle{\mathbb{E}}[\langle{\bf
e}_{k}{\mathcal{H}}^{\ell},{\bf
e}_{k}\rangle]+2RE(1+\rho+\cdots+\rho^{\ell-1})$ $\displaystyle\leq$
$\displaystyle\rho^{\ell}{\mathbb{E}}[\|{\bf
e}_{k}\|^{2}]+2RE\frac{1-\rho^{\ell}}{1-\rho}.$

By the first claim, we obtain

$\displaystyle{\mathbb{E}}[\langle{\bf e}_{k},{\bf
e}_{k+\ell}\rangle]\leq\rho^{\ell}(\|{\bf
e}_{0}\|^{2}+4kR(R+E))+2RE\frac{1-\rho^{\ell}}{1-\rho}.$

Simplifying this leads to the second claim. ∎

######  Remark 4.8.

(1). In the proof of the above lemma, we may not have
${\mathbb{E}}_{i}[A_{Z_{i}}^{+}{\bf c}_{Z_{i}}]=0$. This is in sharp contrast
to Lemma 3.1 (more precisely (3.4)) or the block version using orthogonal
projections [17]. This makes the analysis harder. So in the above upper bound
analysis of $\|{\bf e}_{k+1}\|^{2}$, we used a quite large bound of $\|{\bf
e}_{k}\|^{2}+4R(R+E)$ for simplicity. But when return to the situation of non-
block version, we can remove the assumption ${\mathbb{E}}[\|{\bf e}_{k}\|]\leq
E$ and set $E=0$ directly in the above proof. Consequently, it coincides with
Lemma 3.1.

(2). The assumption ${\mathbb{E}}[\|{\bf e}_{k}\|]\leq E$ in Lemma 4.7 may be
a little strange at first glance. However, the following two reasons indicate
the rationality of this assumption. First, we have the following very coarse
upper bound

${\mathbb{E}}[\|{\bf e}_{k}\|]={\mathbb{E}}[\|{\bf e}_{k-1}-2A_{Z_{i}}^{+}{\bf
c}_{Z_{i}}\|]\leq{\mathbb{E}}[\|{\bf e}_{k-1}\|]+2R\leq\|{\bf e}_{0}\|+2kR.$

Namely, we can replace $E$ with $E_{k}:=|{\bf e}_{0}\|+2kR$ in the above
estimate, which depends on $k$ now. But in practice, the error should decrease
on average. The replacement will affect Theorem 4.9 below by a factor of $N$
for the error term (i.e., Err), which is acceptable when $\|{\bf c}\|\approx
0$. Second, the bounds on the error terms in Lemma 4.7 are affected by $R$,
which further depends on $\|{\bf c}\|$. When $\|{\bf c}\|$ is small, i.e.,
when ${\bf b}$ is close to the column space of $A$, then $R$ can be very
small. As we can see from the error term below, $E$ is multiplied by $R$ too.
So in this case, the overall error can still be very small even if $E$ is
replaced by $\|{\bf e}_{0}\|+2kR$ at the $k$-th step of iteration.

######  Theorem 4.9 (Convergence rate of Algorithm 1 for inconsistent linear
systems).

Assume that $A$ has full column rank. Let ${\bf x}_{*}=A^{+}{\bf b}$ and
$\\{{\bf x}_{0},{\bf x}_{1},{\bf x}_{2},\ldots,{\bf x}_{N-1}\\}$ be a series
of points generated by the procedure (4.3). Assume that ${\mathbb{E}}[\|{\bf
x}_{*}-{\bf x}_{k}\|]\leq E$ for all $k\leq N-1$. Let $R=\frac{\|{\bf
c}\|}{\|A\|_{\tt F}}\max_{i}\|A_{Z_{i}}\|_{\tt F}\|A_{Z_{i}}^{+}\|$. Then

${\mathbb{E}}\left[\,\left\|{\bf x}_{*}-\frac{1}{N}\sum_{k=0}^{N-1}{\bf
x}_{k}\right\|^{2}\,\right]\leq\frac{1}{N}\left(1+\frac{\kappa_{\tt
F}^{2}}{2\gamma^{2}}\right)\|{\bf x}_{*}-{\bf x}_{0}\|^{2}+{\rm Err},$

where ${\rm Err}=O(R(R+E)\kappa_{\tt F}^{2}/\gamma^{2})$ and $\gamma$ is
defined in Lemma 4.2.

###### Proof.

Recall that ${\bf e}_{k}={\bf x}_{k}-{\bf x}_{*}$ and

$\displaystyle{\mathbb{E}}\left[\,\left\|\sum_{k=0}^{N-1}{\bf
e}_{k}\right\|^{2}\,\right]=\sum_{k=0}^{N-1}{\mathbb{E}}[\|{\bf
e}_{k}\|^{2}]+2\sum_{k=0}^{N-2}\sum_{\ell=1}^{N-k-1}{\mathbb{E}}[\langle{\bf
e}_{k},{\bf e}_{k+\ell}\rangle].$

For the first term, by Equation (4.6), we have

$\sum_{k=0}^{N-1}{\mathbb{E}}[\|{\bf e}_{k}\|^{2}]\leq N\|{\bf
e}_{0}\|^{2}+2N^{2}R(R+E).$

For the second term, by Equation (4.7), we have

$\displaystyle\sum_{k=0}^{N-2}\sum_{\ell=1}^{N-k-1}{\mathbb{E}}[\langle{\bf
e}_{k},{\bf e}_{k+\ell}\rangle]$ $\displaystyle\leq$
$\displaystyle\sum_{k=0}^{N-2}\sum_{\ell=1}^{N-k-1}(1-\frac{2\gamma^{2}}{\kappa_{\tt
F}^{2}})^{\ell}\|{\bf e}_{0}\|^{2}$
$\displaystyle+\sum_{k=0}^{N-2}\sum_{\ell=1}^{N-k-1}4kR(R+E)(1-\frac{2\gamma^{2}}{\kappa_{\tt
F}^{2}})^{\ell}$
$\displaystyle+\sum_{k=0}^{N-2}\sum_{\ell=1}^{N-k-1}\frac{RE\kappa_{\tt
F}^{2}}{\gamma^{2}}\left(1-(1-\frac{2\gamma^{2}}{\kappa_{\tt
F}^{2}})^{\ell}\right).$

From the proof of Theorem 3.3, it is easy to estimate that

$\sum_{k=0}^{N-2}\sum_{\ell=1}^{N-k-1}(1-\frac{2\gamma^{2}}{\kappa_{\tt
F}^{2}})^{\ell}\leq\frac{N\kappa_{\tt F}^{2}}{2\gamma^{2}}.$

Similarly, we have

$\sum_{k=0}^{N-2}\sum_{\ell=1}^{N-k-1}4kR(R+E)(1-\frac{2\gamma^{2}}{\kappa_{\tt
F}^{2}})^{\ell}\leq 2R(R+E)\frac{N^{2}\kappa_{\tt F}^{2}}{\gamma^{2}},$

and

$\sum_{k=0}^{N-2}\sum_{\ell=1}^{N-k-1}\frac{RE\kappa_{\tt
F}^{2}}{\gamma^{2}}\left(1-(1-\frac{2\gamma^{2}}{\kappa_{\tt
F}^{2}})^{\ell}\right)\leq\frac{N^{2}RE\kappa_{\tt F}^{2}}{\gamma^{2}}.$

Finally, putting it all together, we obtain the claimed result. ∎

In the end, we consider Algorithm 2. The proofs are similar. Now we have to
use Lemma 4.3 to compute a bound of the operator norm. Basically, we have the
following result. The proof is similar to that of Lemma 4.7.

######  Lemma 4.10.

Assume that ${\mathbb{E}}[\|{\bf e}_{k}\|]\leq E$ for all $k$. Let $q$ be an
integer specified in Algorithm 2. Let

$R=\frac{\sqrt{q}\|A\|\|{\bf c}\|}{\|A\|_{\tt
F}}\max_{Z\subseteq[m]:|Z|=q}\|(A_{Z}A_{Z})^{+}\|^{1/2}.$ (4.8)

then

$\displaystyle{\mathbb{E}}[\|{\bf e}_{k}\|^{2}]$ $\displaystyle\leq$
$\displaystyle\|{\bf e}_{0}\|^{2}+4kR(R+E),$
$\displaystyle{\mathbb{E}}[\langle{\bf e}_{k},{\bf e}_{k+\ell}\rangle]$
$\displaystyle\leq$ $\displaystyle(1-\frac{2q}{\kappa^{2}\kappa_{\tt
F}^{2}})^{\ell}\|{\bf e}_{0}\|^{2}+4kR(R+E)(1-\frac{2q}{\kappa^{2}\kappa_{\tt
F}^{2}})^{\ell}+\frac{RE\kappa_{\tt
F}^{2}}{\gamma}\left(1-(1-\frac{2q}{\kappa^{2}\kappa_{\tt
F}^{2}})^{\ell}\right).$

Finally, similar to the proof of Theorem 4.9, by using Lemma 4.10 we will
obtain the following result.

######  Theorem 4.11 (Convergence rate of Algorithm 2 for inconsistent linear
systems).

Assume that $A$ has full column rank. Let ${\bf x}_{*}=A^{+}{\bf b}$ and
$\\{{\bf x}_{0},{\bf x}_{1},{\bf x}_{2},\ldots,{\bf x}_{N-1}\\}$ be a series
of points generated by the procedure (4.4). Assume that ${\mathbb{E}}[\|{\bf
x}_{*}-{\bf x}_{k}\|]\leq E$ for all $k\leq N-1$. Let $R$ be given in (4.8)
and $q$ be an integer specified in Algorithm 2, then

${\mathbb{E}}\left[\left\|{\bf x}_{*}-\frac{1}{N}\sum_{k=0}^{N-1}{\bf
x}_{k}\right\|^{2}\right]\leq\frac{1}{N}\left(1+\frac{\kappa_{\tt
F}^{2}\kappa^{2}}{2q}\right)\|{\bf x}_{*}-{\bf x}_{0}\|^{2}+{\rm Err},$

where ${\rm Err}=O(R(R+E)\kappa_{\tt F}^{2}\kappa^{2}/q)$.

##  5 Deterministic reflective Kaczmarz algorithm

In [21], the author studied a deterministic version of the algorithm (3.1),
where $i_{k}$ is chosen cyclically, i.e., $i_{k}=(k\mod m)+1$. Some
interesting and new geometric properties about solutions of linear systems
were discovered. Roughly, the two main results of [21] are: (1). Let $\\{{\bf
x}_{0},{\bf x}_{1},{\bf x}_{2},\ldots\\}$ be the points generated by (3.1)
with $i_{k}=(k\mod m)+1$, then the average of $\\{{\bf
x}_{jk}:j=0,1,\ldots,N-1\\}$ is an approximation of the solution ${\bf x}_{*}$
of $A{\bf x}={\bf b}$ that is closest to ${\bf x}_{0}$ when it is consistent.
To ensure a good approximation, $N$ is determined by a quantity, denoted as
$\eta(A)$ in [21], also see Definition 5.15. In [21], it was proved that
$\eta(A)=O(\kappa(A)^{2})$. (2). Moreover, $\\{{\bf
x}_{2jk}:j=0,1,2,\ldots\\}\cup\\{{\bf x}_{(2j+1)k}:j=0,1,2,\ldots\\}$ lie on
two spheres such that the midpoint of the two centers of spheres is ${\bf
x}_{*}$ exactly.

In this section, as an continuation of this research, we investigate more
properties of this algorithm. First, we study the block version of the
algorithm. We will show that all results proved in [21] also hold in the block
version. Second, as a theoretical complement to this algorithm, we further
study more properties about the convergence rate (i.e., $\eta(A)$) and its
connections to least squares.

###  5.1 Block version

Let $A\in\mathbb{R}^{m\times n}$. As discussed above (similar to Algorithm 1),
in the block version we decompose $A$ into several parts according to a row
partition of $A$. More precisely, let $Z=\\{Z_{1},\ldots,Z_{p}\\}$ be a set of
submatrices that form a row partition of $A$, i.e., each $Z_{i}$ is a
submatrix of $A$ generated by some rows of $A$. So there is a permutation
matrix $S\in\mathbb{R}^{m\times m}$ such that

$SA=\begin{pmatrix}Z_{1}\\\ Z_{2}\\\ \vdots\\\ Z_{p}\end{pmatrix}.$ (5.1)

The partition $Z$ will be called a row partition of $A$ for convenience.

######  Remark 5.1.

In the following analysis, for convenience, we shall set $S=I_{m}$ and assume
that the partition is contiguous. All results stated below hold for a general
permutation $S$ and when the partition is not contiguous.

In this section, we mainly consider the following deterministic block
algorithm for solving linear systems:

${\bf x}_{k+1}={\bf x}_{k}-2A_{Z_{i}}^{+}A_{Z_{i}}{\bf
x}_{k}+2A_{Z_{i}}^{+}{\bf b}_{Z_{i}},\quad k=0,1,2,\ldots,$ (5.2)

where $i=(k\mod p)+1$ is chosen cyclically.

From (5.1), when $S=I_{m}$, it is easy to see that

$AA^{{\tt T}}=\begin{pmatrix}\vspace{.1cm}Z_{1}Z_{1}^{{\tt
T}}&Z_{1}Z_{2}^{{\tt T}}&\cdots&Z_{1}Z_{p}^{{\tt T}}\\\
\vspace{.1cm}Z_{2}Z_{1}^{{\tt T}}&Z_{2}Z_{2}^{{\tt T}}&\cdots&Z_{2}Z_{p}^{{\tt
T}}\\\ \vspace{.1cm}\vdots&\vdots&\ddots&\vdots\\\ Z_{p}Z_{1}^{{\tt
T}}&Z_{p}Z_{2}^{{\tt T}}&\cdots&Z_{k}Z_{p}^{{\tt
T}}\end{pmatrix}\in\mathbb{R}^{m\times m}.$

A matrix $X$ is called a block-partition matrix of $A$ with respect to
partition $Z$ if it has the form

$X=\begin{pmatrix}\vspace{.1cm}Z_{1}Z_{1}^{{\tt T}}&0&\cdots&0\\\
\vspace{.1cm}2Z_{2}Z_{1}^{{\tt T}}&Z_{2}Z_{2}^{{\tt T}}&\cdots&0\\\
\vspace{.1cm}\vdots&\vdots&\ddots&\vdots\\\ 2Z_{p}Z_{1}^{{\tt
T}}&2Z_{p}Z_{2}^{{\tt T}}&\cdots&Z_{p}Z_{p}^{{\tt
T}}\end{pmatrix}\in\mathbb{R}^{m\times m}.$ (5.3)

Basically, $X$ is a blocked lower triangular matrix satisfying $X+X^{{\tt
T}}=2AA^{{\tt T}}$. It is easy to see that
$\det(X)=\prod_{i=1}^{p}\det(Z_{i}Z_{i}^{\tt T})$. So $X$ is nonsingular when
all $Z_{i}$ have full row rank.

######  Lemma 5.2.

Let $X$ be a block-partition matrix of $A$ of the form (5.3). We decompose it
as

$X=\begin{pmatrix}\vspace{.1cm}X_{1}&0\\\ Z_{p}X_{2}&Z_{p}Z_{p}^{\tt
T}\end{pmatrix},$

where $X_{1}$ is the top-left corner containing the first $(p-1)\times(p-1)$
blocks, then the pseudoinverse of $X$ satisfies

$X^{+}=\begin{pmatrix}\vspace{.1cm}X_{1}^{+}&0\\\ -(Z_{p}Z_{p}^{\tt
T})^{+}Z_{p}X_{2}X_{1}^{+}&(Z_{p}Z_{p}^{\tt T})^{+}\end{pmatrix}.$

Moreover,

$XX^{+}=\begin{pmatrix}X_{1}X_{1}^{+}&0\\\ 0&(Z_{p}Z_{p}^{\tt
T})(Z_{p}Z_{p}^{\tt T})^{+}\end{pmatrix},\quad
X^{+}X=\begin{pmatrix}X_{1}^{+}X_{1}&0\\\ 0&(Z_{p}Z_{p}^{\tt
T})^{+}(Z_{p}Z_{p}^{\tt T})\end{pmatrix}.$

Consequently, we have $X_{2}(I-X_{1}^{+}X_{1})=0$ and $(I-XX^{+})A=0$.

###### Proof.

We prove this lemma by induction on $k$. The claim is obviously true if $k=1$.
We below prove the general case. It is easy to compute that

$XX^{+}=\begin{pmatrix}\vspace{.1cm}X_{1}X_{1}^{+}&0\\\
[Z_{p}-(Z_{k}Z_{k}^{\tt T})(Z_{k}Z_{k}^{\tt
T})^{+}Z_{p}]X_{2}X_{1}^{+}&(Z_{p}Z_{p}^{\tt T})(Z_{p}Z_{p}^{\tt
T})^{+}\end{pmatrix}=\begin{pmatrix}\vspace{.1cm}X_{1}X_{1}^{+}&0\\\
0&(Z_{p}Z_{p}^{\tt T})(Z_{p}Z_{p}^{\tt T})^{+}\end{pmatrix}.$

In the above, we used the fact that $M-(MM^{\tt T})(MM^{\tt T})^{+}M=0$ for
any matrix $M$.

We also have

$X^{+}X=\begin{pmatrix}\vspace{.1cm}X_{1}^{+}X_{1}&0\\\ (Z_{p}Z_{p}^{\tt
T})^{+}Z_{p}X_{2}(I-X_{1}^{+}X_{1})&(Z_{p}Z_{p}^{\tt T})^{+}(Z_{p}Z_{p}^{\tt
T})\end{pmatrix}=\begin{pmatrix}\vspace{.1cm}X_{1}^{+}X_{1}&0\\\
0&(Z_{p}Z_{p}^{\tt T})^{+}(Z_{p}Z_{p}^{\tt T})\end{pmatrix}$

due to $X_{2}(I-X_{1}^{+}X_{1})=0$. This indeed can be proved by induction at
the same time. To be more exact, we denote

$X_{1}=\begin{pmatrix}\vspace{.1cm}Y_{1}&0\\\ Z_{p-1}Y_{2}&Z_{p-1}Z_{p-1}^{\tt
T}\end{pmatrix}.$

By induction, we have $Y_{2}(I-Y_{1}^{+}Y_{1})=0$ and

$X_{1}^{+}=\begin{pmatrix}\vspace{.1cm}Y_{1}^{+}&0\\\ -(Z_{p-1}Z_{p-1}^{\tt
T})^{+}Z_{p-1}Y_{2}Y_{1}^{+}&(Z_{p-1}Z_{p-1}^{\tt T})^{+}\end{pmatrix}.$

Note that $X_{2}=2(Z_{1}^{\tt T},\cdots,Z_{p-2}^{\tt T},Z_{p-1}^{\tt
T})=(Y_{2},2Z_{p-1}^{\tt T})$, so

$\displaystyle X_{2}(I-X_{1}^{+}X_{1})$ $\displaystyle=$
$\displaystyle(Y_{2},2Z_{p-1}^{\tt
T})\begin{pmatrix}\vspace{.1cm}I-Y_{1}^{+}Y_{1}&0\\\ 0&I-(Z_{p-1}Z_{p-1}^{\tt
T})^{+}\end{pmatrix}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\vspace{.1cm}Y_{2}(I-Y_{1}^{+}Y_{1})\\\
2Z_{p-1}^{\tt T}(I-(Z_{p-1}Z_{p-1}^{\tt T})^{+})\end{pmatrix}=0.$

It is now easy to verify the properties of pseudoinverse for $X$, namely,
$XX^{+}X=X,X^{+}XX^{+}=X^{+}$, and $XX^{+},X^{+}X$ are Hermitian.

We now denote $A=\begin{pmatrix}A_{1}\\\ Z_{k}\end{pmatrix}$, where $A_{1}$
corresponds to $X_{1}$. By induction we have

$XX^{+}A=\begin{pmatrix}\vspace{.1cm}X_{1}X_{1}^{+}A_{1}\\\ (Z_{p}Z_{p}^{\tt
T})(Z_{p}Z_{p}^{\tt
T})^{+}Z_{p}\end{pmatrix}=\begin{pmatrix}\vspace{.1cm}A_{1}\\\
Z_{p}\end{pmatrix}=A.$

This completes the proof. ∎

######  Definition 5.3.

Let $Z=\\{Z_{1},\ldots,Z_{p}\\}$ be a row partition of $A$. Define

${\mathcal{H}}_{Z}(A):={\mathcal{H}}_{Z_{p}}\cdots{\mathcal{H}}_{Z_{2}}{\mathcal{H}}_{Z_{1}}$
(5.4)

as the product of block Householder matrices. When it makes no confusion, we
sometimes simply denote it as ${\mathcal{H}}(A)$. Note that
${\mathcal{H}}_{A}$ and ${\mathcal{H}}(A)$ are different.

######  Lemma 5.4.

Let $Z=\\{Z_{1},\ldots,Z_{p}\\}$ be a contiguous row partition of $A$. Then
${\mathcal{H}}_{Z}(A)=I_{n}-2A^{{\tt T}}X^{+}A$, where $X$ is the block-
partition matrix of $A$ with respect to $Z$ defined by (5.3).

###### Proof.

We prove the result by induction on $p$. For generality, we include $S$ in the
proof. If $p=1$, then $S=I,Z_{1}=A,X=Z_{1}Z_{1}^{\tt T}$, the claim is
obvious. Next, let $Z_{A},Z_{B}$ be partitions of $A,B$ respectively, then by
induction ${\mathcal{H}}_{Z}(A)=I_{n}-2(S_{1}A)^{{\tt T}}X_{A}^{+}(S_{1}A)$
and ${\mathcal{H}}_{Z}(B)=I_{n}-2(S_{2}B)^{{\tt T}}X_{B}^{+}(S_{2}B)$ for some
permutations $S_{1},S_{2}$ and two block-partition matrices $X_{A},X_{B}$.
Denote $C=\begin{pmatrix}A\\\ B\end{pmatrix},S=\begin{pmatrix}S_{1}&\\\
&S_{2}\end{pmatrix}$, then

$2(SC)(SC)^{\tt T}=\begin{pmatrix}\vspace{.1cm}2(S_{1}A)(S_{1}A)^{\tt
T}&2(S_{1}A)(S_{2}B)^{\tt T}\\\ 2(S_{2}B)(S_{1}A)^{\tt
T}&2(S_{2}B)(S_{2}B)^{\tt
T}\end{pmatrix}=\begin{pmatrix}\vspace{.1cm}X_{A}+X_{A}^{{\tt
T}}&2(S_{1}A)(S_{2}B)^{\tt T}\\\ 2(S_{2}B)(S_{1}A)^{\tt T}&X_{B}+X_{B}^{{\tt
T}}\end{pmatrix}.$

By Lemma 5.2, the block-partition matrix, which is denoted as $X_{C}$, of $C$
satisfies

$X_{C}=\begin{pmatrix}\vspace{.1cm}X_{A}&0\\\ 2(S_{2}B)(S_{1}A)^{\tt
T}&X_{B}\end{pmatrix},\quad
X_{C}^{+}=\begin{pmatrix}\vspace{.1cm}X_{A}^{+}&0\\\
-2X_{B}^{+}(S_{2}B)(S_{1}A)^{\tt T}X_{A}^{+}&X_{B}^{+}\end{pmatrix}.$

Now it is easy to check that
${\mathcal{H}}_{Z}(C)={\mathcal{H}}_{Z}(B){\mathcal{H}}_{Z}(A)=I_{n}-2(SC)^{\tt
T}X_{C}^{+}(SC)$. ∎

######  Remark 5.5.

When each $Z_{i}$ contains only one row of $A$, the decomposition of
${\mathcal{H}}_{Z}(A)$ is known as the $WY$ representation for the product of
Householder transformations [20]. In our notation, $W=A^{\tt T}X^{+},Y=A$. The
decomposition we obtained corresponds to the $Y^{\tt T}TY$ decomposition
mentioned in [20], i.e., $T=X^{+}$ in our notation.

######  Definition 5.6.

Fix a partition $Z=\\{Z_{1},\ldots,Z_{p}\\}$ of the rows of $A$. Assume that
the eigenvalues of ${\mathcal{H}}_{Z}(A)$ are $e^{2i\theta_{l}}$ with
$\theta_{l}\in[0,\pi)$ for $l\in\\{1,2,\ldots,n\\}$. Define

$\eta_{Z}(A):=\frac{1}{\min_{\ell:\theta_{\ell}\neq 0}|\sin\theta_{\ell}|}.$
(5.5)

Some properties of $\eta(A)$ when each $Z_{i}$ contains one row of $A$ has
been given in [21]. In the next subsection, we shall study more properties
about $\eta_{Z}(A)$. The proof of the following theorem is the same as that of
[21, Theorem 3.4].

######  Theorem 5.7.

Assume that the linear system $A{\bf x}={\bf b}$ is consistent. Fix a
partition $Z=\\{Z_{1},\ldots,Z_{p}\\}$ of the rows of $A$. Let ${\bf
x}_{0},{\bf x}_{1},\ldots,{\bf x}_{N-1}$ be the vectors generated by the
process (5.2). Then

$\left\|{\bf x}_{*}-\frac{1}{N}\sum_{j=0}^{N-1}{\bf
x}_{jk}\right\|\leq\varepsilon\|{\bf x}_{*}-{\bf x}_{0}\|,$ (5.6)

where ${\bf x}_{*}$ is the solution that has the minimal distance to ${\bf
x}_{0}$ and

$N=\left\lceil\frac{\pi\eta_{Z}(A)}{\varepsilon}\right\rceil.$ (5.7)

Suppose $Z=\\{Z_{1},Z_{2}\\}$ and $Z_{1}$ has the full column rank so that
${\mathcal{H}}_{Z_{1}}=-I$, then the set of $S_{1}=\\{{\bf
x}_{2j}:j=0,1,\ldots,\\}=\\{{\bf x}_{0},{\bf x}_{*}-{\mathcal{H}}_{Z_{2}}({\bf
x}_{0}-{\bf x}_{*})\\}$. The other set is $S_{2}=\\{{\bf
x}_{2j+1}:j=0,1,\ldots,\\}=\\{2{\bf x}_{*}-{\bf x}_{0},{\bf
x}_{*}+{\mathcal{H}}_{Z_{2}}({\bf x}_{0}-{\bf x}_{*})\\}$. Now $S_{1}$ or
$S_{2}$ alone is not enough for us to compute ${\bf x}_{*}$. However, their
union is enough.

As inspired by our previous research [21], we hope the algorithm (5.2)
achieves the following: if we start from ${\bf x}_{0}$, then the algorithm
will always converge to the solution ${\bf x}_{*}$ that is closest to ${\bf
x}_{0}$. As shown in [21], a necessary and sufficient condition for this is
that ${\bf x}_{0}-{\bf x}_{*}$ has no overlap on the eigenspace of
${\mathcal{H}}_{Z}(A)$ corresponding to the eigenvalue 1. We will prove this
in the block version and provide a more rigorous analysis.

######  Proposition 5.8.

Assume that $A{\bf x}={\bf b}$ is consistent. Let $\\{{\bf x}_{0},{\bf
x}_{1},\ldots\\}$ be the vectors generated by the process (5.2), and let ${\bf
x}_{*}$ be the solution that has the minimal distance to ${\bf x}_{0}$. Then
the following statements are equivalent.

  1. 1.

$\lim_{N\rightarrow\infty}\frac{1}{N}\sum_{k=0}^{N-1}{\bf x}_{k}={\bf x}_{*}$.

  2. 2.

${\bf x}_{0}-{\bf x}_{*}$ has no overlap on the eigenspace of
${\mathcal{H}}_{Z}(A)$ of eigenvalue 1.

  3. 3.

For any ${\bf v}$,

${\mathcal{H}}_{Z}(A){\bf v}={\bf v}\Leftrightarrow A{\bf v}=0.$ (5.8)

###### Proof.

$(1\Leftrightarrow 2)$ Suppose ${\bf y}_{0}:={\bf x}_{0}-{\bf x}_{*}={\bf
u}+{\bf v}$, where ${\bf u}$ lies in the eigenspace of ${\mathcal{H}}_{Z}(A)$
of eigenvalues not equal 1, and ${\bf v}$ lies in the eigenspace of
${\mathcal{P}}_{Z,A}$ of eigenvalue 1. In the iteration, we obtain ${\bf
y}_{k}:={\bf x}_{k}-{\bf x}_{*}={\mathcal{H}}_{Z}(A)^{k}{\bf u}+{\bf v}$. As a
result,

$\left\|{\bf x}_{*}-\frac{1}{N}\sum_{k=0}^{N-1}{\bf
x}_{k}\right\|^{2}=\frac{1}{N^{2}}\left\|\sum_{k=0}^{N-1}{\mathcal{H}}_{Z}(A)^{k}{\bf
u}\right\|^{2}+\|{\bf v}\|^{2}\leq\frac{\pi^{2}\eta_{Z}^{2}(A)}{N^{2}}\|{\bf
u}\|^{2}+\|{\bf v}\|^{2},$

where the first estimate is similar to the proof of Theorem 5.7. Thus ${\bf
v}=0$.

$(3\Rightarrow 2)$ Consider the same decomposition as above ${\bf x}_{0}-{\bf
x}_{*}={\bf u}+{\bf v}$. Note that $\|{\bf x}_{0}-({\bf x}_{*}+{\bf
v})\|^{2}=\|{\bf u}\|^{2}\leq\|{\bf u}\|^{2}+\|{\bf v}\|^{2}=\|{\bf
x}_{0}-{\bf x}_{*}\|$. If ${\mathcal{H}}_{Z}(A){\bf v}={\bf v}\Leftrightarrow
A{\bf v}=0$, then ${\bf x}_{*}+{\bf v}$ is also a solution of $A{\bf x}={\bf
b}$. By the property of ${\bf x}_{*}$ we have ${\bf v}=0$.

$(2\Rightarrow 3)$ Suppose ${\mathcal{H}}_{Z}(A){\bf v}={\bf v}$. We apply the
algorithm to ${\bf v}$, we will obtain a solution ${\bf x}_{1}$ such that
${\bf v}-{\bf x}_{1}={\bf u}_{1}+{\bf v}_{1}$. By statement 2, ${\bf
v}_{1}=0$. Similarly, if we start from $2{\bf v}$, then we will obtain ${\bf
x}_{2}$ with $2{\bf v}-{\bf x}_{2}={\bf u}_{2}$. So $A{\bf v}+A({\bf
u}_{1}-{\bf u}_{2})=0$. Since $A{\bf w}=0\Rightarrow{\mathcal{H}}_{Z}(A){\bf
w}={\bf w}$ for any ${\bf w}$, it follows that

$\displaystyle{\bf v}+{\bf u}_{1}-{\bf u}_{2}={\mathcal{H}}_{Z}(A){\bf
v}+{\mathcal{H}}_{Z}(A)({\bf u}_{1}-{\bf u}_{2})={\bf
v}+{\mathcal{H}}_{Z}(A)({\bf u}_{1}-{\bf u}_{2}).$

Hence we have ${\mathcal{H}}_{Z}(A)({\bf u}_{1}-{\bf u}_{2})={\bf u}_{1}-{\bf
u}_{2}$, i.e., ${\bf u}_{1}={\bf u}_{2}$. Therefore, we have $A{\bf v}=0$. ∎

The last condition (5.8) is verifiable before we conduct the algorithm. As a
result, we obtain the following theorem, whose proof is the same as that of
[21, Theorems 3.2, 3.9 and 3.15].

######  Theorem 5.9.

Assume the linear system $A{\bf x}={\bf b}$ is consistent. Fix a partition
$Z=\\{Z_{1},\ldots,Z_{p}\\}$ of rows of $A$. Let

${\mathcal{S}}_{0}:=\\{{\bf x}_{0},{\bf x}_{2k},{\bf
x}_{4k},\ldots\\},\quad{\mathcal{S}}_{1}:=\\{{\bf x}_{k},{\bf x}_{3k},{\bf
x}_{5k},\ldots\\}.$ (5.9)

Let ${\bf c}_{0},{\bf c}_{1}$ be the centers of the minimal sphere supporting
${\mathcal{S}}_{0},{\mathcal{S}}_{1}$ respectively. If (5.8) holds, then ${\bf
c}_{0}+{\bf c}_{1}=2{\bf x}_{*}$.

As a result, we re-obtain the following well-known result. It is obtained by
taking $p=1$ and $Z_{1}=A$ in (5.2).

######  Proposition 5.10.

Assume that $A\in\mathbb{R}^{m\times n}$ has full row rank and ${\bf
b}\in\mathbb{R}^{m}$. Let ${\bf x}_{0}\in\mathbb{R}^{n}$ and ${\bf
x}_{1}={\mathcal{H}}_{A}{\bf x}_{0}+2A(AA^{\tt T})^{+}{\bf b}$, then $A({\bf
x}_{0}+{\bf x}_{1})/2={\bf b}$. Consequently, for any ${\bf x}_{0}$, we have
that ${\bf x}_{0}-A(AA^{\tt T})^{+}A{\bf x}_{0}+A(AA^{\tt T})^{+}{\bf b}$ is a
solution of $A{\bf x}={\bf b}$.

###### Proof.

Since $A$ has full row rank, we have ${\mathcal{H}}_{A}=I_{n}-2A^{+}A$ and
$A{\mathcal{H}}_{A}=-A$. Now it is easy to check the claimed result. ∎

At the end of this part, we give some results on (5.8). The results below are
comparable to [21, Proposition 3.11 and Corollary 3.12].

######  Proposition 5.11.

If $A$ has full row rank, then (5.8) holds.

###### Proof.

By Lemma 5.4, ${\mathcal{H}}_{Z}(A)=I_{n}-2A^{\tt T}X^{-1}A$. If
${\mathcal{H}}_{Z}(A){\bf v}={\bf v}$, then $A^{\tt T}X^{-1}A{\bf v}=0$. Since
$A$ has full row rank, $X$ is invertible. As a result, $A{\bf v}=0$. ∎

Based on this, we below provide an effective approach to check (5.8) for
general matrices.

######  Fact 5.12 (Sylvester’s determinant theorem (also known as Weinstein-
Aronszajn identity)).

For any two matrices $A,B$ we have $\det(I+AB)=\det(I+BA)$.

######  Proposition 5.13.

Suppose $A=\begin{pmatrix}A_{1}\\\ A_{2}\end{pmatrix}\in\mathbb{R}^{m\times
n}$ and ${\rm Rank}(A)={\rm Rank}(A_{1})$. Assume that $A_{1}$ has full row
rank. Let $A_{2}=BA_{1}$ and let $X_{i}$ be the block-partition matrix of
$A_{i}$ for $i=1,2$. Then (5.8) holds if and only if $I-BX_{1}B^{\tt
T}X_{2}^{+}$ is invertible. In particular, if $X_{2}$ is invertible, then
(5.8) holds if and only if $X_{2}-BX_{1}B^{\tt T}$ is invertible.

###### Proof.

By Proposition 5.11, we have ${\mathcal{H}}({A_{1}}){\bf v}={\bf
v}\Leftrightarrow A_{1}{\bf v}=0$. Since $A_{2}$ lies in the row space of
$A_{1}$, we know that ${\mathcal{H}}({A_{1}}){\bf v}={\bf
v}\Rightarrow{\mathcal{H}}({A}){\bf v}={\bf v}$. Indeed,
${\mathcal{H}}({A_{1}}){\bf v}={\bf v}\Rightarrow A_{1}{\bf v}=0\Rightarrow
A_{2}{\bf v}=0\Rightarrow A{\bf v}=0\Rightarrow{\mathcal{H}}({A}){\bf v}={\bf
v}$. Thus (5.8) holds if and only if ${\mathcal{H}}({A}){\bf v}={\bf
v}\Rightarrow{\mathcal{H}}({A_{1}}){\bf v}={\bf v}$.

Since $A_{2}$ has full row rank, $X_{1}$ is invertible. Suppose ${\bf
v}={\mathcal{H}}({A}){\bf v}={\mathcal{H}}({A_{2}}){\mathcal{H}}({A_{1}}){\bf
v}$, then ${\mathcal{H}}({A_{2}}){\bf v}={\mathcal{H}}({A_{1}}){\bf v}$, i.e.,

$(I-2A_{2}^{\tt T}X_{2}^{+}A_{2}){\bf v}=(I-2A_{1}^{\tt T}X_{1}^{-1}A_{1}){\bf
v}.$

So we have

$A_{1}^{\tt T}X_{1}^{-1}A_{1}{\bf v}=A_{2}^{\tt T}X_{2}^{+}A_{2}{\bf
v}=A_{1}^{\tt T}B^{\tt T}X_{2}^{+}BA_{1}{\bf v},$ (5.10)

i.e.,

$A_{1}^{\tt T}(X_{1}^{-1}-B^{\tt T}X_{2}^{+}B)A_{1}{\bf v}=0.$

Since $A_{1}$ has full row rank, we have $0=(X_{1}^{-1}-B^{\tt
T}X_{2}^{+}B)A_{1}{\bf v}=X_{1}^{-1}(I-X_{1}B^{\tt T}X_{2}^{+}B)A_{1}{\bf v}$.
Thus a necessary and sufficient condition for $A_{1}{\bf v}=0$ is
$\det(I-X_{1}B^{\tt T}X_{2}^{+}B)\neq 0$. Indeed, if $\det(I-X_{1}B^{\tt
T}X_{2}^{+}B)=0$, then there is a ${\bf w}\neq 0$ such that $(I-X_{1}B^{\tt
T}X_{2}^{+}B){\bf w}=0$. Since $A_{1}$ has full row rank, there is a ${\bf
v}\neq 0$ such that $A_{1}{\bf v}={\bf w}\neq 0$. A contradiction.

By Sylvester’s determinant theorem (see Fact 5.12), we have

$\det(I-X_{1}B^{\tt T}X_{2}^{+}B)=\det(I-BX_{1}B^{\tt T}X_{2}^{+}).$

If $X_{2}$ is invertible, then it is further equal to $\det(X_{2}-BX_{1}B^{\tt
T})\det(X_{2}^{-1})$. ∎

We below provide more characterisations about the matrix $X_{2}-BX_{1}B^{\tt
T}$.

######  Proposition 5.14.

Making the same assumptions as Proposition 5.13. Let
$Z=\\{Z_{1},\cdots,Z_{p}\\}$ be a row partition of $A_{2}$.

  1. 1.

Assume that (5.8) holds, then $\sum_{i=1}^{p}{\rm Rank}(Z_{i})$ is even. In
particular, if $X_{2}$ is invertible and $d$ is the number of rows of $A_{2}$,
then $d$ is even, or equivalently, if $r={\rm Rank}(A)$, then $m-r$ is even.

  2. 2.

The matrix $X_{2}-BX_{1}B^{\tt T}$ is anti-symmetric, and its $(i,j)$-th block
is

$\begin{cases}\vspace{.1cm}\displaystyle\frac{1}{2}B_{i}(X_{1}^{{\tt
T}}-X_{1})B_{j}^{{\tt T}}&\text{if }i=j,\\\ B_{i}X_{1}^{{\tt T}}B_{j}^{{\tt
T}}&\text{if }i<j,\\\ -B_{i}X_{1}B_{j}^{{\tt T}}&\text{if }i>j,\end{cases}$

where $1\leq i,j\leq p$.

###### Proof.

We shall use the notation appeared in the proof of Proposition 5.13.

(1). Denote $D=\sum_{i=1}^{p}{\rm Rank}(Z_{i})$. Note that

$\det({\mathcal{H}}(A))=\det({\mathcal{H}}(A_{1}))\det({\mathcal{H}}(A_{2}))=(-1)^{D}\det({\mathcal{H}}(A_{1})).$

Since ${\mathcal{H}}(A),{\mathcal{H}}(A_{1})$ are unitary matrices, the
complex eigenvalues appear in conjugate forms. Following the proof of
Proposition 5.13, to ensure that ${\mathcal{H}}(A)$ and ${\mathcal{H}}(A_{1})$
have the same eigenspace of eigenvalue 1, $D$ must be even. When $X_{2}$ is
invertible, we have $d=D$.

(2). Note that $A_{2}=BA_{1}$, or more precisely in block form

$A_{2}=\begin{pmatrix}Z_{1}\\\ \vdots\\\
Z_{p}\end{pmatrix}=\begin{pmatrix}B_{1}A_{1}\\\ \vdots\\\
B_{p}A_{1}\end{pmatrix}.$

So $X_{2}+X_{2}^{\tt T}=2A_{2}A_{2}^{\tt T}=2BA_{1}A_{1}^{\tt T}B^{\tt
T}=B(X_{1}+X_{1}^{\tt T})B^{\tt T}$. The $(i,j)$-th block of $X_{2}$ is

$\begin{cases}\vspace{.1cm}\displaystyle\frac{1}{2}B_{i}(X_{1}+X_{1}^{{\tt
T}})B_{j}^{{\tt T}}&\text{if }i=j,\\\ B_{i}(X_{1}+X_{1}^{{\tt T}})B_{j}^{{\tt
T}}&\text{if }i<j.\end{cases}$

Hence, the $(i,j)$-th block of $X_{2}-BX_{1}B^{\tt T}$ is

$\begin{cases}\vspace{.1cm}\displaystyle\frac{1}{2}B_{i}(X_{1}^{{\tt
T}}-X_{1})B_{j}^{{\tt T}}&\text{if }i=j,\\\ B_{i}X_{1}^{{\tt T}}B_{j}^{{\tt
T}}&\text{if }i<j,\\\ -B_{i}X_{1}B_{j}^{{\tt T}}&\text{if }i>j.\end{cases}$

This means that $X_{2}-BX_{1}B^{\tt T}$ is anti-symmetric. ∎

In practice, we do not have $A_{1},A_{2}$ available in advance. Generally, we
start from a row partition $Z$, then run the algorithm. So what kind of
$Z=\\{Z_{1},\ldots,Z_{p}\\}$ should we choose? From Proposition 5.13, one
choice is to find an $A_{1}$ of full row rank from $Z$. But what if this is
not possible for a given $Z$. In the algorithm, it is required to compute
pseudoinverses of some matrices, so it is preferable to keep the matrices
small, i.e., each partition contains a few rows of $A$. With this intuition,
one natural choice is that each $Z_{i}$ has at most $\text{Rank}(A)$ rows.
From Proposition 5.14, it is necessary to make sure that $m-\text{Rank}(A)$ is
even. In summary, we provide the following simple criterion to run the
algorithm.

> The criterion: Let $Z=\\{Z_{1},\ldots,Z_{p}\\}$ be a row partition. To use
> this algorithm, it is preferable to make sure that each $Z_{i}$ has full-row
> rank and $m-\text{Rank}(A)$ is even.

Both conditions in the above criterion are easy to satisfy. As for the second
condition, we can introduce an extra trivial row if necessary. But this idea
is not always working from the proof of Proposition 5.13. There are some
nonsingularity conditions for it works. From (5.10), a general one is that we
randomly decompose $Z$ into two parts $A_{1}$ and $A_{2}$, then $A_{1}^{\tt
T}X_{1}^{-1}A_{1}-A_{2}^{\tt T}X_{2}^{+}A_{2}$ should be nonsingular.

If the linear system is consistent and $A$ has full column rank, then
${\mathcal{H}}_{A}=-I$. Thus in the iteration (5.2), if $i=1$ and $Z=A$, then
${\bf x}_{k+1}-{\bf x}_{*}=-({\bf x}_{k}-{\bf x}_{*})$, where ${\bf x}_{*}$ is
the solution of $A{\bf x}={\bf b}$. Namely, during the iteration, there are
only two points that are symmetric with respect to the solution. More
generally, given a partition $\\{Z_{1},\ldots,Z_{p}\\}$ of the rows of $A$, if
$Z_{i}$ for some $i$ has full column rank, then ${\mathcal{H}}_{Z_{i}}=-I$.
Even this, we can not delete it in the iteration since we need to make sure
that $m-{\rm Rank}(A)$ is even by Proposition 5.14.

###  5.2 Properties of $\eta_{Z}(A)$

In this subsection, we study more properties of $\eta_{Z}(A)$, which is
relevant to the convergence rate due to Theorem 5.7. We shall mainly consider
the case when each $Z_{i}$ is a row of $A$. Now the definition can be restated
as follows. It is a little different from [21, Definition 1.2], but
essentially the same thing.

######  Definition 5.15.

Let $A$ be an $m\times n$ matrix such that all rows are nonzero. Let
$\mathcal{R}_{i}$ be the reflection generated by the $i$-th row of $A$, and
let $\mathcal{R}_{A}=\mathcal{R}_{n}\cdots\mathcal{R}_{1}$. Denote the
eigenvalues of $\mathcal{R}_{A}$ be $e^{2i\theta_{1}},\ldots,e^{2i\theta_{m}}$
with $\theta_{j}\in[0,\pi)$ for all $j$. Define

$\eta(A):=\frac{1}{\min_{\ell:\theta_{\ell}\neq 0}|\sin(\theta_{\ell})|}.$
(5.11)

In [21, Definition 1.2], $\eta(A)$ is defined by
$\frac{1}{\min_{\ell:\theta_{\ell}\neq 0}|\theta_{\ell}|}.$ Here we feel
$\sin(\theta_{l})$ has more clear geometric meaning. A reason for this change
will be clear shortly. But theoretically, when $\theta_{\ell}$ is small, these
two quantities are close to each other.

The following results are easy to check through the definition of $\eta(A)$.

######  Proposition 5.16.

We have the following properties about $\eta(A):$

  * •

Let $Q$ be an orthogonal matrix, then $\mathcal{R}_{AQ}=Q^{\tt
T}\mathcal{R}_{A}Q$, so $\eta(A)$ is invariant under the right action of
orthogonal matrices.

  * •

Let $D$ be an invertible diagonal matrix, then $\eta(DA)=\eta(A)$, so
$\eta(A)$ is invariant under row scaling.

######  Lemma 5.17 (Lemma 2.2 of [21], also a special case of Lemma 5.4).

Let $W$ be the lower-triangular matrix satisfying $W+W^{\tt T}=2AA^{\tt T}$,
then $\mathcal{R}_{A}=I-2A^{\tt T}W^{-1}A$.

######  Theorem 5.18.

Let $A$ be an $m\times n$ matrix with no zero rows and let its $i$-th row as
$A_{i}$. Denote $c_{ij}=(A_{i}\cdot A_{j})/\|A_{i}\|\|A_{j}\|$. Then the
eigenvalues of $\mathcal{R}_{A}$ are determined by the following equation

$\det\begin{pmatrix}1+x&2c_{12}&\cdots&2c_{1m}\\\
2c_{12}x&1+x&\cdots&2c_{2m}\\\ \vdots&\vdots&\ddots&\vdots\\\
2c_{1m}x&2c_{2m}x&\cdots&1+x\end{pmatrix}=0.$ (5.12)

###### Proof.

Let $(x,{\bf u})$ be an eigenpair of $\mathcal{R}_{A}=I-2A^{\tt T}W^{-1}A$,
then we have $2A^{\tt T}W^{-1}A{\bf u}=(1-x){\bf u}$. Multiplying $A$ on both
sides gives rise to $(W+W^{\tt T})W^{-1}A{\bf u}=(1-x)A{\bf u}$, i.e.,
$(W^{\tt T}W^{-1}+xI)A{\bf u}=0$. If $A{\bf u}\neq 0$, then $\det(W^{\tt
T}W^{-1}+xI)=0$. Equivalently, $\det(W^{\tt T}+xW)=0$. If $A{\bf u}=0$, then
$x=1$. We also have $\det(W^{\tt T}+W)=\det(2AA^{\tt T})=0$. The multiplicity
of this eigenvalue $x=1$ is at least $m-\text{Rank}(A)$. Without considering
multiplicity, the eigenvalues are the same. In matrix form, we have (5.12). ∎

Theorem 5.18 was indeed already proved by Coxeter [5] in 1939. But this paper
is very old and hard to find. A statement of the result without a proof was
given in another paper of Coxeter [6].

When computing condition numbers, we usually consider the eigenvalues of
$AA^{\tt T}$. We now assume that all rows of $A$ are normalised to have norm
1. In matrix form, the characteristic polynomial of $AA^{\tt T}$ is

$\det\begin{pmatrix}1-x&c_{12}&\cdots&c_{1m}\\\ c_{12}&1-x&\cdots&c_{2m}\\\
\vdots&\vdots&\ddots&\vdots\\\ c_{1m}&c_{2m}&\cdots&1-x\end{pmatrix}=0.$
(5.13)

From (5.12) and (5.13) we can see that if $\eta(A)$ is large, i.e., when
$x=e^{2i\theta}\approx 1$ is a solution of (5.12), then
$(1-x)/2\approx\sin^{2}(\theta)$ is an approximate solution of (5.13). The
converse also holds. As a result, $\eta(A)\approx 1/\sigma_{\min}(A)$, where
$\sigma_{\min}$ refers to the minimal nonzero singular value of $A$.

Note that in the general case, the singular values are described by (5.13)
with diagonal entries replaced with
$1-\frac{x}{\|A_{1}\|^{2}},1-\frac{x}{\|A_{2}\|^{2}},\ldots,1-\frac{x}{\|A_{m}\|^{2}}$.
In comparison, $\eta(A)$ is always described by (5.12) because of Proposition
5.16.

In [21, Propositions 3.5 and 3.6], we proved that $\eta(A)=O(\kappa(A)^{2}\log
m)$ for all $A$ and $\eta(A)\approx\kappa(A)$ for random matrices $A$. Here
$\kappa(A)$ is the condition number of $A$. We below give a more quantitative
characterisation between $\eta$ and $\sigma_{\min}^{-1}$.

######  Theorem 5.19.

Let $L$ be the lower triangular matrix such that $L+L^{\tt T}=AA^{\tt T}-I$.
Then we have

$\frac{|\eta^{2}-\sigma^{-2}_{\min}(A)|^{2}}{\eta^{2}\sigma^{-2}_{\min}(A)\Big{(}1+\eta^{2}\Big{)}\Big{(}1+\sigma^{-2}_{\min}(A)\Big{)}}\leq
4\|L\|^{2}.$ (5.14)

###### Proof.

We will use Proposition 2.3. The singular values of $A$ are described by the
generalised eigenvalue problem $AA^{\tt T}{\bf u}=y{\bf u}$. In (5.12), let
$x=1-2y$, then we obtain $\det(AA^{\tt T}-y(I+2L))=0$. So $y$ can be described
by the generalised eigenvalue problem $AA^{\tt T}{\bf u}=y(I+2L){\bf u}$.
Using the notation of Proposition 2.3, we have $Z=(AA^{\tt T},I)$ and
$W=(AA^{\tt T},I+2L)$. By that proposition, we have

$S_{Z}(W)\leq d_{2}(Z,W)\leq\|(I+(AA^{\tt T})^{2})^{-1/2}\|\|2L\|\leq 2\|L\|.$

The first inequality is because $AA^{\tt T}$ is symmetric so that $T$ in
Proposition 2.3 is orthogonal. The second inequality follows from (2.1). The
last inequality is because $\|(I+(AA^{\tt
T})^{2})^{-1/2}\|=(1+\sigma_{\min}(A)^{4})^{-1/2}\leq 1$.

From the definition of $S_{Z}(W)$, we have that for any $y$ satisfying
$AA^{\tt T}{\bf u}=y(I+2L){\bf u}$,

$\frac{|y-\sigma^{2}_{\min}(A)|}{\sqrt{\Big{(}1+|y|^{2}\Big{)}\Big{(}1+\sigma^{2}_{\min}(A)\Big{)}}}\leq
S_{Z}(W)\leq 2\|L\|.$ (5.15)

Now assume that $\eta:=\eta(A)=1/\sin(\theta)$ for some $\theta>0$, i.e.,
there is a $x=e^{2i\theta}$ satisfying (5.12). Then

$y=\frac{1-x}{2}=\sin^{2}(\theta)-i\sin(\theta)\cos(\theta)=\eta^{-2}-i\eta^{-1}\sqrt{1-\eta^{-2}}.$

So
$|y-\sigma^{2}_{\min}(A)|\geq|\eta^{-2}-\sigma^{2}_{\min}(A)|=\frac{|\eta^{2}-\sigma^{-2}_{\min}(A)|}{\eta^{2}\sigma^{-2}_{\min}(A)}$
and $|y|^{2}=\eta^{-2}$. Substituting these into (5.15) yields

$4\|L\|^{2}\geq\frac{|\eta^{2}-\sigma^{-2}_{\min}(A)|^{2}}{\eta^{4}\sigma^{-4}_{\min}(A)\Big{(}1+\eta^{-2}\Big{)}\Big{(}1+\sigma^{2}_{\min}(A)\Big{)}}=\frac{|\eta^{2}-\sigma^{-2}_{\min}(A)|^{2}}{\eta^{2}\sigma^{-2}_{\min}(A)\Big{(}1+\eta^{2}\Big{)}\Big{(}1+\sigma^{-2}_{\min}(A)\Big{)}}.$

This completes the proof. ∎

From the above theorem, we can see that if $\|L\|$ is close to 0, i.e., when
$A$ is close to orthogonal, then $\eta\approx\sigma^{-1}_{\min}(A)$. We can
make this more precise as follows.

######  Definition 5.20 (restricted isometry property [3]).

Let $\delta\in(0,1)$ and $s\in\mathbb{N}$. A matrix $A\in\mathbb{C}^{m\times
n}$ is called to satisfy the $s$-restricted isometry property with restricted
isometry constant $\delta$ if for every submatrix $A_{s}\in\mathbb{C}^{m\times
s}$, we have

$(1-\delta)\|{\bf x}\|^{2}\leq\|A_{s}{\bf x}\|^{2}\leq(1+\delta)\|{\bf
x}\|^{2}$

for every ${\bf x}\in\mathbb{C}^{s}$. Namely,
$\|A_{s}^{\dagger}A_{s}-I_{s}\|_{2\rightarrow 2}\leq\delta$, or equivalently
all the eigenvalues of $A_{s}^{\dagger}A_{s}$ lie in the interval
$[1-\delta,1+\delta]$. Here $A_{s}^{\dagger}$ is the conjugate transpose of
$A_{s}$.

In randomised numerical linear algebra, Johnson-Lindenstrauss lemma shows that
there is a linear map $f:\mathbb{R}^{N}\rightarrow\mathbb{R}^{n}$ satisfying
that for any set $X$ of $m$ points in $\mathbb{R}^{N}$ and integer $n>8(\log
m)/\delta^{2}$, we have
$(1-\delta)\|u-v\|^{2}\leq\|f(u)-f(v)\|^{2}\leq(1+\delta)\|u-v\|^{2}$ for
every $u,v\in X$ with high probability, e.g., see [26]. The map $f$ plays an
important in randomised numerical linear algebra that can reduce a large scale
matrix problem to a small scale but “equivalent” matrix problem with a
tolerable error. Clearly, this map has the restricted isometry property in the
randomised sense.

######  Lemma 5.21 (Theorem 1 of [1]).

Let $M_{\ell}$ be the lower triangular part of $M\in\mathbb{C}^{m\times m}$,
then $\|M_{\ell}\|=O(\|M\|\log m)$ for every $m\geq 2$.

By Theorem 5.19 and Lemma 5.21 with $M=A^{\dagger}A-I$, we have the following
result.

######  Corollary 5.22.

Assume that $A\in\mathbb{R}^{m\times n}$ has $n$-restricted isometry property
with restricted isometry constant $\delta$, then

$\frac{|\eta^{2}-\sigma^{-2}_{\min}(A)|^{2}}{\eta^{2}\sigma^{-2}_{\min}(A)\Big{(}1+\eta^{2}\Big{)}\Big{(}1+\sigma^{-2}_{\min}(A)\Big{)}}\leq
O(\delta^{2}\log^{2}m).$ (5.16)

As shown in Proposition 5.16, $\eta(A)$ is invariant under the right action of
orthogonal matrices. Usually, it is not invariant under left action of
orthogonal matrices. We describe this change as follows.

######  Theorem 5.23.

Assume that $A\in\mathbb{R}^{m\times n}$. Let $P$ be an orthogonal matrix of
dimension $m$, then

$\frac{|\eta(PA)^{2}-\eta(A)^{2}|}{\eta(PA)\eta(A)\sqrt{(1+\eta(A)^{2})(1+\eta(PA)^{2})}}\leq
O\Big{(}\kappa(A)^{2}\|P^{\tt T}L_{PA}P-L_{A}\|^{1/m}\|A\|^{-2/m}\Big{)},$
(5.17)

where $L_{A}$, $L_{PA}$ are lower triangular matrices such that
$L_{A}+L_{A}^{\tt T}=AA^{\tt T}-I$ and $L_{PA}+L_{PA}^{\tt T}=PAA^{\tt
T}P^{\tt T}-I$ respectively.

###### Proof.

From the proof of Theorem 5.19, we know that the eigenvalues of
$\mathcal{R}_{A}$ are described by the generalised eigenvalue problem $AA^{\tt
T}{\bf u}=y(I+2L_{A}){\bf u}$. Similarly, the eigenvalues of
$\mathcal{R}_{PA}$ is described by the generalised eigenvalue problem
$PAA^{\tt T}P^{\tt T}{\bf u}=y(I+2L_{PA}){\bf u}$. Namely, they are described
by $AA^{\tt T}(P^{\tt T}{\bf u})=y(I+2P^{\tt T}L_{PA}P)(P^{\tt T}{\bf u})$. By
Proposition 2.1 with $Z=(AA^{\tt T},I+2L_{A})$ and $W=(AA^{\tt T},I+2P^{\tt
T}L_{PA}P)$, we have

$S_{Z}(W)\leq\frac{\sqrt{\|A\|^{4}+\|I+2L_{A}\|^{2}}}{\gamma(AA^{\tt
T},I+2L_{A})}\left(\frac{\sqrt{2}\|P^{\tt
T}L_{PA}P-L_{A}\|}{\sqrt{\|A\|^{4}+\|I+2L_{A}\|^{2}}}\right)^{1/m}.$

By Remark 2.2, $\gamma(AA^{\tt T},I+2L_{A})\geq\sigma_{\min}(AA^{\tt
T})=\sigma_{\min}(A)^{2}$. By Lemma 5.21, $\|L_{A}\|\leq\|A\|^{2}+1$. So

$\frac{\sqrt{\|A\|^{4}+\|I+2L_{A}\|^{2}}}{\gamma(AA^{\tt T},I+2L_{A})}\leq
O(\kappa_{A}^{2}).$

Similar to the argument used in the proof of Theorem 5.19, we have that

$\frac{|\eta(A)^{-2}-\eta(PA)^{-2}|}{\sqrt{1+\eta(A)^{-2}}\sqrt{1+\eta(PA)^{-2}}}\leq
S_{Z}(W).$

Simplifying this and putting the above all together lead to the claimed
result. ∎

As singular values do not change after orthogonal operations. So from Theorem
5.19 and the triangle inequality, we have another description of the change
described as follows:

$\frac{\Big{|}\eta(PA)^{2}-\eta(A)^{2}\Big{|}^{2}}{\max\Big{(}\eta(A)^{2},\eta(PA)^{2}\Big{)}\sigma^{-2}_{\min}(A)\Big{(}1+\max\Big{(}\eta(A)^{2},\eta(PA)^{2}\Big{)}\Big{)}\Big{(}1+\sigma^{-2}_{\min}(A)\Big{)}}\leq
8(\|L_{A}\|^{2}+\|L_{PA}\|^{2}).$

Similar to Corollary 5.22, for matrices with restricted isometry property, the
above bound on the right hand side can be small, and so $\eta(A)$ is close to
invariant under left actions of orthogonal matrices.

Proposition 5.16 also holds for $\eta_{Z}(A)$. Regarding Theorem 5.18, we
similarly have the following characterisation, whose proof is similar due to
Lemma 5.4.

######  Theorem 5.24.

Let $A$ be an $m\times n$ matrix with row partition
$Z=\\{Z_{1},\ldots,Z_{p}\\}$ such that all $Z_{i}$ has full row rank. Denote
$C_{ij}=(Z_{i}Z_{i}^{\tt T})^{-1/2}(Z_{i}Z_{j}^{\tt T})(Z_{j}Z_{j}^{\tt
T})^{-1/2}$. Then the eigenvalues of ${\mathcal{H}}_{Z}(A)$ are determined by
the following equation

$\det\begin{pmatrix}\vspace{.1cm}(1+x)I&2C_{12}&\cdots&2C_{1p}\\\
\vspace{.1cm}2C_{21}x&(1+x)I&\cdots&2C_{2p}\\\
\vspace{.1cm}\vdots&\vdots&\ddots&\vdots\\\
2C_{p1}x&2C_{p2}x&\cdots&(1+x)I\end{pmatrix}=0.$ (5.18)

Let

$B=\begin{pmatrix}(Z_{1}Z_{1}^{\tt T})^{-1/2}Z_{1}\\\ (Z_{2}Z_{2}^{\tt
T})^{-1/2}Z_{2}\\\ \vdots\\\ (Z_{p}Z_{p}^{\tt T})^{-1/2}Z_{p}\\\
\end{pmatrix},$

then

$BB^{\tt T}-xI=\det\begin{pmatrix}\vspace{.1cm}(1-x)I&C_{12}&\cdots&C_{1p}\\\
\vspace{.1cm}C_{21}&(1-x)I&\cdots&C_{2p}\\\
\vspace{.1cm}\vdots&\vdots&\ddots&\vdots\\\
C_{p1}&C_{p2}&\cdots&(1-x)I\end{pmatrix}.$

From the above analysis, it is not hard to imagine that $\eta_{Z}(A)$ relates
to the minimal singular value of $B$. We will not study this in more detail
here.

###  5.3 Connections of the solution of least squares

In [21], it was shown that the deterministic Kaczmarz method of the scheme
(3.1), where $i_{k}$ is chosen cyclically, cannot be used to solve least
squares. The iterative process indeed solves $A^{\tt T}W^{-1}A{\bf x}=A^{\tt
T}W^{-1}{\bf b}$, where $W$ is the lower triangular matrix such that $W+W^{\tt
T}=AA^{\tt T}-I$. In this section, we analyse in more detail why this is the
case and present a method to fix this. Some notation were introduced in [21,
Section 4].

######  Proposition 5.25.

Let $A\in\mathbb{R}^{m\times n}$ and

$F:=2\sum_{j=1}^{m}\frac{1}{\|A_{j}\|^{2}}\mathcal{R}_{m}\cdots\mathcal{R}_{j+1}A_{j}{\bf
e}_{j}^{\tt T},$

then $F=2A^{\tt T}W^{-1}$, where $W$ is the lower triangular matrix satisfying
$W+W^{\tt T}=2AA^{\tt T}$. Thus $FA=I-\mathcal{R}_{A}$.

###### Proof.

It suffices to prove $FW=2A^{\tt T}$. By definition,

$\displaystyle FW$ $\displaystyle=$ $\displaystyle
2\sum_{j=1}^{m}\frac{1}{\|A_{j}\|^{2}}\mathcal{R}_{m}\cdots\mathcal{R}_{j+1}A_{j}\left(\sum_{k=1}^{j-1}2A_{j}^{\tt
T}A_{k}{\bf e}_{k}^{\tt T}+\|A_{j}\|^{2}{\bf e}_{j}^{\tt T}\right)$
$\displaystyle=$ $\displaystyle
2\sum_{k=1}^{m-1}\left(\sum_{j=k+1}^{m}\frac{2}{\|A_{j}\|^{2}}\mathcal{R}_{m}\cdots\mathcal{R}_{j+1}A_{j}A_{j}^{\tt
T}\right)A_{k}{\bf e}_{k}^{\tt
T}+2\sum_{j=1}^{m}\mathcal{R}_{m}\cdots\mathcal{R}_{j+1}A_{j}{\bf e}_{j}^{\tt
T}$ $\displaystyle=$ $\displaystyle
2\sum_{k=1}^{m-1}\left(\sum_{j=k+1}^{m}\mathcal{R}_{m}\cdots\mathcal{R}_{j+1}(I_{n}-\mathcal{R}_{j})\right)A_{k}{\bf
e}_{k}^{\tt T}+2\sum_{j=1}^{m}\mathcal{R}_{m}\cdots\mathcal{R}_{j+1}A_{j}{\bf
e}_{j}^{\tt T}$ $\displaystyle=$ $\displaystyle
2\sum_{k=1}^{m-1}\left(I_{n}-\mathcal{R}_{m}\cdots\mathcal{R}_{k+1}\right)A_{k}{\bf
e}_{k}^{\tt T}+2\sum_{j=1}^{m}\mathcal{R}_{m}\cdots\mathcal{R}_{j+1}A_{j}{\bf
e}_{j}^{\tt T}$ $\displaystyle=$ $\displaystyle 2A^{\tt T}.$

This completes the proof. ∎

For any $i\in\\{0,1,\ldots,m-1\\}$, define $i$-th shifted product of
reflections as

$\mathcal{R}_{A}^{(i)}:=\mathcal{R}_{i}\cdots\mathcal{R}_{1}\mathcal{R}_{m}\mathcal{R}_{m-1}\cdots\mathcal{R}_{i+1},$
(5.19)

where $\mathcal{R}_{i}\cdots\mathcal{R}_{1}$ is viewed as the identity matrix
when $i=0$. So $\mathcal{R}_{A}=\mathcal{R}_{A}^{(0)}$. we also denote

$F_{A}^{(i)}=\sum_{j=1}^{m}\frac{2}{\|A_{j+i}\|^{2}}\mathcal{R}_{m+i}\mathcal{R}_{m+i-1}\cdots\mathcal{R}_{j+i+1}A_{j+i}{\bf
e}_{j}^{\tt T},$

and the shift matrix as

$Q=\begin{pmatrix}0&1&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&1\\\ 1&0&\cdots&0\\\ \end{pmatrix}_{m\times m}.$

For convenience, we now assume that $m-n$ is even and ${\rm Rank}(A)=n$. For
the linear system $A{\bf x}={\bf b}$, we start from an ${\bf x}_{0}$, then
iterate:

${\bf x}_{k+1}={\bf x}_{k}+2\frac{b_{i_{k}}-A_{i_{k}}^{{\tt T}}{\bf
x}_{k}}{\|A_{i_{k}}\|^{2}}A_{i_{k}}=\mathcal{R}_{i_{k}}{\bf
x}_{k}+2\frac{b_{i_{k}}}{\|A_{i_{k}}\|^{2}}A_{i_{k}},$ (5.20)

where $i_{k}=(k\mod m)+1$. We decompose $\\{{\bf x}_{0},{\bf x}_{1},{\bf
x}_{2},\cdots\\}$ into $m$ parts:

${\mathcal{S}}_{i}=\\{{\bf x}_{i+km}:k=0,1,2,\ldots,\\},\quad
i\in\\{0,1,\ldots,m-1\\}.$

Then each ${\mathcal{S}}_{i}$ lies on a high dimensional sphere, and according
to [21, Theorem 4.2] its centre is the solution of the linear system

$F_{A}^{(i)}Q^{i}A{\bf x}=F_{A}^{(i)}Q^{i}{\bf b}.$ (5.21)

######  Proposition 5.26.

Assume that $m-n$ is even and ${\rm Rank}(A)=n$. Let ${\bf x}_{*}^{(i)}$ be
the solution of $F_{A}^{(i)}Q^{i}A{\bf x}=F_{A}^{(i)}Q^{i}{\bf b}$. Then for
any $i\geq 1$

${\bf x}_{*}^{(i)}=\mathcal{R}_{i}{\bf
x}_{*}^{(i-1)}+\frac{2b_{i}}{\|A_{i}\|^{2}}A_{i}.$ (5.22)

###### Proof.

We first consider the case $i=1$. Now we have

$F^{(1)}=\sum_{j=1}^{m-1}\frac{2}{\|A_{j+1}\|^{2}}\mathcal{R}_{1}\mathcal{R}_{m}\cdots\mathcal{R}_{j+2}A_{j+1}{\bf
e}_{j}^{\tt T}+\frac{2}{\|A_{1}\|^{2}}A_{1}{\bf e}_{m}^{\tt T}.$

We also have
$F^{(1)}_{A}{QA}=I-\mathcal{R}_{QA}=\mathcal{R}_{1}(I-\mathcal{R}_{A})\mathcal{R}_{1}^{\tt
T}$. The equation for $i=1$ is

$(I-\mathcal{R}_{QA}){\bf
x}_{*}^{(1)}=\sum_{j=1}^{m-1}\frac{2b_{j+1}}{\|A_{j+1}\|^{2}}\mathcal{R}_{1}\mathcal{R}_{m}\cdots\mathcal{R}_{j+2}A_{j+1}+\frac{2b_{1}}{\|A_{1}\|^{2}}A_{1}.$

Namely,

$\displaystyle(I-\mathcal{R}_{A})\mathcal{R}_{1}^{\tt T}{\bf x}_{*}^{(1)}$
$\displaystyle=$
$\displaystyle\sum_{j=1}^{m-1}\frac{2b_{j+1}}{\|A_{j+1}\|^{2}}\mathcal{R}_{m}\cdots\mathcal{R}_{j+2}A_{j+1}-\frac{2b_{1}}{\|A_{1}\|^{2}}A_{1}$
$\displaystyle=$ $\displaystyle F^{(0)}{\bf
b}-\frac{2b_{1}}{\|A_{1}\|^{2}}(I-\mathcal{R}_{A})A_{1}.$

So $(I-\mathcal{R}_{A})(\mathcal{R}_{1}^{\tt T}{\bf
x}_{*}^{(1)}+\frac{2b_{1}}{\|A_{1}\|^{2}}A_{1})=F^{(0)}{\bf b}$, i.e., ${\bf
x}_{*}^{(0)}=\mathcal{R}_{1}^{\tt T}{\bf
x}^{(1)}+\frac{2b_{1}}{\|A_{1}\|^{2}}A_{1}$, or equivalently ${\bf
x}_{*}^{(1)}=\mathcal{R}_{1}{\bf
x}_{*}^{(0)}+\frac{2b_{1}}{\|A_{1}\|^{2}}A_{1}$. In this step, we used the
facts that $m-n$ is even and ${\rm Rank}(A)=n$ due to [21, Proposition 4.5].
Similarly, we have ${\bf x}_{*}^{(i)}=\mathcal{R}_{i}{\bf
x}_{*}^{(i-1)}+\frac{2b_{i}}{\|A_{i}\|^{2}}A_{i}$ for any $i\geq 1$. ∎

The connection (5.22) between ${\bf x}_{*}^{(i)}$ and ${\bf x}_{*}^{(i-1)}$ is
similar to the iterative procedure (5.20). The above result also implies that
for any $i\geq 1$

$(I+\mathcal{R}_{i})({\bf x}_{*}^{(i)}-{\bf x}_{*}^{(i-1)})=0,$

where ${\bf x}_{*}^{(m)}:={\bf x}_{*}^{(0)}$. So ${\bf x}_{*}^{(i)}-{\bf
x}_{*}^{(i-1)}=\alpha_{i}|A_{i}\rangle$ for some $\alpha_{i}$, where
$|A_{i}\rangle=\frac{A_{i}}{\|A_{i}\|}$ is the Dirac notation. Therefore,

$\sum_{i=1}^{m}\alpha_{i}|A_{i}\rangle=0.$ (5.23)

By (5.22), we have

$\displaystyle\alpha_{i}=\frac{2b_{i}}{\|A_{i}\|}-\frac{2}{\|A_{i}\|}A_{i}^{\tt
T}{\bf x}_{*}^{(i-1)}.$

As a direct application of Proposition 5.26, we have the following relation.

######  Corollary 5.27.

Let ${\bf x}_{*}^{(i)}$ be defined as in Proposition 5.26, then we have

$\sum_{i=1}^{m}\|A_{i}\|^{2}{\bf
x}_{*}^{(i)}=\sum_{i=1}^{m}\|A_{i}\|^{2}\mathcal{R}_{i}{\bf
x}_{*}^{(i-1)}+2A^{\tt T}{\bf b}.$ (5.24)

In comparison, let ${\bf x}_{\rm LS}$ the solution of the least square problem
$\arg\min\|A{\bf x}-{\bf b}\|$, i.e., $A^{\tt T}A{\bf x}_{\rm LS}=A^{\tt
T}{\bf b}$, then

$\Big{(}\sum_{i=1}^{m}\|A_{i}\|^{2}\Big{)}{\bf x}_{\rm
LS}=\Big{(}\sum_{i=1}^{m}\|A_{i}\|^{2}\mathcal{R}_{i}\Big{)}{\bf x}_{\rm
LS}+2A^{\tt T}{\bf b}.$ (5.25)

The difference of the two equations (5.24), (5.25) explains why we can not use
$\frac{1}{m}\sum_{i=1}^{m}{\bf x}_{*}^{(i)}$ to approximate ${\bf x}_{\rm LS}$
directly, i.e., why the iterative process (3.1) with $i_{k}$ chosen cyclically
does not solve the least squares. But for consistent linear systems, we have
${\bf x}_{*}^{(1)}=\cdots={\bf x}_{*}^{(m)}$. In this case, the above two
equations coincide with each other. Note that ${\bf x}_{*}^{(i)}$ is invariant
under row scaling, while ${\bf x}_{\rm LS}$ does.

We below present more connections between ${\bf x}^{(i)}$ and ${\bf x}_{\rm
LS}$.

######  Proposition 5.28.

Assume $\|A_{i}\|=1$ for all $i$, then we have the following properties:

  1. 1.

$\sum_{i=1}^{m}|A_{i}\rangle\langle A_{i},{\bf x}_{*}^{(i)}-{\bf x}_{\rm
LS}\rangle=0.$

  2. 2.

${\bf x}_{\rm
LS}=(mI_{m}-\sum_{j=1}^{m}\mathcal{R}_{j})^{-1}\sum_{i=1}^{m}(I_{m}-\mathcal{R}_{i}){\bf
x}^{(i-1)}.$

  3. 3.

If we assume that ${\bf x}_{\rm LS}=\sum_{i=1}^{m}\lambda_{i}{\bf
x}_{*}^{(i)}$, then $\lambda_{1},\ldots,\lambda_{m}$ satisfy

$\sum_{i,j=1}^{m}\lambda_{i}\Big{(}\mathcal{R}_{j}{\bf
x}_{*}^{(i)}-\mathcal{R}_{i}{\bf
x}_{*}^{(i-1)}\Big{)}+2\sum_{i=1}^{m}b_{i}(1-m\lambda_{i})A_{i}=0.$

###### Proof.

From (5.24) and (5.25), when $\|A_{i}\|=1$ for all $i$, we have

$\sum_{i=1}^{m}(I-\mathcal{R}_{i}){\bf x}_{*}^{(i-1)}=2A^{\tt T}{\bf
b}=\sum_{i=1}^{m}(I-\mathcal{R}_{i}){\bf x}_{\rm LS}.$

So

$\sum_{i=1}^{m}|A_{i}\rangle\langle A_{i},{\bf x}_{*}^{(i-1)}-{\bf x}_{\rm
LS}\rangle=0.$ (5.26)

This is indeed (5.23), but usually $\alpha_{i}\neq\langle A_{i},{\bf
x}_{*}^{(i-1)}-{\bf x}_{\rm LS}\rangle$. By (5.23) and (5.26), we also have

$\sum_{i=1}^{m}|A_{i}\rangle\langle A_{i},{\bf x}_{*}^{(i)}-{\bf x}_{\rm
LS}\rangle=0.$ (5.27)

From (5.26), we also obtain

${\bf x}_{\rm
LS}=(mI-\sum_{j=1}^{m}\mathcal{R}_{j})^{-1}\sum_{i=1}^{m}(I-\mathcal{R}_{i}){\bf
x}_{*}^{(i-1)}.$

Now assume ${\bf x}_{\rm LS}=\sum_{i=1}^{m}\lambda_{i}{\bf x}_{*}^{(i)}$, then

$m\sum_{i=1}^{m}\lambda_{i}{\bf
x}_{*}^{(i)}=\Big{(}\sum_{j=1}^{m}\mathcal{R}_{j}\Big{)}\sum_{i=1}^{m}\lambda_{i}{\bf
x}_{*}^{(i)}+2\sum_{i=1}^{m}b_{i}A_{i}.$

Namely,

$\sum_{i=1}^{m}\Big{(}m-\sum_{j=1}^{m}\mathcal{R}_{j}\Big{)}\lambda_{i}{\bf
x}_{*}^{(i)}=2A^{\tt T}{\bf b}.$

By (5.22), the LHS is

$m\sum_{i=1}^{m}\lambda_{i}(\mathcal{R}_{i}{\bf x}_{*}^{(i-1)}+2b_{i}A_{i}).$

So we obtain

$\sum_{i,j=1}^{m}\lambda_{i}\Big{(}\mathcal{R}_{j}{\bf
x}_{*}^{(i)}-\mathcal{R}_{i}{\bf
x}_{*}^{(i-1)}\Big{)}+2\sum_{i=1}^{m}b_{i}(1-m\lambda_{i})A_{i}=0.$

This completes the proof. ∎

###  5.4 Some final comments

In matrix form, we know from (5.22) that

$\begin{pmatrix}{\bf x}_{*}^{(i)}\\\
1\end{pmatrix}=\begin{pmatrix}\mathcal{R}_{i}&\frac{2b_{i}}{\|A_{i}\|^{2}}A_{i}\\\
0&1\end{pmatrix}\begin{pmatrix}{\bf x}_{*}^{(i-1)}\\\ 1\end{pmatrix}.$

This is the same procedure as the Kaczmarz method. In particular, when $A$ is
invertible (or has full row rank), then ${\bf x}_{*}^{(i)}={\bf
x}_{*}^{(i-1)}=A^{-1}{\bf b}=:{\bf x}_{*}$, so $\begin{pmatrix}{\bf x}_{*}\\\
1\end{pmatrix}$ lies in the invariant space of

${\mathcal{C}}_{i}:=\begin{pmatrix}\mathcal{R}_{i}&\frac{2b_{i}}{\|A_{i}\|^{2}}A_{i}\\\
0&1\end{pmatrix},\quad i\in[m].$

From the above, for any $\omega_{1},\ldots,\omega_{m}>0$, we also have

$\sum_{i=1}^{m}\omega_{i}{\mathcal{C}}_{i}{\bf x}_{*}=\omega{\bf x}_{*},$
(5.28)

where $\omega=\sum_{i}\omega_{i}$. Recall that Cimmino’s algorithm [4, 2]
reads as follows: Arbitrary choose ${\bf c}_{0}\in\mathbb{R}^{n}$, update

${\bf c}_{k+1}=\sum_{i=1}^{m}\frac{\omega_{i}}{\omega}\left({\bf
c}_{k}+2\frac{b_{i}-A_{i}^{\tt T}{\bf c}_{k}}{\|A_{i}\|^{2}}\right).$

In our notation, if we set ${\bf x}_{k}=({\bf c}_{k},1)^{\tt T}$, then the
above formula becomes

${\bf x}_{k+1}=\sum_{i=1}^{m}\frac{\omega_{i}}{\omega}{\mathcal{C}}_{i}{\bf
x}_{k}.$

This is consistent with (5.28). Cimmino’s algorithm aims to find an invariant
vector of the operator
$\sum_{i=1}^{m}\frac{\omega_{i}}{\omega}{\mathcal{C}}_{i}$.

Another way to understand Kaczmarz algorithm is as follows. Let
$P(t)={\mathcal{C}}_{1}+{\mathcal{C}}_{2}t+\cdots+{\mathcal{C}}_{m}t^{m-1}-mt^{m}$.
We say that $\lambda$ is an eigenvalue of $P(t)$ if there exists a nonzero
vector ${\bf x}$ satisfying $P(\lambda){\bf x}=0$; ${\bf x}$ is then a
corresponding eigenvector. Consider the block companion matrix of $P(t)$:

$C_{P}=\begin{pmatrix}0&I\\\ &0&I\\\ &&\ddots&\ddots\\\ &&&0&I\\\
\frac{{\mathcal{C}}_{1}}{m}&\frac{{\mathcal{C}}_{2}}{m}&\cdots&\frac{{\mathcal{C}}_{m-1}}{m}&\frac{{\mathcal{C}}_{m}}{m}\end{pmatrix}.$

It is easy to check that the eigenpairs $(\lambda,\tilde{{\bf x}})$ of $C_{P}$
must have the form $\tilde{{\bf x}}=({\bf x},\lambda{\bf x},\lambda^{2}{\bf
x},\cdots,\lambda^{m-1}{\bf x})^{\tt T}$ where $P(\lambda){\bf x}=0$. So we
can recover ${\bf x}_{*}$ if $\lambda$ is a solution of
$1+t+\cdots+t^{m-1}-mt^{m}=0$. Particularly, when $\lambda=1$, we have
$\tilde{{\bf x}}_{*}=({\bf x}_{*},{\bf x}_{*},{\bf x}_{*},\cdots,{\bf
x}_{*})^{\tt T}$. If we start from an arbitrary chosen vector $\tilde{{\bf
x}}_{0}=({\bf x}_{0},{\bf x}_{0},{\bf x}_{0},\cdots,{\bf x}_{0})^{\tt T}$ and
repeat the iteration $\tilde{{\bf x}}_{0}\rightarrow C_{P}\tilde{{\bf
x}}_{0}$, then this process converges to $\tilde{{\bf x}}_{*}$. Here the last
entry of ${\bf x}_{0}$ should be set as 1.

##  6 Numerical illustration

In this section, we mainly demonstrate the randomised reflective Kaczmarz
algorithms for solving linear systems. Here we will not compare the efficiency
of them with previous Kaczmarz algorithms, as some of the comparisons have
already been made in [21]. Our main focus here is to demonstrate the
differences of the algorithms presented in this paper. The code, written by
Maple 2024, is available via this link.222code link:
https://drive.google.com/file/d/16Qxeh91V9Isxky6FnJ24uUOibkAVbQPd/view?usp=sharing.
All experiments were conducted on an iMac Pro with processor 3.2 GHz 8-Core
Intel Xeon W.

Figure 1 exhibits the difference between the reflective Kaczmarz algorithm
(3.1) for consistent systems and inconsistent systems. As shown in Figure
1(a), for consistent systems, the points generated by the procedure (3.1) lie
on a sphere exactly. This is also guaranteed theoretically. For inconsistent
systems, as shown in Figure 1(b), the points no longer lie on a sphere but
instead are distributed around the solution in a manner similar to spheres.

(a) Generate a series of points on a sphere around the solution

(b) Generate a series of points converges to the solution

Figure 1: Comparison of reflective Kaczmarz algorithm (3.1) for consistent and
inconsistent linear systems. The red point in the center is the optimal
solution. The blue points are generated by the procedure (3.1). Figure 2:
Comparison of Algorithm 1 and Algorithm 2. Figure 3: Comparison of Algorithm 2
for different choices of $q$ for random linear systems of size $300\times
100$. Figure 4: The runtime (in seconds) of Algorithm 2 for different choices
of $q\in\\{1,2,\ldots,100\\}$ for random linear systems of size $300\times
100$.  Figure 5: A demonstration of ${\bf x}_{*}^{(i)}$, which is denoted as
$s_{i}$ because of some reason of exhibition, and ${\bf x}_{\rm LS}$ in
dimension 2.

In Figure 2, we show the difference between Algorithm 1 and Algorithm 2. The
theoretical results (i.e., Theorems 4.4, 4.6, 4.9 and 4.11) indicate that the
difference is small. This is also verified by numerical tests on random linear
systems of size $300\times 100$. In the plot, the least squares error is given
by $\|{\bf x}_{*}-\frac{1}{N}\sum_{k=0}^{N-1}{\bf x}_{k}\|$. We stopped the
algorithm when this error is smaller than $0.01$. This explains why the two
lines in Figure 2 stop at different numbers of iterations. In comparison,
Algorithm 2 offers more flexibility in implementation, as there is no need to
fix the partition of $A$ in advance.

In Figure 3, we demonstrate the convergence rate of Algorithm 2 for different
choices of $q$. As shown in the figure, to obtain a desired least squares
error, the number of iterations reduces when $q$ increases. This result aligns
with our theoretical expectations, as stated in Theorem 4.6.

Figure 4 shows the overall running time of Algorithm 2 for different choices
of $q$. The minimum runtime is achieved when $q=10$. The runtime decreases
initially as $q$ increases, but there is a sudden increase when $q=13$. After
that the runtime shows a tendency to decrease. All tested examples of size
$300\times 100$ exhibit similar behavior, as shown in Figure 4. This
phenomenon is reasonable because, theoretically, the algorithm converges
faster when more rows are used. At the same time, the runtime of computing
pseudoinverse also increases. So there is a trade-off between the number of
rows used at each step and the runtime for computing the pseudoinverse.

We emphasize again that, in practice, to run Algorithm 1 and Algorithm 2, we
would run the algorithms for a while to obtain a series of points $\\{{\bf
x}_{0},{\bf x}_{1},\ldots,{\bf x}_{N-1}\\}$ and then restart the algorithm
with this average as new initial point. Figures 2, 3 and 4 were indeed
obtained in this way.

Figure 5 illustrates ${\bf x}_{*}^{(i)}$ and ${\bf x}_{\rm LS}$. Each ${\bf
x}_{*}^{(i)}$ is the center of a sphere and ${\bf x}_{*}^{(i)}$ is obtained
from ${\bf x}_{*}^{(i-1)}$ via a reflection and a translation by Proposition
5.26. In Corollary 5.27 and Proposition 5.28, we provided some connections
between ${\bf x}_{*}^{(i)}$ and ${\bf x}_{\rm LS}$. But from Figure 5, it is
hard to imagine how they are connected.

Finally, we make a remark based on our numerical tests. Although the
theoretical results on the convergence rate of block Kaczmarz methods
presented in this paper are similar to those of many previous studies, our
numerical experiments reveal that the restarting strategy causes (block)
reflective Kaczmarz methods to converge much faster than expected. This might
be an advantage of this algorithm over others.

## Acknowledgements

The research is supported by the National Key Research Project of China under
Grant No. 2023YFA1009403.

## References

  * [1] J. R. Angelos, C. C. Cowen, and S. K. Narayan, Triangular truncation and finding the norm of a Hadamard multiplier, Linear Algebra and its Applications, 170 (1992), pp. 117–135.
  * [2] M. Benzi, Gianfranco Cimmino’s contributions to numerical mathematics, Atti del Seminario di Analisi Matematica, Dipartimento di Matematica dell’Universita di Bologna. Volume Speciale: Ciclo di Conferenze in Ricordo di Gianfranco Cimmino, (2004), pp. 87–109.
  * [3] E. J. Candes and T. Tao, Decoding by linear programming, IEEE transactions on Information Theory, 51 (2005), pp. 4203–4215.
  * [4] G. Cimmino, Calcolo approssimato per le soluzioni dei sistemi di equazioni lineari, La Ricerca Scientifica, II (1938), pp. 326–333.
  * [5] H. M. Coxeter, Lösung der aufgabe 245, Jahresbericht der Deutschen Mathematiker Vereinigung, 49 (1939), pp. 4–6.
  * [6] , The product of three reflections, Quarterly of Applied Mathematics, 5 (1947), pp. 217–222.
  * [7] T. Elfving, Block-iterative methods for consistent and inconsistent linear equations, Numerische Mathematik, 35 (1980), pp. 1–12.
  * [8] L. Elsner and J.-g. Sun, Perturbation theorems for the generalized eigenvalue problem, Linear Algebra and its Applications, 48 (1982), pp. 341–357.
  * [9] R. M. Gower and P. Richtárik, Randomized iterative methods for linear systems, SIAM Journal on Matrix Analysis and Applications, 36 (2015), pp. 1660–1690.
  * [10] D. Han, Y. Su, and J. Xie, Randomized Douglas–Rachford Methods for Linear Systems: Improved Accuracy and Efficiency, SIAM Journal on Optimization, 34 (2024), pp. 1045–1070.
  * [11] S. Karczmarz, Angenaherte auflosung von systemen linearer glei-chungen, Bull. Int. Acad. Pol. Sic. Let., Cl. Sci. Math. Nat., (1937), pp. 355–357.
  * [12] J. Liu and S. Wright, An accelerated randomized Kaczmarz algorithm, Mathematics of Computation, 85 (2016), pp. 153–178.
  * [13] A. Ma, D. Needell, and A. Ramdas, Convergence properties of the randomized extended Gauss–Seidel and Kaczmarz methods, SIAM Journal on Matrix Analysis and Applications, 36 (2015), pp. 1590–1604.
  * [14] J. D. Moorman, T. K. Tu, D. Molitor, and D. Needell, Randomized Kaczmarz with averaging, BIT Numerical Mathematics, 61 (2021), pp. 337–359.
  * [15] I. Necoara, Faster randomized block Kaczmarz algorithms, SIAM Journal on Matrix Analysis and Applications, 40 (2019), pp. 1425–1452.
  * [16] D. Needell, Randomized Kaczmarz solver for noisy linear systems, BIT Numerical Mathematics, 50 (2010), pp. 395–403.
  * [17] D. Needell and J. A. Tropp, Paved with good intentions: analysis of a randomized block Kaczmarz method, Linear Algebra and its Applications, 441 (2014), pp. 199–221.
  * [18] D. Needell, R. Ward, and N. Srebro, Stochastic gradient descent, weighted sampling, and the randomized Kaczmarz algorithm, Advances in Neural Information Processing Systems, 27 (2014).
  * [19] D. Needell, R. Zhao, and A. Zouzias, Randomized block Kaczmarz method with projection for solving least squares, Linear Algebra and its Applications, 484 (2015), pp. 322–343.
  * [20] R. Schreiber and C. Van Loan, A storage-efficient WY representation for products of Householder transformations, SIAM Journal on Scientific and Statistical Computing, 10 (1989), pp. 53–57.
  * [21] C. Shao, A deterministic Kaczmarz algorithm for solving linear systems, SIAM Journal on Matrix Analysis and Applications, 44 (2023), pp. 212–239.
  * [22] S. Steinerberger, Surrounding the solution of a linear system of equations from all sides, Quarterly of Applied Mathematics, 79 (2021), pp. 419–429.
  * [23] G. W. Stewart, An Elsner-like perturbation theorem for generalized eigenvalues, Linear Algebra and its Applications, 390 (2004), pp. 1–5.
  * [24] G. W. Stewart and J.-g. Sun, Matrix Perturbation Theory, Academic Press, INC, 1990.
  * [25] T. Strohmer and R. Vershynin, A randomized Kaczmarz algorithm with exponential convergence, Journal of Fourier Analysis and Applications, 15 (2009), pp. 262–278.
  * [26] D. P. Woodruff et al., Sketching as a tool for numerical linear algebra, Foundations and Trends® in Theoretical Computer Science, 10 (2014), pp. 1–157.
  * [27] A. Zouzias and N. M. Freris, Randomized extended Kaczmarz for solving least squares, SIAM Journal on Matrix Analysis and Applications, 34 (2013), pp. 773–793.