# On solving classes of positive-definite quantum linear systems with
quadratically improved runtime in the condition number

Davide Orsucci  Institut für Kommunikation und Navigation, Deutsches Zentrum
für Luft- und Raumfahrt (DLR), Münchener Str. 20, 82234 Weßling, Germany
Vedran Dunjko  Leiden University, Niels Bohrweg 1, 2333 CA Leiden, The
Netherlands

(November 1, 2021)

###### Abstract

Quantum algorithms for solving the Quantum Linear System (QLS) problem are
among the most investigated quantum algorithms of recent times, with potential
applications including the solution of computationally intractable
differential equations and speed-ups in machine learning. A fundamental
parameter governing the efficiency of QLS solvers is $\kappa$, the condition
number of the coefficient matrix $A$, as it has been known since the inception
of the QLS problem that for worst-case instances the runtime scales at least
linearly in $\kappa$ [1]. However, for the case of positive-definite matrices
classical algorithms can solve linear systems with a runtime scaling as
$\sqrt{\kappa}$, a quadratic improvement compared to the the indefinite case.
It is then natural to ask whether QLS solvers may hold an analogous
improvement. In this work we answer the question in the negative, showing that
solving a QLS entails a runtime linear in $\kappa$ also when $A$ is positive
definite. We then identify broad classes of positive-definite QLS where this
lower bound can be circumvented and present two new quantum algorithms
featuring a quadratic speed-up in $\kappa$: the first is based on efficiently
implementing a matrix-block-encoding of $A^{-1}$, the second constructs a
decomposition of the form $A=LL^{\dagger}$ to precondition the system. These
methods are widely applicable and both allow to efficiently solve BQP-complete
problems.

##  1 Introduction

Quantum computation is described using the formalism of linear algebra,
suggesting that quantum methods may be intrinsically well-suited to perform
linear algebraic tasks. Algorithms solving linear systems of equations, in
particular, are at the cornerstone of linear algebra [2, Chapter 2], having
many direct applications and playing a pivotal role in several computational
methods [3]. In the seminal work of Harrow, Hassadim, and Lloyd (HHL) the so-
called Quantum Linear System (QLS) problem was introduced and a quantum
algorithm was presented that allows solving the QLS exponentially faster than
classical algorithms solving classical linear systems [1]. In subsequent
works, several new algorithms have been put forward that solve QLS with
further increased efficiency in comparison to the original HHL algorithm,
improving the runtime dependence on the condition number [4], on the precision
[5] and on the sparsity [6]. Recently, a new approach inspired by adiabatic
quantum computation has introduced a significantly simpler quasi-optimal
solving algorithm [7, 8], significantly narrowing the gap with experimental
implementations [9], making the algorithm compatible with Near-term
Intermediate Scale Quantum (NISQ) devices [10, 11] and leading to the
development of the presently most efficient QLS solvers [12].

A key idea underpinning the possibility of achieving large quantum speed-ups
in linear algebra tasks is the fact that an exponentially large complex vector
can be compactly encoded in the amplitudes of a pure quantum state; e.g., a
$n$-qubit state is described via $2^{n}$ amplitudes. This intuition is indeed
correct for the QLS problem, which has been proven to be BQP-complete [1]: any
quantum computation can be re-formulated as a QLS with only a polynomial
overhead and therefore there exist families of QLS problems that afford super-
polynomial speed-ups compared to classical solution methods (unless
$\textsf{BPP}=\textsf{BQP}$111BPP is the class of decision problems that can
be solved in bounded-error probabilistic polynomial time, BQP are those that
can be solved with bounded-error polynomial time quantum computations. Loosely
speaking, $\textsf{BPP}=\textsf{BQP}$ would mean that the power of quantum
computers is equal to that of classical computers.). While this reduction
shows that almost certainly there exist families of QLS problems that allow an
exponential speed-up compared to all classical methods, the crucial question
is whether there are _natural_ problems that can be directly formulated and
solved as QLS. The ubiquity of linear systems seems to suggests their quantum
variant should be broadly applicable as well, but still it is not guaranteed,
since in the QLS setting further significant constraints have to be met to
obtain exponential speed-ups [13].

A prominent fundamental bottleneck of QLS solvers is that, to efficiently
obtain the solution, it is not sufficient that the coefficient matrix $A$ of
the system $A\textbf{x}=\textbf{b}$ is invertible, but it also have to be
_robustly_ invertible, that is _well-conditioned_ : it is required that the
_condition number_ of the matrix $A$, defined as the ratio between the largest
and the smallest singular values, is small. In fact, solving a QLS necessarily
entails a runtime scaling at least linearly in the condition number (unless
$\textsf{BQP}=\textsf{PSPACE}$222PSPACE is the class of decision problems that
can be solved in classical polynomial space. With a Feymann sum-over-paths
approach one can show that any quantum computation can be classically
simulated in exponential time but with only polynomial space, thus
$\textsf{BQP}\subseteq\textsf{PSPACE}$. It is widely believed that
$\textsf{BQP}\neq\textsf{PSPACE}$.) as was already proven in Ref. [1].
Therefore, an exponential quantum speed-up for QLS solving is achievable only
when the condition number scales polylogarithmically with the system size and,
unfortunately, it seems rather difficult to find natural examples of matrix
families that exhibit such mild growth of the condition number [14]. However,
polynomial quantum speed-ups for linear system solving should be rather
broadly achievable and could still provide a quantum advantage if the degree
of the polynomial is large enough [15]. In this view, obtaining a further
quadratic improvement in the dependence on the condition number for some
restricted classes of matrices, that are however of wide practical interest,
could be of the uttermost importance for obtaining a broader impact of QLS
solving algorithms. A previous publication showed that in the context of
quantum algorithms for solving certain Markov chain problems the use of
specialised QLS solvers for positive-definite matrices provides an improvement
in the condition number dependence [16]. Exploring the general positive-
definite QLS problem is the main focus of this work.

In the rest of the Introduction we motivate why quadratic speed-ups in the
condition number in positive-definite QLS may be expected (Section 1.1),
review some related results present in the literature (Section 1.2) and then
proceed to give high-level overviews of our two algorithms (Section 1.3 and
Section 1.4). In Section 2 we fix the notation and give the main definitions.
In Section 3 we prove that QLS solvers require, even when restricting to
positive-definite matrices, a runtime scaling linearly in the condition
number. We then move to the main results of this work, that is, achieving
improved runtime scaling for solving certain classes of positive-definite QLS:
in Section 4 we show a method based on an efficient implementation of a
matrix-block-encoding of $A^{-1}$ and in Section 5 a method based on
decomposing the coefficient matrix as $A=LL^{\dagger}$ to effectively
precondition the system. Finally, an outlook of possible future research
directions is given in Section 6.

###  1.1 Positive definite linear systems and quadratic speed-ups

In this work, we investigate the efficiency of QLS solving algorithms
specialised to the case where the coefficient matrix $A\in\mathds{C}^{N\times
N}$ is Hermitian and positive definite (PD). We will call this sub-class the
positive-definite Quantum Linear System (PD-QLS) problem.

A first reason to focus specifically on the PD case is that in the classical
setting several problems of practical relevance are formulated as linear
systems involving PD coefficient matrices. A second important reason is that
the few fully worked-out examples in the literature providing strong evidence
of polynomial quantum speed-ups for “natural” QLS problems involve PD
coefficient matrices [14, 17]. In particular, the discretization of a partial
differential equation (PDE) having a positive-definite kernel (such as, e.g.,
Poisson’s equation) results in a large linear system where the coefficient
matrix is positive definite. Montanaro and Pallister perform in Ref. [14] a
detailed analysis of how QLS may be employed to solve the finite-element
formulation of a PDE and show that the linear dependence on the condition
number is the main bottleneck to obtaining large quantum speed-ups. In fact,
the discretization of PDEs for functions defined in $\mathds{R}^{D}$ typically
results in positive-definite linear systems where the condition number scales
as ${\mathcal{O}}\\!\left(N^{2/D}\right)$ [14].

We highlight that one might have reasonably conjectured that it is possible to
have a quadratically better scaling in $\kappa$, the condition number of $A$,
in the PD case. First, note that the runtime lower bound of Ref. [1] is proven
using a special family of matrices that are indefinite (neither positive nor
negative definite) by construction, hence it is not directly applicable to the
PD case. A standard method allows to transform an indefinite linear system
into a PD one, but having a quadratically larger condition number333Namely,
for a given indefinite matrix $A$, the systems $A\textbf{x}=\textbf{b}$ and
$A^{\prime}\textbf{x}=\textbf{b}^{\prime}$, with
$\textbf{b}^{\prime}:=A^{\dagger}\textbf{b}$ and $A^{\prime}:=A^{\dagger}A$,
have the same solution. The matrix $A^{\prime}$ is positive definite and
$\kappa(A^{\prime})=\kappa(A)^{2}$.; hence, the lower bound in Ref. [1]
directly yields a $\sqrt{\kappa}$ runtime lower bound for PD-QLS solvers.
Second, the conjugate gradient (CG) descent method is the most efficient
classical algorithm for solving PD linear systems and requires only
${\mathcal{O}}\big{(}\sqrt{\kappa}\log(1/\varepsilon)\big{)}$ iterations to
converge to the correct solution, up to $\varepsilon$ approximation. Each
iteration consists in the update of some vectors having $N$ entries, where $N$
is the dimension of the linear system, thus resulting in a total runtime in
${\mathcal{O}}\big{(}N\sqrt{\kappa}\log(1/\varepsilon)\big{)}$ [18]. Then, it
might seem plausible that a quantization of the CG method could yield a
quantum algorithm having runtime in
${\mathcal{O}}\big{(}\mathrm{polylog}(N)\sqrt{\kappa}\log(1/\varepsilon)\big{)}$.

This conjectured quadratic speed-up is however not always achievable since, as
we prove in Section 3, PD-QLS solvers have runtime scaling linearly in
$\kappa$ in worst-case problem instances. But this no-go result can be used as
guidance to understand what is preventing us from achieving a better runtime
scaling and, conversely, what additional conditions have to be imposed in
order to achieve a quadratic speed-up in $\kappa$ for PD-QLS solvers.

Method |  Result |  Requirements
---|---|---
Reduction to majority problem |  $\Omega(\kappa)$ query complexity lower bound Proposition 6 |  Access to $A$ via a matrix-block-encoding or via a sparse-oracle access
Block-encoding of $A^{-1}$ |  ${\mathcal{O}}(\sqrt{\kappa})$ query and gate complexity Proposition 12 and Proposition 13 |  Access to normalised matrix-block-encoding of $B=(I-\eta A)$ and $\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|$ is large
Decomposition as $A=LL^{\dagger}$ |  ${\mathcal{O}}(\sqrt{\kappa})$ gate complexity Proposition 16 |  $A$ is the sum of PD local Hamiltonians, b is sparse and $1/\gamma$ in Eq. (112) is small
Table 1: Summary of the main results of this work. The results provided in the
second column of the table hold under the conditions specified in the third
column. The big-$\Omega$ notation is used for runtime lower bounds (in query
complexity), the big-${\mathcal{O}}$ notation for runtime upper bounds
(exhibiting an explicit solving algorithm).

###  1.2 Previous related results

In context of general QLS solvers, i.e. solvers applicable also to indefinite
or non-Hermitian matrices, the best algorithms have a runtime scaling quasi-
linearly in $\kappa$, i.e. scaling as
${\mathcal{O}}\big{(}\kappa\,\mathrm{polylog}(\kappa)\big{)}$, thus almost
saturating the linear lower bound [1]. Note that the original HHL algorithm
has a worse performance, with a runtime scaling as
${\mathcal{O}}(\kappa^{2}/\varepsilon)$ where $\varepsilon$ is the target
precision. The first algorithm to achieve a quasi-linear scaling in $\kappa$
was proposed by Ambainis in Ref. [4], which introduces a technique called
Variable-Time Amplitude Amplification (VTAA) and employs it to optimize to the
HHL algorithm. Subsequently, Childs, Kothari and Somma [5] introduced
polynomial approximations of $A^{-1}$ to exponentially improve the runtime
dependence on the approximation error to
${\mathcal{O}}\big{(}\kappa^{2}\,\mathrm{polylog}(\kappa/\varepsilon)\big{)}$;
they show, furthermore, that the VTAA-based optimization can be used also for
this algorithm, thus yielding a
${\mathcal{O}}\big{(}\kappa\,\mathrm{polylog}(\kappa/\varepsilon)\big{)}$
runtime. Later, Chakraborty et al. showed that also the pseudo-inversion
problem, whereby the matrix $A$ may be non-invertible and even non-square, can
be solved with a runtime in
${\mathcal{O}}\big{(}\kappa\,\sqrt{\gamma}\,\mathrm{polylog}(\kappa/\varepsilon)\big{)}$,
where $\gamma$ parametrises the overlap of b with the subspace where $A$ is
non-singular [19]. Finally, the current state-of-the-art methods for general
QLS solving is given in Ref. [12], which does not rely on VTAA but instead is
based on ideas stemming from adiabatic quantum computation [7], which result
in conceptually simpler algorithms and in a significant improvement of the
polylogarithmic factors.

Furthermore, several specialized quantum algorithms have been introduced with
the scope of more efficiently solving QLS for particular sub-classes of
matrices. A few works, e.g. [17, 20, 21], have investigated the use of
preconditioning to speed-up the solution. The main idea, which is well-known
in classical linear system solving methods, is to look for an invertible
matrix $B$, a so-called _preconditioner_ , such that the matrix $BA$ has a
smaller condition number than $A$, and subsequently solve the equivalent
linear system $BA=B\,\textbf{b}$. An algorithm based on a sparse
preconditioning matrix was introduced in Ref. [17], but it has very little
formal guarantees of performance improvement. Another method based on
circulant preconditioners was presented in Ref. [20], for which it is clearer
how to assess when a runtime improvement can be achieved. Runtime improvements
have been obtained in Ref. [21] applying new preconditioning methods to
Hamiltonians arising in many-body physics. An entirely different approach,
based on hybrid classical-quantum algorithms, has been explored in Ref. [22],
which yields runtime speed-ups for the case where the rows or columns of the
coefficient matrix can be prepared with polylogarithmic-depth quantum
circuits. We also mention the result of Ref. [23], showing that it is possible
to solve QLS for the special class of tridiagonal Toeplitz matrix with a
runtime that scales polylogarithmically in all parameters, condition number
included; note, however, that matrices of this class can be fully specified
with just two real parameters. Finally, the PD-QLS problem has been previously
considered in Ref. [16], where it is suggested that positive-definite systems
could be solved more efficiently than indefinite ones.

###  1.3 Overview of the method based on a matrix-block-encoding of $A^{-1}$

Our first method for solving PD-QLS with improved runtime is based on
implementing as a quantum circuit a unitary matrix ${\mathcal{U}}_{A^{-1}}$
that encodes in a sub-block a matrix proportional to $A^{-1}$, i.e. a so-
called matrix-block-encoding [24, 25]. The solution state
$\left|{A^{-1}\textbf{b}}\right\rangle$ can be subsequently obtained via
matrix-vector multiplication, achieved by applying the circuit encoding
$A^{-1}$ and projecting onto the correct sub-block. This method is analogous
to the one introduced by Childs, Kothari and Somma in Ref. [5], where it is
shown that exponentially precise polynomial approximations of the inverse
function can be constructed, which then allow to implement a matrix-block-
encoding of $A^{-1}$ up to exponentially small error and finally solve the QLS
problem via matrix-vector multiplication. Here, we show that if $A$ is a PD
matrix, an equally good approximation of $A^{-1}$ can be obtained with
polynomials having a quadratically smaller degree, leading to the possibility
of a quadratic speed-up in PD-QLS solving.

More in detail, Ref. [5] considers a Hermitian matrix $A$ whose spectrum is by
assumption contained in the domain
$\mathcal{D}_{\Delta}:=[-1,-\Delta]\cup[+\Delta,+1]$, with $\Delta\leq
1/\kappa$. Then, families of real polynomials $P(x)$ are constructed such that
$|P(x)-1/x|\leq\varepsilon$ on the domain $\mathcal{D}_{\Delta}$ and have a
degree $\ell\in{\mathcal{O}}\big{(}\kappa\log(1/\varepsilon)\big{)}$, see
panel $(a)$ of Figure 1 for an illustrative example. Using either the Linear
Combination of Unitaries (LCU) [26, 27] or the Quantum Signal Processing (QSP)
[28, 24, 25] method, it is possible to implement a matrix-block-encoding of
$P(A)$, up to some rescaling factor $K>0$ such that $P(A)/K$ can be encoded a
sub-block of a unitary matrix, and then we have
$\left|\left|{P(A)-A^{-1}}\right|\right|\,\leq\,\varepsilon$ in operator norm.
The query complexities of LCU and QSP scale at least linearly in the degree of
the polynomial, hence the use of polynomials of low degree is crucial to
construct efficient QLS solvers. Moreover, the normalisation factor $K$ of the
matrix-block-encoding of $P(A)$ scales linearly in $\kappa$ and enters
multiplicatively in the cost of the matrix-vector multiplication step
necessary to produce $\left|{A^{-1}\textbf{b}}\right\rangle$.

In case $A$ is a PD matrix, we can exploit the knowledge that its spectrum is
contained in $\mathcal{D}_{\Delta}^{+}:=[\Delta,1]$ to perform the following
trick: we assume to have access to a _normalised_ matrix-block-encoding of
$B:=I-\eta\,A$, for some constant444We show in Section 4.3 that $B$ can be
efficiently constructed for $\eta=1$ when $A$ is diagonally dominant and for
$\eta=1/J$ when $A$ is the sum of $J$ positive semi-definite local Hamiltonian
terms. $\eta\in(0,2]$, so that the spectrum of $B$ is always contained in the
interval $\mathcal{D}_{B}=[-1,1-\eta\Delta]$. We then construct a polynomial
$P(x)$ that approximates the function $f(x):=1/(1-x)$ on the domain
${\mathcal{D}}_{B}$ up to $\varepsilon$ distance. Using the QSP method to
implement a matrix-block-encoding of the matrix $P(B)$ we then have

$\displaystyle P(B)\approx
f(B)=\frac{1}{I-B}=\frac{1}{I-(I-\eta\,A)}=\frac{1}{\eta}A^{-1}$ (1)

which is the required matrix inverse, up to a $1/\eta$ rescaling factor.
Importantly, our polynomial approximation of $1/(1-x)$ has a degree
$\ell\in{\mathcal{O}}\big{(}\sqrt{\kappa/\eta}\,\log(\kappa/(\eta\varepsilon))\big{)}$,
a quadratically better dependence on $\kappa$ with respect to the
approximation of $1/x$ given in Ref. [5], see panel $(b)$ of Figure 1.
Moreover, the normalisation factor $K$ scales again linearly in $\kappa$,
which is the best dependence achievable.

$(a)$ | $(b)$
---|---
|
Figure 1: In panel $(a)$ are shown two polynomials approximating the function
$1/x$ (in black) chosen from the family of polynomials given in Ref. [5, Lemma
14], using as parameters $b=200,j_{0}=10$ (red curve) and $b=200,j_{0}=20$
(blue curve curve). In panel $(b)$ are shown two polynomials approximating the
function $1/(1-x)$ (in black) chosen from the family given in Eq. (28), using
as parameters $\ell=6,\kappa=15$ (red curve) and $\ell=10,\kappa=9$ (blue
curve). The shaded regions in each panel indicate the intervals where the
polynomial approximation is not accurate.

From a mathematical standpoint, the possibility of a quadratic reduction of
the degree of the approximating polynomial can be interpreted as a consequence
of Bernstein’s inequality [29]. This inequality states that in the class of
real polynomials $P(x)$ of degree $\ell$ such that $|P(x)|\leq 1$ for all
$x\in[-1,+1]$ the derivative $P^{\prime}(x)$ satisfies
$|P^{\prime}(x)|\leq(1-x^{2})^{-1/2}\ell$ for all $-1<x<1$, while we have
$|P^{\prime}(x)|\leq\ell^{2}$ for $x=1$ and $x=-1$ (and these last bounds are
saturated by Chebyschev polynomials). Note that polynomials approximations of
$1/x$ and of $1/(x-1)$ have high derivative close to the singularities,
respectively in $x=0$ and in $x=1$, and because of Bernstein’s inequality the
latter case allows good polynomial approximations having a quadratically lower
degree.

We need next to perform matrix-vector multiplication to obtain the state
$\left|{A^{-1}\textbf{b}}\right\rangle$ and thus solve the QLS; this is
obtained by applying the quantum circuit that encodes the matrix $P(B)\propto
A^{-1}$ to a quantum state of the form
$\left|{0}\right\rangle\left|{\textbf{b}}\right\rangle$ and then post-
selecting the ancilla register to be in $\left|{0}\right\rangle$ or, more
efficiently, using amplitude amplification [30]. The amplitude amplification
step implies a multiplicative overhead of order $1/\kappa$ in the worst case,
yielding a total runtime in
${\mathcal{O}}\big{(}\kappa^{3/2}\log(\kappa/\varepsilon)\big{)}$ for this PD-
QLS solver. However, the efficiency the matrix-vector multiplication depends
on the post-selection success probability and thus on the specific choice for
the vector b. In a best-case scenario the post-selection success probability
is constant and thus the overall runtime of the PD-QLS solver is
${\mathcal{O}}\big{(}\sqrt{\kappa}\log(\kappa/\varepsilon)\big{)}$, while in
the same high-success-probability scenario the QLS solver of Ref. [5] has a
runtime in ${\mathcal{O}}\big{(}\kappa\log(\kappa/\varepsilon)\big{)}$, since
this is the cost of implementing a polynomial approximation of $A^{-1}$ for
indefinite matrices.

We remark that there is a technical assumption that has to be satisfied to
allow the realisation of our method, namely, that it is possible to implement
a _normalised_ matrix block encoding of $B=I-\eta\,A$. This is a non-trivial
task: standard methods could allow us to prepare, for instance, a normalised
matrix-block-encoding of $B/d$ where $d\geq 2$ is the _sparsity_ of $B$, i.e.
the maximum number of non-zero entries in each column of $B$ [5]. Note,
however, that we would then need to implement an approximation of the function
$g(x):=1/(1-xd)$ to obtain $g(B/d)=\frac{1}{\eta}A^{-1}$; the function $g(x)$
has a singularity in $x=1/d$, which is in the interior of the domain
$[-1,+1]$, and thus Bernstein’s inequality precludes us from achieving good
approximations with low-degree polynomials for this function.

We will show that it is possible to implement normalised matrix-block
encodings of $B$ for two special classes of QLS problems: the first is when
$A$ is a diagonally dominant matrix; the second, when $A$ is given as the sum
of local PD Hamiltonian terms, where by “local” we mean that it acts non-
trivially only on a small number of qubits. In these two cases it is therefore
possible to implement our improved PD-QLS solver. We leave the question
whether it is possible in broader generality to prepare normalised block-
encodings as an open research question.

###  1.4 Overview of the method based on the $A=LL^{\dagger}$ decomposition

Our second method for solving PD-QLS with improved runtime is based on finding
a decomposition $A=LL^{\dagger}$, akin to the classical Cholesky decomposition
[31] which exists for all PD matrices, and then use $L^{-1}$ as an efficient
and effective preconditioner. Note, in fact, that
$L^{\dagger}\textbf{x}=\textbf{b}^{\prime}$ for
$\textbf{b}^{\prime}=L^{-1}\textbf{b}$ is a linear system equivalent to the
original one, but the decomposition $A=LL^{\dagger}$ immediately implies
$\kappa(L)=\kappa(L^{\dagger})=\sqrt{\kappa(A)}$ and thus the new system
provably has a quadratically smaller condition number. This decomposition is
similar to a _spectral gap amplification_ of $A$ [32].

The method involves thus two main steps: in the first one we use classical
computation to efficiently obtain a description of a matrix $L$ such that
$LL^{\dagger}=A$ and such that it is possible to efficiently find, using only
classical computation, a description of the vector
$\left|{\textbf{b}^{\prime}}\right\rangle:=\left|{L^{-1}\textbf{b}}\right\rangle$;
in the second one, we use an efficient quantum algorithm, having runtime
quasi-linear in $\kappa$, to solve
$\left|{L^{\dagger}\textbf{x}}\right\rangle=\left|{\textbf{b}^{\prime}}\right\rangle$,
which thus gives
$\left|{\textbf{x}}\right\rangle=\left|{(L^{\dagger})^{-1}L^{-1}\textbf{b}}\right\rangle=\left|{A^{-1}\textbf{b}}\right\rangle$.
Note that the classical descriptions of $L^{\dagger}$ and of
$\left|{\textbf{b}^{\prime}}\right\rangle$ should also allow the efficient
compilation the quantum algorithm used to solve the preconditioned QLS.

The picture is not yet complete, since we actually use a matrix $L$ that is
non-square and thus singular. As a result, the inversion operation has to be
substituted by a pseudo-inversion and the condition number by the _effective_
condition number, the ratio between the largest and the smallest _non-zero_
singular value; when $A$ is invertible the effective condition number of $L$
is quadratically smaller than the condition number of $A$. We also use two
different pseudo-inverses in the classical and in the quantum part of the
computation: in the quantum step the Moore-Penrose pseudo-inverse
$(L^{\dagger})^{+}$ is employed and in the classical preconditioning a
_generalised_ pseudo-inverse $L^{g}$ chosen such that
$(L^{\dagger})^{+}L^{g}=A^{-1}$; thus, they yield the desired solution
$\left|{\textbf{x}}\right\rangle=\left|{A^{-1}\textbf{b}}\right\rangle$ when
applied in sequence. These extensions to non-square matrices and to different
pseudo-inverses are made to give leeway in the design of the classical part of
the algorithm, allowing us to meet the efficiency requirements mentioned
above. Finally, we will employ the QLS solver of [19], which can tackle
pseudo-inversion problems and has a runtime quasi-linear in the effective
condition number.

We show that a fully suitable decomposition of the form $A=LL^{\dagger}$ can
be constructed for the Sum-QLS problem, i.e., when $A$ is a sum of local PD
Hamiltonian terms. In this case, the matrix $L$ is formed by several blocks,
each constructed from a single Hamiltonian term, while the pseudo-inverse
$L^{g}$ is obtained inverting the individual blocks, operations involving only
small matrices and thus classically feasible. We also require that the vector
b is sparse, containing only polynomially many non-zero entries, thus allowing
to efficiently compute the description of
$\left|{\textbf{b}^{\prime}}\right\rangle=\left|{L^{g}\,\textbf{b}}\right\rangle$.

As a final technical caveat, we note that because of the mismatch between the
pseudo-inverse used in the classical and in the quantum part of the algorithm
($L^{g}$ and $(L^{\dagger})^{+}$), the vector $\textbf{b}^{\prime}$ is not
entirely contained in the support of $(L^{\dagger})^{+}$ and thus amplitude
amplification of the component in the correct subspace is required. This
incurs in a multiplicative overhead of a factor $1/\sqrt{\gamma}$, where
$\sqrt{\gamma}>0$ is a known bound for the amplitude of the “good” component
of $\textbf{b}^{\prime}$. This method thus has a provable runtime improvement
over competing QLS solvers whenever $1/\sqrt{\gamma}$ is sufficiently small.

##  2 Notation and definitions

We assume knowledge of the main quantum computation concepts, as given for
instance in Ref. [33]. A quantum computation is described using a Hilbert
space of dimension $2^{n}$ for some $n$, corresponds to a system of $n$
qubits, having a distinguished computational basis.

###  2.1 Linear algebra and asymptotic notation

We consistently denote with $N\in\mathds{N}$ the dimension of the QLS we aim
to solve and we define $n:=\lceil\log_{2}N\rceil$, so that a vector in
$\mathds{C}^{N}$ (possibly padded with zeroes at the end if $N$ is not a power
of two) can be described as pure state of $n$ qubits. For any complex matrix
$A\in\mathds{C}^{N\times M}$ having $N$ rows and $M$ columns we write its
Singular Value Decomposition (SVD) as $A=V\Sigma W^{\dagger}$ where $V$ and
$W$ are unitary matrices of size $N\times N$ and $M\times M$ respectively and
$\Sigma$ is a real non-negative matrix of size $N\times M$ which is uniquely
determined up to reordering of the diagonal entries and contains the singular
values of $A$ on the main diagonal. An Hermitian matrix $A$ is positive
definite if
$\left\langle{\textbf{v}}\right|A\left|{\textbf{v}}\right\rangle>0$ for all
$\left|{\textbf{v}}\right\rangle$ and is is positive semi-definite if
$\left\langle{\textbf{v}}\right|A\left|{\textbf{v}}\right\rangle\geq 0$. The
eigenvalues of $A$ are real, positive (in the definite case) or non-negative
(in the semi-definite case) and coincide with its singular values. For a
general $A$, $\varsigma_{\min},\varsigma_{\max}$ and
$\lambda_{\min},\lambda_{\max}$ denote the minimum and maximum singular values
and eigenvalues, respectively. The Moore-Penrose pseudo-inverse of $A$ is
obtained by applying to the singular values $\varsigma_{i}$ of $A$ the
function $f:\mathds{R}\mapsto\mathds{R}$ defined as $f(x)=1/x$ for $x\neq 0$
and $f(0)=0$. More precisely, for a diagonal matrix $\Sigma$ we define
$\Sigma^{+}=f(\Sigma)$, while for a general matrix $A=V\Sigma W^{\dagger}$ the
pseudo-inverse is given as $A^{+}:=W\Sigma^{+}V^{\dagger}$. Given
$A\in\mathds{C}^{N\times M}$ a _generalised_ pseudo-inverse
$A^{g}\in\mathds{C}^{M\times N}$ is any matrix satisfying the equation
$A\,A^{g}A=A$.

In this work we employ the $\ell^{2}$-norm for vectors
$\left|\left|{\textbf{v}}\right|\right|^{2}:=\sum_{i=1}^{N}\text{v}_{i}^{2}$
and the induced operator norm for matrices
$\left|\left|{A}\right|\right|:=\max_{\textbf{v}\neq
0}\left|\left|{A\textbf{v}}\right|\right|/\left|\left|{\textbf{v}}\right|\right|$.
The condition number of a matrix is given by
$\kappa(A):=\left|\left|{A}\right|\right|||A^{-1}||$. Since we have
$\left|\left|{A}\right|\right|=\varsigma_{\max}(A)$ and
$||A^{-1}||=1/\varsigma_{\min}(A)$, the condition number can be also written
as $\kappa(A)=\frac{\varsigma_{\max}(A)}{\varsigma_{\min}(A)}$. For a singular
matrix $A$ we define the _effective_ condition number to be
$\kappa_{\mathrm{eff}}(A):=\left|\left|{A}\right|\right|||A^{+}||$, which is
equal to the ratio between the largest and the smallest non-zero singular
value of $A$. A Hermitian matrix $A$ is diagonally dominant if $\sum_{j:j\neq
i}|A_{ij}|\leq|A_{ii}|$ for all $i$ and note that $|A_{ii}|=A_{ii}>0$ when $A$
is positive definite.

We use the standard big-${\mathcal{O}}$ and small-$\mathfrak{o}$ notations for
asymptotic scaling, together with the following definitions:
$f(x)\in\Omega(g(x))$ if and only if $g(x)\in{\mathcal{O}}(f(x))$, which is
used to give lower bounds, and
$\Theta(g(x))={\mathcal{O}}(g(x))\,\cap\,\Omega(g(x))$. We also use the
soft-${\mathcal{O}}$ and soft-$\Omega$ notations where
$f(x)\in\widetilde{O}(g(x))$ means
$f(x)\in{\mathcal{O}}\big{(}g(x)\,\text{polylog}\,[g(x)]\big{)}$, and
similarly for $\widetilde{\Omega}(g(x))$, which are used to give more coarse-
grained bounds.

###  2.2 Definition of the Quantum Linear System problem

In this section we introduce the main definitions that are relevant in the
contest of the QLS problem, which is a quantum analogue of the classical
linear algebra problem of solving the system of equations
$A\textbf{x}=\textbf{b}$, having solution $\textbf{x}=A^{-1}\textbf{b}$ when
$A$ is invertible.

We use pure quantum states to encode $N$-dimensional complex vectors. A vector
v enclosed in a bra or in a ket is always assumed to be normalised,
$\big{|}\big{|}\left|{\textbf{v}}\right\rangle\big{|}\big{|}=1$. In particular
we have:

$\displaystyle\left|{\textbf{b}}\right\rangle=\frac{\textbf{b}}{\left|\left|{\textbf{b}}\right|\right|}=\frac{\sum_{i=1}^{N}b_{i}\left|{i}\right\rangle}{\left(\sum_{i=1}^{N}|b_{i}|^{2}\right)^{\\!1/2}}$
(2)
$\displaystyle\left|{A^{-1}\textbf{b}}\right\rangle=\frac{A^{-1}\left|{\textbf{b}}\right\rangle}{\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|}\;.$
(3)

We now give a general definition of the QLS problem. The formulation is
similar to the one provided in Ref. [7] and is intentionally not specifying
the access models employed for the coefficient matrix $A$ and the known-term
vector b, for sake of generality.

######  Definition 1 (Quantum Linear System).

Suppose we have access to a vector
$\textbf{b}\in\mathds{C}^{N}\setminus\\{0\\}$ and to a non-singular matrix
$A\in\mathds{C}^{N\times N}$ (access is given via quantum oracles, or some
kind of implicit or explicit description). We are given two real positive
parameters $\varsigma_{*}$ and $\varsigma^{*}$ such that
$\varsigma_{*}\leq\varsigma_{\min}(A)$ and
$\varsigma_{\max}(A)\leq\varsigma^{*}$, i.e. the singular values of $A$ are
contained in the interval
${\mathcal{D}}_{A}=\big{[}\varsigma_{*},\varsigma^{*}\big{]}$; equivalently,
we are given two parameters $\kappa>1$ and $\alpha>0$ that provide the upper
bounds $\kappa(A)\leq\kappa$ and $\left|\left|{A}\right|\right|\leq\alpha$. We
are also given a target precision $\varepsilon>0$.

The QLS problem then consists in preparing a density matrix
$\rho_{\textbf{x}}$ which is $\varepsilon$-close in trace distance to the
solution vector
$\left|{\textbf{x}}\right\rangle=\left|{A^{-1}\textbf{b}}\right\rangle$; that
is, we require that

$\displaystyle\big{|}\big{|}\,\rho_{\textbf{x}}-\left|{\hskip
0.5pt\textbf{x}\,}\right\rangle\\!\\!\left\langle{\hskip
0.5pt\textbf{x}\,}\right|\big{|}\big{|}_{\mathrm{Tr}}\leq\varepsilon$ (4)

where the trace norm is defined as
$\left|\left|{X}\right|\right|_{\mathrm{Tr}}:=\frac{1}{2}\mathrm{Tr}\left(\sqrt{XX^{\dagger}}\right)$.

######  Definition 2 (Positive-definite Quantum Linear System).

A PD-QLS problem is a QLS problem, as given in Definition 1, where the
coefficient matrix $A$ is Hermitian and positive definite.

We note that a more commonly employed definition requires that the QLS solver
outputs a state $|\widetilde{\textbf{x}}\rangle$ such that
$\big{|}\big{|}\,|\widetilde{\textbf{x}}\rangle-\left|{\textbf{x}}\right\rangle\big{|}\big{|}\leq\varepsilon$.
We prefer to use the trace distance since it is operationally motivated, as it
equals the probability of distinguishing two copies of two quantum states when
optimizing over all possible measurements, see [33, Section 9.2.1]. The trace
distance then also bounds the maximum relative error that can be introduced
when estimating the expectation value
$\mathrm{Tr}(\left|{\textbf{x}}\right\rangle\\!\\!\left\langle{\textbf{x}}\right|M)$
for a given observable $M$ and computing expectation values was the end goal
of Ref. [1]. Finally, for pure normalised states $\left|{\psi}\right\rangle$
and $\left|{\phi}\right\rangle$ the trace distance simplifies as
$d_{\mathrm{Tr}}(\psi,\phi):=\big{|}\big{|}\,\left|{\psi}\right\rangle\\!\\!\left\langle{\psi}\right|-\left|{\phi}\right\rangle\\!\\!\left\langle{\phi}\right|\big{|}\big{|}_{\mathrm{Tr}}=\sqrt{1-|\langle{\psi}|{\phi}\rangle|^{2}}$
and using $|\langle{\psi}|{\phi}\rangle|\geq
1-\frac{1}{2}\big{|}\big{|}\left|{\psi}\right\rangle-\left|{\phi}\right\rangle\big{|}\big{|}^{2}$
we obtain the inequality

$\displaystyle
d_{\mathrm{Tr}}(\psi,\phi)\leq\big{|}\big{|}\left|{\psi}\right\rangle-\left|{\phi}\right\rangle\big{|}\big{|}\,,$
(5)

thus the trace distance between $|\widetilde{\textbf{x}}\rangle$ and
$\left|{\textbf{x}}\right\rangle$ is at least as small as their $\ell^{2}$
distance.

It is customary to assume that $A$ has been rescaled by a factor
$\alpha\geq\left|\left|{A}\right|\right|$, so that we take, without loss of
generality, $\left|\left|{A}\right|\right|\leq 1$ and the only parameter that
needs to be specified is an upper bound to the condition number. This will be
also our convention, unless otherwise specified.

Figure 2: Different access models for QLS algorithms. Panel (a) illustrates
the case where a classical description of $A$ and b (and the target precision
$\varepsilon$) is provided and then used to compile a quantum algorithm
$\mathcal{A}$, solving the QLS for the given $A$ and b. The description of $A$
and b does not need to be fully explicit: it is sufficient that it allows to
efficiently compute the sequence of elementary quantum gates of $\mathcal{A}$.
Panel (b) illustrates the case of a _relativising_ quantum algorithm; in this
case, the sequence of elementary quantum gates only depends on a few
parameters (e.g., $\alpha,\kappa,\varepsilon$ as in Definition 1) while two
fixed sub-routines specify $A$ and b and these sub-routines can be treated as
black-boxes.

###  2.3 Oracles for quantum linear system solving

We now define and discuss a few different access models for $A$ and b, since
the results we present for the PD-QLS solvers crucially depend on which access
model is assumed. A fundamental distinction is between oracular and non-
oracular (a.k.a. relativising and non-relativising) algorithms, see Figure 2.
For instance, in oracular settings it is often possible to establish
unconditional query complexity lower bounds, while in non-oracular settings
non-trivial runtime lower bounds typically can be proven only under some
(reasonable) complexity theory assumptions, such as
$\textsf{P}\neq\textsf{NP}$. Note that most of the literature on QLS assumes
oracular access to $A$ and b [1, 4, 5, 6, 7, 12, 19].

We start defining the access model for b that we assume throughout this paper,
except where differently specified.

######  Definition 3 (State preparation oracle, as in [1]).

Given a vector $\textbf{b}\in\mathds{C}^{N}$ we say that we have quantum
access to a _state preparation oracle_ for b if there is a unitary operator
${\mathcal{U}}_{\textbf{b}}$ such that
${\mathcal{U}}_{\textbf{b}}\left|{0^{n}}\right\rangle=\left|{\textbf{b}^{\prime}}\right\rangle$,
where $\textbf{b}^{\prime}$ is obtained by padding b with zeroes until its
size is a power of $2$.

As already noted by HHL, a general setting where it is possible to efficiently
implement ${\mathcal{U}}_{B}$ is the one presented by Grover and Rudolph [34];
another possibility is that ${\mathcal{U}}_{\textbf{b}}$ is directly encoded
in a qRAM [35], as may be required in quantum machine learning contexts.

Next, we define two models for access to $A$, which we denote as
${\mathcal{P}}_{A}$ and ${\mathcal{U}}_{A}$ and which correspond to a sparse-
matrix-access and matrix-block-encoding, respectively.

######  Definition 4 (Sparse-matrix-access, as in [5]).

Given a Hermitian matrix555Since $A$ is Hermitian, access by rows and by
columns are equivalent. The definition can be extended to non-Hermitian
matrices, but we need to assume access both by rows and by columns. $A$ which
is $d$-sparse (i.e., has at most $d$ non-zero entries in each row and column)
a quantum _sparse-matrix-access_ ${\mathcal{P}}_{A}$ is given by a pair of
oracular functions

$\displaystyle{\mathcal{P}}_{A}:=({\mathcal{P}}_{A}^{pos},{\mathcal{P}}_{A}^{val})$
(6)

where ${\mathcal{P}}_{A}^{pos}$ and ${\mathcal{P}}_{A}^{val}$ specify the
positions of the (potentially) non-zero entries of $A$ and the values of those
entries, respectively, i.e.

$\displaystyle{\mathcal{P}}_{A}^{pos}:\left|{i,\nu}\right\rangle\ \mapsto\
\left|{i,j(i,\nu)}\right\rangle\hskip 15.649pt$
$\displaystyle\qquad\forall~{}i,j\in\\{1,\ldots,N\\}~{}\mathrm{and}~{}\nu\in\\{1,\ldots,d\\}$
(7)
$\displaystyle{\mathcal{P}}_{A}^{val}:\left|{i,j,z}\right\rangle\mapsto\left|{i,j,A_{i,j}\\!\oplus\\!z}\right\rangle$
$\displaystyle\qquad\forall~{}z\in\\{0,1\\}^{*}$ (8)

where $A_{i,j}\\!\oplus\\!z\in\\{0,1\\}^{*}$ denotes a bit string of arbitrary
length that encodes the value $A_{i,j}\in\mathds{C}$.

In order to keep the presentation simple, we assume here and throughout the
paper that numeric representations of complex numbers can be specified exactly
or with a sufficiently high number of digits of precision.

######  Definition 5 (Matrix block encoding, as in [28, 25]).

A unitary operator ${\mathcal{U}}_{A}$ acting on $n+a$ qubits is called an
$(\alpha,a,\varepsilon)$-matrix-block-encoding of a $n$-qubit operator $A$ if
666The circuit ${\mathcal{U}}_{A}$ may also act on other ancilla qubits, which
are in $\left|{0}\right\rangle$ both before and after the application of
${\mathcal{U}}_{A}$. For a given $(\cdot,a,\cdot)$-matrix-block-encoding we
only count the $a$ ancilla qubits that require post-selection to
$\left\langle{0^{a}}\right|$ to obtain the encoding of $A$.

$\displaystyle\big{|}\big{|}\,A-\alpha\,(\left\langle{0^{a}}\right|\otimes
I)\,{\mathcal{U}}_{A}\,(\left|{0^{a}}\right\rangle\otimes
I)\,\big{|}\big{|}\leq\varepsilon$ (9)

which can be expressed also as:

$\displaystyle{\mathcal{U}}_{A}=\left(\begin{array}[]{cc}\widetilde{A}/\alpha&*\\\
&*\end{array}\right)\quad\mathrm{with}~{}\left|\left|{\widetilde{A}-A}\right|\right|\leq\varepsilon\,,$
(12)

where the asterisks ($*$) denote arbitrary matrix blocks of appropriate
dimensions.

We call $\alpha$ the _normalization factor_ of the matrix-block-encoding and
we say in the special case where $\alpha=1$ that ${\mathcal{U}}_{A}$ is a
_normalised_ matrix-block-encoding of $A$.

A technique introduced by Childs allows to implement a $(d,1,0)$-matrix-block-
encoding of $A$, where the normalisation constant $d$ is equal to the sparsity
of $A$, using only a constant number of accesses to ${\mathcal{P}}_{A}$ and
${\mathcal{O}}\big{(}\mathrm{poly}(n)\big{)}$ extra elementary gates, see Ref.
[5] and references therein. In short, we have the reduction:

$\displaystyle{\mathcal{P}}_{A}\ \Longrightarrow\ {\mathcal{U}}_{A}$ (13)

where the arrows means that having access to an oracle of first type allows to
efficiently implement the oracle of the second type.

We also note that other access models to $A$ have been considered in the
literature; for example in Ref. [36] it is assumed that is possible to
efficiently prepare quantum states that are proportional to each one of the
columns of $A$.

###  2.4 Quantum linear systems in non-oracular settings

We consider in Section 5 (and also briefly in Section 4.3) a case that we call
the Sum-of-Hamiltonians QLS (Sum-QLS) problem, which is not formulated as an
oracular algorithm but is based, instead, on a classical description of $A$
and b. In order to obtain efficient Sum-QLS solving algorithms, it is thus
necessary that the descriptions of $A$ and b are compact, depending at most on
${\mathcal{O}}\big{(}\mathrm{poly}(n)\big{)}$ real or complex parameters.

For the known term vector b, we will simply assume that it is a sparse vector
in the computational basis, with at most
${\mathcal{O}}\big{(}\mathrm{poly}(n)\big{)}$ non-zero entries, whose position
is also known. Hence, a fully explicit classical description of b can be
provided and this also enables efficiently implementing a state preparation
circuit ${\mathcal{U}}_{\textbf{b}}$.

For the coefficient matrix $A$ we give a more implicit description: the
entries of $A$ are not specified one-by-one (which would be inefficient, as
the matrix size $N\in\Theta(2^{n})$ is assumed to be very large) but rather
can be computed from only a relatively small set of parameters, scaling
polynomially in the number of qubits. Specifically, we assume that $A$ is
given as the sum of polynomially many local Hamiltonian terms:

$\displaystyle A=\sum_{j=1}^{J}H_{(j)}\qquad\forall j\ H_{(j)}\ \textup{is
positive (semi)-definite,}$ (14)

where the number of terms is $J\in{\mathcal{O}}\big{(}\mathrm{poly}(n)\big{)}$
and each Hamiltonian $H_{(j)}$ acts on a small number of qubits, namely, on at
most ${\mathcal{O}}\big{(}\log(n)\big{)}$ qubits. This case has been
previously considered in Ref. [16].

##  3 Query complexity lower bound

In this section we prove a $\Omega(\kappa)$ lower bound on the runtime of QLS
which is alternative to the ones given by HHL in Ref. [1]. The main innovation
we introduce is that our lower bound applies also to the PD-QLS case, while
the proofs given by HHL only yield a $\widetilde{\Omega}(\sqrt{\kappa})$ lower
bound when specialised to PD matrices. More precisely, we have the following
result.

######  Proposition 6 (Query complexity lower bound).

Consider oracular quantum algorithms that solve the PD-QLS problem as
presented in Definition 2 for different access models to $A$ and b. Namely,
access to b is given via a state preparation oracle
${\mathcal{U}}_{\textbf{b}}$ (Definition 3), while access to $A$ is given
either via a sparse-matrix oracle ${\mathcal{P}}_{A}$ (Definition 4) or via a
matrix-block-encoding ${\mathcal{U}}_{A}$ (Definition 5). Then, PD-QLS solving
algorithms reaching a constant precision $\varepsilon\in{\mathcal{O}}(1)$ have
query complexities
$Q[{\mathcal{U}}_{\textbf{b}}],Q[{\mathcal{U}}_{A}],Q[{\mathcal{P}}_{A}]$ all
in $\Omega\big{(}\min(\kappa,N)\big{)}$.

The proof of these lower bounds is rather technical and can be found in
Appendix A. We now present a weaker result that, however, can be easily proven
as a consequence of the optimality of Grover search [37]; namely, we show that
PD-QLS solving has a linear scaling in $\kappa$ for all $\kappa\leq\sqrt{N}$.

######  Proposition 7.

Consider a PD-QLS problem as presented in Definition 2 and suppose that access
to $A$ is given by a sparse-matrix oracle ${\mathcal{P}}_{A}$ (Definition 4),
with no assumption on the access model for b. Then, a quantum algorithm that
solves the QLS up to any constant precision $\varepsilon\in[0,1)$ must make
$\Omega\big{(}\min(\kappa,\sqrt{N})\big{)}$ accesses to ${\mathcal{P}}_{A}$.

###### Proof.

Consider the search problem of finding an element $j\in\mathcal{S}$, where
${\mathcal{S}}\subseteq\\{1,\ldots,N\\}$ is a subset containing $M$ elements.
The membership of a element $j$ in $\mathcal{S}$ is encoded as a quantum
oracle ${\mathcal{P}}_{\mathcal{S}}$ which flips a ancilla qubit if
$j\in{\mathcal{S}}$ and leaves the ancilla unchanged if
$j\notin{\mathcal{S}}$.

Consider, next, a matrix $A$ that is diagonal (and thus 1-sparse) having
entries

$\displaystyle
A_{j,j}\;=\;\begin{cases}\alpha=\sqrt{\frac{N-M}{N}}&\mathrm{if}~{}j\notin{\mathcal{S}}\\\
\beta=\sqrt{\frac{M}{N}}&\mathrm{if}~{}j\in{\mathcal{S}}\,.\end{cases}$ (15)

The sparse-matrix oracle
${\mathcal{P}}_{A}=({\mathcal{P}}_{A}^{pos},{\mathcal{P}}_{A}^{val})$ for this
diagonal matrix $A$ can be implemented with exactly one access to the
membership oracle ${\mathcal{P}}_{\mathcal{S}}$. In fact,
${\mathcal{P}}_{A}^{pos}$ can be implemented without any access to
${\mathcal{P}}_{\mathcal{S}}$, since it is known that the non-zero entries are
on the diagonal, while ${\mathcal{P}}_{A}^{val}$ can be implemented with a
single access to ${\mathcal{P}}_{\mathcal{S}}$, assuming that $M$ and $N$ are
known: by definition we have
${\mathcal{P}}_{\mathcal{S}}\left|{j,0}\right\rangle=\left|{j,0}\right\rangle$
if $j\notin{\mathcal{S}}$ and
${\mathcal{P}}_{\mathcal{S}}\left|{j,0}\right\rangle=\left|{j,1}\right\rangle$
if $j\in{\mathcal{S}}$, thus it is sufficient to apply on the ancilla system
the transformation $\left|{0}\right\rangle\mapsto\left|{\alpha}\right\rangle$
and $\left|{1}\right\rangle\mapsto\left|{\beta}\right\rangle$, where the
quantum register contains a binary representation of the numbers $\alpha$ and
$\beta$, to implement ${\mathcal{P}}_{A}^{val}$.

Now notice that $A$ is a matrix having condition number
$\kappa=\frac{\alpha}{\beta}=\sqrt{\frac{N-M}{M}}\in\Theta\Big{(}\sqrt{\frac{N}{M}}\,\Big{)}$
assuming $M\leq N/2$. Moreover, $A^{-1}$ is also diagonal, with entries

$\displaystyle(A^{-1})_{j,j}\;=\;\begin{cases}\sqrt{\frac{N}{N-M}}&\text{if}~{}j\notin{\mathcal{S}}\\\
\sqrt{\frac{N}{M}}&\text{if}~{}j\in{\mathcal{S}}\,.\end{cases}$ (16)

Solving the QLS for the known-term vector
$\left|{\textbf{b}}\right\rangle=\left|{\boldsymbol{1}_{N}}\right\rangle=\frac{1}{\sqrt{N}}\sum_{j=1}^{N}\left|{j}\right\rangle$
yields

$\displaystyle\left|{\textbf{x}}\right\rangle\ $ $\displaystyle=\
\frac{A^{-1}\left|{\boldsymbol{1}_{N}}\right\rangle}{\left|\left|{A^{-1}\left|{\boldsymbol{1}_{N}}\right\rangle}\right|\right|}$
(17) $\displaystyle=\
\frac{1}{\sqrt{2N}}\bigg{(}\sqrt{\frac{N}{N-M}}\sum_{j\notin{\mathcal{S}}}\left|{j}\right\rangle+\sqrt{\frac{N}{M}}\sum_{j\in{\mathcal{S}}}\left|{j}\right\rangle\bigg{)}$
(18) $\displaystyle\equiv\
\frac{1}{\sqrt{2}}\Big{(}\left|{j\notin{\mathcal{S}}}\right\rangle+\left|{j\in{\mathcal{S}}}\right\rangle\Big{)},$
(19)

where in the last line we have introduced the normalised vectors
$\left|{j\notin{\mathcal{S}}}\right\rangle:=\frac{1}{\sqrt{N-M}}\sum_{j\notin{\mathcal{S}}}\left|{j}\right\rangle$
and
$\left|{j\in{\mathcal{S}}}\right\rangle:=\frac{1}{\sqrt{M}}\sum_{j\in{\mathcal{S}}}\left|{j}\right\rangle$.
Measuring $\left|{\textbf{x}}\right\rangle$ in the computational basis
therefore solves the search problem with probability $1/2$.

Suppose now that exists a quantum algorithm $\mathcal{A}$ that solves the QLS
problem exactly ($\varepsilon=0$) for PD matrices and that queries
$\mathfrak{o}\big{(}\kappa\big{)}$ times the oracle ${\mathcal{P}}_{A}$.
Applying $\mathcal{A}$ to the diagonal matrix $A$ and
$\left|{\textbf{b}}\right\rangle=\left|{\boldsymbol{1}_{N}}\right\rangle$
would then require only $\mathfrak{o}(\sqrt{N/M}\,)$ calls to
${\mathcal{P}}_{A}$, and hence to ${\mathcal{P}}_{\mathcal{S}}$, in order to
produce $\left|{\textbf{x}}\right\rangle$. This means that using $\mathcal{A}$
we can solve an unstructured search problem using $\mathfrak{o}(\sqrt{N/M}\,)$
queries to ${\mathcal{P}}_{\mathcal{S}}$, violating the optimality of Grover
search.

Next, suppose that the algorithm $\mathcal{A}$ is an approximate PD-QLS
solver, i.e. that it produces a state $\rho_{\textbf{x}}$ such that
$\big{|}\big{|}\,\rho_{\textbf{x}}-\left|{\hskip
0.5pt\textbf{x}\,}\right\rangle\\!\\!\left\langle{\hskip
0.5pt\textbf{x}\,}\right|\big{|}\big{|}_{\mathrm{Tr}}\leq\varepsilon$ for some
constant $\varepsilon<1/2$. Apply a projective measurement to
$\rho_{\textbf{x}}$ that projects it on the space spanned by
$\\{\left|{j}\right\rangle\\}_{j\in{\mathcal{S}}}$ with probability $p$ and
projects it onto the orthogonal subspace with probability $q=1-p$; in the
ideal case ($\varepsilon=0$) we would have $p=q=1/2$. By the operational
definition of the trace distance, the probability distribution $(p,q)$ must be
at most $\varepsilon$-distinguishable from $(1/2,1/2)$, i.e.
$\max\\{|p-1/2|,|q-1/2|\\}\leq\varepsilon$. Thus, the success probability is a
constant $p\geq 1/2-\varepsilon>0$.

Finally, this argument can be extended to any constant precision
$\varepsilon<1$. It is sufficient to define a new diagonal matrix
$\widehat{A}$ by changing the values $\alpha$ and $\beta$ in Eq. (15), so that
the probability $\widehat{p}$ of finding an element $j\in{\mathcal{S}}$ when
measuring
$|\widehat{\textbf{x}}\rangle=|\widehat{A}^{-1}\boldsymbol{1}_{N}\rangle$
satisfies $\widehat{p}>1-\varepsilon$.

∎

Notice that in the proof the vector
$\left|{\textbf{b}}\right\rangle=\left|{\boldsymbol{1}_{N}}\right\rangle$ is
fixed and easy to produce and hence the access model for
$\left|{\textbf{b}}\right\rangle$ plays no role in our reduction. We also
remark that this proof can be straightforwardly modified to prove that the
operation of quantum matrix-vector multiplication (i.e., obtaining a state
proportional to $A\left|{\textbf{b}}\right\rangle$) must also have a linear
cost in $\kappa$. Moreover, since a sparse oracle access ${\mathcal{P}}_{A}$
allows to efficiently implement also a matrix-block encoding
${\mathcal{U}}_{A}$ [5], the same reduction immediately rules out oracular
algorithms that use $\mathfrak{o}(\kappa)$ accesses to ${\mathcal{U}}_{A}$
(for $\kappa\leq\sqrt{N}$). Finally, a simple argument shows that a
$\Omega(\kappa)$ lower bound holds also for the
${\mathcal{U}}_{\textbf{b}}$-query complexity: an initial small difference
between two input states b and $\textbf{b}^{\prime}$ can be magnified
$\kappa$-fold in the corresponding outputs
$\left|{A^{-1}\textbf{b}}\right\rangle$ and
$\left|{A^{-1}\textbf{b}^{\prime}}\right\rangle$, which is impossible unless
one uses at least $\kappa$ accesses to ${\mathcal{U}}_{\textbf{b}}$ [38].

##  4 Method based on low-degree polynomial approximations of $A^{-1}$

We start this section giving a few details on how to use the Quantum Signal
Processing (QSP) method to implement polynomial functions of a matrix that we
can access through a matrix-block-encoding (Section 4.1) and then we provide
the explicit definition of the approximating polynomials of the inverse
function (Section 4.2). Next (Section 4.3) we discuss two cases where we can
implement a matrix-block-encoding of $B=I-\eta\,A$, as required to achieve a
quadratic reduction in the degree of the approximating polynomials. Finally
(Section 4.4) we discuss the cost of matrix-vector multiplication and
summarise the cost of PD-QLS solving via this approach.

###  4.1 The quantum signal processing method

We employ QSP as a tool to implement ${\mathcal{U}}_{A^{-1}}$, a matrix-block-
encoding of $A^{-1}$, assuming that we have access to a normalised matrix-
block-encoding ${\mathcal{U}}_{B}$ of

$\displaystyle B:=I-\eta\,A\,$ (20)

for some $\eta>0$. We assume that the spectrum of $B$ is contained in the
interval $\mathcal{D}_{B}=[-1,1-\eta/\kappa]$, where $\kappa$ is an upper
bound to the condition number of $A$. The QSP method can be stated, already
specialised to our case of interest, as follows [25, Theorem 56].

######  Theorem 8.

Consider a $(\beta,b,\epsilon)$-block-encoding ${\mathcal{U}}_{B}$ of a
Hermitian matrix $B$ and let $P(x)$ be a degree-$\ell$ real polynomial with
$|P(x)|\leq 1/2$ for all $x\in[-1,1]$. Then there is a quantum circuit
${\mathcal{U}}_{P(B/\beta)}$ which is a
$(1,b+2,4\ell\sqrt{\epsilon/\beta})$-block-encoding of $P(B/\beta)$, and
requires $\ell$ applications of $U_{B}$ and $U_{B}^{\dagger}$, a single
application of controlled-$U_{B}$ and ${\mathcal{O}}\big{(}(n+b)\ell\big{)}$
elementary quantum gates. Moreover, the same result holds for polynomials
satisfying $|P(x)|\leq 1$ if $P(x)$ has defined parity, i.e., $P(-x)=P(x)$ or
$P(-x)=-P(x)$.

Importantly, there are explicit classical algorithms that can efficiently
compute a parametrization of $\,{\mathcal{U}}_{P(B/\beta)}$ for any polynomial
$P$ and then compile an explicit quantum circuit that implements it, see Refs.
[39, 40] for the current state-of-the-art.

We now further motivate the need to choose $\beta=1$, i.e., that the matrix-
block-encoding of $B$ is normalised; equivalently, the part proportional to
the identity in the definition (20) must not be rescaled. We remind that, as
presented in Section 1.3, our goal is to implement an polynomial $P(B/\beta)$
approximating $A^{-1}$, that is

$\displaystyle P(B/\beta)\approx
f(B)=\frac{1}{I-B}=\frac{1}{I-(I-\eta\,A)}=\frac{1}{\eta}A^{-1}.$ (21)

Note, however, that in this expression actually we have
$P(x)\approx\frac{1}{1-\beta x}$, a function that has a singularity in
$x=1/\beta\leq 1$. As a consequence of Bernstein’s inequality [29] this
function may have polynomial approximations with quadratically better degree
only when the singularity is in $x=1$, i.e. when we have $\beta=1$.

###  4.2 Polynomial approximation of $1/(1-x)$

In this section we show analytical polynomial approximations of the function
$f(x)=1/(1-x)$ so that we can use the QSP method to implement it for the
matrix $B=I-\eta\,A$ as in Eq. (20). To keep the notation simple we assume
$\eta=1$ and that the spectrum of $A$ is contained in
$\mathcal{D}_{A}=\left[\,\frac{1}{\kappa},\,2\right]$, while we can account
for any value $\eta<1$ simply by rescaling $\kappa$ to $\kappa/\eta$.
Consequently, it is only necessary for the polynomial $P(x)$ to be a good
approximation of our target function in the domain
${\mathcal{D}}_{B}=\left[\,-1,\,1\\!-\\!\frac{1}{\kappa}\,\right]$.

Our starting point will be the polynomial
$\hat{{\mathcal{T}}}_{\ell,\kappa}(x)$, a shifted and rescaled version of
${\mathcal{T}}_{\ell}(x)$, the $\ell$-degree Chebyshev polynomial of the first
kind,

$\displaystyle\hat{{\mathcal{T}}}_{\ell,\kappa}(x):=\frac{{\mathcal{T}}_{\ell}\left(\frac{x+\frac{1}{2\kappa}}{1-\frac{1}{2\kappa}}\right)}{{\mathcal{T}}_{\ell}\left(\frac{1+\frac{1}{2\kappa}}{1-\frac{1}{2\kappa}}\right)}\,,$
(22)

which is the solution of the following minimax problem (see Ref. [3, Theorem
6.25]):

$\displaystyle\hat{{\mathcal{T}}}_{\ell,\kappa}(x)=\underset{\begin{subarray}{c}P\in\mathds{R}_{\ell},\\\
P(1)=1\end{subarray}}{\mathrm{argmin}}\max_{x\in[-1,1-\frac{1}{\kappa}]}\big{|}P(x)\big{|}\,.$
(23)

Chebyschev polynomials satisfy the property that
$|{\mathcal{T}}_{\ell}(x)|\leq 1$ for all $\ell\in\mathds{N}$ and all
$x\in[-1,+1]$, while
${\mathcal{T}}_{\ell}(1+\delta)\geq\frac{1}{2}e^{\ell\sqrt{\delta}}$ for
$0\leq\delta\leq 1/6$, see e.g. [12, Lemma 13]. Using the changes of variables

$\displaystyle
y(x):=\frac{x+\frac{1}{2\kappa}}{1-\frac{1}{2\kappa}}\qquad\quad\delta:=\frac{1+\frac{1}{2\kappa}}{1-\frac{1}{2\kappa}}-1=\frac{1}{\kappa-1/2}$
(24)

we can rewrite the definition in Eq. (22) as

$\displaystyle\hat{{\mathcal{T}}}_{\ell,\kappa}(x)=\frac{{\mathcal{T}}_{\ell}\big{(}y(x)\big{)}}{{\mathcal{T}}_{\ell}\big{(}1+\delta\big{)}}\,.$
(25)

Then, the numerator satisfies
$\big{|}{\mathcal{T}}_{\ell}\big{(}y(x)\big{)}\big{|}\leq 1$ for all
$x\in{\mathcal{D}}_{B}=[-1,1-\frac{1}{\kappa}]$, while the denominator is
${\mathcal{T}}_{\ell}\big{(}1+\delta\big{)}\geq\frac{1}{2}e^{\ell\sqrt{\delta}}$.
This means that it is sufficient to choose
$\ell\geq\frac{1}{\sqrt{\delta}}\log\\!\left(\frac{2}{\epsilon}\right)$, i.e.
$\ell\in\Theta\big{(}\sqrt{\kappa}\log(1/\epsilon)\big{)}$, to obtain
$\big{|}\hat{{\mathcal{T}}}_{\ell,\kappa}(x)\big{|}\leq\epsilon$ on the
interval ${\mathcal{D}}_{B}$.

We then use the following $(2\ell-1)$-degree polynomial as our approximation
of $f(x)=\frac{1}{1-x}$:

$\displaystyle
P_{2\ell-1,\kappa}(x):=\frac{1}{1-x}\left[1-\hat{{\mathcal{T}}}_{\ell,\kappa}(x)\right]^{2}.$
(26)

To see that $P_{2\ell-1,\kappa}(x)$ is indeed a polynomial, note that
$[1-\hat{{\mathcal{T}}}_{\ell,\kappa}(x)]$ has a twofold root in $x=1$, thus
is exactly divisible by $1-x$ and moreover $P_{2\ell-1,\kappa}(1)=0$. This
last property is useful because it allows us to implement the pseudo-inverse
$A^{+}$ for a singular matrix $A$; i.e., in the case in which $A$ has some
eigenvalues that are equal to $0$ (equivalently, $B=I-A$ has eigenvalues equal
to $1$) these will be mapped to $0$ and if moreover all the non-zero
eigenvalues are separated from zero by a gap $1/\kappa$, the the polynomial in
Eq. (26) is a close approximation of the mapping $(1-\lambda)\mapsto
1/\lambda$ for all $\lambda\neq 0$ in the spectrum of $A$. We choose the
degree $\ell\in\Theta\big{(}\sqrt{\kappa}\log(\kappa/\varepsilon)\big{)}$ in
such a way that
$\big{|}\hat{{\mathcal{T}}}_{\ell,\kappa}(x)\big{|}\leq\varepsilon/(3\kappa)$
for all $x\in{\mathcal{D}}_{B}$ and thus we obtain from Eq. (26)

$\displaystyle\left|P_{2\ell-1,\kappa}(x)-\frac{1}{1-x}\right|\;\leq\;\kappa\left(2\,\frac{\varepsilon}{3\kappa}+\frac{\varepsilon^{2}}{9\kappa^{2}}\right)\;\leq\;\varepsilon\qquad\quad\forall
x\in{\mathcal{D}}_{B}\,,$ (27)

that is, we have an $\varepsilon$-close polynomial approximation on the
interval ${\mathcal{D}}_{B}=[-1,1-\frac{1}{\kappa}]$. This directly implies
$\left|\left|{P_{2\ell-1,\kappa}(B)-A^{-1}}\right|\right|\leq\varepsilon$ in
operator norm.

Finally, the polynomial in Eq. (26) has to be normalised so that it becomes
compatible with the QSP method. Therefore we define

$\displaystyle\hat{P}_{2\ell-1,\kappa}(x):=\frac{P_{2\ell-1,\kappa}(x)}{K}\,,\qquad\mathrm{where}~{}K:=\,2\\!\\!\max_{x\in[-1,+1]}|P_{2\ell-1}(x)|\;.$
(28)

With this definition we have $|\hat{P}_{2\ell-1}(x)|\leq\frac{1}{2}$ for
$x\in[-1,+1]$, as required. The normalisation constant satisfies
$K\in\Theta(\kappa)$ for $\ell\in\Omega(\sqrt{\kappa})$, see Appendix B for
the proof of this bound. In conclusion, the QSP method allows us to implement
a $(K,b+2,\varepsilon)$-matrix-block-encoding of $A^{-1}$, where $b$ is the
number of ancilla qubits required in the block-encoding of $B$.

###  4.3 Implementing normalised matrix-block-encodings of $B=I$ \- $\eta A$

In this subsection we show two methods that, under different assumptions,
allow us to efficiently implement a normalised matrix-block-encoding of
$B=I-\eta\,A$. Preliminarily, we remark that this is a non-trivial task, as we
argue with the following three considerations.

First, any Hermitian matrix $M\in\mathds{C}^{N\times N}$ satisfying
$\left|\left|{M}\right|\right|\leq 1$ can be implemented as a normalised sub-
block of a unitary, since an explicit construction is given by [28]

$\displaystyle{\mathcal{U}}_{M}=\left(\begin{array}[]{cc}M&-\sqrt{I-M^{2}}\\\
\sqrt{I-M^{2}}&M\end{array}\right).$ (31)

Implementing the circuit corresponding to ${\mathcal{U}}_{M}$ in general
requires up to ${\mathcal{O}}(N^{2})$ elementary quantum gates [42] and is
thus inefficient; however, efficient constructions are possible in specialised
cases.

Second, standard quantum methods allow us to efficiently implement _sub-
normalised_ matrix-block-encodings of $B$. As a first example, Childs’ quantum
walk operator uses ${\mathcal{O}}(1)$ accesses to a sparse-matrix oracle
${\mathcal{P}}_{B}$ to implement ${\mathcal{U}}_{M}$ for $M=B/d$, where $d$ is
the column-sparsity of $B$ [5]. A second example is to assume that we have
access to ${\mathcal{U}}_{A}$, a block-encoding of $A$ with
$\left|\left|{A}\right|\right|\leq 1$, and then use the LCU lemma [26] to
implement a linear combination of $I$ and $-{\mathcal{U}}_{A}$, yielding a
normalised block-encoding of $p\,I-(1-p)A\equiv p\,B$, for some $p\in(0,1)$.
However, any “black-box” method that aims at amplifying a block-encoding of
$B/\beta$ (with $\beta>1$) to a normalised block-encoding of $B$ is in general
inefficient. This can be proven, for instance, applying the lower bound in
Ref. [25, Theorem 73] to the function $f(x)=\beta\,x$.

Third, it is currently an open problem whether it is possible or not to
implement a normalised matrix-block-encoding ${\mathcal{U}}_{M}$ given access
to a sparse-matrix oracle ${\mathcal{P}}_{M}$ with
$\left|\left|{M}\right|\right|\leq 1$. In absence of general results of this
kind, we then turn to developing specialised methods to efficiently
implementing a normalised block-encoding ${\mathcal{U}}_{B}$, with
$B=I-\eta\,A$, in the cases where $(i)$ $A$ is diagonally dominant or $(ii)$
$A$ is the sum of positive semi-definite local Hamiltonians.

####  4.3.1 Diagonally-dominant coefficient matrix

In this section, we implement a normalised block-encoding of $B=I-A$ employing
the method described in Ref. [25, Lemma 47], which we report here in Lemma 9.
Our construction requires the preparation of some families of states
$\\{\left|{\psi_{i}}\right\rangle\\}_{i},\\{\left|{\phi_{j}}\right\rangle\\}_{j}$
that are well-defined only when $A$ is diagonally-dominant, while attempts at
extending the method to $B=I-\eta\,A$ for non-diagonally-dominant PD matrices
results in non-normalisable states for any $\eta>0$. We also remark that the
problem of solving linear systems involving diagonally dominant PD matrices
(which includes the noteworthy class of Laplacian matrices [43]) is well
studied in classical linear algebra: for these classes of matrices there are
classical algorithms that can solve a linear system substantially faster than
what is possible for more general matrices [44].

######  Lemma 9.

Suppose that we have access to two “state preparation” unitaries $U_{L}$ and
$U_{R}$ (left and right) acting on $a+s$ qubits such that

$\displaystyle U_{L}:\left|{0^{a}}\right\rangle\left|{i}\right\rangle\,$
$\displaystyle\mapsto\left|{\psi_{i}}\right\rangle$ (32) $\displaystyle
U_{R}:\left|{0^{a}}\right\rangle\left|{j}\right\rangle$
$\displaystyle\mapsto\left|{\phi_{j}}\right\rangle,$ (33)

for any $i,j\in\\{1,\ldots,2^{s}\\}$ and for some families of states
$\\{\left|{\psi_{i}}\right\rangle\\}_{i}$ and
$\\{\left|{\phi_{j}}\right\rangle\\}_{j}$. Then, it is immediate to see that
$U_{L}^{\dagger}U_{R}$ is a $(1,a,0)$-matrix-block-encoding of the Gram matrix
$H$ such that $H_{ij}=\langle{\psi_{i}}|{\phi_{j}}\rangle$.

Let $A$ be a Hermitian $d$-sparse diagonally-dominant matrix, i.e.
$\sum_{j\neq i}|A_{ij}|\leq A_{ii}\leq 1$ for all $i$. By Gershgorin theorem
[41], diagonal dominance of a Hermitian matrix is sufficient to guarantee
positive semi-definiteness, i.e. $\lambda_{\min}(A)\geq 0$, and both
$\lambda_{\min}(A)=0$ and $\lambda_{\min}(A)\neq 0$ are possible777If $A$ is
_strictly_ diagonally dominant we have
$\lambda_{\min}(A)\geq\min_{i}\big{\\{}A_{ii}-\sum_{j\neq
i}|A_{ij}|\big{\\}}>0$ and then the condition number is immediately bounded by
$\kappa(A)\leq 1/\lambda_{\min}(A)$ when $\left|\left|{A}\right|\right|\leq
1$.. Consider then the states

$\displaystyle\left|{\psi_{i}}\right\rangle:=\sum_{l\in\mathrm{supp}(A_{i})}\sqrt{\delta_{il}-A_{il}}\left|{l}\right\rangle\;+\;\sqrt{r_{i}}\left|{N+1}\right\rangle$
(34)
$\displaystyle\left|{\phi_{j}}\right\rangle:=\sum_{k\in\mathrm{supp}(A_{j})}\sqrt{\vphantom{A_{il}}\smash[b]{\delta_{jk}-A_{jk}^{*}}}\left|{k}\right\rangle\;+\;\sqrt{r_{j}}\left|{N+1}\right\rangle$
(35)

where $\mathrm{supp}(A_{i})$ are the position of the (at most) $d$ non-zero
entries of the column vector $A_{i}$ and the value $0\leq r_{i}\leq 1$ can be
computed so that the states are normalised, since we have

$\displaystyle|r_{i}|=1-\sum_{l\in\mathrm{supp}(A_{i})}\left|\sqrt{\delta_{il}-A_{il}}\right|^{2}=A_{ii}-\sum_{l\neq
i}\big{|}A_{il}\big{|}\geq 0$ (36)

where we have used $A_{ii}\leq 1$ and the diagonal dominance of $A$. We then
define the following state-preparation unitaries:

$\displaystyle U_{L}:\left|{0^{a}}\right\rangle\left|{i}\right\rangle$
$\displaystyle\;\mapsto\;|0^{b}\rangle\left|{i}\right\rangle\left|{\psi_{i}^{*}}\right\rangle$
(37) $\displaystyle U_{R}:\left|{0^{a}}\right\rangle\left|{j}\right\rangle$
$\displaystyle\;\mapsto\;|0^{b}\rangle\left|{\phi_{j}}\right\rangle\left|{j}\right\rangle,$
(38)

for certain numbers $a$ and $b$ of ancilla qubits initialised in
$\left|{0}\right\rangle$, and where $\left|{\psi_{i}^{*}}\right\rangle$ is the
complex conjugate of the state $\left|{\psi_{i}}\right\rangle$ w.r.t. the
computational basis. Then, $U_{L}^{\dagger}U_{R}$ is a normalised encoding of
the matrix $B=I-A$, as one can verify using $A_{ji}^{*}=A_{ij}$:

$\displaystyle
B_{ij}=\langle{0^{b},i,\psi_{i}^{*}}|{0^{b},\phi_{j},j}\rangle=\underbrace{\sqrt{\vphantom{A_{il}}\smash[b]{\delta_{ji}-A_{ji}^{*}}}}_{\langle{i}|{\phi_{j}}\rangle}\underbrace{\sqrt{\delta_{ij}-A_{ij}}}_{\langle{\psi_{i}^{*}}|{j}\rangle}=\delta_{ij}-A_{ij}\;.$
(39)

The quantum circuit $U_{L}$ can be implemented efficiently (and similarly for
$U_{R}$), as we now show. We use $d$ calls to ${\mathcal{P}}_{A}^{pos}$ to
recover the values $j_{\nu}\equiv j(i,\nu)$ for $\nu\in\\{1,\ldots,d\\}$,
i.e., the positions of all the (potentially) non-zero entries of $A_{i}$. This
corresponds to implementing the following isometry (i.e., a unitary circuit
plus the possibility of adding ancillas):

$\displaystyle\left|{i}\right\rangle\overset{d\times{\mathcal{P}}_{\\!A}^{pos}}{\mapsto}\left|{i}\right\rangle\left|{j_{1},\ldots,j_{d}}\right\rangle.$
(40)

Next, we use $d$ calls to ${\mathcal{P}}_{A}^{val}$ to recover all the values
$A_{i,j}$, i.e.:

$\displaystyle\eqref{eq:comp1}\overset{d\times{\mathcal{P}}_{\\!A}^{val}}{\mapsto}\left|{i}\right\rangle\left|{j_{1},\ldots,j_{d}}\right\rangle\left|{A_{ij_{1}},\ldots,A_{ij_{d}}}\right\rangle.$
(41)

We then use reversible (classical) computation to calculate the numerical
values of all the amplitudes
$\boldsymbol{\psi}^{(i)}:=(\sqrt{-A_{ij_{1}}}^{\,*},\ldots,\sqrt{1-A_{ii}}^{\,*},\ldots,\sqrt{-A_{ij_{d}}}^{\,*},\sqrt{r_{i}}^{\,*})^{T}$:

$\displaystyle\eqref{eq:comp2}\overset{\mathrm{compute}}{\mapsto}\left|{i}\right\rangle\left|{j_{1},\ldots,j_{d}}\right\rangle\left|{A_{ij_{1}},\ldots,A_{ij_{d}}}\right\rangle|\boldsymbol{\psi}^{(i)}\rangle\,.$
(42)

Then, we use a general state preparation quantum circuit which, given a
classical description of the amplitudes of a $(d+1)$-dimensional quantum
state, effectively prepares the corresponding state:

$\displaystyle\eqref{eq:comp3}\overset{\mathrm{prepare}}{\mapsto}\left|{i}\right\rangle\left|{j_{1},\ldots,j_{d}}\right\rangle\left|{A_{ij_{1}},\ldots,A_{ij_{d}}}\right\rangle|\boldsymbol{\psi}^{(i)}\rangle\sum_{\nu=1}^{d+1}\psi_{\nu}^{(i)}\left|{\nu}\right\rangle.$
(43)

Next, we use a single call to ${\mathcal{P}}_{A}^{pos}$ in quantum
superposition to map
$\left|{i,\nu}\right\rangle\mapsto\left|{i,j(i,\nu)}\right\rangle=\left|{i,j_{\nu}}\right\rangle$
and thus we obtain

$\displaystyle\eqref{eq:comp4}\overset{{\mathcal{P}}_{\\!A}^{pos}}{\mapsto}\left|{i}\right\rangle\left|{j_{1},\ldots,j_{d}}\right\rangle\left|{A_{ij_{1}},\ldots,A_{ij_{d}}}\right\rangle|\boldsymbol{\psi}^{(i)}\rangle\sum_{\nu=1}^{d+1}\psi_{\nu}^{(i)}\left|{j_{\nu}}\right\rangle$
(44)

and now note that, adopting the definition $j(i,d+1):=N+1$, on the right-hand
side we have obtained
$\sum_{\nu=1}^{d+1}\psi_{\nu}^{(i)}\left|{j_{\nu}}\right\rangle=\left|{\psi_{i}^{*}}\right\rangle$,
the complex conjugate (w.r.t. the computational basis) of the state defined in
Eq. (34). Finally, we “uncompute” the intermediate registers
$\left|{j_{1},\ldots,j_{d}}\right\rangle$,
$\left|{A_{ij_{1}},\ldots,A_{ij_{d}}}\right\rangle$ and
$|\boldsymbol{\psi}^{(i)}\rangle$, mapping them to
$\left|{0^{b}}\right\rangle$ (where $b$ is the number of left-over ancillas),
performing the steps (40) to (43) in reverse. With a final swapping of
$\left|{i}\right\rangle$ and $\left|{0^{b}}\right\rangle$, we have implemented
the circuit $U_{L}$ given in Eq. (37).

We can now estimate the query and gate cost of implementing $U_{L}$ (and the
cost of implementing $U_{R}$ is the same). Going through the derivation, we
see that $4d+1$ oracle calls to
${\mathcal{P}}_{A}=({\mathcal{P}}_{A}^{pos},{\mathcal{P}}_{A}^{val})$ are
required, that is the query complexity is
$Q[{\mathcal{P}}_{A}]\in{\mathcal{O}}(d)$. Regarding gate complexity, step
(42) requires ${\mathcal{O}}\big{(}d\,\mathrm{poly}(p)\big{)}$ Toffoli gates
(which are universal for classical reversible computation) to perform the
computation up to $p$ digits of precision; moreover, step (43) can be
performed using ${\mathcal{O}}(d)$ control-nots and single-qubit rotations
employing the methods of Ref. [42]. In conclusion, treating the number of
digits of precision as a constant and the single-qubit rotations as exact,
both the query and the gate complexities are in ${\mathcal{O}}(d)$ and we thus
obtain the following result.

######  Proposition 10 (Normalised matrix-block-encoding, diagonally-dominant
case).

Suppose that we have access to a $d$-sparse diagonally-dominant PD matrix
$A\in\mathds{C}^{N\times N}$ via ${\mathcal{P}}_{A}$ as in Definition 4. Then
it is possible to implement a normalised block-encoding of $B=I-A$ with
${\mathcal{O}}(d)$ query and gate complexity, assuming exact single-qubit
rotations and that all arithmetic operations are performed with a constant
number of digits of precision.

We finally remark that, in some cases, it might be possible to implement
$U_{L}$ and $U_{R}$ with a query and gate complexity in
${\mathcal{O}}(\sqrt{d})$ instead of linearly in $d$. First, we may assume
that we have directly access to an oracle ${\mathcal{P}}_{\psi}$ that directly
returns the amplitudes $\psi_{\nu}^{(i)}$ (including the value
$\sqrt{r_{i}}$), instead of needing to compute these vales from the non-zero
entries of $A_{i}$. Second, one can use an algorithm that generalises Grover
search to prepare the state $\left|{\psi_{i}^{*}}\right\rangle$ using
${\mathcal{O}}(\sqrt{d})$ accesses to ${\mathcal{P}}_{\psi}$ [45], and an even
more efficient implementation can be realised using the method of Ref. [46],
which avoids synthesising arithmetical operations and brings about an
improvement of two orders of magnitude over prior works for realistic levels
of precision.

####  4.3.2 Sum of positive semi-definite Hamiltonians

We now consider the case where $A\in\mathds{C}^{N\times N}$ is given by the
sum of positive semi-definite Hamiltonian terms; i.e., we consider the case

$\displaystyle A=\sum_{j=1}^{J}H_{(j)}\qquad\forall j\ H_{(j)}\ \textup{is
positive semi-definite,}$ (45)

which is similar to the case presented in Section 5, but here we allow the
Hamiltonian terms to have eigenvalue zero. We assume that the number $J$ of
Hamiltonian terms scales polynomially in $n=\lceil\log_{2}N\rceil$ and that
each Hamiltonian term $H_{(j)}$ is local, i.e., that it acts upon a small
number of qubits; specifically, we require that the set
${\mathcal{S}}_{j}\subseteq\\{1,\ldots,n\\}$ of qubits upon which $H_{j}$ acts
non-trivially satisfies $|{\mathcal{S}}_{j}|\leq s\in{\mathcal{O}}(\log n)$
for all $j$. Each $H_{(j)}$ can thus be expressed as888 The correct expression
for an operator $H_{(j)}$ acting on a set
${\mathcal{S}}_{j}\subseteq\\{1,\ldots,n\\}$ of $s$ qubits is  $\displaystyle
H_{(j)}=\sum_{\begin{subarray}{c}a_{1}\cdots a_{s}\in\\{0,1\\}^{s}\\\
b_{1}\cdots b_{s}\in\\{0,1\\}^{s}\end{subarray}}h^{(j)}_{a_{1}\cdots
a_{s},b_{1}\cdots
b_{s}}\bigotimes_{r=1}^{n}\;\mathfrak{O}_{r}\qquad\mathrm{with}~{}\mathfrak{O}_{r}=\begin{cases}I=\left|{0}\right\rangle\\!\\!\left\langle{0}\right|\\!+\\!\left|{1}\right\rangle\\!\\!\left\langle{1}\right|&\mathrm{if}~{}r\notin{\mathcal{S}}_{j}\\\
\left|{a_{r}}\right\rangle\\!\left\langle{b_{r}}\right|&\mathrm{if}~{}r\in{\mathcal{S}}_{j}\end{cases}\;.$
(46)

$\displaystyle H_{(j)}=h_{(j)}\otimes I_{\neg{\mathcal{S}}_{j}}$ (47)

where $h_{(j)}$ is a positive-definite matrix of dimension $2^{s}\times
2^{s}$, which can be fully specified with $2^{{\mathcal{O}}(\log
n)}={\mathcal{O}}(\mathrm{poly}\,n)$ parameters. We also impose
$\left|\left|{h_{(j)}}\right|\right|\leq 2$ for all $j$.

Now we define $w_{(j)}:=I-h_{(j)}$ and note that it is a small matrix, of at
most ${\mathcal{O}}(\mathrm{poly}\,n)$ size, with
$\left|\left|{w_{(j)}}\right|\right|\leq 1$. Then, we can rapidly compute,
with classical algorithms, a unitary extension $u_{(j)}$ for each $w_{(j)}$,
each requiring only one ancilla qubit. Specifically, we define

$\displaystyle
u_{(j)}:=\left(\begin{array}[]{cc}w_{(j)}&-\sqrt{I-w_{(j)}^{2}}\\\
\sqrt{I-w_{(j)}^{2}}&w_{(j)}\end{array}\right)$ (50)

and then we (implicitly) define $W_{(j)}\in\mathds{C}^{N\times N}$ and the
unitary $U_{(j)}\in\mathds{C}^{2N\times 2N}$

$\displaystyle W_{(j)}$ $\displaystyle=w_{(j)}\otimes
I_{\neg{\mathcal{S}}_{j}}$ (51) $\displaystyle U_{(j)}$
$\displaystyle=u_{(j)}\otimes I_{\neg{\mathcal{S}}_{j}}$ (52)

where the interpretations are the same as in Eq. (47). Note that each
$U_{(j)}$ can be efficiently compiled as a quantum circuit: it is sufficient
to determine the gate decomposition of the $(s+1)$-qubit matrix $u_{(j)}$ and
then embed the circuit in a $(n+1)$-qubit quantum register. We assume that the
ancilla qubit used for the extension given in Eq. (50) is one and shared
across all circuits $U_{(j)}$.

We then employ the LCU lemma [26] as follows. Given access to the controlled
version of each circuit $U_{(j)}$, it is possible to efficiently implement the
multi-controlled unitary

$\displaystyle
U_{\textsc{Select}}=\sum_{j=1}^{J}\left|{j}\right\rangle\\!\\!\left\langle{j}\right|_{c}\otimes
U_{(j)}$ (53)

where the subscript $c$ denotes the control register. Then, defining a unitary
$\mathrm{Had}$ that acts as
$\mathrm{Had}\left|{0}\right\rangle_{c}=\sum_{j=1}^{J}\frac{1}{\sqrt{J}}\left|{j}\right\rangle_{c}$,
we obtain:

$\displaystyle\big{(}\left\langle{0}\right|_{c}\mathrm{Had}^{\dagger}\otimes
I\big{)}\,U_{\textsc{Select}}\,\big{(}\mathrm{Had}\left|{0}\right\rangle_{c}\otimes
I\big{)}$ $\displaystyle=\sum_{j=1}^{J}\frac{1}{J}U_{(j)}$ (54)

and further post-selecting to the “top-left” corner of each $U_{(j)}$,
according to Eq. (50) we obtain:

$\displaystyle\big{(}\left\langle{0}\right|\otimes
I\big{)}\,\sum_{j=1}^{J}\frac{1}{J}U_{(j)}\,\big{(}\left|{0}\right\rangle\otimes
I\big{)}\,$
$\displaystyle=\sum_{j=1}^{J}\frac{1}{J}W_{(j)}=\sum_{j=1}^{J}\frac{1}{J}\left(I-H_{(j)}\right)=I-\frac{1}{J}A\,.$
(55)

In conclusion, ${\mathcal{U}}_{B}:=(\mathrm{Had}^{\dagger}\otimes
I)\,U_{\textsc{Select}}\,(\mathrm{Had}\otimes I)$ is a normalised matrix-
block-encoding of $B=I-\eta A$, where the factor $\eta=1/J$ scales, by
assumption, polynomially in $n$. The gate complexity of ${\mathcal{U}}_{B}$
can be estimated as follows. Each $u_{(j)}$ requires ${\mathcal{O}}(2^{2s})$
elementary gates to be implemented [42] and thus the gate complexity of
$U_{(j)}$ is also in ${\mathcal{O}}(2^{2s})$, assuming that two-qubit gates
can be applied among arbitrary pairs of qubits. Then, $U_{\textsc{Select}}$
has gate complexity scaling as the sum of the complexities of the individual
$U_{(j)}$ [27] and is thus in ${\mathcal{O}}(J2^{2s})$. The complexity of
$\mathrm{Had}$ is ${\mathcal{O}}(\log J)$, a sub-leading additive term that
can be neglected in the asymptotic gate complexity of ${\mathcal{U}}_{B}$.
Hence we have the following result.

######  Proposition 11 (Normalised matrix-block-encoding, Sum-of-Hamiltonians
case).

Suppose that we have an explicit classical description of positive semi-
definite matrices $h_{(j)}\in\mathds{C}^{2^{s}\times 2^{s}}$ with
$j\in\\{1,\ldots,J\\}$ and consider the positive semi-definite Hamiltonian
terms $H_{(j)}$ each obtained by applying $h_{(j)}$ to a subset of qubits
${\mathcal{S}}_{j}\subseteq\\{1,\ldots,n\\}$, as given in Eq. (47). Consider
then a coefficient matrix $A\in\mathds{C}^{2^{n}\times 2^{n}}$ as given in Eq.
(45). Then it is possible to implement a normalised block-encoding of
$B=I-\frac{1}{J}A$ with a gate complexity in ${\mathcal{O}}(J2^{s})$, assuming
exact single-qubit rotations.

###  4.4 From matrix inversion to solving the quantum linear system problem

Suppose now that we have a matrix-block-encoding of

$\displaystyle\hat{P}_{2\ell-1,\kappa}(B)\approx\frac{1}{K}\frac{1}{I-B}=\frac{A^{-1}}{\eta\,K}$
(56)

where $K\in\Theta(\kappa/\eta)$ is the normalization factor of the matrix-
block-encoding (recall the rescaling of $\kappa$ to $\kappa/\eta$), upper
bounded by ${\mathcal{O}}(\kappa/\eta)$ as we show in Appendix B. This is
equivalent to say that we have implemented the unitary

$\displaystyle{\mathcal{U}}_{A^{-1}}\approx\left(\begin{array}[]{cc}A^{-1}/(\eta\,K)&*\\\
&*\end{array}\right)$ (59)

where the left-upper block corresponds to having $a$ ancilla qubits in
$\left|{0^{a}}\right\rangle$. Then, one can directly solve a QLS by applying
the unitary quantum circuit ${\mathcal{U}}_{A^{-1}}$ to the vector
$\left|{0^{a}}\right\rangle\left|{\textbf{b}}\right\rangle$ and then post-
select the outcome $\left|{0^{a}}\right\rangle$ on the ancilla system;
however, post-selection might introduce large overheads, since we have

$\displaystyle{\mathcal{U}}_{A^{-1}}:\left|{0^{a}}\right\rangle\left|{\textbf{b}}\right\rangle\quad\mapsto\hskip
22.76219pt\frac{1}{\eta\,K}\hskip 19.91692pt$
$\displaystyle\left|{0^{a}}\right\rangle
A^{-1}\left|{\textbf{b}}\right\rangle+\sqrt{1-\frac{1}{\eta^{2}K^{2}}}\,\big{|}\Psi^{\perp}\big{\rangle}$
(60)
$\displaystyle=\,\frac{\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|}{\eta\,K}$
$\displaystyle\left|{0^{a}}\right\rangle\left|{A^{-1}\textbf{b}}\right\rangle+\sqrt{1-\frac{1}{\eta^{2}K^{2}}}\,\big{|}\Psi^{\perp}\big{\rangle}$
(61)

where $\big{|}\Psi^{\perp}\big{\rangle}$ is a state perpendicular to all
states of the form $\left|{0^{a},\psi}\right\rangle$. Therefore, the
probability of successfully obtaining the state
$\left|{A^{-1}\textbf{b}}\right\rangle$ when post-selecting on the ancilla
measurement is

$\displaystyle
p_{\mathrm{succ}}=\frac{\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|^{2}}{\eta^{2}K^{2}}\,.$
(62)

A PD-QLS solver that prepares a matrix-encoding of $A^{-1}$ and obtains
$\left|{A^{-1}\textbf{b}}\right\rangle$ via post-selection requires
${\mathcal{O}}(1/p_{\mathrm{succ}})$ accesses to ${\mathcal{U}}_{\textbf{b}}$
and ${\mathcal{O}}((2\ell-1)/p_{\mathrm{succ}})$ accesses to
${\mathcal{U}}_{B}$. Recall, $2\ell-1$ is the degree of
$\hat{P}_{2\ell-1,\kappa}(B)$ and
$\ell\in\Theta\\!\left(\sqrt{\frac{\kappa}{\eta}}\log\frac{\kappa}{\eta\varepsilon}\right)$.
The query complexities can be quadratically improved to
${\mathcal{O}}(1/\sqrt{p_{\mathrm{succ}}})$ and
${\mathcal{O}}((2\ell-1)/\sqrt{p_{\mathrm{succ}}})$, respectively, using
amplitude amplification [30]. Having implemented an approximate matrix-
encoding of $\widetilde{A}^{-1}$ that is $\varepsilon$-close to $A^{-1}$
(using Definition 5), the output state
$|{\widetilde{A}^{-1}\textbf{b}}\rangle=\widetilde{A}^{-1}\left|{\textbf{b}}\right\rangle/||\widetilde{A}^{-1}\left|{\textbf{b}}\right\rangle\\!||$
satisfies [5, Proposition 9]

$\displaystyle\left|\left|{\big{|}\widetilde{A}^{-1}\textbf{b}\big{\rangle}-\left|{A^{-1}\textbf{b}}\right\rangle}\right|\right|\leq
4\,\varepsilon\,.$ (63)

The same inequality holds in trace distance because of (5) and it is true both
when using post-selection and when using amplitude amplification. In
conclusion, we have the following results.

######  Proposition 12 (Complexity of the PD-QLS solver).

Suppose that we have access to a $(1,b,0)$-matrix-block-encoding
${\mathcal{U}}_{B}$ of $B=I-\eta\,A$, where $\eta\in(0,1]$ and
$A\in\mathds{C}^{N\times N}$ is a PD matrix with eigenvalues contained in the
interval ${\mathcal{D}}_{A}=\big{[}\frac{1}{\kappa},2\big{]}$ for some known
value $\kappa>1$; see e.g. Proposition 10 and Proposition 11 for explicit
constructions.

First, using the QSP method of Theorem 8 and the polynomial approximation
given in Eq. (26), it is possible to implement a $(K,b+2,\varepsilon)$-matrix-
block-encoding of $A^{-1}$, where $K\in\Theta(\kappa/\eta)$, and the method
has a query complexity

$\displaystyle Q[{\mathcal{U}}_{B}]$
$\displaystyle\in{\mathcal{O}}\\!\left(\sqrt{\frac{\kappa}{\eta}}\,\log\frac{\kappa}{\eta\,\varepsilon}\right),$
(64)

where $Q[{\mathcal{U}}_{B}]$ denotes the number of accesses to
${\mathcal{U}}_{B}$. Moreover, the algorithm is gate-efficient, i.e. it
requires ${\mathcal{O}}\big{(}\mathrm{poly}(n)\,Q[{\mathcal{U}}_{B}]\big{)}$
extra elementary quantum gates.

Second, suppose that we want to solve a PD-QLS as in Definition 2, where
access to $A$ is given (indirectly) by ${\mathcal{U}}_{B}$ and access to b via
a state preparation oracle ${\mathcal{U}}_{\textbf{b}}$ as in Definition 3.
Then, the QLS can be solved up to precision ${\mathcal{O}}(\varepsilon)$ using
the $(K,b+2,\varepsilon)$-matrix-block-encoding of $A^{-1}$ and employing
amplitude amplification to perform matrix-vector multiplication with constant
success probability. The total query complexities, in terms of accesses to
${\mathcal{U}}_{B}$ and ${\mathcal{U}}_{\textbf{b}}$, are

$\displaystyle Q[{\mathcal{U}}_{\textbf{b}}]$
$\displaystyle\in{\mathcal{O}}\\!\left(\frac{\kappa}{\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|}\right)$
(65) $\displaystyle Q[{\mathcal{U}}_{B}]$
$\displaystyle\in{\mathcal{O}}\\!\left(\sqrt{\frac{\kappa}{\eta}}\frac{\kappa}{\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|}\log\frac{\kappa}{\eta\,\varepsilon}\right)$
(66)

and the algorithm is gate efficient, that is, the gate complexity is in
${\mathcal{O}}\big{(}Q\,\mathrm{poly}(\log Q,\log N)\big{)}$. A quadratic
speed-up in $\kappa$ (up to polylogarithmic factors) is achieved over general
QLS solvers when
$\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|\in{\mathcal{O}}(\kappa)$.

We now proceed to a worst-case, average case, and best-case scenario analysis
of a PD-QLS solver as given in the previous Proposition.

##### Worst-case scenario:

In the worst case we have
$\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|\in{\mathcal{O}}(1)$
and consequently the post-selection success probability is
$p_{\mathrm{succ}}\in\Omega(1/\kappa^{2})$. One can in alternative use
${\mathcal{O}}(\kappa)$ rounds of amplitude amplification to reach a constant
success probability. Using amplitude amplification, the overall query
complexity is in ${\mathcal{O}}(\kappa^{3/2})$, which is an improvement
compared to the ${\mathcal{O}}(\kappa^{2})$ runtime achieved by the HHL
algorithm, but still falls shorts of the $\widetilde{{\mathcal{O}}}(\kappa)$
runtime that can be achieved using more advanced methods such as Variable-Time
Amplitude Amplification (VTAA) [4] or eigenpath transversal [7].

##### Average-case scenario:

We now look at the distribution of runtimes that arises when using a randomly
chosen vector $\left|{\textbf{b}}\right\rangle$. As observed in the work by
Subaşı and Somma [47, Section III.B] one could model the eigenvalues of a
positive-definite matrix $A$ as uniformly distributed over the interval
$[1/\kappa,1]$ while if $\left|{\textbf{b}}\right\rangle$ is chosen from the
outcome of a random quantum circuit its amplitudes are sampled according to a
Porter-Thomas distribution; then,
$\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|\in{\mathcal{O}}(\sqrt{\kappa})$
almost surely in the regime $1\ll\kappa\ll N$. This implies that a randomly
sampled PD-QLS problem (according to the specified distribution) can be solved
almost surely with a query complexity in ${\mathcal{O}}(\kappa)$, employing
amplitude amplification in the post-selection step. This method then matches
(actually, improves by a polylog($\kappa$) factor) the asymptotic runtime of
more sophisticated methods (such as those that employ VTAA or adiabatic
evolution) when considering “typical” instances of PD-QLS.

##### Best-case scenario:

The largest value that
$\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|$ can reach
is $\kappa$. In such case the post-selection probability is constant and the
overall query complexity is in ${\mathcal{O}}(\sqrt{\kappa})$, even without
employing amplitude-amplification. This is a quadratic improvement for PD-QLS
solving over competing methods working for indefinite QLS, since implementing
a block-encoding of $A^{-1}$ for indefinite matrices already requires
${\mathcal{O}}(\kappa)$ oracle calls to ${\mathcal{U}}_{A}$ [5]. We note these
best-case problems almost never occur under the probabilistic model described
before, but real-world problems have intrinsic structure that could make them
depart from the Porter-Thomas distribution and thus have
$\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|\gg\sqrt{\kappa}$:
for instance, this is the case if $\left|{\textbf{b}}\right\rangle$ has
constant overlap with the eigenvector relative to the largest eigenvalue of
$A^{-1}$. Note, finally, that it is not required that we know in advance how
large the success probability is, since by definition
$\left|{0^{a}}\right\rangle$ heralds the success.

Figure 3: Bi-logarithmic plots of the asymptotic runtimes of general QLS
solvers and of our PD-QLS solvers. The plot on the left represents the scaling
of the query complexity for methods based on direct matrix-vector
multiplication (together with amplitude amplification) in terms of the
variable
$\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|^{-1}\in[\,1,\kappa\,]$.
The plot on the right represents the scaling of the query complexity for
methods based on VTAA in terms of the variable
$\Gamma_{A,\textbf{b}}:=\sqrt{\kappa}\,\frac{\left|\left|{A^{-1/2}\left|{\textbf{b}}\right\rangle}\right|\right|}{\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|}\in[\,1,\sqrt{\kappa}\,]$.
In both plots we ignore poly-logarithmic multiplicative factors and we express
all the variables in units of powers of $\kappa$, assuming $\kappa\gg 1$ as
constant. See the main text for details.

###  4.5 Optimization using Variable-Time Amplitude Amplification

In Appendix C we show how to use VTAA to obtain a PD-QLS solver having
improved asymptotic query complexities. We proceed as in Ref. [5, Section 5]
and Ref. [19, Section 3]: first we reformulate our algorithm as a variable-
stopping-time quantum algorithm and afterwards we apply the VTAA optimisation
to improve its runtime. There is, however, a technical hurdle to overcome: all
previous VTAA-based QLS solvers use a phase estimation subroutine having a
${\mathcal{O}}(\kappa)$ runtime; its use would then preclude us from achieving
a runtime sub-linear in $\kappa$. The main new idea we introduce is to replace
phase estimation with efficient “windowing functions” whose implementation via
QSP requires only $\widetilde{{\mathcal{O}}}(\sqrt{\kappa})$ accesses to
${\mathcal{U}}_{B}$.

More precisely, in Ref. [5] the so-called Gapped Phase Estimation (GPE) method
is introduced, with the purpose of reliably selecting eigenvalues of $A$ that
are larger than some value $\delta$. For these eigenvalues it is possible to
implement an approximate inverse at a reduced cost, scaling as
${\mathcal{O}}(1/\delta)$ instead of ${\mathcal{O}}(\kappa)$. Then, a sequence
of increasingly small values of $\delta$ is considered, until $\delta\leq
1/\kappa$, and VTAA is employed to enhance the success probability. Since GPE
has a query complexity in $\Omega(1/\delta)$ and the required precision is
$\delta\leq 1/\kappa$, its complexity is in $\Omega(\kappa)$.

Instead, we introduce a “windowing function” $W_{\epsilon,\delta}(\lambda)$ to
select eigenvalues of $B=I-\eta A$ that satisfy
$\lambda\in[-1+2\delta,1-2\delta]$ and to reject eigenvalues
$\lambda\in[-1,-1+\delta]\cup[1-\delta,+1]$, except for a small error
probability $\epsilon$. Thus, $W_{\epsilon,\delta}(x)$ is a polynomial
$\epsilon$-close to 1 in the center of the interval $[-1,+1]$ and
$\epsilon$-close to 0 near the edges of the interval, with a steep fall around
the points $\pm(1-1.5\delta)$. The intervals where the function derivative is
large (of order $1/\delta$) are very close to the extrema of the interval
$[-1,+1]$. According to Bernstein’s inequality [29] it is not prohibited that
a windowing function could be implemented with a polynomial having a degree
$\ell\in\widetilde{O}(1/\sqrt{\delta})$, a quadratically smaller degree
compared to case where the large derivative is near the center of the
interval. In Appendix C we then show, with an explicit construction, that
windowing polynomial with degree $\ell\in\widetilde{O}(1/\sqrt{\delta})$
indeed can be implemented. We defer to the Appendix for further details.

The end result is summarized in the following Proposition.

######  Proposition 13 (Complexity of PD-QLS with VTAA).

Consider a PD-QLS where we have access to a normalised matrix-encoding of
$B=I-\eta\,A$ and to a state preparation unitary for b. Then, there is a VTAA-
based solver having target precision $\varepsilon$, constant success
probability, and query complexities given by

$\displaystyle
Q[{\mathcal{U}}_{\textbf{b}}]\in{\mathcal{O}}\\!\left(\sqrt{\log(\kappa)}+\frac{\kappa}{\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|}\right)$
(67) $\displaystyle
Q[{\mathcal{U}}_{B}]\in{\mathcal{O}}\\!\left(\sqrt{\frac{\kappa}{\eta}}\,\Gamma_{A,\textbf{b}}\,\mathrm{polylog}(\kappa,{\tilde{\epsilon}}^{-1},\eta^{-1})\right)$
(68) $\displaystyle\mathrm{with}\quad$
$\displaystyle\Gamma_{A,\textbf{b}}:=\sqrt{\kappa}\,\frac{\left|\left|{A^{-1/2}\left|{\textbf{b}}\right\rangle}\right|\right|}{\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|},\qquad{\tilde{\epsilon}}\in{\mathcal{O}}\\!\left(\frac{\varepsilon}{\kappa\sqrt{\log\kappa}}\right)$
(69) $\displaystyle\mathrm{and}\quad$
$\displaystyle\mathrm{polylog}(\kappa,{\tilde{\epsilon}}^{-1},\eta^{-1})=\log^{2}(\kappa)\log^{7/4}({\tilde{\epsilon}}^{-1})\log^{3/2}(\eta^{-1}).$
(70)

Moreover, the algorithm is gate efficient.

Now we discuss the runtime of this VTAA PD-QLS solver.

  * •

We always have $Q[{\mathcal{U}}_{B}]\geq Q[{\mathcal{U}}_{\textbf{b}}]$, thus
the ${\mathcal{U}}_{B}$-complexity is the dominant factor.

  * •

Compared to Proposition 12, the ${\mathcal{U}}_{\textbf{b}}$ query complexity
increases here only by an additive $\sqrt{\log(\kappa)}$ factor, while the
${\mathcal{U}}_{B}$ complexity typically (i.e., for almost all values of
$\Gamma_{A,\textbf{b}}$) has a polynomial improvement. To prove it, note that
$\left|\left|{A^{-1/2}\left|{\textbf{b}}\right\rangle}\right|\right|^{2}=\left\langle{\textbf{b}}\right|A^{-1}\left|{\textbf{b}}\right\rangle\leq\kappa$
and thus
$\Gamma_{A,\textbf{b}}=\sqrt{\kappa}\,\frac{\left|\left|{A^{-1/2}\left|{\textbf{b}}\right\rangle}\right|\right|}{\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|}\leq\frac{\displaystyle\kappa}{\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|}$.

  * •

Compared to the general VTAA-based QLS solver of Ref. [5], our PD-QLS solver
has a polynomial speed-up for almost all values of $\Gamma_{A,\textbf{b}}$,
see the right plot in Figure 3. This is a consequence of Lemma 20, where we
prove that $\Gamma_{A,\textbf{b}}\in[\,1,\sqrt{\kappa}\,]$.

##  5 Method based on a quadratic reduction of the condition number via matrix
decomposition

In this section we start giving some preliminary considerations on the
approach that we are going to present (Section 5.1) and then give a formal
statement of the Sum-QLS problem that we solve (Section 5.2). Next, we
describe a classical pre-processing step that quadratically reduces the
condition number (Section 5.3) and then present the quantum algorithm solving
the pseudo-inversion problem that originates from the preconditioning (Section
5.4). Finally, we estimate the gate complexity of the resulting Sum-QLS solver
(Section 5.5).

###  5.1 General considerations

We present now the general features of any algorithm that solves PD-QLS
exploiting a decomposition of the form $A=LL^{\dagger}$, as already summarised
in Section 1.4. Note that such $L$ exists for any PD matrix $A$ and that it
may not be unique, especially if we allow $L$ to be non-square. The key
property is that any $L$ such that $A=LL^{\dagger}$ satisfies
$\kappa_{\mathrm{eff}}(L)=\sqrt{\kappa(A)}$, hence a system of the form
$L^{\dagger}\textbf{x}=\textbf{b}^{\prime}$ (for any $\textbf{b}^{\prime}$) is
quadratically better conditioned than the original system
$A\textbf{x}=\textbf{b}$.

The method we introduce is based on finding matrices $L\in\mathds{C}^{N\times
M}$ and $L^{g}\in\mathds{C}^{M\times N}$, with $M>N$, such that the
decomposition $A=LL^{\dagger}$ holds and $L^{g}$ satisfies $LL^{g}=I$, i.e. it
is a right pseudo-inverse; moreover, $L^{g}$ is rectangular with more rows
than columns and thus $L^{g}L\neq I$, i.e. $L^{g}$ cannot be a left pseudo-
inverse. We make then the following observations.

  1. 1.

$L^{g}LL^{\dagger}\textbf{x}=L^{g}\textbf{b}$ is a linear system equivalent to
the original one, having the vector $\textbf{x}=A^{-1}\textbf{b}$ as the
unique solution, but with no guarantee that the condition number of
$L^{g}LL^{\dagger}$ is small.

  2. 2.

$L^{\dagger}\textbf{x}=L^{g}\textbf{b}$ is a over-constrained linear system
that may be inequivalent to the original one (since $L^{g}L\neq I$) and
typically has no proper solution x.

  3. 3.

Finding
$\mathrm{argmin}_{\textbf{x}}\left|\left|{\,L^{\dagger}\textbf{x}-L^{g}\textbf{b}\,}\right|\right|$
is a problem equivalent to the original system. The unique solution is
$\textbf{x}=(L^{\dagger})^{+}L^{g}\textbf{b}$, therefore using999By assumption
$A$ is invertible, hence $L^{\dagger}\in\mathds{C}^{M\times N}$ is full-rank,
and thus using the SVD $L^{\dagger}=W\Sigma^{\dagger}V^{\dagger}$ we get
$(LL^{\dagger})^{-1}L=(V\Sigma\Sigma^{\dagger}V^{\dagger})^{-1}V\Sigma
W^{\dagger}=V(\Sigma^{\dagger})^{+}W^{\dagger}=(L^{\dagger})^{+}$.
$(L^{\dagger})^{+}=(LL^{\dagger})^{-1}L=A^{-1}L$ and $LL^{g}=I$ we get the
required result $\textbf{x}=A^{-1}LL^{g}\textbf{b}=A^{-1}\textbf{b}$.

The goal is thus to convert the linear system $A\textbf{x}=\textbf{b}$ into
the linear regression problem
$\mathrm{argmin}_{\textbf{x}}\left|\left|{\,L^{\dagger}\textbf{x}-L^{g}\textbf{b}\,}\right|\right|$,
having solution $\textbf{x}=(L^{\dagger})^{+}L^{g}\textbf{b}$. This is a non-
trivial task, as we need, given access to $A$ via sparse-matrix oracle or via
some succinct description, to construct a suitable access to $L^{\dagger}$
and, moreover, given access to b, to construct a suitable access to
$\textbf{b}^{\prime}:=L^{g}\textbf{b}$. The latter requirement seems
particularly worrisome, since it involves a pseudo-inversion of the
exponentially large matrix $L$. We show, however, that for the Sum-QLS
problem, whereby $A$ is provided as a sum of local PD terms, one can find a
suitable $L^{g}$ for which a compact classical description can be efficiently
computed.

We also remark that a pseudo-inversion problem can be interpreted as a regular
matrix inversion on the subspace where the matrix $L^{\dagger}$ is full-rank;
thus, solving pseudo-inversion entails a larger runtime compared to solving a
standard QLS, since the appropriate subspace has to be selected via projection
or via amplitude amplification. More formally, the operator
$(L^{\dagger})^{+}\\!:\mathds{C}^{M}\rightarrow\mathds{C}^{N}$ (with $M>N$)
has rank equal to $N$ and thus we have the orthogonal decomposition

$\displaystyle\mathds{C}^{M}=\mathrm{supp}\big{(}(L^{\dagger})^{+}\big{)}+\mathrm{ker}\big{(}(L^{\dagger})^{+}\big{)}\qquad\begin{cases}\dim\mathrm{supp}\big{(}(L^{\dagger})^{+}\big{)}=N\\\
\dim\mathrm{ker}\big{(}(L^{\dagger})^{+}\big{)}=M-N\end{cases}$ (71)

where the support is by definition the subspace orthogonal to the kernel.
Then, calling $\Pi$ and $\Pi^{\perp}=I-\Pi$ the orthogonal projectors on the
support and on the kernel of $(L^{\dagger})^{+}$, respectively101010Using the
identity $(L^{\dagger})^{+}=A^{-1}L$ we obtain
$\mathrm{supp}((L^{\dagger})^{+})=\mathrm{supp}(L)$ and moreover
$\mathrm{ker}((L^{\dagger})^{+})=\mathrm{ker}(L)$., we obtain the identity

$\displaystyle(L^{\dagger})^{+}\left|{\textbf{b}^{\prime}}\right\rangle=(L^{\dagger})^{+}\,\Pi\left|{\textbf{b}^{\prime}}\right\rangle+(L^{\dagger})^{+}\,\Pi^{\perp}\left|{\textbf{b}^{\prime}}\right\rangle=(L^{\dagger})^{+}\,\Pi\left|{\textbf{b}^{\prime}}\right\rangle$
(72)

since by definition we have $(L^{\dagger})^{+}\,\Pi^{\perp}=0$. It is then
evident that only the component $\Pi\left|{\textbf{b}^{\prime}}\right\rangle$,
which is in general a sub-normalised quantum state, plays a role in the
pseudo-inversion algorithm, while the orthogonal component
$\Pi^{\perp}\left|{\textbf{b}^{\prime}}\right\rangle$ can be arbitrary.
Therefore, any quantum pseudo-inversion algorithms implicitly requires the
amplification of the $\Pi\left|{\textbf{b}^{\prime}}\right\rangle$ component,
which therefore entails a gate complexity in ${\mathcal{O}}(1/\sqrt{\gamma})$,
for some known lower bound
$\sqrt{\gamma}\leq\left|\left|{\Pi\left|{\textbf{b}^{\prime}}\right\rangle}\right|\right|$.

###  5.2 Problem statement

In the Sum-QLS we assume that the coefficient matrix $A\in\mathds{C}^{N\times
N}$ is given by an _explicit classical description_ , rather than via oracular
access. This allows us to evade the lower bounds given in Proposition 6, since
those bounds are formulated for relativising (i.e. oracular) algorithms. We
assume, specifically, that $A$ has the form [16]

$\displaystyle A=\sum_{j=1}^{J}H_{(j)}\qquad\forall j\ H_{(j)}\ \textup{is
positive definite.}$ (73)

Here we impose that each $H_{(j)}$ is strictly positive definite (rather than
semi-definite) because of a technical condition that will become clear later;
in essence, the expression in Eq. (112) can diverge if any $H_{(j)}$ is
singular, resulting in an infinite runtime. Each term $H_{(j)}$ is a local
Hamiltonian, i.e. it can be expressed as [see also Eq. (46) in the footnote]

$\displaystyle H_{(j)}=h_{(j)}\otimes I_{\neg{\mathcal{S}}_{j}}$ (74)

where each $h_{(j)}$ is an operator acting on a subset
${\mathcal{S}}_{j}\subseteq\\{1,\ldots,n\\}$ of at most $s$ qubits,
corresponding to a matrix of size at most $2^{s}\times 2^{s}$. We assume that
$J$, the number of Hamiltonian terms, and $s$ are “small”, i.e. we take
$J\in{\mathcal{O}}(\mathrm{poly}\,n)$ and $s\in{\mathcal{O}}(\log n)$, and
that we have a complete classical description of each operator $h_{(j)}$ and
of each subset ${\mathcal{S}}_{j}$. The matrix $A$ is fully specified with at
most $J2^{2s}$ real parameters and with $Jn$ boolean values (defining the sets
${\mathcal{S}}_{j}$), i.e. the number of parameters is
$J2^{2s}+Jn\in{\mathcal{O}}(\mathrm{poly}\,n)$.

We require moreover that the known-term vector b is sparse, containing at most
$d_{\textbf{b}}$ non-zero entries in the computational basis, where
$d_{\textbf{b}}$ also scales polynomially in $n$. This implies that a
preparation circuit ${\mathcal{U}}_{\textbf{b}}$ can be given as an explicit
small quantum circuit. This leads to the following definition for the Sum-QLS
problem.

######  Definition 14 (Sum-of-Hamiltonians Quantum Linear System).

A Sum-QLS problem is a PD-QLS as in Definition 2 with the following
restrictions. The coefficient matrix $A\in\mathds{C}^{N\times N}$, for
$N=2^{n}$, is provided as the sum of PD Hamiltonian terms
$A=\sum_{j=1}^{J}H_{(j)}$, where each $H_{(j)}$ acts on at most $s$ qubits;
each $H_{(j)}$ is fully specified as a PD matrix $h_{(j)}$ of size
$2^{s_{j}}\times 2^{s_{j}}$ with $s_{j}\leq s$, together with the set
${\mathcal{S}}_{j}$ of $s_{j}$ qubits on which $H_{(j)}$ acts upon. The vector
$\textbf{b}\in\mathds{C}^{N}$ is $d_{\textbf{b}}$-sparse and the value and
position of each of the $d_{\textbf{b}}$ non-zero entries is provided.

In Appendix D we show that the Sum-QLS problem is $\mathsf{BQP}$-hard, by
adapting a proof given in HHL [1]. That is, we show that any polynomial-time
quantum computation (in the BQP class) can be re-formulated as a Sum-QLS
problem for some artfully constructed coefficient matrix $A$ and known-term
vector b; therefore no polynomial-time classical probabilistic algorithm (in
the BPP class) can solve the Sum-QLS problem111111In this context, a classical
probabilistic algorithm “solves” a QLS problem if it outputs a value
$n\in\\{1,\ldots,N\\}$ with probability approximately equal to $|x_{n}|^{2}$,
the square of the $n$-th amplitude of the quantum state
$\left|{\textbf{x}}\right\rangle=\left|{A^{-1}\textbf{b}}\right\rangle$.
(unless $\textsf{BPP}=\textsf{BQP}$).

###  5.3 Classical pre-processing step

In this section we describe the classical pre-processing step, providing a
quadratic improvement of the condition number. We decompose each matrix
$h_{(j)}$ as

$\displaystyle h_{(j)}=l_{(j)}l_{(j)}^{\dagger}\,,$ (75)

which can be accomplished for example via the Cholesky decomposition [31], in
which case $l_{(j)}$ is a lower triangular matrix. Notice that each matrix
$h_{(j)}$ is a small matrix of size $2^{s}\times
2^{s}\in{\mathcal{O}}(\mathrm{poly}\,n)$ and thus the Cholesky decomposition
can be performed numerically on a classical computer using
${\mathcal{O}}(\mathrm{poly}\,n)$ operations; specifically, Cholesky
factorisation of a $m\times m$ matrix requires $m^{3}/3$ arithmetic
operations, $m^{3}/6$ additions and $m^{3}/6$ multiplications [31]. The total
number of Hamiltonian terms is $J\in{\mathcal{O}}(\mathrm{poly}\,n)$, implying
that the total runtime for performing the Cholesky decomposition for all
Hamiltonian terms is ${\mathcal{O}}(2^{3s}J)$, which also is polynomial in $n$
under our assumptions.

We save the Cholesky decompositions $l_{(j)}$ in a classical memory, storing
${\mathcal{O}}(J2^{2}s)={\mathcal{O}}(\mathrm{poly}\,n)$ complex values, for
later use. These decompositions implicitly define the operators

$\displaystyle L_{(j)}=l_{(j)}\otimes I_{\neg{\mathcal{S}}_{j}}$ (76)

where each $L_{(j)}$ is of size $2^{n}\times 2^{n}$ and where the
interpretation of this equation is the same as in Eq. (74). We now introduce
the rectangular matrix $L\in\mathds{C}^{N\times JN}$, with $N=2^{n}$, given by

$\displaystyle L:=\left(L_{(1)}\ \bigg{|}\ \cdots\ \bigg{|}\ L_{(J)}\right)$
(77)

and thus we have the required decomposition

$\displaystyle
LL^{\dagger}=\sum_{j=1}^{J}L_{(j)}L_{(j)}^{\dagger}=\sum_{j=1}^{J}H_{(j)}=A\,.$
(78)

We can efficiently implement quantum circuits ${\mathcal{P}}_{L}$
(${\mathcal{P}}_{L^{\dagger}}$) that provide sparse-matrix access to $L$
($L^{\dagger}$). To this end, it is sufficient to convert the classical random
access memory that stores the positions and values of the entries of $L$
($L^{\dagger}$) into a qRAM [35]. The scheme presented in Ref. [48, Section
6.3.5] implements this qRAM with a gate complexity in
${\mathcal{O}}(nJ2^{2s})$ and a circuit depth in
${\mathcal{O}}\big{(}\log(J2^{2s})\big{)}={\mathcal{O}}\big{(}s+\log
J\big{)}$.

Since $h_{(j)}$ is by assumption non-singular, each $l_{(j)}$ is a non-
singular lower-triangular matrix and its inverse $l_{(j)}^{-1}$ is an upper-
triangular matrix which we can efficiently compute and store in a classical
memory using a polynomial amount of space. These matrices implicitly define
operators $L_{(j)}^{-1}=l_{(j)}^{-1}\otimes I_{\neg{\mathcal{S}}_{j}}$ such
that $L_{(j)}L_{(j)}^{-1}=I$. We can then define the matrix

$\displaystyle
L^{g}:=\frac{1}{J}\left(\begin{array}[]{c}L_{(1)}^{-1}\\\\[2.84526pt]
\hline\cr\vdots\\\\[2.84526pt] \hline\cr L_{(J)}^{-1}\end{array}\right).$ (82)

which is a generalised right pseudo-inverse of $L$, i.e. $L^{g}$ satisfies the
equation

$\displaystyle L\,L^{g}=\frac{1}{J}\sum_{j=1}^{J}I=I$ (83)

and thus using $(L^{\dagger})^{+}=A^{-1}L$ we get the required relation
$(L^{\dagger})^{+}L^{g}=A^{-1}L\,L^{g}=A^{-1}$. Using a qRAM the gate
complexity of ${\mathcal{P}}_{L^{g}}$ is equal to that of ${\mathcal{P}}_{L}$
and is thus in ${\mathcal{O}}(nJ2^{2s})$.

We finally introduce the quantum state

$\displaystyle\left|{\textbf{b}^{\prime}}\right\rangle:=\left|{L^{g}\,\textbf{b}}\right\rangle.$
(84)

By assumption, b is sparse and has
$d_{\textbf{b}}\in{\mathcal{O}}(\mathrm{poly}\,n)$ non-zero entries, hence the
vector $\textbf{b}^{\prime}$ has sparsity $d_{\textbf{b}^{\prime}}\leq
d_{\textbf{b}}J2^{s}\in{\mathcal{O}}(\mathrm{poly}\,n)$. It is then possible
efficiently classically compute the positions and the values of all the non-
zero entries of $\textbf{b}^{\prime}$, with a gate complexity in
${\mathcal{O}}(n\,d_{\textbf{b}^{\prime}})$, and thus also compute the
normalisation factor $\left|\left|{\textbf{b}^{\prime}}\right|\right|$, with a
gate complexity in ${\mathcal{O}}(d_{\textbf{b}^{\prime}})$. Using the method
described in Ref. [49], we can then efficiently compile a quantum circuit
${\mathcal{U}}_{\textbf{b}^{\prime}}$ that prepares
$\left|{\textbf{b}^{\prime}}\right\rangle$ and has a gate complexity in
${\mathcal{O}}(n\,d_{\textbf{b}^{\prime}})={\mathcal{O}}(n\,d_{\textbf{b}}J2^{s})$,
assuming that all single qubit rotations are performed exactly.

Employing the classical pre-processing described up to now, we can efficiently
implement quantum circuits ${\mathcal{P}}_{L^{\dagger}}$ and
${\mathcal{U}}_{\textbf{b}^{\prime}}$ that act as a sparse access to
$L^{\dagger}$ and state preparation circuit for $\textbf{b}^{\prime}$,
respectively; importantly, the gate complexities of these unitaries are
independent from $\kappa$. We thus have at hand the necessary tools to
implement the quantum pseudo-inversion algorithm that we present in the
upcoming Section 5.4.

###  5.4 Efficient pseudo-inversion quantum algorithm

In this section we look into a quantum algorithm for the linear regression
problem

$\displaystyle\underset{\textbf{x}}{\mathrm{argmin}}\left|\left|{\,L^{\dagger}\textbf{x}-\textbf{b}^{\prime}\,}\right|\right|$
(85)

having solution
$\left|{\textbf{x}}\right\rangle=\left|{(L^{\dagger})^{+}\textbf{b}^{\prime}}\right\rangle=\left|{A^{-1}\textbf{b}}\right\rangle$.
We then employ the quantum pseudo-inversion algorithm of [19, Corollary 31]
which we report here for completeness.

######  Proposition 15 (Complexity of pseudo-inverse state preparation).

Suppose $\widetilde{\kappa}\geq 2$, $\mathcal{L}\in\mathds{C}^{N\times N}$ is
a Hermitian matrix whose non-zero eigenvalues are contained in the domain
$[-1,-1/\widetilde{\kappa}]\cup[1/\widetilde{\kappa},1]$, and $\varepsilon$ is
the target precision. Assume that we have access to
$\mathcal{U}_{\mathcal{L}}$, an $(\alpha,a,\delta)$-matrix-block-encoding of
$\mathcal{L}$ with
$\delta\in\mathfrak{o}\left(\varepsilon/(\widetilde{\kappa}^{2}\log^{3}\frac{\widetilde{\kappa}}{\varepsilon})\right)$
and $a\in\Omega(\log N)$, and to a state preparation oracle
$\mathcal{U}_{\textbf{v}}$ for a vector v such that
$\left|\left|{\Pi_{\mathcal{L}}\left|{\textbf{v}}\right\rangle}\right|\right|\geq\sqrt{\gamma}$,
where $\Pi_{\mathcal{L}}$ is the orthogonal projector onto the support of
$\mathcal{L}^{+}$ and $\gamma$ is a known positive parameter. Then, there is a
(VTAA-based) quantum algorithm that produces a state $\varepsilon$-close to
$\left|{\mathcal{L}^{+}\textbf{v}}\right\rangle$ and has:

$\displaystyle Q\left[{\mathcal{U}}_{\mathcal{L}}\right]\;$
$\displaystyle\in\;{\mathcal{O}}\\!\left(\frac{\alpha}{\sqrt{\gamma}}\,\widetilde{\kappa}\log^{3}(\widetilde{\kappa})\log^{2}(1/\varepsilon)\right)$
(86) $\displaystyle Q\left[{\mathcal{U}}_{\textbf{v}}\right]\;$
$\displaystyle\in\;{\mathcal{O}}\\!\left(\frac{1}{\sqrt{\gamma}}\,\widetilde{\kappa}\log(\widetilde{\kappa})\right)$
(87)

where $Q[{\mathcal{U}}_{\mathcal{L}}]$ and $Q[{\mathcal{U}}_{\textbf{b}}]$ are
the query complexity in terms of access to ${\mathcal{U}}_{\mathcal{L}}$ and
${\mathcal{U}}_{\textbf{b}}$. The algorithm is gate-efficient, only requiring
${\mathcal{O}}\big{(}a\,Q\left[\mathcal{U}_{\mathcal{L}}\right]\big{)}$ extra
elementary gates.

We will apply Proposition 15 using as coefficient matrix $\mathcal{L}$ the
Hermitian extension of $L^{\dagger}$ and
$\widetilde{\kappa}\equiv\sqrt{\kappa}$. That is, we consider the Hermitian
matrix $\mathcal{L}\in\mathds{C}^{(J+1)N\times(J+1)N}$ and vector
$\textbf{v}\in\mathds{C}^{(J+1)N}$

$\displaystyle\mathcal{L}=\left(\begin{array}[]{cc}0&L^{\dagger}\\\
L&0\end{array}\right)\quad\textbf{v}=\left(\begin{array}[]{c}\textbf{b}^{\prime}\\\
0\end{array}\right)$ (92)

which, after pseudo-inversion, results in the vector

$\displaystyle\textbf{x}=\mathcal{L}^{+}\textbf{v}=\left(\begin{array}[]{cc}0&L^{+}\\\
(L^{\dagger})^{+}&0\end{array}\right)\left(\begin{array}[]{c}\textbf{b}^{\prime}\\\
0\end{array}\right)=\left(\begin{array}[]{c}0\\\
(L^{\dagger})^{+}\,\textbf{b}^{\prime}\end{array}\right)$ (99)

which encodes the quantum state
$\left|{J}\right\rangle\otimes\left|{(L^{\dagger})^{+}\,\textbf{b}^{\prime}}\right\rangle$
and the solution is obtained discarding the $(J+1)$-level ancilla system in
the state $\left|{J}\right\rangle$.

###  5.5 Runtime estimation

In this section we look into methods for estimating of the parameters
$\alpha,\kappa$ and $\gamma$ (i.e., the normalisation of the block-encoding of
$L^{\dagger}$, condition number, and overlap with the support space) that
determine the runtime of the pseudo-inversion algorithm of Proposition 15, and
thus determine the overall complexity of the Sum-QLS solver.

First, we can implement sparse-matrix-accesses ${\mathcal{P}}_{L}$ and
${\mathcal{P}}_{L^{\dagger}}$ using the information stored in a qRAM, as
previously explained. Childs’ walk operator [5] then allows to realise a
block-encoding ${\mathcal{U}}_{\mathcal{L}}$ of $\mathcal{L}$ using
${\mathcal{O}}(1)$ accesses to ${\mathcal{P}}_{L}$ and
${\mathcal{P}}_{L^{\dagger}}$. Note that ${\mathcal{U}}_{\mathcal{L}}$ is also
a block-encoding of $L^{\dagger}$ (after a swap of the position of the block).
The normalisation factor of this block-encoding is
$\alpha=J2^{s}\in{\mathcal{O}}(\mathrm{poly}\,n)$, equal to the sparsity of
$\mathcal{L}$.

Second, we can explicitly bound the condition number of $A$ as follows.
Positive-definiteness of the Hamiltonian terms $H_{(j)}$ implies positive-
definiteness of $A$ and, moreover, the smallest and largest eigenvalues of $A$
satisfy $\lambda_{\min}(A)\geq\sum_{j=1}^{J}\lambda_{\min}(h_{(j)})$ and
$\lambda_{\max}(A)\leq\sum_{j=1}^{J}\lambda_{\max}(h_{(j)})$. Since each
$h_{(j)}$ can be efficiently diagonalised, this means that it is possible to
classically compute these values, which then yield the explicit upper bound

$\displaystyle\kappa(A)\leq\frac{\sum_{j=1}^{J}\lambda_{\max}(h_{(j)})}{\sum_{j=1}^{J}\lambda_{\min}(h_{(j)})}\equiv\kappa\,.$
(100)

Tighter bounds to $\kappa(A)$ could also be obtained via more computationally
intensive numerical methods, e.g. by first summing together groups of
Hamiltonian terms and then diagonalizing each sum of Hamiltonians.

We now move on to lower-bounding the value of the overlap parameter

$\displaystyle\big{|}\big{|}\,\Pi_{\mathcal{L}}\left|{\textbf{v}}\right\rangle\big{|}\big{|}=\big{|}\big{|}\,\Pi_{L}\left|{\textbf{b}^{\prime}}\right\rangle\big{|}\big{|},$
(101)

where $\Pi_{\mathcal{L}}$ are the projectors on the supports of
$\mathcal{L}^{+}$ and of $(L^{\dagger})^{+}$, respectively, and note that the
supports of $L$ and $(L^{\dagger})^{+}$ are equal. Using the identity
$\Pi_{L}=L^{+}L$ we thus obtain

$\displaystyle\Pi_{L}=L^{\dagger}A^{-1}L$
$\displaystyle=\sum_{i,j=1}^{J}\left|{i}\right\rangle\\!\left\langle{j}\right|\otimes
L_{(i)}^{\dagger}A^{-1}L_{(j)}\,.$ (102)

Moreover, we have:

$\displaystyle\left|{\textbf{b}^{\prime}}\right\rangle=\left|{L^{g}\,\textbf{b}}\right\rangle=\frac{1}{\sqrt{{\mathcal{N}}}}\sum_{j=1}^{J}\left|{j}\right\rangle\otimes
L_{(j)}^{-1}\left|{\textbf{b}}\right\rangle$ (103)

where the normalisation factor ${\mathcal{N}}$ is given by

$\displaystyle{\mathcal{N}}=\sum_{j=1}^{J}\left|\left|{L_{(j)}^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|^{2}=\sum_{j=1}^{J}\left\langle{\textbf{b}}\right|H_{(j)}^{-1}\left|{\textbf{b}}\right\rangle.$
(104)

Then we can compute

$\displaystyle\Pi_{L}\left|{\textbf{b}^{\prime}}\right\rangle$
$\displaystyle=\frac{1}{\sqrt{{\mathcal{N}}}}\sum_{i,j=1}^{J}\left|{i}\right\rangle\otimes
L_{(i)}^{\dagger}A^{-1}L_{(j)}L_{(j)}^{-1}\left|{\textbf{b}}\right\rangle$
(105)
$\displaystyle=\frac{J}{\sqrt{{\mathcal{N}}}}\sum_{i=1}^{J}\left|{i}\right\rangle\otimes
L_{(i)}^{\dagger}A^{-1}\left|{\textbf{b}}\right\rangle$ (106)

and finally we obtain

$\displaystyle\left|\left|{\Pi_{L}\left|{\textbf{b}^{\prime}}\right\rangle}\right|\right|^{-1}$
$\displaystyle=\frac{\sqrt{{\mathcal{N}}}}{J}\left[\left\langle{\textbf{b}}\right|A^{-1}{\textstyle\sum_{i=1}^{J}}\big{(}L_{(i)}L_{(i)}^{\dagger}\big{)}A^{-1}\left|{\textbf{b}}\right\rangle\right]^{-1/2}$
(107)
$\displaystyle=\frac{1}{J}\,\sqrt{\frac{\sum_{j=1}^{J}\left\langle{\textbf{b}}\right|H_{(j)}^{-1}\left|{\textbf{b}}\right\rangle}{\left\langle{\textbf{b}}\right|A^{-1}\left|{\textbf{b}}\right\rangle}}\;.$
(108)

We now suppose that a value $\gamma>0$ such that
$\left|\left|{\Pi_{\mathcal{L}}\left|{\textbf{v}}\right\rangle}\right|\right|\geq\sqrt{\gamma}$
is known and we remind that the quantum psuedo-inversion algorithm has a
runtime quasi-linear in $1/\sqrt{\gamma}$. We extensively comment on the
values that $\gamma$ can take in order to understand in which cases the Sum-
QLS solver yields an advantage over competing methods.

  1. 1.

We have
$\left|\left|{\Pi_{L}\left|{\textbf{b}^{\prime}}\right\rangle}\right|\right|\leq
1$, and the inequality is saturated when $H_{(j)}=A/J$ for all $j$.

  2. 2.

The bound
$\sum_{j=1}^{J}\left\langle{\textbf{b}}\right|H_{(j)}^{-1}\left|{\textbf{b}}\right\rangle\leq
J\,\lambda_{*}^{-1}$ holds, where
$\lambda_{*}:=\min_{j}\lambda_{\min}(H_{(j)})$. Assuming that
$\lambda_{*}\in\Omega\big{(}\lambda_{\min}(A)/J\big{)}$, i.e. there is no
Hamiltonian term having a minimum eigenvalue significantly smaller than the
average minimum eigenvalue, we obtain:

$\displaystyle\left|\left|{\Pi_{L}\left|{\textbf{b}^{\prime}}\right\rangle}\right|\right|^{-1}$
$\displaystyle\in{\mathcal{O}}\bigg{(}\frac{1}{J}\,\frac{\sqrt{J\lambda_{*}^{-1}}}{\sqrt{\left\langle{\textbf{b}}\right|A^{-1}\left|{\textbf{b}}\right\rangle}}\bigg{)}={\mathcal{O}}\bigg{(}\frac{\sqrt{\kappa(A)}}{\sqrt{\left\langle{\textbf{b}}\right|A^{-1}\left|{\textbf{b}}\right\rangle}}\bigg{)}\;.$
(109)

  3. 3.

The numerator in Eq. (108) can be explicitly calculated, while the denominator
is in general difficult to compute121212One could use techniques related to
amplitude estimation to bound
$\big{|}\big{|}A^{-1/2}\left|{\textbf{b}}\right\rangle\big{|}\big{|}$, but
this operation could be as difficult as solving the QLS in the first place..
However, assuming $\left|\left|{A}\right|\right|\leq 1$, we have the bounds

$\displaystyle\sqrt{\left\langle{\textbf{b}}\right|A^{-1}\left|{\textbf{b}}\right\rangle}=\left|\left|{A^{-1/2}\left|{\textbf{b}}\right\rangle}\right|\right|\in\big{[}\,1,\sqrt{\kappa}\,\big{]}$
(110)

  4. 4.

Importantly, the expression
$\left\langle{\textbf{b}}\right|A^{-1}\left|{\textbf{b}}\right\rangle$ appears
at the denominator, so that a more “ill-conditioned” vector b results in a
larger overlap and thus in a faster Sum-QLS solver.

  5. 5.

As in the analysis of Section 4.4 we can study the runtime in an average-case
scenario. For randomly chosen $A$ and b (sampled according to suitable
probability distributions) we have
$\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|\in\Theta(\sqrt{\kappa(A)})$
almost surely; under the same assumptions, we also have
$\left|\left|{A^{-1/2}\left|{\textbf{b}}\right\rangle}\right|\right|\in\Theta(\kappa(A)^{1/4})$
almost surely. Inserting this estimation in Eq. (109) we have that

$\displaystyle\left|\left|{\Pi_{L}\left|{\textbf{b}^{\prime}}\right\rangle}\right|\right|^{-1}$
$\displaystyle\;\in\;\Theta\Big{(}\kappa(A)^{1/4}\Big{)}$ (111)

holds almost surely. We conclude that for an average-case Sum-QLS problem the
runtime is in
$\widetilde{{\mathcal{O}}}(\sqrt{\kappa/\gamma})=\widetilde{{\mathcal{O}}}(\kappa^{3/4})$,
if $\sqrt{\gamma}$ is a tight lower bound for
$\left|\left|{\Pi_{L}\left|{\textbf{b}^{\prime}}\right\rangle}\right|\right|$.

Summarising, we have the following result.

######  Proposition 16 (Complexity of the Sum-QLS solver).

Consider the Sum-QLS problem with parameters $N,J,s,\kappa$ and
$d_{\textbf{b}}$ as in Definition 14. There is a classical-quantum algorithm
$\mathcal{A}$ solving the Sum-QLS that has the following features.

The first part of $\mathcal{A}$ consists of an efficient classical pre-
processing algorithm, which outputs a description of the quantum circuits
implementing ${\mathcal{U}}_{L^{\dagger}}$ and
${\mathcal{U}}_{\textbf{b}^{\prime}}$; here ${\mathcal{U}}_{L^{\dagger}}$ is a
$(2^{s}J,1,0)$-matrix-block-encoding of $L^{\dagger}$ [as given in Eq. (77)]
with gate complexity in ${\mathcal{O}}(nJ2^{2s})$ and circuit depth in
${\mathcal{O}}\big{(}s+\log J\big{)}$; while
${\mathcal{U}}_{\textbf{b}^{\prime}}$ is a state preparation unitary for
$\textbf{b}^{\prime}$ [as given in Eq. (82)] with gate complexity in
${\mathcal{O}}(n\,d_{\textbf{b}}\,J2^{s})$. By definition, $L^{\dagger}$ and
$\textbf{b}^{\prime}$ satisfy
$(L^{\dagger})^{+}\textbf{b}^{\prime}=A^{-1}\textbf{b}$.

The second part of $\mathcal{A}$ consists of using the quantum pseudo-
inversion algorithm given in Proposition 15, using the unitaries
${\mathcal{U}}_{L^{\dagger}}$ and ${\mathcal{U}}_{\textbf{b}^{\prime}}$ as
sub-routines, in order to produce a $\varepsilon$-close approximation of the
ideal output
$\left|{\textbf{x}}\right\rangle=\left|{A^{-1}\textbf{b}}\right\rangle$. This
algorithm requires the knowledge of a value $\sqrt{\gamma}$ that satisfies
$\sqrt{\gamma}\leq\left|\left|{\Pi_{L}\textbf{b}^{\prime}}\right|\right|$,
where $\Pi_{L}$ is the projector onto the support of $L$, equivalently:

$\displaystyle\frac{1}{J^{2}}\frac{\sum_{j=1}^{J}\left\langle{\textbf{b}}\right|H_{(j)}^{-1}\left|{\textbf{b}}\right\rangle}{\left\langle{\textbf{b}}\right|A^{-1}\left|{\textbf{b}}\right\rangle}\leq\frac{1}{\gamma}\,.$
(112)

Inserting the previously given expressions for the relevant parameters in Eqs.
(86-87) results in

gate complexity
$\displaystyle\in{\mathcal{O}}\\!\left(\big{(}n\,J2^{2s}\big{)}\,\frac{J2^{s}}{\sqrt{\gamma}}\,\kappa\log^{3}(\kappa)\log^{2}(1/\varepsilon)\,+\,\big{(}n\,d_{\textbf{b}}J2^{s}\big{)}\frac{1}{\sqrt{\gamma}}\,\kappa\log(\kappa)\right)$
(113)
$\displaystyle={\mathcal{O}}\\!\left(\sqrt{\frac{\kappa}{\gamma}}\;\mathrm{poly}\Big{(}n,\log(\kappa/\varepsilon)\Big{)}\right)$
(114)

assuming that $J,2^{s}$ and $d_{\textbf{b}}$ have polynomial dependence on
$n=\log_{2}N$. A quadratic speed-up in $\kappa$ (up to polylogarithmic
factors) is achieved over general QLS solvers when $\gamma\in\Omega(1)$.

We remark that the family of Sum-QLS instances where all the parameters in
Proposition 16 scale polynomially (with the promise, in particular, that Eq.
(112) holds for some $\gamma\in{\mathcal{O}}(\mathrm{poly}\,n)$), can be
solved in polynomial time on a quantum computer, as was already shown in Ref.
[16]. This means that the subset of problems having a polynomial scaling of
the parameters, which we denote Sum-QLSpoly, is contained in BQP. Moreover,
the reduction in Appendix D can map polynomial-sized quantum circuits onto an
instance of Sum-QLSpoly, thus showing that Sum-QLSpoly is also BQP-hard. These
two inclusion then show that Sum-QLSpoly is BQP-complete131313While the
original QLS problem was proven to be BQP-complete already in Ref. [1], our
contribution is to show that adding the constraint that $A$ is the sum of
positive-definite local Hamiltonians does not change the complexity class of
the problem. .

##  6 Discussion and outlook

In this work we have presented two algorithms aiming at solving QLS problems
in the case where the coefficient matrix is positive definite and having (for
certain problem instances) a runtime in ${\mathcal{O}}(\sqrt{\kappa})$, a
quadratic improvement compared to what can be obtained using general QLS
solvers. This improvement has the potential of greatly expanding the classes
of problems where quantum computation can provide a quantum speed-up. For
instance, the discretization of partial differential equations in $D$
dimensions results in PD linear system with $\kappa\in{\mathcal{O}}(N^{2/D})$
[14] and thus having a runtime improvement from ${\mathcal{O}}(\kappa)$ to
${\mathcal{O}}(\sqrt{\kappa})$ is crucial to yield a quantum speed-up in the
physically relevant cases $D=2$ and $D=3$. As a second example, it is possible
estimate the hitting time of a Markov chain solving a QLS for the matrix
$A=I-S$, where $S$ is related to the discriminant matrix of the Markov chain
[16]; since $A$ is positive definite and decomposable as a sum of PD local
Hamiltonian terms [16, Appendix A] our second algorithm could be applicable to
this problem.

In the spirit of finding in the near future real-world applications of quantum
algorithms, we note that there is considerable interest in the possibility of
realising QLS solvers in Noisy Intermediate-scale Quantum (NISQ) devices [10,
11] and we argue that some of our results might be implementable in NISQ
devices too. In particular, the crux of our first algorithm is to find a good
polynomial approximation for $A^{-1}$ with a degree in
${\mathcal{O}}(\sqrt{\kappa})$ and then implement it with the quantum signal
processing method [28, 25]. The quadratic reduction in the degree of the
polynomials renders their realisation more easily compatible with the next
generation of quantum processors.

We note that many further improvements and extensions to our algorithms may be
possible. Regarding the first algorithm (Section 4), it would be important to
extend the classes of matrices for which a normalised matrix-block-encoding of
$B=I-\eta\,A$ can be efficiently implemented to make the method more generally
applicable. Regarding the second algorithm (Section 5) we note that the
specific choice of the generalised pseudo-inverse $L^{g}$ in the classical
step results in a ${\mathcal{O}}(1/\sqrt{\gamma})$ multiplicative overhead in
the runtime, where $\gamma$ is given in Eq. (112). An open question is whether
a different choice of the pseudo-inverse could improve, or eliminate
altogether, this overhead. We also mention the possibility that the
decomposition of $A$ as a sum of local PD Hamiltonians could be computed on-
the-fly by the solver, instead of being given as an external input. We note
that if the sparsity pattern satisfies certain conditions (it is a _chordal_
graph) a decomposition $A=LL^{\dagger}$ that does not increase the sparsity of
$A$ exists, and the characterisation given in Ref. [50, Theorem 2.6] could be
employed to compute it.

We finally mention an open research idea that may be worth investigating. The
eigenpath transversal method has been used in some new algorithms to solve the
QLS problem with time complexity in $\widetilde{{\mathcal{O}}}(\kappa)$ [7, 8,
10, 11, 12]; this method is simpler than the Variable-Time Amplitude
Amplification (VTAA) method and also results in (marginally) improved
runtimes, however, it is not directly applicable to solve a pseudo-inversion
problem. It would be interesting to find a way to adapt the eigenpath
transversal method to make it work also in the case where the coefficient
matrix is singular. As a by-product, it could replace the algorithm given in
Ref. [19] as the sub-routine used in our second algorithm to solve the pseudo-
inversion problem, therefore making it more practical.

## Acknowledgments

This work was supported by the Dutch Research Council (NWO/OCW), as part of
the Quantum Software Consortium programme (project number 024.003.037). Some
ideas present in the paper originated from discussions with Anirban Chowdhury
while DO was visiting the Center for Quantum Information and Control (CQuIC).
The authors acknowledge discussion with Markus Mieth, Anderas Spörl, and thank
András Gilyén for reading the manuscript and for giving us several insightful
comments and suggestions.

## Appendices

##  Appendix A Proof of the query complexity lower bound

In this Appendix we give the proof of the query complexity lower bound
presented in Section 3, which we present here in a more extended form.

######  Proposition 17 (Query complexity lower bound).

Consider oracular quantum algorithms that solve the PD-QLS problem as
presented in Definition 2 for different access models to $A$ and b. Namely,
access to b is given via a state preparation oracle
${\mathcal{U}}_{\textbf{b}}$ (Definition 3), while access to $A$ is given
either via a sparse-matrix oracle ${\mathcal{P}}_{A}$ (Definition 4) or via a
matrix-block-encoding ${\mathcal{U}}_{A}$ (Definition 5). Then, PD-QLS solving
algorithms reaching a constant precision $\varepsilon\in{\mathcal{O}}(1)$ have