A quantum linear system algorithm for dense matrices

TL;DR

This paper introduces a quantum algorithm for solving dense linear systems with a runtime that improves quadratically over previous quantum methods, especially for matrices with bounded spectral norm, leveraging singular value estimation and efficient state preparation.

Contribution

The paper presents a novel quantum linear system algorithm with sparsity-independent runtime scaling for dense matrices, improving efficiency over existing quantum algorithms.

Findings

Achieves $ ilde{O}( au ext{polylog}(n)/ extepsilon)$ runtime for dense matrices.

Provides a quadratic speedup over previous quantum algorithms for dense matrices.

Utilizes a memory architecture for efficient quantum state preparation.

Abstract

Solving linear systems of equations is a frequently encountered problem in machine learning and optimisation. Given a matrix $A$ and a vector $\mathbf b$ the task is to find the vector $\mathbf x$ such that $A \mathbf x = \mathbf b$. We describe a quantum algorithm that achieves a sparsity-independent runtime scaling of $\mathcal{O}(\kappa^2 \|A\|_F \text{polylog}(n)/\epsilon)$, where $n\times n$ is the dimensionality of $A$ with Frobenius norm $\|A\|_F$, $\kappa$ denotes the condition number of $A$, and $\epsilon$ is the desired precision parameter. When applied to a dense matrix with spectral norm bounded by a constant, the runtime of the proposed algorithm is bounded by $\mathcal{O}(\kappa^2\sqrt{n} \text{polylog}(n)/\epsilon)$, which is a quadratic improvement over known quantum linear system algorithms. Our algorithm is built upon a singular value estimation subroutine, which makes use of a memory architecture that allows for efficient preparation of quantum states that correspond to the rows and row Frobenius norms of $A$.

Equations

$$U_{\mathcal{M}}:|i⟩|0⟩\to |i,{\mathrm{A}}_{\mathrm{i}}⟩$$

$$U_{\mathcal{N}}:|0⟩|j⟩\to |{\mathrm{A}}_F,j⟩$$

$$\sum _i{\alpha }_{{\mathrm{v}}_i}|{\mathrm{v}}_i⟩|\overline{{\sigma }_i}⟩$$

$$\mathcal{M}:|\alpha ⟩$$

$$\mathcal{N}:|\alpha ⟩=\sum _{j=1}^n{\alpha }_j|j⟩\to \sum _{j=1}^n{\alpha }_j|{\mathrm{A}}_F,j⟩=|\mathcal{N}\alpha ⟩.$$

$$W|\mathcal{N}{\mathrm{v}}_i⟩=$$

$$=$$

$$\cos {\theta }_i$$

$$=\frac{2{\sigma }_i}{\Vert A{\Vert }_F}⟨{\mathrm{v}}_i|A^{†}|{\mathrm{u}}_i⟩-1$$

$$=\frac{2{\sigma }_i^2}{\Vert A{\Vert }_F^2}-1,$$

$$\sum _i{\beta }_i{|{\mathrm{v}}_i⟩}_A{||{\overline{\lambda }}_i|⟩}_B{||{\overline{\lambda }}_i+\mu |⟩}_C.$$

$$\sum _i(-1{)}^{f_i}{\beta }_i{|{\mathrm{v}}_i⟩}_A{||{\overline{\lambda }}_i|⟩}_B{||{\overline{\lambda }}_i+\mu |⟩}_C{|f_i⟩}_D.$$

$$\sum _i(-1{)}^{f_i}{\beta }_i|{\mathrm{v}}_i⟩\left(\frac{\gamma }{\overline{{\lambda }_i}}|0⟩+\sqrt{1-{\frac{\gamma }{\overline{{\lambda }_i}}}^2}|1⟩\right)$$

$$\begin{gathered} |\mathrm{h}(\lambda )⟩:=\sqrt{1-\mathrm{f}(\lambda {)}^2-\mathrm{g}(\lambda {)}^2}|\mathrm{NO}⟩ \\ +\mathrm{f}(\lambda )|\mathrm{WC}⟩+\mathrm{g}(\lambda )|\mathrm{IC}⟩, \end{gathered}$$

$$|\psi ⟩:=Q|\mathbf{b}⟩|0⟩\to \sum _i{\beta }_i|{\mathbf{v}}_i⟩|\mathrm{h}({\lambda }_i)⟩.$$

$$|\overline{\psi }⟩:=\overline{Q}|\mathbf{b}⟩|0⟩\to \sum _i{\beta }_i|{\mathbf{v}}_i⟩|\mathrm{h}({\overline{\lambda }}_i)⟩.$$

$$\Vert |\mathrm{h}({\lambda }_i)⟩-|\mathrm{h}({\lambda }_j)⟩\Vert \le c\kappa |{\lambda }_i-{\lambda }_j|,$$

$$Re⟨\overline{\psi }|\psi ⟩=\sum _{i=1}^N|{\beta }_i{|}^{2}Re⟨\mathrm{h}({\overline{\lambda }}_i)|\mathrm{h}({\lambda }_i)⟩$$

$$\ge \sum _{i=1}^N|{\beta }_i{|}^2\left(1-\frac{c^2{\kappa }^2{\delta }^2{\Vert A\Vert }_F^2}{2}\right),$$

$$1-Re⟨\overline{\psi }|\psi ⟩\le \sum _{i=1}^N|{\beta }_i{|}^2\left(\frac{c^2{\kappa }^2{\delta }^2{\Vert A\Vert }_F^2}{2}\right)$$

$$\Vert |\overline{\psi }⟩-|\psi ⟩\Vert \le \mathcal{O}\left(\kappa \delta {\Vert A\Vert }_F\right).$$

Full text

Introduction.

A common bottleneck in statistical learning and machine learning algorithms is the inversion of high-dimensional matrices in order to solve linear systems of equations. Examples include covariance matrix inversions in Gaussian processes and support vector machines, as well as data matrix inversions in large scale regression problems Rasmussen and Williams (2004); Bishop (2006).

Recent advances in the field of quantum information processing have provided promising prospects for the efficient solution of high-dimensional linear systems. The breakthrough work of Harrow, Hassidim and Lloyd Harrow et al. (2009) introduced the quantum linear system algorithm (QLSA) that computes the quantum state $\ket{\mathbf{x}}=\ket{A^{-1}\mathbf{b}}$ corresponding to the solution of a linear system $A\mathbf{x}=\mathbf{b}$, where $A\in\mathbb{R}^{n\times n}$ and $\mathbf{x},\mathbf{b}\in\mathbb{R}^{n}$, in time $\mathcal{O}\left(\text{polylog}(n)\right)$ for a sparse and well-conditioned $A$. Unlike the output $A^{-1}\mathbf{b}\in\mathbb{R}^{n}$ of a classical linear system solver, a copy of $\ket{A^{-1}\mathbf{b}}$ does not provide access to the coordinates of $A^{-1}\mathbf{b}$. Nevertheless, it allows us to perform useful computations such as sampling from the solution vector. The QLSA algorithm has inspired several works Rebentrost et al. (2014); Schuld et al. (2016); Wiebe et al. (2012); Wiebe and Granade (2015); Wiebe et al. (2015a, b, 2014, 2016); Zeng and Coecke (2016); Zhao et al. (2015) in the emerging research area of quantum machine learning.

In the classical setting, the best known algorithm for the sampling task performed by the QLSA algorithm requires solving the linear system. The running time for a classical linear system solver scales as $\mathcal{O}\left(n^{\omega}\right)$, where the matrix multiplication exponent $\omega\leq 2.373$ Coppersmith and Winograd (1990); Le Gall (2014). However, as the sub-cubic scaling is difficult to achieve in practice, linear system solvers typically use the Cholesky decomposition and require time $\mathcal{O}\left(n^{3}\right)$ for dense matrices.

The QLSA algorithm Harrow et al. (2009) has running time $\tilde{O}(\kappa^{2}s(A)^{2}/\epsilon)$ where $\kappa$ is the condition number, $s(A)$ is the sparsity or the maximum number of non-zero entries in a row of $A$ and $\epsilon$ is the precision to which the solution is approximated. There have been several improvements to the QLSA algorithm since the original proposal that have improved the running time to linear in $\kappa$ and $s(A)$ and to poly-logarithmic in the precision parameter $\epsilon$ Childs et al. (2015); Ambainis (2012). The work Clader et al. (2013) introduced pre-conditioning for the QLSA algorithm and extended its applicability.

Quantum machine learning is an emerging research area that attempts to harness the power of quantum information processing to obtain speedups for classical machine learning tasks. A number of quantum machine learning algorithms have been proposed Rebentrost et al. (2014); Schuld et al. (2016); Wiebe et al. (2012); Wiebe and Granade (2015); Wiebe et al. (2015a, b, 2014, 2016); Zeng and Coecke (2016); Zhao et al. (2015). Most of these algorithms use a quantum linear system solver as a subroutine. However, as mentioned in Harrow et al. (2009), and later also pointed out in Childs (2009); Aaronson (2015), the QLSA potentially has a few caveats. In particular, the QLSA achieves an exponential speedup over classical algorithms when the matrix $A$ is sparse and well conditioned, due to the sparsity-dependent Hamiltonian simulation subroutine. The potential exponential advantage of QLSA is lost when it is applied to dense matrices, which consistute a large class of interesting applications. Examples include kernel methods,  Wilson et al. (2015), and artificial neural networks, where particularly convolutional neural network architectures rely heavily on subroutines that manipulate large, non-sparse matrices Yangqing (2014); Chetlur and Woolley (2014). Alternative approaches to the quantum linear system problem that avoid sparsity dependence are therefore desirable for a more general application of quantum computing to classical learning problems.

In this letter we present a quantum algorithm for solving linear systems of equations using the quantum singular value estimation (QSVE) algorithm introduced in Kerenidis and Prakash (2017). The proposed algorithm achieves a runtime $\mathcal{O}\left(\kappa^{2}\left\lVert{A}\right\rVert_{F}\cdot\text{polylog}(n)/\epsilon\right)$, where $\kappa$ denotes the condition number of $A$, $\left\lVert{A}\right\rVert_{F}$ is the Frobenius norm and $\epsilon$ is the precision parameter. When the spectral norm $\left\lVert{A}\right\rVert_{*}$ is bounded by a constant, the scaling becomes $\mathcal{O}\left(\kappa^{2}\sqrt{n}\cdot\text{polylog}(n)/\epsilon\right)$, which amounts to a polynomial speed-up over the $\mathcal{O}\left(\kappa^{2}n\cdot\text{polylog}(n)/\epsilon\right)$ scaling achieved by Harrow et al. (2009) when applied to dense matrices.

We start by introducing some preliminaries. For a symmetric matrix $A\in\mathbb{R}^{n\times n}$ with spectral decomposition $A=\sum_{i\in[n]}\lambda_{i}\mathbf{s}_{i}\mathbf{s}_{i}^{\dagger}$, the singular value decomposition is given by $A=\sum_{i}^{r}|\lambda_{i}|\mathbf{s}_{i}\mathbf{s}_{i}^{\dagger}$. We also need the well known quantum phase estimation algorithm:

Theorem 1** (Phase estimation Kitaev (1995)).**

Let unitary $U\ket{v_{j}}=\exp{(i\theta_{j})}\ket{v_{j}}$ with $\theta_{j}\in[-\pi,\pi]$ for $j\in[n]$. There is a quantum algorithm that transforms $\sum_{j\in[n]}\alpha_{j}\ket{v_{j}}\to\sum_{j\in[n]}\alpha_{j}\ket{v_{j}}\ket{\overline{\theta}_{j}}$ such that $|\overline{\theta_{j}}-\theta_{j}|\leq\delta$ for all $j\in[n]$ with probability $1-1/\text{poly}(n)$ in time $\mathcal{O}\left(T_{U}\log{(n)}/\delta\right)$, where $T_{U}$ defines the time to implement $U$.

Quantum singular value estimation can be viewed as an extension of phase estimation to non unitary matrices. It is the main algorithmic primitive required for our linear system solver.

Definition 1** (Quantum singular value estimation).**

Let $A\in\mathbb{R}^{m\times n}$ have singular value decomposition $A=\sum_{i}\sigma_{i}u_{i}v_{i}^{t}$. A quantum singular value estimation algorithm with precision $\delta$ transforms $\sum_{j\in[n]}\alpha_{j}\ket{v_{j}}\to\sum_{j\in[n]}\alpha_{j}\ket{v_{j}}\ket{\overline{\sigma}_{j}}$ such that $|\overline{\sigma_{j}}-\sigma_{j}|\leq\delta$ for all $j\in[n]$ with probability $1-1/\text{poly}(n)$.

A quantum singular value estimation (QSVE) algorithm with running time of $\tilde{O}(\left\lVert{A}\right\rVert_{F}/\delta)$ was presented in Kerenidis and Prakash (2017), where it was used for quantum recommendation systems. An SVE algorithm applied to a symmetric matrix estimates $|\lambda_{i}|$ but does not provide an estimate for $sign(\lambda_{i})$. However, in order to solve linear systems we also need to recover $sign(\lambda_{i})$. We provide a simple procedure for recovering the sign given an SVE algorithm. Our procedure provides a way to construct a quantum linear system solver from a QSVE algorithm in a black box manner.

The main result of this letter is a quantum linear system solver based on the QSVE algorithm Kerenidis and Prakash (2017) that achieves a running time of $\mathcal{O}\left(\kappa^{2}\left\lVert{A}\right\rVert_{F}\cdot\text{polylog}(n)/\epsilon\right)$. We briefly describe the QSVE algorithm in the next section. We then present the quantum linear system solver and provide a complete analysis for the linear system solver as well as a comparison with other approaches to the QLSA in the discussion.

The QSVE algorithm.

The QSVE algorithm requires the ability to efficiently prepare the quantum states corresponding to the rows and columns of matrix $A$. The matrix entries are stored in the following data structure, such that a quantum algorithm with access to this data structure has this ability.

Lemma 1** (Data Structure  Kerenidis and Prakash (2017)).**

Let $A\in\mathbb{R}^{m\times n}$ be a matrix with entries $A_{ij}$ which arrive in an arbitrary order. There exists a data structure with the following properties:

• •

A quantum computer with access to the data structure can perform the following mappings in $\mathcal{O}\left(\text{polylog}(mn)\right)$ time.

[TABLE]

where $\bm{\mathrm{A_{i}}}\in R^{n}$ correspond to the rows of the matrix $A$ and $\bm{\mathrm{A}}_{F}\in\mathbb{R}^{m}$ is a vector whose entries are the $\ell_{2}$ norms of the rows, i.e. $(\bm{\mathrm{A}}_{F})_{i}=\left\lVert{A_{i}}\right\rVert$.

• •

The time required to store a new entry $A_{ij}$ is $\mathcal{O}\left(\text{log}^{2}(mn)\right)$ and data structure size is $O(w\log mn)$ where $w$ is the number of non zero entries in $A$.

A possible realization of this data structure is based on an array of $m$ binary trees, each binary tree contains at most $n$ leaves which store the squared amplitudes of the corresponding matrix entry $|A_{ij}|^{2}$, as well as the sign of $A_{ij}$. An internal node of a tree stores the sum of the elements in the subtree rooted at it. The root of the $i^{th}$ tree then contains $\left\lVert{\mathbf{A}_{i}}\right\rVert^{2},\,i\in[m]$. In order to access the vector of row Frobenius norms, we construct one more binary tree, the $i^{th}$ leaf of which stores $\left\lVert{\mathbf{A}_{i}}\right\rVert^{2}$. A detailed description of such a binary tree memory structure, and the proof of Lemma 1 can be found in Kerenidis and Prakash (2017).

The QSVE algorithm is a quantum walk based algorithm that leverages the connection between the singular values $\sigma_{i}$ of the target matrix $A$ and the principal angles $\theta_{i}$ between certain subspaces associated with $A$. The relation between quantum walks and eigenvalues has been well known in the literature and has been used in several previous results Szegedy (2004); Childs (2010). However, the quantum walk defined by the QSVE algorithm is particularly interesting for linear systems as instead of the sparsity $s(A)$, it depends on the Frobenius norm $\left\lVert{A}\right\rVert_{F}$.

The QSVE algorithm makes use of a factorization $\frac{A}{\|A\|_{F}}=\mathcal{M}^{\dagger}\mathcal{N}$, where $\mathcal{M}\in\mathbb{R}^{mn\times m}$ and $\mathcal{N}\in\mathbb{R}^{mn\times n}$ are isometries. The key idea is that the unitary operator $W$ defined by $W=(2\mathcal{M}\mathcal{M}^{\dagger}-I_{mn})(2\mathcal{N}\mathcal{N}^{\dagger}-I_{mn})$, where $I_{mn}$ is the identity matrix, can be implemented efficiently using the data structure in Lemma 1. Further, $W$ has two dimensional eigenspaces spanned by $\{\mathcal{M}\mathbf{u}_{i},\mathcal{N}\mathbf{v}_{i}\}$ on which it acts as a roation by angle $\theta_{i}$, such that $\cos\frac{\theta_{i}}{2}=\frac{\sigma_{i}}{\|A\|_{F}}$. Note that the sub-space spanned by $\{\mathcal{M}\mathbf{u}_{i},\mathcal{N}\mathbf{v}_{i}\}$ is therefore also spanned by $\{\mathbf{w}_{i}^{+},\mathbf{w}_{i}^{-}\}$, the eigenvectors of $W$ with eigenvalues $\exp(i\theta_{i})$ and $\exp(-i\theta_{i})$ respectively. In particular we may write the following decomposition, $\ket{\mathcal{N}\mathbf{v}_{i}}=\omega_{i}^{+}\ket{\mathbf{w}_{i}^{+}}+\omega_{i}^{-}\ket{\mathbf{w}_{i}^{-}}$, with $|\omega_{i}^{-}|^{2}+|\omega_{i}^{+}|^{2}=1$. Algorithm 1 describes the QSVE algorithm, the analysis is contained in the following lemma.

Lemma 2** (Preparation of the Isometries Kerenidis and Prakash (2017)).**

Let $A\in\mathbb{R}^{m\times n}$ be a matrix with singular value decomposition $A=\sum_{i}\sigma_{i}\bm{\mathrm{u}}_{i}\bm{\mathrm{v}}_{i}^{\dagger}$ stored in the data structure described in Lemma 1. Then there exist matrices $\mathcal{M}\in\mathbb{R}^{mn\times m}$, and $\mathcal{N}\in\mathbb{R}^{mn\times n}$, such that

$\mathcal{M},\mathcal{N}$ are isometries, that is $\mathcal{M}^{\dagger}\mathcal{M}=I_{m}$ and $\mathcal{N}^{\dagger}\mathcal{N}=I_{n}$ such that $A$ can be factorized as $A/\left\lVert{A}\right\rVert_{F}=\mathcal{M}^{\dagger}\mathcal{N}$.

Multiplication by $\mathcal{M},\mathcal{N}$, i.e. the mappings $\ket{\alpha}\rightarrow\ket{\mathcal{M}\alpha}$ and $\ket{\beta}\rightarrow\ket{\mathcal{N}\beta}$ can be performed in time $\mathcal{O}\left(\text{polylog}(mn)\right)$, 2. 2.

The reflections $2\mathcal{M}\mathcal{M}^{\dagger}-I_{mn}$, $2\mathcal{N}\mathcal{N}^{\dagger}-I_{mn}$, and hence the unitary $W=(2\mathcal{M}\mathcal{M}^{\dagger}-I_{mn})(2\mathcal{N}\mathcal{N}^{\dagger}-I_{mn})$ can be implemented in time $\mathcal{O}\left(\text{polylog}(mn)\right)$. 3. 3.

The unitary $W$ acts as rotation by $\theta_{i}$ on the two dimensional invariant subspace $\{\mathcal{M}\mathbf{u}_{i},\mathcal{N}\mathbf{v}_{i}\}$ plane, such that $\sigma_{i}=\cos\frac{\theta_{i}}{2}\|A\|_{F}$, where $\sigma_{i}$ is the $i$-th singular value for $A$.

We outline the ideas involved in the analysis of the QSVE algorithm and refer to Kerenidis and Prakash (2017) for further details. The map $\mathcal{M}$ appends to an arbitrary input state vector $\ket{\alpha}$ a register that encodes the row vectors $\mathbf{A_{i}}$ of $A$, such that

[TABLE]

The map $\mathcal{N}$ similarly appends to an arbitrary input state vector $\ket{\alpha}$ a register that encodes the vector $\bm{\mathrm{A_{F}}}$ whose entries are the $\ell_{2}$ norms $\|\bm{\mathrm{A_{i}}}\|_{F}$ of the rows of $A$,

[TABLE]

The above maps can be efficiently implemented given the memory structure described by Lemma 1.

The factorization of $A$ follows from the amplitude encoding of $\bm{\mathrm{A_{i}}}$ and $\bm{\mathrm{A}}_{F}$. We have $\ket{i,\bm{\mathrm{A_{i}}}}=\frac{1}{\|\bm{\mathrm{A_{i}}}\|}\sum\limits_{j=1}^{n}A_{ij}\ket{i,j}$ and $\ket{\bm{\mathrm{A}}_{F},j}=\frac{1}{\|A\|_{F}}\sum\limits_{i=1}^{m}\|\bm{\mathrm{A_{i}}}\|\ket{i,j}$, implying that $(\mathcal{M}^{\dagger}\mathcal{N})_{ij}=\left<i,\bm{\mathrm{A_{i}}}\left|\bm{\mathrm{A}}_{F},j\right.\right>=\frac{A_{ij}}{\|A\|_{F}}$. Similarly, it follows that $\mathcal{M}$ and $\mathcal{N}$ have orthonormal columns and thus $\mathcal{M}^{\dagger}\mathcal{M}=I_{m}$ and $\mathcal{N}^{\dagger}\mathcal{N}=I_{n}$.

To show the relation between the eigenvalues of $W$ and the singular values of $A$, we consider the following:

[TABLE]

where we used the singular value decomposition $A=\sum_{i}\sigma_{i}\ket{\bm{\mathrm{u}}_{i}}\bra{\bm{\mathrm{v}}_{i}}$, and the fact that the right singular vectors $\{\bm{\mathrm{v}}_{i}\}$ are mutually orthonormal. Note that $W$ rotates $\ket{\mathcal{N}\bm{\mathrm{v}}_{i}}$ in the plane of $\{\mathcal{M}\mathbf{u}_{i},\mathcal{N}\mathbf{v}_{i}\}$ by $\theta_{i}$, such that

[TABLE]

where we have used the fact that $(2\mathcal{M}\mathcal{M}^{\dagger}-I_{mn})$ represents a reflection in $\ket{\mathcal{M}\bm{\mathrm{u}}_{i}}$ and that $A^{\dagger}=\mathcal{N}^{\dagger}\mathcal{M}=\sum_{i}\sigma_{i}\ket{\bm{\mathrm{v}}_{i}}\bra{\bm{\mathrm{u}}_{i}}$. Therefore the angle between $\ket{\mathcal{N}\bm{\mathrm{v}}_{i}}$ and $\ket{\mathcal{M}\bm{\mathrm{u}}_{i}}$ is given by $\frac{\theta_{i}}{2}$, i.e. half of the total rotation angle. Comparing the above expression with the half-angle formula for cosine functions, we obtain the relation $\cos\left(\frac{\theta_{i}}{2}\right)=\frac{\sigma_{i}}{\|A\|_{F}}$.

The two dimensional sub-space spanned by $\{\mathcal{M}\mathbf{u}_{i},\mathcal{N}\mathbf{v}_{i}\}$ is therefore invariant under the action of $W$ which acts on it as a rotation by angle $\theta_{i}$. The corresponding eigenvectors of $W$ have hence eigenvalues $\exp(\pm i\theta_{i})$, and in particular we can perform phase estimation to get an estimate $\pm\overline{\theta_{i}}$ and then compute $\overline{\sigma_{i}}=\cos(\overline{\theta_{i}}/2)\left\lVert{A}\right\rVert_{F}$ to estimates the singular values. We therefore have established the correctness of Algorithm 1 and have the following theorem,

Theorem 2** (Quantum Singular Value Estimation Kerenidis and Prakash (2017)).**

Let $A\in\mathbb{R}^{m\times n}$ be a matrix with singular value decomposition $A=\sum_{i}\sigma_{i}\bm{\mathrm{u}}_{i}\bm{\mathrm{v}}_{i}^{\dagger}$ stored in the data structure in Lemma 1. Further let $\delta>0$ be the precision number. There is an algorithm that runs in $\mathcal{O}\left(\text{polylog}(mn)/\delta\right)$ and performs the mapping $\sum_{i}\alpha_{\bm{\mathrm{v}}_{i}}\ket{\bm{\mathrm{v}}_{i}}\ket{0}\rightarrow\sum_{i}\alpha_{\bm{\mathrm{v}}_{i}}\ket{\bm{\mathrm{v}}_{i}}\ket{\overline{\sigma_{i}}}$ where $\overline{\sigma_{i}}\in\sigma_{i}\pm\delta\left\lVert{A}\right\rVert_{F}$ for all $i$ with probability at least $1-1/\text{poly}(n)$.

The runtime of QSVE is dominated by the phase estimation procedure which returns an $\delta$-close estimate of $\theta_{i}$, s.t. $|\overline{\theta}_{i}-\theta_{i}|\leq 2\delta$, which translates into the estimated singular value via $\overline{\sigma}_{i}=\cos{(\overline{\theta}_{i}/2)}\left\lVert{A}\right\rVert_{F}$. The error in $\sigma_{i}$ can then be bounded from above by $|\overline{\sigma}_{i}-\sigma_{i}|\leq\delta\left\lVert{A}\right\rVert_{F}$. The unitary $W$ can be implemented in time $\mathcal{O}\left(\text{polylog}(mn)\right)$ by Lemma 2, by Theorem 1 the running time for estimating of the singular values with additive error $\delta\left\lVert{A}\right\rVert_{F}$ in $\mathcal{O}\left(\text{polylog}(mn)/\delta\right)$.

Quantum linear system algorithm.

Without loss of generality we can assume that the matrix $A$ is Hermitian as it is well known that the general case can be reduced to the Hermitian case Harrow et al. (2009). A QSVE algorithm immediately yields a linear system solver for positive definite matrices as the estimated singular values and eigenvalues are related via $\overline{\sigma_{i}}=|\overline{\lambda_{i}}|$. In order to solve general linear systems we need to recover the sign of each $\overline{\lambda_{i}}$. We provide a simple algorithm that recovers the signs using the QSVE procedure as a black box incurring only a constant overhead over the QSVE.

The main result of this letter is the following theorem.

Theorem 3**.**

Let $A\in\mathbb{R}^{n\times n}$ be a Hermitian matrix with spectral decomposition $A=\sum_{i}\lambda_{i}\mathbf{\bm{\mathrm{u}}_{i}\bm{\mathrm{u}}_{i}}^{\dagger}$ stored in the data structure in Lemma 1. Further let $\kappa$ be the condition number $A$, and $\left\lVert{A}\right\rVert_{F}$ the Frobenius norm and $\epsilon>0$ be a precision parameter. Then Algorithm 2 has runtime $\mathcal{O}\left(\kappa^{2}\cdot\text{polylog}(n)\cdot\left\lVert{A}\right\rVert_{F}/\epsilon\right)$ that outputs state $\ket{\overline{A^{-1}\mathbf{b}}}$ such that $\left\lVert{\ket{\overline{A^{-1}\mathbf{b}}}-\ket{A^{-1}\mathbf{b}}}\right\rVert\leq\epsilon$.

Proof.

We first argue that Algorithm 2 correctly recovers the sign of the $\lambda_{i}$. The algorithm compares the estimates obtained by performing QSVE for $A$ and for $A^{\prime}=A+\mu I_{n}$, where $\mu$ is a positive scalar to be chosen later. The matrix $A^{\prime}$ has the same eigenvectors as $A$, but has eigenvalues $\lambda_{i}+\mu$. Note that for $\lambda_{i}\geq 0$ we have $|\lambda_{i}+\mu|=|\lambda_{i}|+|\mu|\geq|\lambda_{i}|$, however if $\lambda_{i}\leq-\mu/2$ then $|\lambda_{i}+\mu|\leq|\lambda_{i}|$.

Thus, if the estimates were perfect, then choosing $\mu=2/\kappa$ would recover the sign correctly when the eigenvalues of $A$ lie in the interval $[-1,-1/\kappa]\cup[1/\kappa,1]$. With the choice $\mu=4/\kappa$ and $\epsilon=1/\kappa$ we find that the signs are still correct for all $\lambda_{i}$.

We outline the derivation of the runtime, the analysis of the error bounds appears in the appendix. The running time for QSVE with precision $\epsilon=1/\kappa$ is $\tilde{O}(\kappa\left\lVert{A}\right\rVert_{F})$. Considering the success probability of the post-selection step, we require on average $\mathcal{O}\left(\kappa^{2}\right)$ repetitions of the coherent computation. This can be reduced to $\mathcal{O}\left(\kappa\right)$ using amplitude amplification Brassard et al. (2002). Therefore an upper-bound of the runtime of our algorithm is given by $\mathcal{O}\left(\kappa^{2}\cdot\text{polylog}(n)\left\lVert{A}\right\rVert_{F}/\epsilon\right)$. ∎

Discussion.

The error dependence on the Frobenius norm suggests that our algorithm is most accurate when the $\left\lVert{A}\right\rVert_{F}$ is bounded by some constant, in which case the algorithm returns the output state with a constant $\epsilon$-error in polylogarithmic time even if the matrix is non-sparse. More generally, as in the QLSA algorithm we can assume that the spectral norm $\left\lVert{A}\right\rVert_{*}$ is bounded by a constant, although the Frobenius norm may scale with the dimensionality of the matrix. In such cases we have $\left\lVert{A}\right\rVert_{F}=\mathcal{O}\left(\sqrt{n}\right)$. Hence in such scenarios the proposed algorithm runs in $\mathcal{O}\left(\kappa^{2}\sqrt{n}\cdot\text{polylog}(n)/\epsilon\right)$ and returns the output with a constant $\epsilon$-error.

It was shown in Berry and Childs (2012) that given black-box access to the matrix elements, Hamiltonian simulation with error $\delta_{h}$ can be performed in time $\mathcal{O}\left(n^{2/3}\cdot\text{polylog}(n)/\delta_{h}^{1/3}\right)$. This leads to a linear system algorithm based on Harrow et al. (2009) which scales as $\mathcal{O}\left(\kappa^{2}n^{2/3}\cdot\text{polylog}(n)/\epsilon\right)$, where we have assumed the dominant error comes from phase estimation, and hence the error introduced by the Hamiltonian simulation is neglected. It was also shown numerically that the method of Berry and Childs (2012) attains a typical scaling of $\mathcal{O}\left(\sqrt{n}\cdot\text{polylog}(n)/\delta_{h}^{1/2}\right)$ when applied to randomly selected matrices, leading to a $\mathcal{O}\left(\kappa^{2}\sqrt{n}\cdot\text{polylog}(n)/\epsilon\right)$ linear system algorithm. The work  Berry and Childs (2012) assumes that we have black-box access to the matrix entries, that is quantum queries of the form $\ket{i,j,0}\to\ket{i,j,A_{ij}}$ are allowed. We note that in Lemma 1 we instead assume black-box access to a data structure constructed in linear time from a stream of the matrix entries. Our access model thus differs from the one used in  Berry and Childs (2012) and therefore a direct comparison of the two results is not appropriate. The QSVE-based linear system solver achieves a $O(\sqrt{n})$-scaling with quantum access to the data structure in Lemma 1, and it is an interesting open question if one can achieve a similar scaling in the model with black box access to matrix entries.

We also note that for practical implementations, the constant runtime overhead with respect to a given set of elementary fault-tolerant quantum gates is an important consideration. It has been shown by Scherer et al. Scherer et al. (2017) that current approaches to the QLSA potentially suffer from a large constant overhead, hindering prospects of near-term applications. Whether our proposed QSVE-based algorithm exhibits a more advantageous constant overhead due to the absence of Hamiltonian simulation, remains an open question.

Acknowledgements.

The authors thank Simon Benjamin, Joseph Fitzsimons, Patrick Rebentrost and Nathan Wiebe for helpful comments on the manuscript, and Andrew Childs for his feedback on the earlier version. The authors also acknowledge support from Singapore’s Ministry of Education and National Research Foundation. This material is based on research funded in part by the Singapore National Research Foundation under NRF Award NRF-NRFF2013-01.

I Appendix

In this appendix we establish error bounds on the final state. The analysis below closely follows the analysis of the QLSA algorithm Harrow et al. (2009).

We use the filter functions $\mathrm{f}$ and $\mathrm{g}$  Hansen (1998), which allow us to invert only the the well-conditioned part of the matrix, i.e. the space which is spanned by the eigenspaces with eigenvalues, $\lambda_{i}\geq 1/\kappa$. We define the function $\mathrm{f}(\lambda):=1/(\gamma\kappa\lambda)$ for $|\lambda|\geq 1/\kappa$, $\mathrm{f}(\lambda):=0$ for $\lambda\leq 1/2\kappa$, and $\mathrm{f}(\lambda)$ is a smooth interpolating function $\eta_{1}(\lambda)$ for $1/2\kappa\leq\lambda\leq 1/\kappa$. Similarly, we define $\mathrm{g}(\lambda):=0$ for $|\lambda|\geq 1/\kappa$, $\mathrm{g}(\lambda):=1/2$ for $\lambda\leq 1/2\kappa$, and $\mathrm{g}(\lambda)$ is an interpolating function $\eta_{2}(\lambda)$ for $1/2\kappa\leq\lambda\leq 1/\kappa$. The interpolating functions $\eta_{1},\eta_{2}$ are chosen such that $f^{2}(\lambda)+g^{2}(\lambda)\leq 1$ for all $\lambda\in\mathbb{R}$. The algorithm in the main text corresponds to the choice $\mathrm{g}(\lambda)=0$.

Let $\gamma=\mathcal{O}\left(1/\kappa\right)$ be some constant which assures that the controlled rotation angle is less than $2\pi$ for any eigenvalues. We then define the map

[TABLE]

with $\mathrm{f}^{2}(x)+\mathrm{g}^{2}(x)\leq 1$, where $\ket{\mathrm{NO}}$ indicates that no matrix inversion has taken place, $\ket{\mathrm{IC}}$ means that part of $\ket{\mathbf{b}}$ is in the ill-conditioned subspace of $A$, and $\ket{\mathrm{WC}}$ means that the matrix inversion has taken place and is in the well conditioned subspace of $A$. This allows us to invert only the well conditioned part of the matrix while it flags the ill conditioned ones and interpolates between those two behaviours when $1/(2\kappa)<|\lambda|<1/\kappa$. We therefore only invert eigenvalues which are larger than $1/(2\kappa)$, which motivates the choice of $\mu$ in Algorithm 2.

Let $Q$ be the error-free operation corresponding to the QSVE subroutine followed by the controlled rotation without post-selection, i.e.

[TABLE]

$\overline{Q}$ in contrast describes the same procedure but the phase estimation step is erroneous, i.e.

[TABLE]

We want to bound the error in $\left\lVert{\overline{Q}-Q}\right\rVert$. By choosing a general state $\ket{\mathbf{b}}$, this is equivalent to the bound in $\left\lVert{Q\ket{\mathbf{b}}-\overline{Q}\ket{\mathbf{b}}}\right\rVert:=\left\lVert{\ket{\overline{\psi}}-\ket{\psi}}\right\rVert$. We will make use of the following lemma.

Lemma 3** (Harrow et al. (2009)).**

The map $\lambda\rightarrow\ket{h(\lambda)}$ is $\mathcal{O}\left(\kappa\right)$-Lipschitz, i.e. $\forall\lambda_{i}\neq\lambda_{j}$:

[TABLE]

for some $c\leq\pi/2=\mathcal{O}\left(1\right)$.

As $\left\lVert{\ket{\overline{\psi}}-\ket{\psi}}\right\rVert=\sqrt{2\left(1-Re\left<\overline{\psi}\left|\psi\right.\right>\right)},$ it suffices to lower-bound $Re\left<\overline{\psi}\left|\psi\right.\right>$:

[TABLE]

where we used the error bounds of the QSVE subroutine for the eigenvalue distance, i.e. $|\lambda_{i}-\overline{\lambda}_{i}|\leq\delta\left\lVert{A}\right\rVert_{F}$, which is a consequence of the phase estimation accuracy, and the $\mathcal{O}\left(\kappa\right)$-Lipschitz property in Lemma (3). Since $0\leq Re\left<\overline{\psi}\left|\psi\right.\right>\leq 1$, it follows that

[TABLE]

Using $\sum_{i}|\beta_{i}|^{2}=1$, the distance can be bounded as

[TABLE]

If we require this error to be of $\mathcal{O}\left(\epsilon\right)$, we need to take the phase estimation accuracy to be $\delta=\mathcal{O}\left(\frac{\epsilon}{\kappa\left\lVert{A}\right\rVert_{F}}\right)$. This results in a runtime $\mathcal{O}\left(\kappa\left\lVert{A}\right\rVert_{F}\cdot\text{polylog}(n)/\epsilon\right)$. In order to successfully perform the post-selection step, we need to repeat the algorithm on average $\kappa^{2}$ times. This additional multiplicative factor of $\kappa^{2}$ can be reduced to $\kappa$ using amplitude amplification Brassard et al. (2002). Putting everything together, we have a final runtime of $\mathcal{O}\left(\kappa^{2}\left\lVert{A}\right\rVert_{F}\cdot\text{polylog}(n)/\epsilon\right)$.