# Inertial randomized Kaczmarz algorithms for solving coherent linear systems

Songnian Hea,b, Ziting Wanga and Qiao-Li Donga,b
Tianjin Key Laboratory for Advanced Signal Processing and College of Science,
Civil Aviation University of China, Tianjin 300300, China,  Email:
<EMAIL_ADDRESS><EMAIL_ADDRESS>author. Email:
<EMAIL_ADDRESS>
###### Abstract

In this paper, by regarding the two-subspace Kaczmarz method [20] as an
alternated inertial randomized Kaczmarz algorithm we present a new convergence
rate estimate which is shown to be better than that in [20] under a mild
condition. Furthermore, we accelerate the alternated inertial randomized
Kaczmarz algorithm and introduce a multi-step inertial randomized Kaczmarz
algorithm which is proved to have a faster convergence rate. Numerical
experiments support the theory results and illustrate that the multi-inertial
randomized Kaczmarz algorithm significantly outperform the two-subspace
Kaczmarz method in solving coherent linear systems.

Keywords: System of linear equations, randomized Kaczmarz algorithm, inertial
randomized Kaczmarz algorithm, inertial extrapolation, two-subspace Kaczmarz
method.

##  1  Introduction

Let $\mathbb{R}^{J}$ be the real vector space of $J$-dimensional column
vectors with usual inner product $\langle\cdot,\cdot\rangle$ and norm
$\|\cdot\|$, respectively. We consider the following system of linear
equations:

$Ax=b,$ (1)

where $b\in\mathbb{R}^{I}$, $A\in\mathbb{R}^{I\times J}$ and
$x\in\mathbb{R}^{J}$. Denote by $A_{i},\,i=1,2,\cdots,I$, the row vectors of
$A$, then the linear equations in the system (1) can be expressed as
$A_{i}x=b_{i}$, $i=1,2,\cdots,I$. As is known to all, associated with each
equation $A_{i}x=b_{i}$, there is a hyperplane $H_{i}$ defined by

$H_{i}=\\{x{\in\mathbb{R}^{J}}\,|\,A_{i}x=b_{i}\\},\,\,\,\,\,i=1,2,\cdots,I.$

The orthogonal projection operator $P_{i}$ from $\mathbb{R}^{J}$ onto the
hyperplane $H_{i}$ assigns to each $z\in\mathbb{R}^{J}$ the unique point in
$H_{i}$ closest to $z$, denoted by $P_{i}z,$ (see (7) below for its closed
expression).

Kaczmarz algorithm is one of the most popular methods to solve (1) (see, e.g.,
[15], [6, Chapter 8]), which can be expressed in terms of the operators
$\\{P_{i}\\}_{i=1}^{I}$ as follows:

$x^{k+1}=P_{i(k)}x^{k},$

where $i(k)=(k\mod I)+1$. Because Kaczmarz algorithm uses only single equation
at each iteration step, it is also called a row-action [5, 7] or sequential
method. If the system (1) is consistent, the sequence
$\\{x^{k}\\}_{k=0}^{\infty}$ generated by Kaczmarz algorithm converges to the
solution closest to the initial value $x^{0}$.

The convergence speed of Kaczmarz algorithm heavily depends on the ordering of
the equations. In selecting the equation ordering, the important thing is to
avoid particularly bad ordering, in which case the hyperplanes $H_{i}$ and
$H_{i+1}$ are nearly parallel [13], i.e., $A_{i}$ and $A_{i+1}$ is highly
correlated (its definition is given later). Furthermore, it is very difficult
to estimate the convergence rate of Kaczmarz algorithm.

To overcome this difficulty, Strohmer and Vershynin [24] proposed the
following randomized Kaczmarz algorithm.

Algorithm 1 Randomized Kaczmarz algorithm

Input: $A$, $b$ and $x^{0}$

$k\leftarrow k+1$ Select $i_{k}\in\\{1,2,\ldots,I\\}$ $\\{$Randomly select a
row of $A$ with probability $\frac{\|A_{i_{k}}\|^{2}}{\|A\|_{F}^{2}}\\}$ Set
$x^{k+1}=P_{i_{k}}x^{k}$ $\\{$Perform projection$\\}$

Strohmer and Vershynin proved that Algorithm 1 converges linearly in
expectation when the linear system (1) is consistent and $A$ is a full column
rank matrix with $I\geq J$. Later, Ma and Needell et al. [17] showed that
whether the matrix $A$ satisfies $I\geq J$ with full column rank, or $I<J$
with full row rank, Algorithm 1 converges linearly in expectation to the
minimum norm solution of the linear system (1). Furthermore, Gower and
Richtárik [9] proved that even if the full rank assumption imposed on the
matrix $A$ is removed, the convergence result of Algorithm 1 still holds.
Precisely, their convergence result is as follows.

######  Theorem 1.1.

[9] Assume that the linear system (1) is consistent. Let the initial guess
$x^{0}$ be selected in ${\rm
span}\\{A_{1}^{\top},A_{2}^{\top},\cdots,A_{I}^{\top}\\}$ arbitrarily. Then
the sequence $\\{x^{k}\\}_{k=0}^{\infty}$ generated by Algorithm 1 converges
linearly in expectation to the minimum norm solution $x^{\dagger}$ for the
linear system (1), i.e.,

$\mathbb{E}\left[\|x^{k}-x^{\dagger}\|^{2}\right]\leq\left(1-\frac{\lambda_{\rm
min}(A^{\top}A)}{\|A\|^{2}_{F}}\right)^{k}\|x^{0}-x^{\dagger}\|^{2},\,\,\,\,\forall
k\geq 1,$ (2)

where $\lambda_{\rm min}(A^{\top}A)$ denotes the minimum nonzero eigenvalue of
$A^{\top}A$.

Highly correlated linear systems are often encountered in engineering
practice, for example sampling in geophysics which is the reconstruction of a
bandlimited function from nonuniformly spaced since it can be physically
challenging to take uniform samples. When the sampling is expressed as a
system of linear equations, the sampling points form the rows of a matrix $A$
which are close together and therefore the corresponding rows will be highly
correlated.

The quantity

$\delta(A)=\min_{i\neq j;\,1\leq i,j\leq
I}\left|\left\langle\frac{A_{i}}{\|A_{i}\|},\frac{A_{j}}{\|A_{j}\|}\right\rangle\right|.$

is generally used to measure the coherence between the rows of the matrix $A$.
The larger $\delta(A)$ means the higher coherence between row vectors of $A$.
Hereafter, $\delta(A)$ is called the coherence constant of $A$, and the linear
system (1) is called coherent if $\delta(A)>0$. It is easy to see that when
$\delta(A)$ is large, the convergence speed of Kaczmarz algorithm (even
randomized Kaczmarz algorithm) is relatively slow, regardless of the order of
equations. Therefore, it is of great significance to present accelerated
versions of the randomized Kaczmarz algorithm which can effectively solve the
coherent linear systems.

In [20], Needell and Ward introduced so called two-subspace Kaczmarz method
based on the idea of the orthogonalization, where the matrix $A$ is
standardized which means that each of its rows has unit Euclidean norm before
the implementation of the method.

Algorithm 2 Two-subspace Kaczmarz method

Input: $A$, $b$ and $x^{0}$

$k\leftarrow k+1$ Select $r,s\in\\{1,2,\ldots,I\\}$ $\,\,\\{$Select two
distinct rows of $A$ uniformly at random$\\}$ Set $\mu_{k}\leftarrow\langle
a_{r},a_{s}\rangle$ $\,\,\,\\{\hbox{Compute correlation}\\}$ Set
$y^{k}\leftarrow x^{k-1}+(b_{s}-\langle x^{k-1},a_{s}\rangle)a_{s}$
$\\{\hbox{Perform intermediate projection}\\}$ Set
$v_{k}\leftarrow\frac{\displaystyle a_{r}-\mu_{k}a_{s}}{\displaystyle
1-\mu_{k}^{2}}$ $\\{$Compute vector orthogonal to $a_{s}$ in direction of
$a_{r}$$\\}$ Set $\beta_{k}\leftarrow\frac{\displaystyle
b_{r}-\mu_{k}b_{s}}{\displaystyle 1-\mu_{k}^{2}}$ $\\{$Compute corresponding
measurement$\\}$ $x^{k}\leftarrow y^{k}+(\beta_{k}-\langle
y^{k},v_{k}\rangle)v_{k}$ $\\{$Perform projection$\\}$

The convergence rate estimate was established when $A$ is a full rank
standardized matrix with $I>J$ in [20]. In fact the full rank assumption can
be removed as Theorem 1.1 and the convergence results is given as follows.

######  Theorem 1.2.

Assume that the linear system (1) is consistent. Let the initial guess $x^{0}$
be selected in ${\rm span}\\{A_{1}^{\top},A_{2}^{\top},\cdots,A_{I}^{\top}\\}$
arbitrarily. Let $\\{x^{k}\\}$ be the sequence generated by the two-subspace
Kaczmarz method. Then

$\mathbb{E}\left[\|x^{k}-x^{\dagger}\|^{2}\right]\leq\left[\left(1-\frac{\lambda_{\rm
min}(A^{\top}A)}{\|A\|_{F}^{2}}\right)^{2}-\frac{\lambda_{\rm
min}(A^{\top}A)D}{\|A\|_{F}^{2}}\right]^{k}\|x^{0}-x^{\dagger}\|^{2},$ (3)

where
$D=\min\\{\frac{\delta(A)^{2}(1-\delta(A))}{1+\delta(A)},\frac{\Delta(A)^{2}(1-\Delta(A))}{1+\Delta(A)}\\}$
and $\Delta(A)=\max_{i\neq j;\,1\leq i,j\leq
I}\left|\left\langle\frac{A_{i}}{\|A_{i}\|},\frac{A_{j}}{\|A_{j}\|}\right\rangle\right|.$

Comparing the estimates (2) and (3), it is easy to see that the two-subspace
Kaczmarz method improves the randomized Kaczmarz algorithm if any rows of $A$
are highly correlated.

The inertial extrapolation firstly introduced by Polyak [22] has been widely
used to accelerate the convergence speed of iterative algorithms (see, e.g.,
[8] and [11]) since the cost of each iteration stays basically unchanged. To
introduce the inertial extrapolation, we start by introducing the term “basic
algorithm”. Consider the algorithmic operator
${\bf\Psi}:\mathbb{R}^{J}\rightarrow\mathbb{R}^{I}$ working iteratively by

$x^{k+1}={\bf\Psi}(x^{k}),\,\,\,\,\,\,\forall k\geq 0,$ (4)

which is denoted as the basic algorithm. The inertial version of the basic
algorithm (4) is defined as follows:

$x^{k+1}={\bf\Psi}(x^{k}+\alpha_{k}(x^{k}-x^{k-1})),\,\,\,\,\,\,\forall k\geq
1,$ (5)

where $\alpha_{k}(x^{k}-x^{k-1})$ is referred to as the inertial term and
$\alpha_{k}$ as the inertial parameter. The main feature of inertial type
methods is that the next iterate is defined by making use of the previous two
iterates. The inertial parameter sequence $\\{\alpha_{k}\\}_{k=0}^{\infty}$ is
key to inertial type algorithms and a careful choice is known to accelerate
the convergence rate of the objective function from $O(1/k)$ to $O(1/k^{2})$
or $o(1/k^{2})$ for a large class of first-order algorithms of the convex
optimization (see [1, 4, 21] for details). Since its inception, inertial type
method has received a great deal of attention of many authors, who improved it
in various ways. Mu and Peng [18] firstly introduced alternated inertial
methods applying inertia every other iteration, which was extended to solve
the composite convex optimization problem [14] and variational inequalities
[23]. Recently multi-step inertial methods [16] were proposed based on the
idea of using previous more than two iterates to accelerate the convergence
rate.

The purpose of this paper is three folds. Firstly, we give a new convergence
rate estimate of the two-subspace Kaczmarz method by rewriting it as an
alternated inertial randomized Kaczmarz algorithm and show that it is better
than that in [20] under a mild condition. Secondly, we introduce a multi-step
inertial randomized Kaczmarz algorithm by combining the inertial extrapolation
and the randomized Kaczmarz algorithm and its convergence rate is proved to be
faster than the two-subspace Kaczmarz method. Finally, numerical experiments
support the theory results and illustrate that the proposed algorithm
outperform the two-subspace Kaczmarz method in solving coherent linear
systems.

The paper is structured as follows. In the next section, we introduce some
fundamental concepts and tools, which are needed in the subsequent sections.
In Section 3, we rewrite the two-subspace Kaczmarz method as the alternated
inertial randomized Kaczmarz algorithm and present a new convergence rate. In
Section 4, we propose the multi-step inertial randomized Kaczmarz algorithm as
a further extension of the alternated inertial randomized Kaczmarz algorithm.
Numerical experiments are presented in Section 5 to show the efficiency and
advantage of our proposed algorithm.

##  2 Preliminaries

In this section, some fundamental tools are listed that will be used in the
convergence analysis of the proposed methods.

We will use the following notations:

  1. $\bullet$

$\mathcal{I}=\\{1,2,\cdots,I\\}$ and
$\Lambda=\\{(j,i)\,|\,j,i\in\mathcal{I}\,{\rm with}\,j\neq i\\}$ denote two
index sets, respectively.

  2. $\bullet$

$\lambda_{\min}(A^{\top}A)$ denotes the minimum nonzero eigenvalue of
$A^{\top}A$.

  3. $\bullet$

$\|A\|_{F}$ denotes the Frobenius norm of $A$.

  4. $\bullet$

$\alpha_{F}(A)=\max_{i\in\mathcal{I}}\\{\|A\|^{2}_{F}-\|A_{i}\|^{2}\\}$ and
$\Upsilon=\sum_{(j,i)\in\Lambda}\|A_{j}\|^{2}\|A_{i}\|^{2}$ denote two
constants, respectively.

  5. $\bullet$

$x^{\dagger}$ denotes the minimum norm solution of the problem (1).

  6. $\bullet$

$\mathbb{E}\left[\cdot\right]$ and $\mathbb{E}\left[\cdot\,|\,\cdot\right]$
denote the full expectation and the conditional expectation of a random
variable, respectively.

For a given nonempty closed convex subset $C$ of $\mathbb{R}^{J}$, the
projection operator onto $C$ is the mapping from $\mathbb{R}^{J}$ to $C$ that
assigns to each $x\in\mathbb{R}^{J}$ the unique point in $C$, denoted
$P_{C}x$, closest to $x$ that is,

$P_{C}x=\arg\min_{z\in C}\|z-x\|^{2}.$

A subset $H$ of $\mathbb{R}^{J}$ is called a hyperplane if $H$ is of the form

$H:=\\{x\in\mathbb{R}^{J}\mid\langle x,v\rangle=d\\},$ (6)

where $v$ is a fixed nonzero vector in $\mathbb{R}^{J}$ and $d$ is a number in
$\mathbb{R}$. Obviously, a hyperplane is closed and convex, and we have the
following closed-form formula for projections onto a hyperplane $H$ (see [3,
Chapter 25] or [12]).

######  Lemma 2.1.

Let $H$ be a hyperplane of the form (6). Then, for any $u\in\mathbb{R}^{J}$,
we have

$P_{H}u=u-\frac{\langle u,v\rangle-d}{\|v\|^{2}}v.$ (7)

In what follows, we denote by $P_{i}$ the projection operator onto
$H_{i}=\\{x{\in\mathbb{R}^{J}}\,|\,A_{i}x=b_{i}\\},i=1,2,\cdots,I.$

The following property is a well-known characterization of projections.

######  Lemma 2.2.

([10, Section 3]) Let $C$ be a closed convex subset of $\mathbb{R}^{J}.$ Given
$x\in\mathbb{R}^{J}$ and $z\in C$. Then $z=P_{C}x$ if and only if

$\langle x-z,y-z\rangle\leq 0,\quad\forall\,y\in C.$ (8)

Moreover, if $C$ is a hyperplane, then $z=P_{C}x$ if and only if

$\langle x-z,y-z\rangle=0,\quad\forall\,y\in C.$ (9)

######  Lemma 2.3.

[9] Let $A\in\mathbb{R}^{I\times J}$ be a matrix with row vectors
$A_{1},A_{2},\cdots,A_{I}$. Then the inequality

$\|Ax\|^{2}\geq\lambda_{\min}(A^{\top}A)\|x\|^{2}$ (10)

holds for all $x\in{\rm
span}\\{A_{1}^{\top},A_{2}^{\top},\cdots,A_{I}^{\top}\\}$.__

Throughout the rest of this paper, all discussions are based on the following
two basic assumptions:

  1. (H1)

The system of linear equations (1) is consistent, that is, its solution set,
denoted by $S$, is nonempty.

  2. (H2)

Any two equations in (1) are different, that is, the same equation is not
allowed to appear twice in (1).

##  3 Alternated inertial randomized Kaczmarz algorithm

In this section, we will give a better convergence rate estimate of the two-
subspace Kaczmarz method than (3) by regarding it as an alternated inertial
randomized Kaczmarz algorithm.

To this end, we first introduce an alternated inertial randomized Kaczmarz
algorithm as follows.

Algorithm 3 Alternated inertial randomized Kaczmarz algorithm

Input: $A$ and $b$

Step 0: Take $x^{0}\in{\rm
span}\\{A_{1}^{\top},A_{2}^{\top},\cdots,A_{I}^{\top}\\}$ arbitrarily and set
$k=1$.

Step 1: Select an index $(j_{k},i_{k})\in\Lambda$ with probability
$p_{(j_{k},i_{k})}=\frac{\|A_{j_{k}}\|^{2}\|A_{i_{k}}\|^{2}}{\Upsilon}$.

Step 2: Set  $y^{k}=P_{j_{k}}x^{k-1},$ (11) Set
$\beta_{k}=\frac{(A_{i_{k}}y^{k}-b_{i_{k}})\langle
A_{j_{k}},A_{i_{k}}\rangle}{\|A_{j_{k}}\|^{2}\|A_{i_{k}}\|^{2}-\langle
A_{j_{k}},A_{i_{k}}\rangle^{2}},$ (12) Calculate
$x^{k}=P_{i_{k}}(y^{k}+\beta_{k}A_{j_{k}}^{\top}).$ (13)

Step 3: Set $k:=k+1$ and return to Step 1.

######  Remark 3.1.

We give three important characteristics of Algorithm 3.

(i) $\beta_{k}$ given by (12) is well-defined for all $k\geq 0$. Indeed, since
$j_{k}\neq i_{k}$ can be guaranteed by Algorithm 3, it is easy to see that
$\|A_{j_{k}}\|^{2}\|A_{i_{k}}\|^{2}-\langle A_{j_{k}},A_{i_{k}}\rangle^{2}\neq
0$ under the assumption (H2).

(ii) Algorithm 3 is essentially an alternated inertial algorithm. In fact,
from (11) and (7), we get

$y^{k}-x^{k-1}=-\frac{A_{j_{k}}x^{k-1}-b_{j_{k}}}{\|A_{j_{k}}\|^{2}}A_{j_{k}}^{\top}.$
(14)

Obviously, in the case where $y^{k}\neq x^{k-1}$, we have

$y^{k}+\beta_{k}A_{j_{k}}^{\top}=y^{k}+\beta^{\prime}_{k}(y^{k}-x^{k-1}),$
(15)

where

$\beta^{\prime}_{k}=-\frac{(A_{i_{k}}y^{k}-b_{i_{k}})\|A_{j_{k}}\|^{2}\langle
A_{j_{k}},A_{i_{k}}\rangle}{(A_{j_{k}}x^{k-1}-b_{j_{k}})\left[\|A_{j_{k}}\|^{2}\|A_{i_{k}}\|^{2}-\langle
A_{j_{k}},A_{i_{k}}\rangle^{2}\right]}.$ (16)

Consequently, (13) can be rewritten as the form

$x^{k}=P_{i_{k}}(y^{k}+\beta^{\prime}_{k}(y^{k}-x^{k-1})).$ (17)

The formula (17) means that there is an inertial extrapolation in $x^{k}$.
Since there is no inertial extrapolation in $y^{k}$, Algorithm 3 is an
alternated inertial algorithm.

(iii) When $A$ is standardized, it is easy to verify that Algorithm 3 equals
to the two-subspace Kaczmarz method. The reason that we prefer the form of
Algorithm 3 instead of the two-subspace Kaczmarz method in this section is
that the former is convenient for us to make full use of the properties of the
projection operator to analyze the convergence rate more precisely, which can
be seen from Lemma 3.1 and its proof below.

We now give the following fundamental result, which is crucial to estimate the
convergence rate of Algorithm 3.

######  Lemma 3.1.

Let $\\{x^{k}\\}_{k=0}^{\infty}$ be the sequence generated by Algorithm 3.
Then for any $x^{*}\in S$ and $k\geq 1$, the following statements hold:

  1. (i)

$\|y^{k}-x^{*}\|^{2}=\|x^{k-1}-x^{*}\|^{2}-\frac{\displaystyle|A_{j_{k}}x^{k-1}-b_{j_{k}}|^{2}}{\displaystyle\|A_{j_{k}}\|^{2}};$

  2. (ii)

$\|y^{k}-x^{k-1}\|=\frac{\displaystyle|A_{j_{k}}x^{k-1}-b_{j_{k}}|}{\displaystyle\|A_{j_{k}}\|}$;

  3. (iii)

$\left\langle A^{\top}_{j_{k}},y^{k}-x^{*}\right\rangle=0$;

  4. (iv)

The inertial parameter sequence $\\{\beta_{k}\\}_{k=1}^{\infty}$ in (12) is
optimal in the sense that $\beta_{k}$ makes $\|x^{k}-x^{*}\|^{2}$ obtain its
minimum value;

  5. (v)

$\|x^{k}-x^{*}\|^{2}=\|y^{k}-x^{*}\|^{2}-\frac{\displaystyle|A_{i_{k}}y^{k}-b_{i_{k}}|^{2}}{\displaystyle\|A_{i_{k}}\|^{2}}\frac{\displaystyle
1}{\displaystyle 1-\left\langle\frac{\displaystyle
A_{j_{k}}}{\displaystyle\|A_{j_{k}}\|},\frac{\displaystyle
A_{i_{k}}}{\displaystyle\|A_{i_{k}}\|}\right\rangle^{2}};$

  6. (vi)

$\|x^{k}-y^{k}\|^{2}=\frac{\displaystyle|A_{i_{k}}y^{k}-b_{i_{k}}|^{2}}{\displaystyle\|A_{i_{k}}\|^{2}}\frac{\displaystyle
1}{\displaystyle 1-\left\langle\frac{\displaystyle
A_{j_{k}}}{\displaystyle\|A_{j_{k}}\|},\frac{\displaystyle
A_{i_{k}}}{\displaystyle\|A_{i_{k}}\|}\right\rangle^{2}}.$

###### Proof.

(i) Noting the fact that $A_{j_{k}}x^{*}=b_{j_{k}}$, we get

$A_{j_{k}}(x^{k-1}-x^{*})=A_{j_{k}}x^{k-1}-A_{j_{k}}x^{*}=A_{j_{k}}x^{k-1}-b_{j_{k}},$

and consequently, this together with (11) and (7) implies that

$\displaystyle\|y^{k}-x^{*}\|^{2}$ (18) $\displaystyle=$
$\displaystyle\|P_{j_{k}}x^{k-1}-x^{*}\|^{2}$ $\displaystyle=$
$\displaystyle\|x^{k-1}-x^{*}-\frac{A_{j_{k}}x^{k-1}-b_{j_{k}}}{\|A_{j_{k}}\|^{2}}A_{j_{k}}^{\top}\|^{2}$
$\displaystyle=$ $\displaystyle\|x^{k-1}-x^{*}\|^{2}-2\left\langle
x^{k-1}-x^{*},\frac{A_{j_{k}}x^{k-1}-b_{j_{k}}}{\|A_{j_{k}}\|^{2}}A_{j_{k}}^{\top}\right\rangle+\frac{|A_{j_{k}}x^{k-1}-b_{j_{k}}|^{2}}{\|A_{j_{k}}\|^{2}}$
$\displaystyle=$
$\displaystyle\|x^{k-1}-x^{*}\|^{2}-\frac{|A_{j_{k}}x^{k-1}-b_{j_{k}}|^{2}}{\|A_{j_{k}}\|^{2}}.$

(ii) From (14), we have

$\|y^{k}-x^{k-1}\|=\frac{|A_{j_{k}}x^{k-1}-b_{j_{k}}|}{\|A_{j_{k}}\|}.$

(iii) Noting $A_{j_{k}}x^{*}=b_{j_{k}}$, this conclusion is a direct
consequence of (14).

(iv) Obviously, in order to make Algorithm 3 obtain the optimal convergence
rate, for the current $y^{k}$, $\beta_{k}$ should be chosen such that
$\|x^{k}-x^{*}\|^{2}$ reaches its minimum value. To show that the inertial
parameter sequence $\\{\beta_{k}\\}_{k=1}^{\infty}$ in (12) is optimal, it
suffices to verify that

$\beta_{k}=\arg\min\limits_{\beta\in\mathbb{R}}\|P_{i_{k}}(y^{k}+\beta
A_{j_{k}}^{\top})-x^{*}\|^{2}$

holds for all $k\geq 1$.

Similar to the derivation of (18), we get

$\displaystyle\|P_{i_{k}}(y^{k}+\beta A_{j_{k}}^{\top})-x^{*}\|^{2}$
$\displaystyle=$ $\displaystyle\|y^{k}-x^{*}+\beta
A_{j_{k}}^{\top}\|^{2}-\frac{|A_{i_{k}}(y^{k}+\beta
A_{j_{k}}^{\top})-b_{i_{k}}|^{2}}{\|A_{i_{k}}\|^{2}}.$

Using (iii), one has

$\displaystyle\|y^{k}-x^{*}+\beta A_{j_{k}}^{\top}\|^{2}=$
$\displaystyle\|y^{k}-x^{*}\|^{2}+\beta^{2}\|A_{j_{k}}\|^{2}+2\beta\langle
y^{k}-x^{*},A_{j_{k}}^{\top}\rangle$ $\displaystyle=$
$\displaystyle\|y^{k}-x^{*}\|^{2}+\beta^{2}\|A_{j_{k}}\|^{2}.$

On the other hand,

$\displaystyle|A_{i_{k}}(y^{k}+\beta A_{j_{k}}^{\top})-b_{i_{k}}|^{2}=$
$\displaystyle|A_{i_{k}}y^{k}-b_{i_{k}}|^{2}+\beta^{2}\langle
A_{j_{k}},A_{i_{k}}\rangle^{2}$ $\displaystyle+2\beta\langle
A_{j_{k}},A_{i_{k}}\rangle(A_{i_{k}}y^{k}-b_{i_{k}}).$

Combining above three equalities, we have

$\displaystyle\|P_{i_{k}}(y^{k}+\beta A_{j_{k}}^{\top})-x^{*}\|^{2}$ (19)
$\displaystyle=$
$\displaystyle\|y^{k}-x^{*}\|^{2}-\frac{|A_{i_{k}}y^{k}-b_{i_{k}}|^{2}}{\|A_{i_{k}}\|^{2}}$
$\displaystyle+\beta^{2}\left(\|A_{j_{k}}\|^{2}-\frac{\langle
A_{j_{k}},A_{i_{k}}\rangle^{2}}{\|A_{i_{k}}\|^{2}}\right)-2\beta\frac{\langle
A_{j_{k}},A_{i_{k}}\rangle(A_{i_{k}}y^{k}-b_{i_{k}})}{\|A_{i_{k}}\|^{2}}.$

Set

$\displaystyle h(\beta):=$
$\displaystyle\beta^{2}\left(\|A_{j_{k}}\|^{2}-\frac{\langle
A_{j_{k}},A_{i_{k}}\rangle^{2}}{\|A_{i_{k}}\|^{2}}\right)$ (20)
$\displaystyle-2\beta\frac{\langle
A_{j_{k}},A_{i_{k}}\rangle(A_{i_{k}}y^{k}-b_{i_{k}})}{\|A_{i_{k}}\|^{2}},\,\,\,\forall\beta\in\mathbb{R}.$

For the current iterate $y^{k}$, noting that

$\|y^{k}-x^{*}\|^{2}-\frac{|A_{i_{k}}y^{k}-b_{i_{k}}|^{2}}{\|A_{i_{k}}\|^{2}}$

is a constant, it is easy to see from (19) and (20) that,
$\|P_{i_{k}}(y^{k}+\beta A_{j_{k}}^{\top})-x^{*}\|^{2}$ reaching its minimum
value is equivalent to the function $h(\beta)$ reaching its minimum value.
Note that $h(\beta)$ is a quadratic function of $\beta$, therefore, minimizing
$h(\beta)$ yields that the $k$th optimal inertial parameter is given by (12).

(v) The conclusion can be obtained directly by setting $\beta=\beta_{k}$ in
(19).

(vi) Using (12) and (13), we have

$x^{k}=y^{k}+\frac{(A_{i_{k}}y^{k}-b_{i_{k}})\left[\langle
A_{j_{k}},A_{i_{k}}\rangle
A_{j_{k}}^{\top}-\|A_{j_{k}}\|^{2}A_{i_{k}}^{\top}\right]}{\|A_{j_{k}}\|^{2}\|A_{i_{k}}\|^{2}-\langle
A_{j_{k}},A_{i_{k}}\rangle^{2}}.$ (21)

The desired conclusion follows from (21). ∎

######  Remark 3.2.

(i) Lemma 3.1 (iv) provides a reference for us to find $\beta_{k}$, whose idea
is just the selection strategy of inertial parameters in some inertial type
algorithms.

(ii) The accelerating convergence effect of Algorithm 3 can be clearly shown
by (v) and (vi) of Lemma 3.1. In fact, the iterative sequence generated by
Algorithm 1 only holds

$\|x^{k}-x^{*}\|^{2}=\|y^{k}-x^{*}\|^{2}-\frac{|A_{i_{k}}y^{k}-b_{i_{k}}|^{2}}{\|A_{i_{k}}\|^{2}},$

and

$\|x^{k}-y^{k}\|^{2}=\frac{|A_{i_{k}}y^{k}-b_{i_{k}}|^{2}}{\|A_{i_{k}}\|^{2}}.$

While, due to the inertia extrapolation, the iterative sequence generated by
Algorithm 3 has smaller $\|x^{k}-x^{*}\|^{2}$ and bigger
$\|x^{k}-y^{k}\|^{2}$. Indeed, from (v) and (vi) of Lemma 3.1, we get

$\|x^{k}-x^{*}\|^{2}\leq\|y^{k}-x^{*}\|^{2}-\frac{|A_{i_{k}}y^{k}-b_{i_{k}}|^{2}}{\|A_{i_{k}}\|^{2}}\frac{1}{1-\delta(A)^{2}},$

and

$\|x^{k}-y^{k}\|^{2}\geq\frac{|A_{i_{k}}y^{k}-b_{i_{k}}|^{2}}{\|A_{i_{k}}\|^{2}}\frac{1}{1-\delta(A)^{2}}.$

Obviously, the larger $\delta(A)$, the better the acceleration effect of
Algorithm 3.

Finally we analyze the convergence rate of Algorithm 3, i.e, the two-subspace
Kaczmarz method as follows.

######  Theorem 3.1.

Choose $x^{0}\in{\rm span}\\{A_{1}^{\top},A_{2}^{\top},\cdots,A_{I}^{\top}\\}$
arbitrarily and let $\left\\{x^{k}\right\\}_{k=0}^{\infty}$ be a sequence
generated by Algorithm 3. Then $\left\\{x^{k}\right\\}_{k=0}^{\infty}$
converges linearly in expectation to the minimum norm solution $x^{\dagger}$.
Precisely, for all $k\geq 1$, there holds the convergence rate estimate:

$\mathbb{E}\left[\|x^{k}-x^{\dagger}\|^{2}\right]\leq\left(1-\frac{\lambda_{\min}(A^{\top}A)}{\alpha_{F}(A)(1-\delta^{2}(A))}\right)^{k}\left(1-\frac{\lambda_{\min}(A^{\top}A)}{\|A\|_{F}^{2}}\right)^{k}\|x^{0}-x^{\dagger}\|^{2}.$
(22)

###### Proof.

Without losing generality, we assume that the system of equations (1) has been
standardized, that is, $\|A_{i}\|^{2}=1$ holds for each $i\in\mathcal{I}$. In
this case, it is very easy to see that $\|A\|_{F}^{2}=I$, $\alpha_{F}(A)=I-1$,
and $\Upsilon=I(I-1)$ hold. Consequently, Step 1 of Algorithm 3 is to select
$(j_{k},i_{k})\in\Lambda$ uniformly at random, i.e., with probability
$\frac{1}{I(I-1)}$.

Note that $A_{i}x^{\dagger}=b_{i}$ holds for all $i\in\mathcal{I}$, we have
from (i) and (v) of Lemma 3.1 that

$\|y^{k}-x^{\dagger}\|^{2}=\|x^{k-1}-x^{\dagger}\|^{2}-\displaystyle|A_{j_{k}}(x^{k-1}-x^{\dagger})|^{2}$
(23)

and

$\|x^{k}-x^{\dagger}\|^{2}\leq\|y^{k}-x^{\dagger}\|^{2}-\frac{\displaystyle|A_{i_{k}}(y^{k}-x^{\dagger})|^{2}}{(1-\delta^{2}(A))}$
(24)

hold for all $k\geq 1$, respectively. By (24) and the definition of
conditional expectation, we have

$\displaystyle\mathbb{E}\left[\|x^{k}-x^{\dagger}\|^{2}\,|\,x^{k-1}\right]$
(25) $\displaystyle\leq$
$\displaystyle\frac{1}{I(I-1)}\sum_{(j_{k},i_{k})\in\Lambda}\left\\{\|y^{k}-x^{\dagger}\|^{2}-\frac{\displaystyle|A_{i_{k}}(y^{k}-x^{\dagger})|^{2}}{(1-\delta^{2}(A))}\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{I(I-1)}\sum_{j_{k}\in\mathcal{I}}\sum_{i_{k}\in\mathcal{I}-\\{j_{k}\\}}\left\\{\|y^{k}-x^{\dagger}\|^{2}-\frac{\displaystyle|A_{i_{k}}(y^{k}-x^{\dagger})|^{2}}{(1-\delta^{2}(A))}\right\\}$
$\displaystyle=$
$\displaystyle\frac{1}{I}\sum_{j_{k}\in\mathcal{I}}\|y^{k}-x^{\dagger}\|^{2}-\frac{1}{I(I-1)(1-\delta^{2}(A))}\sum_{j_{k}\in\mathcal{I}}\sum_{i_{k}\in\mathcal{I}-\\{j_{k}\\}}|A_{i_{k}}(y^{k}-x^{\dagger})|^{2}.$

Noting that $y^{k}-x^{\dagger}\in{\rm
span}\\{A_{1}^{\top},A_{2}^{\top},\cdots,A_{I}^{\top}\\}$ and
$A_{j_{k}}(y^{k}-x^{\dagger})=0$ due to $y^{k}=P_{j_{k}}x^{k-1}$, it follows
from Lemma 2.3 that

$\displaystyle\sum_{i_{k}\in\mathcal{I}-\\{j_{k}\\}}|A_{i_{k}}(y^{k}-x^{\dagger})|^{2}=$
$\displaystyle\sum_{i_{k}\in\mathcal{I}}|A_{i_{k}}(y^{k}-x^{\dagger})|^{2}$
(26) $\displaystyle=$
$\displaystyle\|A(y^{k}-x^{\dagger})\|^{2}\geq\lambda_{\min}(A^{\top}A)\|y^{k}-x^{\dagger}\|^{2}.$

Combining (25) and (26), we have

$\mathbb{E}\left[\|x^{k}-x^{\dagger}\|^{2}\,|\,x^{k-1}\right]\leq\left(1-\frac{\lambda_{\min}(A^{\top}A)}{(I-1)(1-\delta^{2}(A))}\right)\frac{1}{I}\sum_{j_{k}\in\mathcal{I}}\|y^{k}-x^{\dagger}\|^{2}.$
(27)

Noting $x^{k-1}-x^{\dagger}\in{\rm
span}\\{A_{1}^{\top},A_{2}^{\top},\cdots,A_{I}^{\top}\\}$, we have from (23)
and Lemma 2.3 that

$\displaystyle\frac{1}{I}\sum_{j_{k}\in\mathcal{I}}\|y^{k}-x^{\dagger}\|^{2}=$
$\displaystyle\|x^{k-1}-x^{\dagger}\|^{2}-\frac{1}{I}\sum_{j_{k}\in\mathcal{I}}|A_{j_{k}}(x^{k-1}-x^{\dagger})|^{2}$
(28) $\displaystyle=$
$\displaystyle\|x^{k-1}-x^{\dagger}\|^{2}-\frac{1}{I}\|A(x^{k-1}-x^{\dagger})\|^{2}$
$\displaystyle\leq$
$\displaystyle\left(1-\frac{\lambda_{\min}(A^{\top}A)}{I}\right)\|x^{k-1}-x^{\dagger}\|^{2}.$

Substituting (28) into (27), we have

$\mathbb{E}\left[\|x^{k}-x^{\dagger}\|^{2}\,|\,x^{k-1}\right]\leq\left(1-\frac{\lambda_{\min}(A^{\top}A)}{(I-1)(1-\delta^{2}(A))}\right)\left(1-\frac{\lambda_{\min}(A^{\top}A)}{I}\right)\|x^{k-1}-x^{\dagger}\|^{2}.$
(29)

Hence, taking the full expectation at both sides of (29), we obtain

$\mathbb{E}\left[\|x^{k}-x^{\dagger}\|^{2}\right]\leq\left(1-\frac{\lambda_{\min}(A^{\top}A)}{(I-1)(1-\delta^{2}(A))}\right)\left(1-\frac{\lambda_{\min}(A^{\top}A)}{I}\right)\mathbb{E}\left[\|x^{k-1}-x^{\dagger}\|^{2}\right].$
(30)

Noting $\|A\|_{F}^{2}=I$ and $\alpha_{F}(A)=I-1$, (22) follows by applying
mathematical induction to (30). This completes the proof. ∎

Although Algorithm 3 is equivalent to the two-subspace Kaczmarz method when
$A$ is standardized, the proof of Theorem 3.1 is completely different from
that of Theorem 1.2, and the former is simpler. Furthermore, the convergence
rate estimate (22) and (3) are different and obviously, (22) is better than
(3) in at least two aspects:

  1. (i)

(22) is only related to $\delta(A)$, while (3) is related to both $\delta(A)$
and $\Delta(A)$;

  2. (ii)

When $\delta(A)=0$, the acceleration effect of the two-subspace Kaczmarz
method can not be reflected by (3) due to $D=0$, but (22) clearly shows the
acceleration effect of Algorithm 3 since $\alpha_{F}(A)<\|A\|_{F}^{2}$.

Further, we can prove that (22) is better than (3) under a mild condition.

######  Theorem 3.2.

If $2\leq r\leq I\leq
1+\frac{9\left(1+\sqrt{1-\frac{8}{9r}}\right)^{2}r^{2}(r-1)}{4},$ then (22) is
better than (3), where $r$ is the rank of $A$.

###### Proof.

For convenience, we assume that the system of equations (1) has been
standardized, i.e., $\|A_{i}\|^{2}=1$ for all $i\in\mathcal{I}$. Also,
$\lambda_{\rm min}(A^{\top}A)$ and $\delta(A)$ are simply expressed as
$\lambda_{\rm min}$ and $\delta$, respectively. Hence (22) and (3) become

$\mathbb{E}\left[\|x^{k}-x^{\dagger}\|^{2}\right]\leq\left(1-\frac{\lambda_{\min}}{(I-1)(1-\delta^{2})}\right)^{k}\left(1-\frac{\lambda_{\min}}{I}\right)^{k}\|x^{0}-x^{\dagger}\|^{2}.$
(31)

and

$\mathbb{E}\left[\|x^{k}-x^{\dagger}\|^{2}\right]\leq\left[\left(1-\frac{\lambda_{\rm
min}}{I}\right)^{2}-\frac{\lambda_{\rm
min}D}{I}\right]^{k}\|x^{0}-x^{\dagger}\|^{2},$ (32)

respectively. Note that

$\left(1-\frac{\lambda_{\min}}{(I-1)(1-\delta^{2})}\right)\left(1-\frac{\lambda_{\min}}{I}\right)=\left(1-\frac{\lambda_{\rm
min}}{I}\right)^{2}-\frac{\lambda_{\rm min}}{I}\left(1-\frac{\lambda_{\rm
min}}{I}\right)\frac{1+(I-1)\delta^{2}}{(I-1)(1-\delta^{2})},$

to prove the desired conclusion, from (31) and (32), we need to show that

$\left(1-\frac{\lambda_{\rm
min}}{I}\right)\frac{1+(I-1)\delta^{2}}{(I-1)(1-\delta^{2})}\geq D$ (33)

holds for all $\delta\in[0,1].$ Since $r\lambda_{\rm min}\leq I$ and
$D\leq\frac{\delta^{2}(1-\delta)}{1+\delta}$, to show (33), it suffices to
verify

$\left(1-\frac{1}{r}\right)\frac{1+(I-1)\delta^{2}}{(I-1)(1-\delta^{2})}\geq\frac{\delta^{2}(1-\delta)}{1+\delta},\,\,\forall\delta\in[0,1],$

that is,

$\left(1-\frac{1}{r}\right)\left(\frac{1}{I-1}+\delta^{2}\right)\geq\delta^{2}(1-\delta)^{2},\,\,\forall\delta\in[0,1].$
(34)

Define a real function:

$f(\delta)=\left(1-\frac{1}{r}\right)\left(\frac{1}{I-1}+\delta^{2}\right)-\delta^{2}(1-\delta)^{2},\,\,\forall\delta\in[0,1].$
(35)

It is easy to verify that $f$ obtains its minimum value at
$\hat{\delta}=\frac{2}{3r(1+\sqrt{1-\frac{8}{9r}})}$ and this implies that
proving (34) is equivalent to proving $f(\hat{\delta})\geq 0$. To prove
$f(\hat{\delta})\geq 0$, from (35), it suffices to verify that

$\left(1-\frac{1}{r}\right)\left(\frac{1}{I-1}+\hat{\delta}^{2}\right)\geq\hat{\delta}^{2}$
(36)

holds. Indeed, it is very easy to verify that (36) is equivalent to the
condition $2\leq r\leq I\leq
1+\frac{9(1+\sqrt{1-\frac{8}{9r}})^{2}r^{2}(r-1)}{4}.$ This completes the
proof. ∎

######  Remark 3.3.

A great deal of matrices satisfy the condition $2\leq r\leq I\leq
1+\frac{9(1+\sqrt{1-\frac{8}{9r}})^{2}r^{2}(r-1)}{4}$, for example, the row
full rank matrix, and column full rank matrix with $I\leq
1+\frac{9(1+\sqrt{1-\frac{8}{9J}})^{2}J^{2}(J-1)}{4}$.

##  4 Multi-step inertial randomized Kaczmarz algorithm

The alternated inertial randomized Kaczmarz algorithm is computationally
efficient and simple to implement, but this method still could be accelerated.
In this section, we accelerate the alternated inertial randomized Kaczmarz
algorithm based on the idea in multi-step inertial algorithms and thus name it
as the multi-step inertial randomized Kaczmarz algorithm.

Algorithm 4 Multi-step inertial randomized Kaczmarz algorithm

Input: $A$ and $b$

Step 0: Take $x^{0}\in{\rm
span}\\{A_{1}^{\top},A_{2}^{\top},\cdots,A_{I}^{\top}\\}$ arbitrarily and set
$w^{0}=x^{0}$. Select an index $i_{0}\in\\{1,2,\cdots,I\\}$ with probability
$p_{i_{0}}=\frac{\|A_{i_{0}}\|^{2}}{\|A\|^{2}_{F}}$. Calculate
$x^{1}=P_{i_{0}}w^{0}.$ (37) Set $k=1$.

Step 1: Select an index $i_{k}\in\mathcal{I}-\\{i_{k-1}\\}$ with probability
$p_{i_{k}}=\frac{\|A_{i_{k}}\|^{2}}{\|A\|^{2}_{F}-\|A_{i_{k-1}}\|^{2}}$,

Step 2: Set  $\gamma_{k}=\frac{(A_{i_{k}}x^{k}-b_{i_{k}})\langle
A_{i_{k-1}},A_{i_{k}}\rangle}{\|A_{i_{k-1}}\|^{2}\|A_{i_{k}}\|^{2}-\langle
A_{i_{k-1}},A_{i_{k}}\rangle^{2}}.$ (38) Set
$w^{k}=x^{k}+\gamma_{k}A_{i_{k-1}}^{\top}.$ (39) Calculate
$x^{k+1}=P_{i_{k}}w^{k}.$ (40)

Step 3: Set $k:=k+1$ and return to Step 1.

######  Remark 4.1.

There are two illustrations to Algorithm 4.

(i) Algorithm 4 is essentially a multi-step inertial algorithm. For each
$k\geq 1$, noting $x^{k}=P_{i_{k-1}}w^{k-1}$, we get

$x^{k}-w^{k-1}=-\frac{A_{i_{k-1}}w^{k-1}-b_{i_{k-1}}}{\|A_{i_{k-1}}\|^{2}}A_{i_{k-1}}^{\top}.$
(41)

It is easy to see from (39) and (41) that, in the case $x^{k}\neq w^{k-1}$,
(39) can be rewritten as the form:

$w^{k}=x^{k}+\gamma_{k}^{\prime}(x^{k}-w^{k-1}),$ (42)

where

$\gamma_{k}^{\prime}=-\frac{(A_{i_{k}}x^{k}-b_{i_{k}})\langle
A_{i_{k-1}},A_{i_{k}}\rangle\|A_{i_{k-1}}\|^{2}}{(A_{i_{k-1}}w^{k-1}-b_{i_{k-1}})[\|A_{i_{k-1}}\|^{2}\|A_{i_{k}}\|^{2}-\langle
A_{i_{k-1}},A_{i_{k}}\rangle^{2}]}.$ (43)

Consequently, (40) is of the form:

$x^{k+1}=P_{i_{k}}(x^{k}+\gamma_{k}^{\prime}(x^{k}-w^{k-1})).$ (44)

It is easy to see from (42) and (44) that constructing $x^{k+1}$ involves
$x^{0},x^{1},\cdots,x^{k}$, hence Algorithm 4 is a multi-step inertial
algorithm (see, e.g. [8]).

(ii) Assume that $A$ is standardized. Set $\mu_{k}=\langle
A_{i_{k-1}},A_{i_{k}}\rangle$ and $r^{k}=A_{i_{k}}x^{k}-b_{i_{k}}$. Then
(38)-(40) of Algorithm 4 can be simply expressed as: for the current
$x^{k}\,(k\geq 0)$, update $x^{k+1}$ by the formula

$x^{k+1}=x^{k}+\frac{r^{k}\mu_{k}}{1-\mu_{k}^{2}}A^{\top}_{i_{k-1}}-\frac{r^{k}}{1-\mu_{k}^{2}}A^{\top}_{i_{k}}.$
(45)

(45) is very suitable for the implementation of Algorithm 4.

######  Lemma 4.1.

The inertial parameter sequence $\\{\gamma_{k}\\}_{k=0}^{\infty}$ in (38) is
optimal in the sense that, for each $k\geq 1$, $\gamma_{k}$ makes
$\|x^{k+1}-x^{*}\|^{2}$ obtain its minimum value for arbitrary $x^{*}\in S$.

###### Proof.

For each $k\geq 0$, it concludes from (37), (40), (7) and the fact
$A_{i_{k}}x^{*}=b_{i_{k}}$ that

$\displaystyle\|x^{k+1}-x^{*}\|^{2}$ (46) $\displaystyle=$
$\displaystyle\|P_{i_{k}}w^{k}-x^{*}\|^{2}$ $\displaystyle=$
$\displaystyle\left\|w^{k}-x^{*}-\frac{A_{i_{k}}w^{k}-b_{i_{k}}}{\|A_{i_{k}}\|^{2}}A_{i_{k}}^{\top}\right\|^{2}$
$\displaystyle=$ $\displaystyle\|w^{k}-x^{*}\|^{2}-2\left\langle
w^{k}-x^{*},\frac{A_{i_{k}}w^{k}-b_{i_{k}}}{\|A_{i_{k}}\|^{2}}A_{i_{k}}^{\top}\right\rangle+\frac{|A_{i_{k}}w^{k}-b_{i_{k}}|^{2}}{\|A_{i_{k}}\|^{2}}$
$\displaystyle=$
$\displaystyle\|w^{k}-x^{*}\|^{2}-\frac{|A_{i_{k}}w^{k}-b_{i_{k}}|^{2}}{\|A_{i_{k}}\|^{2}}.$

For any $k\geq 1$, noting $A_{i_{k-1}}x^{*}=b_{i_{k-1}}$, we have from (41)
that

$\langle x^{k}-x^{*},A_{i_{k-1}}^{\top}\rangle=0.$ (47)

Consequently, we get

$\|w^{k}-x^{*}\|^{2}=\|x^{k}-x^{*}+\gamma_{k}A_{i_{k-1}}^{\top}\|^{2}=\|x^{k}-x^{*}\|^{2}+\gamma_{k}^{2}\|A_{i_{k-1}}\|^{2},\,\,\forall
k\geq 1.$ (48)

By direct calculation, we get

$\displaystyle|A_{i_{k}}w^{k}-b_{i_{k}}|^{2}$ (49) $\displaystyle=$
$\displaystyle|A_{i_{k}}(x^{k}+\gamma_{k}A_{i_{k-1}}^{\top})-b_{i_{k}}|^{2}$
$\displaystyle=$
$\displaystyle|A_{i_{k}}x^{k}-b_{i_{k}}|^{2}+\gamma_{k}^{2}\langle
A_{i_{k-1}},A_{i_{k}}\rangle^{2}+2\gamma_{k}\langle
A_{i_{k-1}},A_{i_{k}}\rangle(A_{i_{k}}x^{k}-b_{i_{k}}),\,\,\,\forall k\geq 1.$

Combining (46), (48) and (49), we obtain

$\displaystyle\|x^{k+1}-x^{*}\|^{2}$ (50) $\displaystyle=$
$\displaystyle\|x^{k}-x^{*}\|^{2}-\frac{|A_{i_{k}}x^{k}-b_{i_{k}}|^{2}}{\|A_{i_{k}}\|^{2}}+\gamma_{k}^{2}\left(\|A_{i_{k-1}}\|^{2}-\frac{\langle
A_{i_{k-1}},A_{i_{k}}\rangle^{2}}{\|A_{i_{k}}\|^{2}}\right)$
$\displaystyle-2\gamma_{k}\frac{\langle
A_{i_{k-1}},A_{i_{k}}\rangle(A_{i_{k}}x^{k}-b_{i_{k}})}{\|A_{i_{k}}\|^{2}},\,\,\,\forall
k\geq 1.$

By an argument similar to (iv) of Lemma 3.1, it is easy to see that
$\gamma_{k}$ in (38) such that $\|x^{k+1}-x^{*}\|^{2}$ reaches its minimum
value. ∎

######  Lemma 4.2.

Let $\\{x^{k}\\}_{k=0}^{\infty}$ be the sequence generated by Algorithm 4.
Then for any $k\geq 1$ and $x^{*}\in S$, the following statements hold:

  1. (i)

$\|x^{k+1}-x^{*}\|^{2}=\|x^{k}-x^{*}\|^{2}-\frac{\displaystyle|A_{i_{k}}x^{k}-b_{i_{k}}|^{2}}{\displaystyle\|A_{i_{k}}\|^{2}}\frac{\displaystyle
1}{\displaystyle 1-\left\langle\frac{\displaystyle
A_{i_{k-1}}}{\displaystyle\|A_{i_{k-1}}\|},\frac{\displaystyle
A_{i_{k}}}{\displaystyle\|A_{i_{k}}\|}\right\rangle^{2}}.$ (51)

  2. (ii)

$\|x^{k+1}-x^{k}\|^{2}=\frac{\displaystyle|A_{i_{k}}x^{k}-b_{i_{k}}|^{2}}{\displaystyle\|A_{i_{k}}\|^{2}}\frac{\displaystyle
1}{\displaystyle 1-\left\langle\frac{\displaystyle
A_{i_{k-1}}}{\displaystyle\|A_{i_{k-1}}\|},\frac{\displaystyle
A_{i_{k}}}{\displaystyle\|A_{i_{k}}\|}\right\rangle^{2}}.$ (52)

###### Proof.

(i) Obviously, (51) follows from substituting (38) into (50).

(ii) From (39) and (40), we have

$x^{k+1}-x^{k}=\gamma_{k}\left\\{A_{i_{k-1}}^{\top}-\frac{\langle
A_{i_{k-1}},A_{i_{k}}\rangle}{\|A_{i_{k}}\|^{2}}A_{i_{k}}^{\top}\right\\}-\frac{A_{i_{k}}x^{k}-b_{i_{k}}}{\|A_{i_{k}}\|^{2}}A_{i_{k}}^{\top}.$

Noting that

$\left\langle A_{i_{k-1}}^{\top}-\frac{\langle
A_{i_{k-1}},A_{i_{k}}\rangle}{\|A_{i_{k}}\|^{2}}A_{i_{k}}^{\top},A_{i_{k}}^{\top}\right\rangle=0,$

we have

$\|x^{k+1}-x^{k}\|^{2}=\frac{|A_{i_{k}}x^{k}-b_{i_{k}}|^{2}}{\|A_{i_{k}}\|^{2}}+\gamma_{k}^{2}\left(\|A_{i_{k-1}}\|^{2}-\frac{\langle
A_{i_{k-1}},A_{i_{k}}\rangle^{2}}{\|A_{i_{k}}\|^{2}}\right)$ (53)

Hence, (52) yields by substituting (38) into (53). ∎

######  Theorem 4.1.

Choose $x^{0}\in{\rm span}\\{A_{1}^{\top},A_{2}^{\top},\cdots,A_{I}^{\top}\\}$
arbitrarily and let $\left\\{x^{k}\right\\}_{k=0}^{\infty}$ be a sequence
generated by Algorithm 4. Then $\left\\{x^{k}\right\\}_{k=0}^{\infty}$
converges linearly in expectation to the minimum norm solution $x^{\dagger}$.
Precisely, for all $k\geq 1$, there holds the following statement:

$\mathbb{E}\left[\|x^{k+1}-x^{\dagger}\|^{2}\right]\leq\left(1-\frac{\lambda_{\min}(A^{\top}A)}{\alpha_{F}(A)(1-\delta^{2}(A))}\right)^{k}\left(1-\frac{\lambda_{\min}(A^{\top}A)}{\|A\|_{F}^{2}}\right)\|x^{0}-x^{\dagger}\|^{2}.$
(54)

###### Proof.

Note that $w^{0}=x^{0}$, we have from (37), (46) and (7) that

$\|x^{1}-x^{\dagger}\|^{2}=\|x^{0}-x^{\dagger}\|^{2}-\frac{|A_{i_{0}}x^{0}-b_{i_{0}}|^{2}}{\|A_{i_{0}}\|^{2}}.$

Due to $x^{0}-x^{\dagger}\in{\rm
span}\\{A_{1}^{\top},A_{2}^{\top},\cdots,A_{I}^{\top}\\}$, this together with
Lemma 2.3 leads to

$\displaystyle\mathbb{E}\left[\|x^{1}-x^{\dagger}\|^{2}\right]$
$\displaystyle=\|x^{0}-x^{\dagger}\|^{2}-\sum_{i_{0}\in\mathcal{I}}\frac{\|A_{i_{0}}\|^{2}}{\|A\|_{F}^{2}}\frac{|A_{i_{0}}x^{0}-b_{i_{0}}|^{2}}{\|A_{i_{0}}\|^{2}}$
(55)
$\displaystyle=\|x^{0}-x^{\dagger}\|^{2}-\sum_{i_{0}\in\mathcal{I}}\frac{\|A_{i_{0}}\|^{2}}{\|A\|_{F}^{2}}\frac{|A_{i_{0}}(x^{0}-x^{\dagger})|^{2}}{\|A_{i_{0}}\|^{2}}$
$\displaystyle\leq\|x^{0}-x^{\dagger}\|^{2}-\frac{1}{\|A\|_{F}^{2}}\|A(x^{0}-x^{\dagger})\|^{2}$
$\displaystyle\leq\left(1-\frac{\lambda_{\min}(A^{\top}A)}{\|A\|_{F}^{2}}\right)\|x^{0}-x^{\dagger}\|^{2}.$

For any $k\geq 1$, from (51), we have

$\|x^{k+1}-x^{\dagger}\|^{2}\leq\|x^{k}-x^{\dagger}\|^{2}-\frac{\displaystyle|A_{i_{k}}x^{k}-b_{i_{k}}|^{2}}{\displaystyle\|A_{i_{k}}\|^{2}}\frac{\displaystyle
1}{\displaystyle 1-\delta^{2}(A)}.$ (56)

From Step 1 of Algorithm 4 and (56), we get

$\displaystyle\mathbb{E}\left[\|x^{k+1}-x^{\dagger}\|^{2}\,|\,x^{k}\right]$
(57) $\displaystyle\leq$
$\displaystyle\|x^{k}-x^{\dagger}\|^{2}-\frac{\displaystyle 1}{\displaystyle
1-\delta^{2}(A)}\sum_{i_{k}\in\mathcal{I}-\\{i_{k-1}\\}}\frac{\|A_{i_{k}}\|^{2}}{\|A\|^{2}_{F}-\|A_{i_{k-1}}\|^{2}}\frac{\displaystyle|A_{i_{k}}x^{k}-b_{i_{k}}|^{2}}{\displaystyle\|A_{i_{k}}\|^{2}}$
$\displaystyle\leq$
$\displaystyle\|x^{k}-x^{\dagger}\|^{2}-\frac{\displaystyle
1}{\displaystyle\alpha_{F}(A)(1-\delta^{2}(A))}\sum_{i_{k}\in\mathcal{I}-\\{i_{k-1}\\}}|A_{i_{k}}x^{k}-b_{i_{k}}|^{2}.$

On the other hand, from (47), it yields $A_{i_{k-1}}x^{k}-b_{i_{k-1}}=0$. This
together with Lemma 2.3 leads to

$\displaystyle\sum_{i_{k}\in\mathcal{I}-\\{i_{k-1}\\}}|A_{i_{k}}x^{k}-b_{i_{k}}|^{2}$
(58) $\displaystyle=$
$\displaystyle\sum_{i_{k}\in\mathcal{I}}|A_{i_{k}}x^{k}-b_{i_{k}}|^{2}=\sum_{i_{k}\in\mathcal{I}}|A_{i_{k}}(x^{k}-x^{\dagger})|^{2}\geq\lambda_{\min}(A^{\top}A)\|x^{k}-x^{\dagger}\|^{2}.$

Combining (57) and (58), it yields

$\mathbb{E}\left[\|x^{k+1}-x^{\dagger}\|^{2}\,|\,x^{k}\right]\leq\left(1-\frac{\lambda_{\min}(A^{\top}A)}{\alpha_{F}(A)(\displaystyle
1-\delta^{2}(A))}\right)\|x^{k}-x^{\dagger}\|^{2}.$ (59)

Hence, taking the full expectation at both sides of (60), we obtain

$\mathbb{E}\left[\|x^{k+1}-x^{\dagger}\|^{2}\right]\leq\left(1-\frac{\lambda_{\min}(A^{\top}A)}{\alpha_{F}(A)(\displaystyle
1-\delta^{2}(A))}\right)\mathbb{E}\left[\|x^{k}-x^{\dagger}\|^{2}\right].$
(60)

Using mathematical induction, (54) follows from (55) and (60). ∎

######  Remark 4.2.

Since each iteration of Algorithm 3 utilizes two rows of the matrix $A$, one
needs to equate a single iteration of Algorithm 3 with two iterations of
Algorithm 4 for fair comparison of the convergence rate. It is obvious that
the estimate (54) is better than the estimate (22) in this sense.

##  5 Numerical Results

In this section, we perform several experiments to compare the convergence
speed of the multi-step inertial randomized Kaczmarz algorithm (denoted as
MIRK) with that of the two-subspace Kaczmarz method (denoted as TSK).

We use “IT” and “CPU” to represent the arithmetical averages of the required
numbers of iteration steps required and the elapsed CPU times taken to run the
algorithms 50 times, respectively. The speed-up of MIRK against TSK is defined
by

$\hbox{speed-up}=\frac{\hbox{CPU of TSK}}{\hbox{CPU of MIRK}}.$

To test methods, we construct various types of $I\times J$ matrices $A$ and
set the entries of $A$ to be independent identically distributed uniform
random variables on some interval $[c,1]$. Changing the value of $c$ will
appropriately change the coherence of $A$.

In the implementation of the methods, we use the function $unifrnd$ to
randomly generate a solution vector $x^{*}\in\mathbb{R}^{J}$, and take the
right-hand side $b\in\mathbb{R}^{I}$ as $Ax^{*}$. This ensures that the
randomly generated matrices constitute a solvable system of equations $Ax=b$.
For each matrix construction, two algorithms are run with the initial vector
$x^{0}=0$ and the same fixed matrix, and terminated once the relative solution
error, defined by

$\hbox{RSE}=\frac{\|x^{k}-x^{\dagger}\|^{2}}{\|x^{\dagger}\|^{2}}$

at the current $x^{k}$, satisfies $\hbox{RSE}\leq 10^{-6}$. In addition, all
experiments are carried out using MATLAB (version R2018a) on a personal
computer with 2.30 GHz central processing unit (Intel(R) Core(TM) i5-6300HQ
CPU), 12.00GB memory, and Windows operating system (Windows 10).

Table 1:  IT and CPU of TSK and MIRK for $1000\times 3000$ matrices with different $c$.  $c$ | 0.9 | 0.5 | 0.1 | -0.4
---|---|---|---|---
TSK | CPU | 3.5000 | 3.3125 | 3.1719 | 2.2188
IT | $2.7362\times 10^{4}$ | $2.6849\times 10^{4}$ | $2.6548\times 10^{4}$ | $2.3490\times 10^{4}$
MIRK | CPU | 2.4844 | 2.4688 | 2.4063 | 2.2031
IT | $3.7174\times 10^{4}$ | $3.7282\times 10^{4}$ | $3.6742\times 10^{4}$ | $3.3572\times 10^{4}$
speed-up | 1.4088 | 1.3417 | 1.3182 | 1.0071
Table 2:  IT and CPU of TSK and MIRK for $2000\times 1000$ matrices with different $c$.  $c$ | 0.9 | 0.5 | 0.1 | -0.2
---|---|---|---|---
TSK | CPU | 3.4531 | 3.3594 | 2.9375 | 2.7031
IT | $5.0883\times 10^{4}$ | $4.7767\times 10^{4}$ | $4.5534\times 10^{4}$ | $6.6282\times 10^{4}$
MIRK | CPU | 2.7969 | 2.7344 | 2.6563 | 2.6250
IT | $6.8314\times 10^{4}$ | $6.7658\times 10^{4}$ | $6.7208\times 10^{4}$ | $7.9606\times 10^{4}$
speed-up | 1.2346 | 1.2286 | 1.1059 | 1.0298

(a) (b) (c) (b)

Figure 1: RSE versus CPU with $1000\times 3000$ matrices for TSK and MIRK,
when $c=0.9$(a), $c=0.5$(b), $c=0.1$(c) and $c=-0.4$(d).

(a) (b) (c) (b)

Figure 2: RSE versus CPU with $2000\times 1000$ matrices for TSK and MIRK,
when $c=0.9$(a), $c=0.5$(b), $c=0.1$(c) and $c=-0.2$(d).

(a) (b)

Figure 3: CPU versus $c$ with $1000\times 3000$ matrices(a) and $2000\times
1000$ matrices(b) for TSK and MIRK.

(a) (b)

Figure 4:  CPU versus $J$(a) and $I$(b) for TSK and MIRK when $I$ is fixed or
$J$ is fixed.

We firstly take the matrix sizes as $1000\times 3000$ and $2000\times 1000$,
and change the values of $c$. The numbers of iteration steps and the computing
times for both MIRK and TSK methods are reported in Tables 1-2, from which we
see that the MIRK algorithm always outperforms the TSK method for big $c$ in
terms of CPU time, whether $A$ is fat or thin. According to the data, the
speed-up becomes larger as $c$ goes close to 1. When the matrix is fat, the
maximum speed-up can be 1.4088. When the matrix is thin, the maximum speed-up
can be 1.2346. The data in Tables 1-2 is displayed directly in Figures 1-2
which depict the curves of the relative solution error versus CPU time for the
different rows and columns of matrix $A$ and different values of $c$. It can
be concluded from the figures that when $c=0.9$, the two methods have the
largest difference in the decay of the relative solution error and the MIRK
algorithm is obviously better than the TSK method. The improvement effect of
the MIRK algorithm over the TSK method is weakened when $c$ successively
decreases or even takes a negative value. It is observed that the drop curves
almost coincide in Figure 1(d) ($c=-0.4$) and Figure 2(d) ($c=-0.2$), which
means that the improvement effect disappears. These can also be seen in Figure
3, which illustrats that the bigger $c$ is, the more obvious improvement
effect of the MIRK algorithm over the TSK method.

Table 3:  IT and CPU of TSK and MIRK for $I$-by-$J$ matrices with $I=1000$ and different $J$.  $J$ | 2000 | 3000 | 4000 | 5000
---|---|---|---|---
TSK | CPU | 4.2719 | 3.4356 | 8.9900 | 10.5266
IT | $5.0309\times 10^{4}$ | $2.7114\times 10^{4}$ | $2.1190\times 10^{4}$ | $1.8685\times 10^{4}$
MIRK | CPU | 3.0731 | 2.5031 | 6.7650 | 7.8109
IT | $6.7512\times 10^{4}$ | $3.7046\times 10^{4}$ | $2.9676\times 10^{4}$ | $2.6438\times 10^{4}$
speed-up | 1.3901 | 1.3725 | 1.3289 | 1.3477
Table 4:  IT and CPU of TSK and MIRK for $I$-by-$J$ matrices with $J=1000$ and different $I$.  $I$ | 2000 | 3000 | 4000 | 5000
---|---|---|---|---
TSK | CPU | 3.3872 | 2.1891 | 8.8016 | 8.6222
IT | $5.0554\times 10^{4}$ | $2.7714\times 10^{4}$ | $2.1325\times 10^{4}$ | $1.8773\times 10^{4}$
MIRK | CPU | 2.7056 | 1.8053 | 7.4206 | 7.8966
IT | $6.7936\times 10^{4}$ | $3.7014\times 10^{4}$ | $2.9572\times 10^{4}$ | $2.6326\times 10^{4}$
speed-up | 1.2519 | 1.2126 | 1.1861 | 1.0919

(a) (b)

Figure 5:  CPU versus $J$(a) and $I$(b) for TSK and MIRK when $I$ is fixed or
$J$ is fixed.

We also fix $c=0.9$ and observe the performance of two algorithms for
different row $I$ and and column $J$ of $A$. We set $I=1000$ and change the
values of $J$ in Table 3, which shows that in the fat case the speed-up can
attain at most 1.3901 and least 1.3289. We fix $J=1000$ and change the values
of $I$ in Table 4 which illustrates that in thin case, the speed-up can be at
most 1.2519 and at least 1.0919. It is concluded that the improvement effect
is better for the fat case and decreases as $A$ becomes fatter or thinner. The
above observations are visually plotted in Figure 4 which shows that the
curves of the CPU time versus $J$ or $I$ of the MIRK algorithm are always
lower than those of the TSK method no matter what the size of the matrix is.

##  6 Some concluding remarks

In this paper, we firstly give a new convergence rate of the two-subsapce
Kaczmarz method by regarding it as an alternated inertial randomized Kaczmarz
algorithm. Secondly, we provide a multi-step inertial randomized Kaczmarz
algorithm which accelerates the alternated inertial randomized Kaczmarz
algorithm. The preliminary numerical examples are presented to support the
theory results.

It is worth to investigate the combination of the inertial extrapolation with
variants of Kaczmarz method, for example the greedy randomized Kaczmarz method
[2] and the randomized block Kaczmarz method [19].

## References

  * [1] Attouch, H., Peypouquet, J., The rate of convergence of Nesterov’s acclerated forward-backward method is actually faster than $1/k^{2}$, SIAM J. Optim. 26(3), (2016), 1824–1834.
  * [2] Bai, Z.Z., Wu, W.T., On greedy randomized Kaczmarz method for solving large sparse linear systems, SIAM J. Sci. Comput., 40, (2018), A592-A606.
  * [3] Bauschke, H.H., Combettes, P.L., Convex analysis and monotone operator theory in Hilbert spaces. Springer 2010.
  * [4] Beck, A., Teboulle, M., A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM J. Imaging Sci. 2(1)(2009), 183–202.
  * [5] Browne, J., DePierro, A., A row-action alternative to the EM algorithm for maximizing likelihoods in emission tomography, IEEE Trans. Med. Imag. 15, (1996), 687–699.
  * [6] Byrne, C.L., Iterative Optimization in Inverse Problems. CRC Press, Boca Raton (2014).
  * [7] Censor, Y., Row-action methods for huge and sparse systems and their applications. SIAM Review. 23, (1981), 444–464.
  * [8] Dong, Q.L., Huang, J., Li, X.H., Cho, Y.J., Rassias, Th.M., MiKM: Multi-step inertial Krasnosel’skiǐ–Mann algorithm and its applications, J. Global Optim. 73(4), (2019), 801–824.
  * [9] Gower, R. M., Richtárik, P., Stochastic dual ascent for solving linear systems, arXiv: 1512.06890(2015), 28pages.
  * [10] Goebel, K., Reich, S., Uniform convexity, hyperbolic geometry and non-expansive mappings. M. Dekker (1984).
  * [11] Goldstein, T., O’Donoghue, B., Setzer, S., Baraniuk, R., Fast alternating direction optimization methods, SIAM J. Imaging Sc., 7(3), (2014), 1588–1623.
  * [12] He, S.N., Yang, C.P., Duan, P.C., Realization of the hybrid method for Mann iterations. Applied Mathematics and Computation. 217, (2010), 4239–4247.
  * [13] Herman, G.T., Meyer, L., Algebraic reconstruction technique can be made computationlly efficient, IEEE T. Med. Imaging, 12, (1993), 600–609.
  * [14] Iutzeler, F., Malick, J., On the proximal gradient algorithm with alternated inertia, J. Optimiz. Theory App. 176, (2018), 688–710.
  * [15] Kaczmarz, S., Angenäherte Auflösung von Systemen linearer Gleichungen, Bull. Int. Acad. Polon. Sci. Lett. A, 35, (1937), 355–357.
  * [16] Liang, J.W., Convergence Rates of First–Order Operator Splitting Methods. Optimization and Control [math.OC]. Normandie Université; GREYC CNRS UMR 6072, 2016. English.
  * [17] Ma, A., Needell, D., Ramdas, A., Convergence properties of the randomized extended Gauss-Seidel and Kaczmarz methods, SIAM J. Matrix Anal. Appl., 36, (2015), 1590–1604.
  * [18] Mu, Z., Peng, Y., A note on the inertial proximal point method, Stat. Optim. Inf. Comput. 3, (2015), 241–248.
  * [19] Needell, D., Tropp, J.A., Paved with good intentions: analysis of a randomized block Kaczmarz method, Linear Algebra Appl. 441, (2014), 199–221.
  * [20] Needell, D., Ward, R., Two-subspace projection method for coherent overdetermined systems, J. Fourier Anal. Appl., 19, (2013), 256–269.
  * [21] Nesterov, Y.E., A method for solving the convex programming problem with convergence rate $O(1/k^{2})$, Dokl. Akad. Nauk SSSR, 269, (1983), 543–547 (in Russian).
  * [22] Polyak, B.T., Some methods of speeding up the convergence of iteration methods, U.S.S.R. Comput. Math. Math. Phys. 4, (1964), 1–17.
  * [23] Shehu, Y., Iyiola, O.S., Projection methods with alternating inertial steps for variational inequalities: weak and linear convergence, Appl. Numer. Math. 157, (2020), 315–337.
  * [24] Strohmer, T., Vershynin, R., A randomized Kaczmarz algorithm with exponential convergence, J. Fourier Anal. Appl., 15, (2009), 262–278.