#  Faster quantum-inspired algorithms for solving linear systems

Changpeng Shao<EMAIL_ADDRESS>School of Mathematics, University
of Bristol, UK  Ashley Montanaro<EMAIL_ADDRESS>School of
Mathematics, University of Bristol, UK  Phasecraft Ltd. UK

###### Abstract

We establish an improved classical algorithm for solving linear systems in a
model analogous to the QRAM that is used by quantum linear solvers. Precisely,
for the linear system $A{\bf x}={\bf b}$, we show that there is a classical
algorithm that outputs a data structure for ${\bf x}$ allowing sampling and
querying to the entries, where ${\bf x}$ is such that $\|{\bf x}-A^{-1}{\bf
b}\|\leq\epsilon\|A^{-1}{\bf b}\|$. This output can be viewed as a classical
analogue to the output of quantum linear solvers. The complexity of our
algorithm is $\widetilde{O}(\kappa_{F}^{6}\kappa^{2}/\epsilon^{2})$, where
$\kappa_{F}=\|A\|_{F}\|A^{-1}\|$ and $\kappa=\|A\|\|A^{-1}\|$. This improves
the previous best algorithm [Gilyén, Song and Tang, arXiv:2009.07268] of
complexity $\widetilde{O}(\kappa_{F}^{6}\kappa^{6}/\epsilon^{4})$. Our
algorithm is based on the randomized Kaczmarz method, which is a particular
case of stochastic gradient descent. We also find that when $A$ is row sparse,
this method already returns an approximate solution ${\bf x}$ in time
$\widetilde{O}(\kappa_{F}^{2})$, while the best quantum algorithm known
returns $|{\bf x}\rangle$ in time $\widetilde{O}(\kappa_{F})$ when $A$ is
stored in the QRAM data structure. As a result, assuming access to QRAM and if
$A$ is row sparse, the speedup based on current quantum algorithms is
quadratic.

##  1 Introduction

Since the discovery of the Harrow-Hassidim-Lloyd algorithm [16] for solving
linear systems, many quantum algorithms for solving machine learning problems
were proposed, e.g. [20, 19, 21, 24, 31, 2]. Most of them claimed that quantum
computers could achieve exponential speedups over classical algorithms.
However, as shown by Tang [36] and some subsequent works [6, 10, 17, 35, 7,
11, 5] (currently known as quantum-inspired algorithms), for many previously
discovered quantum machine learning algorithms, the quantum computer actually
only achieves a polynomial speedup in the dimension. Tang’s work makes us
reconsider the advantages of quantum computers in machine learning. For
instance, one may wonder what is the best possible separation between quantum
and classical computers for solving a given machine learning problem?

Here we focus on the solving of linear systems, which is a fundamental problem
in machine learning. For a linear system $A{\bf x}={\bf b}$, the currently
best quantum algorithm seems to be the one building on the technique of block-
encoding. For simplicity, assume that the singular values of $A$ lie in
$[1/\kappa,1]$, where $\kappa$ is the condition number of $A$. Let
$\alpha\in\mathbb{R}^{+}$, if $A/\alpha$ is efficiently encoded as a block of
a unitary, the complexity to obtain $|A^{-1}{\bf b}\rangle$ up to error
$\epsilon$ is $\widetilde{O}(\alpha\kappa)$ [4], where the $\widetilde{O}$
notation hides logarithmic factors in all parameters (including
$\kappa,\epsilon$ and the size of $A$). For example, when we have the QRAM
data structure for $A$, we can take $\alpha=\|A\|_{F}$.111In fact,
$\alpha=\mu(A)$ for some $\mu(A)\leq\|A\|_{F}$ (see [12]). The definition of
$\mu(A)$ is quite complicated. In this paper, for convenience of the
comparison we just use $\|A\|_{F}$.  When $A$ is $s$-sparse and given in the
sparse-access input model, we can take $\alpha=s\|A\|_{\max}$. In the
classical case, the quantum-inspired algorithm given in [6] costs
$\widetilde{O}(\kappa_{F}^{6}\kappa^{22}/\epsilon^{6})$ to output a classical
analogue to $|A^{-1}{\bf b}\rangle$, where
$\kappa_{F}=\|A\|_{F}\|A^{-1}\|$,222In the literature, $\kappa_{F}$ does not
have a standard name. It is named the scaled condition number in [34, 15]. and
$\epsilon$ is the accuracy. This result was recently improved to
$\widetilde{O}(\kappa_{F}^{6}\kappa^{6}/\epsilon^{4})$ [11]. Although the
quantum-inspired algorithms are polylog in the dimension, the heavy dependence
on $\kappa_{F},\kappa,\epsilon$ may make them hard to use in practice [1]. So
it is necessary to improve the efficiency of quantum-inspired algorithms.
Meanwhile, this can help us better understand the quantum advantages in
certain problems.

QRAM is a widely used data structure in many quantum machine learning
algorithms (e.g. [20, 19, 21]). It allows us efficiently to prepare a quantum
state corresponding to the data and to perform quantum linear algebraic
operations. The quantum-inspired algorithms usually assume a data structure
which is a classical analogue to the QRAM. For a vector, this data structure
allows us to do the sampling and query operations to the entries. The sampling
operation can be viewed as measurements of the quantum state.

As for the problem of solving linear systems, in the QRAM model, the input of
quantum algorithms is the QRAM data structures for $A$ and ${\bf b}$, and the
output is $|{\bf x}\rangle$ such that $\||{\bf x}\rangle-|A^{-1}{\bf
b}\rangle\|\leq\epsilon$. In the classical setting, the input of quantum-
inspired algorithms is the classical analogue of QRAM data structures for $A$
and ${\bf b}$, and the output is a data structure for ${\bf x}$ such that
$\|{\bf x}-A^{-1}{\bf b}\|\leq\epsilon\|A^{-1}{\bf b}\|$. This makes the
comparison of quantum/quantum-inspired algorithms reasonable. There are also
some other quantum models, such as the sparse-access input model. This allows
efficient quantum algorithms for solving sparse linear systems [8]. However,
because of the similarity between QRAM and the model used in quantum-inspired
algorithms, in this paper we will mainly focus on the QRAM model for quantum
linear solvers.

Randomized iterative methods have received a lot of attention recently in
solving linear systems because of their connection to the stochastic gradient
descent method. The randomized Kaczmarz method [34] is a typical example. It
was first proposed by Kaczmarz in 1937 [18], and rediscovered by Gordon,
Bender and Herman in 1970 [14]. After that, there appear many generalizations
[15, 26, 27, 32]. Kaczmarz method is a particular case of the stochastic
gradient descent algorithm for solving linear systems. In its simplest
version, at each step of iteration, it randomly chooses a row of the matrix
with probability corresponding to the norm of that row, then finds the closest
solution that satisfies this linear constraint. Geometrically, the next point
is the orthogonal projection of the previous one onto the hyperplane defined
by the chosen linear constraint. The simple iterative structure and clear
geometric meaning greatly simplifies the analysis of this method and makes it
popular in applications like computed tomography [3]. Usually, the randomized
iterative methods need to do some importance sampling according to the norms
of the rows. Sampling access is made available by the data structure for
quantum-inspired algorithms and the QRAM for quantum linear solvers. So in
this paper, for the fairness of comparison we also assume a similar data
structure for the randomized iterative methods.

###  1.1 Main results

In this paper, we propose two quantum-inspired algorithms for solving linear
systems. One is based on the randomized Kaczmarz method [34], and the other
one is based on the randomized coordinate descent method [22]. The second
algorithm only works for symmetric positive-definite (SPD) matrices. Our
results are summarized in Table 1.

| Complexity | Reference | Assumptions
---|---|---|---
Quantum algorithm | $\widetilde{O}(\kappa_{F})$ |  [4], 2018 |
Randomized classical algorithm | $\widetilde{O}(s\kappa_{F}^{2})$ |  [34], 2009 | row sparse
$\widetilde{O}(s{{\rm Tr}(A)\|A^{-1}\|})$ |  [22], 2010 | sparse, SPD
Quantum-inspired algorithm | $\widetilde{O}(\kappa_{F}^{6}\kappa^{6}/\epsilon^{4})$ |  [11], 2020 |
$\widetilde{O}(\kappa_{F}^{6}\kappa^{2}/\epsilon^{2})$ | Theorem 12 |
$\widetilde{O}({{\rm Tr}(A)^{3}\|A^{-1}\|^{3}\kappa}/{\epsilon^{2}})$ | Theorem 13 | SPD
Table 1: Comparison of different algorithms for solving the linear system
$A{\bf x}={\bf b}$, where
$\kappa_{F}=\|A\|_{F}\|A^{-1}\|,\kappa=\|A\|\|A^{-1}\|$ and $s$ is the row
sparsity. SPD = symmetric positive definite.

In the table, the quantum algorithm is based on the QRAM model. The randomized
classical and quantum-inspired algorithms assume a data structure that is a
classical analogue to the QRAM (see Section 3). So all those algorithms are
based on similar data structures, which make the comparison in the table fair
enough.

  * •

When $A$ is dense, the complexity of our algorithm is
$\widetilde{O}(\kappa_{F}^{6}\kappa^{2}/\epsilon^{2})$. If $A$ is additionally
SPD, the complexity becomes $\widetilde{O}({{\rm
Tr}(A)^{3}\|A^{-1}\|^{3}\kappa}/{\epsilon^{2}})$. This reduces the complexity
of the previous quantum-inspired algorithm given by Gilyén et al. [11] by at
least a 4th power in terms of the dependence on $\kappa$, and improves the
dependence on $1/\epsilon$ from quartic to quadratic. One result not shown in
the table is about the description of the solution. In [11], the authors
define an $s$-sparse description of a vector ${\bf x}$ as an $s$-sparse vector
${\bf y}$ such that ${\bf x}=A^{\dagger}{\bf y}$. The sparsity of ${\bf y}$
determines the runtime of querying one entry of ${\bf x}$. In [11], they
obtain an $\widetilde{O}(\kappa_{F}^{2}\kappa^{2}/\epsilon^{2})$-sparse
description of the approximate solution in cost
$\widetilde{O}(\kappa_{F}^{6}\kappa^{2}/\epsilon^{4})$. In comparison, we
obtain an $\widetilde{O}(\kappa_{F}^{2})$-sparse description in cost
$\widetilde{O}(\kappa_{F}^{6}\kappa^{2}/\epsilon^{2})$ by Theorem 12.

  * •

When $A$ is row sparse, i.e., each row is sparse, the randomized Kaczmarz
method can find an approximate sparse solution in time
$\widetilde{O}(s\kappa_{F}^{2})$, where $s$ is the row sparsity. Moreover, the
output is a vector, which is stronger than the outputs of quantum and quantum-
inspired algorithms. Assuming $s=\widetilde{O}(1)$, this result means that the
quantum speedup of solving row sparse linear systems is quadratic if using the
QRAM data structure. It also means that to explore large quantum speedups in
solving sparse linear systems, we may need to move our attentions to other
quantum models. For example, it is known that when $A$ is row and column
sparse, the quantum algorithm has complexity
$\widetilde{O}(s\|A\|_{\max}\|A^{-1}\|)$ in the sparse access input model [4],
which is better than the quantum result we show in the table.

  * •

When $A$ is sparse and SPD, the randomized coordinate descent method can find
an approximate sparse solution (also a sparse vector) in time
$\widetilde{O}(s{\rm Tr}(A)\|A^{-1}\|)$. In this special case, it may happen
that quantum computers do not achieve a speedup in solving linear systems. For
example, suppose $A$ is a density matrix with
$\|A\|=\Theta(1),s=\widetilde{O}(1)$, then the randomized coordinate descent
method costs $\widetilde{O}(\kappa)$, which is also the cost of the quantum
linear solver. As shown in [30], to solve a sparse SPD linear system in a
quantum computer, the dependence on $\kappa$ is at least linear. However, for
some special SPD linear systems, quantum algorithms can achieve quadratic
speedups in terms of $\kappa$.333In [30], the authors mentioned two special
cases. Case 1: Access to a block-encoding of $I-\eta A$ and $\|A^{-1}{\bf
b}\|\in O(\kappa)$, where $\eta\in(0,1]$. Case 2: $A$ is a sum of SPD local
Hamiltonians, ${\bf b}$ is sparse and a parameter $\gamma\in O(1)$ that
quantifies the overlap of ${\bf b}$ with the subspace where $A$ is non-
singular. However, the special cases mentioned in [30] seem not directly
related to the case we mentioned above.

  * •

In Table 1, the complexity of the quantum algorithm is the cost to obtain the
quantum state corresponding to the solution. For quantum-inspired algorithms,
the table shows the complexity to obtain the data structure for the solution
that allows sampling and querying. This output can be viewed as a classical
analogue to the output of the quantum algorithm. If we want to estimate the
norm of the solution, then there is an extra factor $1/\epsilon$ for the
quantum algorithm [4] and an extra factor $1/\epsilon^{2}$ for the quantum-
inspired algorithms [6] (also see Lemma 7). However, for the randomized
classical algorithms, the complexity does not change because their outputs are
sparse vectors whose sparsity has the same order as the time complexity. So to
estimate the norm of the solution, the quantum and quantum-inspired algorithms
are worse than the randomized classical algorithms in terms of the precision
in the row sparse case.

  * •

For solving a row sparse linear system, the table shows that the complexities
of quantum and classical randomized algorithms are logarithmic in the
precision $\epsilon$. However, if we write down the complexity more precisely,
the classical randomized algorithms may have better dependence on $\epsilon$,
at least compared with the currently best quantum algorithm. Indeed, for the
classical randomized algorithms, the dependence on $\epsilon$ is
$O(\log(1/\epsilon))$, while it is $O(\log^{2}(1/\epsilon))$ for the quantum
algorithm [4]. If the linear system is also column sparse, the dependence on
$\epsilon$ can be reduced to $O(\log(1/\epsilon))$ in the quantum case [23].

###  1.2 Summary of the techniques

Our main algorithm is based on a generalization of Kaczmarz method. As briefly
introduced above (also see Section 2 for more details), the Kaczmarz method
aims to solve one linear constraint at each step of the iteration, which can
be viewed as a particular case of the stochastic gradient method. To design an
efficient quantum-inspired algorithm, we concentrate on the dual form of the
Kaczmarz method. Namely, introducing ${\bf y}$ such that ${\bf
x}=A^{\dagger}{\bf y}$, then focusing on the iteration on ${\bf y}$. The
linear system $A{\bf x}={\bf b}$ can be reformulated as the least squares
problem (LSP) ${\bf x}=\arg\min_{{\bf x}}\|A{\bf x}-{\bf b}\|$. It turns out
that ${\bf y}$ converges to the optimal solution of the dual problem of the
LSP. Our idea is to apply the iterative method to find a sparse approximate
solution ${\bf y}$ of the dual problem. The sparsity is guaranteed by the
iterative method. Now ${\bf x}=A^{\dagger}{\bf y}$ gives an approximate
solution of the LSP. Since ${\bf y}$ is sparse, it is easy to query an entry
of ${\bf x}$. To sample an entry of ${\bf x}$, we can use the rejection
sampling idea.

The technical part of this paper is the complexity analysis of the rejection
sampling method. However, it turns out that the main cost of the above
algorithm comes from the computation of ${\bf y}$, whose complexity analysis
is straightforward. It is time-consuming if we directly use the iterative
method because it relates to the calculation of inner products. To alleviate
this, we modify the iteration by using the sampling idea. We prove that this
modified method has the same convergence rate as the original one by choosing
some parameters appropriately.

Notation. For a matrix $A$, we use $A^{\dagger}$ to denote its conjugate
transpose, $A^{-1}$ to denote its Moore-Penrose inverse. The scaled condition
number is $\kappa_{F}=\|A\|_{F}\|A^{-1}\|$, and the condition number is
$\kappa=\|A\|\|A^{-1}\|$, where $\|\cdot\|_{F}$ is the Frobenius norm and
$\|\cdot\|$ is the spectral norm. The row sparsity of $A$ is defined as the
maximum of the number of nonzero zero entries of any row. For a unit vector
${\bf v}$, we sometimes denote it as $|{\bf v}\rangle$, and denote ${\bf
v}^{\dagger}$ as $\langle{\bf v}|$. For two vectors ${\bf a},{\bf b}$, we use
$\langle{\bf a}|{\bf b}\rangle$ or ${\bf a}\cdot{\bf b}$ to denote their inner
product. If $m$ is an integer, we define $[m]:=\\{1,2,\ldots,m\\}$. For matrix
$A$, we denote its $i$-th row as $A_{i*}$, its $j$-th column as $A_{*j}$ and
the $(i,j)$-th entry as $A_{i,j}$. For a vector ${\bf v}$, we denote the
$i$-th component as $v_{i}$. We use $\\{{\bf e}_{1},\ldots,{\bf e}_{m}\\}$ to
enote the standard basis of $\mathbb{C}^{m}$.

##  2 Randomized Kaczmarz algorithm

Let $A$ be an $m\times n$ matrix, ${\bf b}$ be an $m\times 1$ vector. Solving
the linear system $A{\bf x}={\bf b}$ can be described as a quadratic
optimization problem, i.e., the so-called least squares problem (LSP)

$\arg\min_{{\bf x}\in\mathbb{C}^{n}}\quad\|A{\bf x}-{\bf b}\|^{2}.$ (1)

Classically, there are many algorithms for solving the LSP [13]. The best
known classical algorithm might be the conjugate-gradient (CG) method.
Compared to CG, the Kaczmarz algorithm has received much attention recently
due to its simplicity. For instance, at each step of iteration, CG uses the
whole information of the matrix $A$ and $O(mn)$ operations, while Kaczmarz
method only uses one row of $A$ and $O(n)$ operations. So Kaczmarz method
becomes effective when the entire matrix (data set) cannot be loaded into
memory.

######  Definition 1 (Randomized Kaczmarz algorithm [34]).

Let ${\bf x}_{0}$ be an arbitrary initial approximation to the solution of
(1). For $k=0,1,\ldots$, compute

${\bf x}_{k+1}={\bf x}_{k}+\frac{b_{r_{k}}-\langle A_{r_{k}*}|{\bf
x}_{k}\rangle}{\|A_{r_{k}*}\|^{2}}A_{r_{k}*},$ (2)

where $A_{r_{k}*}$ is the $r_{k}$-th row of $A$, and $r_{k}$ is chosen from
$[m]$ randomly with probability proportional to $\|A_{r_{k}*}\|^{2}$.

The Kaczmarz algorithm has a clear geometric meaning, see Figure 1.
Geometrically, ${\bf x}_{k+1}$ is the orthogonal projection of ${\bf x}_{k}$
onto the hyperplane $A_{r_{k}*}\cdot{\bf x}=b_{r_{k}}$. So at each step of
iteration, this method finds the closest vector that satisfies the chosen
linear constraint.

$A_{1*}\cdot{\bf x}=b_{1}$$A_{2*}\cdot{\bf x}=b_{2}$${\bf x}_{0}$${\bf
x}_{1}$${\bf x}_{2}$${\bf x}_{3}$${\bf x}_{4}$${\bf x}_{5}$${\bf x}_{6}$${\bf
x}_{*}$ Figure 1: Illustration of Kaczmarz algorithm when $m=2$, where ${\bf
x}_{*}$ is the optimal solution.

From the definition, we can simplify the notation by defining

$\tilde{{\bf b}}={\rm diag}(\|A_{i*}\|^{-1}:i\in[m]){\bf
b},\quad\widetilde{A}={\rm diag}(\|A_{i*}\|^{-1}:i\in[m])A,$ (3)

so that $\tilde{b}_{i}=b_{i}/\|A_{i*}\|$ and
$\widetilde{A}_{i*}=A_{i*}/\|A_{i*}\|$. Then we can rewrite (2) as

${\bf x}_{k+1}={\bf
x}_{k}+\left(\tilde{b}_{r_{k}}-\langle\widetilde{A}_{r_{k}*}|{\bf
x}_{k}\rangle\right)\widetilde{A}_{r_{k}*}=\left(I-|\widetilde{A}_{r_{k}*}\rangle\langle\widetilde{A}_{r_{k}*}|\right){\bf
x}_{k}+\tilde{b}_{r_{k}}\widetilde{A}_{r_{k}*}.$ (4)

In [33], the Kaczmarz algorithm has been quantized based on the description
(4). As an orthogonal projector, it is not hard to extend
$I-|\widetilde{A}_{i*}\rangle\langle\widetilde{A}_{i*}|$ to a block of a
unitary operator. For example, a unitary of the form
$I\otimes(I-|\widetilde{A}_{i*}\rangle\langle\widetilde{A}_{i*}|)+X\otimes|\widetilde{A}_{i*}\rangle\langle\widetilde{A}_{i*}|$,
where $X$ is Pauli-X. This unitary is closely related to the Grover diffusion
operator when we apply it to $|-\rangle|{\bf x}\rangle$.

######  Lemma 2 (Theorem 2 of [34]).

Let ${\bf x}_{*}$ be the solution of (1). Then the randomized Kaczmarz
algorithm converges to ${\bf x}_{*}$ in expectation, with the average error

${\mathbb{E}}[\|{\bf x}_{T}-{\bf
x}_{*}\|^{2}]\leq(1-\kappa_{F}^{-2})^{T}\|{\bf x}_{0}-{\bf x}_{*}\|^{2},$ (5)

where $\kappa_{F}=\|A\|_{F}\|A^{-1}\|$.

The above result holds for all ${\bf x}_{0}$, so throughout this paper we
shall assume ${\bf x}_{0}=0$ for simplicity. In the above lemma, by Markov’s
inequality, with high probability 0.99, we have $\|{\bf x}_{T}-{\bf
x}_{*}\|^{2}\leq 100(1-\kappa_{F}^{-2})^{T}\|{\bf x}_{*}\|^{2}$. To make sure
the error is bounded by $\epsilon^{2}\|{\bf x}_{*}\|^{2}$ with high
probability, it suffices to choose

$T=O(\kappa_{F}^{2}\log(1/\epsilon)).$ (6)

From (2), we know that there exist $y_{k,0},\ldots,y_{k,k-1}$ such that ${\bf
x}_{k}=\sum_{j=0}^{k-1}y_{k,j}A_{r(j)*},$ that is ${\bf x}_{k}=A^{\dagger}{\bf
y}_{k}$ for some vector ${\bf y}_{k}$. One obvious fact is that ${\bf y}_{k}$
is at most $k$-sparse. In [11], ${\bf y}_{k}$ is called the sparse description
of ${\bf x}_{k}$. So this description comes easily in the Kaczmarz method.
Actually, ${\bf y}_{k}$ converges to the optimal solution of the dual problem
of the LSP (1). More precisely, an alternative formulation of LSP is to find
the least-norm solution of the linear system

$\min_{{\bf x}\in\mathbb{C}^{n}}\quad\frac{1}{2}\|{\bf x}\|^{2},\quad{\rm
s.t.}\quad A{\bf x}={\bf b}.$

The dual problem takes the form:

$\min_{{\bf y}\in\mathbb{C}^{m}}\quad g({\bf y}):=\frac{1}{2}\|A^{\dagger}{\bf
y}\|^{2}-{\bf b}\cdot{\bf y}.$ (7)

The randomized Kaczmarz algorithm (2) is equivalent to one step of the
stochastic gradient descent method [28] applied to
$\frac{1}{2}(A_{r_{k}*}\cdot{\bf x}-b_{r_{k}})^{2}$ with stepsize
$1/\|A_{r_{k}*}\|^{2}$. It is also equivalent to one step of the randomized
coordinate descent method [29] applied to the dual problem (7). Namely, when
we apply this method to the $r_{k}$-th component, the gradient is
$\nabla_{r_{k}}g=\langle A_{r_{k}*}|A^{\dagger}|{\bf y}_{k}\rangle-b_{r_{k}}$.
Setting the stepsize as $1/\|A_{r_{k}*}\|^{2}$ leads to the updating

${\bf y}_{k+1}={\bf y}_{k}+\frac{b_{r_{k}}-\langle A_{r_{k}*}|A^{\dagger}|{\bf
y}_{k}\rangle}{\|A_{r_{k}*}\|^{2}}{\bf e}_{r_{k}},$ (8)

where $\\{{\bf e}_{1},\ldots,{\bf e}_{m}\\}$ is the standard basis of
$\mathbb{C}^{m}$. We can recover the original Kaczmarz iteration (2) by
multiplying $A^{\dagger}$ on both sides of (8).

The Kaczmarz method can be generalized to use multiple rows at each step of
iteration [15]. Instead of orthogonally projecting to a hyperplane determined
by one linear constraint, we can orthogonally project to a vector space
defined by several linear constraints. This idea is similar to the minibatch
stochastic gradient descent. However, this generalization needs to compute the
pseudoinverse of certain matrices. A simple pseudoinverse-free variant is as
follows [26]:

${\bf x}_{k+1}={\bf
x}_{k}+\frac{\alpha}{2q}\sum_{i\in{\mathcal{T}}_{k}}\frac{b_{i}-\langle
A_{i*}|{\bf x}_{k}\rangle}{\|A_{i*}\|^{2}}A_{i*},$ (9)

where ${\mathcal{T}}_{k}$ is a random set of $q$ row indices sampled with
replacement, $\alpha$ is known as the relaxation parameter. Each index
$i\in[m]$ is put into ${\mathcal{T}}_{k}$ with probability proportional to
$\|A_{i*}\|^{2}$. For instance, let $\alpha=2$, then (9) can be rewritten in a
parallel form

${\bf x}_{k+1}=\frac{1}{q}\sum_{i\in{\mathcal{T}}_{k}}\left({\bf
x}_{k}+\frac{b_{i}-\langle A_{i*}|{\bf
x}_{k}\rangle}{\|A_{i*}\|^{2}}A_{i*}\right).$

It was proved in [26, Corollary 3] that if $\alpha=q=\|A\|_{F}^{2}/\|A\|^{2}$
(also see a related proof in Appendix B),

${\mathbb{E}}[\|{\bf x}_{T}-{\bf
x}_{*}\|^{2}]\leq\left(1-\kappa^{-2}\right)^{T}\|{\bf x}_{0}-{\bf
x}_{*}\|^{2}.$ (10)

So the number of iterations required to achieve error $\epsilon$ via (9) is
$O(\kappa^{2}\log(1/\epsilon))$.

##  3 Preliminaries of quantum-inspired algorithms

This section briefly recalls some definitions and results about quantum-
inspired algorithms that will be used in this paper. The main reference is
[6].

######  Definition 3 (Query access).

For a vector ${\bf v}\in\mathbb{C}^{n}$, we have $Q({\bf v})$, query access to
${\bf v}$ if for all $i\in[n]$, we can query $v_{i}$. Likewise, for a matrix
$A\in\mathbb{C}^{m\times n}$, we have $Q(A)$ if for all
$(i,j)\in[m]\times[n]$, we can query $A_{ij}$. Let ${\bf q}({\bf v})$ (or
${\bf q}(A)$) denote the (time) cost of such a query.

######  Definition 4 (Sampling and query access to a vector).

For a vector ${\bf v}\in\mathbb{C}^{n}$, we have $SQ({\bf v})$, sampling and
query access to ${\bf v}$, if we can:

  1. 1.

Sample: obtain independent samples $i\in[n]$ following the distribution
$\mathcal{D}_{{\bf v}}\in\mathbb{R}^{n}$, where $\mathcal{D}_{{\bf
v}}(i)=|v_{i}|^{2}/\|{\bf v}\|^{2}$.

  2. 2.

Query: query entries of ${\bf v}$ as in $Q({\bf v})$;

  3. 3.

Norm: query $\|{\bf v}\|$.

Let ${\bf q}({\bf v}),{\bf s}({\bf v})$, and ${\bf n}({\bf v})$ denote the
cost of querying entries, sampling indices, and querying the norm
respectively. Further define ${\bf s}{\bf q}({\bf v}):={\bf q}({\bf v})+{\bf
s}({\bf v})+{\bf n}({\bf v})$.

######  Definition 5 (Sampling and query access to a matrix).

For any matrix $A\in\mathbb{C}^{m\times n}$, we have $SQ(A)$ if we have
$SQ(A_{i*})$ for all $i\in[m]$ and $SQ({\bf a})$, where ${\bf a}$ is the
vector of row norms. We define ${\bf s}{\bf q}(A)$ as the cost to obtain
$SQ(A)$.

The sampling and query access defined above can be viewed as a classical
analogue to the QRAM data structure, which is widely used in many quantum
machine learning (QML) algorithms. In QML, this data structure can be used to
efficiently prepare the quantum states of vectors, which are usually the
inputs of the QML algorithms. For matrices, it allows us to find an efficient
block-encoding so that many quantum linear algebraic techniques become
effective. However, the QRAM (and Definition 5) is a quite strong data
structure because it contains a lot of information (e.g. $\|A_{i*}\|$ and
$\|A\|_{F}$) about $A$, which is usually not easy to obtain. For instance, it
usually needs $\Theta(\|A\|_{0})$ operations to build this data structure,
which $\|A\|_{0}$ is the number of nonzero entries of $A$. For quantum-
inspired algorithms, we need to use certain sampling techniques to make sure
that the complexity has a polylog dependence on the dimension, so it seems
quite natural to use the above data structure. In comparison, there are also
some other efficient quantum models, e.g., sparse-access model, in which the
QML algorithms are still efficient.

The sampling operation on ${\bf v}$ can be viewed as a classical analogue of
the measurement of the quantum state $|{\bf v}\rangle=\|{\bf
v}\|^{-1}\sum_{i\in[n]}v_{i}|i\rangle$ in the computational basis. For some
problems, we do not have the state $|{\bf v}\rangle$ exactly. Instead, what we
have is a state of the form $\sin(\theta)|{\bf
v}\rangle|0\rangle+\cos(\theta)|{\bf w}\rangle|1\rangle$. Then the probability
of obtaining $|i,0\rangle$ is $\frac{v_{i}^{2}}{\|{\bf
v}\|^{2}}\sin^{2}(\theta)$. This leads to the definition of oversampling for
quantum-inspired algorithms.

######  Definition 6 (Oversampling and query access).

For ${\bf v}\in\mathbb{C}^{n}$ and $\phi\geq 1$, we have $SQ_{\phi}({\bf v})$,
$\phi$-oversampling and query access to ${\bf v}$, if we have $Q({\bf v})$ and
$SQ(\tilde{{\bf v}})$ for a vector $\tilde{{\bf v}}\in\mathbb{C}^{n}$
satisfying $\|\tilde{{\bf v}}\|^{2}=\phi\|{\bf v}\|^{2}$ and
$|\tilde{v}_{i}|^{2}\geq|v_{i}|^{2}$ for all $i\in[n]$. We denote ${\bf
s}_{\phi}(v)={\bf s}(\tilde{{\bf v}}),{\bf q}_{\phi}(v)={\bf q}(\tilde{{\bf
v}}),{\bf n}_{\phi}(v)={\bf n}(\tilde{{\bf v}})$, and ${\bf s}{\bf
q}_{\phi}({\bf v})={\bf s}_{\phi}(v)+{\bf q}_{\phi}(v)+{\bf q}(v)+{\bf
n}_{\phi}(v)$.

For example, if we view ${\bf v}$ as the vector $({\bf v},0,\ldots,0)$ and
$\tilde{{\bf v}}$ as the vector $({\bf v},\|{\bf v}\|\cot(\theta){\bf w})$ so
that they have the same dimension, then it is easy to see that $\|\tilde{{\bf
v}}\|^{2}=\|{\bf v}\|^{2}/\sin^{2}(\theta)$ and $\phi=1/\sin^{2}(\theta)$. The
condition $|\tilde{v}_{i}|^{2}\geq|v_{i}|^{2}$ is also satisfied. But the
above definition is stronger than this intuition.

######  Lemma 7 (Lemma 2.9 of [6]).

Given $SQ_{\phi}({\bf v})$ and $\delta\in(0,1]$ we can sample from
$\mathcal{D}_{{\bf v}}$ with success probability at least $1-\delta$ and cost
$O(\phi{\bf s}{\bf q}_{\phi}({\bf v})\log(1/\delta))$. We can also estimate
$\|{\bf v}\|$ up to relative error $\epsilon$ in time $O(\epsilon^{-2}\phi{\bf
s}{\bf q}_{\phi}({\bf v})\log(1/\delta))$.

######  Lemma 8 (Lemma 2.10 of [6]).

Given $SQ_{\varphi_{i}}({\bf v}_{i})$ and $\lambda_{i}\in\mathbb{C}$ for
$i\in[k]$, denote ${\bf v}=\sum_{i}\lambda_{i}{\bf v}_{i}$, then we have
$SQ_{\phi}({\bf v})$ for $\phi=k\sum_{i}\varphi_{i}\|\lambda_{i}{\bf
v}_{i}\|^{2}/\|{\bf v}\|^{2}$ and ${\bf s}{\bf q}_{\phi}({\bf v})=\max_{i}{\bf
s}_{\varphi_{i}}({\bf v}_{i})+\sum_{i}{\bf q}({\bf v}_{i})$.

The first lemma tells us how to construct the original sampling access from an
over-sampled one. By Definition 6, given $SQ_{\phi}({\bf v})$, we have $Q({\bf
v})$. Lemma 7 tells us how to obtain $S({\bf v})$ from $SQ_{\phi}({\bf v})$.
So this lemma actually tells us how to obtain $SQ({\bf v})$ from
$SQ_{\phi}({\bf v})$. The second lemma is about how to construct a linear
combination of the data structures. Equivalently, it is about constructing the
data structure for $M{\bf c}$ when we have the data structures of $M$ and
${\bf c}$, where $M$ is a matrix and ${\bf c}$ is a vector.

In our algorithms below, ${\bf v}_{i}$ will be a row of $A$, $k$ will refer to
the sparsity of a certain known vector and ${\bf v}$ will be the approximate
solution to the linear system. So in Lemma 8, ${\bf s}{\bf q}_{\phi}({\bf
v})=O(k{\bf s}{\bf q}(A))$. Combing the above two lemmas, given $SQ({\bf
v}_{i})$ and $\lambda_{i}$, the cost to sample and query ${\bf v}$ is
$\widetilde{O}(k\phi{\bf s}{\bf q}(A))$, and the cost to estimate its norm is
$\widetilde{O}(k\phi{\bf s}{\bf q}(A)/\epsilon^{2})$. So a main difficulty in
the the complexity analysis is the estimation of $\phi$.

##  4 The algorithm for solving linear systems

###  4.1 Solving general linear systems

In this section, based on the Kaczmarz method, we give a quantum-inspired
classical algorithm for solving linear systems. The basic idea is similar to
[11]. However, our result has a lower complexity and the analysis is much
simpler. The main idea is as follows: we first use the iteration (8) to
compute ${\bf y}_{T}$, then output $SQ({\bf x}_{T})$ with ${\bf
x}_{T}=A^{\dagger}{\bf y}_{T}$ by Lemmas 7 and 8.

When $A$ is dense, the calculation of $\langle A_{r_{k}*}|A^{\dagger}|{\bf
y}_{k}\rangle$ in (8) is expensive. Note that $\langle
A_{r_{k}*}|A^{\dagger}|{\bf y}_{k}\rangle=\sum_{j=1}^{n}A_{r_{k},j}\langle
A_{*j}|{\bf y}_{k}\rangle$. As an inner product, we can use the sampling idea
(e.g. Monte Carlo) to approximate it. Now let $j_{1},\ldots,j_{d}$ be drawn
randomly from $[n]$ under the distribution that ${\rm
Pr}[j]=\|A_{*j}\|^{2}/\|A\|_{F}^{2}$. This distribution is designed to put
more weight on the most important rows [9, 34]. Let
$S=\\{j_{1},\ldots,j_{d}\\}$ be the sequence of the chosen indices. We modify
the updating of ${\bf y}$ to:

${\bf y}_{k+1}={\bf
y}_{k}+\frac{1}{\|A_{r_{k}*}\|}\left(\tilde{b}_{r_{k}}-\frac{1}{d}\sum_{j\in
S}\tilde{A}_{r_{k},j}\langle A_{*j}|{\bf
y}_{k}\rangle\frac{\|A\|_{F}^{2}}{\|A_{*j}\|^{2}}\right){\bf e}_{r_{k}}.$ (11)

For conciseness, define the matrix $D$ as

$D=\sum_{j\in S}\frac{\|A\|_{F}^{2}}{d\|A_{*j}\|^{2}}|j\rangle\langle j|.$
(12)

Then

$\displaystyle{\bf y}_{k+1}$ $\displaystyle=$ $\displaystyle{\bf
y}_{k}+\frac{\tilde{b}_{r_{k}}-\langle\tilde{A}_{r_{k}*}|DA^{\dagger}|{\bf
y}_{k}\rangle}{\|A_{r_{k}*}\|}{\bf e}_{r_{k}}$ (13) $\displaystyle=$
$\displaystyle{\bf
y}_{k}+\frac{\tilde{b}_{r_{k}}-\langle\tilde{A}_{r_{k}*}|A^{\dagger}|{\bf
y}_{k}\rangle}{\|A_{r_{k}*}\|}{\bf
e}_{r_{k}}+\frac{\langle\tilde{A}_{r_{k}*}|(I-D)A^{\dagger}|{\bf
y}_{k}\rangle}{\|A_{r_{k}*}\|}{\bf e}_{r_{k}}.$ (14)

Compared to the original iteration (8), the procedure (14) contains a
perturbation term. To analyze its convergence rate, we need to determine $d$.
Before that, we first analyze the properties of the following random variable:

$\mu_{k}:=\langle\tilde{A}_{r_{k}*}|(I-D)A^{\dagger}|{\bf
y}_{k}\rangle=\langle\tilde{A}_{r_{k}*}|I-D|{\bf
x}_{k}\rangle=\langle\tilde{A}_{r_{k}*}|{\bf
x}_{k}\rangle-\frac{1}{d}\sum_{j\in
S}\tilde{A}_{r_{k},j}x_{k,j}\frac{\|A\|_{F}^{2}}{\|A_{*j}\|^{2}}.$ (15)

######  Lemma 9.

The mean and the variance of $\mu_{k}$ satisfy
${\mathbb{E}}_{D}[\mu_{k}]=0,{\rm
Var}_{D}[\mu_{k}]\leq\frac{\|A\|_{F}^{2}\|{\bf
x}_{k}\|^{2}}{d\min_{j\in[n]}\|A_{*j}\|^{2}}.$

###### Proof.

The results follow from the standard calculation.

$\displaystyle{\mathbb{E}}_{D}[\mu_{k}]$ $\displaystyle=$
$\displaystyle\langle\tilde{A}_{r_{k}*}|{\bf
x}_{k}\rangle-{\mathbb{E}}_{D}\left[\frac{1}{d}\sum_{j\in
S}\tilde{A}_{r_{k},j}x_{k,j}\frac{\|A\|_{F}^{2}}{\|A_{*j}\|^{2}}\right]$
$\displaystyle=$ $\displaystyle\langle\tilde{A}_{r_{k}*}|{\bf
x}_{k}\rangle-{\mathbb{E}}_{D}\left[\tilde{A}_{r_{k},j}x_{k,j}\frac{\|A\|_{F}^{2}}{\|A_{*j}\|^{2}}\right]$
$\displaystyle=$ $\displaystyle\langle\tilde{A}_{r_{k}*}|{\bf
x}_{k}\rangle-\sum_{j=1}^{n}\tilde{A}_{r_{k},j}x_{k,j}\frac{\|A\|_{F}^{2}}{\|A_{*j}\|^{2}}\frac{\|A_{*j}\|^{2}}{\|A\|_{F}^{2}}$
$\displaystyle=$ $\displaystyle 0.$ $\displaystyle{\rm Var}_{D}[\mu_{k}]$
$\displaystyle=$ $\displaystyle{\rm Var}_{D}\left[\frac{1}{d}\sum_{j\in
S}\tilde{A}_{r_{k},j}x_{k,j}\frac{\|A\|_{F}^{2}}{\|A_{*j}\|^{2}}\right]$
$\displaystyle=$ $\displaystyle\sum_{j\in S}{\rm
Var}_{D}\left[\frac{1}{d}\tilde{A}_{r_{k},j}x_{k,j}\frac{\|A\|_{F}^{2}}{\|A_{*j}\|^{2}}\right]$
$\displaystyle\leq$
$\displaystyle\frac{1}{d}{\mathbb{E}}_{D}\left[\left(\tilde{A}_{r_{k},j}x_{k,j}\frac{\|A\|_{F}^{2}}{\|A_{*j}\|^{2}}\right)^{2}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{d}\sum_{j=1}^{n}\left(\tilde{A}_{r_{k},j}x_{k,j}\frac{\|A\|_{F}^{2}}{\|A_{*j}\|^{2}}\right)^{2}\frac{\|A_{*j}\|^{2}}{\|A\|_{F}^{2}}$
$\displaystyle=$
$\displaystyle\frac{1}{d}\sum_{j=1}^{n}\tilde{A}_{r_{k},j}^{2}x_{k,j}^{2}\frac{\|A\|_{F}^{2}}{\|A_{*j}\|^{2}}$
$\displaystyle\leq$
$\displaystyle\frac{1}{d}\frac{\|A\|_{F}^{2}}{\min_{j\in[n]}\|A_{*j}\|^{2}}\|{\bf
x}_{k}\|^{2}.$

The last inequality is caused by the fact that $\tilde{A}_{r_{k}*}$ is
normalized so that $\tilde{A}_{r_{k},j}^{2}\leq 1$. ∎

From (14), the updating of ${\bf x}_{k}=A^{\dagger}{\bf y}_{k}$ then becomes

${\bf x}_{k+1}={\bf x}_{k}+(\tilde{b}_{r_{k}}-\langle\tilde{A}_{r_{k}*}|{\bf
x}_{k}\rangle)\tilde{A}_{r_{k}*}+\langle\tilde{A}_{r_{k}*}|I-D|{\bf
x}_{k}\rangle\tilde{A}_{r_{k}*}.$ (16)

Next, we prove a similar result to Lemma 2. The idea of our proof is similar
to that of Lemma 2 from [34], which nicely utilizes the geometrical structure
of the Kaczmarz method.

######  Proposition 10.

Choose $d=\frac{4\|A\|_{F}^{2}\kappa_{F}^{2}(\log
2/\epsilon^{2})}{\epsilon^{2}\min_{j\in[n]}\|A_{*j}\|^{2}}$, then after
$T=\kappa_{F}^{2}(\log 2/\epsilon^{2})$ steps of iteration of the procedure
(16), we have ${\mathbb{E}}[\|{\bf x}_{T}-{\bf
x}_{*}\|^{2}]\leq\epsilon^{2}\|{\bf x}_{*}\|^{2}$.

###### Proof.

For any vector ${\bf v}$ such that $A{\bf v}\neq 0$ we have $\|{\bf
v}\|^{2}\leq\|A^{-1}\|^{2}\|A{\bf v}\|^{2}$. This implies

$\frac{\|{\bf v}\|^{2}}{\kappa_{F}^{2}}=\frac{\|{\bf
v}\|^{2}}{\|A\|_{F}^{2}\|A^{-1}\|^{2}}\leq\sum_{i=1}^{m}\langle\tilde{A}_{i*}|{\bf
v}\rangle^{2}\frac{\|A_{i*}\|^{2}}{\|A\|_{F}^{2}}={\mathbb{E}}[\langle\tilde{A}_{i*}|{\bf
v}\rangle^{2}].$ (17)

Here the random variable is defined by

${\rm Pr}[i]=\frac{\|A_{i*}\|^{2}}{\|A\|_{F}^{2}}.$

For simplicity, we set the initial vector ${\bf x}_{0}=0$. In procedure (16),
we denote $\tilde{{\bf x}}_{k+1}={\bf
x}_{k}+(\tilde{b}_{r_{k}}-\langle\tilde{A}_{r_{k}*}|{\bf
x}_{k}\rangle)\tilde{A}_{r_{k}*}$, it is the orthogonal projection of ${\bf
x}_{k}$ onto the hyperplane $\tilde{A}_{r_{k}*}\cdot{\bf
x}=\tilde{b}_{r_{k}}$. Thus we can rewrite the iteration as

${\bf x}_{k+1}=\tilde{{\bf x}}_{k+1}+\mu_{k}\tilde{A}_{r_{k}*}.$

The orthogonality implies (see Figure 2)

$\displaystyle\|{\bf x}_{k+1}-{\bf x}_{*}\|^{2}$ $\displaystyle=$
$\displaystyle\|\tilde{{\bf x}}_{k+1}-{\bf x}_{*}\|^{2}+\mu_{k}^{2}$
$\displaystyle=$ $\displaystyle\|{\bf x}_{k}-{\bf x}_{*}\|^{2}-\|{\bf
x}_{k}-\tilde{{\bf x}}_{k+1}\|^{2}+\mu_{k}^{2}$ $\displaystyle=$
$\displaystyle\|{\bf x}_{k}-{\bf x}_{*}\|^{2}-\langle{\bf x}_{k}-{\bf
x}_{*}|\tilde{A}_{r_{k}*}\rangle^{2}+\mu_{k}^{2}$ $\displaystyle=$
$\displaystyle\left(1-\langle\frac{{\bf x}_{k}-{\bf x}_{*}}{\|{\bf x}_{k}-{\bf
x}_{*}\|}|\tilde{A}_{r_{k}*}\rangle^{2}\right)\|{\bf x}_{k}-{\bf
x}_{*}\|^{2}+\mu_{k}^{2}.$

${\bf x}_{*}$$\tilde{A}_{r_{k}*}\cdot{\bf x}=\tilde{b}_{r_{k}}$${\bf
x}_{k}$$\tilde{{\bf x}}_{k+1}$${\bf x}_{k+1}$$\mu_{k}$ Figure 2: Illustration
of the orthogonality in the proof of Proposition 10.

Taking the expectation yields

${\mathbb{E}}[\|{\bf x}_{k+1}-{\bf
x}_{*}\|^{2}]=\left(1-{\mathbb{E}}\left[\langle\frac{{\bf x}_{k}-{\bf
x}_{*}}{\|{\bf x}_{k}-{\bf
x}_{*}\|}|\tilde{A}_{r_{k}*}\rangle^{2}\right]\right){\mathbb{E}}[\|{\bf
x}_{k}-{\bf x}_{*}\|^{2}]+{\mathbb{E}}[\mu_{k}^{2}].$

In the above, the expectation is taken first on $D$ by fixing $k,r_{k}$, then
on $r_{k}$ by fixing $k$, and finally on $k$. Also if $A({\bf x}_{k}-{\bf
x}_{*})=0$, then ${\bf x}_{k}$ is already the optimal solution, so we assume
that this does not happen before convergence.

By equation (17),

$1-{\mathbb{E}}\left[\langle\frac{{\bf x}_{k}-{\bf x}_{*}}{\|{\bf x}_{k}-{\bf
x}_{*}\|}|\tilde{A}_{r_{k}*}\rangle^{2}\right]\leq 1-\kappa_{F}^{-2}.$

By Lemma 9, we obtain

$\displaystyle{\mathbb{E}}[\|{\bf x}_{k+1}-{\bf x}_{*}\|^{2}]$
$\displaystyle\leq$
$\displaystyle\left(1-\kappa_{F}^{-2}\right){\mathbb{E}}[\|{\bf x}_{k}-{\bf
x}_{*}\|^{2}]+\frac{\|A\|_{F}^{2}}{d\min_{j\in[n]}\|A_{*j}\|^{2}}{\mathbb{E}}[\|{\bf
x}_{k}\|^{2}]$ $\displaystyle\leq$
$\displaystyle\left(1-\kappa_{F}^{-2}\right){\mathbb{E}}[\|{\bf x}_{k}-{\bf
x}_{*}\|^{2}]+\frac{2\|A\|_{F}^{2}}{d\min_{j\in[n]}\|A_{*j}\|^{2}}(\|{\bf
x}_{*}\|^{2}+{\mathbb{E}}[\|{\bf x}_{k}-{\bf x}_{*}\|^{2}])$ $\displaystyle=$
$\displaystyle\left(1-\kappa_{F}^{-2}+\frac{2\|A\|_{F}^{2}}{d\min_{j\in[n]}\|A_{*j}\|^{2}}\right){\mathbb{E}}[\|{\bf
x}_{k}-{\bf
x}_{*}\|^{2}]+\frac{2\|A\|_{F}^{2}}{d\min_{j\in[n]}\|A_{*j}\|^{2}}\|{\bf
x}_{*}\|^{2}.$

Now we choose $d=\frac{4\|A\|_{F}^{2}\kappa_{F}^{2}(\log
2/\epsilon^{2})}{\epsilon^{2}\min_{j\in[n]}\|A_{*j}\|^{2}}$, then the above
estimation leads to

${\mathbb{E}}[\|{\bf x}_{k+1}-{\bf
x}_{*}\|^{2}]\leq\left(1-\kappa_{F}^{-2}+\frac{1}{2}\kappa_{F}^{-2}\epsilon^{2}(\log
2/\epsilon^{2})^{-1}\right){\mathbb{E}}[\|{\bf x}_{k}-{\bf
x}_{*}\|^{2}]+\frac{\epsilon^{2}}{2\kappa_{F}^{2}(\log 2/\epsilon^{2})}\|{\bf
x}_{*}\|^{2}.$

Setting $\rho:=1-\kappa_{F}^{-2}+\frac{1}{2}\kappa_{F}^{-2}\epsilon^{2}(\log
2/\epsilon^{2})^{-1}$ which is less than 1, we finally obtain

$\displaystyle{\mathbb{E}}[\|{\bf x}_{T}-{\bf x}_{*}\|^{2}]$
$\displaystyle\leq$ $\displaystyle\rho^{T}{\mathbb{E}}[\|{\bf x}_{0}-{\bf
x}_{*}\|^{2}]+\frac{\epsilon^{2}}{2\kappa_{F}^{2}(\log 2/\epsilon^{2})}\|{\bf
x}_{*}\|^{2}\sum_{i=0}^{T-1}\rho^{i}$ $\displaystyle\leq$
$\displaystyle\rho^{T}\|{\bf
x}_{*}\|^{2}+\frac{T\epsilon^{2}}{2\kappa_{F}^{2}(\log 2/\epsilon^{2})}\|{\bf
x}_{*}\|^{2}.$

Let $T=\kappa_{F}^{2}(\log 2/\epsilon^{2})$, then

$\displaystyle\rho^{T}$ $\displaystyle=$
$\displaystyle\exp(\kappa_{F}^{2}(\log
2/\epsilon^{2})\log(1-\kappa_{F}^{-2}+\kappa_{F}^{-2}\epsilon^{2}(\log
1/\epsilon)^{-1}))$ $\displaystyle\approx$ $\displaystyle\exp(-(\log
2/\epsilon^{2})(1-\epsilon^{2}(\log 1/\epsilon)^{-1})))$
$\displaystyle\approx$ $\displaystyle\epsilon^{2}/2.$

Therefore, ${\mathbb{E}}[\|{\bf x}_{T}-{\bf
x}_{*}\|^{2}]\leq\epsilon^{2}\|{\bf x}_{*}\|^{2}$ as claimed. ∎

######  Theorem 11.

Assume that ${\bf x}_{*}=\arg\min\|A{\bf x}-{\bf b}\|$. Given $SQ(A),Q({\bf
b})$, there is an algorithm that returns $SQ({\bf x})$ such that $\|{\bf
x}-{\bf x}_{*}\|\leq\epsilon\|{\bf x}_{*}\|$ with probability at least 0.99
and

${\bf s}{\bf q}({\bf
x})=\widetilde{O}\left(\kappa_{F}^{6}\left(\frac{\|A\|_{F}^{2}({\bf s}{\bf
q}(A)+{\bf q}({\bf
b}))}{\epsilon^{2}\min_{j\in[n]}\|A_{*j}\|^{2}}+\kappa^{2}{\bf s}{\bf
q}(A)\right)\right)=\widetilde{O}(\kappa_{F}^{8}({\bf s}{\bf q}(A)+{\bf
q}({\bf b}))/\epsilon^{2}).$ (18)

###### Proof.

Let $d=\frac{4\|A\|_{F}^{2}\kappa_{F}^{2}(\log
2/\epsilon^{2})}{\epsilon^{2}\min_{j\in[n]}\|A_{*j}\|^{2}}$. By the updating
formula (11), we know that ${\bf y}_{k}$ is at most $k$-sparse. So in the
$k$-th step of the iteration, the calculation of the summation in (11) costs
$O(kd)$. At each step of the iteration, we sample a row of $A$ and query the
corresponding entry of ${\bf b}$. So the overall cost of $k$-th step is
$O(kd({\bf s}{\bf q}(A)+{\bf q}({\bf b})))$. After convergence, it costs
$O(T^{2}d({\bf s}{\bf q}(A)+{\bf q}({\bf b})))$ in total, where
$T=O(\kappa_{F}^{2}\log(1/\epsilon))$.

Note that ${\bf x}_{T}=A^{\dagger}{\bf y}_{T}$ and ${\bf y}_{T}$ is
$T$-sparse. By Lemma 8 (taking ${\bf v}_{i}$ to be rows of $A$), we can obtain
$SQ_{\phi}({\bf x}_{T})$ in time $O(T{\bf s}{\bf q}(A))$, where

$\phi=T\frac{\sum_{i}y_{T,i}^{2}\|A_{i*}\|^{2}}{\|{\bf x}_{T}\|^{2}}.$

In Appendix A, we shall show that $\phi=O(T^{2}\kappa^{2})$. By Lemma 7, we
can obtain $SQ({\bf x}_{T})$ in time $O(\phi T{\bf s}{\bf
q}(A))=O(T^{3}\kappa^{2}{\bf s}{\bf q}(A)).$ Therefore, the total cost is
$O(T^{2}d({\bf s}{\bf q}(A)+{\bf q}({\bf b}))+T^{3}\kappa^{2}{\bf s}{\bf
q}(A))$. Substituting the values of $d,T$ into it gives rise to the claimed
result. ∎

We next use the Kaczmarz method with averaging (9) to reduce the dependence on
$\kappa_{F}$ in the above theorem. Similar to (16), we change the updating
rule into

${\bf x}_{k+1}={\bf
x}_{k}+\frac{1}{2}\sum_{i\in{\mathcal{T}}_{k}}(\tilde{b}_{i}-\langle\tilde{A}_{i*}|{\bf
x}_{k}\rangle)\tilde{A}_{i*}+\frac{1}{2}\sum_{i\in{\mathcal{T}}_{k}}\langle\tilde{A}_{i*}|I-D_{i}|{\bf
x}_{k}\rangle\tilde{A}_{i*},$ (19)

The definition of $D_{i}$ is similar to (12), which depends on $i$ now. In
Appendix B, we shall prove that if

$d=\frac{\|A\|_{F}^{4}\kappa^{2}\log(1/\epsilon^{2})}{\epsilon^{2}\|A\|^{2}\min_{j\in[n]}\|A_{*j}\|^{2}}=\widetilde{O}(\kappa_{F}^{4}/\epsilon^{2}),$
(20)

then after $T=O(\kappa^{2}\log(1/\epsilon))$ steps of iteration, we have
${\mathbb{E}}[\|{\bf x}_{T}-{\bf x}_{*}\|^{2}]\leq\epsilon^{2}\|{\bf
x}_{0}-{\bf x}_{*}\|^{2}$. Together with the estimation of $\phi$ in Appendix
C, we have the following improved result.

######  Theorem 12.

Assume that ${\bf x}_{*}=\arg\min\|A{\bf x}-{\bf b}\|$. Given $SQ(A),Q({\bf
b})$, there is an algorithm that returns $SQ({\bf x})$ such that $\|{\bf
x}-{\bf x}_{*}\|\leq\epsilon\|{\bf x}_{*}\|$ with probability at least 0.99
and

${\bf s}{\bf q}({\bf x})=\widetilde{O}\left({\kappa_{F}^{6}\kappa^{2}({\bf
s}{\bf q}(A)+{\bf q}({\bf b}))}/{\epsilon^{2}}\right).$ (21)

###### Proof.

From the updating formula (19), the updating of ${\bf y}_{k}$ is

${\bf y}_{k+1}={\bf
y}_{k}+\frac{1}{2}\sum_{i\in{\mathcal{T}}_{k}}\frac{\tilde{b}_{i}-\langle\tilde{A}_{i*}|D_{i}A^{\dagger}|{\bf
y}_{k}\rangle}{\|A_{i*}\|}{\bf e}_{i}.$

So ${\bf y}_{k}$ is at most $(k\kappa_{F}^{2}/\kappa^{2})$-sparse. Thus the
cost in step $k$ is $O(kd\kappa_{F}^{2}({\bf s}{\bf q}(A)+{\bf q}({\bf
b}))/\kappa^{2})=\widetilde{O}(k\kappa_{F}^{6}({\bf s}{\bf q}(A)+{\bf q}({\bf
b}))/\kappa^{2}\epsilon^{2})$. After convergence, the total cost is
$\widetilde{O}(T^{2}\kappa_{F}^{6}({\bf s}{\bf q}(A)+{\bf q}({\bf
b}))/\kappa^{2}\epsilon^{2})=\widetilde{O}(\kappa_{F}^{6}\kappa^{2}({\bf
s}{\bf q}(A)+{\bf q}({\bf b}))/\epsilon^{2})$. By the estimation in Appendix
C, we have $\phi=O(\kappa_{F}^{2}\kappa^{4})$, so the cost to obtain $SQ({\bf
x}_{T})$ is $O(\phi T\kappa_{F}^{2}{\bf s}{\bf
q}(A)/\kappa^{2})=\widetilde{O}(\kappa_{F}^{4}\kappa^{4}{\bf s}{\bf q}(A))$.
Hence, the overall cost is $\widetilde{O}(\kappa_{F}^{6}\kappa^{2}({\bf s}{\bf
q}(A)+{\bf q}({\bf b}))/\epsilon^{2})$. ∎

###  4.2 Solving symmetric positive definite linear systems

Now suppose $A$ is a symmetric positive definite matrix. Then a simple
approach to solve $A{\bf x}={\bf b}$ is the randomized coordinate descent
iteration [15]

${\bf x}_{k+1}={\bf x}_{k}-\frac{\langle A_{r_{k}*}|{\bf x}_{k}\rangle-
b_{r_{k}}}{A_{r_{k},r_{k}}}{\bf e}_{r_{k}},$ (22)

where the probability of choosing $r_{k}$ is $A_{r_{k},r_{k}}/{\rm Tr}(A)$. If
we apply the above iteration to the normal equation $A^{\dagger}A{\bf
x}=A^{\dagger}{\bf b}$, we will obtain the Kaczmarz method. As for the
iteration (22), it is not hard to show that

${\mathbb{E}}[\|{\bf x}_{k}-{\bf
x}_{*}\|^{2}]\leq\left(1-\frac{1}{\|A^{-1}\|{\rm Tr(A)}}\right)^{k}\|{\bf
x}_{0}-{\bf x}_{*}\|^{2}.$

In fact, assume $A{\bf x}_{*}={\bf b}$, then (22) implies that

${\bf x}_{k+1}-{\bf x}_{*}=\left(I-\frac{|r_{k}\rangle\langle
r_{k}|A}{A_{r_{k},r_{k}}}\right)({\bf x}_{k}-{\bf x}_{*}).$

Thus

$\displaystyle{\mathbb{E}}[\|{\bf x}_{k+1}-{\bf x}_{*}\|^{2}]$
$\displaystyle\leq$
$\displaystyle\left\|{\mathbb{E}}\left[I-\frac{|r_{k}\rangle\langle
r_{k}|A}{A_{r_{k},r_{k}}}\right]\right\|{\mathbb{E}}[\|{\bf x}_{k}-{\bf
x}_{*}\|^{2}]$ $\displaystyle=$ $\displaystyle\left\|I-\frac{A}{{\rm
Tr}(A)}\right\|{\mathbb{E}}[\|{\bf x}_{k}-{\bf x}_{*}\|^{2}]$ $\displaystyle=$
$\displaystyle\left(1-\frac{1}{\|A^{-1}\|{\rm
Tr(A)}}\right){\mathbb{E}}[\|{\bf x}_{k}-{\bf x}_{*}\|^{2}].$

Similar to the idea introduced above to handle Kaczmarz method, we can use the
randomized coordinate descent method to solve dense linear systems $A{\bf
x}={\bf b}$ such that $A$ is symmetric positive definite. The only difference
is that $\|A\|_{F}^{2}$ becomes ${\rm Tr}(A)$, and $\|A^{-1}\|^{2}$ becomes
$\|A^{-1}\|$. Moreover, the randomized coordinate descent (22) can also be
parallelized into

${\bf x}_{k+1}={\bf
x}_{k}-\frac{\alpha}{2q}\sum_{i\in{\mathcal{T}}_{k}}\frac{\langle A_{i*}|{\bf
x}_{k}\rangle-b_{i}}{A_{i,i}}{\bf e}_{i}.$ (23)

If we take $\alpha=q={\rm Tr}(A)/\|A\|$, then it converges after
$O(\kappa\log(1/\epsilon))$ iterations. The proof is the same as that of (9).
Therefore, by a similar argument to the proof of Theorem 12, we have

######  Theorem 13.

Let $A$ be symmetric positive definite. Assume that ${\bf
x}_{*}=\arg\min\|A{\bf x}-{\bf b}\|$. Given $SQ(A),Q({\bf b})$, there is an
algorithm that returns $SQ({\bf x})$ such that $\|{\bf x}-{\bf
x}_{*}\|\leq\epsilon\|{\bf x}_{*}\|$ with probability at least 0.99 and

${\bf s}{\bf q}({\bf x})=\widetilde{O}\left(\frac{{\rm
Tr}(A)^{3}\|A^{-1}\|^{3}\kappa({\bf s}{\bf q}(A)+{\bf q}({\bf
b}))}{\epsilon^{2}}\right).$ (24)

The above result is actually more general than Theorem 12. If we consider the
normal equation $A^{\dagger}A{\bf x}=A^{\dagger}{\bf b}$, Theorem 13 implies
Theorem 12.

##  5 Discussion

###  5.1 Solving row sparse linear systems

When $A$ is row sparse, the inner product in the Kaczmarz method (2) is not
expensive to calculate. As a result, we can directly apply the Kaczmarz method
to solve the linear system $A{\bf x}={\bf b}$. The result is summarized below,
and the proof is straightforward.

######  Theorem 14.

Assume that $A$ has row sparsity $s$. Let ${\bf x}_{*}=\arg\min\|A{\bf x}-{\bf
b}\|$. Given $SQ(A),SQ({\bf b})$, then there is an algorithm that returns an
$O(s\kappa_{F}^{2}\log(1/\epsilon))$-sparse vector ${\bf x}$ such that $\|{\bf
x}-{\bf x}_{*}\|\leq\epsilon\|{\bf x}_{*}\|$ with probability at least 0.99 in
time

$O(s\kappa_{F}^{2}\log(1/\epsilon)).$ (25)

###### Proof.

In the Kaczmarz iteration (2), since $A$ has row sparsity $s$, it follows that
the calculation of the inner product $\langle A_{r_{k}*}|{\bf x}_{k}\rangle$
costs $O(s)$ operations. This is also the cost of the $k$-th step of
iteration. After convergence, the total cost is $O(Ts)$, where
$T=O(\kappa_{F}^{2}\log(1/\epsilon))$. As for the sparsity of ${\bf x}_{k}$,
we can choose the initial approximation ${\bf x}_{0}=0$, then ${\bf x}_{k}$ is
at most $sk$-sparse. ∎

Similarly, based on the randomized coordinate descent method (22), we have

######  Theorem 15.

Suppose that $A$ is symmetric positive definite and $s$-sparse. Let ${\bf
x}_{*}=\arg\min\|A{\bf x}-{\bf b}\|$. Given $SQ(A),SQ({\bf b})$, then there is
an algorithm that returns an $O({\rm Tr}(A)\|A^{-1}\|\log(1/\epsilon))$-sparse
vector ${\bf x}$ such that $\|{\bf x}-{\bf x}_{*}\|\leq\epsilon\|{\bf
x}_{*}\|$ with probability at least 0.99 in time

$O(s{\rm Tr}(A)\|A^{-1}\|\log(1/\epsilon)).$ (26)

Differently from quantum and quantum-inspired algorithms, the outputs of the
above two randomized algorithms are vectors. It is known that if $A$ has row
and column sparsity $s$ and given in the sparse access input model, the
currently best known quantum algorithm for matrix inversion costs
$O(s\|A\|_{\max}\|A^{-1}\|\log^{2}(1/\epsilon))$, where $\|A\|_{\max}$ refers
to the maximal entry in the sense of absolute value [4]. The comparison
between this quantum algorithm and the randomized algorithm stated in Theorem
14 may be not fair enough because they are building on different models – the
quantum algorithm uses the sparse access model, while the later one assumes
access to $SQ(A)$. However, if we only concentrate on the final complexity,
the quantum linear solver achieves large speedups if
$\|A\|_{F}\gg\|A\|_{\max}$. This can often happen, especially when the rank of
$A$ is large (e.g. $A$ is a sparse unitary). However, when the rank of $A$ is
polylog in the dimension, then $\|A\|_{F}=\widetilde{O}(s\|A\|_{\max})$ if $A$
is row and column sparse. In this case, the quantum computer only achieves a
quadratic speedup. For a row sparse but column dense matrix $A$, if it is
stored in the QRAM data structure, the complexity to solve $A{\bf x}={\bf b}$
using a quantum computer is $O(\kappa_{F}\log^{2}(1/\epsilon))$ [4]. This only
gives a quadratic speedup in terms of $\kappa_{F}$ over Theorem 14.

By Theorem 15, when $A$ is symmetric positive definite, the randomized
coordinate descent method finds an approximate solution in time $O({\rm
Tr}(A)\|A^{-1}\|\log(1/\epsilon))$. For simplicity, let us consider the
special case that $A$ is a sparse density operator, then ${\rm Tr(A)}=1$.
Under certain conditions [12], we can efficiently find a unitary such that $A$
is a sub-block. This means that the complexity of the quantum linear solver
for this particular linear system is $O(\|A\|\|A^{-1}\|\log^{2}(1/\epsilon))$,
while the algorithm based on Theorem 15 costs $O(\|A^{-1}\|\log(1/\epsilon))$.
So if $\|A\|=\Theta(1)$, it is possible that there is no quantum speedup based
on the current quantum linear solvers.

In conclusion, quantum linear solvers can build on many different models, e.g.
sparse-access input model, QRAM model, etc. Whenever we have an efficient
block-encoding, quantum computers can solve the corresponding linear system
efficiently. In comparison, so far there is only one model for quantum-
inspired algorithms. The above two theorems indicate that to explore large
quantum speedups for sparse problems, it is better to focus on other models
rather than the QRAM model. This also leads to the problem of exploring other
models for quantum-inspired algorithms.

###  5.2 Reducing the dependence on $\|A\|_{F}$

The Frobenius norm can be large for high-rank matrices. For example, when we
use the discretization method (e.g. Euler method) to solve a linear system of
ordinary differential equations, the obtained linear system contains the
identity as a sub-matrix. So the Frobenius norm can be as large as the
dimension of the problem. The performance of quantum-inspired algorithms is
highly affected by the Frobenius norm, so currently they are mainly effective
for low-rank systems. To enlarge the application scope, it is necessary to
consider the problem of how to reduce or even remove the the dependence on the
Frobenius norm in the complexity.

To the best of our knowledge, there are no generalizations of the Kaczmarz
method such that the total cost is independent of $\|A\|_{F}$. However, it can
be generalized such that $\|{\mathbb{E}}[{\bf x}_{k}-{\bf
x}_{*}]\|\leq(1-\kappa^{-r})^{k}\|{\bf x}_{0}-{\bf x}_{*}\|$ and each step of
iteration is not expensive, where $r\in\\{1,2\\}$. For example, a simple
generalization is

${\bf x}_{k+1}={\bf
x}_{k}+\frac{\|A\|_{F}^{2}}{\|A\|^{2}}(\tilde{b}_{r_{k}}-\langle\widetilde{A}_{r_{k}*}|{\bf
x}_{k}\rangle)\widetilde{A}_{r_{k}*}.$

Generally, the relaxation parameter ${\|A\|_{F}^{2}}/{\|A\|^{2}}\gg 2$, so
this iterative scheme is not convergent,444To make sure that it is convergent,
the relaxation parameter should be chosen from $(0,2)$. Usually, the optimal
choice is 1. that is we usually do not have ${\mathbb{E}}[\|{\bf x}_{k}-{\bf
x}_{*}\|^{2}]\|\leq\rho^{k}\|{\bf x}_{0}-{\bf x}_{*}\|^{2}$ for some $\rho<1$.
However, we do have $\|{\mathbb{E}}[{\bf x}_{k}-{\bf
x}_{*}]\|\leq(1-\kappa^{-2})^{k}\|{\bf x}_{0}-{\bf x}_{*}\|$. The proof is
straightforward. If we set $T=O(\kappa^{2}\log(1/\epsilon))$, then
$\|{\mathbb{E}}[{\bf x}_{T}-{\bf x}_{*}]\|\leq\epsilon\|{\bf x}_{0}-{\bf
x}_{*}\|$. There are also some other generalizations [27, 32, 25] such that
$T=O(\kappa\log(1/\epsilon))$, which is optimal for this kind of iterative
schemes. One idea to apply this result is to generate $L$ samples ${\bf
x}_{T1},\ldots,{\bf x}_{TL}$ such that $\|{\mathbb{E}}[{\bf
x}_{T}]-\frac{1}{L}\sum_{l=1}^{L}{\bf x}_{Tl}\|\leq\epsilon\|{\mathbb{E}}[{\bf
x}_{T}]\|$. The matrix Bernstein inequality might be helpful to determine $L$.
However, to apply this method, we have to find tight bounds for $\|{\bf
x}_{T}-{\bf x}_{*}\|$ and ${\mathbb{E}}[\|{\bf x}_{T}-{\bf x}_{*}\|^{2}]$,
which seems not that easy.

##  6 Outlook

In this paper, for solving linear systems, we reduced the separation between
quantum and quantum-inspired algorithms from
$\kappa_{F}:\kappa_{F}^{6}\kappa^{6}/\epsilon^{4}$ to
$\kappa_{F}:\kappa_{F}^{6}\kappa^{2}/\epsilon^{2}$. It is interesting to
consider whether the parameters of our algorithm can be improved further. For
the class of iterative schemes including the Kaczmarz iteration, the optimal
convergence rate equals the condition number $\kappa$ (as opposed to
$\kappa_{F}^{2}$ as in the Kaczmarz method), which is achieved by the
conjugate gradient method. However, there is currently no such result for the
Kaczmarz iterative schemes themselves. As discussed in Section 5.2, in [27,
32, 25], it was proved that there are generalizations of the Kaczmarz method
so that $\|{\mathbb{E}}[{\bf x}_{T}-{\bf
x}_{*}]\|\leq(1-\kappa^{-1})^{T}\|{\bf x}_{0}-{\bf x}_{*}\|$. This convergence
result is a little weak but may be helpful to reduce the separation in terms
of $\kappa_{F}$ further.

## Acknowledgement

We would like to thank Ryan Mann for useful discussions. This paper was
supported by the QuantERA ERA-NET Cofund in Quantum Technologies implemented
within the European Union’s Horizon 2020 Programme (QuantAlgo project), and
EPSRC grants EP/L021005/1 and EP/R043957/1. This project has received funding
from the European Research Council (ERC) under the European Union’s Horizon
2020 research and innovation programme (grant agreement No. 817581). No new
data were created during this study.

##  Appendix A Estimation of $\phi$

Using the notation (15), we can decompose the updating rule (14) of ${\bf y}$
as follows

${\bf y}_{k+1}={\bf y}_{k}+{\bf z}_{k}+{\bf z}_{k}^{\prime},$ (27)

where

$\displaystyle{\bf
z}_{k}=\frac{\tilde{b}_{r_{k}}-\langle\tilde{A}_{r_{k}*}|A^{\dagger}|{\bf
y}_{k}\rangle}{\|A_{r_{k}*}\|}{\bf e}_{r_{k}},\quad{\bf
z}_{k}^{\prime}=\frac{\mu_{k}}{\|A_{r_{k}*}\|}{\bf e}_{r_{k}}.$

Denote

$Z=\|A{\bf x}_{*}-{\bf b}\|=\min_{{\bf x}}\|A{\bf x}-{\bf
b}\|,\quad\Lambda={\rm diag}(\|A_{i*}\|^{2}:i\in[m]).$ (28)

In the following, for any two vectors ${\bf a},{\bf b}$, we define
$\langle{\bf a}|{\bf b}\rangle_{\Lambda}=\langle{\bf a}|\Lambda|{\bf
b}\rangle$. To bound $\phi$, it suffices to bound $\|{\bf
y}_{T}\|_{\Lambda}^{2}$. From (27), it is plausible that $\|{\bf
y}_{k+1}\|_{\Lambda}^{2}\approx\|{\bf y}_{k}\|_{\Lambda}^{2}+\|{\bf
z}_{k}\|_{\Lambda}^{2}+\|{\bf z}_{k}^{\prime}\|_{\Lambda}^{2}$. In the
following, we shall prove that in fact this holds up to a constant factor on
average. We shall choose the initial vector as 0 for simplicity, i.e., ${\bf
y}_{0}=0$.

First, we consider the case $Z=0$, that is ${\bf b}=A{\bf x}_{*}$. In the
following, we first fix $k,r_{k}$ and compute the mean value over the random
variable $D$, then we compute the mean value over $r_{k}$ by fixing $k$.
Finally, we calculate the mean value over the random variable $k$.

From now on, we assume that
$d=\frac{4\|A\|_{F}^{2}T}{\epsilon^{2}\min_{j\in[n]}\|A_{*j}\|^{2}}$ and
$T=O(\kappa_{F}^{2}\log(1/\epsilon))$. By Lemma 9,

$\displaystyle{\mathbb{E}}_{D}[\|{\bf z}_{k}^{\prime}\|_{\Lambda}^{2}]$
$\displaystyle\leq$ $\displaystyle\frac{\|A\|_{F}^{2}\|{\bf
x}_{k}\|^{2}}{d\min_{j\in[n]}\|A_{*j}\|^{2}}=\frac{\epsilon^{2}\|{\bf
x}_{k}\|^{2}}{4T}\leq\frac{\epsilon^{2}(\|{\bf x}_{k}-{\bf x}_{*}\|^{2}+\|{\bf
x}_{*}\|^{2})}{2T},$ $\displaystyle{\mathbb{E}}_{D}[\langle{\bf z}_{k}|{\bf
z}_{k}^{\prime}\rangle_{\Lambda}]$ $\displaystyle=$
$\displaystyle\left(\tilde{b}_{r_{k}}-\langle\tilde{A}_{r_{k}*}|{\bf
x}_{k}\rangle\right){\mathbb{E}}_{D}[\mu_{k}]=0,$
$\displaystyle{\mathbb{E}}_{D}[\langle{\bf y}_{k}|{\bf
z}_{k}^{\prime}\rangle_{\Lambda}]$ $\displaystyle=$
$\displaystyle\|A_{r_{k}*}\|\langle{\bf y}_{k}|{\bf
e}_{r_{k}}\rangle{\mathbb{E}}_{D}[\mu_{k}]=0.$

About the norm of ${\bf z}_{k}$, we have

$\displaystyle{\mathbb{E}}[\|{\bf z}_{k}\|_{\Lambda}^{2}]$ $\displaystyle=$
$\displaystyle\sum_{r_{k}=1}^{m}\frac{{\|A_{r_{k}*}\|}^{2}}{\|A\|_{F}^{2}}\frac{(\tilde{b}_{r_{k}}-\langle\tilde{A}_{r_{k}*}|A^{\dagger}|{\bf
y}_{k}\rangle)^{2}}{\|A_{r_{k}*}\|^{2}}\|A_{r_{k}*}\|^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{\|A\|_{F}^{2}}\sum_{r_{k}=1}^{m}(b_{r_{k}}-\langle
A_{r_{k}}|{\bf x}_{k}\rangle)^{2}$ $\displaystyle\leq$ $\displaystyle
2\frac{\|{\bf b}\|^{2}+\|A\|_{F}^{2}\|{\bf x}_{k}\|^{2}}{\|A\|_{F}^{2}}\leq
2\frac{\|{\bf b}\|^{2}}{\|A\|_{F}^{2}}+4\|{\bf x}_{*}\|^{2}+4\|{\bf
x}_{*}-{\bf x}_{k}\|^{2}.$

As for the inner product between ${\bf y}_{k}$ and ${\bf z}_{k}$, we have the
following estimate

$\displaystyle{\mathbb{E}}[\langle{\bf y}_{k}|{\bf z}_{k}\rangle_{\Lambda}]$
$\displaystyle=$
$\displaystyle\sum_{r_{k}=1}^{m}\frac{{\|A_{r_{k}*}\|}^{2}}{\|A\|_{F}^{2}}\frac{(\tilde{b}_{r_{k}}-\langle\tilde{A}_{r_{k}*}|A^{\dagger}|{\bf
y}_{k}\rangle)}{\|A_{r_{k}*}\|}\langle{\bf y}_{k}|{\bf
e}_{r_{k}}\rangle\|A_{r_{k}*}\|^{2}$ $\displaystyle=$
$\displaystyle\sum_{r_{k}=1}^{m}\frac{{\|A_{r_{k}*}\|}^{2}}{\|A\|_{F}^{2}}{(b_{r_{k}}-\langle
A_{r_{k}*}|A^{\dagger}|{\bf y}_{k}\rangle)}\langle{\bf y}_{k}|{\bf
e}_{r_{k}}\rangle$ $\displaystyle=$ $\displaystyle\frac{\langle{\bf b}|{\bf
y}_{k}\rangle_{\Lambda}-\|{\bf
x}_{k}\|_{\Lambda}^{2}}{\|A\|_{F}^{2}}=\frac{\langle{\bf
x}_{*}|A^{\dagger}|{\bf y}_{k}\rangle_{\Lambda}-\|{\bf
x}_{k}\|_{\Lambda}^{2}}{\|A\|_{F}^{2}}$ $\displaystyle=$
$\displaystyle\frac{\langle{\bf x}_{*}|{\bf x}_{k}\rangle_{\Lambda}-\|{\bf
x}_{k}\|_{\Lambda}^{2}}{\|A\|_{F}^{2}}\leq\|{\bf x}_{*}\|\|{\bf
x}_{k}\|+\|{\bf x}_{k}\|^{2}$ $\displaystyle\leq$ $\displaystyle 3\|{\bf
x}_{*}\|^{2}+\|{\bf x}_{*}\|\|{\bf x}_{k}-{\bf x}_{*}\|+2\|{\bf x}_{k}-{\bf
x}_{*}\|^{2}$ $\displaystyle\leq$ $\displaystyle\frac{7}{2}\|{\bf
x}_{*}\|^{2}+\frac{5}{2}\|{\bf x}_{k}-{\bf x}_{*}\|^{2}.$

In the above, we used the fact that for any two vectors ${\bf a},{\bf b}$ we
have $|\langle{\bf a}|{\bf b}\rangle_{\Lambda}|\leq\|\Lambda\|^{2}|\langle{\bf
a}|{\bf b}\rangle|$, and $\|\Lambda\|^{2}\leq\|A\|_{F}^{2}$.

Hence, we have

$\displaystyle{\mathbb{E}}[\|{\bf y}_{k+1}\|_{\Lambda}^{2}]$ $\displaystyle=$
$\displaystyle{\mathbb{E}}[\|{\bf y}_{k}\|_{\Lambda}^{2}]+{\mathbb{E}}[\|{\bf
z}_{k}\|_{\Lambda}^{2}]+{\mathbb{E}}[\|{\bf
z}_{k}^{\prime}\|_{\Lambda}^{2}]+2{\mathbb{E}}[\langle{\bf y}_{k}|{\bf
z}_{k}\rangle_{\Lambda}]$ $\displaystyle\leq$
$\displaystyle{\mathbb{E}}[\|{\bf y}_{k}\|_{\Lambda}^{2}]+2\frac{\|{\bf
b}\|^{2}}{\|A\|_{F}^{2}}+(9+\frac{\epsilon^{2}}{2T}){\mathbb{E}}[\|{\bf
x}_{k}-{\bf x}_{*}\|^{2}]+(11+\frac{\epsilon^{2}}{2T})\|{\bf x}_{*}\|^{2}$
$\displaystyle\leq$ $\displaystyle{\mathbb{E}}[\|{\bf
y}_{k}\|_{\Lambda}^{2}]+2\frac{\|{\bf
b}\|^{2}}{\|A\|_{F}^{2}}+(20+\frac{\epsilon^{2}}{T})\|{\bf x}_{*}\|^{2},$

where we use that ${\mathbb{E}}[\|{\bf x}_{k}-{\bf x}_{*}\|^{2}]\leq\|{\bf
x}_{*}\|^{2}$ by Lemma 2. Therefore,

$\displaystyle{\mathbb{E}}[\|{\bf y}_{T}\|_{\Lambda}^{2}]\leq
T\left(2\frac{\|{\bf
b}\|^{2}}{\|A\|_{F}^{2}}+(20+\frac{\epsilon^{2}}{T})\|{\bf
x}_{*}\|^{2}\right).$

This means that, with high probability,

$\phi=T\frac{\|{\bf y}_{T}\|_{\Lambda}^{2}}{\|{\bf x}_{T}\|^{2}}\leq
T^{2}\left(2\frac{\|{\bf b}\|^{2}}{\|A\|_{F}^{2}\|{\bf
x}_{T}\|^{2}}+(20+\frac{\epsilon^{2}}{T})\frac{\|{\bf x}_{*}\|^{2}}{\|{\bf
x}_{T}\|^{2}}\right)=O(T^{2}).$

When $Z\neq 0$, then ${\bf b}=A{\bf x}_{*}+{\bf c}$ for some vector ${\bf c}$
of norm $Z$ that is not in the column space of $A$. This can happen when $A$
is not full rank. Since ${\bf c}$ is independent of $A$, we cannot bound it in
terms of $A$. In this case, the only change is ${\mathbb{E}}[\langle{\bf
y}_{k}|{\bf z}_{k}\rangle_{\Lambda}]$, which is now bounded by

$\displaystyle{\mathbb{E}}[\langle{\bf y}_{k}|{\bf z}_{k}\rangle_{\Lambda}]$
$\displaystyle\leq$ $\displaystyle\frac{\langle{\bf c}|{\mathbb{E}}[{\bf
y}_{k}]\rangle_{\Lambda}}{\|A\|_{F}^{2}}+\frac{7}{2}\|{\bf
x}_{*}\|^{2}+\frac{5}{2}\|{\bf x}_{k}-{\bf x}_{*}\|^{2}$ $\displaystyle\leq$
$\displaystyle\frac{\|A\|^{2}}{\|A\|_{F}^{2}}\|{\bf c}\|\|{\mathbb{E}}[{\bf
y}_{k}]\|+\frac{7}{2}\|{\bf x}_{*}\|^{2}+\frac{5}{2}\|{\bf x}_{k}-{\bf
x}_{*}\|^{2}$ $\displaystyle=$
$\displaystyle\frac{\kappa^{2}Z}{\kappa_{F}^{2}}\|{\mathbb{E}}[{\bf
y}_{k}]\|+\frac{7}{2}\|{\bf x}_{*}\|^{2}+\frac{5}{2}\|{\bf x}_{k}-{\bf
x}_{*}\|^{2}.$

At the end, we obtain

$\displaystyle{\mathbb{E}}[\|{\bf
y}_{k+1}\|_{\Lambda}^{2}]\leq{\mathbb{E}}[\|{\bf
y}_{k}\|_{\Lambda}^{2}]+2\frac{\|{\bf
b}\|^{2}}{\|A\|_{F}^{2}}+(20+\frac{\epsilon^{2}}{T})\|{\bf
x}_{*}\|^{2}+\frac{2\kappa^{2}Z}{\kappa_{F}^{2}}\|{\mathbb{E}}[{\bf
y}_{k}]\|.$

From (27), we know that

${\mathbb{E}}[{\bf
y}_{k+1}]=\left(I-\frac{AA^{\dagger}}{\|A\|_{F}^{2}}\right){\mathbb{E}}[{\bf
y}_{k}]+\frac{{\bf b}}{\|A\|_{F}^{2}}.$

This means

$\|{\mathbb{E}}[{\bf
y}_{k}]\|=\left\|\sum_{i=0}^{k-1}\left(I-\frac{AA^{\dagger}}{\|A\|_{F}^{2}}\right)^{i}\frac{{\bf
b}}{\|A\|_{F}^{2}}\right\|\leq\sum_{i=0}^{k-1}(1-\kappa_{F}^{-2})^{i}\frac{\|{\bf
b}\|}{\|A\|_{F}^{2}}\leq\frac{\kappa_{F}^{2}\|{\bf b}\|}{\|A\|_{F}^{2}}.$

Therefore,

$\displaystyle{\mathbb{E}}[\|{\bf y}_{T}\|_{\Lambda}^{2}]$ $\displaystyle\leq$
$\displaystyle T\left(\frac{2\|{\bf
b}\|^{2}}{\|A\|_{F}^{2}}+(20+\frac{\epsilon^{2}}{T})\|{\bf
x}_{*}\|^{2}+\frac{2\kappa^{2}\|{\bf b}\|Z}{\|A\|_{F}^{2}}\right)$
$\displaystyle=$ $\displaystyle T\left(\frac{2\|A{\bf
x}_{*}\|^{2}}{\|A\|_{F}^{2}}+(20+\frac{\epsilon^{2}}{T})\|{\bf
x}_{*}\|^{2}+\frac{2\kappa^{2}\|{\bf b}\|Z+2Z^{2}}{\|A\|_{F}^{2}}\right).$

Finally, by Markov’s inequality, with high probability we have

$\phi=T\frac{\|{\bf y}_{T}\|_{\Lambda}^{2}}{\|{\bf
x}_{T}\|^{2}}=O\left(T^{2}+T^{2}\frac{\kappa^{2}\|{\bf
b}\|Z+Z^{2}}{\|A\|_{F}^{2}\|{\bf x}_{*}\|^{2}}\right),$

where we used the fact that ${\bf x}_{T}\approx{\bf x}_{*}$ and $\|A{\bf
x}_{*}\|^{2}/\|A\|_{F}^{2}\leq\|{\bf x}_{*}\|^{2}$. Note that $\|A\|_{F}\|{\bf
x}_{*}\|=\|A\|_{F}\|A^{-1}{\bf b}\|\geq\|A\|_{F}\|{\bf b}\|/\|A\|\geq Z$, so
the second term is bounded by $O(T^{2}\kappa^{4}Z/\kappa_{F}^{2}\|{\bf b}\|)$.
In the worst case, $\phi=O(T^{2}\kappa^{2})$ because
$\kappa^{2}Z\leq\kappa_{F}^{2}\|{\bf b}\|$.

##  Appendix B Estimation of the convergence rate

For simplicity, we assume that $A{\bf x}_{*}={\bf b}$. Following from the
analysis of [26], the convergence rate does not change too much if the LSP is
not consistent. By the updating formula (19), we have

$\displaystyle{\bf x}_{k+1}-{\bf x}_{*}$ $\displaystyle=$ $\displaystyle{\bf
x}_{k}-{\bf
x}_{*}+\frac{1}{2}\sum_{i\in{\mathcal{T}}_{k}}(\tilde{b}_{i}-\langle\tilde{A}_{i*}|{\bf
x}_{k}\rangle)|\tilde{A}_{i*}\rangle+\frac{1}{2}\sum_{i\in{\mathcal{T}}_{k}}\langle\tilde{A}_{i*}|I-D_{i}|{\bf
x}_{k}\rangle|\tilde{A}_{i*}\rangle$ $\displaystyle=$
$\displaystyle\left(I-\frac{1}{2}\sum_{i\in{\mathcal{T}}_{k}}|\tilde{A}_{i*}\rangle\langle\tilde{A}_{i*}|\right)({\bf
x}_{k}-{\bf
x}_{*})+\frac{1}{2}\sum_{i\in{\mathcal{T}}_{k}}\langle\tilde{A}_{i*}|I-D_{i}|{\bf
x}_{k}\rangle|\tilde{A}_{i*}\rangle.$

Below, we try to bound ${\mathbb{E}}[\|{\bf x}_{k+1}-{\bf x}_{*}\|^{2}]$. We
will follow the notation of (15). But here to avoid any confusion, we denote

$\mu_{ik}:=\langle\tilde{A}_{i*}|I-D_{i}|{\bf x}_{k}\rangle.$

In the following, the result of Lemma 9 will be used, and
$|{\mathcal{T}}_{k}|=\|A\|_{F}^{2}/\|A\|^{2}$.

First, we have

$\displaystyle{\mathbb{E}}_{D}\left[\left\|\sum_{i\in{\mathcal{T}}_{k}}\langle\tilde{A}_{i*}|I-D_{i}|{\bf
x}_{k}\rangle|\tilde{A}_{i*}\rangle\right\|^{2}\right]$ $\displaystyle=$
$\displaystyle\sum_{i,j\in{\mathcal{T}}_{k}}\langle\tilde{A}_{i*}|\tilde{A}_{j*}\rangle{\mathbb{E}}_{D_{i},D_{j}}\left[\mu_{ik}\mu_{jk}\right]$
$\displaystyle=$
$\displaystyle\sum_{i\in{\mathcal{T}}_{k}}{\mathbb{E}}_{D_{i}}[\mu_{ik}^{2}]+\sum_{i\neq
j}\langle\tilde{A}_{i*}|\tilde{A}_{j*}\rangle{\mathbb{E}}_{D_{i}}\left[\mu_{ik}\right]{\mathbb{E}}_{D_{j}}[\mu_{jk}]$
$\displaystyle\leq$ $\displaystyle\frac{\|A\|_{F}^{4}\|{\bf
x}_{k}\|^{2}}{d\|A\|^{2}\min_{l\in[n]}\|A_{*l}\|^{2}}.$

By Lemma 9, we have

$\displaystyle{\mathbb{E}}_{D}\left[\left\langle\left(I-\frac{1}{2}\sum_{i\in{\mathcal{T}}_{k}}|\tilde{A}_{i*}\rangle\langle\tilde{A}_{i*}|\right)({\bf
x}_{k}-{\bf
x}_{*})\Bigg{|}\sum_{i\in{\mathcal{T}}_{k}}\mu_{ik}|\tilde{A}_{i*}\rangle\right\rangle\right]$
$\displaystyle=$
$\displaystyle\left\langle\left(I-\frac{1}{2}\sum_{i\in{\mathcal{T}}_{k}}|\tilde{A}_{i*}\rangle\langle\tilde{A}_{i*}|\right)({\bf
x}_{k}-{\bf
x}_{*})\Bigg{|}\sum_{i\in{\mathcal{T}}_{k}}{\mathbb{E}}_{D}[\mu_{ik}]|\tilde{A}_{i*}\rangle\right\rangle=0.$

Hence, after computing the mean value over $D$, we have

$\displaystyle\|{\bf x}_{k+1}-{\bf
x}_{*}\|^{2}\leq\left\|\left(I-\frac{1}{2}\sum_{i\in{\mathcal{T}}_{k}}|\tilde{A}_{i*}\rangle\langle\tilde{A}_{i*}|\right)({\bf
x}_{k}-{\bf x}_{*})\right\|^{2}+\frac{\|A\|_{F}^{4}\|{\bf
x}_{k}\|^{2}}{4d\|A\|^{2}\min_{l\in[n]}\|A_{*l}\|^{2}}.$

For the first term, we now compute its mean value over the random variable
${\mathcal{T}}_{k}$

$\displaystyle{\mathbb{E}}\left[\left\|\left(I-\frac{1}{2}\sum_{i\in{\mathcal{T}}_{k}}|\tilde{A}_{i*}\rangle\langle\tilde{A}_{i*}|\right)({\bf
x}_{k}-{\bf x}_{*})\right\|^{2}\right]$ $\displaystyle=$
$\displaystyle\langle{\bf x}_{k}-{\bf
x}_{*}|{\mathbb{E}}\left[\left(I-\frac{1}{2}\sum_{i\in{\mathcal{T}}_{k}}|\tilde{A}_{i*}\rangle\langle\tilde{A}_{i*}|\right)^{2}\right]|{\bf
x}_{k}-{\bf x}_{*}\rangle.$

And it can be shown that

$\displaystyle{\mathbb{E}}\left[\left(I-\frac{1}{2}\sum_{i\in{\mathcal{T}}_{k}}|\tilde{A}_{i*}\rangle\langle\tilde{A}_{i*}|\right)^{2}\right]$
$\displaystyle=$
$\displaystyle{\mathbb{E}}\left[I-\frac{3}{4}\sum_{i\in{\mathcal{T}}_{k}}|\tilde{A}_{i*}\rangle\langle\tilde{A}_{i*}|+\frac{1}{4}\sum_{i,j\in{\mathcal{T}}_{k},i\neq
j}|\tilde{A}_{i*}\rangle\langle\tilde{A}_{i*}|\tilde{A}_{j*}\rangle\langle\tilde{A}_{j*}|\right]$
$\displaystyle=$ $\displaystyle
I-\frac{3q}{4}\frac{A^{\dagger}A}{\|A\|_{F}^{2}}+\frac{1}{4}(q^{2}-q)\left(\frac{A^{\dagger}A}{\|A\|_{F}^{2}}\right)^{2}$
$\displaystyle\preceq$
$\displaystyle\left(1-\frac{3q}{4}\frac{\sigma_{\min}(A^{\dagger}A)}{\|A\|_{F}^{2}}+\frac{1}{4}(q^{2}-q)\left(\frac{\sigma_{\min}(A^{\dagger}A)}{\|A\|_{F}^{2}}\right)^{2}\right)I$
$\displaystyle\preceq$ $\displaystyle\left(1-\frac{1}{2\kappa^{2}}\right)I.$

In the last step, we used the result $q=\|A\|_{F}^{2}/\|A\|^{2}$ so that the
second term is $3/4\kappa^{2}$ and the third term is less than
$1/4\kappa^{4}\leq 1/4\kappa^{2}$.

Therefore,

$\displaystyle{\mathbb{E}}[\|{\bf x}_{k+1}-{\bf
x}_{*}\|^{2}]\leq(1-\frac{1}{2\kappa^{2}}){\mathbb{E}}[\|{\bf x}_{k}-{\bf
x}_{*}\|^{2}]+\frac{\|A\|_{F}^{4}{\mathbb{E}}[\|{\bf
x}_{k}\|^{2}]}{4d\|A\|^{2}\min_{l\in[n]}\|A_{*l}\|^{2}}.$

Now set $T=O(\kappa^{2}\log(2/\epsilon^{2}))$, and

$d=\frac{\|A\|_{F}^{4}T}{\epsilon^{2}\|A\|^{2}\min_{l\in[n]}\|A_{*l}\|^{2}}.$
(29)

Then

$\displaystyle{\mathbb{E}}[\|{\bf x}_{k+1}-{\bf x}_{*}\|^{2}]$
$\displaystyle\leq$ $\displaystyle(1-\frac{1}{2\kappa^{2}}){\mathbb{E}}[\|{\bf
x}_{k}-{\bf x}_{*}\|^{2}]+\frac{\epsilon^{2}}{4T}{\mathbb{E}}[\|{\bf
x}_{k}\|^{2}]$ $\displaystyle\leq$
$\displaystyle(1-\frac{1}{2\kappa^{2}}+\frac{\epsilon^{2}}{2T}){\mathbb{E}}[\|{\bf
x}_{k}-{\bf x}_{*}\|^{2}]+\frac{\epsilon^{2}}{2T}\|{\bf x}_{*}\|^{2}.$

Finally, we obtain

$\displaystyle{\mathbb{E}}[\|{\bf x}_{T}-{\bf
x}_{*}\|^{2}]\lesssim(1-\frac{1}{2\kappa^{2}}+\frac{\epsilon^{2}}{2T})^{T}\|{\bf
x}_{0}-{\bf x}_{*}\|^{2}+\frac{\epsilon^{2}}{2}\|{\bf
x}_{*}\|^{2}\leq\epsilon^{2}\|{\bf x}_{*}\|^{2}.$

##  Appendix C Estimation of $\phi$ for Kaczmarz method with averaging

The calculation here is similar to that in Appendix A. For simplicity, denote

${\bf y}_{k+1}={\bf y}_{k}+\frac{1}{2}{\bf w}_{k}+\frac{1}{2}{\bf
w}_{k}^{\prime},$

where

$\displaystyle{\bf w}_{k}:=\sum_{i\in{\mathcal{T}}_{k}}\frac{b_{i}-\langle
A_{i*}|{\bf x}_{k}\rangle}{\|A_{i*}\|^{2}}{\bf e}_{i*},\quad{\bf
w}_{k}^{\prime}:=\sum_{i\in{\mathcal{T}}_{k}}\frac{\mu_{ik}}{\|A_{i*}\|}{\bf
e}_{i*}.$

In this section, $d$ is given in formula (29). From a similar estimation in
Appendix A, we have ${\mathbb{E}}_{D}[\langle{\bf w}_{k}|{\bf
w}_{k}^{\prime}\rangle_{\Lambda}]={\mathbb{E}}_{D}[\langle{\bf y}_{k}|{\bf
w}_{k}^{\prime}\rangle_{\Lambda}]=0$ and

$\displaystyle{\mathbb{E}}_{D}[\|{\bf
w}_{k}^{\prime}\|_{\Lambda}^{2}]\leq\frac{\|A\|_{F}^{2}}{\|A\|^{2}}\frac{\|A\|_{F}^{2}\|{\bf
x}_{k}\|^{2}}{d\min_{j\in[n]}\|A_{*j}\|^{2}}=\frac{\epsilon^{2}\|{\bf
x}_{k}\|^{2}}{4T}\leq\frac{\epsilon^{2}(\|{\bf x}_{k}-{\bf x}_{*}\|^{2}+\|{\bf
x}_{*}\|^{2})}{2T}.$

As for $\|{\bf w}_{k}\|_{\Lambda}^{2}$, we still have

$\displaystyle{\mathbb{E}}[\|{\bf w}_{k}\|_{\Lambda}^{2}]$ $\displaystyle=$
$\displaystyle{\mathbb{E}}\left[\sum_{i,j\in{\mathcal{T}}_{k}}\frac{b_{i}-\langle
A_{i*}|{\bf x}_{k}\rangle}{\|A_{i*}\|^{2}}\frac{b_{j}-\langle A_{j*}|{\bf
x}_{k}\rangle}{\|A_{j*}\|^{2}}\langle{\bf e}_{i*}|{\bf
e}_{j*}\rangle\|A_{i*}\|^{2}\right]$ $\displaystyle=$
$\displaystyle\frac{\|A\|_{F}^{2}}{\|A\|^{2}}{\mathbb{E}}\left[\frac{(b_{i}-\langle
A_{i*}|{\bf x}_{k}\rangle)^{2}}{\|A_{i*}\|^{2}}\right]$ $\displaystyle=$
$\displaystyle\frac{\|{\bf b}-A{\bf x}_{k}\|^{2}}{\|A\|^{2}}$
$\displaystyle\leq$ $\displaystyle\frac{2\|{\bf b}\|^{2}}{\|A\|^{2}}+2\|{\bf
x}_{k}\|^{2}.$

By the estimation of ${\mathbb{E}}[\langle{\bf y}_{k}|{\bf
z}_{k}\rangle_{\Lambda}]$ in Appendix A and noting that
$\|\Lambda\|\leq\|A\|$, we also have

$\displaystyle{\mathbb{E}}[\langle{\bf y}_{k}|{\bf
w}_{k}\rangle_{\Lambda}]=\frac{\|A\|_{F}^{2}}{\|A\|^{2}}{\mathbb{E}}[\langle{\bf
y}_{k}|\frac{b_{i}-\langle A_{i*}|{\bf x}_{k}\rangle}{\|A_{i*}\|^{2}}{\bf
e}_{i*}\rangle]\leq\frac{7}{2}\|{\bf x}_{*}\|^{2}+\frac{5}{2}\|{\bf
x}_{k}-{\bf x}_{*}\|^{2}.$

All the estimations above do not change. The constant 1/2 in the decomposition
of ${\bf y}_{k+1}$ does not affect the upper bound, so when $Z=0$ we have

$\displaystyle{\mathbb{E}}[\|{\bf y}_{T}\|_{\Lambda}^{2}]\leq
T\left(2\frac{\|{\bf
b}\|^{2}}{\|A\|_{F}^{2}}+(20+\frac{\epsilon^{2}}{T})\|{\bf
x}_{*}\|^{2}\right),$

where $T=\widetilde{O}(\kappa^{2})$. Therefore

$\phi=T\frac{\kappa_{F}^{2}}{\kappa^{2}}\frac{\|{\bf
y}_{T}\|_{\Lambda}^{2}}{\|{\bf x}_{T}\|^{2}}\leq
T^{2}\frac{\kappa_{F}^{2}}{\kappa^{2}}\left(2\frac{\|{\bf
b}\|^{2}}{\|A\|_{F}^{2}\|{\bf
x}_{T}\|^{2}}+(20+\frac{\epsilon^{2}}{T})\frac{\|{\bf x}_{*}\|^{2}}{\|{\bf
x}_{T}\|^{2}}\right)=O(\kappa_{F}^{2}\kappa^{2}).$

When $Z\neq 0$, we similarly have $\phi=O(\kappa_{F}^{2}\kappa^{4})$.

## References

  * [1] Juan Miguel Arrazola, Alain Delgado, Bhaskar Roy Bardhan, and Seth Lloyd.  Quantum-inspired algorithms in practice.  Quantum, 4:307, 2020.  arXiv:1905.10415.
  * [2] Jacob Biamonte, Peter Wittek, Nicola Pancotti, Patrick Rebentrost, Nathan Wiebe, and Seth Lloyd.  Quantum machine learning.  Nature, 549(7671):195–202, 2017.
  * [3] Yair Censor, Stavros Andrea Zenios, et al.  Parallel optimization: Theory, algorithms, and applications.  Oxford University Press on Demand, 1997.
  * [4] Shantanav Chakraborty, András Gilyén, and Stacey Jeffery.  The Power of Block-Encoded Matrix Powers: Improved Regression Techniques via Faster Hamiltonian Simulation.  In 46th International Colloquium on Automata, Languages, and Programming (ICALP 2019), pages 33:1–33:14, 2019.  arXiv:1804.01973.
  * [5] Nadiia Chepurko, Kenneth L Clarkson, Lior Horesh, and David P Woodruff.  Quantum-inspired algorithms from randomized numerical linear algebra.  arXiv preprint arXiv:2011.04125, 2020.
  * [6] Nai-Hui Chia, András Gilyén, Tongyang Li, Han-Hsuan Lin, Ewin Tang, and Chunhao Wang.  Sampling-based sublinear low-rank matrix arithmetic framework for dequantizing quantum machine learning.  In Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing (STOC 2020), pages 387–400, 2020.  arXiv:1910.06151.
  * [7] Nai-Hui Chia, Tongyang Li, Han-Hsuan Lin, and Chunhao Wang.  Quantum-inspired classical sublinear-time algorithm for solving low-rank semidefinite programming via sampling approaches.  2019\.  arXiv:1901.03254.
  * [8] Andrew M Childs, Robin Kothari, and Rolando D Somma.  Quantum algorithm for systems of linear equations with exponentially improved dependence on precision.  SIAM Journal on Computing, 46(6):1920–1950, 2017.
  * [9] Petros Drineas, Ravi Kannan, and Michael W Mahoney.  Fast Monte Carlo algorithms for matrices II: Computing a low-rank approximation to a matrix.  SIAM Journal on computing, 36(1):158–183, 2006.
  * [10] András Gilyén, Seth Lloyd, and Ewin Tang.  Quantum-inspired low-rank stochastic regression with logarithmic dependence on the dimension.  2018\.  arXiv:1811.04909.
  * [11] András Gilyén, Zhao Song, and Ewin Tang.  An improved quantum-inspired algorithm for linear regression.  2020\.  arXiv:2009.07268.
  * [12] András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe.  Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics.  In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, pages 193–204, 2019.
  * [13] Gene H. Golub and Charles F. van Loan.  Matrix Computations.  The Johns Hopkins University Press, Baltimore, fourth edition, 2013.
  * [14] Richard Gordon, Robert Bender, and Gabor T Herman.  Algebraic reconstruction techniques (art) for three-dimensional electron microscopy and x-ray photography.  Journal of theoretical Biology, 29(3):471–481, 1970.
  * [15] Robert M Gower and Peter Richtárik.  Randomized iterative methods for linear systems.  SIAM Journal on Matrix Analysis and Applications, 36(4):1660–1690, 2015.
  * [16] Aram W Harrow, Avinatan Hassidim, and Seth Lloyd.  Quantum algorithm for linear systems of equations.  Phys Rev Lett, 103(15):150502, 2009.  arXiv:0811.3171.
  * [17] Dhawal Jethwani, Franccois Le Gall, and Sanjay K Singh.  Quantum-inspired classical algorithms for singular value transformation.  2019\.  arXiv:1910.05699.
  * [18] S Karczmarz.  Angenaherte auflosung von systemen linearer glei-chungen.  Bull. Int. Acad. Pol. Sic. Let., Cl. Sci. Math. Nat., pages 355–357, 1937.
  * [19] Iordanis Kerenidis, Jonas Landman, Alessandro Luongo, and Anupam Prakash.  q-means: A quantum algorithm for unsupervised machine learning.  In Advances in Neural Information Processing Systems, pages 4134–4144, 2019.
  * [20] Iordanis Kerenidis and Anupam Prakash.  Quantum Recommendation Systems.  In 8th Innovations in Theoretical Computer Science Conference (ITCS 2017), pages 49:1–49:21, 2017.  arXiv:1603.08675.
  * [21] Iordanis Kerenidis and Anupam Prakash.  A quantum interior point method for LPs and SDPs.  ACM Transactions on Quantum Computing, 1(1):1–32, 2020.
  * [22] Dennis Leventhal and Adrian S Lewis.  Randomized methods for linear constraints: convergence rates and conditioning.  Mathematics of Operations Research, 35(3):641–654, 2010.
  * [23] Lin Lin and Yu Tong.  Optimal polynomial based quantum eigenstate filtering with application to solving quantum linear systems.  Quantum, 4:361, 2020.
  * [24] Seth Lloyd, Masoud Mohseni, and Patrick Rebentrost.  Quantum principal component analysis.  Nature Physics, 10(9):631–633, 2014.  arXiv:1307.0401.
  * [25] Nicolas Loizou and Peter Richtárik.  Momentum and stochastic momentum for stochastic gradient, newton, proximal point and subspace descent methods.  Computational Optimization and Applications, 77(3):653–710, 2020\.
  * [26] Jacob D Moorman, Thomas K Tu, Denali Molitor, and Deanna Needell.  Randomized Kaczmarz with averaging.  2020\.  arXiv:2002.04126.
  * [27] Ion Necoara.  Faster randomized block Kaczmarz algorithms.  SIAM J Matrix Anal Appl, 40(4):1425–1452, 2019.  arXiv:1902.09946.
  * [28] Arkadi Nemirovski, Anatoli Juditsky, Guanghui Lan, and Alexander Shapiro.  Robust stochastic approximation approach to stochastic programming.  SIAM J Optim, 19(4):1574–1609, 2009.
  * [29] Yu Nesterov.  Efficiency of coordinate descent methods on huge-scale optimization problems.  SIAM J Optim, 22(2):341–362, 2012.
  * [30] Davide Orsucci and Vedran Dunjko.  On solving classes of positive-definite quantum linear systems with quadratically improved runtime in the condition number.  arXiv preprint arXiv:2101.11868, 2021.
  * [31] Patrick Rebentrost, Masoud Mohseni, and Seth Lloyd.  Quantum support vector machine for big data classification.  Phys Rev Lett, 113(13):130503, 2014.  arXiv:1307.0471.
  * [32] Peter Richtárik and Martin Takác.  Stochastic reformulations of linear systems: algorithms and convergence theory.  SIAM Journal on Matrix Analysis and Applications, 41(2):487–524, 2020.
  * [33] Changpeng Shao and Hua Xiang.  Row and column iteration methods to solve linear systems on a quantum computer.  Phys Rev A, 101(2):022322, 2020.  arXiv:1905.11686.
  * [34] Thomas Strohmer and Roman Vershynin.  A randomized kaczmarz algorithm with exponential convergence.  J Fourier Anal Appl, 15(2):262, 2009.  arXiv:math/0702226.
  * [35] Ewin Tang.  Quantum-inspired classical algorithms for principal component analysis and supervised clustering.  2018\.  arXiv:1811.00414.
  * [36] Ewin Tang.  A quantum-inspired classical algorithm for recommendation systems.  In Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing (STOC 2019), pages 217–228, 2019.  arXiv:1807.04271.