11institutetext: Lionel Tondji 22institutetext: Institute for Analysis and
Algebra, TU Braunschweig, 38092 Braunschweig, Germany,
22email<EMAIL_ADDRESS>33institutetext: Dirk A. Lorenz
44institutetext: Institute for Analysis and Algebra, TU Braunschweig, 38092
Braunschweig, Germany,
44email<EMAIL_ADDRESS>
# Faster Randomized Block Sparse Kaczmarz by Averaging††thanks: Funding: The
work of the authors has been supported by the ITN-ETN project TraDE-OPT funded
by the European Union’s Horizon 2020 research and innovation programme under
the Marie Skłodowska-Curie grant agreement No 861137. This work represents
only the author’s view and the European Commission is not responsible for any
use that may be made of the information it contains.

Lionel Tondji  Dirk A. Lorenz

(Received: date / Accepted: date)

###### Abstract

The standard randomized sparse Kaczmarz (RSK) method is an algorithm to
compute sparse solutions of linear systems of equations and uses sequential
updates, and thus, does not take advantage of parallel computations. In this
work, we introduce a parallel (mini batch) version of RSK based on averaging
several Kaczmarz steps. Naturally, this method allows for parallelization and
we show that it can also leverage large over-relaxation. We prove linear
expected convergence and show that, given that parallel computations can be
exploited, the method provably provides faster convergence than the standard
method. This method can also be viewed as a variant of the linearized Bregman
algorithm, a randomized dual block coordinate descent update, a stochastic
mirror descent update, or a relaxed version of RSK and we recover the standard
RSK method when the batch size is equal to one. We also provide estimates for
inconsistent systems and show that the iterates converges to an error in the
order of the noise level. Finally, numerical examples illustrate the benefits
of the new algorithm.

###### Keywords:

Randomized Kaczmarz Sparse solutions Parallel methods

###### MSC:

65F10 68W20 68W10 90C25

##  1 Introduction

In this work we are concerned with the fundamental problem of approximating
sparse solutions of large scale linear systems of the form

(1) $\mathbf{A}x=b$

with matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$ and right hand side
$b\in\mathbb{R}^{m}$. Linear systems like (1) arise in several fields of
engineering and physics problems, such as sensor networks tropp2011improved ,
signal processing hanke1990acceleration , partial differential equations
olshanskii2014iterative , filtering khan2007distributed , computerized
tomography hounsfield1973computerized , optimal control
patrascu2017nonasymptotic , inverse problems jiao2017preasymptotic ;
rabelo2022kaczmarzillposed  and machine learning, to name just a few. When
$\mathbf{A}$ is too large to fit in memory, direct methods for solving
equation (1) are not feasible and iterative methods are preferred. As long as
one can afford full matrix vector products or the system matrix fits in
memory, Krylov methods including the conjugate gradient (CG) algorithms
stiefel1952methods  are the industrial standard. On the other hand, randomized
methods such as the randomized (block) Kaczmarz Kac37 ; SV09  and coordinate
descent method nesterov2012efficiency  are effective if a single matrix vector
product is too expensive and in some situations are even more efficient than
CG method (see, e.g., SV09  for an example). Linear convergence of the
Kaczmarz method has been shown in the randomized case in SV09  and
popa2018convrate  analyzes convergence rates in the deterministic case.
Moreover block Kaczmarz methods moorman2021randomized ; necoara2019faster ;
needell2014paved ; P15 ; richtarik2020stochastic  have received much attention
for their high efficiency for solving (1) and distributed implementations. In
this work we propose and analyze the randomized sparse Kaczmarz method SL19
and show that parallel computations and averaging as in moorman2021randomized
lead to faster convergence.

###  1.1 Related Work

####  Randomized Kaczmarz.

In the large data regime, the randomized Kaczmarz (RK) is a popular iterative
method for solving linear systems. In each iteration111We use subscript
indices for components of a vector, columns or rows of a matrix, and also as
iteration indices. But the meaning should always be clear from the context., a
row vector $a_{i}^{T}$ of $\mathbf{A}$ is chosen at random from the system (1)
and the current iterate $x_{k}$ is projected onto the solution space of that
equation to obtain $x_{k+1}$. Geometrically, at each iteration

$x_{k+1}=\operatorname*{argmin}_{x\in\mathbb{R}^{n}}\,\|x-x_{k}\|_{2}^{2}\quad\text{s.t.}\quad\langle
a_{i},x\rangle=b_{i}.$

It has been observed that the convergence of RK method can be accelerated by
introducing relaxation. In a relaxed variant of RK, a step is taken in the
direction of this projection with the size of the step depending on a
relaxation parameter. Explicitly, the relaxed RK update is given by

(2) $\displaystyle x_{k+1}=x_{k}-w_{k,i}\tfrac{\langle a_{i}\,,\,x_{k}\rangle-
b_{i}}{\|a_{i}\|_{2}^{2}}\cdot a_{i}\,,$

with initial values $x_{0}=0$, where the $w_{k,i},\,i\in\\{1,\dots,m\\}$ are
relaxation parameters. Note that this update rule requires low cost per
iteration and storage of order $\mathcal{O}(n)$. For consistent systems the
relaxation parameters must satisfy

(3)
$0<\liminf_{k\rightarrow\infty}w_{k,i}\leq\limsup_{k\rightarrow\infty}w_{k,i}<2$

to ensure convergence 10.1145/359340.359351 . Fixing the relaxation parameters
$w_{k,i}=1$ for all iterations $k$ and indices $i$ leads to the standard RK
method. In necoara2019faster , a block Kaczmarz variant under the name
randomized block Kaczmarz (RBK) has been analyzed. Linear convergence in
expectation was shown for consistent systems of equations, with a rate
depending on the geometric properties of the matrix, its submatrices, and on
the size of the blocks. The convergence rate given in necoara2019faster
depends on the block size and the stochastic conditioning parameter of the
most ill-conditioned block of the partition when a partition is used and on
the most ill-conditioned block of the entire matrix $\mathbf{A}$ when the
indices are sampled i.i.d. with replacement. The paper needell2014paved
considers more general sampling strategies such as sampling from a partition
of the rows of the matrix. In du2020randomized  the authors investigate an
extension of the randomized averaged block Kaczmarz method that can solve
least squares problems. A parallel version of RK where a weighted average of
independent updates is used was studied in moorman2021randomized . They showed
that as the number of threads increases, the rate of convergence improves and
the convergence horizon for inconsistent systems decreases. Another more
general class of block methods are sketch-and-project methods
richtarik2020stochastic ; gower2019adaptive . For a linear system
$\mathbf{A}x=b$, sketch-and-project methods iteratively project the current
iterate onto the solution space of a sketched subsystem
$\mathbf{S}^{T}\mathbf{A}x=\mathbf{S}^{T}b$. In particular, RK is a sketch-
and-project method with $\mathbf{S}$ being rows of the identity matrix.

#### Randomized Sparse Kaczmarz.

Recently, a new variant of the standard RK method namely the randomized sparse
Kaczmarz method (RSK) SL19 ; LWSM14 ; P15  with almost the same low cost and
storage requirements has shown good performance in approximating sparse
solutions of large consistent linear systems. It uses two variables
$x_{k}^{*}$ and $x_{k}$ and the relaxed RSK update is given by

(4) $\displaystyle x_{k+1}^{*}$ $\displaystyle=x_{k}^{*}-w_{k,i}\tfrac{\langle
a_{i}\,,\,x_{k}\rangle-b_{i}}{\|a_{i}\|_{2}^{2}}\cdot a_{i}\,,$ $\displaystyle
x_{k+1}$ $\displaystyle=S_{\lambda}(x_{k+1}^{*})$

with initial values $x_{0}=x_{0}^{*}=0$, $\lambda>0$, and the soft shrinkage
operator

$S_{\lambda}(x)=\max\\{|x|-\lambda,0\\}\cdot\operatorname{sign}(x)\,.$

Fixing the relaxation parameters $w_{k,i}=1$ for all iterations $k$ and
indices $i$ lead to the standard RSK method. For consistent systems the
iterates of the standard RSK method converge in expectation to the solution of
the regularized _Basis Pursuit Problem_

(5)
$\min_{x\in\mathbb{R}^{n}}\lambda\cdot\|x\|_{1}+\tfrac{1}{2}\cdot\|x\|_{2}^{2}\quad\mbox{s.t.}\quad\mathbf{A}x=b\,.$

The advantage of (5) is the strong convexity property of the objective
function. It is important to note that (see friedlander2008exact ) for a large
but finite parameter $\lambda>0$, the solution of (5) gives a solution of

(6) $\min_{x\in\mathbb{R}^{n}}\|x\|_{1}\quad\mbox{s.t.}\quad\mathbf{A}x=b\,,$

which is the famous _Basis Pursuit Problem_ chen1998basispursuit . The
$\ell_{1}$-norm has been used in many applications, with the goal of obtaining
sparse or even sparsest solutions of underdetermined systems of linear
equations and least-squares problems which is the basis of the theory of
compressed sensing candes2006compressive ; donoho2006compressed . The sparsest
solution is given by minimizing the so-called zero-norm, $\|x\|_{0}$ counting
the number of nonzero components in $x$. However, it is computationally
intractable even for the simplest instances due to its combinatorial nature
(it is even strongly NP-complete, see problem [MP5] in the seminal reference
garey1979NP ). In candes2006robust ; donoho2005sparse  reasonable conditions
are given under which a solution of (6) is a sparsest solution. In
schopfer2022extended , an extension of the RSK with linear expected
convergence has been proposed for solving sparse least squares and impulsive
noise problems while requiring only one additional column of the system matrix
in each iteration.

#### Block Sparse Kaczmarz Methods.

In this setting, a subset of rows $\mathbf{A}_{\tau_{k}}$ is used at each
iteration, with $\tau_{k}\subseteq\left\\{1,\dots,m\right\\}$ and
$|\tau_{k}|>1$ where $|\tau_{k}|$ denote the cardinality of the set of indices
$\tau_{k}$. We usually have two approaches. The first variant is simply a
_block generalization_ of the basic sparse Kaczmarz and the update is given by

(7) $\displaystyle x_{k+1}^{*}$
$\displaystyle=x_{k}^{*}-w_{k,i}\tfrac{\mathbf{A}_{\tau_{k}}^{T}(\mathbf{A}_{\tau_{k}}x_{k}-b_{\tau_{k}})}{\|\mathbf{A}_{\tau_{k}}\|^{2}_{2}}\,,$
$\displaystyle x_{k+1}$ $\displaystyle=S_{\lambda}(x_{k+1}^{*})$

Such block variants are considered, e.g., in necoara2019faster ;
needell2014paved  (with $\lambda=0$) and in LSW14 ; P15  for $\lambda\geq 0$
and we refer to these iterative process as _block sparse Kaczmarz_ method.
When $|\tau_{k}|=m$, we refer to (7) as the _linearized Bregman method_
cai2009convergence ; yin2010analysis ; LSW14  and when $|\tau_{k}|=1$ as the
_randomized sparse Kaczmarz method_ LWSM14 ; SL19 . The main drawback of (7)
is that it is not adequate for distributed implementations. The second variant
of block sparse Kaczmarz can take advantage of distributed computing: each
iteration takes $\eta$ steps of the relaxed randomized sparse Kaczmarz,
independently in parallel, averages the results and applies the soft shrinkage
to form the next iterate. This leads to the following iteration:

(8) $\displaystyle x_{k+1}^{*}$
$\displaystyle=x_{k}^{*}-\frac{1}{\eta}\sum_{i\in\tau_{k}}w_{i}\tfrac{\langle
a_{i}\,,\,x_{k}\rangle-b_{i}}{\|a_{i}\|_{2}^{2}}\cdot a_{i}\,,$ $\displaystyle
x_{k+1}$ $\displaystyle=S_{\lambda}(x_{k+1}^{*})$

with initial values $x_{0}=x_{0}^{*}=0$, where
$\tau_{k}\subseteq\left\\{1,\dots,m\right\\}$ denotes a random set of $\eta$
row indices sampled with replacement (and the $i$-th row is chosen with
probability $p_{i}$) and $w_{i}$ represents the weight corresponding to the
$i$-th row. Such block variants are considered, e.g., in
richtarik2020stochastic ; gower2019adaptive ; moorman2021randomized ;
necoara2019faster ; miao2022greedy  (with $\lambda=0$). Method (8) is the main
method presented and analyzed in this paper and we refer to it as the
_randomized sparse Kaczmarz with averaging_ (RSKA) with more details in
Algorithm 1. Note that the update (8) is easy to implement on distributed
computing units, and it is comparable in terms of cost per iteration to the
basic sparse Kaczmarz update i.e., of order $\mathcal{O}(\eta n)$. If
$\tau_{k}$ is a set of one index, we recover the relaxed RSK method and if in
addition the weights are chosen as $w_{i}=1$ for
$i\in\left\\{1,\dots,m\right\\}$, we recover the standard RSK method.

Algorithm 1 Randomized Sparse Kaczmarz with Averaging (RSKA)

1:starting points $x_{0}=x_{0}^{*}=0\in\mathbb{R}^{n}$, matrix
$\mathbf{A}\in\mathbb{R}^{m\times n}$ with rows
$0\not=a_{i}^{T}\in\mathbb{R}^{n}$ and vector $b\in\mathbb{R}^{m}$, batch size
$\eta$, weights $\\{w_{i}\\}_{i=1}^{m}\in\mathbb{R}$ and probabilities $p_{i}$

2:(approximate) solution of
$\min_{x\in\mathbb{R}^{n}}\lambda\|x\|_{1}+\frac{1}{2}\|x\|_{2}^{2}$ s.t.
$\mathbf{A}x=b$

3:initialize $k=0$

4:repeat

5: $\tau_{k}\leftarrow\eta$ indices sampled from $\left\\{1,\dots,m\right\\}$

6: Compute $\delta_{k}=\tfrac{1}{\eta}\sum_{i\in\tau_{k}}w_{i}\tfrac{\langle
a_{i},x_{k}\rangle-b_{i}}{\|a_{i}\|^{2}_{2}}\cdot a_{i}$ $\triangleright$ In
parallel

7: Update $x_{k+1}^{*}=x_{k}^{*}-\delta_{k}$

8: Update $x_{k+1}=S_{\lambda}(x^{*}_{k+1})$

9: Increment $k=k+1$

10:until a stopping criterion is satisfied

###  1.2  Contribution and Organization.

To the best of our knowledge, the proposed block variant (8) have not yet been
proposed and analyzed for the randomized sparse Kaczmarz method. In this work,
we make the following contributions:

  * •

We propose a mini batch version termed _RSKA_ of the randomized sparse
Kaczmarz method, a general algorithm that unifies a variety of other methods
such as the randomized Kaczmarz and the randomized sparse Kaczmarz with their
relaxed variants. It is theoretically well-motivated, can exploit parallel
computation and converges linearly in expectation.

  * •

We prove that our proposal leads to faster convergence than its standard
counterpart. We also validate this empirically and we provide implementations
of our algorithm in Python.

The remainder of the paper is organized as follows. Section 2 provides a brief
overview on convexity and Bregman distances. In section 3 we give several
interpretations of our method. Section 4 provides convergence guarantees for
our proposed method. In Section 5, numerical experiments demonstrate the
effectiveness of RSKA and provides several insights regarding its behavior and
its hyper-parameters. Finally, Section 6 draws some conclusions.

###  1.3 Notation

In this section we introduce notation that will be used throughout. The first
$m$ integers are denoted by $[m]\overset{\text{def}}{=}\\{1,2,\dots,m\\}$.
Given a symmetric positive definite matrix $\mathbf{B},$ we equip the space
$\mathbb{R}^{n}$ with the Euclidean inner product defined by

$\langle x,y\rangle_{\mathbf{B}}\overset{\text{def}}{=}\langle
x,\mathbf{B}y\rangle=\sum_{i,j\in[n]}x_{i}\mathbf{B}_{ij}y_{j},\quad
x,y\in\mathbb{R}^{n}$

We also define the induced norm:
$\|.\|_{\mathbf{B}}^{2}\overset{\text{def}}{=}\langle\cdot,\cdot\rangle_{\mathbf{B}}$
and use the short-hand notation $\|.\|$ to mean $\|.\|_{\mathbf{I}}$ to denote
the standard 2-norm. Let $\mathbf{A}\in\mathbb{R}^{m\times n}$ be a real
matrix. By $\mathbf{Range}(\mathbf{A}),\|\mathbf{A}\|_{F}$ and $a_{i}^{T}$ we
denote its range space, Frobenius norm and $i$th row respectively, and by
$A^{{\dagger}}$ we denote the (Moore-Penrose) pseudo inverse. The indicator
function of a set $C$ is denoted by

$\delta_{C}(x)\overset{\text{def}}{=}\begin{cases}0,&\text{if }x\in C\\\
+\infty,&\text{if }x\not\in C\end{cases}$

By $e_{i}$ we denote the $i$th column of the identity matrix
$\mathbf{I}_{n}\in\mathbb{R}^{n\times n}$. For a random vector $x_{i}$ that
depends on a random index $i\in[q]$ (where $i$ is chosen with probability
$p_{i}$) we denote $\mathbb{E}_{i\sim
p}[x_{i}]\overset{\text{def}}{=}\sum_{i\in[q]}p_{i}x_{i}$ and we will just
write $\mathbb{E}[x_{i}]$ when the probability distribution is clear from the
context. Let $\sigma_{i}(\mathbf{A})$ be the $i$th singular value of $A$ (when
ordered decreasingly) and $\sigma_{\min}(\mathbf{A})$ and
$\sigma_{\max}(\mathbf{A})$ be the smallest and largest singular values of
$\mathbf{A}$, respectively. They are given by

(9)
$\sigma_{\min}(\mathbf{A})\overset{\text{def}}{=}\operatorname*{\min}_{x\in\mathbb{R}^{n},\
x\neq
0}\,\frac{\|\mathbf{A}x\|_{2}}{\|x\|_{2}}\quad\text{and}\quad\sigma_{\max}(\mathbf{A})\overset{\text{def}}{=}\sigma_{1}(\mathbf{A})\overset{\text{def}}{=}\operatorname*{\max}_{x\in\mathbb{R}^{n},\
x\neq 0}\,\frac{\|\mathbf{A}x\|_{2}}{\|x\|_{2}}$

Finally, a result we will need, if $\mathbf{A}$ is a symmetric positive semi-
definite matrix the largest singular value of $\mathbf{A}$ (the $L_{2}$
induced matrix norm or the spectral norm ) can be defined instead as

(10)
$\|\mathbf{A}\|_{2}\overset{\text{def}}{=}\sigma_{\max}(\mathbf{A})=\max_{x\in\mathbb{R}^{n},\
x\neq
0}\,\frac{|\langle\mathbf{A}x,x\rangle|}{\|x\|_{2}^{2}}=\max_{x\in\mathbb{R}^{n},\
x\neq 0}\,\frac{\|\mathbf{A}x\|_{2}}{\|x\|_{2}}.$

Thus clearly

$\frac{|\langle
x,x\rangle_{\mathbf{A}}|}{\|x\|_{2}^{2}}\leq\sigma_{\max}(\mathbf{A}).$

Let
$\tilde{\sigma}_{\min}(\mathbf{A})\overset{\text{def}}{=}\min\\{\sigma_{\min}(\mathbf{A}_{J})\mid
J\subseteq[n],\mathbf{A}_{J}\neq 0\\}$ where $\mathbf{A}_{J}$ denotes the
submatrix of $\mathbf{A}$ that is built up by the columns indexed by $J$ and
$|x|_{\mathrm{min}}\overset{\text{def}}{=}\min\\{|x_{j}|\mid x_{j}\neq 0\\}$.

##  2 Basic notions

At first we recall some well known concepts and properties of convex functions
and Bregman distances and later give upper-bounds of singular values of sum of
matrices. Let $f:\mathbb{R}^{n}\to\mathbb{R}$ be convex (note that we assume
that $f$ is finite everywhere, hence also continuous). The _subdifferential_
of $f$ is defined by

$\partial f(x)\overset{\text{def}}{=}\\{x^{*}\in\mathbb{R}^{n}\,|\,f(y)\geq
f(x)+\langle x^{*}\,,\,y-x\rangle\>\mbox{for all $y\in\mathbb{R}^{n}$}\\}$

at any $x\in\mathbb{R}^{n}$ is nonempty, compact and convex.

The function $f:\mathbb{R}^{n}\to\mathbb{R}$ is said to be _$\alpha$ -strongly
convex_, if for all $x,y\in\mathbb{R}^{n}$ and subgradients $x^{*}\in\partial
f(x)$ we have

$f(y)\geq f(x)+\langle
x^{*}\,,\,y-x\rangle+\tfrac{\alpha}{2}\cdot\|y-x\|_{2}^{2}\,.$

If $f$ is $\alpha$-strongly convex then $f$ is coercive, i.e.

$\lim_{\|x\|_{2}\to\infty}f(x)=\infty\,,$

and its _Fenchel conjugate_ $f^{*}:\mathbb{R}^{n}\to\mathbb{R}$ given by

$f^{*}(x^{*})\overset{\text{def}}{=}\sup_{y\in\mathbb{R}^{n}}\langle
x^{*}\,,\,y\rangle-f(y)$

is also convex, finite everywhere and coercive.

Additionally, $f^{*}$ is differentiable with a _Lipschitz-continuous gradient_
with constant $L_{f^{*}}=\frac{1}{\alpha}$, i.e. for all
$x^{*},y^{*}\in\mathbb{R}^{n}$ we have

$\|\nabla f^{*}(x^{*})-\nabla f^{*}(y^{*})\|_{2}\leq
L_{f^{*}}\cdot\|x^{*}-y^{*}\|_{2}\,,$

which implies the estimate

(11) $f^{*}(y^{*})\leq f^{*}(x^{*})-\langle\nabla
f^{*}(x^{*})\,,\,y^{*}-x^{*}\rangle+\tfrac{L_{f^{*}}}{2}\cdot\|x^{*}-y^{*}\|_{2}^{2}\,.$

######  Example 1.

The objective function

(12)
$f(x)\overset{\text{def}}{=}\lambda\cdot\|x\|_{1}+\tfrac{1}{2}\cdot\|x\|_{2}^{2}$

is strongly convex with constant $\alpha=1$ and its conjugate function can be
computed with the soft shrinkage operator as

$f^{*}(x^{*})=\tfrac{1}{2}\cdot\|S_{\lambda}(x^{*})\|_{2}^{2}\quad\mbox{with}\quad\nabla
f^{*}(x^{*})=S_{\lambda}(x^{*})\,.$

######  Definition 1.

The _Bregman distance_ $D_{f}^{x^{*}}(x,y)$ between $x,y\in\mathbb{R}^{n}$
with respect to $f$ and a subgradient $x^{*}\in\partial f(x)$ is defined as

$D_{f}^{x^{*}}(x,y)\overset{\text{def}}{=}f(y)-f(x)-\langle
x^{*}\,,\,y-x\rangle\,.$

Fenchel’s equality states that $f(x)+f^{*}(x^{*})=\langle x\,,\,x^{*}\rangle$
if $x^{*}\in\partial f(x)$ and implies that the Bregman distance can be
written as

$D_{f}^{x^{*}}(x,y)=f^{*}(x^{*})-\langle x^{*}\,,\,y\rangle+f(y)\,.$

######  Example 2 (cf. SL19 ).

For $f(x)=\frac{1}{2}\cdot\|x\|_{2}^{2}$ we just have $\partial f(x)=\\{x\\}$
and $D_{f}^{x^{*}}(x,y)=\frac{1}{2}\|x-y\|_{2}^{2}$. For
$f(x)=\lambda\cdot\|x\|_{1}+\tfrac{1}{2}\cdot\|x\|_{2}^{2}$ and any
$x^{*}=x+\lambda\cdot s\in\partial f(x)$ we have

$D_{f}^{x^{*}}(x,y)=\frac{1}{2}\cdot\|x-y\|_{2}^{2}+\lambda\cdot(\|y\|_{1}-\langle
s\,,\,y\rangle)\,.$

The following properties are crucial for the convergence analysis of the
randomized algorithms. They immediately follow from the definition of the
Bregman distance and the assumption of strong convexity of $f$, cf. LSW14 .
For all $x,y,z\in\mathbb{R}^{n}$ and $x^{*}\in\partial f(x)$,
$y^{*}\in\partial f(y)$, $z^{*}\in\partial f(z)$ we have

(13) $\frac{\alpha}{2}\|x-y\|_{2}^{2}\leq D_{f}^{x^{*}}(x,y)\leq\langle
x^{*}-y^{*}\,,\,x-y\rangle\leq\|x^{*}-y^{*}\|_{2}\cdot\|x-y\|_{2}$ (14)
$D_{f}^{x^{*}}(x,y)+D_{f}^{y^{*}}(y,z)-D_{f}^{x^{*}}(x,z)=\langle
x^{*}-y^{*}\,,\,z-y\rangle$

Note that if $f$ is differentiable with a Lipschitz-continuous gradient, then
we also have the (better) upper estimate $D_{f}^{x^{*}}(x,y)\leq
L_{f}\cdot\|x-y\|_{2}^{2}$, but in general this need not be the case.

The following Theorem, will be use in the convergence analysis more precisely
in Lemma (6).

######  Theorem 2.1 ((horn1994topics, , Theorem 3.3.16(c))).

Let $\mathbf{A},\mathbf{B}\in\mathbb{R}^{m\times n}$ be given and let
$p=\min\\{m,n\\}$. Then it holds for the decreasingly ordered singular values
of $\mathbf{A},\mathbf{B},\mathbf{A+B}$ that

$|\sigma_{i}(\mathbf{A+B})-\sigma_{i}(\mathbf{A})|\leq\sigma_{1}(\mathbf{B}),\quad\text{for}\quad
i\in[p].$

In particular we have

$\sigma_{i}(\mathbf{A}+\mathbf{B})\geq\sigma_{i}(\mathbf{A})-\sigma_{1}(\mathbf{B}),\quad\text{for}\quad
i\in[p].$

##  3 Interpretations

We can view the randomized sparse Kaczmarz with averaging algorithm as an
optimization method for solving a specific primal or dual optimization
problem. More precisely, the RSKA algorithm is a particular case of the
following.

###  3.1 Randomized Block/Parallel Coordinate Descent

Considering the regularized Basis Pursuit Problem as primal problem

(15)
$\min_{x\in\mathbb{R}^{n}}\lambda\cdot\|x\|_{1}+\tfrac{1}{2}\cdot\|x\|_{2}^{2}\quad\mbox{s.t.}\quad\mathbf{A}x=b\,.$

The dual of optimization problem (15) takes the form of a quadratic program:

(16)
$\min_{y\in\mathbb{R}^{m}}\tfrac{1}{2}\cdot\|S_{\lambda}(\mathbf{A}^{T}y)\|_{2}^{2}-\langle
b\,,\,y\rangle\,.$

where the primal variable $x$ and the dual variable $y$ are related through
the relation $x=S_{\lambda}(\mathbf{A}^{T}y)$. Let us define the primal and
dual objective functions

$f(x)=\lambda\cdot\|x\|_{1}+\tfrac{1}{2}\cdot\|x\|_{2}^{2}+\delta_{\\{0\\}}(b-\mathbf{A}x)$

and

$g(y)=\tfrac{1}{2}\cdot\|S_{\lambda}(\mathbf{A}^{T}y)\|_{2}^{2}-\langle
b\,,\,y\rangle,$

respectively. One iteration of the RSKA algorithm can be viewed as one step of
the randomized block coordinate descent (RBCD) applied to the dual problem
(16) when the weights $w_{i}$ are chosen in a particular form P15 . More
formally a negative gradient step in the random $i$-th component of $y$ having
$\nabla_{i_{k}}g(y)=a_{i_{k}}^{T}S_{\lambda}(\mathbf{A}^{T}y)-b_{i_{k}}$ with
step size $t_{k}=\tfrac{1}{\|a_{i_{k}}\|^{2}_{2}}$ yields

$\displaystyle y_{k+1}$
$\displaystyle=y_{k}-\tfrac{1}{\|a_{i_{k}}\|^{2}_{2}}\nabla_{i_{k}}g(y_{k})e_{i_{k}}$

We easily recover (4) by simply multiplying this update with $\mathbf{A}^{T}$
and using the relation between the primal and dual variables given by
$x=S_{\lambda}(\mathbf{A}^{T}y)$. Consider the dual problem (16). If we choose
the particular weights
$w_{k,i}=\tfrac{\|a_{i}\|^{2}_{2}}{\sum_{i\in\tau_{k}}\|a_{i}\|^{2}_{2}}$, the
block coordinate descent method applied to $g$ from (16) reads

$\displaystyle x_{k+1}^{*}$
$\displaystyle=x_{k}^{*}-\frac{1}{\eta\sum_{i\in\tau_{k}}\|a_{i}\|^{2}_{2}}\sum_{i\in\tau_{k}}(\langle
a_{i}\,,\,x_{k}\rangle-b_{i})\cdot a_{i}\,,$ $\displaystyle x_{k+1}$
$\displaystyle=S_{\lambda}(x_{k+1}^{*})$

Parallel coordinate descent richtarik2016parallel  applied to the dual problem
(16) with learning rate $t_{k}=\frac{1}{\eta\|a_{i}\|^{2}_{2}}$ yields

(17)
$y_{k+1}=y_{k}-\sum_{i\in\tau_{k}}\frac{1}{\eta\|a_{i}\|^{2}_{2}}\nabla_{i}g(y_{k})e_{i}$

In Update (17) only coordinate $i\in\tau_{k}$ are updated in $y_{k+1}$ and the
remaining coordinate are unchanged. Multiplying this update with
$\mathbf{A}^{T}$ and using the relation between the primal and dual variables,
we recover (4) with particular weights $w_{i}=1$. However, for general weights
$w_{k}$, the RSKA algorithm cannot be interpreted in these ways, and thus our
scheme is more general. It is important to note that in P15  for the
randomized block sparse Kaczmarz method of type $(\ref{eq:bsk})$ sublinear
convergence rates have been obtained by identifying the iteration as a
randomized block coordinate gradient descent method applied to the objective
function $g$ of the unconstrained dual of $f$. However, the rates given in P15
are in terms of the dual objective function $g$, and not of the primal
iterates only, although, as mentioned there in the conclusions, the
experimental results indicate that such rates also hold for the primal
iterates.

###  3.2 Stochastic Mirror Descent with Stochastic Polyak Stepsize

The stochastic mirror descent (SMD) method and its variants beck2003mirror ;
lan2012validation ; nemirovski2009robust  is one of the most widely used
family of algorithms in stochastic optimization for non-smooth, Lipschitz
continuous– convex and non-convex functions. Starting with the orginal work of
nemirovskij1983problem , SMD has been studied in the context of convex
programming nemirovski:hal-00976649 , saddle-point problems
mertikopoulos2018optimistic , and monotone variational inequalities
mertikopoulos2018stochastic . Now we draw a connection of Algorithm 1 to the
stochastic mirror descent method using stochastic Polyak stepsizes. We
consider a set of sketching matrices $\mathbf{S}_{i}\in\mathbb{R}^{m\times s}$
($i\in[m]$), define
$\mathbf{Z}_{i}\overset{\text{def}}{=}\mathbf{S}_{i}(\mathbf{S}_{i}^{T}\mathbf{A}\mathbf{A}^{T}\mathbf{S}_{i})^{\dagger}\mathbf{S}_{i}^{T}$
and consider the stochastic convex quadratic $h_{\mathbf{S}_{i}}$

(18)
$h_{\mathbf{S}_{i}}(x)\overset{\text{def}}{=}\frac{1}{2}\|\mathbf{A}x-b\|_{\mathbf{Z}_{i}}^{2}=\frac{1}{2}(\mathbf{A}x-b)^{T}\mathbf{Z_{i}}(\mathbf{A}x-b),$

(recall that $\mathbf{A}\in\mathbb{R}^{m\times n}$, $b\in\mathbb{R}^{m}$).

The general sketched mirror descent update with learning rate $t_{k}$ and a
mirror map $f$ works as follows: For a given iterate $x_{k}$ draw a sketching
matrix $\mathbf{S}_{i_{k}}$ (actually, one draws an index $i_{k}$) at random
and update

(19)
$x_{k+1}=\operatorname*{arg\,min}_{x\in\mathbb{R}^{n}}\,\Bigg{\\{}\langle\nabla
h_{\mathbf{S}_{i_{k}}}(x_{k}),x-x_{k}\rangle+\frac{1}{t_{k}}D^{x_{k}^{*}}_{f}(x_{k},x)\Bigg{\\}},\quad
x_{k}^{*}\in\partial f(x_{k}),$

which yields to the following update:

(20) $\displaystyle x_{k+1}^{*}$ $\displaystyle=x_{k}^{*}-t_{k}\nabla
h_{\mathbf{S}_{i_{k}}}(x_{k})\,,$ $\displaystyle x_{k+1}$
$\displaystyle=\nabla f^{*}(x_{k+1}^{*}).$

where $f^{*}$ denote the Fenchel conjugate of $f$.

One can show that

(21) $\displaystyle\footnotesize D^{x_{k+1}^{*}}_{f}(x_{k+1},\hat{x})\leq
D^{x_{k}^{*}}_{f}(x_{k},\hat{x})-t_{k}\langle\nabla
h_{\mathbf{S}_{i_{k}}}(x_{k}),x_{k}-\hat{x}\rangle+\frac{t_{k}^{2}}{2}\|\nabla
h_{\mathbf{S}_{i_{k}}}(x_{k})\|^{2}$

(actually, this also follows from Lemma 2 below with $\varPhi(x)=\langle\nabla
h_{\mathbf{S}_{i_{k}}}(x_{k}),x-x_{k}\rangle$, and a strongly convex function
$f$). If we select $t_{k}$ such that the RHS of inequality (21) is minimized,
we obtain

(22) $t_{k}=\frac{\langle\nabla
h_{\mathbf{S}_{i_{k}}}(x_{k}),x_{k}-\hat{x}\rangle}{\|\nabla
h_{\mathbf{S}_{i_{k}}}(x_{k})\|^{2}}.$

Since $\nabla h_{\mathbf{S}}(x)=\mathbf{A}^{T}\mathbf{Z}(\mathbf{A}x-b)$ we
get (cf. richtarik2020stochastic )

(23)
$h_{\mathbf{S}_{i_{k}}}(x)-h_{\mathbf{S}_{i_{k}}}(\hat{x})\overset{h_{\mathbf{S}_{i_{k}}}(\hat{x})=0}{=}h_{\mathbf{S}_{i_{k}}}(x)=\tfrac{1}{2}\|\nabla
h_{\mathbf{S}_{i_{k}}}(x)\|^{2}=\tfrac{1}{2}\langle\nabla
h_{\mathbf{S}_{i_{k}}}(x_{k}),x-\hat{x}\rangle$

we get that the optimal step size is in fact simply

(24)
$t_{k}=\frac{2\big{[}h_{\mathbf{S}_{i_{k}}}(x_{k})-h_{\mathbf{S}_{i_{k}}}(\hat{x})\big{]}}{\|\nabla
h_{\mathbf{S}_{i_{k}}}(x_{k})\|^{2}}.$

This quotient is known as _stochastic mirror Polyak stepsize_ d2021stochastic
, (note that in this particular case, we always get $t_{k}=1$). The randomized
Kaczmarz (4) is equivalent to one step of the stochastic mirror descent (19)
with the mirror Polyak stepsize $t_{k}$ given in (22), whereas the mirror
Polyak stepsize in it general form d2021stochastic ; loizou2021stochastic  is
given by
$t_{k}=\frac{h_{\mathbf{S}_{i_{k}}}(x_{k})-h_{\mathbf{S}_{i_{k}}}(\hat{x})}{c\|\nabla
h_{\mathbf{S}}(x_{k})\|^{2}_{2}}.$ The parameter $0<c\in\mathbb{R}$ in the
step size is an important quantity which can be set theoretically based on the
properties of the function under study. In loizou2021stochastic  it is
suggested that, for optimal convergence, one should select $c=1/2$ for
strongly convex functions and $c=0.2$ for non-convex functions. Morever, the
general RSKA method (see Algorithm 1) falls into the general sketched-and-
project framework in the context of mirror descent with
$\mathbf{Z}\overset{\text{def}}{=}\frac{1}{\eta}\sum_{j\in\tau}w_{j}\mathbf{S}_{j}(\mathbf{S}_{j}^{T}\mathbf{A}\mathbf{A}^{T}\mathbf{S}_{j})^{\dagger}\mathbf{S}_{j}^{T}$
with $\mathbf{S}_{j}=e_{j}$, $t_{k}=1$ and
$f(x)=\lambda\cdot\|x\|_{1}+\tfrac{1}{2}\cdot\|x\|_{2}^{2}$. In fact, the
sketch-and-project form of updates (8) (resp. 20) is given by:

(25) $\displaystyle x_{k+1}^{*}$
$\displaystyle=x_{k}^{*}-\mathbf{A}^{T}\mathbf{Z}(\mathbf{A}x_{k}-b),$
$\displaystyle x_{k+1}$ $\displaystyle=\nabla f^{*}(x_{k+1}^{*}),$

with some random $\tau_{i}\subset[m]$ with cardinality $\eta$.

##  4 Convergence analysis

In this section we show expected linear convergence for the randomized sparse
Kaczmarz with averaging method. Before that we give the update satisfy by the
iterate $x_{k}^{*}$ in Algorithm (1). From line 5 of Algorithm (1), it holds
that:

$\displaystyle x_{k+1}^{*}$
$\displaystyle=x_{k}^{*}-\frac{1}{\eta}\sum_{i\in\tau_{k}}w_{i}\tfrac{\langle
a_{i}\,,\,x_{k}\rangle-b_{i}}{\|a_{i}\|_{2}^{2}}\cdot a_{i}$
$\displaystyle=x^{*}_{k}-\tfrac{1}{\eta}\sum_{i\in\tau_{k}}w_{i}(e_{i}^{T}\mathbf{A})^{T}\cdot\tfrac{e_{i}^{T}(\mathbf{A}x_{k}-b)}{\|a_{i}\|^{2}_{2}}$
$\displaystyle=x^{*}_{k}-\mathbf{A}^{T}\tfrac{1}{\eta}\sum_{i\in\tau_{k}}w_{i}\cdot\tfrac{e_{i}e_{i}^{T}}{\|a_{i}\|^{2}_{2}}(\mathbf{A}x_{k}-b)$

To simplify notation, we define the following matrices.

######  Definition 2.

Let $\mathbf{Diag}(d_{1},d_{2},\dots,d_{m})$ denote the diagonal matrix with
$d_{1},d_{2},\dots,d_{m}$ on the diagonal. We define the following matrices:

  * •

Weighted sampling matrix:

$\displaystyle\mathbf{M}_{k}=\tfrac{1}{\eta}\sum_{i\in\tau_{k}}w_{i}\tfrac{e_{i}e_{i}^{T}}{\|a_{i}\|^{2}_{2}}$

  * •

Normalization matrix:

$\displaystyle\mathbf{D}=\mathbf{Diag}(\|a_{1}\|,\|a_{2}\|,\dots,\|a_{m}\|)$

so that the matrix $\mathbf{D}^{-1}\mathbf{A}$ has rows with unit norm.

  * •

Probability matrix:

$\displaystyle\mathbf{P}=\mathbf{Diag}(p_{1},p_{2},\dots,p_{m})$

where $p_{j}=\mathbb{P}(i=j)$.

  * •

Weight matrix:

$\displaystyle\mathbf{W}=\mathbf{Diag}(w_{1},w_{2},\dots,w_{m})$

where $w_{i}$ represents the weight corresponding to the $i$-th row.

The following lemma and proof are taken from (moorman2021randomized, , Lemma
1) and we include the full proof for completeness. The lemma gives the first
and second moment of the random matrices $\mathbf{M}_{k}$,
$\mathbf{A}^{T}\mathbf{M}_{k}$ respectively and will be use in our convergence
analysis.

######  Lemma 1.

Let $\mathbf{M}_{k},\mathbf{P},\mathbf{W}$ and $\mathbf{D}$ be defined as in
Definition 2. Then

$\displaystyle\mathbb{E}_{k}\left[\mathbf{M}_{k}\right]$
$\displaystyle=\mathbf{P}\mathbf{W}\mathbf{D}^{-2}\quad\text{and}$
$\displaystyle\mathbb{E}_{k}\left[(\mathbf{A}^{T}\mathbf{M}_{k})^{T}\cdot(\mathbf{A}^{T}\mathbf{M}_{k})\right]$
$\displaystyle=\tfrac{1}{\eta}\mathbf{PW^{2}D^{-2}}+(1-\tfrac{1}{\eta})\mathbf{PWD^{-2}AA^{T}PWD^{-2}}.$

###### Proof.

Let $\mathbb{E}_{i}[\cdot]$ denote $\mathbb{E}_{i\sim p}[\cdot]$. From the
definition of the weighted sampling matrix $\mathbf{M}_{k}$ as the weighted
average of the i.i.d. sampling matrices
$\tfrac{e_{i}e_{i}^{T}}{\|a_{i}\|^{2}_{2}},$ we see that

$\displaystyle\mathbb{E}_{k}\left[\mathbf{M}_{k}\right]=\mathbb{E}_{k}\left[\tfrac{1}{\eta}\sum_{i\in\tau_{k}}w_{i}\tfrac{e_{i}e_{i}^{T}}{\|a_{i}\|^{2}_{2}}\right]=\mathbb{E}_{i}\left[w_{i}\tfrac{e_{i}e_{i}^{T}}{\|a_{i}\|^{2}_{2}}\right]=\sum_{i=1}^{m}p_{i}w_{i}\tfrac{e_{i}e_{i}^{T}}{\|a_{i}\|^{2}_{2}}=\mathbf{P}\mathbf{W}\mathbf{D}^{-2}.$

In the same way, we have

$\displaystyle\mathbb{E}_{k}\left[(\mathbf{A}^{T}\mathbf{M}_{k})^{T}\cdot(\mathbf{A}^{T}\mathbf{M}_{k})\right]$
$\displaystyle=\mathbb{E}_{k}\left[(\tfrac{1}{\eta}\sum_{i\in\tau_{k}}w_{i}\tfrac{e_{i}e_{i}^{T}}{\|a_{i}\|^{2}_{2}})\mathbf{A}\cdot\mathbf{A}^{T}\cdot(\tfrac{1}{\eta}\sum_{j\in\tau_{k}}w_{j}\tfrac{e_{j}e_{j}^{T}}{\|a_{j}\|^{2}_{2}})\right]$
$\displaystyle=\mathbb{E}_{k}\left[(\tfrac{1}{\eta}\sum_{i\in\tau_{k}}w_{i}\tfrac{e_{i}a_{i}^{T}}{\|a_{i}\|^{2}_{2}})\cdot(\tfrac{1}{\eta}\sum_{j\in\tau_{k}}w_{j}\tfrac{a_{j}e_{j}^{T}}{\|a_{j}\|^{2}_{2}})\right]$
$\displaystyle=\tfrac{1}{\eta}\mathbb{E}_{i}\left[(w_{i}\tfrac{e_{i}a_{i}^{T}}{\|a_{i}\|^{2}_{2}})\cdot(w_{i}\tfrac{a_{i}e_{i}^{T}}{\|a_{i}\|^{2}_{2}})\right]+(1-\tfrac{1}{\eta})\mathbb{E}_{i}\left[w_{i}\tfrac{e_{i}a_{i}^{T}}{\|a_{i}\|^{2}_{2}}\right]\mathbb{E}_{i}\left[w_{i}\tfrac{a_{i}e_{i}^{T}}{\|a_{i}\|^{2}_{2}}\right]$
$\displaystyle=\tfrac{1}{\eta}\mathbb{E}_{i}\left[w_{i}^{2}\tfrac{e_{i}e_{i}^{T}}{\|a_{i}\|^{2}_{2}}\right]+(1-\tfrac{1}{\eta})\mathbb{E}_{i}\left[w_{i}\tfrac{e_{i}e_{i}^{T}}{\|a_{i}\|^{2}_{2}}\right]\mathbf{A}\mathbf{A}^{T}\mathbb{E}_{i}\left[w_{i}\tfrac{e_{i}e_{i}^{T}}{\|a_{i}\|^{2}_{2}}\right]$
$\displaystyle=\tfrac{1}{\eta}\mathbf{P}\mathbf{W}^{2}\mathbf{D}^{-2}+(1-\tfrac{1}{\eta})\mathbf{P}\mathbf{W}\mathbf{D}^{-2}\mathbf{A}\mathbf{A}^{T}\mathbf{P}\mathbf{W}\mathbf{D}^{-2},$

by separating the cases where $i=j$ from those where $i\neq j$ and utilizing
the independence of the indices sampled in $\tau_{k}$ . ∎

We now present convergence results for the proposed method. We start our
analysis by characterizing the error bound between two consecutive iterates
and the error bound between the Bregman distance of the iterates, the solution
and the residual in the following lemma.

######  Lemma 2.

Let $f,\varPhi:\mathbb{R}^{n}\rightarrow\mathbb{R}\cup\\{+\infty\\}$ be
convex, where $\operatorname{dom}(f)=\mathbb{R}^{n}$ and
$\operatorname{dom}(\varPhi)\neq\emptyset$. Let
$\mathcal{X}\subseteq\mathrm{dom}(\varPhi)$ be nonempty and convex,
$x_{k}\in\mathbb{R}^{n}$, $x_{k}^{*}\in\partial f(x_{k})$. Assume that

$x_{k+1}\in\operatorname*{arg\,min}_{x\in\mathcal{X}}\,\Bigg{\\{}\varPhi(x)+D^{x_{k}^{*}}_{f}(x_{k},x)\Bigg{\\}}.$

Then there exist subgradient $x_{k+1}^{*}\in\partial f(x_{k+1})$ such that it
holds

$\displaystyle\varPhi(y)+D^{x_{k}^{*}}_{f}(x_{k},y)\geq\varPhi(x_{k+1})+D^{x_{k}^{*}}_{f}(x_{k},x_{k+1})+D^{x_{k+1}^{*}}_{f}(x_{k+1},y)$

for any $y\in\mathcal{X}$.

###### Proof.

Let denote by $J(x)=\varPhi(x)+D^{x_{k}^{*}}_{f}(x_{k},x)$. Since $J$ and
$\mathcal{X}$ are convex and $x_{k+1}$ minimizes $J$ over $\mathcal{X}$, there
exists a subgradient $d\in\partial J(x_{k+1})$ such that

$\langle d,y-x_{k+1}\rangle\geq 0\,\,\quad\forall\,y\in\mathcal{X}.$

Since $f$ is finite everywhere, we have
$\operatorname{dom}(D_{f}(\cdot,u))=\mathbb{R}^{n}$ for all
$u\in\mathbb{R}^{n}$. Since $\operatorname{dom}(\varPhi)$ is nonempty and
convex, $\varPhi$ has nonempty relative interior. So the subgradient sum rule
applies and we obtain that

$\partial J(x_{k+1})=\partial\varPhi(x_{k+1})+(\partial
f(x_{k+1})-x_{k}^{*}).$

Hence, there exist subgradients $g\in\partial\varPhi(x_{k+1})$,
$x_{k+1}^{*}\in\partial f(x_{k+1})$ such that

$\langle g+(x_{k+1}^{*}-x_{k}^{*}),\ y-x_{k+1}\rangle\geq
0\,\,\quad\forall\,y\in\mathcal{X}.$

Therefore using the property of the subgradient and (14), we have for all
$y\in\mathcal{X}$

$\displaystyle\varPhi(y)\geq$ $\displaystyle\ \varPhi(x_{k+1})+\langle
g,y-x_{k+1}\rangle$ $\displaystyle\geq$ $\displaystyle\
\varPhi(x_{k+1})+\langle x_{k}^{*}-x_{k+1}^{*},\ y-x_{k+1}\rangle$
$\displaystyle\overset{\text{\tiny\eqref{eq:E}}}{=}$ $\displaystyle\
\varPhi(x_{k+1})+D^{x_{k}^{*}}_{f}(x_{k},x_{k+1})-D^{x_{k}^{*}}_{f}(x_{k},y)+D^{x_{k+1}^{*}}_{f}(x_{k+1},y).$

∎

The following lemma provides an error bound for the Bregman distance.

######  Lemma 3 (SL19 ).

Let $\tilde{\sigma}_{\min}(\mathbf{A})$ and $|\hat{x}|_{\mathrm{min}}$ be
defined as in subsection 1.3. Then for any $x\in\mathbb{R}^{n}$ with $\partial
f(x)\cap\mathbf{Range}(\mathbf{A}^{T})\neq 0$ and for all
$\hat{x}=\mathbf{A}^{T}y\in\partial f(x)\cap\mathbf{Range}(\mathbf{A}^{T})$,
we have

(26)
$D_{f}^{x_{k}^{*}}(x_{k},\hat{x})\leq\gamma\cdot\|\mathbf{A}x_{k}-b\|_{2}^{2}$

where

(27)
$\gamma=\frac{1}{\tilde{\sigma}^{2}_{\text{min}}(\mathbf{A})}\frac{|\hat{x}|_{\mathrm{min}}+2\lambda}{|\hat{x}|_{\mathrm{min}}}$

To effectively use Lemma 3 we need the following assumption which characterize
the coupling between the weight matrix and the probability matrix.

######  Assumption 1.

The weight matrix $\mathbf{W}$ and the probability matrix $\mathbf{P}$ are
linked by the following coupling

$\mathbf{P}\mathbf{W}\mathbf{D}^{-2}=\frac{\alpha}{\|\mathbf{A}\|_{F}^{2}}\mathbf{I}$

some scalar relaxation parameter $\alpha>0$.

Assumption (1) has been used in moorman2021randomized  for the inconsistent
case (i.e. $b-\mathbf{A}\hat{x}\neq 0$) and for $\lambda=0$. We include a
motivation for completeness. As the batch size $\eta$ goes to $\infty$ (recall
that we sample with replacement) we have

$\lim_{\eta\to\infty}\mathbf{M}_{k}=\mathbb{E}_{i\sim
p}\bigg{[}w_{i}\tfrac{e_{i}e_{i}^{T}}{\|a_{i}\|^{2}_{2}}\bigg{]}=\mathbf{P}\mathbf{W}\mathbf{D}^{-2}$

Therefore the averaged RSK update of Eq. (8) approaches the deterministic
update:

$\displaystyle x_{k+1}$
$\displaystyle=(\mathbf{I}-\mathbf{A}^{T}\mathbf{P}\mathbf{W}\mathbf{D}^{-2}\mathbf{A})x_{k}+\mathbf{A}^{T}\mathbf{P}\mathbf{W}\mathbf{D}^{-2}b$
$\displaystyle x_{k+1}-\hat{x}$
$\displaystyle=(\mathbf{I}-\mathbf{A}^{T}\mathbf{P}\mathbf{W}\mathbf{D}^{-2}\mathbf{A})(x_{k}-\hat{x})+\mathbf{A}^{T}\mathbf{P}\mathbf{W}\mathbf{D}^{-2}(b-\mathbf{A}\hat{x})$

In order to have that $(x_{k+1}-\hat{x})$ goes to zero in the limit we should
require that this limiting error update has the zero vector as a fixed point,
i.e.

$0=\mathbf{A}^{T}\mathbf{P}\mathbf{W}\mathbf{D}^{-2}(b-\mathbf{A}\hat{x})$

This is guaranteed if $\mathbf{P}\mathbf{W}\mathbf{D}^{-2}=\beta\mathbf{I}$.
But for $\lambda\neq 0$, we do not have a bound of the form
$D_{f}^{x_{k}^{*}}(x_{k},\hat{x})\leq\gamma\cdot\|Ax_{k}-b\|^{2}_{\mathbf{P}\mathbf{W}\mathbf{D}^{-2}}$
that is the reason why for that case we need to assume
$\mathbf{P}\mathbf{W}\mathbf{D}^{-2}=\beta\mathbf{I}$.

Since in this case
$\mathbf{P}\mathbf{W}^{2}\mathbf{D}^{-2}=\mathbf{P}\mathbf{W}\mathbf{D}^{-2}\mathbf{W}$,
Lemma 1 becomes:

######  Lemma 4.

Let $\mathbf{M}_{k},\mathbf{P},\mathbf{W}$ and $\mathbf{D}$ be defined as in
Definition 2 and let Assumption 1 hold. Then

$\mathbb{E}_{k}\left[\mathbf{M}_{k}\right]=\frac{\alpha\mathbf{I}}{\|\mathbf{A}\|^{2}_{F}}$

and

$\mathbb{E}_{k}\left[\mathbf{M_{k}^{T}AA^{T}M_{k}}\right]=\tfrac{1}{\eta}\frac{\alpha\mathbf{W}}{\|\mathbf{A}\|^{2}_{F}}+\alpha^{2}(1-\tfrac{1}{\eta})\frac{\mathbf{AA^{T}}}{\|\mathbf{A}\|^{4}_{F}}$

######  Lemma 5.

Under Assumption 1, for the iterates $x_{k}$ of Algorithm 1, it holds that:

$\displaystyle\mathbb{E}_{k}\left[D_{f}^{x_{k+1}^{*}}(x_{k+1},\hat{x})\right]\leq
D_{f}^{x_{k}^{*}}(x_{k},\hat{x})-\frac{\alpha}{\|\mathbf{A}\|^{2}_{F}}(1-\sigma_{\max}(\mathbf{T}))\|\mathbf{A}x_{k}-b\|^{2}_{2}.$

with

$\mathbf{T}=\tfrac{1}{2\eta}\mathbf{W}+\tfrac{\alpha}{2}(1-\tfrac{1}{\eta})\frac{\mathbf{A}\mathbf{A}^{T}}{\|\mathbf{A}\|^{2}_{F}}$

###### Proof.

Using Lemma 2 with $f(x)=\lambda\|x\|_{1}+\frac{1}{2}\|x\|_{2}^{2}$ and
$\varPhi(x)=\langle\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b),x-x_{k}\rangle$,
$y=\hat{x}$, it holds that:

$\displaystyle D^{x_{k+1}^{*}}_{f}(x_{k+1},\hat{x})$ $\displaystyle\leq
D^{x_{k}^{*}}_{f}(x_{k},\hat{x})+\varPhi(\hat{x})-\varPhi(x_{k+1})-D^{x_{k}^{*}}_{f}(x_{k},x_{k+1})$
$\displaystyle=D^{x_{k}^{*}}_{f}(x_{k},\hat{x})-\langle\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b),x_{k}-\hat{x}\rangle+\langle\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b),x_{k}-x_{k+1}\rangle-D^{x_{k}^{*}}_{f}(x_{k},x_{k+1})$
$\displaystyle\leq
D^{x_{k}^{*}}_{f}(x_{k},\hat{x})-\langle\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b),x_{k}-\hat{x}\rangle+\|\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b)\|\cdot\|x_{k}-x_{k+1}\|-\tfrac{1}{2}\|x_{k}-x_{k+1}\|^{2}$
$\displaystyle\leq
D^{x_{k}^{*}}_{f}(x_{k},\hat{x})-\langle\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b),x_{k}-\hat{x}\rangle+\tfrac{1}{2}\|\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b)\|^{2}$

We have:

$\displaystyle\mathbb{E}_{k}\left[\Bigl{\langle}\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b),x_{k}-\hat{x}\Bigr{\rangle}\right]$
$\displaystyle=\mathbb{E}_{k}\left[\Bigl{\langle}\mathbf{M}_{k}\cdot(\mathbf{A}x_{k}-b),\mathbf{A}x_{k}-b\Bigr{\rangle}\right]$
$\displaystyle\overset{\text{\tiny Lemma\leavevmode\nobreak\
\ref{lemma1}}}{=}\Bigl{\langle}\mathbf{A}x_{k}-b,\mathbf{A}x_{k}-b\Bigr{\rangle}_{\mathbf{PWD}^{-2}}\overset{\text{\tiny
Lemma\leavevmode\nobreak\
\ref{lemma2}}}{=}\frac{\alpha}{\|\mathbf{A}\|^{2}_{F}}\|\mathbf{A}x_{k}-b\|^{2}_{2}$

and with
$\mathbf{T}=\tfrac{1}{2\eta}\mathbf{W}+\tfrac{\alpha}{2}(1-\tfrac{1}{\eta})\frac{\mathbf{A}\mathbf{A}^{T}}{\|\mathbf{A}\|^{2}_{F}}$
we get

$\displaystyle\footnotesize\mathbb{E}_{k}\left[\Bigl{\|}\mathbf{A}^{T}\cdot\mathbf{M}_{k}\cdot(\mathbf{A}x_{k}-b)\Bigr{\|}_{2}^{2}\right]$
$\displaystyle=\Bigl{\langle}\mathbf{A}x_{k}-b,\mathbb{E}_{k}\left[\mathbf{M}_{k}^{T}\mathbf{A}\mathbf{A}^{T}\mathbf{M}_{k}\right]\cdot(\mathbf{A}x_{k}-b)\Bigr{\rangle}$
$\displaystyle\overset{\text{\tiny Lemma\leavevmode\nobreak\
\ref{lemma2}}}{=}\Bigl{\langle}\mathbf{A}x_{k}-b,\Bigg{(}\tfrac{1}{\eta}\frac{\alpha\mathbf{W}}{\|\mathbf{A}\|^{2}_{F}}+\alpha^{2}(1-\tfrac{1}{\eta})\frac{\mathbf{A}\mathbf{A}^{T}}{\|\mathbf{A}\|^{4}_{F}}\Bigg{)}\cdot(\mathbf{A}x_{k}-b)\Bigr{\rangle}$
$\displaystyle=\tfrac{2\alpha}{\|\mathbf{A}\|^{2}_{F}}\Bigl{\langle}\mathbf{A}x_{k}-b,\mathbf{A}x_{k}-b\Bigr{\rangle}_{\mathbf{T}}$
$\displaystyle\leq\tfrac{2\alpha}{\|\mathbf{A}\|^{2}_{F}}\sigma_{\max}(\mathbf{T})\|\mathbf{A}x_{k}-b\|^{2}_{2}.$

Thus combining everything together gives us

$\displaystyle\mathbb{E}_{k}\left[D_{f}^{x_{k+1}^{*}}(x_{k+1},\hat{x})\right]$
$\displaystyle\leq
D_{f}^{x_{k}^{*}}(x_{k},\hat{x})-\mathbb{E}_{k}\left[\Bigl{\langle}\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b),x_{k}-\hat{x}\Bigr{\rangle}\right]+\frac{1}{2}\mathbb{E}_{k}\left[\Bigl{\|}\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b)\Bigr{\|}_{2}^{2}\right]$
$\displaystyle\leq
D_{f}^{x_{k}^{*}}(x_{k},\hat{x})-\frac{\alpha}{\|\mathbf{A}\|^{2}_{F}}\|\mathbf{A}x_{k}-b\|^{2}_{2}+\tfrac{1}{2}\tfrac{2\alpha}{\|\mathbf{A}\|^{2}_{F}}\sigma_{\max}(\mathbf{T})\|\mathbf{A}x_{k}-b\|^{2}_{2}$
$\displaystyle\leq
D_{f}^{x_{k}^{*}}(x_{k},\hat{x})-\frac{\alpha}{\|\mathbf{A}\|^{2}_{F}}(1-\sigma_{\max}(\mathbf{T}))\|\mathbf{A}x_{k}-b\|^{2}_{2}.$

∎

The following lemma gives an upper and a lower bound for the largest singular
value of $\mathbf{T}$ which will be use in the convergence of RSKA iterates.

######  Lemma 6.

Let

$\mathbf{T}=\tfrac{1}{2\eta}\mathbf{W}+\tfrac{\alpha}{2}(1-\tfrac{1}{\eta})\frac{\mathbf{A}\mathbf{A}^{T}}{\|\mathbf{A}\|^{2}_{F}}$

Then the largest singular value of $\mathbf{T}$ satisfies:

$\displaystyle\tfrac{1}{2\eta}\bigg{(}\sigma_{\max}(\mathbf{W})-\tfrac{\alpha}{\|\mathbf{A}\|^{2}_{F}}(\eta-1)\sigma_{\max}^{2}(\mathbf{A})\bigg{)}$
$\displaystyle\leq\sigma_{\max}(\mathbf{T})$
$\displaystyle\leq\tfrac{1}{2\eta}\bigg{(}\sigma_{\max}(\mathbf{W})+\tfrac{\alpha}{\|\mathbf{A}\|^{2}_{F}}(\eta-1)\sigma_{\max}^{2}(\mathbf{A})\bigg{)}$

In addition, If $\mathbf{W}=\alpha\mathbf{I}$, then $\mathbf{T}$ is positive
semi-definite and

$\sigma_{\max}(\mathbf{T})=\tfrac{1}{2\eta}\bigg{(}\alpha+\tfrac{\alpha}{\|\mathbf{A}\|^{2}_{F}}(\eta-1)\sigma_{\max}^{2}(\mathbf{A})\bigg{)}$

###### Proof.

The first part of the proof follows easily from Theorem 2.1. If
$\mathbf{W}=\alpha\mathbf{I}$, we have: If $\lambda$ is an eigenvalue of
$\mathbf{A}\mathbf{A}^{T}$ then
$\tfrac{\alpha}{\eta}+\tfrac{\alpha}{2}(1-\tfrac{1}{\eta})\tfrac{\lambda}{\|\mathbf{A}\|_{F}^{2}}$
is an eigenvalue of $\mathbf{T}$. From this we deduce the last equality as
well as that $\mathbf{T}$ is positive semi-definite. ∎

###  4.1 General Convergence Result

In this part, we present general convergence results for the iterates from
RSKA method.

######  Theorem 4.1 (Noiseless case).

Consider $\eta>1$, $\gamma$ as defined in equation (27) (Lemma (3)), let
Assumption 1 hold and assume that

(28) $\displaystyle
0<\alpha<2\frac{(\eta-\tfrac{1}{2}\sigma_{\max}(\mathbf{W}))\|\mathbf{A}\|_{F}^{2}}{\sigma_{\max}(\mathbf{A})^{2}(\eta-1)}.$

Then the random iterates $x_{k}$ produced by Algorithm 1 converge in
expectation with a linear rate to the unique solution $\hat{x}$ of
$\operatorname*{\min}_{\mathbf{A}x=b}\lambda\|x\|_{1}+\frac{1}{2}\|x\|_{2}^{2}$,
more precisely, with

(29) $\displaystyle
q=1-\frac{1}{\gamma}\cdot\frac{L(\alpha)}{\|\mathbf{A}\|^{2}_{F}}\,\in(0,1),$

and

$\displaystyle
L(\alpha)=\alpha-\tfrac{\alpha}{2\eta}\bigg{(}\tfrac{\alpha}{\|\mathbf{A}\|^{2}_{F}}(\eta-1)\sigma^{2}_{\max}(\mathbf{A})+\sigma_{\max}(\mathbf{W})\bigg{)},$

it holds that

$\displaystyle\mathbb{E}\left[D_{f}^{x_{k+1}^{*}}(x_{k+1},\hat{x})\right]$
$\displaystyle\leq
q\cdot\mathbb{E}\left[D_{f}^{x_{k}^{*}}(x_{k},\hat{x})\right]\,$
$\displaystyle\mathbb{E}\bigg{[}\|x_{k}-\hat{x}\|_{2}^{2}\bigg{]}$
$\displaystyle\leq 2\cdot q^{k}\cdot f(\hat{x}).$

###### Proof.

Combining Lemma 5 with equation (26) gives

$\displaystyle\mathbb{E}_{k}\left[D_{f}^{x_{k+1}^{*}}(x_{k+1},\hat{x})\right]$
$\displaystyle\leq
D_{f}^{x_{k}^{*}}(x_{k},\hat{x})-\tfrac{\alpha}{\|\mathbf{A}\|_{F}^{2}\gamma}(1-\sigma_{\max}(\mathbf{T}))D_{f}^{x_{k}^{*}}(x_{k},\hat{x})$
$\displaystyle\leq\left(1-\tfrac{\alpha}{\|\mathbf{A}\|_{F}^{2}\gamma}(1-\sigma_{\max}(\mathbf{T}))\right)D_{f}^{x_{k}^{*}}(x_{k},\hat{x}).$

Using the rule of total expectation we get

$\displaystyle\mathbb{E}\left[D_{f}^{x_{k+1}^{*}}(x_{k+1},\hat{x})\right]$
$\displaystyle\leq\left(1-\tfrac{1}{\gamma}\cdot\tfrac{\alpha\cdot(1-\sigma_{\max}(\mathbf{T}))}{\|\mathbf{A}\|_{F}^{2}}\right)\mathbb{E}\left[D_{f}^{x_{k}^{*}}(x_{k},\hat{x})\right].$

To get a rate $q\in(0,1)$, we need that $(1-\sigma_{\max}(\mathbf{T}))>0$,
i.e. $\sigma_{\max}(\mathbf{T})<1$ which holds true from (28). From Lemma 6,
it follows that $L(\alpha)\leq\alpha(1-\sigma_{\max}(\mathbf{T}))$ and thus we
get $\mathbb{E}\left[D_{f}^{x_{k+1}^{*}}(x_{k+1},\hat{x})\right]\leq
q\cdot\mathbb{E}\left[D_{f}^{x_{k}^{*}}(x_{k},\hat{x})\right]$, with
$q=1-\frac{1}{\gamma}\cdot\frac{L(\alpha)}{\|\mathbf{A}\|^{2}_{F}}$. The
inequality in terms of $\|x_{k}-\hat{x}\|_{2}^{2}$ is obtained by using the
first inequality of equation (13) since $f$ is $1$-strongly convex. ∎

From (29) we see that we want to choose $\alpha$ such that $L(\alpha)$ is as
large as possible:

######  Corollary 4.2.

Let Assumption 1 hold true. Then the relaxation parameter $\alpha$ and the
constant $L$ which yields the fastest convergence rate guarantee in Theorem
4.1 are as follows:

  1. (a)

General Weights: If $\eta>\max(1,\sigma_{\max}(\mathbf{W})/2)$ then

$\displaystyle\alpha^{*}=\frac{\|\mathbf{A}\|^{2}_{F}}{\sigma^{2}_{\max}(\mathbf{A})(\eta-1)}(\eta-\tfrac{\sigma_{\max}(\mathbf{W})}{2}),$

and

$\displaystyle
L(\alpha^{*})=\frac{\|\mathbf{A}\|^{2}_{F}}{8\sigma^{2}_{max}(\mathbf{A})}\cdot\frac{(2\eta-\sigma_{max}(\mathbf{W}))^{2}}{\eta(\eta-1)}.$

  2. (b)

Uniform Weights i.e $\mathbf{W}=\alpha\mathbf{I}$:

$\alpha^{*}=\frac{\eta}{1+(\eta-1)\tfrac{\sigma^{2}_{\max}(\mathbf{A})}{\|\mathbf{A}\|^{2}_{F}}}$

and

$L(\alpha^{*})=\frac{\eta}{2+2(\eta-1)\tfrac{\sigma^{2}_{\max}(\mathbf{A})}{\|\mathbf{A}\|^{2}_{F}}}.$

###### Proof.

In the case (a) of general weights, in order to get the tightest lower bound,
we maximized the concave function $L(\alpha)$ and obtain

$\displaystyle\alpha=\|\mathbf{A}\|_{F}^{2}(\eta-\sigma_{\max}(\mathbf{W})/2)/((\eta-1)\sigma_{\max}(\mathbf{A})^{2}$

which fulfills (28) and thus gives the best $q$ is Theorem 4.1. Since we need
$\alpha\geq 0$ we have to assume $\eta\geq\sigma_{\max}(\mathbf{W})/2$. This
gives $\alpha^{*}$ and plugging this into $L(\alpha)$ give us $L(\alpha^{*})$.
In the case (b) of uniform weights, we maximized $L(\alpha)$ for all $\eta$
with $\mathbf{W}=\alpha\mathbf{I}$. ∎

A few remarks about the interpretation of the above theorem are in order:

######  Remark 1 (Overrelaxation).

  * •

When a single thread $\eta=1$ is used in the case (b) of uniform weight, we
see that our optimal relaxation parameter is $\alpha^{*}=1$. Whereas, when
multiple threads $\eta>1$ are used, we see

$1<\alpha^{*}\leq\eta,$

i.e. the method allows for large over-relaxation when the number of threads is
high and we will see in Section 5.2 that this does indeed lead to faster
convergence.

  * •

Our relaxation parameter $\alpha^{*}$ in the case of uniform weights is the
same as the relaxation parameter $\alpha^{RT}$ suggested in
richtarik2020stochastic  although they do not treat the sparse case and only
consider uniform weights.

  * •

Finally, for $\eta=1$ in the case of general weights we have, due to the
coupling in Assumption 1, that the weights fulfill
$w_{i}=\tfrac{\alpha\|a_{i}\|^{2}}{p_{i}\|\mathbf{A}\|_{F}^{2}}$. We could
estimate $\sigma_{\max}(\mathbf{W})\leq\tfrac{\alpha}{\min_{i}p_{i}}$ but this
would be a quite crude estimate. By choosing the classical probabilities
$p_{i}=\|a_{i}\|^{2}/\|\mathbf{A}\|_{F}^{2}$ we would get
$\sigma_{\max}(\mathbf{W})=\alpha$ and get that the relaxation parameter need
to fulfill $\alpha\in(0,2)$ and that $\alpha^{*}=1$ is the optimal relaxation
parameter.

######  Remark 2 (Relation to standard randomized sparse Kaczmarz).

In the case $\mathbf{W}=\mathbf{I},\eta=1,\alpha=1$, we have
$\mathbf{T}=\tfrac{1}{2}\mathbf{I}$ and we recover the rate of the standard
RSK. It holds

$\displaystyle\mathbb{E}_{k}\left[D_{f}^{x_{k+1}^{*}}(x_{k+1},\hat{x})\right]\leq
D_{f}^{x_{k}^{*}}(x_{k},\hat{x})-\frac{1}{2\|\mathbf{A}\|^{2}_{F}}\|\mathbf{A}x_{k}-b\|^{2}_{2}$

which is obtained in SL19  and it is shown that this leads to a linear
convergence rate in expectation,

(30) $\mathbb{E}\bigg{[}\|x_{k}-\hat{x}\|_{2}^{2}\bigg{]}\leq
2\cdot\bigg{(}1-\frac{1}{2\gamma\|\mathbf{A}\|^{2}_{F}}\bigg{)}^{k}\cdot
f(\hat{x}).$

which implies that we reach accuracy
$\mathbb{E}\bigg{[}\|x_{k}-\hat{x}\|_{2}^{2}\bigg{]}\leq 2\cdot\epsilon\cdot
f(\hat{x})$ in at most $k\geq 2\gamma\|A\|^{2}_{F}\log(\tfrac{1}{\epsilon})$
iterations.

######  Remark 3 (Convergence rate for the linearized Bregman method).

Similarly to Lemma 5 we can show that the linearized Bregman algorithm
cai2009linearized ; yin2010analysis ; LSW14

(31) $\displaystyle x_{k+1}^{*}$
$\displaystyle=x_{k}^{*}-\frac{A^{T}(Ax_{k}-b)}{\|A\|^{2}_{2}},$
$\displaystyle x_{k+1}$ $\displaystyle=S_{\lambda}(x_{k+1}^{*})$

has linear convergence rate given by

(32) $\big{\|}x_{k}-\hat{x}\big{\|}_{2}^{2}\leq
2\cdot\bigg{(}1-\frac{1}{2\gamma\|\mathbf{A}\|^{2}_{2}}\bigg{)}^{k}\cdot
f(\hat{x}).$

Although we suspect that this result (32) is not new we could not find it in
the literature.

Inspired by richtarik2020stochastic , the following remarks apply to the
uniform weight case.

######  Remark 4 (Mini-batch vs. full-batch).

Let
$H(\eta)=\frac{1}{L(\alpha^{*})}=\frac{2}{\eta}+2(1-\frac{1}{\eta})\tfrac{\sigma^{2}_{\max}(\mathbf{A})}{\|\mathbf{A}\|^{2}_{F}}$
be the inverse of $L(\alpha^{*})$. Recall from Theorem 4.1 that
$L(\alpha^{*})$ influences the convergence rate: The larger $L(\alpha^{*})$,
the faster the convergence.

Since $\frac{\|\mathbf{A}\|^{2}_{F}}{\sigma^{2}_{\max}(\mathbf{A})}\geq 1$,
$H$ is a nonincreasing function of $\eta$ and we have $H(1)=2$,
$H(\infty)\overset{\text{def}}{=}\lim_{\eta\rightarrow\infty}H(\eta)=\tfrac{2\sigma^{2}_{\max}(\mathbf{A})}{\|\mathbf{A}\|^{2}_{F}}$.
In the asymptotic regime $\eta\rightarrow\infty,$ Algorithm 1 becomes
linearized Bregman algorithm for minimizing (15), and $H(\infty)$ is the rate
of linearized Bregman cf. (32). This shows that the averaging method
interpolates between the basic method and the linearized Bregman. By
increasing $\eta$, the quantity
$\frac{H(1)}{H(\infty)}=\frac{\|\mathbf{A}\|^{2}_{F}}{\sigma^{2}_{\max}(\mathbf{A})}$
controls the maximum (guaranteed) speedup in the iteration complexity
achievable. Comparing (29) with the convergence rate (30) of the basic sparse
Kaczmarz method, we get an improvement of $2L(\alpha^{*})>1$ which shows that
for the RSKA algorithm, we can get a speed-up even of order approximately
$\frac{\|\mathbf{A}\|^{2}_{F}}{\sigma^{2}_{\max}(\mathbf{A})}$ compared to the
rate of the basic sparse Kaczmarz algorithm (the also the comparison in Table
1).

For $\eta\geq\frac{\|\mathbf{A}\|^{2}_{F}}{\sigma^{2}_{\max}(\mathbf{A})},$ we
get $H(\eta)\leq 2H(\infty)$ , which is the performance of the full batch (up
to a factor of 2). This means that it does not make sense to use a minibatch
size larger than
$\frac{\|\mathbf{A}\|^{2}_{F}}{\sigma^{2}_{\max}(\mathbf{A})}$. Moreover,
notice that $H(\eta)\geq\frac{1}{\eta}H(1)$ for all $\eta$, show that the
number of iterations does not decrease linearly in the minibatch size $\eta$.
From a total complexity perspective cf. Table 1, this also means that in a
computational regime where processing $\eta$ basic method updates costs $\eta$
times as much as processing a single update, the decrease in iteration
complexity cannot compensate for the increase in cost per iteration, which
means that the choice $\eta=1$ is optimal. On the other hand, if a parallel
processor is available, a larger $\eta$ will be better.

| RSK | RSKA | linBreg
---|---|---|---
iteration complexity | $\mathcal{O}\bigg{(}2\gamma\|\mathbf{A}\|^{2}_{F}\log(\tfrac{1}{\epsilon})\bigg{)}$ |  $\mathcal{O}\bigg{(}\gamma\frac{\|\mathbf{A}\|^{2}_{F}}{L(\alpha^{*})}\log(\tfrac{1}{\epsilon})\bigg{)}$ | $\mathcal{O}\bigg{(}2\gamma\|\mathbf{A}\|^{2}_{2}\log(\tfrac{1}{\epsilon})\bigg{)}$
cost per iteration | $\mathcal{O}(n)$ |  $\mathcal{O}(\eta n)$ | $\mathcal{O}(mn)$
Table 1: Complexity of different methods.

###  4.2 Noisy right hand sides

In the noisy case we consider a consistent linear system $\mathbf{A}x=b$, but
assume that instead of $b$ we only have access to some vector $b^{\delta}$
with $\|b-b^{\delta}\|_{2}\leq\delta$. In the context if Kaczmarz methods,
such errors in the right hand side have been considered in
needell2010randomized ; zouzias2013randomized . Since the system
$\mathbf{A}x=b^{\delta}$ is most likely inconsistent, RSKA will not solve the
optimization problem (5) but hopefully the iterates still come close to the
solution. If we assume a bound $\|b-b^{\delta}\|_{2}\leq\delta$ for $b$ with
$\mathbf{A}\hat{x}=b$ we get the following result on the convergence of the
method:

######  Theorem 4.3 (Noisy case).

Assume that instead of exact data $b\in\mathbf{Range}(\mathbf{A})$ only a
noisy right hand side $b^{\delta}\in\mathbb{R}^{m}$ with
$\|b^{\delta}-b\|_{2}\leq\delta$ is given. Consider $\eta>1$, $\epsilon>0$,
$\gamma$ as defined in equation (27) (Lemma (3)), let Assumption 1 hold and
assume that

(33) $\displaystyle
0<\alpha<2\frac{((1-\epsilon)\eta-\tfrac{1}{2}\sigma_{\max}(\mathbf{W}))\|A\|_{F}^{2}}{\sigma_{\max}(\mathbf{A})^{2}(\eta-1)}.$

If the iterates $x_{k}$ of the RSKA method from Algorithm 1 are computed with
$b$ replaced by $b^{\delta}$, then, with the contraction factor $a$, we have :

(34) $\displaystyle
a=1-\frac{\alpha}{\gamma}\cdot\frac{(1-\epsilon-\sigma_{\max}(\mathbf{T}))}{\|\mathbf{A}\|^{2}_{F}}\,\in(0,1),$

and the expected rate of convergence is

$\displaystyle\mathbb{E}\left[D_{f}^{x_{k+1}^{*}}(x_{k+1},\hat{x})\right]$
$\displaystyle\leq
a\cdot\mathbb{E}\left[D_{f}^{x_{k}^{*}}(x_{k},\hat{x})\right]+c\delta^{2}\,$
$\displaystyle\mathbb{E}\bigg{[}\|x_{k}-\hat{x}\|_{2}^{2}\bigg{]}$
$\displaystyle\leq 2\cdot a^{k}\cdot f(\hat{x})+\frac{c}{1-a}\delta^{2}.$

where
$c=\tfrac{\alpha}{\|A\|^{2}_{F}}\bigg{(}\sigma_{\max}(\mathbf{T})+\tfrac{1}{\epsilon}\sigma_{\max}^{2}(\mathbf{T}^{{}^{\prime}})\bigg{)},\quad\mathbf{T}^{{}^{\prime}}=\mathbf{T}-\tfrac{1}{2}\mathbf{I}$

###### Proof.

Assuming that a noisy observed data $b^{\delta}\in\mathbb{R}^{m}$ instead of
$b$ with $\|b^{\delta}-b\|_{2}\leq\delta$ is given, where $b=A\hat{x}$. The
update in this case is given by:

(35) $\displaystyle x^{*}_{k+1}$
$\displaystyle=x^{*}_{k}-\tfrac{1}{\eta}\sum_{i\in\tau_{k}}w_{i}\tfrac{\langle
a_{i},x_{k}\rangle-b_{i}^{\delta}}{\|a_{i}\|^{2}_{2}}\cdot a_{i},$
$\displaystyle x_{k+1}$ $\displaystyle=S_{\lambda}(x^{*}_{k+1})$

where $\eta=|\tau_{k}|$,which in terms of matrix multiplication is equal to:

(36) $\displaystyle x^{*}_{k+1}$
$\displaystyle=x^{*}_{k}-\mathbf{A}^{T}\cdot\mathbf{M}_{k}\cdot(\mathbf{A}x_{k}-b^{\delta}),$
$\displaystyle x_{k+1}$ $\displaystyle=S_{\lambda}(x^{*}_{k+1})$

We introduce the abbreviation

$x_{k}^{\delta}\overset{\text{def}}{=}\hat{x}+\mathbf{A}^{T}\cdot\mathbf{M}_{k}\cdot(b-b^{\delta})$

and use Lemma 2 with $f(x)=\lambda\|x\|_{1}+\frac{1}{2}\|x\|_{2}^{2}$ and
$\varPhi(x)=\langle\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b^{\delta}),x-x_{k}\rangle$,
and $y=x_{k}^{\delta}$ and get

$\displaystyle D^{x_{k+1}^{*}}_{f}(x_{k+1},x_{k}^{\delta})$ $\displaystyle\leq
D^{x_{k}^{*}}_{f}(x_{k},x_{k}^{\delta})+\varPhi(x_{k}^{\delta})-\varPhi(x_{k+1})-D^{x_{k}^{*}}_{f}(x_{k},x_{k+1})$
$\displaystyle=D^{x_{k}^{*}}_{f}(x_{k},x_{k}^{\delta})-\langle\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b^{\delta}),x_{k}-x_{k}^{\delta}\rangle$
$\displaystyle\qquad+\langle\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b^{\delta}),x_{k}-x_{k+1}\rangle-D^{x_{k}^{*}}_{f}(x_{k},x_{k+1})$
$\displaystyle\leq
D^{x_{k}^{*}}_{f}(x_{k},x_{k}^{\delta})-\langle\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b^{\delta}),x_{k}-x_{k}^{\delta}\rangle$
$\displaystyle\qquad+\|\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b^{\delta})\|\cdot\|x_{k}-x_{k+1}\|-\tfrac{1}{2}\|x_{k}-x_{k+1}\|^{2}$
$\displaystyle\leq
D^{x_{k}^{*}}_{f}(x_{k},x_{k}^{\delta})-\langle\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b^{\delta}),x_{k}-x_{k}^{\delta}\rangle+\tfrac{1}{2}\|\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b^{\delta})\|^{2}\vspace{-0.8cm}$

so that

(37) $\displaystyle D_{f}^{x_{k+1}^{*}}(x_{k+1},x_{k}^{\delta})$
$\displaystyle\leq
D_{f}^{x_{k}^{*}}(x_{k},x_{k}^{\delta})-\langle\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b^{\delta}),x_{k}-x_{k}^{\delta}\rangle+\frac{1}{2}\Bigl{\|}\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b^{\delta})\Bigr{\|}_{2}^{2}\vspace{-0.8cm}$

Unfolding the expression of $D_{f}^{x_{k+1}^{*}}(x_{k+1},x_{k}^{\delta})$ and
$D_{f}^{x_{k}^{*}}(x_{k},x_{k}^{\delta}),$ we get:

$\displaystyle\vspace{-0.8cm}D_{f}^{x_{k+1}^{*}}(x_{k+1},\hat{x})$
$\displaystyle\leq
D_{f}^{x_{k}^{*}}(x_{k},\hat{x})-\langle\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b^{\delta}),x_{k}-\hat{x}\rangle+\frac{1}{2}\Bigl{\|}\mathbf{A}^{T}\mathbf{M}_{k}(\mathbf{A}x_{k}-b^{\delta})\Bigr{\|}_{2}^{2}\vspace{-0.8cm}$

On the other hand,

$\displaystyle\mathbb{E}_{k}\left[\langle\mathbf{M}_{k}\cdot(\mathbf{A}x_{k}-b^{\delta}),\mathbf{A}x_{k}-b\rangle\right]$
$\displaystyle=\frac{\alpha}{\|\mathbf{A}\|^{2}_{F}}\|\mathbf{A}x_{k}-b\|^{2}_{2}+\frac{\alpha}{\|\mathbf{A}\|^{2}_{F}}\langle
b-b^{\delta},\mathbf{A}x_{k}-b\rangle\vspace{-0.8cm}$

and

$\displaystyle\mathbb{E}_{k}\left[\Bigl{\|}\mathbf{A}^{T}\cdot\mathbf{M}_{k}\cdot(\mathbf{A}x_{k}-b^{\delta})\Bigr{\|}_{2}^{2}\right]$
$\displaystyle=\langle\mathbf{A}x_{k}-b,\mathbb{E}_{k}\left[\mathbf{M}_{k}^{T}\mathbf{A}\mathbf{A}^{T}\mathbf{M}_{k}\right]\cdot(\mathbf{A}x_{k}-b)\rangle$
$\displaystyle+\langle
b-b^{\delta},\mathbb{E}_{k}\left[\mathbf{M}_{k}^{T}\mathbf{A}\mathbf{A}^{T}\mathbf{M}_{k}\right]\cdot(b-b^{\delta})\rangle$
$\displaystyle+2\langle\mathbf{A}x_{k}-b,\mathbb{E}_{k}\left[\mathbf{M}_{k}^{T}\mathbf{A}\mathbf{A}^{T}\mathbf{M}_{k}\right]\cdot(b-b^{\delta})\rangle$

In summary we have:

(38)
$\displaystyle\begin{split}&-\mathbb{E}_{k}\left[\langle\mathbf{A}^{T}\cdot\mathbf{M}_{k}\cdot(\mathbf{A}x_{k}-b^{\delta}),x_{k}-\hat{x}\rangle\right]+\frac{1}{2}\mathbb{E}_{k}\left[\Bigl{\|}\mathbf{A}^{T}\cdot\mathbf{M}_{k}\cdot(\mathbf{A}x_{k}-b^{\delta})\Bigr{\|}_{2}^{2}\right]\\\
&=-\frac{\alpha}{\|\mathbf{A}\|^{2}_{F}}\|\mathbf{A}x_{k}-b\|^{2}_{2}-\frac{\alpha}{\|\mathbf{A}\|^{2}_{F}}\langle
b-b^{\delta},\mathbf{A}x_{k}-b\rangle+\tfrac{1}{2}\langle
b-b^{\delta},\mathbb{E}_{k}\left[\mathbf{M}_{k}^{T}\mathbf{A}\mathbf{A}^{T}\mathbf{M}_{k}\right](b-b^{\delta})\rangle\\\
&+\tfrac{1}{2}\langle\mathbf{A}x_{k}-b,\mathbb{E}_{k}\left[\mathbf{M}_{k}^{T}\mathbf{A}\mathbf{A}^{T}\mathbf{M}_{k}\right](\mathbf{A}x_{k}-b)\rangle+\langle\mathbf{A}x_{k}-b,\mathbb{E}_{k}\left[\mathbf{M}_{k}^{T}\mathbf{A}\mathbf{A}^{T}\mathbf{M}_{k}\right]\cdot(b-b^{\delta})\rangle\\\
&\leq-\frac{\alpha}{\|\mathbf{A}\|^{2}_{F}}\|\mathbf{A}x_{k}-b\|^{2}_{2}+\tfrac{\alpha}{\|\mathbf{A}\|^{2}_{F}}\sigma_{\max}(\mathbf{T})\|\mathbf{A}x_{k}-b\|^{2}_{2}+\tfrac{\alpha}{\|\mathbf{A}\|^{2}_{F}}\sigma_{\max}(\mathbf{T})\|b-b^{\delta}\|^{2}_{2}\\\
&+\tfrac{2\alpha}{\|\mathbf{A}\|^{2}_{F}}\langle\mathbf{A}x_{k}-b,\Bigg{(}\tfrac{1}{2\eta}W-\tfrac{1}{2}I+\tfrac{\alpha}{2}(1-\tfrac{1}{\eta})\frac{AA^{T}}{\|\mathbf{A}\|^{2}_{F}}\Bigg{)}\cdot(b-b^{\delta})\rangle\\\
&\leq-\frac{\alpha}{\|\mathbf{A}\|^{2}_{F}}\|\mathbf{A}x_{k}-b\|^{2}_{2}+\tfrac{\alpha}{\|\mathbf{A}\|^{2}_{F}}\sigma_{\max}(\mathbf{T})\|\mathbf{A}x_{k}-b\|^{2}_{2}+\tfrac{\alpha}{\|\mathbf{A}\|^{2}_{F}}\sigma_{\max}(\mathbf{T})\|b-b^{\delta}\|^{2}_{2}\\\
&+\tfrac{\alpha\epsilon}{\|\mathbf{A}\|^{2}_{F}}\|\mathbf{A}x_{k}-b\|^{2}_{2}+\tfrac{\alpha}{\epsilon\|\mathbf{A}\|^{2}_{F}}\Bigg{\|}\Bigg{(}\tfrac{1}{2\eta}W-\tfrac{1}{2}I+\tfrac{\alpha}{2}(1-\tfrac{1}{\eta})\frac{\mathbf{A}\mathbf{A}^{T}}{\|\mathbf{A}\|^{2}_{F}}\Bigg{)}\cdot(b-b^{\delta})\Bigg{\|}^{2}_{2}\\\
&=-\frac{\alpha}{\|\mathbf{A}\|^{2}_{F}}(1-\epsilon-\sigma_{\max}(\mathbf{T}))\|\mathbf{A}x_{k}-b\|^{2}_{2}+\frac{\alpha}{\|\mathbf{A}\|^{2}_{F}}(\sigma_{\max}(\mathbf{T})+\frac{1}{\epsilon}\sigma_{\max}^{2}(\mathbf{T}^{{}^{\prime}}))\|b-b^{\delta}\|^{2}_{2}.\end{split}$

To get an improvement after each iteration, we need
$(1-\epsilon-\sigma_{\max}(\mathbf{T}))>0$, i.e.
$\sigma_{\max}(\mathbf{T})<1-\epsilon$ which hold true because of (33).
Combining (37) and (38) and using the error bound from Lemma 3 we arrive at

$\displaystyle\mathbb{E}_{k}\left[D_{f}^{x_{k+1}^{*}}(x_{k+1},\hat{x})\right]$
$\displaystyle\leq\left(1-\tfrac{\alpha}{\gamma\|\mathbf{A}\|_{F}^{2}}(1-\epsilon-\sigma_{\max}(\mathbf{T}))\right)D_{f}^{x_{k}^{*}}(x_{k},\hat{x})$
$\displaystyle\qquad+\tfrac{\alpha}{\|\mathbf{A}\|^{2}_{F}}(\sigma_{\max}(\mathbf{T})+\tfrac{1}{\epsilon}\sigma^{2}_{\max}(\mathbf{T}^{\prime}))\delta^{2},$

which concludes the proof. The inequality in terms of the norm is obtained by
using the first inequality of (13) since $f$ is $1$-strongly convex. ∎

######  Remark 5.

Note that for $\mathbf{W}=\mathbf{I},\eta=1,\alpha=1$, we have
$\mathbf{T}=\tfrac{1}{2}\mathbf{I}$ meaning $\mathbf{T}^{{}^{\prime}}=0$ and
sending $\epsilon$ to zero give us $\frac{c}{1-a}=\gamma$ which recover the
rate of the standard RSK in the noisy case showed in SL19 .

##  5 Numerical Experiments

We present several experiments to demonstrate the effectiveness of Algorithm 1
under various conditions. In particular, we study the effects of the
relaxation parameter $\alpha$, the number of threads $\eta$, the sparsity
parameter $\lambda$, the weight matrix $\mathbf{W}$, and the probability
matrix $\mathbf{P}$. The simulations were performed in Python on an Intel Core
i7 computer with 16GB RAM. We start by comparing several variants of RSKA
algorithms with randomized Kaczmarz (RK) and randomized sparse Kaczmarz (RSK).
We consider the following RSKA variants:

  1. (a)

v1: RSKA with $\mathbf{W}=\mathbf{I}$, i.e. $\alpha=1$, with a coupling such
that
$\mathbf{P}\mathbf{W}\mathbf{D}^{-2}=\frac{\alpha\mathbf{I}}{\|\mathbf{A}\|^{2}_{F}}$
i.e $p_{i}=\frac{\|a_{i}\|^{2}_{2}}{\|\mathbf{A}\|^{2}_{F}}$.

  2. (b)

v2: RSKA with a uniform weight matrix $\mathbf{W}=\alpha^{*}\mathbf{I}$ i.e
$p_{i}=\frac{\|a_{i}\|^{2}_{2}}{\|\mathbf{A}\|^{2}_{F}}$, where $\alpha^{*}$
is from Corollary 4.2.

  3. (c)

v3: RSKA with a general diagonal weight matrix $\mathbf{W}$ with diagonal
entries $w_{i}$ sample i.i.d. from the uniform distribution on $[0,1]$ and
$p_{i}=\frac{\|a_{i}\|^{2}_{2}}{\|\mathbf{A}\|^{2}_{F}}$.

  4. (d)

v4: RSKA with a general diagonal weight matrix $\mathbf{W}$ with diagonal
entries $w_{i}$ sample i.i.d. from the uniform distribution on $[0,1]$ and set
$p_{i}$ proportional to $\tfrac{\|a_{i}\|^{2}}{w_{i}}$, i.e. we have
$\mathbf{P}\mathbf{W}\mathbf{D}^{-2}=\alpha\mathbf{I}$ with
$\alpha=\big{(}\sum_{j}\tfrac{\|a_{j}\|^{2}}{w_{j}}\big{)}^{-1}$.

The rationale behind these choices are: Version v1 just chooses no weight and
standard probabilities. In v2 we still choose standard probabilities but use a
coupling of weights and probabilities with the uniform weight $\alpha^{*}$ of
which we have seen in Corollary 4.2 that is has a certain optimality property.
In v3 we still use standard probabilites but do not use a coupling of
probabilties and weights. Since we do not have any indications on how to
choose weights in any optimal way, we just used random weights here. In
contrast to v3, we do use a coupling of weights and probabilities in v4, but
again, since we do not have results on how to choose optimal weights, we
choose random ones.

Note that RSK and RK are both special cases of RSKA, namely

  1. (a)

RSKA with $\mathbf{W}=\mathbf{I}$, i.e. $\alpha=1$, with a coupling such that
$\mathbf{P}\mathbf{W}\mathbf{D}^{-2}=\frac{\alpha\mathbf{I}}{\|\mathbf{A}\|^{2}_{F}}$
i.e $p_{i}=\frac{\|a_{i}\|^{2}_{2}}{\|\mathbf{A}\|^{2}_{F}}$ and $\eta=1$,

  2. (b)

RSKA with $\mathbf{W}=\mathbf{I}$, i.e. $\alpha=1$, with a coupling such that
$\mathbf{P}\mathbf{W}\mathbf{D}^{-2}=\frac{\alpha\mathbf{I}}{\|\mathbf{A}\|^{2}_{F}}$
i.e $p_{i}=\frac{\|a_{i}\|^{2}_{2}}{\|\mathbf{A}\|^{2}_{F}}$, $\eta=1$ and
$\lambda=0$.

Synthetic data for the experiments is generated as follows: All elements of
the data matrix $\mathbf{A}\in\mathbb{R}^{m\times n}$ are chosen independent
and identically distributed from the standard normal distribution
$\mathcal{N}(0,1)$. We constructed overdetermined, square, and underdetermined
linear systems. To construct sparse solutions $\hat{x}\in\mathbb{R}^{n}$ we
choose $s$ indices from $\left\\{1,\dots,n\right\\}$ at random and placed
zeros at these positions. and the corresponding right hand sides are
$b=\mathbf{A}\hat{x}\in\mathbb{R}^{m}$ while the respective noisy right hand
sides are $b^{\delta}$ and are obtained by adding Gaussian noise (see Section
5.4 below). We generally chose $\eta=1+\frac{1}{10}\text{min}(m,n)$ for RSKA,
unless something else is indicated. Note that in the overdetermined case with
no noise there will be a unique solution $\hat{x}$ since with probability $1$
the matrices $\mathbf{A}$ have full rank, and so all methods are expected to
converge to the same solution $\hat{x}$ in this case.

For each experiment, we run independent trials each starting with the initial
iterate $x_{0}=0$. We measure performance by plotting the relative residual
error $\|\mathbf{A}x-b\|/\|b\|$ and the error $\|x_{k}-\hat{x}\|/\|\hat{x}\|$
against the number of full iterations. Thick line shows mean over the total
number of trials and shaded area which represent the standard deviation over
the trials are plotted when appropriate.

Figure 1 shows the result for a five times overdetermined and consistent
system without noise where the value $\lambda=1$ was used for RSK and RSKA.
Note that the usual RK and RSKA variants performs consistently well over all
trials, while the performance of RSK differs drastically between different
instances. Moreover we observe experimentally that choosing the theoretically
optimal overrelaxation parameter $\alpha^{*}$ from Corollary 4.2 for the RSKA
v2 method give us faster convergence.

Figure 2 and Figure 3 shows the results respectively for a two times resp.
five times underdetermined and consistent system without noise where the
values $\lambda=1$, resp. $\lambda=3$ was used for RSK and RSKA. Methods like
RSK and RSKA take advantage of the fact that the vectors $\hat{x}$ are very
sparse. Moreover, in Figure 2, the RSK and RK methods do not reduce the
residual as fast as the RSKA method. However, since the problem is
underdetermined, the RK method does not converge to a sparse solution and
hence, the error does not converge to zero. Figure 4 and Figure 5 show results
for further values of $m$ and $n$ and show similar behavior as Figure 2.

(a) Relative Residual

(b) Error

Figure 1: A comparison of randomized Kaczmarz (blue), randomized sparse
Kaczmarz (green) and RSKA method (black), $m=100,n=20,$ sparsity $s=10$,
$\eta=11$, $\lambda=1$, no noise and 20 runs. Thick line shows mean over all
trials.

(a) Relative Residual

(b) Error

Figure 2: A comparison of randomized Kaczmarz (blue), randomized sparse
Kaczmarz (green) and RSKA method (black), $m=100,n=200,$ sparsity $s=10$,
$\eta=11$, $\lambda=1$, no noise and 10 runs. Thick line shows mean over all
trials, shaded area represent the standard deviation.

(a) Relative Residual

(b) Error

Figure 3: A comparison of randomized Kaczmarz (blue), randomized sparse
Kaczmarz (green) and RSKA method (black), $m=100,n=500,$ sparsity $s=10$,
$\eta=11$, $\lambda=3$, no noise and 10 runs. Thick line shows mean over all
trials, shaded area represent the standard deviation.

(a) Relative Residual

(b) Error

Figure 4: A comparison of randomized Kaczmarz (blue), randomized sparse
Kaczmarz (green) and RSKA method (black), $m=200,n=600,$ sparsity $s=10$,
$\eta=21$, $\lambda=3$, no noise and 5 runs. Thick line shows mean over all
trials, shaded area represent the standard deviation.

(a) Relative Residual

(b) Error

Figure 5: A comparison of randomized Kaczmarz (blue), randomized sparse
Kaczmarz (green) and RSKA method (black), $m=300,n=300,$ sparsity $s=10$,
$\eta=31$, $\lambda=3$, no noise and 5 runs. Thick line shows mean over all
trials, shaded area represent the standard deviation.

###  5.1 The effect of the number of threads $\eta$

In Figure 6, 7, 8, and 9, we see the effects of the number of threads $\eta$
in the error of Algorithm 1 for the variant $v2$. We used under- and
overdetermined and consistent system, with $\lambda\in\\{0.01,3\\}$ . For the
small $\lambda=0.01$, the RSKA behaves almost like the standard randomized
Kaczmarz with averaging from moorman2021randomized , while for the larger
value $\lambda=3$, we see the typical behavior for the sparse Kaczmarz method
which stagnates from time to time and switches to faster improvement inbetween
LSW14 . As the number of threads $\eta$ increases, we see a corresponding
decrease in the number of iterations needed to reach a certain accuracy,
however at some point increasing $\eta$ does not improve the method in
accordance with Remark 1. For smaller values of $\eta$ we roughly see an
$\eta$-times speedup in the number of iterations. Thus it is clear that the
averaging will pay off as soon as the updates in RSKA can be done in parallel.

(a) Relative Residual

(b) Error

Figure 6: The effect of the number of threads $\eta$ on the error and relative
residual versus iteration for Algorithm 1 on the variant $v2$. $m=200,n=600,$
sparsity $s=10$, $\lambda=0.01$, no noise.

(a) Relative Residual

(b) Error

Figure 7: The effect of the number of threads $\eta$ on the error and relative
residual versus iteration for Algorithm 1 on the variant $v2$. $m=200,n=600,$
sparsity $s=10$, $\lambda=3$, no noise.

(a) Relative Residual

(b) Error

Figure 8: The effect of the number of threads $\eta$ on the error and relative
residual versus iteration for Algorithm 1 on the variant $v2$. $m=100,n=10,$
sparsity $s=10$, $\lambda=0.01$, no noise.

(a) Relative Residual

(b) Error

Figure 9: The effect of the number of threads $\eta$ on the error and relative
residual versus iteration for Algorithm 1 on the variant $v2$. $m=100,n=10,$
sparsity $s=10$, $\lambda=3$, no noise.

###  5.2 The effect of the relaxation parameter $\alpha$

In Figure 10, we observe the effect on the convergence rate as we vary the
relaxation parameter $\alpha$. We used an underdetermined and consistent
system with $\lambda=0.1$, $\eta=8$ with variant $v2$. In fact increasing
$\alpha$ allow us to get smaller error, however the method can ultimately
diverge for some larger values of the relaxation parameter $\alpha$. This is
observed in Figure 11 which plot the relative residual and the error after 100
iterations on a smaller example for various relaxations parameters $\alpha$
and batch sizes $\eta$. The theoretically optimal parameter $\alpha^{*}$ from
Corollary 4.2 is indicated as a dot. We used an overdetermined and consistent
system with $\lambda=1$ with variant $v2$. The plots confirms that $\alpha$
can not be chosen too large, i.e there exists an $\eta$-depended upper bound
for $\alpha$ which leads to convergence (cf. Theorem 4.1). However we also
observe that a larger relaxation parameter than $\alpha^{*}$ from Corollary
4.2 leads to even faster convergence.

(a) Relative Residual

(b) Error

Figure 10: The effect of the relaxation parameter $\alpha$ on the error and
relative residual versus iteration for Algorithm 1 on the variant $v2$.
$m=200,n=600,$ sparsity $s=10$, $\eta=8$, $\lambda=0.1$, no noise. Figure 11:
The effect of the relaxation parameter $\alpha$ for various values of $\eta$
on the relative residual after $100$ iterations of Algorithm 1 with the
variant $v2$. $m=100,n=10,$ sparsity $s=10$, $\lambda=6$, no noise. Circle
markers are estimates of the optimal relaxation parameter using Corollary 4.2.

###  5.3 The effect of the sparsity parameter $\lambda$

In Figure 12, we see the effects of the sparsity parameter $\lambda$ on the
approximation error of Algorithm 1 and the randomized Kaczmarz method (RK). We
used an underdetermined and consistent system with $\eta=21$ with variant $v2$
of RSKA. We observed that as you increase the sparsity parameter $\lambda$,
the relative residual and the reconstruction error get worse and RSK is more
affected by this behavior whereas RSKA keep it performance along different
$\lambda$. The first row of Figure 12 correspond to experiment where we run
1000 iterations and in the second row the methods are run for
$(1+\lambda)1000$ iterations for all $\lambda$.

(a) Relative Residual

(b) Error

(c) Relative Residual

(d) Error

Figure 12: A comparison of randomized sparse Kaczmarz (blue) and RSKA-$v2$
method (orange) in terms of the sparsity parameter $\lambda$. $m=200,n=600,$
sparsity $s=10$, $\eta=21$, no noise. In the first row $(a)$ and $(b)$,
methods are run for 1000 iterations for all $\lambda$ whereas in the second
row $(c)$ and $(d)$, methods are run for $(1+\lambda)1000$ iterations for all
$\lambda$.

###  5.4 Noisy case

In this part we are interested in the effectiveness of the RSKA method on
inconsistent systems. We construct a sparse $\hat{x}$ with normally
distributed non-zero entries and set $b=\mathbf{A}\hat{x},\quad
b^{\epsilon}=b+\epsilon$ where $\epsilon$ is a random vector uniformly
distributed on a sphere with radius $\ell$ that correspond to the relative
noise level such that $\|b-b^{\epsilon}\|_{2}\leq\ell$. Figures 13, 14 and 15
show the results for noisy right hand sides, all with $\lambda=1$ for RSK and
RSKA. Figure 14 uses a five times overdetermined system with $10\%$ relative
noise, Fig. 15 has the same noise level and a five times underdetermined
system. In the underdetermined case, all methods consistently stagnate at a
residual level which is comparable to the noise level, however, in all
settings, RSKA variants achieves faster convergence than RSK which in turn is
faster than RK. Regarding the reconstruction error, RSKA and RSK achieve
reconstructions with an error in the size of the noise level, while RSKA
achieves an even lower reconstruction error.

(a) Relative Residual

(b) Error

Figure 13: A comparison of randomized Kaczmarz (blue), randomized sparse
Kaczmarz (green) and RSKA method (black), $m=100,n=100,$ sparsity $s=10$,
$\eta=11$, $\lambda=1$, noise level $l=0.1$ and 5 runs. Thick line shows mean
over all trials, shaded area represent the standard deviation.

(a) Relative Residual

(b) Error

Figure 14: A comparison of randomized Kaczmarz (blue), randomized sparse
Kaczmarz (green) and RSKA method (black), $m=500,n=100,$ sparsity $s=10$,
$\eta=11$, $\lambda=1$, noise level $l=0.1$ and 5 runs. Thick line shows mean
over all trials, shaded area represent the standard deviation.

(a) Relative Residual

(b) Error

Figure 15: A comparison of randomized Kaczmarz (blue), randomized sparse
Kaczmarz (green) and RSKA method (black), $m=100,n=500,$ sparsity $s=10$,
$\eta=11$, $\lambda=1$, noise level $l=0.1$ and 5 runs. Thick line shows mean
over all trials, shaded area represent the standard deviation.

##  6 Conclusion

We proved that the iterates of the randomized sparse Kaczmarz with averaging
method (Algorithm 1) are expected to converge linearly for consistent linear
systems. Moreover we show that the iterates reach an error threshold in the
order of the noise-level in the noisy case. We gave a general error bound in
terms of the sparsity parameter $\lambda$, the number of threads $\eta$ and a
relaxation parameter $\alpha$. Numerical experiments show that the method
performs consistently well over a range of values of $\lambda$ (which is
different for the version without averaging), very good reconstruction quality
as $\lambda$ increases, confirm the theoretical results, and demonstrate the
benefit of using Algorithm 1 to recover sparse solutions of linear systems,
even in the noisy case. We demonstrate that the rate of convergence for
Algorithm 1 improve both in theory and practice as the number of threads
$\eta$ increases. Moreover, we derive an optimal value for the relaxation
parameter $\alpha$ which gives the fastest convergence speed and our numerical
experiments indicate that this optimal value for $\alpha$ (in v2 and v4 of the
algorithms in Section 5) does indeed provide fast convergence.

#### Conflict of interest.

The authors declare that they have no conflict of interest.

#### Availability of data and materials

The data for the numerical example are randomly generated matrices.

#### Code availability.

The computational code is only prototypal, but it is available from the
authors upon request.

## References

  * (1) Beck, A., Teboulle, M.: Mirror descent and nonlinear projected subgradient methods for convex optimization.  Operations Research Letters 31(3), 167–175 (2003)
  * (2) Cai, J.F., Osher, S., Shen, Z.: Convergence of the linearized Bregman iteration for $\ell_{1}$-norm minimization.  Mathematics of Computation 78(268), 2127–2136 (2009)
  * (3) Cai, J.F., Osher, S., Shen, Z.: Linearized Bregman iterations for compressed sensing.  Mathematics of computation 78(267), 1515–1536 (2009)
  * (4) Candès, E.J.: Compressive sampling.  In: International Congress of Mathematicians. Vol. III, pp. 1433–1452. Eur. Math. Soc., Zürich (2006)
  * (5) Candès, E.J., Romberg, J., Tao, T.: Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information.  IEEE Transactions on information theory 52(2), 489–509 (2006)
  * (6) Chen, S.S., Donoho, D.L., Saunders, M.A.: Atomic decomposition by basis pursuit.  SIAM J. Sci. Comput. 20(1), 33–61 (1998).  DOI 10.1137/S1064827596304010.  URL https://doi.org/10.1137/S1064827596304010
  * (7) Donoho, D.L.: Compressed sensing.  IEEE Trans. Inform. Theory 52(4), 1289–1306 (2006).  DOI 10.1109/TIT.2006.871582.  URL https://doi.org/10.1109/TIT.2006.871582
  * (8) Donoho, D.L., Tanner, J.: Sparse nonnegative solution of underdetermined linear equations by linear programming.  Proceedings of the national academy of sciences 102(27), 9446–9451 (2005)
  * (9) D’Orazio, R., Loizou, N., Laradji, I., Mitliagkas, I.: Stochastic mirror descent: Convergence analysis and adaptive variants via the mirror stochastic polyak stepsize.  arXiv preprint arXiv:2110.15412 (2021)
  * (10) Du, K., Si, W.T., Sun, X.H.: Randomized extended average block kaczmarz for solving least squares.  SIAM Journal on Scientific Computing 42(6), A3541–A3559 (2020)
  * (11) Friedlander, M.P., Tseng, P.: Exact regularization of convex programs.  SIAM Journal on Optimization 18(4), 1326–1350 (2008)
  * (12) Garey, M.R., Johnson, D.S.: Computers and intractability.  A Series of Books in the Mathematical Sciences. W. H. Freeman and Co., San Francisco, Calif. (1979).  A guide to the theory of NP-completeness
  * (13) Gower, R.M., Molitor, D., Moorman, J., Needell, D.: On adaptive sketch-and-project for solving linear systems.  SIAM Journal on Matrix Analysis and Applications 42(2), 954–989 (2021).  DOI 10.1137/19M1285846.  URL https://doi.org/10.1137/19M1285846
  * (14) Hanke, M., Niethammer, W.: On the acceleration of Kaczmarz’s method for inconsistent linear systems.  Linear Algebra and its Applications 130, 83–98 (1990)
  * (15) Herman, G.T., Lent, A., Lutz, P.H.: Relaxation methods for image reconstruction.  Commun. ACM 21(2), 152–158 (1978).  DOI 10.1145/359340.359351.  URL https://doi.org/10.1145/359340.359351
  * (16) Horn, R.A., Horn, R.A., Johnson, C.R.: Topics in matrix analysis.  Cambridge university press (1994)
  * (17) Hounsfield, G.N.: Computerized transverse axial scanning (tomography): Part 1. description of system.  The British journal of radiology 46(552), 1016–1022 (1973)
  * (18) Jiao, Y., Jin, B., Lu, X.: Preasymptotic convergence of randomized Kaczmarz method.  Inverse Problems 33(12), 125,012, 21 (2017).  DOI 10.1088/1361-6420/aa8e82.  URL https://doi.org/10.1088/1361-6420/aa8e82
  * (19) Kaczmarz, S.: Angenäherte Auflösung von Systemen linearer Gleichungen.  Bull. Internat. Acad. Polon. Sci. Lettres A pp. 355–357 (1937)
  * (20) Khan, U.A., Moura, J.M.: Distributed Kalman filters in sensor networks: Bipartite fusion graphs.  In: 2007 IEEE/SP 14th Workshop on Statistical Signal Processing, pp. 700–704. IEEE (2007)
  * (21) Lan, G., Nemirovski, A., Shapiro, A.: Validation analysis of mirror descent stochastic approximation method.  Mathematical programming 134(2), 425–458 (2012)
  * (22) Loizou, N., Vaswani, S., Laradji, I.H., Lacoste-Julien, S.: Stochastic polyak step-size for sgd: An adaptive learning rate for fast convergence.  In: International Conference on Artificial Intelligence and Statistics, pp. 1306–1314. PMLR (2021)
  * (23) Lorenz, D.A., Schöpfer, F., Wenger, S.: The linearized Bregman method via split feasibility problems: Analysis and generalizations.  SIAM J. Imaging Sciences 7(2), 1237–1262 (2014)
  * (24) Lorenz, D.A., Wenger, S., Schöpfer, F., Magnor, M.: A sparse Kaczmarz solver and a linearized Bregman method for online compressed sensing.  In: 2014 IEEE international conference on image processing (ICIP), pp. 1347–1351. IEEE (2014)
  * (25) Mertikopoulos, P., Lecouat, B., Zenati, H., Foo, C.S., Chandrasekhar, V., Piliouras, G.: Optimistic mirror descent in saddle-point problems: Going the extra(-gradient) mile.  In: International Conference on Learning Representations (2019).  URL https://openreview.net/forum?id=Bkg8jjC9KQ
  * (26) Mertikopoulos, P., Staudigl, M.: Stochastic mirror descent dynamics and their convergence in monotone variational inequalities.  Journal of optimization theory and applications 179(3), 838–867 (2018)
  * (27) Miao, C.Q., Wu, W.T.: On greedy randomized average block Kaczmarz method for solving large linear systems.  Journal of Computational and Applied Mathematics 413, 114,372 (2022)
  * (28) Moorman, J.D., Tu, T.K., Molitor, D., Needell, D.: Randomized Kaczmarz with averaging.  BIT Numerical Mathematics 61(1), 337–359 (2021)
  * (29) Necoara, I.: Faster randomized block Kaczmarz algorithms.  SIAM Journal on Matrix Analysis and Applications 40(4), 1425–1452 (2019)
  * (30) Needell, D.: Randomized Kaczmarz solver for noisy linear systems.  BIT Numerical Mathematics 50(2), 395–403 (2010)
  * (31) Needell, D., Tropp, J.A.: Paved with good intentions: analysis of a randomized block Kaczmarz method.  Linear Algebra and its Applications 441, 199–221 (2014)
  * (32) Nemirovski, A., Juditsky, A., Lan, G., Shapiro, A.: Robust stochastic approximation approach to stochastic programming.  SIAM Journal on optimization 19(4), 1574–1609 (2009)
  * (33) Nemirovski, A.S., Juditsky, A.B., Lan, G., Shapiro, A.: Robust Stochastic Approximation Approach to Stochastic Programming.  SIAM Journal on Optimization 19(4), 1574–1609 (2009).  DOI 10.1137/070704277.  URL https://hal.archives-ouvertes.fr/hal-00976649
  * (34) Nemirovskij, A.S., Yudin, D.B.: Problem complexity and method efficiency in optimization.  Wiley-Interscience (1983)
  * (35) Nesterov, Y.: Efficiency of coordinate descent methods on huge-scale optimization problems.  SIAM Journal on Optimization 22(2), 341–362 (2012)
  * (36) Olshanskii, M.A., Tyrtyshnikov, E.E.: Iterative methods for linear systems: theory and applications.  SIAM (2014)
  * (37) Patrascu, A., Necoara, I.: Nonasymptotic convergence of stochastic proximal point methods for constrained convex optimization.  The Journal of Machine Learning Research 18(1), 7204–7245 (2017)
  * (38) Petra, S.: Randomized sparse block Kaczmarz as randomized dual block-coordinate descent.  Analele Stiintifice Ale Universitatii Ovidius Constanta-Seria Matematica 23(3), 129–149 (2015)
  * (39) Popa, C.: Convergence rates for Kaczmarz-type algorithms.  Numer. Algorithms 79(1), 1–17 (2018).  DOI 10.1007/s11075-017-0425-7.  URL https://doi.org/10.1007/s11075-017-0425-7
  * (40) Rabelo, J.C., Saporito, Y.F., Leitão, A.: On stochastic Kaczmarz type methods for solving large scale systems of ill-posed equations.  Inverse Problems 38(2), Paper No. 025,003, 23 (2022).  DOI 10.1088/1361-6420/ac3f80.  URL https://doi.org/10.1088/1361-6420/ac3f80
  * (41) Richtárik, P., Takáč, M.: Parallel coordinate descent methods for big data optimization.  Mathematical Programming 156(1), 433–484 (2016)
  * (42) Richtárik, P., Takác, M.: Stochastic reformulations of linear systems: algorithms and convergence theory.  SIAM Journal on Matrix Analysis and Applications 41(2), 487–524 (2020)
  * (43) Schöpfer, F., Lorenz, D.A.: Linear convergence of the randomized sparse Kaczmarz method.  Mathematical Programming 173(1), 509–536 (2019).  URL https://link.springer.com/article/10.1007/s10107-017-1229-1
  * (44) Schöpfer, F., Lorenz, D.A., Tondji, L., Winkler, M.: Extended randomized kaczmarz method for sparse least squares and impulsive noise problems.  Lineare Algebra and Applications 652, 132–154 (2022)
  * (45) Stiefel, E.: Methods of conjugate gradients for solving linear systems.  J. Res. Nat. Bur. Standards 49, 409–435 (1952)
  * (46) Strohmer, T., Vershynin, R.: A randomized Kaczmarz algorithm with exponential convergence.  Journal of Fourier Analysis and Applications 15(2), 262–278 (2009)
  * (47) Tropp, J.A.: Improved analysis of the subsampled randomized Hadamard transform.  Advances in Adaptive Data Analysis 3(01n02), 115–126 (2011)
  * (48) Yin, W.: Analysis and generalizations of the linearized Bregman method.  SIAM Journal on Imaging Sciences 3(4), 856–877 (2010)
  * (49) Zouzias, A., Freris, N.M.: Randomized extended Kaczmarz for solving least squares.  SIAM Journal on Matrix Analysis and Applications 34(2), 773–793 (2013)