query complexities
$Q[{\mathcal{U}}_{\textbf{b}}],Q[{\mathcal{U}}_{A}],Q[{\mathcal{P}}_{A}]$ all
in $\Omega\big{(}\min(\kappa,N)\big{)}$. More precisely, we have:

  1. 1.

$Q[{\mathcal{U}}_{b}]\in\Omega\big{(}\min(\kappa,N)\big{)}$, independently
from the access model for $A$;

  2. 2.

$Q[{\mathcal{U}}_{A}]\in\Omega\big{(}\min(\kappa,N)\big{)}$, independently
from the access model for b, when ${\mathcal{U}}_{A}$ is a normalised matrix-
block-encoding;

  3. 3.

$Q[{\mathcal{P}}_{A}]\in\Omega\big{(}\min(\kappa,N)\big{)}$, independently
from the access model for b, when $A$ is a matrix with constant sparsity.

In the main text we prove a weaker result using a reduction to the quantum
search problem, which has a query complexity in $\Omega(\sqrt{N/M})$, where
$M\in[N]:=\\{1,\ldots,N\\}$ is the number of marked elements; here, we use
instead a reduction to a “promise majority” problem, which has a query
complexity in $\Omega(N/M)$, where $M\in[N]$ is the margin of the majority. As
a result, we can prove that solving a PD-QLS has linear scaling of the query
complexity in the condition number for all $\kappa\in{\mathcal{O}}(N)$. To
prove Proposition 17, we first introduce the PromiseMajorityM problem as
follows.

######  Definition 18 (PromiseMajorityM).

Given a vector $y\in\\{0,1\\}^{N}$, a value $M\in[N]$ (we assume for
simplicity that $N+M$ is even) and given the promise that we either have

  * (Case 0)

$y_{i}=0$ for $N/2+M/2$ of the entries

  * (Case 1)

$y_{i}=1$ for $N/2+M/2$ of the entries

the PromiseMajorityM problem consists in determining which of the two is the
case.

We assume that we have access to $y$ via a quantum oracle ${\mathcal{P}}_{y}$
that acts as ${\mathcal{P}}_{y}\left|{i,z}\right\rangle=\left|{i,z\oplus
y_{i}}\right\rangle$ for all $i\in[N]$ and for $z\in\\{0,1\\}$. We also remind
that the two-sided bounded-error quantum query complexity $\mathcal{Q}_{2}$ of
a boolean function is defined as the minimum number of accesses to the input
of the function (i.e., to ${\mathcal{P}}_{y}$) that are necessary to correctly
output the value of the function with probability at least $2/3$, both for the
positive and for the negative instances. Then we have the following Lemma:

######  Lemma 19.

The two-sided bounded-error quantum query complexity $\mathcal{Q}_{2}$ of
PromiseMajorityM, in terms of accesses to ${\mathcal{P}}_{y}$, is
$\mathcal{Q}_{2}(\textsc{PromiseMajority}_{M})\in\Omega(N/M)$.

###### Proof.

This follows immediately from Ref. [51, Corollary 1.2]. ∎

We now show that (relativising) PD-QLS solving algorithms can be used to
compute PromiseMajorityM; the lower bound on the query complexity of
PromiseMajorityM directly translates into a lower bound on the query
complexity of the PD-QLS solvers. We will prove separately the three cases of
Proposition 17, with each proof building upon the previous ones.

###### Proof.

_Case 1._

We assume that $y\in\\{0,1\\}^{N}$ is in the domain of PromiseMajorityM, i.e.
$y$ either contains exactly $N/2+M/2$ zeros or $N/2+M/2$ ones, and we define
$\textbf{b}\in\mathds{C}^{N+1}$:

$\displaystyle\left\\{\begin{array}[]{l}b_{i}=(-1)^{y_{i}}\qquad\qquad\mathrm{for}\
i\in[N]\\\\[5.69054pt] b_{N+1}=\sqrt{N+M}\end{array}\right.$ (A.3)

where the value $b_{N+1}$ is fixed to provide a “phase reference” and avoid
ambiguity on the global sign. We have
$\left|{\textbf{b}}\right\rangle=\textbf{b}/\sqrt{2N+M}$ and
$\left|{\textbf{b}}\right\rangle$ can be implemented by first preparing a
state proportional to $(1,\ldots,1,\sqrt{N+M})$ and then applying the correct
phases; this can be done with one oracle call to each of ${\mathcal{P}}_{y}$
and ${\mathcal{P}}_{y}^{\dagger}$, via the transformations

$\displaystyle\left|{i}\right\rangle\ \overset{\mathrm{ancilla}}{\mapsto}\
\left|{i,0}\right\rangle\ \overset{{\mathcal{P}}_{y}}{\mapsto}\
\left|{i,y_{i}}\right\rangle\ \overset{I\otimes Z}{\mapsto}\
(-1)^{y_{i}}\left|{i,y_{i}}\right\rangle\
\overset{{\mathcal{P}}_{y}^{\dagger}}{\mapsto}\
(-1)^{y_{i}}\left|{i,0}\right\rangle\ \overset{\mathrm{discard}}{\mapsto}\
(-1)^{y_{i}}\left|{i}\right\rangle$ (A.4)

and extended by linearity to superpositions. Next, we introduce the vector
$\boldsymbol{1}_{N}:=(1,\ldots,1)^{T}$ containing $N$ ones and then define

$\displaystyle K^{\prime}\equiv K\oplus
0:=\frac{1}{N}\boldsymbol{1}_{N}\boldsymbol{1}_{N}^{\,T}\oplus 0$ (A.5)

as a matrix of size $(N+1)\times(N+1)$, so that $K^{\prime\,2}=K^{\prime}$,
and finally

$\displaystyle A:=I-(1-\epsilon)K^{\prime}$ (A.6)

where $\epsilon$ is a small parameter that we will define shortly. The matrix
$A$ can be used as a coefficient matrix for a PD-QLS solver since $A$ is
positive definite and $\left|\left|{A}\right|\right|=1$. Moreover, the
condition number of $A$ is exactly $\kappa(A)=1/\epsilon$.

Next we have:

$\displaystyle A^{-1}$
$\displaystyle=I+\sum_{t=1}^{\infty}\big{(}(1-\epsilon)K^{\prime}\big{)}^{t}$
(A.7) $\displaystyle=I+\frac{1-\epsilon}{\epsilon}K^{\prime},$ (A.8)

where the summation converges since
$\left|\left|{(1-\epsilon)K}\right|\right|<1$. Let’s apply $A^{-1}$ to b:

$\displaystyle A^{-1}\textbf{b}$
$\displaystyle=\begin{cases}\textbf{b}+\frac{1-\epsilon}{\epsilon}\frac{M}{N}\,\boldsymbol{1}_{N}^{\prime}&\mathrm{if~{}}y\mathrm{~{}has~{}a~{}majority~{}of~{}0}\\\
\textbf{b}-\frac{1-\epsilon}{\epsilon}\frac{M}{N}\,\boldsymbol{1}_{N}^{\prime}&\mathrm{if~{}}y\mathrm{~{}has~{}a~{}majority~{}of~{}1}\end{cases}$
(A.9)
$\displaystyle=\begin{cases}\textbf{b}+\boldsymbol{1}_{N}^{\prime}&\qquad\ \
\,\mathrm{if~{}}y\mathrm{~{}has~{}a~{}majority~{}of~{}0}\\\
\textbf{b}-\boldsymbol{1}_{N}^{\prime}&\qquad\ \
\,\mathrm{if~{}}y\mathrm{~{}has~{}a~{}majority~{}of~{}1}\end{cases}$ (A.10)

where we have introduced $\boldsymbol{1}_{N}^{\prime}:=(1,\ldots,1,0)^{T}$ and
we have chosen $\epsilon$ so that $\frac{1-\epsilon}{\epsilon}\frac{M}{N}=1$,
giving $\kappa(A)=\frac{1}{\epsilon}=\frac{N+M}{M}$. Introducing a boolean
value $f=\textsc{PromiseMajority}_{M}(y)$, i.e. $f\in\\{0,1\\}$ is equal to
the majority of $y$, we can rewrite the vector $A^{-1}\textbf{b}$ entry-wise
as

$\displaystyle\left\\{\begin{array}[]{llc}{[A^{-1}\textbf{b}]_{i}}&=(-1)^{f}\cdot
2&\qquad\mathrm{if~{}}i\mathrm{~{}is~{}such~{}that~{}}y_{i}=f\\\
{[A^{-1}\textbf{b}]_{i}}&=0&\qquad\mathrm{if~{}}i\mathrm{~{}is~{}such~{}that~{}}y_{i}\neq
f\\\ {[A^{-1}\textbf{b}]_{N+1}}&=\sqrt{N+M}&\end{array}\right.$ (A.14)

Then we have
$\left|{A^{-1}\textbf{b}}\right\rangle=A^{-1}\textbf{b}/\left|\left|{A^{-1}\textbf{b}}\right|\right|$,
with

$\displaystyle\left|\left|{A^{-1}\textbf{b}}\right|\right|^{2}=2^{2}\,\frac{N+M}{2}+\sqrt{N+M}^{\,2}=3(N+M)$
(A.15)

so that we get

$\displaystyle\left|{A^{-1}\textbf{b}}\right\rangle=\sqrt{\frac{1}{3}}\left|{N+1}\right\rangle\,+\,(-1)^{f}\sqrt{\frac{2}{3}}\sum_{i:y_{i}=f}\sqrt{\frac{2}{N+M}}\left|{i}\right\rangle.$
(A.16)

We then perform a projective measurement where one of the possible measurement
outcomes is

$\displaystyle\left|{``+"}\right\rangle:=\sqrt{\frac{1}{2}}\left|{N+1}\right\rangle\,+\,\sqrt{\frac{1}{2}}\sum_{i=1}^{N}\frac{1}{\sqrt{N}}\left|{i}\right\rangle.$
(A.17)

The cases $f=0$ and $f=1$ in Eq. (A.16) can be distinguished with constant
advantage, since:

$\displaystyle\langle{``+"}|{A^{-1}\textbf{b}}\rangle$
$\displaystyle=\sqrt{\frac{1}{6}}\,+\,(-1)^{f}\frac{N+M}{2}\sqrt{\frac{1}{3}}\sqrt{\frac{2}{N(N+M)}}$
(A.18) $\displaystyle=\frac{1}{\sqrt{6}}\left(1+(-1)^{f}\sqrt{1+M/N}\right).$
(A.19)

Note that the two cases can still be distinguished with constant probability
if we replace
$\left|{\textbf{x}}\right\rangle=\left|{A^{-1}\textbf{b}}\right\rangle$ with
any approximation $\rho_{\textbf{x}}$ which is sufficiently close to it.

To summarise, suppose we have to solve a PromiseMajorityM problem and that we
can exploit as a subroutine an oracular quantum algorithm $\mathcal{A}$ that,
given access to ${\mathcal{U}}_{\textbf{b}}$, prepares the state
$\left|{A^{-1}\textbf{b}}\right\rangle$ with sufficiently high precision;
suppose moreover that $\mathcal{A}$ has a query complexity
$Q[{\mathcal{U}}_{\textbf{b}}]=g(\kappa)$, for some function
$g:\mathds{R}^{+}\rightarrow\mathds{N}$. Then, $\mathcal{A}$ can be used to
prepare the state in Eq. (A.16) and solve PromiseMajorityM with constant
distinguishing advantage. The ${\mathcal{P}}_{y}$-query complexity of
$\mathcal{A}$ is
$Q[{\mathcal{P}}_{y}]=2\,Q[{\mathcal{U}}_{\textbf{b}}]=2\,g\big{(}\kappa(A)\big{)}=2\,g\\!\left(\frac{N+M}{M}\right)$.
The lower bound $\mathcal{Q}_{2}(\textsc{PromiseMajority}_{M})\in\Omega(N/M)$
then directly implies $g(\kappa)\in\Omega\big{(}\min(\kappa,N)\big{)}$.

∎

###### Proof.

_Case 2._

We modify the construction given in the previous proof and encode the input
$y$ in the entries of the coefficient matrix, with the goal of showing that
$Q[{\mathcal{U}}_{A}]\in\Omega\big{(}\min(\kappa,N)\big{)}$. To this end, we
define the vector $\textbf{u}\in\mathds{R}^{N+1}$ and a diagonal matrix
$D\in\mathds{R}^{(N+1)\times(N+1)}$

$\displaystyle\left\\{\begin{array}[]{ll}u_{i}=1&\mathrm{for}\
i\in[N]\\\\[5.69054pt] u_{N+1}=\sqrt{N+M}&\end{array}\right.\hskip
28.45274pt\left\\{\begin{array}[]{ll}D_{i,i}=(-1)^{y_{i}}&\mathrm{for}\
i\in[N]\\\\[5.69054pt] D_{N+1,N+1}=1&\end{array}\right.$ (A.24)

and notice the vector b in Eq. (A.3) satisfies $\textbf{b}=D\textbf{u}$. We
also define the coefficient matrix $A^{\prime}$

$\displaystyle A^{\prime}:=DAD$ (A.25)

where $A$ is given in Eq. (A.6). Note that $D$ is unitary and self-inverse,
hence $A^{\prime}$ is positive definite, $\kappa(A^{\prime})=\kappa(A)$, and
moreover $A^{\prime-1}=D^{-1}A^{-1}D^{-1}=DA^{-1}D$.

It is possible to implement exactly (i.e., ideally with zero error) a
normalised matrix block of $A^{\prime}$ using at most $4$ calls to
${\mathcal{P}}_{y}$. First, we consider the unitary $\mathrm{Had}^{\prime}$
that prepares the state
$\left|{\boldsymbol{1}^{\prime}}\right\rangle=\left|{(1,1,\ldots,1,0)^{T}}\right\rangle$,
that is
$\mathrm{Had}\left|{0}\right\rangle=\left|{\boldsymbol{1}^{\prime}}\right\rangle$.
Then, the matrix

$\displaystyle{\mathcal{U}}_{A}:=\left(\begin{array}[]{cc}\\!\mathrm{Had}\\!&0\\\
0&\\!\mathrm{Had}\\!\end{array}\right)\left(\begin{array}[]{cc}I-(1-\epsilon)\left|{0}\right\rangle\\!\\!\left\langle{0}\right|&-\sqrt{1-(1-\epsilon)^{2}}\left|{0}\right\rangle\\!\\!\left\langle{0}\right|\\\
\sqrt{1-(1-\epsilon)^{2}}\left|{0}\right\rangle\\!\\!\left\langle{0}\right|&I-(1-\epsilon)\left|{0}\right\rangle\\!\\!\left\langle{0}\right|\end{array}\right)\left(\begin{array}[]{cc}\\!\mathrm{Had}^{\dagger}\\!&0\\\
0&\\!\mathrm{Had}^{\dagger}\\!\end{array}\right)$ (A.32)

is a normalised matrix-block-encoding of $A=I-(1-\epsilon)K^{\prime}$. Note
that the matrix in the centre can be interpreted as the
$\left|{0}\right\rangle\\!\\!\left\langle{0}\right|$-controlled version of the
Pauli-$X$ rotation $e^{i\theta X}$ (with $\cos\theta=-(1-\epsilon)$) and is
thus efficiently implementable. The operations in Eq. (A.4) correspond to a
unitary quantum circuit that can be written as $D\oplus{\mathcal{U}}$, for
some unitary ${\mathcal{U}}$, and finally we obtain that
${\mathcal{U}}_{A^{\prime}}:=(D\oplus{\mathcal{U}})\,{\mathcal{U}}_{A}\,(D\oplus{\mathcal{U}}^{\dagger})$
is a matrix-block-encoding of $A^{\prime}$.

We now consider the linear system $A^{\prime}\textbf{x}=\textbf{u}$ and a
quantum algorithm that prepares the corresponding solution state

$\displaystyle\left|{A^{\prime-1}\,\textbf{u}}\right\rangle=\left|{DA^{-1}D\,\textbf{u}}\right\rangle=D\left|{A^{-1}\,\textbf{b}}\right\rangle.$
(A.33)

The state $D\left|{A^{-1}\,\textbf{b}}\right\rangle$ can be transformed into
$\left|{A^{-1}\,\textbf{b}}\right\rangle$ using the steps given in Eq. (A.4),
which only requires two extra accesses to ${\mathcal{P}}_{y}$. This state
allows to solve PromiseMajorityM with constant probability and thus the same
considerations made in the preceding proof yield the result
$Q[{\mathcal{U}}_{A}]\in\Omega\big{(}\min(\kappa,N)\big{)}$.

∎

###### Proof.

_Case 3._

We start proving again a lower bound on the query complexity
$Q[{\mathcal{U}}_{\textbf{b}}]$, as in Case 1., but for a PD-QLS where the
coefficient matrix $A$ is sparse; then, we use the method used in the proof of
Case 2. to convert it into a lower bound on $Q[{\mathcal{P}}_{A}]$.

Given $y\in\\{0,1\\}^{N}$ satisfying the PromiseMajorityM condition, we
introduce the known term vector $\textbf{b}\in\mathds{R}^{N+1}$

$\displaystyle\left\\{\begin{array}[]{l}b_{i}=(-1)^{y_{i}}\qquad\qquad\qquad\mathrm{for}\
i\in[N]\\\\[5.69054pt] b_{N+1}=\sqrt{N}\,c_{0}\end{array}\right.$ (A.36)

where $c_{0}$ is a positive constant that we will fix later, and we have
$\left|{\textbf{b}}\right\rangle=\textbf{b}\Big{/}\sqrt{N(1+c_{0}^{2})}$.

Next, we define $B^{\prime}\in\mathds{R}^{(N+1)\times(N+1)}$ as

$\displaystyle B^{\prime}=B\oplus 0$ (A.37)

where $B\in\mathds{R}^{N\times N}$ is a symmetric sparse matrix, having $d$
non-zero entries in each row and column which are all equal to $\frac{1}{d}$,
for some constant $d$. Since each row and column of $B$ sums to one, $B$ can
be interpreted as the adjacency matrix of a Markov chain on a $d$-sparse
graph. We require that the graph corresponding to $B$ is not bipartite and
that the spectral gap of $B$ is large (i.e., the Markov chain is ergodic and
rapidly mixing). These properties guarantee that $B^{t}$ quickly converges to
$K=\frac{1}{N}\boldsymbol{1}_{N}\boldsymbol{1}_{N}^{\,T}$ for
$t\rightarrow\infty$. Since $B$ is symmetric, its spectrum is real and,
because of the Perron-Froebenius theorem [52], the spectrum is contained in
the interval $[-1,+1]$ and includes an eigenvalue $\lambda=1$ with
multiplicity one; the fact that $B^{t}$ converges to
$\frac{1}{N}\boldsymbol{1}_{N}\boldsymbol{1}_{N}^{\,T}$ implies that $-1$
cannot be an eigenvalue of $B$. We then define the spectral gap $\delta(B)$ as
the positive parameter $\delta(B):=\min_{\lambda\neq 1}\\{1-|\lambda|\\}$,
where the minimum is taken over the eigenvalues of $B$.

There are families of so-called _expander graphs_ such that both the sparsity
and the spectral gap are constant, see [53, Chapter 21]. We assume that $B$
belongs to one of these expander families and thus, in particular, there is a
(known) positive constant $c_{1}$ such that

$\displaystyle\frac{1}{\delta(B)}\leq c_{1}$ (A.38)

for all sizes $N\in\mathds{N}$. Moreover, $\boldsymbol{1}_{N}$ is the unique
$+1$ eigenvector and hence, from the spectral decomposition of $B$, we can
write

$\displaystyle
B=\frac{1}{N}\boldsymbol{1}_{N}\boldsymbol{1}_{N}^{\,T}+\sum_{\lambda\neq
1}\lambda\,\textbf{v}_{\lambda}\textbf{v}_{\lambda}^{\,{\dagger}}=K+R\,.$
(A.39)

Here $K=\frac{1}{N}\boldsymbol{1}_{N}\boldsymbol{1}_{N}^{\,T}$,
$\textbf{v}_{\lambda}$ are normalised eigenvectors of $B$, and thus
$R\in\mathds{R}^{N\times N}$ is a matrix of rank $N-1$ with
$\left|\left|{R}\right|\right|=1-\delta(B)$ and $KR=RK=0$. Then we define:

$\displaystyle A:=I-(1-\epsilon)B^{\prime}$ (A.40)

for some $\epsilon>0$ that we will define shortly. The matrix $A$ can be used
as a coefficient matrix in a PD-QLS solver since it is positive definite and
with norm one. Moreover, $A$ is $(d+1)$-sparse (where $d$ is the sparsity of
$B$, which is constant) and the condition number of $A$ is exactly
$\kappa(A)=1/\epsilon$.

Next we have

$\displaystyle A^{-1}$
$\displaystyle=I+\sum_{t=1}^{\infty}\big{(}(1-\epsilon)B^{\prime}\big{)}^{t}\equiv
I+\mathcal{B}$ (A.41)

and then defining
$\mathcal{K}:=\sum_{t=1}^{\infty}\big{(}(1-\epsilon)K^{\prime}\big{)}^{t}=\frac{1-\epsilon}{\epsilon}K^{\prime}$
we get

$\displaystyle\left|\left|{\mathcal{B}-\mathcal{K}}\right|\right|$
$\displaystyle\leq\sum_{t=1}^{\infty}\left|\left|{\big{[}(1-\epsilon)B\big{]}^{t}-\big{[}(1-\epsilon)K\big{]}^{t}}\right|\right|$
(A.42)
$\displaystyle=\sum_{t=1}^{\infty}(1-\epsilon)^{t}\left|\left|{(K+R)^{t}-K^{t}}\right|\right|$
(A.43)
$\displaystyle=\sum_{t=1}^{\infty}(1-\epsilon)^{t}\left|\left|{(K+R^{t})-K}\right|\right|$
(A.44) $\displaystyle\leq\sum_{t=1}^{\infty}(1-\delta(B))^{t}$ (A.45)
$\displaystyle\leq\frac{1}{\delta(B)}\leq c_{1}$ (A.46)

where we have used $\left|\left|{R}\right|\right|=1-\delta(B)$, $K^{2}=K$ and
$KR=RK=0$.

Applying $A^{-1}=I+\mathcal{B}$ to b we thus obtain:

$\displaystyle A^{-1}\textbf{b}$
$\displaystyle=\big{[}(I+\mathcal{B}-\mathcal{K})+\mathcal{K}\big{]}\textbf{b}$
(A.47)
$\displaystyle=(I+\mathcal{B}-\mathcal{K})\textbf{b}^{\prime}+\begin{cases}b_{N+1}\textbf{e}_{N+1}+\frac{1-\epsilon}{\epsilon}\frac{M}{N}\,\boldsymbol{1}_{N}^{\prime}&\
\mathrm{if~{}}y\mathrm{~{}has~{}a~{}majority~{}of~{}0}\\\
b_{N+1}\textbf{e}_{N+1}-\frac{1-\epsilon}{\epsilon}\frac{M}{N}\,\boldsymbol{1}_{N}^{\prime}&\
\mathrm{if~{}}y\mathrm{~{}has~{}a~{}majority~{}of~{}1}\end{cases}$ (A.48)
$\displaystyle=(I+\mathcal{B}-\mathcal{K})\textbf{b}^{\prime}+\begin{cases}b_{N+1}\textbf{e}_{N+1}+c_{0}\boldsymbol{1}_{N}^{\prime}&\qquad\;\mathrm{if~{}}y\mathrm{~{}has~{}a~{}majority~{}of~{}0}\\\
b_{N+1}\textbf{e}_{N+1}-c_{0}\boldsymbol{1}_{N}^{\prime}&\qquad\;\mathrm{if~{}}y\mathrm{~{}has~{}a~{}majority~{}of~{}1}\end{cases}$
(A.49)

where $\textbf{b}^{\prime}$ is a vector equal to the first $N$ entries of b
(and $b^{\prime}_{N+1}=0$), $\boldsymbol{1}_{N}^{\prime}:=(1,\ldots,1,0)^{T}$
while $\textbf{e}_{N+1}$ is the vector having a one in position $N+1$;
moreover, we choose $\epsilon$ such that
$\frac{1-\epsilon}{\epsilon}\frac{M}{N}=c_{0}$, where $c_{0}$ was introduced
in the definition of b in Eq. (A.36). Then, fixing the constant $c_{0}$ as
$c_{0}=100c_{1}$ and using the triangle inequality, we obtain the following
upper bound

$\displaystyle\mathcal{\sqrt{N}}=\left|\left|{A^{-1}\textbf{b}}\right|\right|$
$\displaystyle\leq\Big{(}\overbrace{c_{0}^{2}N}^{|b_{N+1}|^{2}}+\\!\\!\\!\overbrace{c_{0}^{2}N}^{\left|\left|{c_{0}\boldsymbol{1}_{N}^{\prime}}\right|\right|^{2}}\Big{)}^{1/2}+\overbrace{\sqrt{N}\,[1+c_{1}]}^{\geq\,\left|\left|{(I+\mathcal{B}-\mathcal{K})\textbf{b}^{\prime}}\right|\right|}$
(A.50)
$\displaystyle\leq\left(\sqrt{c_{0}^{2}+c_{0}^{2}}+2\,c_{1}\right)\sqrt{N}\ =\
\left(100\sqrt{2}+2\right)c_{1}\,\sqrt{N}$ (A.51) $\displaystyle\leq
144\,c_{1}\,\sqrt{N},$ (A.52)

where we have assumed $c_{1}\geq 1$. We also have the lower bound

$\displaystyle\sqrt{\mathcal{N}}=\left|\left|{A^{-1}\textbf{b}}\right|\right|$
$\displaystyle\geq\Big{(}\overbrace{c_{0}^{2}N}^{|b_{N+1}|^{2}}+\\!\\!\\!\overbrace{c_{0}^{2}N}^{\left|\left|{c_{0}\boldsymbol{1}_{N}^{\prime}}\right|\right|^{2}}\Big{)}^{1/2}-\overbrace{\sqrt{N}\,[1+c_{1}]}^{\geq\,\left|\left|{(I+\mathcal{B}-\mathcal{K})\textbf{b}^{\prime}}\right|\right|}$
(A.53) $\displaystyle\geq\left(100\sqrt{2}-2\right)c_{1}\,\sqrt{N}$ (A.54)
$\displaystyle\geq 139\,c_{1}\,\sqrt{N}.$ (A.55)

Then we have, using $f=\textsc{PromiseMajority}_{M}(y)$,

$\displaystyle\left|{A^{-1}\textbf{b}}\right\rangle=\frac{A^{-1}\textbf{b}}{\left|\left|{A^{-1}\textbf{b}}\right|\right|}$
$\displaystyle=\frac{c_{0}\,\sqrt{N}}{\sqrt{{\mathcal{N}}}}\left|{N+1}\right\rangle\,+\,(-1)^{f}\frac{c_{0}\,\sqrt{N}}{\sqrt{{\mathcal{N}}}}\sum_{i=1}^{N}\frac{1}{\sqrt{N}}\left|{i}\right\rangle+\left|{\psi}\right\rangle$
(A.56)
$\displaystyle=\frac{100}{144}\left|{N+1}\right\rangle\,+\,(-1)^{f}\frac{100}{144}\sum_{i=1}^{N}\frac{1}{\sqrt{N}}\left|{i}\right\rangle+\left|{\psi^{\prime}}\right\rangle,$
(A.57)

where
$\left|{\psi}\right\rangle:=(\mathcal{B}-\mathcal{K}+I)\textbf{b}^{\prime}/\mathcal{\sqrt{N}}$
is a sub-normalised perturbation vector with
$\big{|}\big{|}\left|{\psi}\right\rangle\big{|}\big{|}\leq\frac{2}{139}$,
while subtracting line (A.57) from line (A.56) one obtains

$\displaystyle\Big{|}\Big{|}\,\left|{\psi^{\prime}}\right\rangle-\left|{\psi}\right\rangle\,\Big{|}\Big{|}\leq\sqrt{2}\left(\frac{100}{139}-\frac{100}{144}\right)\leq
0.04$ (A.58)

and then we have:

$\displaystyle\Big{|}\Big{|}\,\left|{\psi^{\prime}}\right\rangle\,\Big{|}\Big{|}\leq
0.04+\frac{2}{139}\leq 0.06\,.$ (A.59)

The cases $f=0$ and $f=1$ in Eq. (A.57) can be distinguished with constant
advantage using the swap test with the state $\left|{``+"}\right\rangle$
defined in Eq. (A.17), since we have:

$\displaystyle~{}~{}\;\langle{``+"}|{A^{-1}\textbf{b}}\rangle\,=\,\frac{1}{\sqrt{2}}\frac{100}{144}\,+\,(-1)^{f}\frac{1}{\sqrt{2}}\frac{100}{144}\,+\,\langle{``+"}|{\,\psi^{\prime}\,}\rangle$
(A.60) $\displaystyle\Longleftrightarrow\quad$
$\displaystyle\begin{cases}\big{|}\,\langle{``+"}|{A^{-1}\textbf{b}}\rangle\,\big{|}\,\geq\,0.92&\mathrm{if~{}}y\mathrm{~{}has~{}a~{}majority~{}of~{}0}\\\
\big{|}\,\langle{``+"}|{A^{-1}\textbf{b}}\rangle\,\big{|}\,\leq\,0.06&\mathrm{if~{}}y\mathrm{~{}has~{}a~{}majority~{}of~{}1}\,.\end{cases}$
(A.61)

Next, we proceed as in the proof of Case 2. and define an equivalent PD-QLS
where the vector $y$ is encoded in the entries of the coefficient matrix. We
thus define the vector $\textbf{u}\in\mathds{R}^{N+1}$ and the diagonal matrix
$D\in\mathds{R}^{(N+1)\times(N+1)}$

$\displaystyle\left\\{\begin{array}[]{ll}u_{i}=1&\mathrm{for}\
i\in[N]\\\\[5.69054pt] u_{N+1}=\sqrt{N}\,c_{0}&\end{array}\right.\hskip
28.45274pt\left\\{\begin{array}[]{ll}D_{i,i}=(-1)^{y_{i}}&\mathrm{for}\
i\in[N]\\\\[5.69054pt] D_{N+1,N+1}=1&\end{array}\right.$ (A.66)

and thus the identity $\textbf{b}=D\textbf{u}$ holds. We then introduce
$A^{\prime}$, given by

$\displaystyle A^{\prime}:=DAD$ (A.67)

where $A$ is as in Eq. (A.40). The matrix $D$ is unitary and self-inverse,
hence $\kappa(A^{\prime})=\kappa(A)$, and
$A^{\prime-1}=D^{-1}A^{-1}D^{-1}=DA^{-1}D$.

Notice that the position of the non-zero entries of $A^{\prime}$ are the same
as in $A$ and thus independent from $y$, while the sign of a non-zero entry
$A^{\prime}_{i,j}=(-1)^{y_{i}+y_{j}}$ can be obtained querying
${\mathcal{P}}_{y}$ once with input $\left|{i}\right\rangle$ and once with
input $\left|{j}\right\rangle$. These queries can be performed in quantum
superposition and thus two accesses to ${\mathcal{P}}_{y}$ are sufficient to
implement a quantum sparse-matrix-access ${\mathcal{P}}_{A}^{\prime}$.
Therefore it is possible to prepare, with the same
${\mathcal{P}}_{y}$-complexity as discussed previously, the state

$\displaystyle\left|{A^{\prime-1}\,\textbf{u}}\right\rangle=\left|{DA^{-1}D\,\textbf{u}}\right\rangle=D\left|{A^{-1}\,\textbf{b}}\right\rangle\,.$
(A.68)

Finally, we can obtain $\left|{A^{-1}\,\textbf{b}}\right\rangle$ using two
extra accesses to ${\mathcal{P}}_{y}$ by applying the transformations given in
Eq. (A.4). This proves that a quantum algorithm that solves the PD-QLS having
access to ${\mathcal{P}}_{A}^{\prime}$ necessarily has a query complexity
$Q[{\mathcal{P}}_{A^{\prime}}]\in\Omega\big{(}\min(\kappa,N)\big{)}$.

∎

##  Appendix B Scaling of the normalisation factor of the matrix-block-
encoding

In this Appendix, we consider the polynomial

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

as was defined in Eq. (28) and where we have

$\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)}\;.$
(B.2)

We prove that the normalisation factor
$K:=2\,\max_{x\in[-1,+1]}P_{2\ell-1,\kappa}(x)$ satisfies
$K\in\Theta(\kappa)$, provided that $\ell\geq c\sqrt{\kappa}$ for some
constant $c$ that we will determine later.

Notice that by construction $P_{2\ell-1,\kappa}(1)=0$ and the polynomial is
positive for $x\in[-1,+1)$. To study the properties of the local maxima of
$P_{2\ell-1,\kappa}(x)$ in the interval $[-1,+1]$ we compute the derivative of
$P(x)$ using the property $\frac{\partial}{\partial
x}{\mathcal{T}}_{\ell}(x)=\ell\,{\mathcal{U}}_{\ell-1}(x)$, where
${\mathcal{U}}_{\ell}(x)\in\mathds{R}_{\ell}[x]$ is a Chebyshev polynomial of
the second kind. We have:

$\displaystyle\frac{\partial P_{2\ell-1,\kappa}(x)}{\partial
x}=\frac{\left[1-\hat{{\mathcal{T}}}_{\ell,\kappa}(x)\right]^{2}}{(1-x)^{2}}\;-\;2\ell\,\frac{1-\hat{{\mathcal{T}}}_{\ell,\kappa}(x)}{1-x}\frac{{\mathcal{U}}_{\ell-1}\left(\frac{x+\frac{1}{2\kappa}}{1-\frac{1}{2\kappa}}\right)}{\left(1-\frac{1}{2\kappa}\right){\mathcal{T}}_{\ell}\left(\frac{1+\frac{1}{2\kappa}}{1-\frac{1}{2\kappa}}\right)}\;.$
(B.3)

We set the derivative equal to $0$ and simplify the expression assuming
$1-x\neq 0$ and $1-\hat{{\mathcal{T}}}_{\ell,\kappa}(x)\neq 0$:

$\displaystyle\left(1-\frac{1}{2\kappa}\right)\left[{\mathcal{T}}_{\ell}\\!\left(\frac{1+\frac{1}{2\kappa}}{1-\frac{1}{2\kappa}}\right)-{\mathcal{T}}_{\ell}\\!\left(\frac{x+\frac{1}{2\kappa}}{1-\frac{1}{2\kappa}}\right)\right]\;-\;2\ell\,(1-x)\,{\mathcal{U}}_{\ell-1}\\!\left(\frac{x+\frac{1}{2\kappa}}{1-\frac{1}{2\kappa}}\right)=0$
(B.4)

Then we use the change of variables $y(x)$ and $\delta(\kappa)$ given in Eqs.
(24) and their inverses $\kappa(\delta)=\frac{1}{\delta}+\frac{1}{2}$,
$x(y)=\frac{y-\delta/2}{1+\delta/2}$ to rewrite the previous equation as

$\displaystyle\left(\frac{1}{1+\delta/2}\right)\big{[}{\mathcal{T}}_{\ell}(1+\delta)-{\mathcal{T}}_{\ell}(y)\big{]}\;-\;2\ell\,\left(\frac{1+\delta-y}{1+\delta/2}\right)\,{\mathcal{U}}_{\ell-1}(y)=0$
(B.5)

which is equivalent to:

$\displaystyle\boxed{{\mathcal{T}}_{\ell}(1+\delta)-{\mathcal{T}}_{\ell}(y)\;=\;2\ell\,(1+\delta-y)\,{\mathcal{U}}_{\ell-1}(y)}$
(B.6)

for $x\in[-1,+1]$ or, equivalently, $y\in[-1,1+\delta]$.

The polynomial $P_{2\ell-1,\kappa}(x)$ is by construction non-negative on the
domain $x\in[-1,+1]$ and its derivative in $x=1-\frac{1}{\kappa}$
(corresponding to $y=1$) is positive, provided that
$\ell\in\Omega(1/\sqrt{\delta})$. In fact, the derivative of
$P_{2\ell-1,\kappa}(x)$ is positive if and only if the left hand side of Eq.
(B.6) is larger than the right hand side; using
${\mathcal{T}}_{\ell}(1+\delta)\geq\frac{1}{2}e^{\ell\sqrt{\delta}}$ for
$0\leq\delta\leq 3-2\sqrt{2}$, ${\mathcal{T}}_{\ell}(1)=1$,
${\mathcal{U}}_{\ell-1}(1)=\ell$ we obtain from (B.6) the inequality
$\frac{1}{2}e^{\ell\sqrt{\delta}}-1\overset{!}{>}2\,\ell^{2}\delta$, which is
satisfied for $\ell\geq 4.36/\sqrt{\delta}$. Since we have, moreover,
$P_{2\ell-1,\kappa}(1)=0$, the function $P_{2\ell-1,\kappa}(x)$ does not have
any local maximum in $\left[\,-1,1-\frac{1}{\kappa},\right]$ and must have one
or more local maxima $x_{*}\in\left(\,1\\!-\\!\frac{1}{\kappa},+1\,\right]$;
equivalently, Eq. (B.6) must have at least one solution
$y_{*}\in(1,1+\delta]$.

Any local maximum $y_{*}=y(x^{*})$ satisfies:

$\displaystyle{\mathcal{T}}_{\ell}(y_{*})={\mathcal{T}}_{\ell}(1+\delta)-2\ell\,(1+\delta-
y_{*})\,{\mathcal{U}}_{\ell-1}(y_{*})$ (B.7)

and substituting ${\mathcal{T}}_{\ell}(y_{*})$ in the definition (B.1) gives

$\displaystyle P_{2\ell-1,\kappa}(x_{*})$
$\displaystyle=\frac{1+\delta/2}{1+\delta-
y_{*}}\left[1-\frac{{\mathcal{T}}_{\ell}(y_{*})}{{\mathcal{T}}_{\ell}(1+\delta)}\right]^{2}$
(B.8) $\displaystyle=4\ell^{2}(1+\delta/2)(1+\delta-
y_{*})\,\frac{{\mathcal{U}}_{\ell-1}(y_{*})^{2}}{{\mathcal{T}}_{\ell}(1+\delta)^{2}}\;.$
(B.9)

From Eq. (B.6) we directly have

$\displaystyle{\mathcal{U}}_{\ell-1}(y_{*})\leq\frac{{\mathcal{T}}_{\ell}(1+\delta)}{2\ell(1+\delta-
y_{*})}$ (B.10)

and inserting this inequality in Eq. (B.9) we have that any local maximum
$x_{*}$ satisfies

$\displaystyle P_{2\ell-1,\kappa}(x_{*})$
$\displaystyle\leq\frac{1+\delta/2}{1+\delta-y_{*}}$ (B.11)
$\displaystyle\leq\frac{3/2}{1+\delta-y_{*}}\;.$ (B.12)

In summary, it will be sufficient to show $1+\delta-y_{*}\in\Omega(\delta)$ to
prove that the normalisation constant satisfies
$K=2\,\max_{\\{x_{*}\\}}|P_{2\ell-1,\kappa}(x_{*})|\in{\mathcal{O}}(\kappa)$,
where the maximisation is over the set of (potentially multiple) local maxima
$x_{*}\in\left[\,1\\!-\\!\frac{1}{\kappa},+1\,\right]$.

To this end, we rewrite Eq. (B.6) as:

$\displaystyle\frac{{\mathcal{T}}_{\ell}(1+\delta)-{\mathcal{T}}_{\ell}(y_{*})}{(1+\delta)-y_{*}}\;=\;2\,\frac{\partial{\mathcal{T}}_{\ell}}{\partial
y}(y_{*})\;.$ (B.13)

Since both the first and the second derivative of ${\mathcal{T}}_{\ell}(y)$
are positive for $y\geq 1$ (i.e., the derivative of ${\mathcal{T}}_{\ell}(y)$
is monotonically increasing), this equation can be satisfied only if141414In
this appendix we employ the notation $\overset{!}{\leq}$ to mark inequalities
that we still need to prove and with $\leq$ an inequality that has already
been proven.

$\displaystyle 2\frac{\partial{\mathcal{T}}_{\ell}}{\partial
y}(y_{*})\overset{!}{\leq}\frac{\partial{\mathcal{T}}_{\ell}}{\partial
y}(1+\delta)$ (B.14)

or equivalently, defining $1+\delta_{*}:=y_{*}$

$\displaystyle
2\,{\mathcal{U}}_{\ell-1}(1+\delta_{*})\,\overset{!}{\leq}\,{\mathcal{U}}_{\ell-1}(1+\delta)\;.$
(B.15)

A Chebyshev polynomial of the second kind can be written as

$\displaystyle{\mathcal{U}}_{\ell-1}(y)=\frac{(y+\sqrt{y^{2}-1})^{\ell}-(y+\sqrt{y^{2}-1})^{-\ell}}{2\sqrt{y^{2}-1}}$
(B.16)

hence we have

$\displaystyle{\mathcal{U}}_{\ell-1}(1+\delta_{*})$
$\displaystyle\leq\frac{(1+\delta_{*}+\sqrt{2\delta_{*}+\delta_{*}^{2}})^{\ell}}{2\sqrt{2\delta_{*}+\delta_{*}^{2}}}$
(B.17)
$\displaystyle\leq\frac{(1+1.1\sqrt{2\,\delta_{*}})^{\ell}}{2\sqrt{2\delta_{*}}}$
(B.18)

for sufficiently small $\delta_{*}$ (the constant $1.1$ is somewhat
arbitrary), while we have

$\displaystyle{\mathcal{U}}_{\ell-1}(1+\delta)$
$\displaystyle=\frac{(1+\delta+\sqrt{2\delta+\delta^{2}})^{\ell}\big{[}1-(1+\delta+\sqrt{2\delta+\delta^{2}})^{-2\ell}\big{]}}{2\sqrt{2\delta+\delta^{2}}}$
(B.19) $\displaystyle\geq\frac{(1+\sqrt{2\delta})^{\ell}\times
0.8}{2.4\sqrt{2\delta}}$ (B.20)
$\displaystyle=\frac{2}{3}\frac{(1+\sqrt{2\delta})^{\ell}}{2\sqrt{2\delta}}$
(B.21)
$\displaystyle\geq\frac{2}{3}\frac{(1+\sqrt{2\delta})^{\ell}}{2\sqrt{3\delta_{*}}}$
(B.22)

which holds for $\ell\geq 1/\sqrt{2\delta}$ and in the last step we have
assumed $\delta_{*}\geq\frac{2}{3}\delta$ (notice that in the opposite case
$\delta_{*}<\frac{2}{3}\delta$ we would already have $1+\delta-
y_{*}=\delta-\delta_{*}>\frac{1}{3}\delta$). Thus, the inequality in (B.15) is
implied by

$\displaystyle
2\,\frac{(1+1.1\sqrt{2\,\delta_{*}})^{\ell}}{2\sqrt{2\delta_{*}}}\;$
$\displaystyle\overset{!}{\leq}\;\frac{2}{3}\frac{(1+\sqrt{2\delta})^{\ell}}{2\sqrt{3\delta_{*}}}$
(B.23) $\displaystyle\Longleftrightarrow\hskip
14.22636pt(1+1.1\sqrt{2\,\delta_{*}})^{\ell}\;$
$\displaystyle\overset{!}{\leq}\;\frac{\sqrt{2}}{3\sqrt{3}}\,(1+\sqrt{2\delta})^{\ell}$
(B.24) $\displaystyle\Longleftrightarrow\hskip
19.91692pt1+1.1\sqrt{2\,\delta_{*}}\ \ \;$
$\displaystyle\overset{!}{\leq}\;e^{-\frac{1}{\ell}\log(3\sqrt{3/2})}\,(1+\sqrt{2\delta})\;.$
(B.25)

From the inequality $e^{-x}\geq 1-x$ for $x\geq 0$ we see that the inequality
above is implied by

$\displaystyle 1.1\sqrt{2\,\delta_{*}}\;$
$\displaystyle\overset{!}{\leq}\;\big{(}1-1.31/\ell\big{)}\,(1+\sqrt{2\delta})-1$
(B.26) $\displaystyle=\;\sqrt{2\delta}-\frac{1}{\ell}1.31\,(1+\sqrt{2\delta})$
(B.27)

with $1.31\geq\log(3\sqrt{3/2})$. We now make the assumption that
$\frac{1}{\ell}1.31\,(1+\sqrt{2\delta})\leq\frac{1}{10}\sqrt{2\delta}$, which
is implied by $\ell\geq 13.1+9.27/\sqrt{\delta}$ and is compatible with the
requirement $\ell\in\Omega(\delta^{-1/2})$. Thus it is sufficient to impose:

$\displaystyle 1.1\sqrt{2\,\delta_{*}}\;$
$\displaystyle\overset{!}{\leq}\;\frac{9}{10}\sqrt{2\delta}$ (B.28)
$\displaystyle\Longleftrightarrow\qquad\qquad\delta_{*}\;$
$\displaystyle\overset{!}{\leq}\;\left(\frac{9}{11}\right)^{2}\delta\;.$
(B.29)

We also remark that $(9/11)^{2}\geq 2/3$, so that this last inequality is
compatible with the assumption $\delta_{*}\geq\frac{2}{3}\delta$ we made
earlier.

Plugging the bound in Eq. (B.29) into (B.11), together with the definition of
$\delta(\kappa)$, results in the explicit bound $K\leq 6.05\,\kappa$, i.e.
$K\in{\mathcal{O}}(\kappa)$, provided that $\ell\geq
13.1+9.27\sqrt{\kappa-1/2}$, i.e. $\ell\in\Omega(\sqrt{\kappa})$.

##  Appendix C VTAA optimization improving the runtime values of the first
algorithm

In this Appendix we use VTAA to speed-up the asymptotic runtime of the
algorithm based on polynomial approximations of $1/(1-x)$. Preliminarily, we
introduce a Lemma showing that this VTAA algorithm, as summarised in
Proposition 13, does provide an asymptotic query complexity improvement
(compared to general QLS solvers) for most values of $\Gamma_{A,\textbf{b}}$,
see the right plot in Figure 3.

######  Lemma 20.

Defining
$\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|}$
we have $\Gamma_{A,\textbf{b}}\in[1,\sqrt{\kappa}\,]$, under the usual
assumption that the spectrum of $A$ is contained in $[1/\kappa,1]$.

###### Proof.

We expand $\left|{\textbf{b}}\right\rangle$ in a basis of eigenvectors of $A$,
that is
$\left|{\textbf{b}}\right\rangle=\sum_{\lambda}\beta_{\lambda}\left|{\lambda}\right\rangle$
with $A\left|{\lambda}\right\rangle=\lambda\left|{\lambda}\right\rangle$,
$\langle{\lambda}|{\lambda^{\prime}}\rangle=\delta_{\lambda,\lambda^{\prime}}$,
and $\sum_{\lambda}|\beta_{\lambda}|^{2}=1$. Then, we can write

$\displaystyle\frac{\left|\left|{A^{-1/2}\left|{\textbf{b}}\right\rangle}\right|\right|^{2}}{\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|^{2}}=\frac{\sum_{\lambda}|\beta_{\lambda}|^{2}f_{\lambda}}{\sum_{\lambda}|\beta_{\lambda}|^{2}g_{\lambda}}\
$ $\displaystyle\in\
\big{[}\min_{\lambda}\\{f_{\lambda}/g_{\lambda}\\}\,,\,\max_{\lambda}\\{f_{\lambda}/g_{\lambda}\\}\big{]}\
\subseteq\ \big{[}\kappa^{-1},1\big{]}$ (C.1)

for $f_{\lambda}=\lambda^{-1}$, $g_{\lambda}=\lambda^{-2}$ and thus
$f_{\lambda}/g_{\lambda}=\lambda\in[\kappa^{-1},1]$. We take the square root
of both extrema and multiply by $\sqrt{\kappa}$ to conclude.

∎

###  C.1 Variable-time amplitude amplification review

In this Section we review the general VTAA method, mainly following Ref. [19,
Section 3].

######  Definition 21 (Variable-stopping-time quantum algorithm).

A _variable-stopping-time quantum algorithm_
${\mathcal{A}}={\mathcal{A}}_{m}\cdot\ldots\cdot{\mathcal{A}}_{1}\cdot{\mathcal{A}}_{0}$
is given by the application of $m+1$ sub-algorithms ${\mathcal{A}}_{j}$ in
sequence, acting on the Hilbert space
$\mathcal{H}=\mathcal{H}_{C}\otimes\mathcal{H}_{F}\otimes\mathcal{H}_{S}$
where $S$ is a “system register” of arbitrary size, $F$ is a “flag qubit” that
heralds success and $\mathcal{H}_{C}=\bigotimes_{j=0}^{m}\mathcal{H}_{C_{j}}$
is a “clock register” containing $m+1$ qubits $C_{0},C_{1},\ldots,C_{m}$.
$\mathcal{H}$ is initially prepared in the all-zero state
$\left|{0}\right\rangle_{\mathrm{all}}=\left|{0,0,0}\right\rangle_{C,F,S}$.
Each ${\mathcal{A}}_{j}$ acts on
$\mathcal{H}_{C_{j}}\otimes\mathcal{H}_{F}\otimes\mathcal{H}_{S}$ and
$\left|{1}\right\rangle_{C_{j}}$ indicates that the algorithm stops after the
application of ${\mathcal{A}}_{j}$; i.e., each ${\mathcal{A}}_{j}$ is a
controlled algorithm that acts if and only if the previous $j$ qubits in the
clock register are in the state $\left|{0}\right\rangle^{\otimes
j}\in\bigotimes_{i=0}^{j-1}\mathcal{H}_{C_{i}}$. We assume that all branches
of the computation end by step $m$. The successful branches of the algorithm
are those where the flag is in the state $\left|{1}\right\rangle_{F}$ and thus
we define

$\displaystyle
p_{\mathrm{succ}}:=\big{|}\big{|}\,\Pi_{\checkmark}{\mathcal{A}}\left|{0}\right\rangle_{\mathrm{all}}\big{|}\big{|}^{2}\qquad\mathrm{and}\qquad\left|{\psi_{\mathrm{succ}}}\right\rangle:=\frac{\Pi_{\checkmark}{\mathcal{A}}\left|{0}\right\rangle_{\mathrm{all}}}{\big{|}\big{|}\,\Pi_{\checkmark}{\mathcal{A}}\left|{0}\right\rangle_{\mathrm{all}}\big{|}\big{|}}$
(C.2)

where
$\Pi_{\checkmark}:=I_{C}\otimes\left|{1}\right\rangle\\!\\!\left\langle{1}\right|_{F}\otimes
I_{S}$.

By construction, a variable-stopping-time quantum algorithm produces a quantum
state of the form
${\mathcal{A}}\left|{0}\right\rangle_{\mathrm{all}}=\sum_{j=0}^{m}\alpha_{j}\left|{1_{j}}\right\rangle_{C}\left|{\Psi_{j}}\right\rangle_{F,S}$,
where

$\displaystyle\left|{1_{j}}\right\rangle_{C}:=\left|{0}\right\rangle^{\otimes
m-j}\left|{1}\right\rangle\left|{0}\right\rangle^{\otimes j}$ (C.3)

is a state having the qubit $C_{j}$ in $\left|{1}\right\rangle$ and the other
clock qubits in $\left|{0}\right\rangle$, while
$\left|{\Psi_{j}}\right\rangle_{F,S}$ are normalised quantum states in
$\mathcal{H}_{F}\otimes\mathcal{H}_{S}$ with
$\sum_{j=0}^{m}|\alpha_{j}|^{2}=1$.

######  Definition 22 (Stopping times).

We introduce a sequence of _stopping times_ $t_{\mathrm{min}}\equiv
t_{0}<t_{1}<t_{2}<\ldots<t_{m}\equiv t_{\mathrm{max}}$, where each
$t_{j}\in\mathds{N}$ is the complexity of the sub-algorithm
${\mathcal{A}}_{\leq j}:={\mathcal{A}}_{j}\cdot\ldots\cdot{\mathcal{A}}_{0}$.
The probability $p_{j}$ of stopping at time $t_{j}$ (i.e., after the execution
of algorithm ${\mathcal{A}}_{\leq j}$) and the $\ell_{2}$-_average stopping
time_ are defined as

$\displaystyle p_{j}:=\big{|}\big{|}\,\Pi_{C_{j}}{\mathcal{A}}_{\leq
j}\left|{0}\right\rangle_{\mathrm{all}}\big{|}\big{|}^{2}\qquad\quad
t_{\mathrm{avg}}:=\sqrt{\sum_{j=0}^{m}p_{j}\,t_{j}^{2}}$ (C.4)

with
$\Pi_{C_{j}}:=\left|{1_{j}}\right\rangle\\!\\!\left\langle{1_{j}}\right|_{C}\otimes
I_{F}\otimes I_{S}$. We define $\Pi_{\mathrm{stop}\leq
j}:=\sum_{i=0}^{j}\Pi_{C_{i}}$ and

$\displaystyle\Pi_{\mathrm{bad}}^{j}$ $\displaystyle:=\Pi_{\mathrm{stop}\leq
j}\cdot\big{(}I_{C}\otimes\left|{0}\right\rangle\\!\\!\left\langle{0}\right|_{F}\otimes
I_{S}\big{)}$ (C.5) $\displaystyle\Pi_{\mathrm{mg}}^{j}$
$\displaystyle:=I-\Pi_{\mathrm{bad}}^{j}$ (C.6)

which project onto “bad” and “maybe good” subspaces after the application of
${\mathcal{A}}_{\leq j}$. The “maybe good” subspace at step $j+1$ is contained
in the “maybe good” subspace at step $j$ and, by construction,
$\Pi_{\mathrm{mg}}^{m}{\mathcal{A}}\left|{0}\right\rangle_{\mathrm{all}}=\Pi_{\checkmark}{\mathcal{A}}\left|{0}\right\rangle_{\mathrm{all}}$.

In the definition above the $t_{j}$ can represent any complexity measure
(e.g., query complexity, gate complexity, circuit depth). In this Appendix we
only compute the query complexity $Q$ of the sub-algorithm
${\mathcal{A}}_{\leq j}$, but we remark that all the algorithms we consider
here are gate-efficient, i.e., the gate complexity is in
${\mathcal{O}}\big{(}Q\,\mathrm{poly}(\log Q,\log N)\big{)}$.

######  Definition 23 (Variable-Time Amplitude Amplification).

Given a variable-stopping-time algorithm
${\mathcal{A}}={\mathcal{A}}_{m}\cdot\ldots\cdot{\mathcal{A}}_{0}$ as in
Definition 21 and given a sequence of $m+1$ non-negative integers
$(k_{0},k_{1},\ldots,k_{m})$, we recursively define a variable-time
amplification
${\mathcal{A}}^{\prime}={\mathcal{A}}_{m}^{\prime}\cdot\ldots\cdot{\mathcal{A}}_{0}^{\prime}$
as follows. Setting ${\mathcal{A}}_{-1}:=I$, each ${\mathcal{A}}_{j}^{\prime}$
implements a standard $k_{j}$-step amplitude amplification that uses
${\mathcal{A}}_{j}{\mathcal{A}}_{j-1}^{\prime}$ and its inverse a total of
$2k_{j}+1$ times, where the input state is
$|{\psi_{\mathrm{in}}^{j}}\rangle={\mathcal{A}}_{j}{\mathcal{A}}_{j-1}^{\prime}\left|{0}\right\rangle_{\mathrm{all}}$
and the target state is
$\Pi_{\mathrm{mg}}^{j}|{\psi_{\mathrm{in}}^{j}}\rangle$. That is,
${\mathcal{A}}_{j}^{\prime}$ begins preparing the input state

$\displaystyle~{}|{\psi_{\mathrm{in}}^{j}}\rangle=\
{\mathcal{A}}_{j}{\mathcal{A}}_{j-1}^{\prime}\left|{0}\right\rangle_{\mathrm{all}}\,=\,\sin(\theta_{j})\,|{\psi_{\mathrm{mg}}^{j}}\rangle+\cos(\theta_{j})\,|{\psi_{\mathrm{bad}}^{j}}\rangle$
(C.7) $\displaystyle\mathrm{with}~{}$
$\displaystyle\begin{cases}|{\psi_{\mathrm{mg}}^{j}}\rangle\,\propto\;\Pi_{\mathrm{mg}}^{j}|{\psi_{\mathrm{in}}^{j}}\rangle\\\
|{\psi_{\mathrm{bad}}^{j}}\rangle\propto\;\Pi_{\mathrm{bad}}^{j}|{\psi_{\mathrm{in}}^{j}}\rangle\end{cases}\qquad\theta_{j}:=\arcsin\\!\big{(}\left|\left|{\Pi_{\mathrm{mg}}^{j}{\mathcal{A}}_{j}{\mathcal{A}}_{j-1}^{\prime}\left|{0}\right\rangle_{\mathrm{all}}}\right|\right|\big{)}\;\in\;[0,\pi/2]$
(C.8)

and then uses $k_{j}$ accesses to the reflections
$R_{\mathrm{out}}^{j}:=I-2\,\Pi_{\mathrm{mg}}^{j}$ and
$R_{\mathrm{in}}^{j}:=2\,|{\psi_{\mathrm{in}}^{j}}\rangle\langle{\psi_{\mathrm{in}}^{j}}|-I$,
where $R_{\mathrm{in}}^{j}$ is implemented using
${\mathcal{A}}_{j}{\mathcal{A}}_{j-1}^{\prime}$ and
$({\mathcal{A}}_{j}{\mathcal{A}}_{j-1}^{\prime})^{\dagger}$ once, to obtain
the output state

$\displaystyle|{\psi_{\mathrm{out}}^{j}}\rangle=\ $
$\displaystyle{\mathcal{A}}_{j}^{\prime}\left|{0}\right\rangle_{\mathrm{all}}\,=\,\sin\\!\big{[}(2k_{j}+1)\theta_{j}\big{]}\,|{\psi_{\mathrm{mg}}^{j}}\rangle+\cos\\!\big{[}(2k_{j}+1)\theta_{j}\big{]}\,|{\psi_{\mathrm{bad}}^{j}}\rangle\,.$
(C.9)

Note that if we set $k_{0}=\ldots=k_{m}=0$ we have
${\mathcal{A}}^{\prime}\equiv{\mathcal{A}}$. Also, by the recursive structure
of VTAA the first sub-algorithm ${\mathcal{A}}_{1}$ is used a total of
$\prod_{j=0}^{m}(2k_{j}+1)$ times in ${\mathcal{A}}^{\prime}$, which grows
exponentially in $m$ if $k_{j}\geq 1$. Nonetheless, VTAA can provide a speed-
up when the amplification parameters $(k_{0},k_{1},\ldots,k_{m})$ are chosen
appropriately. More precisely, the following result can be derived from Ref.
[19, Lemma 22].

######  Proposition 24 (Result of VTAA).

Using the notation of the previous definitions, let ${\mathcal{A}}^{\prime}$
be a variable-time amplification such that each ${\mathcal{A}}^{\prime}_{j}$
uses $k_{j}$ steps of amplitude amplification, where

$\displaystyle\frac{\pi}{8\,\theta_{j}}-\frac{1}{2}\,\leq\,k_{j}\,\leq\,\frac{\pi}{4\,\theta_{j}}-\frac{1}{2}$
(C.10)

Then, with the definitions in Eq. (C.2), ${\mathcal{A}}^{\prime}$ outputs the
state

$\displaystyle{\mathcal{A}}^{\prime}\left|{0}\right\rangle_{\mathrm{all}}\,=\,\sqrt{p_{\mathrm{succ}}^{\prime}}\,|{\psi_{\mathrm{succ}}}\rangle_{C,F,S}\,+\,\sqrt{1-p_{\mathrm{succ}}^{\prime}}\,|{\psi^{\perp}}\rangle_{C,F,S}$
(C.11)

where success is heralded by $\left|{1}\right\rangle_{F}$, the success
probability satisfies $p_{\mathrm{succ}}^{\prime}\in\Theta(1)$, and the global
query complexity is

$\displaystyle
Q^{\prime}\in{\mathcal{O}}\\!\left(t_{\mathrm{max}}\sqrt{m}+\frac{t_{\mathrm{avg}}}{\sqrt{p_{\mathrm{succ}}}}\sqrt{m\,\log(t_{\mathrm{max}}/t_{\mathrm{min}})}\right).$
(C.12)

We can compare Eq. (C.12) with standard amplitude amplification, which has a
query complexity
$Q\in{\mathcal{O}}(t_{\mathrm{max}}/\sqrt{p_{\mathrm{succ}}})$. In our
algorithm, $m$ will scale logarithmically in $t_{\mathrm{max}}$, so the total
runtime is
$\widetilde{{\mathcal{O}}}(t_{\mathrm{max}}+t_{\mathrm{avg}}/\sqrt{p_{\mathrm{succ}}})$.
Hence, VTAA can provide a speed-up when the average runtime is much shorter
than the maximum runtime.

We remark that a sequence of good values $(k_{0},\ldots,k_{m})$ such that
conditions in Eq. (C.10) are satisfied can be obtained efficiently by means of
an iterative classical-quantum pre-processing algorithm. Specifically, suppose
that some values $(k_{0},\ldots,k_{j})$ satisfying (C.10) have been already
found; then, one can compile the corresponding VTAA algorithm
${\mathcal{A}}_{j}^{\prime}$ and use phase estimation on the state
$|{\psi_{\mathrm{in}}^{j+1}}\rangle={\mathcal{A}}_{j+1}{\mathcal{A}}_{j}^{\prime}\left|{0}\right\rangle_{\mathrm{all}}$
to obtain an estimate of $\sin(\theta_{j+1})$ up to constant multiplicative
precision; this allows to find a value $k_{j+1}$ which satisfies (C.10) with
probability $1-p_{\mathrm{fail}}$ with a cost
${\mathcal{O}}\big{(}\frac{1}{\theta_{j+1}}\log(1/p_{\mathrm{fail}})\big{)}$
[30], which is asymptotically equal to the query complexity of
${\mathcal{A}}_{j+1}^{\prime}$, apart from a multiplicative
$\log(1/p_{\mathrm{fail}})$ overhead. Once $k_{j+1}$ has been precomputed one
can directly compile ${\mathcal{A}}_{j+1}^{\prime}$, without the need of
performing phase estimation “online”. The total failure probability is
$p_{\mathrm{fail}}^{\mathrm{tot}}\leq m\,p_{\mathrm{fail}}$ and the query cost
of the hybrid classical-quantum pre-processing algorithm has only a
${\mathcal{O}}\big{(}\log(1/p_{\mathrm{fail}})\big{)}$ multiplicative overhead
compared to the expression in Eq. (C.12).

###  C.2 Efficient implementation of windowing polynomials

In this Section we introduce the “windowing functions” that we use to replace
the Gapped Phase Estimation (GPE) subroutine, introduced in Ref. [5]. Using
GPE gives an additive ${\mathcal{O}}(\kappa)$ runtime overhead, which is too
costly for our purposes.

######  Lemma 25 (Efficient windowing polynomials).

Given $\epsilon,\delta\in(0,1/2]$, there exists an even polynomial
$W_{\epsilon,\delta}(x)=W_{\epsilon,\delta}(-x)$ of degree
$\ell\in{\mathcal{O}}\Big{(}\frac{1}{\sqrt{\delta}}\,\mathrm{polylog}(\delta^{-1},\epsilon^{-1})\Big{)}$,
where
$\mathrm{polylog}(\delta^{-1},\epsilon^{-1})\equiv\log^{1/4}(\epsilon^{-1})\big{[}\log(\epsilon^{-1})+\log(\delta^{-1})\big{]}$,
satisfying the inequalities

$\displaystyle
W_{\epsilon,\delta}(x)\in\begin{cases}[1-\epsilon,1]&\mathrm{if}~{}x\in[0,1-2\delta]\\\
[-1,+1]&\mathrm{if}~{}x\in(1-2\delta,1-\delta)\\\
[-\epsilon,+\epsilon]&\mathrm{if}~{}x\in[1-\delta,1]\;.\end{cases}$ (C.13)

The windowing function $W_{\epsilon,\delta}(x)$ can be computed efficiently
with classical algorithms.

A family of polynomials satisfying the inequalities in (C.13) is already given
in Ref. [25, Lemma 29], but they have a degree
$\ell\in\widetilde{{\mathcal{O}}}(\delta^{-1})$. We seek to achieve the same
result with a quadratically smaller degree,
$\ell\in\widetilde{{\mathcal{O}}}(\delta^{-1/2})$. Also, it is crucial that
the polynomial approximation has even parity: using QSP one can implement
matrix polynomials when the polynomial $P$ satisfies $|P(x)|\leq 1$ for
$|x|\leq 1$ and $P$ has definite parity, while for a polynomial $P^{\prime}$
without definite parity the more restrictive constraint $|P^{\prime}(x)|\leq
1/2$ has to be satisfied (see Theorem 8).

###### Proof.

The proof proceeds in two steps: first, we find an analytic function $f(x)$
that satisfies the inequalities in (C.13) within error $\epsilon/2$, and then
show that the Chebyschev expansion converges very quickly to it, so that
choosing a degree
$\ell\in\widetilde{{\mathcal{O}}}\big{(}\delta^{-1/2}\big{)}$ is sufficient to
be within $\epsilon/2$-distance from $f(x)$ for all $x$.

_First part:_ We introduce the normal distribution cumulative function

$\displaystyle\Phi(x):=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{x}e^{-\frac{t^{2}}{2}}\,dt\qquad
x\in\mathds{R}$ (C.14)

normalised so that $0\leq\Phi(x)\leq 1$. We then define

$\displaystyle\mathcal{W}_{\sigma,\delta}(x):=\Phi\\!\left(\frac{x+1-1.5\,\delta}{\sigma}\right)\Phi\\!\left(\frac{-x+1-1.5\,\delta}{\sigma}\right)$
(C.15)

where $\sigma>0$ has to be chosen so that
$\mathcal{W}_{\sigma,\delta}(1-2\delta)\geq 1-\epsilon/2$ and
$\mathcal{W}_{\sigma,\delta}(1-\delta)\leq\epsilon/2$. Using
$\Phi(-x)=1-\Phi(x)$ and the monotonicity of $\Phi$, both inequalities are
implied by $\Phi\\!\left(-\frac{0.5\,\delta}{\sigma}\right)\leq\epsilon/4$.
Using the bound
$\Phi(-x)\leq\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{-x}\frac{-t}{x}e^{-t^{2}/2}\,dt=\frac{e^{-x^{2}/2}}{\sqrt{2\pi}x}\leq\frac{e^{-x^{2}/2}}{\sqrt{2\pi}}$
for $x\geq 1$, it is sufficient to choose
$\sigma\in\Theta\big{(}\delta/\sqrt{\log(\epsilon^{-1})}\big{)}$ to obtain the
desired result.

_Second part:_ We consider the Chebyschev series associated to
$\mathcal{W}_{\sigma,\delta}(x)$. Note that the function
$\mathcal{W}_{\sigma,\delta}(x)$ has even parity, hence also the associated
Chebyschev series has even parity [54]. Truncating the Chebyschev series at
degree $\ell$, we get a sequence of polynomials
$[S_{\ell}\mathcal{W}_{\sigma,\delta}](x)$ converging uniformly to
$\mathcal{W}_{\sigma,\delta}(x)$ for $\ell\rightarrow\infty$. More precisely,
we have the following result [54, Theorem 5.16].

######  Lemma 26.

Suppose $f(x)$ can be extended to an analytic function on the ellipse
$E_{r}\subseteq\mathds{C}$

$\displaystyle E_{r}=\Big{\\{}\ a\cos\theta+i\,b\sin\theta\ \Big{|}\
a=\tfrac{1}{2}(r+r^{-1}),\ b=\tfrac{1}{2}(r-r^{-1}),\ \theta\in[0,2\pi)\
\Big{\\}}$ (C.16)

for some $r>1$; then, the $\ell$-th degree truncation of the Chebyschev series
is a polynomial $[S_{\ell}f](x)$ that satisfies for all $x\in[-1,+1]$

$\displaystyle\Big{|}f(x)-[S_{\ell}f](x)\Big{|}\leq\frac{M}{r^{\ell}(r-1)}\qquad
with~{}M=\sup\big{\\{}|f(z)|:z\in E_{r}\big{\\}}.$ (C.17)

We apply this theorem to $f(x)=\mathcal{W}_{\sigma,\delta}(x)$, which is
analytic in the whole complex plane. The main technical hurdle is to upper
bound $M$. We choose $r=\frac{1+\sqrt{\sigma}}{1-\sqrt{\sigma}}\geq
1+2\sqrt{\sigma}$, so that we have $a=\frac{1+\sigma}{1-\sigma}$ and
$b=\frac{2\sqrt{\sigma}}{1-\sigma}$, where $a$ and $b$ are the semi-axis of
the ellipse $E_{r}$ along the real and imaginary axis, respectively. We have,
for $x,y\in\mathds{R}$:

$\displaystyle|\Phi(x+iy)|$
$\displaystyle=\tfrac{1}{\sqrt{2\pi}}\Big{|}\int_{-\infty}^{x}e^{-\frac{t^{2}}{2}}dt+\int_{x}^{x+iy}e^{-\frac{t^{2}}{2}}dt\Big{|}$
(C.18) $\displaystyle\leq
1+\tfrac{1}{\sqrt{2\pi}}\Big{|}\int_{0}^{iy}e^{-\frac{(x+i\tau)^{2}}{2}}d\tau\Big{|}$
(C.19) $\displaystyle\leq
1+\tfrac{1}{\sqrt{2\pi}}\,e^{-x^{2}/2}\int_{0}^{|y|}e^{\frac{\tau^{2}}{2}}d\tau$
(C.20) $\displaystyle\leq
1+\tfrac{1}{\sqrt{2\pi}}\,e^{-x^{2}/2}\,e^{y^{2}/2}\,|y|$ (C.21)
$\displaystyle\leq 1+\tfrac{1}{\sqrt{2\pi}}e^{(-x^{2}+y^{2})/2}$ (C.22)

where the last inequality holds for $|y|\leq 1$. We now upper bound

$\displaystyle
M_{\mathrm{half}}:=\sup\left\\{\left|\Phi\\!\left(\tfrac{z+1-1.5\,\delta)}{\sigma}\right)\right|:z\in
E_{r}\right\\}$ (C.23)

so that we will have $M=\sup_{E_{r}}|\mathcal{W}_{\sigma,\delta}(z)|\leq
M_{\mathrm{half}}^{2}$. From (C.22) we only need to maximize the expression
$1+\frac{1}{\sqrt{2\pi}}e^{(-x^{2}+y^{2})/2}$ for

$\displaystyle x=\frac{a}{\sigma}\cos\theta+\frac{1-1.5\delta}{\sigma}\qquad
y=\frac{b}{\sigma}\sin\theta\,.$ (C.24)

Equivalently, we can maximize $-x^{2}+y^{2}$. Substituting $\cos\theta=u$ and
$\sin^{2}\theta=1-u^{2}$ we get

$\displaystyle-x^{2}+y^{2}$
$\displaystyle=\frac{1}{\sigma^{2}}\big{[}-(au+1-1.5\delta)^{2}+b^{2}(1-u^{2})\big{]}$
(C.25)
$\displaystyle\leq\frac{b^{2}}{2\sigma^{2}}\left[1-\frac{(1-1.5\delta)^{2}}{a^{2}+b^{2}}\right].$
(C.26)

We now have $b\in{\mathcal{O}}(\sqrt{\sigma})$,
$\frac{1}{a^{2}+b^{2}}=\frac{1-2\sigma+\sigma^{2}}{1+6\sigma+\sigma^{2}}=1-{\mathcal{O}}(\sigma)$
and thus

$\displaystyle\mathrm{sup}\,\\{-x^{2}+y^{2}\\}\ $
$\displaystyle\in{\mathcal{O}}\Big{(}\frac{\sigma}{\sigma^{2}}\Big{)}\Big{(}1-[1-{\mathcal{O}}(\delta)][1-{\mathcal{O}}(\sigma)]\Big{)}$
(C.27)
$\displaystyle={\mathcal{O}}\big{(}\sigma^{-1}\big{)}{\mathcal{O}}\big{(}\delta+\sigma\big{)}={\mathcal{O}}\big{(}\sqrt{\log\epsilon^{-1}}\,\big{)}$
(C.28)

where we have used
$\delta\in\Theta\big{(}\sigma\sqrt{\log\epsilon^{-1}}\big{)}$. We therefore
have $M_{\mathrm{half}}\in
e^{{\mathcal{O}}(\sqrt{\log\epsilon^{-1}}\,)}\subset{\mathcal{O}}\big{(}\epsilon^{-1}\big{)}$
and thus also $M\leq M_{\mathrm{half}}^{2}\in{\mathcal{O}}(\epsilon^{-2}\,)$.
We thus set
$\ell\in\Theta\big{(}\frac{1}{\sqrt{\sigma}}\log(\epsilon^{-3}\sigma^{-1/2})\big{)}$
and compute

$\displaystyle\Big{|}\mathcal{W}_{\sigma,\delta}(x)-[S_{\ell}\mathcal{W}_{\sigma,\delta}](x)\Big{|}$
$\displaystyle\leq\frac{M}{r^{\ell+1}(r-1)}$ (C.29)
$\displaystyle\leq\frac{C\,\epsilon^{-2}}{(1+2\sqrt{\sigma})^{\frac{D}{\sqrt{\sigma}}\log(\epsilon^{-3}\sigma^{-1/2})}\,2\sqrt{\sigma}}$
(C.30)
$\displaystyle\leq\frac{C\,\epsilon^{-2}}{e^{\log(\epsilon^{-3}\sigma^{-1/2})}\,2\sqrt{\sigma}}=(\epsilon^{3}\,\sigma^{1/2})\;\frac{C\,\epsilon^{-2}}{2\sqrt{\sigma}}$
(C.31) $\displaystyle\leq\epsilon/5$ (C.32)

where $C$ and $D$ are positive constants. We recall that
$\sigma\in\Theta\big{(}\delta/\sqrt{\log(\epsilon^{-1})}\big{)}$ and thus

$\displaystyle\ell$
$\displaystyle\in\Theta\Big{(}\frac{1}{\sqrt{\sigma}}\big{[}3\log(\epsilon^{-1})+\frac{1}{2}\log(\sigma^{-1})\big{]}\Big{)}$
(C.33)
$\displaystyle=\Theta\Big{(}\frac{1}{\sqrt{\delta}}\log^{1/4}(\epsilon^{-1})\big{[}\log(\epsilon^{-1})+\log(\delta^{-1})\big{]}\Big{)}$
(C.34)

is sufficient to guarantee that $[S_{\ell}\mathcal{W}_{\sigma,\delta}](x)$ is
within $\epsilon/5$ distance from $\mathcal{W}_{\sigma,\delta}(x)$.

Finally, we renormalise $[S_{\ell}\mathcal{W}_{\sigma,\delta}](x)$ dividing it
by the maximum value attained for $x\in[-1,1]$; by construction, this value is
smaller than $1+\epsilon/5$, thus the maximum distance from
$\mathcal{W}_{\sigma,\delta}$ after normalisation is bounded by
$1-\frac{1-\epsilon/5}{1+\epsilon/5}<\epsilon/2$.

∎

###  C.3 Polynomial approximations on increasingly larger domains

Now we explain how the windowing functions can be used to select domain where
a polynomial of low degree can be a good approximation of the function
$1/(1-x)$, if it is applied to eigenvalues contained in such domain.

As usual, the spectrum of $A$ is contained in
${\mathcal{D}}_{A}=\left[\,\frac{1}{\kappa},\,1\right]$ and, correspondingly,
the spectrum of $B=I-\eta\,A$ is contained in
${\mathcal{D}}_{B}=\left[1-\eta,1-\frac{\eta}{\kappa}\,\right]\subseteq\left[0,1-\frac{\eta}{\kappa}\,\right]$.
We set $m:=\lceil\log_{2}\kappa\rceil+1$, $\delta_{j}:=\eta\,2^{-j}$ for
$j=\\{1,\ldots,m\\}$ and fix ${\tilde{\epsilon}}>0$, a parameter related to
the target precision $\varepsilon$. According to Eq. (27) and Eq. (C.13) we
can find polynomials $P_{{\tilde{\epsilon}},\delta_{j}}(x)$ and
$W_{{\tilde{\epsilon}},\delta_{j}}(x)$ such that

$\displaystyle\Big{|}P_{{\tilde{\epsilon}},\delta_{j}}(x)-\frac{1}{1-x}\Big{|}\leq{\tilde{\epsilon}}\qquad\forall
x\in[-1,1-\delta_{j}]$ (C.35) $\displaystyle
W_{{\tilde{\epsilon}},\delta_{j}}(x)\qquad\quad\mathrm{satisfies~{}Eq.}~{}\eqref{eq:W_inequalities}~{}\mathrm{for~{}each~{}}\delta_{j}$
(C.36)

with degrees in
${\mathcal{O}}\big{(}\delta^{-1/2}\log({\tilde{\epsilon}}^{-1}\delta^{-1})\big{)}$
and
${\mathcal{O}}\big{(}\delta^{-1/2}\log^{1/4}({\tilde{\epsilon}}^{-1})\log({\tilde{\epsilon}}^{-1}\delta^{-1})\big{)}$,
respectively. We also define a normalisation factor

$\displaystyle
K:=2\max_{j}\max_{x\in[-1,1]}|P_{{\tilde{\epsilon}},\delta_{j}}(x)|$ (C.37)

which can be set to coincide with the factor $K$ defined in Eq. (28) and
introduce the shorthands

$\displaystyle
P_{j}(x):=P_{{\tilde{\epsilon}},\delta_{j}}(x)/K\quad\mathrm{and}\quad
W_{j}(x):=W_{{\tilde{\epsilon}},\delta_{j}}(x).$ (C.38)

The windowing function $W_{j}(x)$ is used to select an interval
$[-1+2\delta_{j},1-2\delta_{j}]$ where $P_{j}(x)$ is a good approximation of
the inverse. For any eigenvalue $\lambda$ of $B$ we have
$P_{j}(\lambda)\approx\frac{1}{K}\frac{1}{1-\lambda}$ when $\lambda\leq
1-\delta_{j}$ and $W_{j}(\lambda)\approx 0$ when $\lambda\geq 1-\delta_{j}$.
More precisely, we have

$\displaystyle W_{j}(B)P_{j}(B)$
$\displaystyle=W_{j}(B)\frac{1}{K}\frac{1}{I-B}+\Delta_{j}^{\tilde{\epsilon}}$
(C.39)
$\displaystyle=W_{j}(B)\,\frac{A^{-1}}{\eta\,K}+\Delta_{j}^{\tilde{\epsilon}}$
(C.40)

where $\Delta_{j}^{\tilde{\epsilon}}$ are arbitrary matrices having operator
norm smaller or equal to ${\tilde{\epsilon}}$. Note that $W_{j}(B)$ and
$P_{j}(B)$ commute, since they are both polynomial functions of $B$. Moreover,
we can set without loss of generality $W_{m}(B)\equiv I$, since we already
have $P_{m}(B)\approx A^{-1}/(\eta\,K)$ on the entire domain of $A$. Finally,
we have $K\in\Theta(\kappa/\eta)$, as derived in the main text.

###  C.4 Variable-stopping-time PD-QLS solver: definition

We now proceed to reformulate our PD-QLS solver as a variable-stopping-time
algorithm. For notation convenience, from now on we often drop the dependency
of a polynomial from its variable and simply write $P$ instead of $P(x)$.

######  Algorithm 27 (Variable-stopping-time linear system solving).

We preliminarily define two variable-stopping-time algorithms ${\mathcal{A}}$
and ${\mathcal{B}}$. ${\mathcal{A}}$ is the core PD-QLS solving function and
${\mathcal{B}}$ is used to un-compute the clock registers at the end of the
process.

The algorithm
${\mathcal{A}}={\mathcal{A}}_{m}\cdot\ldots\cdot{\mathcal{A}}_{1}\cdot{\mathcal{A}}_{0}$,
with $m:=\lceil\log_{2}\kappa\rceil+1$, acts on the Hilbert space
$\mathcal{H}=\mathcal{H}_{C}\otimes\mathcal{H}_{F}\otimes\mathcal{H}_{S}$
(respectively clock register, flag qubit, system register with
$\dim\mathcal{H}_{S}=N$), plus ancillary registers if needed. We set
${\mathcal{A}}_{0}=I_{C}\otimes I_{F}\otimes{\mathcal{U}}_{\textbf{b}}$ (i.e.,
${\mathcal{A}}_{0}$ prepares $\left|{\textbf{b}}\right\rangle$ in register
$S$) while the other sub-algorithms ${\mathcal{A}}_{j}$ use oracular access
only to ${\mathcal{U}}_{B}$. Each ${\mathcal{A}}_{j}$, for $j\geq 1$, consists
of the following two steps:

  1. 1.

Conditioned on the qubits $C_{1},\ldots,C_{j-1}$ being in
$\left|{0}\right\rangle^{\otimes j-1}$, use QSP and a single Pauli $X$ gate to
implement the unitary

$\displaystyle{\mathcal{U}}_{j}^{\prime}:=\left(\begin{array}[]{cc}\sqrt{1-W_{j}^{2}}&-W_{j}\\\
W_{j}&\sqrt{1-W_{j}^{2}}\end{array}\right)\quad\mathrm{acting~{}on}\
\mathcal{H}_{C_{j}}\otimes\mathcal{H}_{S},$ (C.43)

with the window function $W_{j}(B)$ given in Eq. (C.38).

  2. 2.

Conditioned on $C_{j}$ being in $\left|{1}\right\rangle_{C_{j}}$, use QSP and
a single Pauli $X$ gate to implement the unitary

$\displaystyle{\mathcal{U}}_{j}^{\prime\prime}:=\left(\begin{array}[]{cc}\sqrt{I-P_{j}^{2}}&-P_{j}\\\
P_{j}&\sqrt{I-P_{j}^{2}}\end{array}\right)\quad\mathrm{acting~{}on}\
\mathcal{H}_{F}\otimes\mathcal{H}_{S},$ (C.46)

with the polynomial approximations $P_{j}(B)$ given in Eq. (C.38).

The algorithm
${\mathcal{B}}={\mathcal{B}}_{m}\cdot\ldots\cdot{\mathcal{B}}_{1}\cdot{\mathcal{B}}_{0}$
acts on the same Hilbert space as ${\mathcal{A}}$. The initialisation step is
skipped, i.e., we set ${\mathcal{B}}_{0}=I$, and $\forall j\geq 1$ we define
${\mathcal{B}}_{j}$ as the application of the unitary
${\mathcal{U}}_{j}^{\prime}$ given in Eq. (C.43) controlled on
$C_{1},\ldots,C_{j-1}$ being in $\left|{0}\right\rangle^{\otimes j-1}$. The
unitaries ${\mathcal{U}}_{j}^{\prime\prime}$ are not applied.

The complete PD-QLS solving algorithm is then defined as follows.

  1. 1.

In a initial pre-processing step (that needs to be run only once), a sequence
of integers $(k_{0},k_{1},\ldots,k_{m})$ that satisfy Eq. (C.10) is determined
using a phase estimation algorithm.

  2. 2.

The main part of the algorithm consists in running ${\mathcal{A}}^{\prime}$,
which is a VTAA version of ${\mathcal{A}}$ where each
${\mathcal{A}}_{j}^{\prime}$ implements a $k_{j}$-step amplification of
${\mathcal{A}}_{j}{\mathcal{A}}_{j-1}^{\prime}$.

  3. 3.

The final post-processing consists in applying ${\mathcal{B}}^{\dagger}$ to
the output of ${\mathcal{A}}^{\prime}$ in order to un-compute the clock
register. We then measure the flag qubit and when the result is
$\left|{1}\right\rangle_{F}$ (which happens with constant probability) we
output the $S$ register (which contains a state close to
$\left|{A^{-1}\textbf{b}}\right\rangle$).

###  C.5 Variable-stopping-time PD-QLS solver: correctness

We will now analyse the correctness of Algorithm 27 and we begin introducing a
Lemma.

######  Lemma 28.

The state after the application of ${\mathcal{A}}_{k}$ defined in Algorithm 27
is given by

$\displaystyle{\mathcal{A}}_{\leq k}\left|{0}\right\rangle_{\mathrm{all}}=$
$\displaystyle\sum_{j=1}^{k}\left|{1_{j}}\right\rangle_{C}\left(\left|{0}\right\rangle_{F}M_{j-1}W_{j}\sqrt{I-P_{j}^{2}}\left|{\textbf{b}}\right\rangle_{S}+\left|{1}\right\rangle_{F}M_{j-1}W_{j}P_{j}\left|{\textbf{b}}\right\rangle_{S}\right)$
(C.47) $\displaystyle+M_{k}\left|{0,0,\textbf{b}}\right\rangle_{C,F,S}$ (C.48)

with $M_{0}:=I$, $M_{j}:=\prod_{i=1}^{j}\sqrt{I-W_{i}^{2}}$ and, using
$W_{m}=I$, at the last step we have $M_{m}=0$. Moreover, the algorithm
${\mathcal{B}}$ defined in Algorithm 27 acts as

$\displaystyle{\mathcal{B}}\left|{0,f,\phi}\right\rangle_{C,F,S}=$
$\displaystyle\sum_{j=1}^{k}\left|{1_{j}}\right\rangle_{C}\left|{f}\right\rangle_{F}M_{j-1}W_{j}\left|{\phi}\right\rangle_{S}$
(C.49)

for any input state $\left|{\phi}\right\rangle_{S}$ and $f\in\\{0,1\\}$.

###### Proof.

We proceed by induction over $k$. For $k=0$ we have ${\mathcal{A}}_{\leq
0}\left|{0}\right\rangle_{\mathrm{all}}={\mathcal{A}}_{0}\left|{0}\right\rangle_{\mathrm{all}}=\left|{0,0,\textbf{b}}\right\rangle_{C,F,S}$.
Using that $P_{j},M_{j},W_{j}$ are functions of $B$ and thus commute, we have
at step $k+1$

$\displaystyle\eqref{eq:induction1}\
\overset{{\mathcal{U}}_{k+1}^{\prime}}{\mapsto}\ $
$\displaystyle\sum_{j=1}^{k}\left|{1_{j}}\right\rangle_{C}\left(\left|{0}\right\rangle_{F}M_{j-1}W_{j}\sqrt{I-P_{j}^{2}}\left|{\textbf{b}}\right\rangle_{S}+\left|{1}\right\rangle_{F}M_{j-1}W_{j}P_{j}\left|{\textbf{b}}\right\rangle_{S}\right)$
(C.50)
$\displaystyle+\left|{1_{k+1}}\right\rangle_{C}\left|{0}\right\rangle_{F}M_{k}W_{k+1}\left|{\textbf{b}}\right\rangle_{F}+\left|{0,0}\right\rangle_{C,F}M_{k}\sqrt{1-W_{k+1}^{2}}\left|{\textbf{b}}\right\rangle_{C,F,S}$
(C.51) $\displaystyle\eqref{eq:induction2}\
\overset{{\mathcal{U}}_{k+1}^{\prime\prime}}{\mapsto}\ $
$\displaystyle\sum_{j=1}^{k}\left|{1_{j}}\right\rangle_{C}\left(\left|{0}\right\rangle_{F}M_{j-1}W_{j}\sqrt{I-P_{j}^{2}}\left|{\textbf{b}}\right\rangle_{S}+\left|{1}\right\rangle_{F}M_{j-1}W_{j}P_{j}\left|{\textbf{b}}\right\rangle_{S}\right)$
(C.52)
$\displaystyle+\left|{1_{k+1}}\right\rangle_{C}\left(\left|{0}\right\rangle_{F}M_{k}W_{k+1}\sqrt{I-P_{k+1}^{2}}\left|{\textbf{b}}\right\rangle_{S}+\left|{1}\right\rangle_{F}M_{k}W_{k+1}P_{k+1}\left|{\textbf{b}}\right\rangle_{S}\right)$
(C.53) $\displaystyle+M_{k+1}\left|{0,0,\textbf{b}}\right\rangle_{C,F,S}.$
(C.54)

Equation (C.49) can be verified similarly. ∎

Next, we consider the output state of ${\mathcal{A}}^{\prime}$, the VTAA
version of ${\mathcal{A}}$, which features an amplification of the amplitude
of the $\left|{1}\right\rangle_{F}$ component, i.e., it is a state of of the
form

$\displaystyle{\mathcal{A}}^{\prime}\left|{0}\right\rangle_{\mathrm{all}}$
$\displaystyle=\sqrt{p_{\mathrm{succ}}^{\prime}}\,|{1}\rangle_{F}|{\psi_{\mathrm{succ}}}\rangle_{C,S}\,+\,\sqrt{1-p_{\mathrm{succ}}^{\prime}}\,|{0}\rangle_{F}|{\psi_{\mathrm{fail}}}\rangle_{C,S}$
(C.55)

where the success probability is constant,
$p_{\mathrm{succ}}^{\prime}\in\Theta(1)$, and where we have

$\displaystyle|{\psi_{\mathrm{succ}}}\rangle_{C,S}$
$\displaystyle=\frac{1}{\sqrt{{\mathcal{N}}}}\sum_{j=1}^{m}\left|{1_{j}}\right\rangle_{C}M_{j-1}W_{j}P_{j}\left|{\textbf{b}}\right\rangle_{S}$
(C.56) $\displaystyle{\mathcal{N}}$
$\displaystyle=\sum_{j=1}^{m}\big{|}\big{|}M_{j-1}W_{j}P_{j}\left|{\textbf{b}}\right\rangle\big{|}\big{|}^{2}\,.$
(C.57)

We now claim that, for an appropriate choice of the algorithm parameters, the
error in $\left|{\psi_{\mathrm{succ}}}\right\rangle_{C,S}$ is upper bounded by
${\mathcal{O}}(\varepsilon)$, and we will prove this claim in the next
section. Thus we have:

$\displaystyle|{\psi_{\mathrm{succ}}}\rangle_{C,S}$
$\displaystyle=\sum_{j=1}^{m}\left|{1_{j}}\right\rangle_{C}M_{j-1}W_{j}\left|{A^{-1}\textbf{b}}\right\rangle_{S}+{\mathcal{O}}(\varepsilon)$
(C.58)

where ${\mathcal{O}}(x)$ here denotes an arbitrary vector with 2-norm bounded
by ${\mathcal{O}}(x)$. Note that Eq. (C.49) for
$\left|{\phi}\right\rangle_{S}=\left|{A^{-1}\textbf{b}}\right\rangle_{S}$
implies that
$\sum_{j}\left|{1_{j}}\right\rangle_{C}M_{j-1}W_{j}\left|{A^{-1}\textbf{b}}\right\rangle_{S}$
is a normalised state.

The final step of the PD-QLS algorithm consists in applying
${\mathcal{B}}^{\dagger}$ to the state in Eq. (C.55). Using again Eq. (C.49)
we obtain:

$\displaystyle{\mathcal{B}}^{\dagger}{\mathcal{A}}^{\prime}\left|{0}\right\rangle_{\mathrm{all}}=\sqrt{p_{\mathrm{succ}}^{\prime}}\left|{0,1,A^{-1}\textbf{b}}\right\rangle_{C,F,S}+\sqrt{1-p_{\mathrm{succ}}^{\prime}}\left|{0,0,\psi^{\prime}}\right\rangle_{C,F,S}+{\mathcal{O}}(\varepsilon)\,.$
(C.59)

Measuring the flag in $\left|{1}\right\rangle_{F}$ then results with constant
success probability in a vector that has ${\mathcal{O}}({\tilde{\epsilon}})$
distance from the ideal output $\left|{A^{-1}\textbf{b}}\right\rangle$.

###  C.6 Variable-stopping-time PD-QLS solver: error bound

In order to derive Eq. (C.58) we use Eq. (C.40) and write

$\displaystyle|{\psi_{\mathrm{succ}}}\rangle_{C,S}$
$\displaystyle=\frac{1}{\sqrt{{\mathcal{N}}}}\sum_{j=1}^{m}\left|{1_{j}}\right\rangle_{C}M_{j-1}W_{j}\frac{A^{-1}}{\eta\,K}\left|{\textbf{b}}\right\rangle_{S}+\frac{1}{\sqrt{{\mathcal{N}}}}\sum_{j=1}^{m}\left|{1_{j}}\right\rangle_{C}\Delta_{j}^{\tilde{\epsilon}}\left|{\textbf{b}}\right\rangle_{S}$
(C.60)
$\displaystyle=\frac{\sqrt{\widetilde{{\mathcal{N}}}}}{\sqrt{{\mathcal{N}}}}\sum_{j=1}^{m}\left|{1_{j}}\right\rangle_{C}M_{j-1}W_{j}\left|{A^{-1}\textbf{b}}\right\rangle_{S}+{\mathcal{O}}({\tilde{\epsilon}}\sqrt{m/{\mathcal{N}}})$
(C.61) $\displaystyle=\hskip
22.76219pt\sum_{j=1}^{m}\left|{1_{j}}\right\rangle_{C}M_{j-1}W_{j}\left|{A^{-1}\textbf{b}}\right\rangle_{S}+{\mathcal{O}}({\tilde{\epsilon}}\sqrt{m/{\mathcal{N}}})$
(C.62)

and $\widetilde{{\mathcal{N}}}$ is defined below. To prove the last step,
recall that
$\sum_{j}\left|{1_{j}}\right\rangle_{C}M_{j-1}W_{j}\left|{A^{-1}\textbf{b}}\right\rangle_{S}$
is a normalised state, hence Eq. (C.61) implies
$\big{|}\sqrt{\widetilde{{\mathcal{N}}}/{\mathcal{N}}}-1\big{|}\in{\mathcal{O}}({\tilde{\epsilon}}\sqrt{m/{\mathcal{N}}})$.
Now we estimate ${\mathcal{N}}$ by using
$\big{|}\sqrt{\widetilde{{\mathcal{N}}}}-\sqrt{{\mathcal{N}}}\big{|}\in{\mathcal{O}}({\tilde{\epsilon}}\sqrt{m})$
and computing

$\displaystyle\widetilde{{\mathcal{N}}}$
$\displaystyle:=\sum_{j=1}^{m}\left|\left|{M_{j-1}W_{j}\frac{A^{-1}}{\eta\,K}\left|{\textbf{b}}\right\rangle}\right|\right|^{2}$
(C.63)
$\displaystyle=\frac{||A^{-1}\left|{\textbf{b}}\right\rangle||^{2}}{\eta^{2}K^{2}}\sum_{j=1}^{m}\left|\left|{M_{j-1}W_{j}\left|{A^{-1}\textbf{b}}\right\rangle}\right|\right|^{2}$
(C.64)
$\displaystyle=\frac{||A^{-1}\left|{\textbf{b}}\right\rangle||^{2}}{\eta^{2}K^{2}}.$
(C.65)

As proven in the main text, $K\in{\mathcal{O}}(\kappa/\eta)$. Thus we have
$\widetilde{{\mathcal{N}}}\in\Omega(1/\kappa^{2})$ and choosing

$\displaystyle{\tilde{\epsilon}}\in{\mathcal{O}}\\!\left(\frac{\varepsilon}{\kappa\,\sqrt{\log\kappa}}\right)$
(C.66)

we also have ${\mathcal{N}}\in\Omega(1/\kappa^{2})$ and
${\tilde{\epsilon}}\sqrt{m/{\mathcal{N}}}\in{\mathcal{O}}(\varepsilon)$. This
condition is sufficient to guarantee that the output state
$\left|{\psi_{\mathrm{succ}}}\right\rangle$ is within $\varepsilon$-distance
from the ideal output.

###  C.7 Variable-stopping-time PD-QLS solver: query complexity

We can now provide a upper bounds the query complexity of the VTAA algorithm
using Eq. (C.12):

$\displaystyle
Q^{\prime}\in{\mathcal{O}}\\!\left(t_{\mathrm{max}}\sqrt{m}+\frac{t_{\mathrm{avg}}}{\sqrt{p_{\mathrm{succ}}}}\sqrt{m\,\log(t_{\mathrm{max}}/t_{\mathrm{min}})}\right).$

We first compute the ${\mathcal{U}}_{\textbf{b}}$-complexity. Since the non-
amplified algorithm ${\mathcal{A}}$ accesses ${\mathcal{U}}_{\textbf{b}}$ only
in the first step, we have
$t_{\mathrm{min}}=t_{\mathrm{max}}=t_{\mathrm{avg}}=1$. Using
$m\,{\tilde{\epsilon}}^{2}\in{\mathcal{O}}\big{(}\varepsilon^{2}/\kappa^{2}\big{)}$
we get

$\displaystyle p_{\mathrm{succ}}$ $\displaystyle={\mathcal{N}}\ \in\
{\mathcal{O}}\\!\left(\widetilde{{\mathcal{N}}}+m\,{\tilde{\epsilon}}^{2}\right)={\mathcal{O}}\\!\left(\frac{||A^{-1}\left|{\textbf{b}}\right\rangle||^{2}}{\kappa^{2}}\right)$
(C.67)

and the using $m=\lceil\log_{2}\kappa\rceil+1$ we have

$\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).$
(C.68)

We now compute the ${\mathcal{U}}_{B}$-complexity. We have:

$\displaystyle t_{j}$ $\displaystyle=\sum_{i=1}^{j}[\deg(P_{i})+\deg(W_{i})]\
\in\
{\mathcal{O}}\Big{(}\delta_{j}^{-1/2}\log^{1/4}({\tilde{\epsilon}}^{-1})\log({\tilde{\epsilon}}^{-1}\delta_{j}^{-1})\Big{)}$
(C.69) $\displaystyle t_{\mathrm{max}}$
$\displaystyle=\sum_{i=1}^{m}[\deg(P_{i})+\deg(W_{i})]\ \in\
{\mathcal{O}}\Big{(}\sqrt{\kappa/\eta}\,\log^{1/4}({\tilde{\epsilon}}^{-1})\log({\tilde{\epsilon}}^{-1}\kappa/\eta)\Big{)}$
(C.70) $\displaystyle t_{\mathrm{min}}$ $\displaystyle=\hskip
14.22636pt\deg(P_{1})+\deg(W_{1})\ \ \in\
{\mathcal{O}}\Big{(}\sqrt{1/\eta}\,\log^{1/4}({\tilde{\epsilon}}^{-1})\log({\tilde{\epsilon}}^{-1}/\eta)\Big{)}.$
(C.71)

Writing
$\left|{\textbf{b}}\right\rangle\equiv\sum_{\lambda}\beta_{\lambda}\left|{\lambda}\right\rangle$,
where $\left|{\lambda}\right\rangle$ denotes an eigenvector of $B$ relative to
the eigenvalue $\lambda$ and $\sum_{\lambda}|\beta_{\lambda}|^{2}=1$, the
probability of stopping at time $t_{j}$ is

$\displaystyle p_{j}$
$\displaystyle=\big{|}\big{|}\,\Pi_{C_{j}}{\mathcal{A}}_{\leq
j}\left|{0}\right\rangle_{\mathrm{all}}\big{|}\big{|}^{2}=\big{|}\big{|}\,M_{j-1}W_{j}\left|{\textbf{b}}\right\rangle\big{|}\big{|}^{2}$
(C.72)
$\displaystyle=\sum_{\lambda}|M_{j-1}(\lambda)W_{j}(\lambda)|^{2}\,|\beta_{\lambda}|^{2}$
(C.73) $\displaystyle\in\
{\mathcal{O}}\Big{(}{\tilde{\epsilon}}+\sum\nolimits_{\lambda\in(1-4\delta_{j},1-\delta_{j}]}|\beta_{\lambda}|^{2}\Big{)}.$
(C.74)

Here, we have used $W_{j}(\lambda)\leq{\tilde{\epsilon}}$ for $\lambda\geq
1-\delta_{j}$, while for $\lambda<1-2\delta_{j-1}=1-4\delta_{j}$ we have

$\displaystyle
M_{j-1}(\lambda)\leq\sqrt{1-W_{j-1}^{2}(\lambda)}\leq\sqrt{1-(1-{\tilde{\epsilon}})^{2}}\leq\sqrt{2{\tilde{\varepsilon}}}\,,$
(C.75)

i.e., the expression $|M_{j-1}(\lambda)W_{j}(\lambda)|^{2}$ is non-negligible
only for $\lambda\in(1-4\delta_{j},1-\delta_{j}]$. Now we estimate the
$\ell_{2}$-average runtime
$t_{\mathrm{avg}}=\sqrt{\sum_{j=1}^{m}p_{j}t_{j}^{2}}$ with

$\displaystyle t_{\mathrm{avg}}^{2}$
$\displaystyle\in{\mathcal{O}}\\!\left(\sum_{j=1}^{m}\Big{[}{\tilde{\epsilon}}+\sum\nolimits_{\lambda\in(1-4\delta_{j},1-\delta_{j}]}|\beta_{\lambda}|^{2}\Big{]}\Big{[}\delta_{j}^{-1/2}\log^{1/4}({\tilde{\epsilon}}^{-1})\log({\tilde{\epsilon}}^{-1}\delta_{j}^{-1})\Big{]}^{2}\right)$
(C.76)
$\displaystyle\subseteq{\mathcal{O}}\\!\left({\tilde{\epsilon}}\,t_{\mathrm{max}}^{2}+\sum\nolimits_{\lambda}\frac{|\beta_{\lambda}|^{2}}{1-\lambda}\,\log^{1/2}({\tilde{\epsilon}}^{-1})\log^{2}({\tilde{\epsilon}}^{-1}\kappa/\eta)\right)$
(C.77)

where we have used that for any positive function $f(\lambda)$

$\displaystyle\sum_{j=1}^{m}\sum_{\lambda\in(1-4\delta_{j},1-\delta_{j}]}f(\lambda)\,\frac{1}{\delta_{j}}\
\in\ \Theta\bigg{(}\sum_{\lambda}f(\lambda)\,\frac{1}{1-\lambda}\bigg{)}.$
(C.78)

We then write
$\sum_{\lambda}\frac{|\beta_{\lambda}|^{2}}{1-\lambda}=\left|\left|{\sum_{\lambda}\frac{\beta_{\lambda}}{\sqrt{I-B}}\left|{\lambda}\right\rangle}\right|\right|^{2}=\left|\left|{(\eta\,A)^{-1/2}\left|{\textbf{b}}\right\rangle}\right|\right|^{2}$
and obtain

$\displaystyle t_{\mathrm{avg}}$
$\displaystyle\in{\mathcal{O}}\\!\left(\frac{1}{\sqrt{\eta}}\big{|}\big{|}A^{-1/2}\left|{\textbf{b}}\right\rangle\big{|}\big{|}\,\log^{1/4}({\tilde{\epsilon}}^{-1})\log({\tilde{\epsilon}}^{-1}\kappa/\eta)\right)$
(C.79)

where we used
$\left|\left|{A^{-1/2}\left|{\textbf{b}}\right\rangle}\right|\right|\geq 1$
and thus
$\sqrt{{\tilde{\epsilon}}}\,t_{\mathrm{max}}\in{\mathcal{O}}\left(\frac{1}{\sqrt{\eta}\log^{1/4}\kappa}\log^{1/4}({\tilde{\epsilon}}^{-1})\log({\tilde{\epsilon}}^{-1}\kappa/\eta)\right)$
is a sub-leading contribution to $t_{\mathrm{avg}}$. Then, using
$1/\sqrt{p_{\mathrm{succ}}}\in{\mathcal{O}}\big{(}\kappa/\left|\left|{A^{-1}\left|{\textbf{b}}\right\rangle}\right|\right|\big{)}$
and

$\displaystyle\log(t_{\mathrm{max}}/t_{\mathrm{min}})\ \in\
{\mathcal{O}}\big{(}\log(\kappa)+\log\log({\tilde{\epsilon}}^{-1}\kappa/\eta)\big{)}\
\subseteq\
{\mathcal{O}}\big{(}\log({\tilde{\epsilon}}^{-1}\kappa/\eta)\big{)}$ (C.80)

we finally obtain

$\displaystyle Q[{\mathcal{U}}_{B}]$
$\displaystyle\in{\mathcal{O}}\Bigg{(}\sqrt{\frac{\kappa}{\eta}}\,\underbrace{\sqrt{\kappa}\frac{\left|\left|{A^{-1/2}\left|{\textbf{b}}\right\rangle}\right|\right|}{\big{|}\big{|}A^{-1}\left|{\textbf{b}}\right\rangle\big{|}\big{|}}}_{:=\Gamma_{A,\textbf{b}}}\underbrace{\log^{1/4}({\tilde{\epsilon}}^{-1})\log({\tilde{\epsilon}}^{-1}\kappa/\eta)\sqrt{\log(\kappa)\log({\tilde{\epsilon}}^{-1}\kappa/\eta)}}_{\mathrm{polylog}(\kappa,{\tilde{\epsilon}}^{-1},\eta^{-1})}\Bigg{)}.$
(C.81)

##  Appendix D Proof that the Sum-QLS problem is BQP-hard

In this Appendix we prove that the Sum-QLS problem is $\mathsf{BQP}$-hard,
therefore no efficient classical algorithm can solve the problem (unless
$\textsf{BPP}=\textsf{BQP}$). The initial part of the proof is equal to the
reduction presented by HHL: it is possible to construct a QLS problem
$M\textbf{x}=\textbf{e}_{1}$ (where $\textbf{e}_{1}$ is a canonical vector
with a one in the first position) which is BQP-hard for a class of sparse
indefinite matrices $M$ that are easily constructible. This can be converted
to the equivalent QLS problem $M\textbf{x}=M^{\dagger}\textbf{e}_{1}$, where
the matrix $A=M^{\dagger}M$ is by construction PD and
$\left|{M^{\dagger}\textbf{e}_{1}}\right\rangle$ is easy to prepare. What
remains to be proven is that $A$ admits an explicit Sum-QLS structure, i.e.
that one can construct a decomposition as a sum of PD local Hamiltonian terms
and that the overlap with the support (bounded by the parameter $\gamma$)
scales polynomially.

Preliminarily, we introduce a couple of useful definitions: the specific
quantum circuit model (universal for quantum computation) that we aim to
“simulate” as a Sum-QLS, and the Feynman-Kitaev clock construction.

##### Quantum circuit model:

We describe an arbitrary quantum computation $\mathcal{C}$ on $n$ qubits as
the application of $T$ elementary gates $\\{U_{0},U_{1}\ldots,U_{T-1}\\}$ to
an initial state $\left|{0}\right\rangle^{\otimes n}$ and the output of the
computation is the quantum state $U_{T-1}\cdots
U_{1}U_{0}\left|{0}\right\rangle^{\otimes n}$. We assume that the elementary
gate set consists of gates acting either on one or two qubits, e.g. arbitrary
single-qubit rotations and control-nots. Hence, the $t$-th elementary gate is
described by a unitary $U_{t}\in\mathds{C}^{N\times N}$ with $N=2^{n}$ which
can be written as $U_{t}=u_{t}\otimes I_{\neg{\mathcal{S}}_{t}}$, where
$u_{t}$ acts on a subset ${\mathcal{S}}_{t}$ of the qubits and is either a
$2\times 2$ ($|{\mathcal{S}}_{t}|=1$) or a $4\times 4$
($|{\mathcal{S}}_{t}|=2)$ unitary matrix.

##### Feynman-Kitaev clock:

We introduce the following unitary matrix acting on
$\mathds{C}^{3T}\otimes\mathds{C}^{n}$, which is based on the _Feynman-Kitaev
clock_ construction (see e.g. [55] and references therein):

$\displaystyle U:=$
$\displaystyle\sum_{t=0}^{3T-1}\left|{t\\!+\\!1}\right\rangle\\!\left\langle{t}\right|_{c}\otimes
U_{t}$ (D.1) $\displaystyle=$ $\displaystyle\
\sum_{t=0}^{T-1}\,\left|{t\\!+\\!1}\right\rangle\\!\left\langle{t}\right|_{c}\otimes
U_{t}\,+\,\left|{T\\!+\\!t\\!+\\!1}\right\rangle\\!\left\langle{T\\!+\\!t}\right|_{c}\otimes
I\,+\,\left|{2T\\!+\\!t\\!+\\!1}\right\rangle\\!\left\langle{2T\\!+\\!t}\right|_{c}\otimes
U_{T-t-1}^{\dagger}$ (D.2)

where we have defined $U_{t}=I$ for $T\leq t\leq 2T-1$ and
$U_{t}=U^{\dagger}_{3T-t-1}$ for $2T\leq t\leq 3T-1$ and the sums in the clock
register (denoted by the subscript $c$) are taken modulo $3T$.

######  Proposition 29 (Sum-QLS is BQP-hard).

Let $\mathcal{C}$ be a given $n$-qubit $T$-gate quantum circuit. Then, there
is a Sum-QLS problem (as in Definition 14) with

  1. 1.

$n^{\prime}\in{\mathcal{O}}(n+\log T)$
[$N^{\prime}\in{\mathcal{O}}(2^{n^{\prime}})$ is the size of $A$]

  2. 2.

$\kappa\in{\mathcal{O}}(T^{2})$ [condition number]

  3. 3.

$J\in{\mathcal{O}}(T)$ [number of PD Hamiltonian terms]

  4. 4.

$s\in{\mathcal{O}}(\log T)$ [locality of the Hamiltonian terms]

  5. 5.

$d_{\textbf{b}}\in{\mathcal{O}}(1)$ [sparsity of the known-term vector b]

  6. 6.

$\gamma\in\Omega(T^{-2})$ [overlap with the support, Eq. (112)]

that is equivalent to $\mathcal{C}$, i.e., solving this Sum-QLS problem (up to
a small constant error $\varepsilon$) allows to obtain the output state of
$\mathcal{C}$ with constant probability and constant precision. According to
Eq. (113), the algorithm presented in Section 5 can solve this Sum-QLS problem
with a gate complexity in ${\mathcal{O}}\big{(}\mathrm{poly}(n,T)\big{)}$.
When $T\in{\mathcal{O}}(\mathrm{poly}\,n)$, the gate complexity of the Sum-QLS
solver is also in ${\mathcal{O}}(\mathrm{poly}\,n)$, i.e., Sum-QLS problem is
BQP-hard. Moreover, the Sum-QLSpoly, defined as the subclass of problems where
the six parameters listed above all scale polynomially in $n$, is BQP-
complete.

###### Proof.

Starting from the Feynman-Kitaev unitary $U$ encoding a quantum circuit
$\mathcal{C}$ as in Eq. (D.1), we introduce the matrix

$\displaystyle M:=I-U\,e^{-1/T}$ (D.3)

which can be written in block form (where each block has size $2^{n}\times
2^{n}$) as

$\displaystyle M=$
$\displaystyle\left(\begin{array}[]{ccccc}\ddots&~{}~{}\ddots{}{}{}{}&0&&\\\
\ddots&I&-U_{t}e^{-1/T}&0&\\\ &0&I&-U_{t+1}e^{-1/T}&\\\
&&0&I&~{}~{}~{}\ddots{}{}{}\\\ -U_{3T-1}e^{-1/T}&&&\ddots&\ddots\\\
\end{array}\right)$ (D.9)

and moreover we introduce:

$\displaystyle A$
$\displaystyle:=M^{\dagger}M=(1+e^{-2/T})I-e^{-1/T}(U+U^{\dagger})$ (D.10)

where $A$ is by construction Hermitian positive definite and the singular
values of $M$ are the square root of the eigenvalues of $A$. From the
eigenvalue inequality $-2\leq\lambda(U+U^{\dagger})\leq 2$ we obtain

$\displaystyle\begin{cases}\lambda_{\max}(A)\leq 4\\\
\lambda_{\min}(A)\geq\big{(}1-e^{-1/T}\big{)}^{2}\geq\frac{(1-e^{-1})^{2}}{T^{2}}\end{cases}\Longrightarrow\quad\kappa(A)\leq\frac{4}{(1-e^{-1})^{2}}\,T^{2}$
(D.11)

hence $\kappa(M)=\sqrt{\kappa(A)}\in{\mathcal{O}}(T)$. Using the identity
$U^{3T}=I$, the inverse of $M$ can be expanded as

$\displaystyle M^{-1}$
$\displaystyle=\sum_{t^{\prime}=0}^{\infty}U^{t^{\prime}}e^{-t^{\prime}/T}=\frac{e^{3}}{e^{3}-1}\sum_{t=0}^{3T-1}U^{t^{\prime}}e^{-t^{\prime}/T}$
(D.12)
$\displaystyle=\frac{e^{3}}{e^{3}-1}\sum_{t^{\prime}=0}^{3T-1}e^{-t^{\prime}/T}\sum_{t=0}^{3T-1}\left|{t\\!+\\!t^{\prime}}\right\rangle\\!\left\langle{t}\right|\otimes
U_{[t:t+t^{\prime}]}$ (D.13)

where the notation $U_{[a:b]}$ indicates the ordered product of unitary
operators from $a$ to $b-1$, $U_{[a:b]}:=U_{b-1}U_{b-2}\cdots U_{a}$, with
$U_{[t:t]}=I$.

As done already by HHL in Ref. [1] step, we consider the QLS problem

$\displaystyle M\textbf{x}=\textbf{e}_{1}$ (D.14)

where $\textbf{e}_{1}\in\mathds{C}^{3T2^{n}}$ is the unit vector with a one in
the first position (i.e., it is the vector representation of the state
$\left|{0}\right\rangle_{c}\left|{0}\right\rangle^{\otimes n}$). Using the
state
$\left|{\textbf{x}}\right\rangle=\left|{M^{-1}\textbf{e}_{1}}\right\rangle$
the output state of the quantum computation is obtained whenever upon
measurement of the clock register a time $t\in[T:2T-1]$ is returned, which
happens with probability $e^{-2}/(1+e^{-2}+e^{-4})\geq 0.11$. Therefore, any
quantum circuit $\mathcal{C}$ having $T$ gates can be restated as a QLS and
solved (using e.g. the HHL algorithm) with a gate complexity scaling as
$\mathrm{poly}(n,T)$. If the circuit $\mathcal{C}$ has
$T\in{\mathcal{O}}(\mathrm{poly}\,n)$, then also the QLS solver has gate
complexity in ${\mathcal{O}}(\mathrm{poly}\,n)$, showing that the QLS problem
is BQP-hard.

We now consider the PD linear system

$\displaystyle A\textbf{x}=\textbf{b}\qquad\mathrm{for}\ A=M^{\dagger}M\,,\
\;\textbf{b}$ $\displaystyle=M^{\dagger}\textbf{e}_{1}$ (D.15)

which is equivalent to the system in Eq. (D.14) and can be cast as a Sum-QLS
problem.

Notice, first, that b is a sparse vector, with $d_{\textbf{b}}\leq 3$: if the
elementary gate set consists of control-nots and single-qubit rotations, then
the first row of $M$ has at most $3$ non-zero entries.

The matrix $A$ has size $3T2^{n}\times 3T2^{n}$ and can be written as

$\displaystyle A=\sum_{t=0}^{3T-1}H_{(t)}\;.$ (D.16)

where each Hamiltonian terms $H_{(t)}$ is positive definite and only acts on
the clock register plus either one or two extra qubits (i.e., the same number
of qubits on which the gates from the elementary set act). We start writing
$A$ in block form

$\displaystyle A=\left(\begin{array}[]{ccccc}\ddots&\ddots&0&&\\\
\ddots&I(1+e^{-2/T})&-U_{t}e^{-1/T}&0&\\\
~{}~{}~{}0{}{}{}&-U_{t}^{\dagger}e^{-1/T}&I(1+e^{-2/T})&-U_{t+1}e^{-1/T}&~{}~{}~{}0{}{}{}\\\
&0&-U_{t+1}^{\dagger}e^{-1/T}&I(1+e^{-2/T})&\ddots\\\ &&0&\ddots&\ddots\\\
\end{array}\right)$ (D.22)

and then introduce, for each $t\in\\{0,\ldots,3T-1\\}$, the Hamiltonian term

$\displaystyle
H_{(t)}\;=\;I\,\delta\;+\;e^{-1/T}\left(\begin{array}[]{c|cc|c}~{}0{}&0&0&~{}0{}\\\
\hline\cr 0&I&-U_{t}&0\\\ 0&-U_{t}^{\dagger}&I&0\\\ \hline\cr 0&0&0&0\\\
\end{array}\right)\begin{array}[]{l}\hskip 1.0pt\\}~{}\mathrm{size}~{}t\cdot
2^{n}\\\\[6.25958pt] \Big{\\}}~{}\mathrm{size}~{}2\cdot 2^{n}\\\\[7.11317pt]
\hskip 1.0pt\\}~{}\mathrm{size}~{}(3T-t-2)\cdot 2^{n}\end{array}$ (D.30)

where we have defined

$\displaystyle\delta\;:=\;\frac{1}{3T}\big{(}1+e^{-2/T}-2e^{-1/T}\big{)}\;\geq\;\frac{(1-e^{-1})^{2}}{3T^{3}}\;\geq\;\frac{1}{7.51\,T^{3}}\;,$
(D.31)

and thus the decomposition in Eq. (D.22) holds. Note that the matrix
$\left(\begin{smallmatrix}I&-U_{t}\\\
-U_{t}^{\dagger}&I\end{smallmatrix}\right)$ has eigenvalues $\lambda=0$ and
$\lambda=2$, hence is positive semi-definite, while each $H_{(t)}$ has the
smallest eigenvalue equal to $\delta$. Moreover, each $H_{(t)}$ is a local
Hamiltonian. Specifically, we have $H_{(t)}=h_{(t)}\otimes
I_{\neg{\mathcal{S}}_{t}}$, where $h_{(t)}$ acts on a set ${\mathcal{S}}_{t}$
consisting of $s=\lceil\log_{2}3T\rceil+q$ qubits, where $q=1$ if $U_{t}$
corresponds to a single qubit gate and $q=2$ if it corresponds to a two-qubit
gate. Explicitly, we have

$\displaystyle
h_{(t)}\;=\;\left(\begin{array}[]{c|cc|c}I\,\delta&0&0&~{}0{}\\\ \hline\cr
0&I(e^{-1/T}\\!+\\!\delta)\\!&-u_{t}e^{-1/T}&0\\\
0&-u_{t}^{\dagger}e^{-1/T}&\\!I(e^{-1/T}\\!+\\!\delta)&0\\\ \hline\cr
0&0&0&I\,\delta\\\ \end{array}\right)\begin{array}[]{l}\hskip
1.0pt\\}~{}\mathrm{size}~{}t\cdot 2^{q}\\\\[6.25958pt]
\Big{\\}}~{}\mathrm{size}~{}2\cdot 2^{q}\\\\[7.11317pt] \hskip
1.0pt\\}~{}\mathrm{size}~{}(3T-t-2)\cdot 2^{q}\end{array}$ (D.39)

where each $u_{t}$ is a matrix specifying the $t$-th elementary quantum gate,
with $U_{t}=u_{t}\otimes I_{\neg{\mathcal{S}}_{t}}$.

We now need to estimate a lower bound $\gamma$ for the overlap parameter as in
Eq. (112). We have

$\displaystyle\left\langle{\textbf{b}}\right|A^{-1}\left|{\textbf{b}}\right\rangle$
$\displaystyle=\left\langle{M^{\dagger}\textbf{e}_{1}}\right|M^{-1}M^{-{\dagger}}\left|{M^{\dagger}\textbf{e}_{1}}\right\rangle$
(D.40)
$\displaystyle=\frac{\left\langle{\textbf{e}_{1}}\right|MM^{-1}M^{-{\dagger}}M^{\dagger}\left|{\textbf{e}_{1}}\right\rangle}{\left|\left|{M^{\dagger}\left|{\textbf{e}_{1}}\right\rangle}\right|\right|^{2}}$
(D.41)
$\displaystyle=\frac{1}{\left|\left|{M^{\dagger}\left|{\textbf{e}_{1}}\right\rangle}\right|\right|^{2}}=\frac{1}{1+e^{-2/T}}\geq
1/2$ (D.42)

while from (D.30) and (D.31) we obtain

$\displaystyle\left\langle{\textbf{b}}\right|H_{(j)}^{-1}\left|{\textbf{b}}\right\rangle\,\leq\,\delta^{-1}\,\leq\,7.51\,T^{3}\,.$
(D.43)

Using $J=3T$ we finally get

$\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}{(3T)^{2}}\,\frac{3T\times
7.51\,T^{3}}{1/2}\,\leq\,5.01\,T^{2}\,=:\,\frac{1}{\gamma}\,.$ (D.44)

In conclusion, given the sequence of gates $u_{t}$ and the set of qubits
${\mathcal{S}}_{t}$ to which they are applied in the quantum circuit
$\mathcal{C}$, we can efficiently compute the values and the positions of the
non-zero entries of $H_{(t)}$, as given in Eq. (D.30). This then explicitly
describes $A$ as a sum of PD local Hamiltonian terms, as required by the
definition of the Sum-QLS problem. Going through the derivation, we see that
the relevant parameters scale as stated in points (1) to (6) in the statement
of the Proposition. Using Eq. (113), it follows that the algorithm presented
in Section 5 solves this problem with a gate complexity in
${\mathcal{O}}\big{(}\mathrm{poly}(n,T)\big{)}={\mathcal{O}}(\mathrm{poly}\,n)$,
in the case where the original circuit is itself a polynomial time quantum
circuit (i.e. $T\in{\mathcal{O}}(\mathrm{poly}\,n)$). Finally, note that the
$\varepsilon$-error in precision of Sum-QLS solver is amplified at most by a
constant factor when post-selecting the clock register to show a time
$t\in[T:2T-1]$. Thus, it is possible to obtain the output of $\mathcal{C}$
with an error that is bounded by a constant, as required per the definition of
the BQP class. This finally proves that the Sum-QLS problem is BQP-hard; at
the same time, it proves that Sum-QLSpoly, the sub-class of problems having
polynomially scaling parameters, is solvable in quantum polynomial time and is
thus BQP-complete.

∎

## References

  * [1] A. W. Harrow, A. Hassidim, and S. Lloyd, _Quantum algorithm for linear systems of equations_ , Physical Review Letters 103, 150502 (2009) [arXiv:0811.3171].
  * [2] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, _Numerical Recipes: The Art of Scientific Computing, Third Edition_ , Cambridge University Press (2007).
  * [3] Y. Saad, _Iterative methods for sparse linear systems_ , SIAM (2003).
  * [4] A. Ambainis, _Variable time amplitude amplification and quantum algorithms for linear algebra problems_ , 29th Symposium on Theoretical Aspects of Computer Science 14, 636–647 (2012) [arXiv:1010.4458].
  * [5] A. M. Childs, R. Kothari, and R. D. Somma, _Quantum linear systems algorithm with exponentially improved dependence on precision_ , SIAM J. Comput. 46, 1920–1950 (2017) [arXiv:1511.02306].
  * [6] L. Wossnig, Z. Zhao, and A. Prakash, _Quantum linear system algorithm for dense matrices_ , Physical Review Letters 120, 050502 (2018) [arXiv:1704.06174].
  * [7] Y. Subaşı, R. D. Somma, and D. Orsucci. _Quantum algorithms for linear systems of equations inspired by adiabatic quantum computing_ , Physical Review Letters 122, 060504 (2019) [arXiv:1805.10549].
  * [8] Dong An and Lin Lin, _Quantum linear system solver based on time-optimal adiabatic quantum computing and quantum approximate optimization algorithm_ , arXiv:1909.05500 (2019).
  * [9] J. Wen, X. Kong, S. Wei, B. Wang, T. Xin, and G. Long, _Experimental realization of quantum algorithms for a linear system inspired by adiabatic quantum computing_ , Physical Review A 99, 012320 (2019) [arXiv:1806.03295].
  * [10] C. Bravo-Prieto, R. LaRose, M. Cerezo, Y. Subaşı, L. Cincio, and P. J. Coles, _Variational quantum linear solver: A hybrid algorithm for linear systems_ , arXiv:1909.05820 (2019).
  * [11] H. Y. Huang, K. Bharti, and P. Rebentrost, _Near-term quantum algorithms for linear systems of equations_ , arXiv:1909.07344.
  * [12] L. Lin and Y. Tong, _Optimal quantum eigenstate filtering with application to solving quantum linear systems_ , Quantum 4, 361 (2020) [arXiv:1910.14596].
  * [13] S. Aaronson, _Read the fine print_ , Nature Physics 11, 291-293 (2015) [citeseerx].
  * [14] A. Montanaro and S. Pallister, _Quantum algorithms and the finite element method_ , Physical Review A 93, 032324 (2016) [arXiv:1512.05903].
  * [15] R. Babbush, J. McClean, C. Gidney, S. Boixo, and H. Neven, _Focus beyond quadratic speedups for error-corrected quantum advantage_ , PRX Quantum 2 (2021) [arXiv:2011.04149].
  * [16] A. N. Chowdhury and R. D. Somma, _Quantum algorithms for Gibbs sampling and hitting-time estimation_ , Quant. Inf. Comp. 17, 0041–0064 (2017) [arXiv:1603.02940].
  * [17] B. D. Clader, B. C. Jacobs, and C. R. Sprouse, _Preconditioned quantum linear system algorithm_ , Physical Review Letters 110, 250504 (2013) [arXiv:1301.2340].
  * [18] J. R. Shewchuk, _An introduction to the conjugate gradient method without the agonizing pain_ , Carnegie Mellon University (1994).
  * [19] S. Chakraborty, A. Gilyén, and S. Jeffery, _The power of block-encoded matrix powers: improved regression techniques via faster hamiltonian simulation_ , 46th International Colloquium on Automata, Languages, and Programming (ICALP 2019) [arXiv:1804.01973].
  * [20] C. Shao and H. Xiang, _Quantum circulant preconditioner for linear system of equations_ , Physical Review A 98, 062321 [arXiv:1807.04563].
  * [21] Y. Tong, D. An, N. Wiebe, and L. Lin. _Fast inversion, preconditioned quantum linear system solvers, fast Green’s-function computation, and fast evaluation of matrix functions_ , Physical Review A 104, 032422 [arXiv:2008.13295].
  * [22] B. Wu, M. Ray, L. Zhao, X. Sun, and P. Rebentrost, _Quantum-classical algorithms for skewed linear systems with optimized Hadamard test_ , Physical Review A 103, 042422 [arXiv:2009.13288].
  * [23] A. C. Vazquez, R. Hiptmair, and S. Woerner, _Enhancing the Quantum Linear Systems Algorithm using Richardson Extrapolation_ , arXiv:2009.04484 (2020).
  * [24] G. H. Low and I. L. Chuang, _Optimal Hamiltonian simulation by quantum signal processing_ , Physical Review Letters 118, 010501 (2017) [arXiv:1606.02685].
  * [25] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, _Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics_ , 51st Annual ACM SIGACT Symposium on Theory of Computing, 193–204 (2019) [arXiv:1806.01838].
  * [26] A. M. Childs and N. Wiebe, _Hamiltonian Simulation Using Linear Combinations of Unitary Operations_ , Quantum Information & Computation [arXiv:1202.5822].
  * [27] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari and R. D. Somma, _Simulating Hamiltonian dynamics with a truncated Taylor series_ , Physical Review Letters 114, 090502 (2015) [arXiv:1412.4687].
  * [28] G. H. Low and I. L. Chuang, _Hamiltonian simulation by qubitization_ , Quantum 3, 163 (2019).
  * [29] A. C. Schaeffer, _Inequalities of A. Markoff and S. Bernstein for polynomials and related functions_ , Bulletin of the American Mathematical Society 47, 565–579(1941).
  * [30] G. Brassard, P. Hoyer, M. Mosca, and A. Tapp, _Quantum Amplitude Amplification and Estimation_ , Contemporary Mathematics 305, 53–74 (2002) [arXiv:0005055].
  * [31] See e.g. the Cholesky decomposition page on Wikipedia.
  * [32] R. D. Somma and S. Boixo, _Spectral gap amplification_ , SIAM Journal on Computing 42, 593-610 (2013) [arXiv:1110.2494].
  * [33] M. A. Nielsen, and I. Chuang, _Quantum computation and quantum information_ , Cambridge University Press (2000).
  * [34] L. Grover and T. Rudolph, _Creating superpositions that correspond to efficiently integrable probability distributions_ , arXiv:0208112 (2002).
  * [35] V. Giovannetti, S. Lloyd, and L. Maccone, _Quantum random access memory_ , Physical Review Letters 100, 160501 (2008) [arXiv:0708.1879].
  * [36] I. Kerenidis and A. Prakash, _Quantum recommendation systems_ , arXiv:1603.08675.
  * [37] M. Boyer,G. Brassard, P. Høyer and A. Tapp, _Tight bounds on quantum searching_ , Progress of Physics 46, 493–505 (1998) [arXiv:9605034].
  * [38] András Gilyén, private communication.
  * [39] R. Chao, D. Ding, A. Gilyén, C. Huang and M. Szegedy, _Finding angles for quantum signal processing with machine precision_ , arXiv:2003.02831 (2020).
  * [40] Y. Dong, X. Meng, K. B. Whaley and L. Lin, _Efficient phase factor evaluation in quantum signal processing_ , Phys. Rev. A 103, 042419 (2021) [arXiv:2002.11649].
  * [41] See e.g. the Gershgorin circle theorem page on Wikipedia.
  * [42] V. V. Shende, S. S. Bullock, and I. L. Markov, _Synthesis of quantum-logic circuits_ , IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 25, 1000–1010 (2006) [arXiv:0406176].
  * [43] R. Merris, _Laplacian matrices of graphs: a survey_ , Linear algebra and its applications 197, 143–176 (1994).
  * [44] D. A. Spielman, _Algorithms, graph theory, and linear equations in Laplacian matrices_ , Proceedings of the International Congress of Mathematicians 2010, 2698–2722 (2010).
  * [45] L. K. Grover, _Synthesis of quantum superpositions by quantum computation_ , Physical Review Letters 85, 1334 (2000).
  * [46] Y. R. Sanders, G. H. Low, A. Scherer and D. W. Berry, _Black-box quantum state preparation without arithmetic_ , Physical Review Letters 122, 020502 (2019) [arXiv:1807.03206].
  * [47] R. D. Somma and Y. Subaşı, _Quantum state verification in the quantum linear systems problem_ , PRX Quantum 2, 010315 (2021) [arXiv:2007.15698].
  * [48] A. Gilyén, _Quantum walk based search methods and algorithmic applications_ , Doctoral dissertation, Eötvös Loránd University (2014).
  * [49] E. Malvetti, R. Iten, and R. Colbeck, _Quantum Circuits for Sparse Isometries_ arXiv:2006.00016 (2020).
  * [50] X. Jiang, _Minimum rank positive semidefinite matrix completion with chordal sparsity pattern_ , Doctoral dissertation, UCLA (2017).
  * [51] A. Nayak and F. Wu, _The quantum query complexity of approximating the median and related statistics_ , Proceedings of the 31st annual ACM symposium on Theory of computing, 384–393 (1999) [arXiv:9804066].
  * [52] S. U. Pillai, T. Suel, and S. Cha, _The Perron-Frobenius theorem: some of its applications_ , IEEE Signal Processing Magazine 22, 62–75 (2005).
  * [53] S. Arora and B. Barak, _Computational complexity: a modern approach_ , Cambridge University Press (2009).
  * [54] J. C. Mason, and D. C. Handscomb, _Chebyshev polynomials_ , CRC press (2002).
  * [55] J. Bausch and E. Crosson, _Analysis and limitations of modified circuit-to-Hamiltonian constructions_ , Quantum 2, 94 (2018).