# A Catalyst Framework for the Quantum Linear System Problem via the Proximal
Point Algorithm

Junhyung Lyle Kim Nai-Hui Chia Anastasios Kyrillidis
Department of Computer Science, Rice University
{jlylekim, nc67<EMAIL_ADDRESS>
###### Abstract

Solving systems of linear equations is a fundamental problem, but it can be
computationally intensive for classical algorithms in high dimensions.
Existing quantum algorithms can achieve exponential speedups for the quantum
linear system problem (QLSP) in terms of the problem dimension, but even such
a theoretical advantage is bottlenecked by the condition number of the
coefficient matrix. In this work, we propose a new quantum algorithm for QLSP
inspired by the classical proximal point algorithm (PPA). Our proposed method
can be viewed as a meta-algorithm that allows inverting a modified matrix via
an existing QLSP_solver, thereby directly approximating the solution vector
instead of approximating the inverse of the coefficient matrix. By carefully
choosing the step size $\eta$, the proposed algorithm can effectively
precondition the linear system to mitigate the dependence on condition numbers
that hindered the applicability of previous approaches.

##  1 Introduction

Background. Solving systems of linear equations is a fundamental problem with
many applications spanning science and engineering. Mathematically, for a
given (Hermitian) matrix $A\in\mathbb{C}^{N\times N}$ and vector
$b\in\mathbb{C}^{N}$, the goal is to find the $N$-dimensional vector
$x^{\star}$ that satisfies $Ax^{\star}=b$. While classical algorithms, such as
Gaussian elimination [14, 22], conjugate gradient method [21], $LU$
factorization [37, 38], $QR$ factorization [13, 26], and iterative Krylov
subspace methods [25, 32], can solve this problem, their complexity scales (at
worst) cubically with $N$, the dimension of $A$, motivating the development of
quantum algorithms that could potentially achieve speedups.

Indeed, for the quantum linear system problem (QLSP) in Definition 1, Harrow,
Hassidim, and Lloyd (a.k.a. the HHL algorithm) [20] showed that the dependence
on the problem dimension exponentially reduces to
$\mathcal{O}\left(\text{poly}\log(N)\right)$, with query complexity of
$\mathcal{O}(\kappa^{2}/\varepsilon)$, under some (quantum) access model for
$A$ and $b$ (c.f., Definitions 2 and 3). Here, $\kappa$ is the condition
number of $A$, defined as the ratio of the largest to the smallest singular
value of $A$. Subsequent works, such as the work by Ambainis [2] and by
Childs, Kothari, and Somma [9] (a.k.a. the CKS algorithm) improve the
dependence on $\kappa$ and $\varepsilon$; see Table 1. The best quantum
algorithm for QLSP is based on the discrete adiabatic theorem [11], achieving
the query complexity of $\mathcal{O}(\kappa\cdot\log(1/\varepsilon))$,
matching the lower bound [33].111Note that this also matches the iteration
complexity of the (classical) conjugate gradient method [21].

Challenges in existing methodologies. A common limitation of the existing
quantum algorithms is that the dependence on the condition number $\kappa$
must be small to achieve the quantum advantage. To put more context, for any
quantum algorithm in Table 1 to achieve an exponential advantage over
classical algorithms, $\kappa$ needs to be in the order of
$\text{poly}\log(N)$, where $N$ is the dimension of $A$. For instance, for
$N=10^{3},$ $\kappa$ needs to be around $4$ to exhibit the exponential
advantage. However, condition numbers are often large in real-world problems
[34].

Moreover, it was proven [33, Proposition 6] that for QLSP, even when $A$ is
positive-definite, the dependence on condition number cannot be improved from
$\mathcal{O}(\kappa)$. This is in contrast to the classical algorithms, such
as the conjugate gradient method, which achieves $\mathcal{O}(\sqrt{\kappa})$
reduction in complexity for the linear systems of equations with $A\succ 0$.
Such observation reinforces the importance of alleviating the dependence on
the condition number $\kappa$ for quantum algorithms, which is our aim.

Our contributions. We present a novel meta-algorithm for solving the QLSP,
based on the proximal point algorithm (PPA) [36, 18]; see Algorithm 1.
Notably, in contrast to some existing methods that approximate $A^{-1}$ [20,
2, 9, 16, 33], our algorithm directly approximates $x^{\star}=A^{-1}b$ through
an iterative process based on PPA. Classically, PPA is known to improve the
“conditioning” of the problem at hand, compared to the gradient descent [41,
42, 1]; it can also be accelerated [19, 23].

An approximate proximal point algorithm for classical convex optimization has
been proposed under the name of “catalyst” [27] in the machine learning
community. Our proposed method operates similarly and can be viewed as a
generic acceleration scheme for QLSP where one can plug in different
QLSP_solver –e.g., HHL, CKS– to achieve generic (constant-level) acceleration.
In Figure 1, we illustrate the case where the best quantum algorithm for QLSP
[11], based on the discrete adiabatic theorem, is utilized as a subroutine for
Algorithm 1. The improvement in the query complexitycompared to the original
algorithm to achieve a fixed accuracy increases as $\kappa$ increases.

Figure 1: Query complexity scaling with respect to the condition number
$\kappa$. Here, the best quantum algorithm [11] is used as QLSP_solver for the
subroutine in Algorithm 1. Simply by “wrapping” the QLSP_solver, one can
achieve much better scaling with respect to the condition number $\kappa$.

Intuitively, by the definition of PPA detailed in Section 3, our proposed
method allows to invert a modified matrix $I+\eta A$ and arrive at the same
solution $A^{-1}b$, as shown below:

$x_{t+1}=(I+\eta A)^{-1}(x_{t}+\eta
b)\underset{\eta\rightarrow\infty}{=}A^{-1}b.$

Yet, the key feature and the main distinction here is the introduction of a
tunable (step size) parameter $\eta$ that allows pre-conditioning the linear
system. By carefully choosing $\eta$, we can invert the modified matrix
$I+\eta A$ that is better conditioned than $A$, thereby mitigating the
dependence on $\kappa$ of the existing quantum algorithms; see also Remark 3
and Figure 1.

Our meta-algorithmic framework complements and provides advantages over prior
works on QLSP. Most importantly, it alleviates the strict requirements on the
condition number, enabling quantum speedups for a broader class of problems
where $\kappa$ may be large. Moreover, the $\eta$ parameter can be tuned to
balance the convergence rate and precision demands, providing greater
flexibility in optimizing the overall algorithmic complexity.

In summary, our contributions and findings are as follows:

  * •

We propose a meta-algorithmic framework for the quantum linear system problem
(QLSP) based on the proximal point algorithm. Unlike existing quantum
algorithms for QLSP, which rely on different unitary approximations of
$A^{-1}$, our proposed method allows to invert a modified matrix with a
smaller condition number (c.f., Remark 3 and Lemma 1).

  * •

To our knowledge, this is the first framework for QLSP where a tunable
parameter $\eta$ allows the user to control the trade-off between the runtime
(query complexity) and the approximation error (solution quality).
Importantly, there exist choices of $\eta>0$ that allow to decrease the
runtime while maintaining the same error level (c.f., Theorem 3).

  * •

Our proposed method allows to achieve significant constant-level improvements
in the query complexity, even compared to the best quantum algorithm for QLSP
[11], simply by using it as a subroutine of Algorithm 1. This is possible as
the improvement “grows faster” than the overhead (c.f., Figure 2 and Theorem
6).

##  2 Problem Setup and Related Work

Notation. Matrices are represented with uppercase letters as in
$A\in\mathbb{C}^{N\times N}$; vectors are represented with lowercase letters
as in $b\in\mathbb{C}^{N}$, and are distinguished from scalars based on the
context. The condition number of a matrix $A$, denoted as $\kappa$, is the
ratio of the largest to the smallest singular value of $A$. We denote
$\|\cdot\|$ as the Euclidean $\ell_{2}$-norm. Qubit is the fundamental unit in
quantum computing, analogous to the bit in classical computing. The state of a
qubit is represented using the bra-ket notation, where a single qubit state
$|\psi\rangle\in\mathbb{C}^{2}$ can be expressed as a linear combination of
the basis states $|0\rangle=[0\;1]^{\top}\in\mathbb{C}^{2}$ and
$|1\rangle=[1\;0]^{\top}\in\mathbb{C}^{2}$, as in
$|\psi\rangle=\alpha|0\rangle+\beta|1\rangle$; here,
$\alpha,\beta\in\mathbb{C}$ are called amplitudes and encode the probability
of the qubit collapsing to either, such that $|\alpha|^{2}+|\beta|^{2}=1$.
Thus, $|\psi\rangle$ represents a normalized column vector by definition.
$|\cdot\rangle$ is a column vector (called bra), and its conjugate transpose
(called ket), denoted by $\langle\cdot|,$ is defined as
$\langle\cdot|=|\cdot\rangle^{*}.$ Generalizing the above, an $n$-qubit state
is a unit vector in $n$-qubit Hilbert space, defined as the Kronecker product
of $n$ single qubit states, i.e.,
$\mathcal{H}=\otimes_{i=1}^{n}\mathbb{C}^{2}\cong\mathbb{C}^{2^{n}}$. It is
customary to write $2^{n}=N$. Quantum states can be manipulated using quantum
gates, represented by unitary matrices that act on the state vectors. For
example, a single-qubit gate $U\in\mathbb{C}^{2\times 2}$ acting on a qubit
state $|\psi\rangle$ transforms it to $U|\psi\rangle$, altering the state’s
probability amplitudes according to the specific operation represented by $U$.

###  2.1 The Quantum Linear System Problem

In the quantum setting, the goal of the quantum linear system problem (QLSP)
is to prepare a quantum state proportional to the vector $x^{\star}$. That is,
we want to output
$|x^{\star}\rangle:=\frac{\sum_{i}x^{\star}_{i}|i\rangle}{\|\sum_{i}x^{\star}_{i}|i\rangle\|}$
where the vector $x^{\star}=[x^{\star}_{1},\dots,x^{\star}_{N}]^{\top}$
satisfies $Ax^{\star}=b$. Formally, we define the QLSP problem as below.

######  Definition 1 (Quantum Linear System Problem [9]).

Let $A$ be an $N\times N$ Hermitian matrix satisfying $\|A\|=1$ with condition
number $\kappa$ and at most $s$ nonzero entries in any row or column. Let $b$
be an $N$-dimensional vector, and let $x^{\star}:=A^{-1}b$. We define the
quantum states $|b\rangle$ and $|x^{\star}\rangle$ as in:

$\displaystyle|b\rangle:=\frac{\sum_{i=0}^{N-1}b_{i}|i\rangle}{\left\|\sum_{i=0}^{N-1}b_{i}|i\rangle\right\|}\quad\text{and}\quad|x^{\star}\rangle:=\frac{\sum_{i=0}^{N-1}x^{\star}_{i}|i\rangle}{\Big{\|}\sum_{i=0}^{N-1}x^{\star}_{i}|i\rangle\Big{\|}}.$
(1)

Given access to $A$ via $\mathcal{P}_{A}$ in Definition 3 or $\mathcal{U}_{A}$
in Definition 4, and access to the state $|b\rangle$ via $\mathcal{P}_{B}$ in
Definition 2, the goal of QLSP is to output a state $|\tilde{x}\rangle$ such
that $\||\tilde{x}\rangle-|x^{\star}\rangle\|\leqslant\varepsilon$.

As in previous works [20, 9, 2, 16, 33, 11], we assume that access to $A$ and
$b$ is provided by black-box subroutines that we detail below. We start with
the state preparation oracle for the vector $b\in\mathbb{C}^{N}$.

######  Definition 2 (State preparation oracle [20]).

Given a vector $b\in\mathbb{C}^{N}$, there exists a procedure
$\mathcal{P}_{B}$ that prepares the state
$|b\rangle:=\frac{\sum_{i}b_{i}|i\rangle}{\|\sum_{i}b_{i}|i\rangle\|}$ in time
$\mathcal{O}(\text{poly}\log(N))$.

We assume two encoding models for $A$: the sparse-matrix-access in Definition
3 denoted by $\mathcal{P}_{A}$, and the matrix-block-encoding model [15, 30]
in Definition 4, denoted by $\mathcal{U}_{A}$, respectively.

######  Definition 3 (Sparse matrix access [9]).

Given a $N\times N$ Hermitian matrix $A$ with operator norm $\|A\|\leqslant 1$
and at most $s$ nonzero entries in any row or column, $\mathcal{P}_{A}$ allows
the following mapping:

$\displaystyle|j,\ell\rangle$ $\displaystyle\mapsto|j,\nu(j,\ell)\rangle$
$\displaystyle\forall j\in[N]\text{ and }\ell\in[s],$ (2)
$\displaystyle|j,k,z\rangle|0\rangle$ $\displaystyle\mapsto|j,k,z\oplus
A_{jk}\rangle$ $\displaystyle\forall j,k\in[N],$ (3)

where $\nu:[N]\times[s]\rightarrow[N]$ in (2) computes the row index of the
$\ell$th nonzero entry of the $j$th column, and the third register of (3)
holds a bit string representing the entry of $A_{jk}$.

######  Definition 4 (Matrix block-encoding [15]).

A unitary operator $\mathcal{U}_{A}$ acting on $n+c$ qubits is called an
$(\alpha,c,\varepsilon)$-matrix-block-encoding of a $n$-qubit operator $A$ if

$\|A-\alpha((\langle 0^{c}|\otimes
I)~{}\mathcal{U}_{A}~{}(|0^{c}\rangle\otimes I))\|\leqslant\varepsilon.$

The above can also be expressed as follows:

$\mathcal{U}_{A}=\begin{bmatrix}\nicefrac{{\tilde{A}}}{{\alpha}}&*\\\
*&*\end{bmatrix}\quad\text{with}\quad\|\tilde{A}-A\|\leqslant\varepsilon,$

where $*$’s denote arbitrary matrix blocks with appropriate dimensions.

In [9], it was shown that a $(s,1,0)$-matrix-block-encoding of $A$ is possible
using a constant number of calls to $\mathcal{P}_{A}$ in Definition 3 (and
$\mathcal{O}(\text{poly}(n))$ extra elementary gates). In short,
$\mathcal{P}_{A}$ in Definition 3 implies efficient implementation of
$\mathcal{U}_{A}$ in Definition 4 (c.f., [15, Lemma 48]). We present both for
completeness as different works rely on different access models; however, our
proposed meta-algorithm can provide generic acceleration for any QLSP solver,
regardless of the encoding method.

QLSP_solver | Query Complexity | Key Technique/Result
---|---|---
HHL [20] | $\mathcal{O}\left(\kappa^{2}/\varepsilon\right)$ | First quantum algorithm for QLSP
Ambainis [2] | $\mathcal{O}\left(\kappa\log^{3}(\kappa)/\varepsilon^{3}\right)$ | Variable Time Amplitude Amplification
CKS [9] | $\mathcal{O}\left(\kappa\cdot\text{poly}\log(\kappa/\varepsilon)\right)$ | (Truncated) Chebyshev bases via LCU
Subaşı et al. [40] | $\mathcal{O}\left((\kappa\log\kappa)/\varepsilon)\right)$ | Adiabatic Randomization Method
An & Lin [3] | $\mathcal{O}\left(\kappa\cdot\text{poly}\log(\kappa/\varepsilon)\right)$ | Time-Optimal Adiabatic Method
Lin & Tong [28] | $\mathcal{O}\left(\kappa\cdot\log(\kappa/\varepsilon)\right)$ | Zeno Eigenstate Filtering
Costa, et al. [11] | $\mathcal{O}\left(\kappa\cdot\log(1/\varepsilon)\right)$ | Discrete Adiabatic Theorem
Ours |  $\kappa\rightarrow\frac{\kappa(1+\eta)}{\kappa+\eta}$ for all above | Proximal Point Algorithm
Table 1: Query complexities and key results used in related works on QLSP. Our
proposed framework in Algorithm 1 allows to improve the dependence on the
condition number $\kappa$ for any QLSP_solver.

###  2.2 Related Work

Quantum algorithms. We summarize the related quantum algorithms for QLSP and
their query complexities in Table 1; all QLSP_solvers share the exponential
improvement on the input dimension, $\mathcal{O}(\text{poly}\log(N))$. The HHL
algorithm [20] utilizes quantum subroutines including $(i)$ Hamiltonian
simulation [12, 29, 10] that applies the unitary operator $e^{iAt}$ to
$|b\rangle$ for a superposition of different times $t$, $(ii)$ phase
estimation [24] that allows to decompose $|b\rangle$ into the eigenbasis of
$A$ and to find its corresponding eigenvalues, and $(iii)$ amplitude
amplification [6, 17, 7] that allows to implement the final state with
amplitudes the same with the elements of $x^{\star}$. Subsequently, [2]
achieved a quadratic improvement on the condition number at the cost of worse
error dependence; the main technical contribution was to improve the amplitude
amplification, the previous bottleneck. CKS [9] significantly improved the
suboptimality by the linear combination of unitaries. We review these
subroutines in the appendix.

[11] is the state-of-the-art QLS algorithm based on the adiabatic framework,
which was spearheaded by [40] and improved in [3, 28]. These are significantly
different from the aforementioned HHL-based approaches. Importantly, Algorithm
1 is oblivious to such differences and provides generic acceleration.

Lower bounds. Along with the proposal of the first quantum algorithm for QLSP,
[20] also proved the lower bound of $\Omega(\kappa)$ queries to the entries of
the matrix is needed for general linear systems. In [33, Proposition 6], this
lower bound was surprisingly extended to the case of positive-definite
systems. This is in contrast to the classical optimization literature, where
methods such as the conjugate gradient method [21] achieve
$\sqrt{\kappa}$-acceleration for positive-definite systems.

##  3 The Proximal Point Algorithm for QLSP

We now introduce our proposed methodology, summarized in Algorithm 1. At a
high level, our method can be viewed as a meta-algorithm, where one can plug
in any existing QLSP_solver as a subroutine to achieve generic acceleration.
We first review the proximal point algorithm (PPA), a classical optimization
method on which Algorithm 1 is based.

###  3.1 The Proximal Point Algorithm

We take a step back from the QLSP in Definition 1 and introduce the proximal
point algorithm (PPA), a fundamental optimization method in convex
optimization [36, 18, 35, 4]. PPA is an iterative algorithm that proceeds by
minimizing the original function plus an additional quadratic term, as in:

$x_{t+1}=\arg\min_{x}\Big{\\{}f(x)+\frac{1}{2\eta}\|x-x_{t}\|_{2}^{2}\Big{\\}}.$
(4)

As a result, it changes the “conditioning” of the problem; if $f(\cdot)$ is
convex, the optimization problem in (4) can be strongly convex [1]. By the
first-order optimality condition [5, Eq. (4.22)], (4) can be written in the
following form, also known as the implicit gradient descent (IGD):

$\displaystyle x_{t+1}$ $\displaystyle=x_{t}-\eta\nabla f(x_{t+1}).$ (5)

(5) is an implicit method and generally cannot be implemented. However, the
case we are interested in is the quadratic minimization problem:

$\min_{x}f(x)=\frac{1}{2}x^{\top}Ax-b^{\top}x.$ (6)

Then, we have the closed-form update for (5) as follows:

$\displaystyle x_{t+1}$ $\displaystyle=x_{t}-\eta(Ax_{t+1}-b)=x_{t}-\eta
A(x_{t+1}-x^{\star}).$

In particular, by rearranging and unfolding, we have

$\displaystyle x_{t+1}$ $\displaystyle=(I+\eta A)^{-1}(x_{t}+\eta
b)=\dots=(I+\eta A)^{-(t+1)}x_{0}+\eta b\sum_{k=1}^{t+1}(I+\eta A)^{-k}.$ (7)

The above expression sheds some light on how applying PPA can differ from
simple inversion: $A^{-1}b$. In particular, PPA enables to invert a modified
matrix $I+\eta A$ based on $\eta$; see also Remark 3 below.

Further, since $b=Ax^{\star}$, we can equivalently express the series of
operations in (7) as follows:

$\displaystyle x_{t+1}-x^{\star}$ $\displaystyle=(I+\eta
A)^{-1}(x_{t}-x^{\star})=\cdots=(I+\eta A)^{-t}(x_{0}-x^{\star}).$ (8)

The above expression helps compute the number of iterations required for PPA,
given $\eta>0$, for finding $\varepsilon$-approximate solution, as we detail
in Section 4.

###  3.2 Meta-Algorithm for QLSP via PPA

We present our proposed method, which is extremely simple as summarized in
Algorithm 1.

Algorithm 1 Proximal Point Algorithm for the Quantum Linear System Problem

1:Input: An oracle $\mathcal{P}_{\eta,A}$ that prepares sparse-access or
block-encoding of $\frac{I+\eta A}{\|I+\eta A\|}$ (c.f., Definition 3 and
Definition 4); a state preparation oracle $\mathcal{P}_{b,x_{0},\eta}$ that
prepares $|x_{0}+\eta b\rangle$ (c.f., Definition 2); and a tunable step size
$\eta>0$.

2:Subroutine: Invoke any QLSP_solver such that $\big{|}\psi_{I+\eta
A,x_{0}+\eta b}\big{\rangle}\approx\Big{|}\big{(}\frac{I+\eta A}{\|I+\eta
A\|}\big{)}^{-1}(x_{0}+\eta b)\Big{\rangle}$.

3:Output: Normalized quantum state $\big{|}\psi_{I+\eta A,x_{0}+\eta
b}\big{\rangle}$.

4:Benefit: Improved dependence on condition number as in Remark 3.

Line 2 is the cornerstone of the algorithm where one can employ any
QLSP_solver –like HHL [20], CKS [9] or the recent work based on discrete
adiabatic approach [11]– to the (normalized) matrix $(I+\eta A)/\|I+\eta A\|$,
enabled by the PPA approach. In other words, line 3 can be seen as the output
of applying $\texttt{QLSP\\_solver}(\frac{I+\eta A}{\|I+\eta A\|},|x_{0}+\eta
b\rangle)$. We make some remarks on the input.

######  Remark 1 (Access to $(I+\eta A)/\|I+\eta A\|$).

Our algorithm necessitates an oracle $P_{\eta,A}$, which can provide either
sparse access to $(I+\eta A)/\|I+\eta A\|$ (as in Definition 3) or its block
encoding (as in Definition 4). For sparse access to $(I+\eta A)/\|I+\eta A\|$,
direct access is feasible from sparse access to $A$, when all diagonal entries
of $A$ are non-zero. Specifically, the access described in (3) is identical to
that of $A$, while the access in (2) modifies $A_{jj}$ to $(1+\eta
A_{jj})/\|I+\eta A\|$. If $A$ contains zero diagonal entries, sparse access to
$(I+\eta A)/\|I+\eta A\|$ requires only two additional uses of (3).

For block encoding, $(I+\eta A)/\|I+\eta A\|$ can be achieved through a linear
combination of block-encoded matrices, as demonstrated in [15, Lemma 52]. This
method allows straightforward adaptation of existing data structures that
facilitate sparse-access or block-encoding of $A$ to also support $(I+\eta
A)/\|I+\eta A\|$.

In addition, original data or data structures that provide sparse access or
block-encoding of $A$ can be easily modified for access to $(I+\eta
A)/\|I+\eta A\|$. For instance, suppose that we are given a matrix as a sum of
multiple small matrices as in the local Hamiltonian problem or Hamiltonian
simulation (see [31] for an introduction). Then the sparse access and block
encoding of both $A$ and $(I+\eta A)/\|I+\eta A\|$ can be efficiently derived
from the sum of small matrices.

######  Remark 2 (Access to $|x_{0}+\eta b\rangle$).

We prepare the initial state $|x_{0}+\eta b\rangle$, instead of $|b\rangle$,
reflecting the PPA update in (7). For this step, we assume there exists an
oracle $\mathcal{P}_{b,x_{0},\eta}$ such that $|x_{0}+\eta b\rangle$ can be
efficiently prepared similarly to Definition 2; otherwise, we can simply
initialize $x_{0}$ with zero vector.

Since we are inverting the (normalized) modified matrix $\frac{I+\eta
A}{\|I+\eta A\|}$, the spectrum changes as follows.

######  Lemma 1.

Let $A$ be an $N\times N$ Hermitian matrix satisfying $\|A\|=1$ with condition
number $\kappa$. Then, the condition number of the modified matrix in
Algorithm 1, $\frac{I+\eta A}{\|I+\eta A\|}$, is given by
$\hat{\kappa}=\frac{\kappa(1+\eta)}{\kappa+\eta}$.

Notice that the modified condition number $\hat{\kappa}$ depends both on
$\kappa$ and the step size parameter $\eta$ for PPA. As a result, $\eta$ plays
a crucial role in the overall performance of Algorithm 1. The introduction of
the tunable parameter $\eta$ in the context of QLSP is a main property that
differentiates Algorithm 1 from other quantum algorithms. We summarize the
trade-off of $\eta$ in the following remark.

######  Remark 3.

The modified condition number,
$\hat{\kappa}=\frac{\kappa(1+\eta)}{\kappa+\eta}$, in conjunction with the PPA
convergence in (8) introduces a trade-off based on the tunable parameter
$\eta.$

  * •

Large $\eta$ regime: setting large $\eta$ allows PPA to converge fast, as can
be seen in (15). On the other hand, the benefit of the modified condition
number diminishes and recovers the original $\kappa$:

$\displaystyle\hat{\kappa}=\frac{\kappa(1+\eta)}{\kappa+\eta}\underset{\eta\rightarrow\infty}{\longrightarrow}\kappa,$

  * •

Small $\eta$ regime: setting small $\eta$ slows down the convergence rate of
PPA and requires more number of iterations, as can be seen in (15). On the
other hand, the modified condition number becomes increasingly better
conditioned, as in:

$\displaystyle\hat{\kappa}=\frac{\kappa(1+\eta)}{\kappa+\eta}\underset{\eta\rightarrow
0}{\longrightarrow}1.$

##  4 Theoretical Analysis

As shown in Algorithm 1 as well as the PPA iteration explained in (7), the
main distinction of our proposed method from the existing QLSP algorithms is
that we invert the modified matrix $(I+\eta A)/\|I+\eta A\|$ instead of the
original matrix $A.$ Precisely, our goal is to bound the following:

$\bigg{\|}\underbrace{\big{|}\psi_{I+\eta A,x_{0}+\eta
b}\big{\rangle}-{\color[rgb]{.5,.5,.5}\big{|}\big{(}\tfrac{I+\eta A}{\|I+\eta
A\|}\big{)}^{-1}(x_{0}+\eta b)\big{\rangle}}}_{\texttt{QLSP\\_solver}\text{
error }(\text{c.f.,
Proposition}~{}\ref{prop:qlsp_error})}+\underbrace{{\color[rgb]{.5,.5,.5}\big{|}\big{(}\tfrac{I+\eta
A}{\|I+\eta A\|}\big{)}^{-1}(x_{0}+\eta
b)\big{\rangle}}-|A^{-1}b\rangle}_{\text{PPA error }(\text{c.f.,
Proposition}~{}\ref{prop:ppa-numiter})}\bigg{\|}\leqslant\varepsilon,$ (9)

where $\big{|}\psi_{I+\eta A,x_{0}+\eta b}\big{\rangle}$ is the output of
Algorithm 1, and $|A^{-1}b\rangle=\frac{A^{-1}b}{\|A^{-1}b\|}$ is the target
quantum state of QLSP based on Definition 1. In between, the term
$\big{|}\big{(}\tfrac{I+\eta A}{\|I+\eta A\|}\big{)}^{-1}(x_{0}+\eta
b)\big{\rangle}$ is added and subtracted, reflecting the modified inversion
due to PPA. Specifically, the first pair of terms quantifies the error coming
from the inexactness of any QLSP_solver used as a subroutine of Algorithm 1.
The second pair of terms quantifies the error coming from PPA in estimating
$A^{-1}b$, as can be seen in (8). In the following subsections, we analyze
each term carefully. Due to space limitations, we defer all proofs to the
appendix at the end of the paper.

Figure 2:  How the “improvement” term, the “overhead” term, and their sum
“total” allows improvement compared to the baselines from [9] and [11] (c.f.,
Theorem 6).  Right panel: using CKS [9] as baseline (c.f., Theorem 4); Left
panel: using [11] as baseline (c.f., Theorem 5). The key insight is that the
rate of “improvement” is faster than the rate of “overhead,” which plateaus
quickly thanks to the logarithm as can be seen in Theorem 6.

###  4.1 Inverting the Modified Matrix $(I+\eta A)/\|I+\eta A\|$

In the first part, we apply any QLSP_solver to invert the modified matrix
$(I+\eta A)/\|I+\eta A\|$. Importantly, we first need to encode it into a
quantum computer, and this requires normalization so that the resulting
encoded matrix is unitary, which is necessary for any quantum computer
operation (c.f., Remarks 1 and 2). As in this part we are simply invoking
existing QLSP_solver, we just have to control the error $\varepsilon_{1}$ and
its contribution to the final accuracy $\varepsilon$, as summarized below.

######  Proposition 1 (QLSP_solver error).

Assume access to the oracle $\mathcal{P}_{\eta,A}$ in Algorithm 1, similarly
to the Definition 3 or 4 (c.f., Remark 1). Further, assume access to the
oracle $\mathcal{P}_{b,x_{0},\eta}$ similarly to the Definition 2 (c.f.,
Remark 2), respectively. Then, there exists quantum algorithms, such as the
ones in Table 1, satisfying

$\bigg{\|}\Big{|}\psi_{I+\eta A,x_{0}+\eta
b}\Big{\rangle}-\Big{|}\big{(}\tfrac{I+\eta A}{\|I+\eta
A\|}\big{)}^{-1}(x_{0}+\eta b)\Big{\rangle}\bigg{\|}\leqslant\varepsilon_{1}.$
(10)

######  Example 1 (CKS polynomial [9]).

CKS performed the following polynomial approximation of $A^{-1}$:

$A^{-1}\approx P(A)=\sum_{i}\alpha_{i}\mathcal{T}_{i}(A),$ (11)

where $\mathcal{T}_{i}$ denotes the Chevyshev polynomials (of the first kind).
Then, given that $\|P(A)-A^{-1}\|\leqslant\varepsilon_{1}<\frac{1}{2}$, [9,
Proposition 9] proves the following:

$\displaystyle\left\|\frac{P(A)|\phi\rangle}{\big{\|}P(A)|\phi\rangle\big{\|}}-\frac{A^{-1}|\phi\rangle}{\big{\|}A^{-1}|\phi\rangle\big{\|}}\right\|\leqslant
4\varepsilon_{1}.$ (12)

Adapting Example 1 to Algorithm 1. The approximation in (11) is only possible
when the norm of the matrix to be inverted is upper-bounded by 1. Hence, we
can use the following approximation:

$\displaystyle\|I+\eta
A\|^{-1}\sum_{i}\alpha_{i}\mathcal{T}_{i}\left(\frac{I+\eta A}{\|I+\eta
A\|}\right)=\|I+\eta A\|^{-1}\cdot P\left(\frac{I+\eta A}{\|I+\eta
A\|}\right),$ (13)

which allows similar step as (12) to hold, and further allows the dependence
on $\kappa$ for the CKS algorithm to be alleviated as summarized in Lemma 1;
see also the illustration in Figure 2 (right).

###  4.2 Approximating $x^{\star}=A^{-1}b$ via PPA

For the PPA error (second pair of terms in (9)), the step size $\eta$ needs to
be properly set up so that a single step of PPA, i.e., $x_{1}=(I+\eta
A)^{-1}(x_{0}+\eta b)$ is close enough to $A^{-1}b=x^{\star}.$ To achieve
that, observe from (8):

$\displaystyle\|x_{t+1}-x^{\star}\|$ $\displaystyle\leqslant\|(I+\eta
A)^{-t}\|\cdot\|x_{0}-x^{\star}\|\leqslant\frac{1}{(1+\eta\sigma_{\min})^{t}}\|x_{0}-x^{\star}\|,$
(14)

where $\sigma_{\min}$ is the smallest singular value of $A$. We desire the RHS
to be less than $\varepsilon_{2}$. Denoting $\|x_{0}-x^{\star}\|:=d$ and using
$\sigma_{\min}=1/\kappa$ (c.f., Definition 1), we can compute the lower bound
on the number of iterations $t$ as:

$\displaystyle\frac{d}{(1+\eta/\kappa)^{t}}$
$\displaystyle\leqslant\varepsilon_{2}\implies\frac{\log(\frac{d}{\varepsilon_{2}})}{\log(1+\frac{\eta}{\kappa})}\leqslant
t.$ (15)

Based on the above analysis, we can compute the number of iterations $t$
required to have $\varepsilon$-optimal solution of the quadratic problem in
(6). Here, (15) is well defined in the sense that the lower bound on $t$ is
positive, as long as $\eta>0$. In other words, if the LHS of (15) is less than
1, that means PPA converges to $\varepsilon_{2}$-approximate solution in one
step, with a proper step size.

Again, our goal is to achieve (9). Hence, we have to characterize how (14)
results in the proximity in corresponding normalized quantum states. We
utilize Lemma 2 below.

######  Lemma 2.

Let $x$ and $y$ be two vectors. Suppose $\|x-y\|\leqslant\epsilon$ for some
small positive scalar $\epsilon$. Then, the distance between the normalized
vectors satisfies the following:

$\left\|\frac{x}{\|x\|}-\frac{y}{\|y\|}\right\|\leqslant\frac{\epsilon}{\sqrt{\|x\|\cdot\|y\|}}.$

Now, we characterize the error of the normalized quantum state of as an output
of PPA.

######  Proposition 2.

Running the PPA in (7) for a single iteration with
$\eta=\kappa\left(\frac{d}{\varepsilon_{2}}-1\right)$, where
$d:=\|x_{0}-A^{-1}b\|$, results in the normalized quantum state satisfying:

$\displaystyle\left\|\big{|}\big{(}\tfrac{I+\eta A}{\|I+\eta
A\|}\big{)}^{-1}(x_{0}+\eta
b)\big{\rangle}-\big{|}A^{-1}b\big{\rangle}\right\|\leqslant\frac{\varepsilon_{2}}{\Psi},$
(16)

where $\Psi:=\sqrt{\|(I+\eta A)^{-1}(x_{0}+\eta b)\|\cdot\|A^{-1}b\|}$.

###  4.3 Overall Complexity and Improvement

Equipped with Propsitions 1 and 2, we arrive at the following theorem.

######  Theorem 3 (Main result).

Consider solving the QLSP problem in Definition 1 with Algorithm 1, of which
the approximation error $\varepsilon$ can be decomposed as (9), recalled
below:

$\bigg{\|}\underbrace{\big{|}\psi_{I+\eta A,x_{0}+\eta
b}\big{\rangle}-{\color[rgb]{.5,.5,.5}\big{|}\big{(}\tfrac{I+\eta A}{\|I+\eta
A\|}\big{)}^{-1}(x_{0}+\eta b)\big{\rangle}}}_{\texttt{QLSP\\_solver}\text{
error}~{}\leqslant~{}\varepsilon_{1}~{}(\text{c.f.,
Proposition}~{}\ref{prop:qlsp_error})}+\underbrace{{\color[rgb]{.5,.5,.5}\big{|}\big{(}\tfrac{I+\eta
A}{\|I+\eta A\|}\big{)}^{-1}(x_{0}+\eta
b)\big{\rangle}}-|A^{-1}b\rangle}_{\text{PPA
error}~{}\leqslant~{}\varepsilon_{2}/\Psi~{}(\text{c.f.,
Proposition}~{}\ref{prop:ppa-numiter})}\bigg{\|}\leqslant\varepsilon.$

Suppose the existence of a QLSP_solver satisfying the assumptions of
Proposition 1 such that (10) is satisfied with
$\varepsilon_{1}=\varepsilon/c$, for $c>1$. Further, suppose the assumptions
of Proposition 2 hold, i.e., single-run PPA with
$\eta=\kappa\left(\frac{d}{\varepsilon_{2}}-1\right)$ is implemented with
accuracy $\varepsilon_{2}=\left(1-\frac{1}{c}\right)\varepsilon\cdot\Psi$.
Then, the output of Algorithm 1 satisfies:

$\displaystyle\bigg{\|}\big{|}\psi_{I+\eta A,x_{0}+\eta
b}\big{\rangle}-\big{|}A^{-1}b\big{\rangle}\bigg{\|}\leqslant\varepsilon.$

Moreover, the dependence on the condition number of QLSP_solver changes from
$\kappa$ to $\frac{\kappa(1+\eta)}{\kappa+\eta}$.

We now further interpret the analysis from the previous subsection and compare
it with other QLSP algorithms (e.g., CKS [9]). We first recall the CKS query
complexity below.

######  Theorem 4 (CKS complexity, Theorem 5 in [9]).

The QLSP in Definition 1 can be solved to $\varepsilon$-accuracy by a quantum
algorithm that makes
$\mathcal{O}\left(\kappa\cdot\emph{poly}\log\left(\frac{\kappa}{\varepsilon}\right)\right)$
queries to the oracles $\mathcal{P}_{A}$ in Definition 3 and $\mathcal{P}_{B}$
in Definition 2.

It was an open question whether there exists a quantum algorithm that can
match the lower bound:
$\Omega\left(\kappa\cdot\log\left(\frac{1}{\varepsilon}\right)\right)$ [33]. A
recent quantum algorithm based on the complicated discrete adiabatic theorem
[11] was shown to match this lower bound. We recall their result below.

######  Theorem 5 (Optimal complexity, Theorem 19 in [11]).

The QLSP in Definition 1 can be solved to $\varepsilon$-accuracy by a quantum
algorithm that makes
$\mathcal{O}\left(\kappa\cdot\log\left(\frac{1}{\varepsilon}\right)\right)$
queries to the oracles $\mathcal{U}_{A}$ in Definition 4 and $\mathcal{P}_{B}$
in Definition 2.

A natural direction to utilize Algorithm 1 is to use the best QLSP_solver
[11], which has the query complexity
$\mathcal{O}\left(\kappa\cdot\log\left(\frac{1}{\varepsilon}\right)\right)$
from Theorem 5. We summarize this in the next theorem.

######  Theorem 6 (Improving the optimal complexity).

Consider running Algorithm 1 with [11] as the candidate for QLSP_solver, which
has the original complexity of
$\mathcal{O}\left(\kappa\cdot\log\left(\frac{1}{\varepsilon}\right)\right)$
(c.f., Theorem 5). The modified complexity of Algorithm 1 via Theorem 3 can be
written and decomposed to:

$\displaystyle\hat{\kappa}\cdot\log\left(\frac{c}{\varepsilon}\right)=\underbrace{\frac{\kappa(1+\eta)}{\kappa+\eta}\cdot\log\left(\frac{1}{\varepsilon}\right)}_{\text{Improvement}}+\underbrace{\frac{\kappa(1+\eta)}{\kappa+\eta}\cdot\log(c)}_{\text{Overhead}},$
(17)

where the “improvement” comes from $\hat{\kappa}\leqslant\kappa$, and the
“overhead” is due to the weight of $\varepsilon_{1}=\varepsilon/c$ (and
subsequently $\varepsilon_{2}=\left(1-\frac{1}{c}\right)\varepsilon\cdot\Psi$)
in Theorem 3. Further, with
$\eta=\kappa\left(\frac{d}{\varepsilon_{2}}-1\right)$, it follows

$\displaystyle\hat{\kappa}=\frac{\kappa(1+\eta)}{\kappa+\eta}=\kappa-\frac{(c-1)(\kappa-1)\Psi\varepsilon}{c\cdot
d}\leqslant\kappa.$ (18)

We illustrate Theorem 6 in Figure 2 (left); thanks to the flexibility of
Algorithm 1 and Theorem 3, similar analysis can done with CKS [9] as the
baseline, as illustrated in Figure 2 (right).

Intuitively, based on the decomposition in (17), we can see that the constant
$c$ –which controls the weight of $\varepsilon_{1}$ and $\varepsilon_{2}$ in
Theorem 3– enters a logarithmic term in the “overhead.” On the contrary, the
(additive) improvement in $\kappa$ is proportional to the term
$\frac{(c-1)}{c}$, as can be seen in (18).

##  5 Conclusion and Discussion

In this work, we proposed a novel quantum algorithm for solving the quantum
linear systems problem (QLSP), based on the proximal poin algorithm (PPA).
Specifically, we showed that implementing a single-step PPA is possible by
utilizing existing QLSP solvers. We designed a meta-algorithm where any QLSP
solver can be utilized as a subroutine to improve the dependence on the
condition number. Even the current best quantum algorithm [11] can be
significantly accelerated via Algorithm 1, especially when the problem is ill-
conditioned.

Limitations. A main limitation of this work is that only constant-level
improvement is attained; yet, due to the existing lower bound [20, 33],
asymptotic improvement is not possible. Another limitation is that only a
single-step PPA is implemented in Algorithm 1 at the moment. This necessitates
$\eta$ to be fairly large, as can be seen in (7) or (15). Then, in conjunction
with Remark 3 and (18), the benefit from the modified condition number
$\hat{\kappa}=\frac{\kappa(1+\eta)}{\kappa+\eta}$ diminishes. To address this,
two interesting future directions can be considered.

Implementing multi-step PPA. A natural idea is to implement multi-step PPA.
Based on (7), two-step PPA can be written as :

$\displaystyle x_{t+1}$ $\displaystyle=(I+\eta A)^{-2}(x_{t-1}+\eta b)+(I+\eta
A)^{-1}\eta b$ $\displaystyle=\big{(}(I+\eta A)^{-2}+(I+\eta
A)^{-1}\big{)}\eta b,$

where in the second step we set $x_{t-1}$ with zero vector (i.e.,
initialization). As can be seen, this requires implementing different powers
(and their addition) of the modified matrix $I+\eta A$, which is possible via
[15, Lemmas 52 and 53], and may provide overall improvement over Algorithm 1.

Warm starting. The modified condition number $\hat{\kappa}$ of $(I+\eta
A)/\|I+\eta A\|$ relies on the initial point of PPA, $x_{0}$. In particular,
$\hat{\kappa}$ can be expressed as in (18) from Theorem 6, where
$d:=\|x_{0}-x^{\star}\|$. Therefore, a better initialization $x_{0}$ can
result in a bigger improvement in the overall complexity.

To achieve that, for instance, suppose that we have some knowledge about the
solution $x^{\star}$, such that it is not far from a specific point or a set
of vectors. Then we can use that set of vectors as our starting point.
Alternatively, we can also run other QLSP solvers or quantum-inspired linear
system solvers (e.g., [8, 39]) with constant accuracy, and use the solution as
the initialization.

## References

  * Ahn and Sra [2022] Ahn, K. and S. Sra (2022).  Understanding nesterov’s acceleration via proximal point method.  In Symposium on Simplicity in Algorithms (SOSA), pp. 117–130. SIAM.
  * Ambainis [2012] Ambainis, A. (2012).  Variable time amplitude amplification and quantum algorithms for linear algebra problems.  In STACS’12 (29th Symposium on Theoretical Aspects of Computer Science), Volume 14, pp. 636–647. LIPIcs.
  * An and Lin [2022] An, D. and L. Lin (2022).  Quantum linear system solver based on time-optimal adiabatic quantum computing and quantum approximate optimization algorithm.  ACM Transactions on Quantum Computing 3(2), 1–28.
  * Bauschke and Combettes [2019] Bauschke, H. and P. Combettes (2019).  Convex analysis and monotone operator theory in hilbert spaces, corrected printing.
  * Boyd and Vandenberghe [2004] Boyd, S. P. and L. Vandenberghe (2004).  Convex optimization.  Cambridge university press.
  * Brassard and Hoyer [1997] Brassard, G. and P. Hoyer (1997).  An exact quantum polynomial-time algorithm for Simon’s problem.  In Proceedings of the Fifth Israeli Symposium on Theory of Computing and Systems, pp. 12–23. IEEE.
  * Brassard et al. [2002] Brassard, G., P. Hoyer, M. Mosca, and A. Tapp (2002).  Quantum amplitude amplification and estimation.  Contemporary Mathematics 305, 53–74.
  * Chia et al. [2022] Chia, N.-H., A. P. Gilyén, T. Li, H.-H. Lin, E. Tang, and C. Wang (2022).  Sampling-based sublinear low-rank matrix arithmetic framework for dequantizing quantum machine learning.  Journal of the ACM 69(5), 1–72.
  * Childs et al. [2017] Childs, A. M., R. Kothari, and R. D. Somma (2017).  Quantum algorithm for systems of linear equations with exponentially improved dependence on precision.  SIAM Journal on Computing 46(6), 1920–1950.
  * Childs et al. [2018] Childs, A. M., D. Maslov, Y. Nam, N. J. Ross, and Y. Su (2018).  Toward the first quantum simulation with quantum speedup.  Proceedings of the National Academy of Sciences 115(38), 9456–9461.
  * Costa et al. [2022] Costa, P. C., D. An, Y. R. Sanders, Y. Su, R. Babbush, and D. W. Berry (2022).  Optimal scaling quantum linear-systems solver via discrete adiabatic theorem.  PRX Quantum 3(4), 040303.
  * Feynman [1982] Feynman, R. P. (1982).  Simulating physics with computers.  International Journal of Theoretical Physics 21(6/7).
  * Francis [1961] Francis, J. G. (1961).  The QR transformation a unitary analogue to the LR transformation—part 1.  The Computer Journal 4(3), 265–271.
  * Gauss [1877] Gauss, C. F. (1877).  Theoria motus corporum coelestium in sectionibus conicis solem ambientium, Volume 7.  FA Perthes.
  * Gilyén et al. [2019] Gilyén, A., Y. Su, G. H. Low, and N. Wiebe (2019).  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, pp. 193–204.
  * Gribling et al. [2021] Gribling, S., I. Kerenidis, and D. Szilágyi (2021).  Improving quantum linear system solvers via a gradient descent perspective.  arXiv preprint arXiv:2109.04248.
  * Grover [1998] Grover, L. K. (1998).  Quantum computers can search rapidly by using almost any transformation.  Physical Review Letters 80(19), 4329.
  * Güler [1991] Güler, O. (1991).  On the convergence of the proximal point algorithm for convex minimization.  SIAM journal on control and optimization 29(2), 403–419.
  * Güler [1992] Güler, O. (1992).  New proximal point algorithms for convex minimization.  SIAM Journal on Optimization 2(4), 649–664.
  * Harrow et al. [2009] Harrow, A. W., A. Hassidim, and S. Lloyd (2009).  Quantum algorithm for linear systems of equations.  Physical review letters 103(15), 150502.
  * Hestenes and Stiefel [1952] Hestenes, M. R. and E. Stiefel (1952).  Methods of conjugate gradients for solving.  Journal of research of the National Bureau of Standards 49(6), 409.
  * Higham [2011] Higham, N. J. (2011).  Gaussian elimination.  Wiley Interdisciplinary Reviews: Computational Statistics 3(3), 230–238.
  * Kim et al. [2022] Kim, J. L., P. Toulis, and A. Kyrillidis (2022).  Convergence and stability of the stochastic proximal point algorithm with momentum.  In Learning for Dynamics and Control Conference, pp. 1034–1047. PMLR.
  * Kitaev [1995] Kitaev, A. Y. (1995).  Quantum measurements and the abelian stabilizer problem.  arXiv preprint quant-ph/9511026.
  * Krylov [1931] Krylov, A. (1931).  On the numerical solution of equation by which are determined in technical problems the frequencies of small vibrations of material systems.  News Acad. Sci. USSR 7, 491–539.
  * Kublanovskaya [1962] Kublanovskaya, V. N. (1962).  On some algorithms for the solution of the complete eigenvalue problem.  USSR Computational Mathematics and Mathematical Physics 1(3), 637–657.
  * Lin et al. [2015] Lin, H., J. Mairal, and Z. Harchaoui (2015).  A universal catalyst for first-order optimization.  Advances in neural information processing systems 28.
  * Lin and Tong [2020] Lin, L. and Y. Tong (2020).  Optimal polynomial based quantum eigenstate filtering with application to solving quantum linear systems.  Quantum 4, 361.
  * Lloyd [1996] Lloyd, S. (1996).  Universal quantum simulators.  Science 273(5278), 1073–1078.
  * Low and Chuang [2019] Low, G. H. and I. L. Chuang (2019).  Hamiltonian simulation by qubitization.  Quantum 3, 163.
  * Nielsen and Chuang [2001] Nielsen, M. A. and I. L. Chuang (2001).  Quantum computation and quantum information, Volume 2.  Cambridge university press Cambridge.
  * Nocedal and Wright [1999] Nocedal, J. and S. J. Wright (1999).  Numerical optimization.  Springer.
  * Orsucci and Dunjko [2021] Orsucci, D. and V. Dunjko (2021).  On solving classes of positive-definite quantum linear systems with quadratically improved runtime in the condition number.  Quantum 5, 573.
  * Papyan [2020] Papyan, V. (2020).  Traces of class/cross-class structure pervade deep learning spectra.  The Journal of Machine Learning Research 21(1), 10197–10260.
  * Parikh et al. [2014] Parikh, N., S. Boyd, et al. (2014).  Proximal algorithms.  Foundations and trends® in Optimization 1(3), 127–239.
  * Rockafellar [1976] Rockafellar, R. T. (1976).  Monotone operators and the proximal point algorithm.  SIAM journal on control and optimization 14(5), 877–898.
  * Schwarzenberg-Czerny [1995] Schwarzenberg-Czerny, A. (1995).  On matrix factorization and efficient least squares solution.  Astronomy and Astrophysics Supplement, v. 110, p. 405 110, 405.
  * Shabat et al. [2018] Shabat, G., Y. Shmueli, Y. Aizenbud, and A. Averbuch (2018).  Randomized LU decomposition.  Applied and Computational Harmonic Analysis 44(2), 246–272.
  * Shao and Montanaro [2022] Shao, C. and A. Montanaro (2022).  Faster quantum-inspired algorithms for solving linear systems.  ACM Transactions on Quantum Computing 3(4), 1–23.
  * Subaşı et al. [2019] Subaşı, Y., R. D. Somma, and D. Orsucci (2019).  Quantum algorithms for systems of linear equations inspired by adiabatic quantum computing.  Physical review letters 122(6), 060504.
  * Toulis et al. [2014] Toulis, P., E. Airoldi, and J. Rennie (2014).  Statistical analysis of stochastic gradient methods for generalized linear models.  In International Conference on Machine Learning, pp. 667–675. PMLR.
  * Toulis and Airoldi [2017] Toulis, P. and E. M. Airoldi (2017).  Asymptotic and finite-sample properties of estimators based on stochastic gradients.  The Annals of Statistics 45(4), 1694–1727.

## Appendix

##  Appendix A Missing Proofs

###  A.1 Proof of Lemma 1

###### Proof.

By the assumption on QLSP in Definition 1, the singular values of $A$ is
contained in the interval $\big{[}\frac{1}{\kappa},1\big{]}$. Thus, for the
normalized modified matrix, i.e., $\frac{I+\eta A}{\|I+\eta A\|}$, the
singular values are contained in the interval
$\big{[}\frac{\kappa+\eta}{\kappa(1+\eta)},1\big{]}$, leading to the modified
condition number $\hat{\kappa}=\frac{\kappa(1+\eta)}{\kappa+\eta}$. ∎

###  A.2 Proof of Lemma 2

###### Proof.

We need to bound:

$\displaystyle\left\|\hat{x}-\hat{y}\right\|=\left\|\frac{x}{\|x\|}-\frac{y}{\|y\|}\right\|.$

Consider the inner product form:

$\displaystyle\left\|\frac{x}{\|x\|}-\frac{y}{\|y\|}\right\|^{2}=2-2\left\langle\frac{x}{\|x\|},\frac{y}{\|y\|}\right\rangle=2-2\frac{\langle
x,y\rangle}{\|x\|\|y\|}.$

The inner product can be bounded using $\|x-y\|\leqslant\epsilon$:

$\displaystyle\langle x,y\rangle$
$\displaystyle=\frac{1}{2}\left(\|x\|^{2}+\|y\|^{2}-\|x-y\|^{2}\right)$
$\displaystyle\geqslant\frac{1}{2}\left(\|x\|^{2}+\|y\|^{2}-\epsilon^{2}\right).$

Now, substituting this into the normalized inner product:

$\displaystyle\frac{\langle
x,y\rangle}{\|x\|\|y\|}\geqslant\frac{\frac{1}{2}\left(\|x\|^{2}+\|y\|^{2}-\epsilon^{2}\right)}{\|x\|\|y\|}=\frac{1}{2}\left(\frac{\|x\|}{\|y\|}+\frac{\|y\|}{\|x\|}-\frac{\epsilon^{2}}{\|x\|\|y\|}\right).$

Using $\frac{\|x\|}{\|y\|}+\frac{\|y\|}{\|x\|}\geqslant 2$ from AM-GM
inequality, we get:

$\displaystyle\frac{\langle x,y\rangle}{\|x\|\|y\|}\geqslant
1-\frac{\epsilon^{2}}{2\|x\|\|y\|}.$

Substituting the above back into the distance expression, we get:

$\displaystyle\left\|\frac{x}{\|x\|}-\frac{y}{\|y\|}\right\|^{2}\leqslant
2-2\left(1-\frac{\epsilon^{2}}{2\|x\|\|y\|}\right)=\frac{\epsilon^{2}}{\|x\|\|y\|},$

which implies the desired result:

$\displaystyle\left\|\frac{x}{\|x\|}-\frac{y}{\|y\|}\right\|\leqslant\frac{\epsilon}{\sqrt{\|x\|\|y\|}}.$

∎

###  A.3 Proof of Proposition 1

###### Proof.

Since the input matrix $\frac{I+\eta A}{\|I+\eta A\|}$ is normalized, and due
to the reasoning in Remark 1, we can have sparse-acess or block-encoding of
$\frac{I+\eta A}{\|I+\eta A\|}$. Then, the result of Proposition 1 simply
follows from the result of each related work in Table 1. ∎

###  A.4 Proof of Proposition 2

###### Proof.

As explained in Section 4.2 of the main text, we can compute the number of
iterations for PPA using (15), recalled below, which is based on (14).

$\displaystyle\frac{d}{(1+\eta/\kappa)^{t}}$
$\displaystyle\leqslant\varepsilon_{2}\implies\frac{\log(\frac{d}{\varepsilon_{2}})}{\log(1+\frac{\eta}{\kappa})}\leqslant
t.$

In this work, since we are implementing a single-step PPA, we want $t=1$ from
(15). To achieve that, we can have:

$\displaystyle\log\left(\frac{d}{\varepsilon_{2}}\right)$
$\displaystyle=\log\left(1+\frac{\eta}{\kappa}\right)\Rightarrow$
$\displaystyle\eta$
$\displaystyle=\kappa\left(\frac{d}{\varepsilon_{2}}-1\right).$

∎

###  A.5 Proof of Theorem 3

###### Proof.

Recall the decomposition:

$\bigg{\|}\underbrace{\big{|}\psi_{I+\eta A,x_{0}+\eta
b}\big{\rangle}-{\color[rgb]{.5,.5,.5}\big{|}\big{(}\tfrac{I+\eta A}{\|I+\eta
A\|}\big{)}^{-1}(x_{0}+\eta b)\big{\rangle}}}_{\texttt{QLSP\\_solver}\text{
error}~{}\leqslant~{}\varepsilon_{1}}+\underbrace{{\color[rgb]{.5,.5,.5}\big{|}\big{(}\tfrac{I+\eta
A}{\|I+\eta A\|}\big{)}^{-1}(x_{0}+\eta
b)\big{\rangle}}-|A^{-1}b\rangle}_{\text{PPA
error}~{}\leqslant~{}\varepsilon_{2}}\bigg{\|}\leqslant\varepsilon.$

The QLSP_solver error can be provided by Proposition 1, with the accuracy of
$\varepsilon_{1}=\varepsilon/c$, for $c>1$. Similarly, $\eta$ is chosen based
on Proposition 2, with
$\varepsilon_{2}=\left(1-\frac{1}{c}\right)\varepsilon\cdot\Psi$. Then, by
adding the result from each proposition, we have

$\displaystyle\varepsilon_{1}+\frac{\varepsilon_{2}}{\Psi}=\frac{\varepsilon}{c}+\left(1-\frac{1}{c}\right)\varepsilon=\varepsilon.$

Note that since QLSP_solver is only called once in Algorithm 1, the final
accuracy $\varepsilon$ is the same for the output of Algorithm 1 and the
subroutine QLSP_solver being called. ∎

###  A.6 Proof of Theorem 6

(17) simply follows from plugging in the modified condition number
$\hat{\kappa}$, along with $\varepsilon_{1}=\varepsilon/c$. That is,

$\displaystyle\mathcal{O}\left(\hat{\kappa}\cdot\log\left(\frac{1}{\varepsilon_{1}}\right)\right)$
$\displaystyle=\mathcal{O}\left(\frac{\kappa(1+\eta)}{\kappa+\eta}\cdot\log\left(\frac{c}{\varepsilon}\right)\right)$
$\displaystyle=\mathcal{O}\left(\frac{\kappa(1+\eta)}{\kappa+\eta}\cdot\log\left(\frac{1}{\varepsilon}\right)+\frac{\kappa(1+\eta)}{\kappa+\eta}\cdot\log(c)\right).$

Note that similar analysis can be done for CKS [9], starting from
$\mathcal{O}\left(\hat{\kappa}\cdot\log\left(\frac{\hat{\kappa}}{\varepsilon_{1}}\right)\right)$.

For (18), it follows from the below steps:

$\displaystyle\frac{\kappa(1+\eta)}{\kappa+\eta}$
$\displaystyle=\frac{1+\eta}{1+\eta/\kappa}$
$\displaystyle=\frac{1+\frac{\kappa
dc}{(c-1)\varepsilon\Psi}-\kappa}{\frac{\kappa dc}{(c-1)\varepsilon\Psi}}$
$\displaystyle=\kappa-\frac{\kappa-1}{\frac{dc}{(c-1)\varepsilon\Psi}}$
$\displaystyle=\kappa-\frac{(c-1)\varepsilon\Psi(\kappa-1)}{dc}.$