Quantum algorithm for linear differential equations with exponentially improved dependence on precision

TL;DR

This paper introduces a quantum algorithm for solving linear differential equations that significantly reduces complexity dependence on precision, leveraging recent quantum linear systems techniques.

Contribution

It presents a novel quantum algorithm with exponential improvement in precision dependence for solving linear differential equations.

Findings

Complexity is polynomial in log(1/error)

No additional stability hypotheses needed

Builds on recent quantum linear systems advances

Abstract

We present a quantum algorithm for systems of (possibly inhomogeneous) linear ordinary differential equations with constant coefficients. The algorithm produces a quantum state that is proportional to the solution at a desired final time. The complexity of the algorithm is polynomial in the logarithm of the inverse error, an exponential improvement over previous quantum algorithms for this problem. Our result builds upon recent advances in quantum linear systems algorithms by encoding the simulation into a sparse, well-conditioned linear system that approximates evolution according to the propagator using a Taylor series. Unlike with finite difference methods, our approach does not require additional hypotheses to ensure numerical stability.

Equations

$$\frac{d\vec{x}}{dt}=A\vec{x}+\vec{b}$$

$$\vec{x}(t)=\operatorname{exp}\mathopen{}\mathclose{{}\left(At}\right)\vec{x}(0)+(\operatorname{exp}\mathopen{}\mathclose{{}\left(At}\right)-I)A^{-1}\vec{b}.$$

$$T_{k}(z)\vcentcolon=\displaystyle\sum\limits_{j=0}^{k}\dfrac{z^{j}}{j!}\approx\operatorname{exp}\mathopen{}\mathclose{{}\left(z}\right)$$

$$S_{k}(z)\vcentcolon=\displaystyle\sum_{j=1}^{k}\frac{z^{j-1}}{j!}\approx(\operatorname{exp}\mathopen{}\mathclose{{}\left(z}\right)-1)z^{-1}$$

$$\vec{x}(h)\approx T_k(Ah)\vec{x}(0)+S_k(Ah)h\vec{b}.$$

$$C_{m,k,p}(A)$$

$$\displaystyle\quad-\displaystyle\sum\limits_{i=0}^{m-1}\displaystyle\sum\limits_{j=0}^{k}\mathopen{}\mathclose{{}\left|(i+1)(k+1)\rangle\langle i(k+1)+j}\right|\otimes I-\displaystyle\sum\limits_{j=d-p+1}^{d}\mathopen{}\mathclose{{}\left|j\rangle\langle j-1}\right|\otimes I,$$

$$C_{m,k,p}(Ah)\mathopen{}\mathclose{{}\left|x}\right\rangle=\mathopen{}\mathclose{{}\left|0}\right\rangle\mathopen{}\mathclose{{}\left|x_{\mathrm{in}}}\right\rangle+h\displaystyle\sum_{i=0}^{m-1}|i(k+1)+1\rangle\mathopen{}\mathclose{{}\left|b}\right\rangle,$$

$$C_{2,3,2}(Ah)|x⟩=\text{setcounter}MaxMatrixCols11\left(\begin{array}{ccccccccccc}I& & & & & & & & & & \\ -Ah& I& & & & & & & & & \\ & -Ah/2& I& & & & & & & & \\ & & -Ah/3& I& & & & & & & \\ -I& -I& -I& -I& I& & & & & & \\ & & & & -Ah& I& & & & & \\ & & & & & -Ah/2& I& & & & \\ & & & & & & -Ah/3& I& & & \\ & & & & -I& -I& -I& -I& I& & \\ & & & & & & & & -I& I& \\ & & & & & & & & & -I& I\end{array}\right)|x⟩=\left(\begin{array}{c}|x_{\mathrm{in}}⟩\\ h|b⟩\\ 0\\ 0\\ 0\\ h|b⟩\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right).$$

$$\mathopen{}\mathclose{{}\left|x}\right\rangle=C_{m,k,p}(Ah)^{-1}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left|0}\right\rangle\mathopen{}\mathclose{{}\left|x_{\mathrm{in}}}\right\rangle+h\displaystyle\sum_{i=0}^{m-1}|i(k+1)+1\rangle\mathopen{}\mathclose{{}\left|b}\right\rangle}\right],$$

$$\mathopen{}\mathclose{{}\left|x}\right\rangle=\displaystyle\sum\limits_{i=0}^{m-1}\sum\limits_{j=0}^{k}\mathopen{}\mathclose{{}\left|i(k+1)+j}\right\rangle\mathopen{}\mathclose{{}\left|x_{i,j}}\right\rangle+\displaystyle\sum\limits_{j=0}^{p}\mathopen{}\mathclose{{}\left|m(k+1)+j}\right\rangle\mathopen{}\mathclose{{}\left|x_{m,j}}\right\rangle$$

$$\displaystyle\mathopen{}\mathclose{{}\left|x_{0,0}}\right\rangle$$

$$\displaystyle\mathopen{}\mathclose{{}\left|x_{i,0}}\right\rangle$$

$$\displaystyle\mathopen{}\mathclose{{}\left|x_{i,1}}\right\rangle$$

$$\displaystyle\mathopen{}\mathclose{{}\left|x_{i,j}}\right\rangle$$

$$\displaystyle\mathopen{}\mathclose{{}\left|x_{m,j}}\right\rangle$$

$$\displaystyle\mathopen{}\mathclose{{}\left|x_{0,0}}\right\rangle$$

$$\displaystyle\mathopen{}\mathclose{{}\left|x_{0,j}}\right\rangle$$

$$\displaystyle\mathopen{}\mathclose{{}\left|x_{1,0}}\right\rangle$$

$$\displaystyle\approx\operatorname{exp}\mathopen{}\mathclose{{}\left(Ah}\right)\mathopen{}\mathclose{{}\left|x_{\mathrm{in}}}\right\rangle+(\operatorname{exp}\mathopen{}\mathclose{{}\left(Ah}\right)-I)A^{-1}\mathopen{}\mathclose{{}\left|b}\right\rangle,$$

$$\displaystyle\mathopen{}\mathclose{{}\left|x_{1,j}}\right\rangle$$

$$\displaystyle\mathopen{}\mathclose{{}\left|x_{2,0}}\right\rangle$$

$$\displaystyle\approx\operatorname{exp}\mathopen{}\mathclose{{}\left(2Ah}\right)\mathopen{}\mathclose{{}\left|x_{\mathrm{in}}}\right\rangle+(\operatorname{exp}\mathopen{}\mathclose{{}\left(2Ah}\right)-I)A^{-1}\mathopen{}\mathclose{{}\left|b}\right\rangle,$$

$$⋮$$

$$\displaystyle\mathopen{}\mathclose{{}\left|x_{m-1,0}}\right\rangle$$

$$\displaystyle\approx\operatorname{exp}\mathopen{}\mathclose{{}\left(Ah(m-1)}\right)\mathopen{}\mathclose{{}\left|x_{\mathrm{in}}}\right\rangle+(\operatorname{exp}\mathopen{}\mathclose{{}\left(Ah(m-1)}\right)-I)A^{-1}\mathopen{}\mathclose{{}\left|b}\right\rangle,$$

$$\displaystyle\mathopen{}\mathclose{{}\left|x_{m-1,j}}\right\rangle$$

$$\displaystyle\mathopen{}\mathclose{{}\left|x_{m,0}}\right\rangle$$

$$\displaystyle\approx\operatorname{exp}\mathopen{}\mathclose{{}\left(Ahm}\right)\mathopen{}\mathclose{{}\left|x_{\mathrm{in}}}\right\rangle+(\operatorname{exp}\mathopen{}\mathclose{{}\left(Ahm}\right)-I)A^{-1}\mathopen{}\mathclose{{}\left|b}\right\rangle,$$

$$\displaystyle\mathopen{}\mathclose{{}\left|x_{m,j}}\right\rangle$$

$$\displaystyle\approx\operatorname{exp}\mathopen{}\mathclose{{}\left(Ahm}\right)\mathopen{}\mathclose{{}\left|x_{\mathrm{in}}}\right\rangle+(\operatorname{exp}\mathopen{}\mathclose{{}\left(Ahm}\right)-I)A^{-1}\mathopen{}\mathclose{{}\left|b}\right\rangle,$$

$$\frac{d\vec{x}}{dt}=A\vec{x}+\vec{b}$$

$$\vec{x}(t)=\operatorname{exp}\mathopen{}\mathclose{{}\left(At}\right)\vec{x}_{\mathrm{in}}+(\operatorname{exp}\mathopen{}\mathclose{{}\left(At}\right)-I)A^{-1}\vec{b}$$

$$\mathopen{}\mathclose{{}\left\|C_{m,k,p}(\lambda)^{-1}\mathopen{}\mathclose{{}\left|l}\right\rangle}\right\|\leq\sqrt{1.04eI_{0}(2)(m+p)}$$

$$\mathopen{}\mathclose{{}\left|\mathopen{}\mathclose{{}\left\langle n}\right|C_{m,k,p}(\lambda)^{-1}\mathopen{}\mathclose{{}\left|l}\right\rangle}\right|\leq\sqrt{1.04e}.$$

$$T_{b,k}(z):=\sum _{j=b}^k\frac{b!z^{j-b}}{j!}$$

$$C_{m,k,p}(\lambda)\mathopen{}\mathclose{{}\left|x}\right\rangle=\mathopen{}\mathclose{{}\left|l}\right\rangle$$

$$\mathopen{}\mathclose{{}\left|x}\right\rangle=\displaystyle\sum\limits_{i=0}^{m-1}\sum\limits_{j=0}^{k}x_{i,j}\mathopen{}\mathclose{{}\left|i(k+1)+j}\right\rangle+\displaystyle\sum\limits_{j=0}^{p}x_{m,j}\mathopen{}\mathclose{{}\left|m(k+1)+j}\right\rangle$$

$$x_{i,0}-\sum _{j=0}^k{x}_{i-1,j}$$

$$x_{i,j}-(\lambda /j)x_{i,j-1}$$

Full text

1 Introduction

One of the original motivations for developing a quantum computer was to efficiently simulate Hamiltonian dynamics, i.e., differential equations of the form $\frac{d\vec{x}}{dt}=A\vec{x}$ where $A$ is anti-Hermitian. Given a suitable description of $A$, a copy of the initial quantum state \mathopen{}\mathclose{{}\left|x(0)}\right\rangle, and an evolution time $T$, the goal is to produce a quantum state that is $\epsilon$-close to the final state \mathopen{}\mathclose{{}\left|x(T)}\right\rangle. The first algorithms for this problem had complexity polynomial in $1/\epsilon$ [13, 2, 10, 4]. Subsequent work gave an algorithm with complexity \operatorname{poly}\mathopen{}\mathclose{{}\left(\operatorname{log}\mathopen{}\mathclose{{}\left(1/\epsilon}\right)}\right)—an exponential improvement—which is optimal in a black-box model [5]. More recent work has streamlined these algorithms and improved their dependence on other parameters [6, 7, 15, 8, 14, 16].

While Hamiltonian simulation has been a focus of quantum algorithms research, the more general problem of simulating linear differential equations of the form $\frac{d\vec{x}}{dt}=A\vec{x}+\vec{b}$ for arbitrary (sparse) $A$ is less well studied. Reference [3] solves this problem using a quantum linear systems algorithm (QLSA) to implement linear multistep methods, which represent the differential equations with a system of linear equations by discretizing time. The complexity of this approach is \operatorname{poly}\mathopen{}\mathclose{{}\left(1/\epsilon}\right). Considering the recent improvements to the complexity of Hamiltonian simulation, it is natural to ask whether linear differential equations can be solved more efficiently as a function of $\epsilon$.

Hamiltonian simulation is a central component of the QLSA, and the techniques underlying \operatorname{poly}\mathopen{}\mathclose{{}\left(\operatorname{log}\mathopen{}\mathclose{{}\left(1/\epsilon}\right)}\right) Hamiltonian simulation have been adapted to give a QLSA with complexity \operatorname{poly}\mathopen{}\mathclose{{}\left(\operatorname{log}\mathopen{}\mathclose{{}\left(1/\epsilon}\right)}\right) [11]. However, even if this improved QLSA is used to implement the algorithm of Ref. [3], the overall complexity is still \operatorname{poly}\mathopen{}\mathclose{{}\left(1/\epsilon}\right), since the multistep method itself is a significant source of error.

In a similar vein, the QLSA of Ref. [12] has \operatorname{poly}\mathopen{}\mathclose{{}\left(1/\epsilon}\right) complexity even when using a Hamiltonian simulation algorithm with \operatorname{poly}\mathopen{}\mathclose{{}\left(\operatorname{log}\mathopen{}\mathclose{{}\left(1/\epsilon}\right)}\right) complexity, simply because phase estimation has complexity \operatorname{poly}\mathopen{}\mathclose{{}\left(1/\epsilon}\right). Reference [11] provides a QLSA with \operatorname{poly}\mathopen{}\mathclose{{}\left(\operatorname{log}\mathopen{}\mathclose{{}\left(1/\epsilon}\right)}\right) complexity by avoiding phase estimation and instead directly inverting the linear system using a linear combination of unitaries (LCU). Thus, one might consider realizing the solution of $\frac{d\vec{x}}{dt}=A\vec{x}+\vec{b}$ as a linear combination of unitaries. Unfortunately, in the general case where $A$ is not anti-Hermitian, the best implementation of this approach that we are aware of has an exponentially small success probability.

In this paper, we circumvent these limitations and present a quantum algorithm for linear differential equations with complexity \operatorname{poly}\mathopen{}\mathclose{{}\left(\operatorname{log}\mathopen{}\mathclose{{}\left(1/\epsilon}\right)}\right), an exponential improvement over Ref. [3]. As in Ref. [3], our approach applies the QLSA. However, instead of using a linear multistep method, we encode a truncation of the Taylor series of \operatorname{exp}\mathopen{}\mathclose{{}\left(At}\right), the propagator for the differential equation, into a linear system. Since it effectively implements a linear combination of operations, our approach is conceptually similar to quantum simulation via linear combinations of unitaries, but we achieve significantly better performance by constructing this linear combination stepwise through a system of linear equations. This alternative to direct application of LCU methods might be advantageous for other quantum algorithms.

In addition to scaling well with the simulation error, our algorithm has favorable performance as a function of other parameters. The complexity is nearly linear in the evolution time, which is a quadratic improvement over Ref. [3] and is nearly optimal [4]. The complexity is also nearly linear in the sparsity of $A$ and in a parameter characterizing the decay of the solution vector. The latter dependence is necessary since producing a normalized version of a subnormalized solution vector is equivalent to postselection, which is computationally intractable [1], as discussed further in Section 8. Along similar lines, we assume that the eigenvalues of $A$ have non-positive real part since it is intractable to simulate exponentially growing solutions. (This improves upon Ref. [3], where the eigenvalues $\lambda$ of $A$ must satisfy $|\arg(-\lambda)|\leq\alpha$ for some constant $\alpha$ depending on the stability of the multistep method.) For a precise statement of the main result, see Theorem 9.

This paper is organized as follows. Section 2 describes how we encode the solution of a system of differential equations into a system of linear equations. The following three sections analyze properties of this system: Section 3 bounds its condition number, Section 4 analyzes how well it approximates the differential equation, and Section 5 shows that a measurement of its solution vector provides a solution of the differential equation with appreciable probability. In Section 6 we explain how to prepare the state that is input to the QLSA using black boxes for the initial condition and inhomogeneous term of the differential equation. We formally state and prove our main result in Section 7. Finally, we conclude in Section 8 with a discussion of the result and some open problems.

2 Constructing the Linear System

As in Ref. [3] we consider a differential equation of the form

[TABLE]

where $A$ and $\vec{b}$ are time-independent. This has the exact solution

[TABLE]

Define

[TABLE]

and

[TABLE]

where the approximations hold for large $k$. Then for short evolution time $h$ (namely h\leq 1/\mathopen{}\mathclose{{}\left\|A}\right\|, where \mathopen{}\mathclose{{}\left\|\cdot}\right\| denotes the spectral norm) and large $k$, we can approximate the solution by

[TABLE]

This approximate solution can be used in turn as an initial condition for another step of evolution, and we can repeat this procedure as desired for a total number of steps $m$.

We encode this procedure in a linear system using the following family of matrices.

Definition 1**.**

Let $A$ be an $N\times N$ matrix, and let $m,k,p\in\mathbb{Z}^{+}$. Define

[TABLE]

where $d\vcentcolon=m(k+1)+p$, and $I$ is the $N\times N$ identity matrix.

Now consider the linear system

[TABLE]

where \mathopen{}\mathclose{{}\left|x_{\mathrm{in}}}\right\rangle,\mathopen{}\mathclose{{}\left|b}\right\rangle\in\mathbb{C}^{N} and $h\in\mathbb{R}^{+}$. The first register labels a natural block structure for $C_{m,k,p}$. For example, the system C_{2,3,2}(Ah)|x\rangle=\mathopen{}\mathclose{{}\left|0}\right\rangle\mathopen{}\mathclose{{}\left|x_{\mathrm{in}}}\right\rangle+h\sum_{i=0}^{1}|4i+1\rangle\mathopen{}\mathclose{{}\left|b}\right\rangle is as follows:

[TABLE]

After performing $m$ steps of the evolution approximated with a Taylor series of order $k$, the solution is kept constant for $p$ steps. This ensures a significant probability of obtaining the solution at the final time, similarly as in Ref. [3]. Note that $C_{m,k,p}(Ah)$ is nonsingular, since it is a lower-triangular matrix with nonzero diagonal entries. If $N=1$, then $C_{m,k,p}(Ah)$ is a $(d+1)\times(d+1)$ matrix.

The solution of Eq. (7) is

[TABLE]

which can be written as

[TABLE]

for some \mathopen{}\mathclose{{}\left|x_{i,j}}\right\rangle\in\mathbb{C}^{N}. By the definition of $C_{m,k,p}(Ah)$, these \mathopen{}\mathclose{{}\left|x_{i,j}}\right\rangles satisfy

[TABLE]

From these equations, we obtain

[TABLE]

In these approximations, we assume $k$ is sufficiently large that we can neglect the truncation errors \mathopen{}\mathclose{{}\left\|T_{k}(Ah)-\operatorname{exp}\mathopen{}\mathclose{{}\left(Ah}\right)}\right\| and \mathopen{}\mathclose{{}\left\|S_{k}(Ah)h-(\operatorname{exp}\mathopen{}\mathclose{{}\left(Ah}\right)-I)A^{-1}}\right\| (we make this more precise in Section 4). Note that \mathopen{}\mathclose{{}\left|x}\right\rangle (defined by Eq. (10)) includes a piece that can be interpreted as the history state of the evolution

[TABLE]

with the initial condition $\vec{x}(0)=\vec{x}_{\mathrm{in}}$. More precisely, \mathopen{}\mathclose{{}\left|x_{i,0}}\right\rangle is a good approximation of the system’s state at time $ih$, for any $i\in\{0,1,\dots,m\}$. Furthermore, \mathopen{}\mathclose{{}\left|x_{m,0}}\right\rangle=\mathopen{}\mathclose{{}\left|x_{m,1}}\right\rangle=\dots=\mathopen{}\mathclose{{}\left|x_{m,p}}\right\rangle is a good approximation of

[TABLE]

for $t=mh$.

If we measured the first register of the state \mathopen{}\mathclose{{}\left|x}\right\rangle/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x}\right\rangle}\right\| in the standard basis, we would obtain the state \mathopen{}\mathclose{{}\left|x_{i,j}}\right\rangle/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{i,j}}\right\rangle}\right\| for random $i,j$. By choosing a large $p$, we can ensure there is a high probability of obtaining \mathopen{}\mathclose{{}\left|x_{m,0}}\right\rangle/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{m,0}}\right\rangle}\right\|, \mathopen{}\mathclose{{}\left|x_{m,1}}\right\rangle/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{m,1}}\right\rangle}\right\|, $\dots$, or \mathopen{}\mathclose{{}\left|x_{m,p}}\right\rangle/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{m,p}}\right\rangle}\right\| (as we show in Section 5). Then this probability can be raised to \Omega\mathopen{}\mathclose{{}\left(1}\right) by using amplitude amplification (or by classical repetition). This is how we prepare a state close to \vec{x}(t)/\mathopen{}\mathclose{{}\left\|\vec{x}(t)}\right\| for $t=mh$.

Note that the state we generate is not of the same form as the history state in Ref. [3], which only encodes $\vec{x}(t)$ at intermediate times. The solution of our linear system not only encodes $\vec{x}(t)$ at intermediate times (via $|x_{i,0}\rangle$) but also encodes $((Ah)^{j}/j!)\vec{x}(t)+((Ah)^{j-1}/j!)h\vec{b}$ at intermediate times (via $|x_{i,j}\rangle$).

To analyze the performance of this approach to solving differential equations, we establish three properties of this system of linear equations. First, since the complexity of the best known QLSAs grows linearly with condition number (and sublinear complexity is impossible unless $\textsf{BQP}=\textsf{PSPACE}$ [12]), we analyze the condition number of $C_{m,k,p}$ (Section 3). Second, we show that the solution of the linear system includes a piece that is close to the solution of the associated differential equation (Section 4). Third, we show that this piece can be obtained from a measurement that succeeds with appreciable probability (Section 5).

3 Condition Number

In this section, we upper bound the condition number of the matrix $C_{m,k,p}(A)$ under mild assumptions about $A$. We begin with a technical lemma that upper bounds the norms of the columns of the inverse of this matrix.

Lemma 2**.**

Let $\lambda\in\mathbb{C}$ such that \mathopen{}\mathclose{{}\left|\lambda}\right|\leq 1 and $\mathrm{Re}(\lambda)\leq 0$. Let $m,k,p\in\mathbb{Z}^{+}$ such that $k\geq 5$ and $(k+1)!\geq 2m$, and let $d=m(k+1)+p$. Then for any $n,l\in\{0,1,\dots,d\}$,

[TABLE]

with $I_{0}(2)<2.28$ a modified Bessel function of the first kind, and

[TABLE]

Proof.

Recall the definitions of $T_{k}(z)$ and $S_{k}(z)$ in Eqs. (3) and (4), respectively. We also define

[TABLE]

for $b\leq k$.

Fix any $l\in\{0,1,\dots,d\}$. Suppose the solution of the linear system

[TABLE]

is

[TABLE]

for some $x_{i,j}\in\mathbb{C}$. By the definition of $C_{m,k,p}(\lambda)$, the $x_{i,j}$s should satisfy

[TABLE]

where $\delta_{i,j}=1$ if $i=j$, and [math] otherwise.

We consider the cases $0\leq l<m(k+1)$ and $m(k+1)\leq l\leq d$ separately.

• •

Case 1: $0\leq l<m(k+1)$. Suppose $l=a(k+1)+b$ for some $0\leq a<m$ and $0\leq b\leq k$. In this case, Eq. (34) implies

[TABLE]

Since \mathopen{}\mathclose{{}\left|\lambda}\right|\leq 1, for any $b\leq j\leq k$, we have

[TABLE]

Furthermore, since \mathopen{}\mathclose{{}\left|\lambda}\right|\leq 1 and $\mathrm{Re}(\lambda)\leq 0$, by Lemma 10 and Lemma 11 in Appendix A, we have

[TABLE]

[TABLE]

Consequently, we have

[TABLE]

and

[TABLE]

It follows that

[TABLE]

Using these facts, we obtain

[TABLE]

where in the fifth step we use the facts

[TABLE]

• •

Case 2: $m(k+1)\leq l\leq d$. Suppose $l=m(k+1)+b$ for some $0\leq b\leq p$. In this case, Eq. (34) implies

[TABLE]

It follows that

[TABLE]

In both of the above cases, we have \mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x}\right\rangle}\right\|\leq\sqrt{1.04eI_{0}(2)(m+p)} and \mathopen{}\mathclose{{}\left|\mathopen{}\mathclose{{}\left\langle n|x}\right\rangle}\right|\leq\sqrt{1.04e} for any $n\in\{0,1,\dots,d\}$, as claimed. ∎

Now we are ready to upper bound the norm of the inverse of the matrix.

Lemma 3**.**

Let $A=VDV^{-1}$ be a diagonalizable matrix, where D=\operatorname{diag}\mathopen{}\mathclose{{}\left(\lambda_{0},\lambda_{1},\dots,\lambda_{N-1}}\right) satisfies \mathopen{}\mathclose{{}\left|\lambda_{i}}\right|\leq 1 and $\mathrm{Re}(\lambda_{i})\leq 0$ for $i\in\{0,1,\dots,N-1\}$. Let $m,k,p\in\mathbb{Z}^{+}$ such that $k\geq 5$ and $(k+1)!\geq 2m$. Then

[TABLE]

where \kappa_{V}=\mathopen{}\mathclose{{}\left\|V}\right\|\cdot\mathopen{}\mathclose{{}\left\|V^{-1}}\right\| is the condition number of $V$.

Proof.

For convenience, we will drop the subscripts $m,k,p$, and use $C(\cdot)$ to denote $C_{m,k,p}(\cdot)$. We diagonalize $C(A)$ as

[TABLE]

where \tilde{V}\vcentcolon=\sum_{j=0}^{d}\mathopen{}\mathclose{{}\left|j\rangle\langle j}\right|\otimes V has condition number $\kappa_{\tilde{V}}=\kappa_{V}$. Then we may upper bound \mathopen{}\mathclose{{}\left\|C(A)^{-1}}\right\| in terms of \mathopen{}\mathclose{{}\left\|C(D)^{-1}}\right\| as

[TABLE]

We have

[TABLE]

where we maximize over all states \mathopen{}\mathclose{{}\left|\psi}\right\rangle\in\mathbb{C}^{(d+1)N} The state \mathopen{}\mathclose{{}\left|\psi}\right\rangle can be written as \mathopen{}\mathclose{{}\left|\psi}\right\rangle=\sum_{l=0}^{d}\mathopen{}\mathclose{{}\left|l}\right\rangle\mathopen{}\mathclose{{}\left|\psi_{l}}\right\rangle for some \mathopen{}\mathclose{{}\left|\psi_{l}}\right\rangle\in\mathbb{C}^{N}. Then we have

[TABLE]

Now let \mathopen{}\mathclose{{}\left|\psi_{l}}\right\rangle=\sum_{j=0}^{N-1}\psi_{j,l}\mathopen{}\mathclose{{}\left|j}\right\rangle for some $\psi_{j,l}\in\mathbb{C}$. With D=\sum_{j=0}^{N-1}\lambda_{j}\mathopen{}\mathclose{{}\left|j\rangle\langle j}\right|, we have C(D)=\sum_{j=0}^{N-1}C(\lambda_{j})\otimes\mathopen{}\mathclose{{}\left|j\rangle\langle j}\right|. Hence, by Lemma 28, we obtain

[TABLE]

Using this expression in Eq. (60) gives

[TABLE]

Using this result in Eq. (59) yields

[TABLE]

Combining this expression with Eq. (58) then gives Eq. (56), as claimed. ∎

It remains to upper bound the norm of the matrix.

Lemma 4**.**

Let $A$ be an $N\times N$ matrix such that \mathopen{}\mathclose{{}\left\|A}\right\|\leq 1. Let $m,k,p\in\mathbb{Z}^{+}$, and $k\geq 5$. Then

[TABLE]

Proof.

Observe that $C\vcentcolon=C_{m,k,p}(A)$ can be written as the sum of three matrices:

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

where $d=m(k+1)+p$. One can easily check that \mathopen{}\mathclose{{}\left\|C_{1}}\right\|=1, \mathopen{}\mathclose{{}\left\|C_{2}}\right\|=\sqrt{k+1}, and \mathopen{}\mathclose{{}\left\|C_{3}}\right\|=\max\{\mathopen{}\mathclose{{}\left\|A}\right\|,1\}=1 (this is trivial for $C_{1}$, and follows directly from a calculation of $C_{2}C_{2}^{\dagger}$ and $C_{3}C_{3}^{\dagger}$ for the other cases). Consequently,

[TABLE]

as claimed. ∎

Combining Lemma 3 and Lemma 64, we obtain the following upper bound on the condition number of $C_{m,k,p}(A)$:

Theorem 5**.**

Let $A=VDV^{-1}$ be a diagonalizable matrix such that \mathopen{}\mathclose{{}\left\|A}\right\|\leq 1, D=\operatorname{diag}\mathopen{}\mathclose{{}\left(\lambda_{0},\lambda_{1},\dots,\lambda_{N-1}}\right) and $\mathrm{Re}(\lambda_{i})\leq 0$, for $i\in\{0,1,\dots,N-1\}$. Let $m,k,p\in\mathbb{Z}^{+}$ such that $k\geq 5$ and $(k+1)!\geq 2m$. Let $C\vcentcolon=C_{m,k,p}(A)$, and let \kappa_{C}=\mathopen{}\mathclose{{}\left\|C}\right\|\cdot\mathopen{}\mathclose{{}\left\|C^{-1}}\right\| be the condition number of $C$. Then

[TABLE]

where \kappa_{V}=\mathopen{}\mathclose{{}\left\|V}\right\|\cdot\mathopen{}\mathclose{{}\left\|V^{-1}}\right\| is the condition number of $V$.

4 Solution Error

In this section, we prove that the solution of the linear system defined by Eq. (7) encodes a good approximation of the solution of the differential equation defined by Eq. (25) with the initial condition $\vec{x}(0)=\vec{x}_{\mathrm{in}}$.

Theorem 6**.**

Let $A=VDV^{-1}$ be a diagonalizable matrix, where D=\operatorname{diag}\mathopen{}\mathclose{{}\left(\lambda_{0},\lambda_{1},\dots,\lambda_{N-1}}\right) satisfies $\mathrm{Re}(\lambda_{i})\leq 0$ for $i\in\{0,1,\dots,N-1\}$. Let $h\in\mathbb{R}^{+}$ such that \mathopen{}\mathclose{{}\left\|Ah}\right\|\leq 1. Let \mathopen{}\mathclose{{}\left|{x}_{\mathrm{in}}}\right\rangle,\mathopen{}\mathclose{{}\left|b}\right\rangle\in\mathbb{C}^{N}, and let \mathopen{}\mathclose{{}\left|x(t)}\right\rangle be defined by Eq. (26). Let $m,k,p\in\mathbb{Z}^{+}$ such that $k\geq 5$ and $(k+1)!\geq 2m$. Let \mathopen{}\mathclose{{}\left|x_{i,j}}\right\rangle be defined by Eqs. (9) and (10). Then for any $j\in\{0,1,\dots,m\}$,

[TABLE]

where \kappa_{V}=\mathopen{}\mathclose{{}\left\|V}\right\|\cdot\mathopen{}\mathclose{{}\left\|V^{-1}}\right\| is the condition number of $V$.

Proof.

Note that \mathopen{}\mathclose{{}\left|x(jh)}\right\rangle (the solution of the differential equation) satisfies the recurrence relation

[TABLE]

while \mathopen{}\mathclose{{}\left|x_{j,0}}\right\rangle (in the solution of the associated linear system) satisfies the recurrence relation

[TABLE]

Recall that $T_{k}(\lambda)=\sum_{j=0}^{k}\frac{\lambda^{j}}{j!}$ and $S_{k}(\lambda)=\sum_{j=1}^{k}\frac{\lambda^{j-1}}{j!}$. In addtion, we have \mathopen{}\mathclose{{}\left|x(0)}\right\rangle=\mathopen{}\mathclose{{}\left|x_{0,0}}\right\rangle=\mathopen{}\mathclose{{}\left|x_{\mathrm{in}}}\right\rangle.

Define \mathopen{}\mathclose{{}\left|y(t)}\right\rangle\vcentcolon=V^{-1}\mathopen{}\mathclose{{}\left|x(t)}\right\rangle and \mathopen{}\mathclose{{}\left|y_{i,j}}\right\rangle=V^{-1}\mathopen{}\mathclose{{}\left|x_{i,j}}\right\rangle. We will give an upper bound on \delta_{j}\vcentcolon=\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|y(jh)}\right\rangle-\mathopen{}\mathclose{{}\left|y_{j,0}}\right\rangle}\right\| and convert it into an upper bound on \epsilon_{j}\vcentcolon=\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x(jh)}\right\rangle-\mathopen{}\mathclose{{}\left|x_{j,0}}\right\rangle}\right\|. Since $A=VDV^{-1}$, we have \operatorname{exp}\mathopen{}\mathclose{{}\left(Ah}\right)=V\operatorname{exp}\mathopen{}\mathclose{{}\left(Dh}\right)V^{-1}, $T_{k}(Ah)=VT_{k}(Dh)V^{-1}$, and $S_{k}(Ah)=VS_{k}(Dh)V^{-1}$. Then Eq. (72) implies

[TABLE]

where \mathopen{}\mathclose{{}\left|c}\right\rangle=V^{-1}\mathopen{}\mathclose{{}\left|b}\right\rangle. Meanwhile, Eq. (73) implies

[TABLE]

In addition, we have \mathopen{}\mathclose{{}\left|y(0)}\right\rangle=\mathopen{}\mathclose{{}\left|y_{0,0}}\right\rangle=\mathopen{}\mathclose{{}\left|y_{\mathrm{in}}}\right\rangle\vcentcolon=V^{-1}\mathopen{}\mathclose{{}\left|x_{\mathrm{in}}}\right\rangle.

Now since $\mathrm{Re}(\lambda_{i})\leq 0$ for any $i\in\{0,1,\dots,N-1\}$, we have \mathopen{}\mathclose{{}\left\|\operatorname{exp}\mathopen{}\mathclose{{}\left(Dh}\right)}\right\|\leq 1. Moreover, since \mathopen{}\mathclose{{}\left\|Ah}\right\|\leq 1, we have \mathopen{}\mathclose{{}\left|\lambda_{i}h}\right|\leq 1 for any $i\in\{0,1,\dots,N-1\}$. Then since $\mathrm{Re}(\lambda_{i}h)\leq 0$ for any $i\in\{0,1,\dots,N-1\}$, by Lemma 10 in Appendix A, we get

[TABLE]

and by Lemma 12 in Appendix A, we get

[TABLE]

The error \delta_{j}=\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|y(jh)}\right\rangle-\mathopen{}\mathclose{{}\left|y_{j,0}}\right\rangle}\right\| can be bounded as follows. Note that $\delta_{0}=0$, and using the triangle inequality, we have

[TABLE]

This implies

[TABLE]

Next, we give an upper bound on \max_{0\leq i\leq m}\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|y_{i,0}}\right\rangle}\right\|. Fix any $i\in\{0,1,\dots,m\}$. Then we have

[TABLE]

where the bra \mathopen{}\mathclose{{}\left\langle i(k+1)}\right| acts on the first register, and

[TABLE]

Therefore, by the triangle inequality,

[TABLE]

To bound both terms, we consider the general expression \mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left\langle i(k+1)}\right|C_{m,k,p}(D)^{-1}\mathopen{}\mathclose{{}\left|s}\right\rangle\mathopen{}\mathclose{{}\left|\beta}\right\rangle}\right\| for integer $s$. Let \mathopen{}\mathclose{{}\left|\beta}\right\rangle=\sum_{l=0}^{N-1}\beta_{l}\mathopen{}\mathclose{{}\left|l}\right\rangle for some $\beta_{l}\in\mathbb{C}$. Then since C_{m,k,p}(D)=\sum_{l=0}^{N-1}C_{m,k,p}(\lambda_{l})\otimes\mathopen{}\mathclose{{}\left|l\rangle\langle l}\right|, we get

[TABLE]

Then, by Lemma 28, we have

[TABLE]

Hence we obtain

[TABLE]

Using these two facts and the triangle inequality, Eq. (82) implies

[TABLE]

Since this holds for any $i\in\{0,1,\dots,m\}$, we get

[TABLE]

Now using Eqs. (79) and (88), we get

[TABLE]

Finally, recall that \mathopen{}\mathclose{{}\left|x(t)}\right\rangle=V\mathopen{}\mathclose{{}\left|y(t)}\right\rangle and \mathopen{}\mathclose{{}\left|x_{i,j}}\right\rangle=V\mathopen{}\mathclose{{}\left|y_{i,j}}\right\rangle. Thus we have

[TABLE]

for any $j\in\{0,1,\dots,m\}$, as claimed. ∎

Theorem 6 implies that by choosing (k+1)!\geq 3\kappa_{V}m\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{\mathrm{in}}}\right\rangle}\right\|+mh\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|b}\right\rangle}\right\|}\right)/\epsilon, we can ensure \mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{j,0}}\right\rangle-\mathopen{}\mathclose{{}\left|x(jh)}\right\rangle}\right\|\leq\epsilon for any $j\in\{0,1,\dots,m\}$, so that \mathopen{}\mathclose{{}\left|x_{j,0}}\right\rangle is close to \mathopen{}\mathclose{{}\left|x(jh)}\right\rangle for any $j$. Furthermore, by Lemma 142 in Appendix B, if \mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{j,0}}\right\rangle-\mathopen{}\mathclose{{}\left|x(jh)}\right\rangle}\right\|\leq\epsilon and \mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x(jh)}\right\rangle}\right\|\geq\alpha, then \mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{j,0}}\right\rangle/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{j,0}}\right\rangle}\right\|-\mathopen{}\mathclose{{}\left|x(jh)}\right\rangle/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x(jh)}\right\rangle}\right\|}\right\|\leq 2\epsilon/\alpha. So, provided that \mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x(jh)}\right\rangle}\right\| is large enough, the normalized state \mathopen{}\mathclose{{}\left|x_{j,0}}\right\rangle/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{j,0}}\right\rangle}\right\| is also close to the normalized state \mathopen{}\mathclose{{}\left|x(jh)}\right\rangle/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x(jh)}\right\rangle}\right\|.

5 Success Probability

In the previous section, we have shown that \mathopen{}\mathclose{{}\left|x_{m,0}}\right\rangle/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{m,0}}\right\rangle}\right\|=\mathopen{}\mathclose{{}\left|x_{m,1}}\right\rangle/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{m,1}}\right\rangle}\right\|=\dots=\mathopen{}\mathclose{{}\left|x_{m,p}}\right\rangle/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{m,p}}\right\rangle}\right\| is a good approximation of \mathopen{}\mathclose{{}\left|x(mh)}\right\rangle/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x(mh)}\right\rangle}\right\|, provided the truncation order $k$ is sufficiently large. In this section, we show that such a state can be obtained with non-negligible probability by measuring the first register of \mathopen{}\mathclose{{}\left|x}\right\rangle/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x}\right\rangle}\right\| in the standard basis, provided the padding parameter $p$ is sufficiently large.

Theorem 7**.**

Let $A=VDV^{-1}$ be a diagonalizable matrix, where D=\operatorname{diag}\mathopen{}\mathclose{{}\left(\lambda_{0},\lambda_{1},\dots,\lambda_{N-1}}\right) satisfies $\mathrm{Re}(\lambda_{i})\leq 0$ for $i\in\{0,1,\dots,N-1\}$. Let $h\in\mathbb{R}^{+}$ such that \mathopen{}\mathclose{{}\left\|Ah}\right\|\leq 1. Let \mathopen{}\mathclose{{}\left|{x}_{\mathrm{in}}}\right\rangle,\mathopen{}\mathclose{{}\left|b}\right\rangle\in\mathbb{C}^{N}, and let \mathopen{}\mathclose{{}\left|x(t)}\right\rangle be defined by Eq. (26). Let $m,k,p\in\mathbb{Z}^{+}$ such that (k+1)!\geq 70\kappa_{V}m(\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{\mathrm{in}}}\right\rangle}\right\|+mh\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|b}\right\rangle}\right\|)/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x(mh)}\right\rangle}\right\|, where \kappa_{V}=\mathopen{}\mathclose{{}\left\|V}\right\|\cdot\mathopen{}\mathclose{{}\left\|V^{-1}}\right\| is the condition number of $V$. Let g={\max_{t\in[0,mh]}\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x(t)}\right\rangle}\right\|}/{\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x(mh)}\right\rangle}\right\|}. Let \mathopen{}\mathclose{{}\left|x}\right\rangle be defined by Eq. (9) and let \mathopen{}\mathclose{{}\left|x_{i,j}}\right\rangle be defined by Eq. (10). Then for any $j\in\{0,1,\dots,p\}$,

[TABLE]

Proof.

Recall that \mathopen{}\mathclose{{}\left|x_{m,j}}\right\rangle=\mathopen{}\mathclose{{}\left|x_{m,0}}\right\rangle for all $j\in\{1,2,\dots,p\}$. Thus it is sufficient to prove Eq. (91) for $j=0$.

Define

[TABLE]

and

[TABLE]

Then \mathopen{}\mathclose{{}\left|x}\right\rangle=\mathopen{}\mathclose{{}\left|x_{\rm good}}\right\rangle+\mathopen{}\mathclose{{}\left|x_{\rm bad}}\right\rangle and \mathopen{}\mathclose{{}\left\langle x_{\rm good}|x_{\rm bad}}\right\rangle=0, so

[TABLE]

Next we give a lower bound on \mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{m,0}}\right\rangle}\right\| and an upper bound on \mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{\rm bad}}\right\rangle}\right\|. Let q=\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x(mh)}\right\rangle}\right\|. Then by the definition of $g$, we have \mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x(ih)}\right\rangle}\right\|\leq gq for $i\in\{0,1,\dots,m-1\}$. Meanwhile, by Theorem 6 and the choice of $k$, we have

[TABLE]

As a result, we get

[TABLE]

and

[TABLE]

Now for any $i\in\{0,1,\dots,m-1\}$, since \mathopen{}\mathclose{{}\left|x_{i,j}}\right\rangle=({Ah}/{j})\mathopen{}\mathclose{{}\left|x_{i,j-1}}\right\rangle, for $j\in\{2,3,\dots,k\}$, we get

[TABLE]

Then, since \mathopen{}\mathclose{{}\left\|Ah}\right\|\leq 1, we get

[TABLE]

Next, using the fact \mathopen{}\mathclose{{}\left|x_{i+1,0}}\right\rangle=\mathopen{}\mathclose{{}\left|x_{i,0}}\right\rangle+\sum_{j=1}^{k}\mathopen{}\mathclose{{}\left|x_{i,j}}\right\rangle and the triangle inequality, we get

[TABLE]

which implies

[TABLE]

Then it follows from Eq. (99) that

[TABLE]

Now using Eqs. (96), (101), and (102), we obtain

[TABLE]

where the last step follows from

[TABLE]

as in Eq. (50). Thus, combining Eqs. (94), (97), and (103) we get

[TABLE]

as claimed. ∎

Theorem 91 implies that by choosing $p=m$, we can make the probability at least $1/78g^{2}$ for obtaining the state \mathopen{}\mathclose{{}\left|x_{m,j}}\right\rangle/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{m,j}}\right\rangle}\right\| when measuring the first register of \mathopen{}\mathclose{{}\left|x}\right\rangle/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x}\right\rangle}\right\| in the standard basis. This probability can be increased to \Omega\mathopen{}\mathclose{{}\left(1}\right) by amplitude amplification, which uses O\mathopen{}\mathclose{{}\left(g}\right) repetitions of the above procedure.

6 State Preparation

To apply the QLSA to a linear system of the form $M|x\rangle=|y\rangle$, we must also be able to prepare the state $|y\rangle$. To quantify the complexity of this subroutine, Ref. [3] assumes that $\vec{x}_{\mathrm{in}}$ and $\vec{b}$ are sparse vectors whose entries are given by oracles. Instead of assuming sparsity, here we simply assume that we have controlled oracles that produce states proportional to $\vec{x}_{\mathrm{in}}$ and $\vec{b}$, respectively.

The following lemma shows how to use these oracles to produce the state appearing on the right-hand side of our linear system. We write \mathopen{}\mathclose{{}\left|\bar{\varphi}}\right\rangle to denote the normalized version of \mathopen{}\mathclose{{}\left|\varphi}\right\rangle, i.e., \mathopen{}\mathclose{{}\left|\bar{\varphi}}\right\rangle\vcentcolon=\mathopen{}\mathclose{{}\left|\varphi}\right\rangle/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|\varphi}\right\rangle}\right\|, for any \mathopen{}\mathclose{{}\left|\varphi}\right\rangle.

Lemma 8**.**

Let $\mathcal{O}_{x}$ be a unitary that maps \mathopen{}\mathclose{{}\left|1}\right\rangle\mathopen{}\mathclose{{}\left|\varphi}\right\rangle to \mathopen{}\mathclose{{}\left|1}\right\rangle\mathopen{}\mathclose{{}\left|\varphi}\right\rangle for any \mathopen{}\mathclose{{}\left|\varphi}\right\rangle and maps \mathopen{}\mathclose{{}\left|0}\right\rangle\mathopen{}\mathclose{{}\left|0}\right\rangle to \mathopen{}\mathclose{{}\left|0}\right\rangle\mathopen{}\mathclose{{}\left|\bar{x}_{\mathrm{in}}}\right\rangle, where \bar{x}_{\mathrm{in}}=\vec{x}_{\mathrm{in}}/\mathopen{}\mathclose{{}\left\|\vec{x}_{\mathrm{in}}}\right\|. Let $\mathcal{O}_{b}$ be a unitary that maps \mathopen{}\mathclose{{}\left|0}\right\rangle\mathopen{}\mathclose{{}\left|\varphi}\right\rangle to \mathopen{}\mathclose{{}\left|0}\right\rangle\mathopen{}\mathclose{{}\left|\varphi}\right\rangle for any \mathopen{}\mathclose{{}\left|\varphi}\right\rangle and maps \mathopen{}\mathclose{{}\left|1}\right\rangle\mathopen{}\mathclose{{}\left|0}\right\rangle to \mathopen{}\mathclose{{}\left|1}\right\rangle\mathopen{}\mathclose{{}\left|\bar{b}}\right\rangle, where $\bar{b}=\vec{b}/\|\vec{b}\|$ . Suppose we know \mathopen{}\mathclose{{}\left\|\vec{x}_{\mathrm{in}}}\right\| and $\|\vec{b}\|$. Then the state proportional to

[TABLE]

can be produced with a constant number of calls to $\mathcal{O}_{x}$ and $\mathcal{O}_{b}$, and \operatorname{poly}\mathopen{}\mathclose{{}\left(\operatorname{log}\mathopen{}\mathclose{{}\left(mk}\right)}\right) elementary gates.

Proof.

Consider the initial state $|0\rangle|0\rangle$, where the first register is the $(d+1)$-dimensional register corresponding to the block-level indexing of $C_{m,k,p}$, and the second register is the $N$-dimensional register that stores $\vec{b}$ and $\vec{x}_{\mathrm{in}}$. We perform the unitary

[TABLE]

on the first register to get the state

[TABLE]

Next, we apply the unitaries $\mathcal{O}_{x}$ and $\mathcal{O}_{b}$ (in arbitrary order), and obtain the state

[TABLE]

Finally, we apply a unitary that maps \mathopen{}\mathclose{{}\left|0}\right\rangle to \mathopen{}\mathclose{{}\left|0}\right\rangle and maps \mathopen{}\mathclose{{}\left|1}\right\rangle to $\frac{1}{\sqrt{m}}\sum_{j=0}^{m-1}|j(k+1)+1\rangle$ on the first register to get the state we need. This can be done by standard techniques using \operatorname{poly}\mathopen{}\mathclose{{}\left(\operatorname{log}\mathopen{}\mathclose{{}\left(mk}\right)}\right) elementary gates. ∎

7 Main Result

In this section, we formally state and prove our main result on the quantum algorithm for linear ordinary differential equations.

Theorem 9**.**

Suppose $A=VDV^{-1}$ is an $N\times N$ diagonalizable matrix, where D=\operatorname{diag}\mathopen{}\mathclose{{}\left(\lambda_{0},\lambda_{1},\dots,\lambda_{N-1}}\right) satisfies $\mathrm{Re}(\lambda_{j})\leq 0$ for any $j\in\{0,1,\dots,N-1\}$. In addition, suppose $A$ has at most $s$ nonzero entries in any row and column, and we have an oracle $\mathcal{O}_{A}$ that computes these entries. Suppose $\vec{x}_{\mathrm{in}}$ and $\vec{b}$ are $N$-dimensional vectors with known norms and that we have two controlled oracles, $\mathcal{O}_{x}$ and $\mathcal{O}_{b}$, that prepare the states proportional to $\vec{x}_{\mathrm{in}}$ and $\vec{b}$, respectively. Let $\vec{x}$ evolve according to the differential equation

[TABLE]

with the initial condition $\vec{x}(0)=\vec{x}_{\mathrm{in}}$. Let $T>0$ and

[TABLE]

Then there exists a quantum algorithm that produces a state $\epsilon$-close to \vec{x}(T)/\mathopen{}\mathclose{{}\left\|\vec{x}(T)}\right\| in $\ell^{2}$ norm, succeeding with probability \Omega\mathopen{}\mathclose{{}\left(1}\right), with a flag indicating success, using

[TABLE]

queries to $\mathcal{O}_{A}$, $\mathcal{O}_{x}$, and $\mathcal{O}_{b}$, where \kappa_{V}=\mathopen{}\mathclose{{}\left\|V}\right\|\cdot\mathopen{}\mathclose{{}\left\|V^{-1}}\right\| is the condition number of $V$ and \beta=(\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{\mathrm{in}}}\right\rangle}\right\|+T\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|b}\right\rangle}\right\|)/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x(T)}\right\rangle}\right\|. The gate complexity of this algorithm is larger than its query complexity by a factor of \operatorname{poly}\mathopen{}\mathclose{{}\left(\operatorname{log}\mathopen{}\mathclose{{}\left(\kappa_{V}sg\beta T\mathopen{}\mathclose{{}\left\|A}\right\|N/\epsilon}\right)}\right).

Proof.

Recall that we use \mathopen{}\mathclose{{}\left|\bar{\varphi}}\right\rangle to denote the normalized version of \mathopen{}\mathclose{{}\left|\varphi}\right\rangle.

Statement of the Algorithm

Let h=T/\lceil T\mathopen{}\mathclose{{}\left\|A}\right\|\rceil, m=p=T/h=\lceil T\mathopen{}\mathclose{{}\left\|A}\right\|\rceil, $\delta=\epsilon/(25\sqrt{m}g)$, $\epsilon\leq 1/2$, and

[TABLE]

with

[TABLE]

This choice of $k$ ensures that $(k+1)!\geq\Omega$. Moreover, since $\Omega\geq 70$, $k\geq 5$. Hence this is an appropriate choice of $k$ for the conditions of Theorem 91. We build the linear system

[TABLE]

We use the QLSA from Ref. [11] to solve this linear system and obtain a state \mathopen{}\mathclose{{}\left|\bar{x}^{\prime}}\right\rangle such that \mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|\bar{x}}\right\rangle-\mathopen{}\mathclose{{}\left|\bar{x}^{\prime}}\right\rangle}\right\|\leq\delta. Then we measure the first register of \mathopen{}\mathclose{{}\left|\bar{x}^{\prime}}\right\rangle in the standard basis, and conditioned on the outcome being in

[TABLE]

we output the state of the second register. We will show that the probability of this event happening is \Omega\mathopen{}\mathclose{{}\left(1/g^{2}}\right). Using amplitude amplification [9], we can raise this probability to \Omega\mathopen{}\mathclose{{}\left(1}\right) with O\mathopen{}\mathclose{{}\left(g}\right) repetitions of the above procedure.

Proof of Correctness

Let $d=m(k+1)+p$. Let \mathopen{}\mathclose{{}\left|x_{i,j}}\right\rangle be defined by Eq. (10), and let \mathopen{}\mathclose{{}\left|x_{l}}\right\rangle=\mathopen{}\mathclose{{}\left|x_{i,j}}\right\rangle for $l=i(k+1)+j$. Then we have

[TABLE]

Note that \mathopen{}\mathclose{{}\left|x_{l}}\right\rangle=\mathopen{}\mathclose{{}\left|x_{m,0}}\right\rangle for any $l\in S$. Then by Theorem 6 and our choice of parameters, we have for any $l\in S$,

[TABLE]

Now let $\alpha_{l}$ be such that

[TABLE]

Then for any $l\in S$, we have (by Theorem 91 and our choice of parameters)

[TABLE]

By Eq. (118) and Lemma 142 in Appendix B, we have

[TABLE]

Now suppose the QLSA outputs the state

[TABLE]

which satisfies

[TABLE]

Then for any $l\in S$, by Lemma 14 in Appendix B, we have

[TABLE]

Then it follows from Eqs. (121) and (124) that for any $l\in S$,

[TABLE]

Furthermore, by Lemma 15 in Appendix B, we have

[TABLE]

Therefore, if we measure the first register of \mathopen{}\mathclose{{}\left|\bar{x}^{\prime}}\right\rangle in the standard basis, the probability of getting some outcome $l\in S$ is

[TABLE]

and when this happens the state of the second register becomes \mathopen{}\mathclose{{}\left|\bar{x}_{l}^{\prime}}\right\rangle for some $l\in S$, which is $\epsilon$-close to the desired \mathopen{}\mathclose{{}\left|\bar{x}(T)}\right\rangle in $\ell^{2}$ norm. The success probability can be raised to \Omega\mathopen{}\mathclose{{}\left(1}\right) by using $O(g)$ rounds of amplitude amplification.

Analysis of the Complexity

The matrix $C_{m,k,p}(A)$ is a $(d+1)N\times(d+1)N$ matrix with O\mathopen{}\mathclose{{}\left(ks}\right) nonzero entries in any row or column. By Theorem 5 and our choice of parameters, the condition number of $C_{m,k,p}(A)$ is O\mathopen{}\mathclose{{}\left(\kappa_{V}km}\right). Consequently, by Theorem 5 of Ref. [11], the QLSA produces the state \mathopen{}\mathclose{{}\left|x^{\prime}}\right\rangle with

[TABLE]

queries to the oracles $\mathcal{O}_{A}$, $\mathcal{O}_{x}$, and $\mathcal{O}_{b}$, where \beta=\mathopen{}\mathclose{{}\left(\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x_{\mathrm{in}}}\right\rangle+T\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|b}\right\rangle}\right\|}\right\|}\right)/\mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|x(T)}\right\rangle}\right\|, and its gate complexity is larger by a factor of \operatorname{poly}\mathopen{}\mathclose{{}\left(\operatorname{log}\mathopen{}\mathclose{{}\left(\kappa_{V}kmsN/\delta}\right)}\right)=\operatorname{poly}\mathopen{}\mathclose{{}\left(\operatorname{log}\mathopen{}\mathclose{{}\left(\kappa_{V}sg\beta T\mathopen{}\mathclose{{}\left\|A}\right\|N/\epsilon}\right)}\right). Since amplitude amplification requires only O\mathopen{}\mathclose{{}\left(g}\right) repetitions of this procedure, the query complexity of our algorithm is O\mathopen{}\mathclose{{}\left(\kappa_{V}sgT\mathopen{}\mathclose{{}\left\|A}\right\|\cdot\operatorname{poly}\mathopen{}\mathclose{{}\left(\operatorname{log}\mathopen{}\mathclose{{}\left(\kappa_{V}sg\beta T\mathopen{}\mathclose{{}\left\|A}\right\|/\epsilon}\right)}\right)}\right) and its gate complexity is larger by a factor of \operatorname{poly}\mathopen{}\mathclose{{}\left(\operatorname{log}\mathopen{}\mathclose{{}\left(\kappa_{V}sg\beta T\mathopen{}\mathclose{{}\left\|A}\right\|N/\epsilon}\right)}\right), as claimed. ∎

8 Discussion

In this paper, we have presented an algorithm for solving (possibly inhomogeneous) linear ordinary differential equations with constant coefficients, with exponentially improved performance over the algorithm of Ref. [3].

The complexity of our algorithm depends on the parameter $g$ defined in Eq. (111), which characterizes the extent to which the final solution vector decays relative to the solution vector at earlier times. The success probability obtained from a single solution of the linear system (as analyzed in Section 5) decays with $g$, so we boost the success probability using amplitude amplification. Dramatically improved dependence on $g$ is unlikely since a simulation of evolution in which the state decays can be used to implement postselection, following a procedure like that of Proposition 5 of Ref. [1]. In particular, the ability to postselect on an exponentially small amplitude would imply $\textsf{BQP}=\textsf{PP}$, which is considered implausible. Note that the algorithm of Ref. [3] has essentially the same limitation (although there it was stated in terms of an assumed upper bound on $g$).

Aside from its exponentially improved dependence on $\epsilon$, our approach has other advantages over Ref. [3]. First, since linear multistep methods may be unstable, Ref. [3] requires careful selection of a stable method. In contrast, the propagator \operatorname{exp}\mathopen{}\mathclose{{}\left(At}\right) exactly evolves the state forward in time, so we need not concern ourselves with numerical stability. Furthermore, our algorithm has nearly linear scaling with respect to the evolution time $T$, which is a quadratic improvement over Ref. [3]. Since Hamiltonian simulation is a special case of our algorithm with $g=1$, the no-fast-forwarding theorem [4] implies that the dependence of our algorithm on $T$ is optimal up to logarithmic factors. We also obtain polynomial improvements over Ref. [3] for the dependence of the complexity on the sparsity $s$, norm \mathopen{}\mathclose{{}\left\|A}\right\|, and condition number $\kappa_{V}$ of the transformation that diagonalizes $A$.

Although Theorem 9 assumes that $A$ is diagonalizable, our algorithm can also produce approximate solutions for non-diagonalizable $A$. This is because diagonalizable matrices are dense within the set of all complex matrices: for any non-diagonalizable matrix $A$ and any $\delta>0$, there is a diagonalizable matrix $B$ such that \mathopen{}\mathclose{{}\left\|A-B}\right\|<\delta. Using such a $B$ in place of $A$, we can simulate Eq. (1) approximately. This approach can yield a matrix $B$ whose diagonalizing transformation has a condition number polynomial in $1/\delta$, so the complexity would no longer be \operatorname{poly}\mathopen{}\mathclose{{}\left(\operatorname{log}\mathopen{}\mathclose{{}\left(1/\epsilon}\right)}\right).

Our work raises some natural open problems. For simplicity, we have assumed that the matrix $A$ and the inhomogeneity $\vec{b}$ are independent of $t$. More generally, can one solve differential equations with time-dependent coefficients with complexity \operatorname{poly}\mathopen{}\mathclose{{}\left(\operatorname{log}\mathopen{}\mathclose{{}\left(1/\epsilon}\right)}\right)? This is possible in the case of quantum simulation, i.e., when $A(t)$ is anti-Hermitian and $\vec{b}=0$ [5]. However, some aspects of our analysis appear to fail in the time-dependent case.

Finally, while our algorithm has optimal dependence on $\epsilon$ and nearly optimal dependence on $T$, the joint dependence on these parameters might be improved. Recent work gave an algorithm for Hamiltonian simulation with complexity \smash{O(t+\frac{\operatorname{log}\mathopen{}\mathclose{{}\left(1/\epsilon}\right)}{\operatorname{log}\operatorname{log}\mathopen{}\mathclose{{}\left(1/\epsilon}\right)})}, providing an optimal tradeoff. Can similar complexity be attained for more general differential equations?

Acknowledgements

This research was supported by the Canadian Institute for Advanced Research, the National Science Foundation (grant number 1526380), and IARPA (contract number D15PC00242). In addition, DWB is funded by an Australian Research Council Discovery Project (DP160102426).

Appendix A Lemmas About Taylor Series

In this appendix, we establish some basic lemmas about the approximation of functions by truncated Taylor series.

Lemma 10**.**

Let $k\in\mathbb{Z}^{+}$, $z\in\mathbb{C}$, \mathopen{}\mathclose{{}\left|z}\right|\leq 1, and $\mathrm{Re}(z)\leq 0$. Define $T_{k}(z)\vcentcolon=\sum_{j=0}^{k}\frac{z^{j}}{j!}$. Then \mathopen{}\mathclose{{}\left|T_{k}(z)-\operatorname{exp}\mathopen{}\mathclose{{}\left(z}\right)}\right|\leq 1/(k+1)! and \mathopen{}\mathclose{{}\left|T_{k}(z)}\right|\leq 1+1/(k+1)!.

Proof.

Given $\varphi$ such that $z=e^{i\varphi}|z|$, define a function f(x):=\operatorname{exp}\mathopen{}\mathclose{{}\left(e^{i\varphi}x}\right) for $x\in\mathbb{R}$. Then the integral form of the remainder for Taylor’s theorem gives

[TABLE]

so

[TABLE]

Because $\mathrm{Re}(z)\leq 0$, |\operatorname{exp}\mathopen{}\mathclose{{}\left(e^{i\varphi}t}\right)|\leq 1. Taking $x=|z|$, we therefore have

[TABLE]

Using this bound,

[TABLE]

and the lemma follows. ∎

Lemma 11**.**

Let $b\in\mathbb{Z}^{*}$, $k\in\mathbb{Z}^{+}$, $k\geq 5$, and $b\leq k$. Let $z\in\mathbb{C}$, \mathopen{}\mathclose{{}\left|z}\right|\leq 1, and $\mathrm{Re}(z)\leq 0$. Define $T_{b,k}(z)\vcentcolon=\sum_{j=b}^{k}\frac{b!z^{j-b}}{j!}$. Then \mathopen{}\mathclose{{}\left|T_{b,k}(z)}\right|\leq\sqrt{1.04}.

Proof.

Again we can take $z=e^{i\varphi}|z|$ and f(x):=\operatorname{exp}\mathopen{}\mathclose{{}\left(e^{i\varphi}x}\right) for $x\in\mathbb{R}$. If $b=0$, then $T_{b,k}(z)=T_{k}(z)$, so from Lemma 10 we have \mathopen{}\mathclose{{}\left|T_{b,k}(z)}\right|\leq 1+1/(k+1)!<\sqrt{1.04}. If $b=k$, then \mathopen{}\mathclose{{}\left|T_{b,k}(z)}\right|=1. If $b=k-1$, then

[TABLE]

Then we obtain

[TABLE]

since $\mathrm{Re}(z)\leq 0$ and \mathopen{}\mathclose{{}\left|z}\right|\leq 1. Then, for $k\geq 5$, we obtain \mathopen{}\mathclose{{}\left|T_{b,k}(z)}\right|\leq\sqrt{1.04}.

Otherwise, $1\leq b<k-1$, and using the integral form of the remainder gives

[TABLE]

Then limiting the sum to $k$ gives

[TABLE]

Therefore, we obtain a formula for $T_{b,k}(z)$ as

[TABLE]

Again |\operatorname{exp}\mathopen{}\mathclose{{}\left(e^{i\varphi}t}\right)|\leq 1, so taking the absolute value gives us

[TABLE]

where in the last line we have used $k\geq 5$. Hence for $0\leq b\leq k$, we obtain \mathopen{}\mathclose{{}\left|T_{b,k}(z)}\right|\leq\sqrt{1.04} as required. ∎

Lemma 12**.**

Let $k\in\mathbb{Z}^{+}$. Let $z\in\mathbb{C}$, $|z|\leq 1$. Define $S_{k}(z)\vcentcolon=\sum_{j=1}^{k}\frac{z^{j-1}}{j!}$. Then \mathopen{}\mathclose{{}\left|S_{k}(z)-(\operatorname{exp}\mathopen{}\mathclose{{}\left(z}\right)-1)z^{-1}}\right|\leq 1/(k+1)!.

Proof.

Again we can take $z=e^{i\varphi}|z|$ and f(x):=\operatorname{exp}\mathopen{}\mathclose{{}\left(e^{i\varphi}x}\right) for $x\in\mathbb{R}$. Then using Eq. (130),

[TABLE]

Dividing by $e^{i\varphi}x$ gives

[TABLE]

Taking the absolute value then gives

[TABLE]

as claimed. ∎

Appendix B Lemmas About Quantum States

In this appendix, we prove some technical lemmas about approximation of quantum states.

Lemma 13**.**

Let \mathopen{}\mathclose{{}\left|\psi}\right\rangle and \mathopen{}\mathclose{{}\left|\varphi}\right\rangle be two vectors such that \mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|\psi}\right\rangle}\right\|\geq\alpha>0 and \mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left|\psi}\right\rangle-\mathopen{}\mathclose{{}\left|\varphi}\right\rangle}\right\|\leq\beta. Then

[TABLE]

Proof.

Using the triangle inequality, we get

[TABLE]

as claimed. ∎

Lemma 14**.**

Let $|\psi\rangle=\alpha|0\rangle|\psi_{0}\rangle+\sqrt{1-\alpha^{2}}|1\rangle|\psi_{1}\rangle$ and $|\varphi\rangle=\beta|0\rangle|\varphi_{0}\rangle+\sqrt{1-\beta^{2}}|1\rangle|\varphi_{1}\rangle$, where \mathopen{}\mathclose{{}\left|\psi_{0}}\right\rangle, \mathopen{}\mathclose{{}\left|\psi_{1}}\right\rangle, \mathopen{}\mathclose{{}\left|\varphi_{0}}\right\rangle, \mathopen{}\mathclose{{}\left|\varphi_{1}}\right\rangle are unit vectors, and $\alpha,\beta\in[0,1]$. Suppose \mathopen{}\mathclose{{}\left\||\psi\rangle-|\varphi\rangle}\right\|\leq\delta<\alpha. Then \mathopen{}\mathclose{{}\left\||\varphi_{0}\rangle-|\psi_{0}\rangle}\right\|\leq\frac{2\delta}{\alpha-\delta}.

Proof.

First note that \mathopen{}\mathclose{{}\left\||\psi\rangle-|\varphi\rangle}\right\|\leq\delta implies

[TABLE]

In addition, by the triangle inequality, we have

[TABLE]

Then Eqs. (144) and (145) imply

[TABLE]

Then by the triangle inequality, we get

[TABLE]

as claimed. ∎

Lemma 15**.**

Let $|\psi\rangle=\alpha|0\rangle|\psi_{0}\rangle+\sqrt{1-\alpha^{2}}|1\rangle|\psi_{1}\rangle$ and $|\varphi\rangle=\beta|0\rangle|\varphi_{0}\rangle+\sqrt{1-\beta^{2}}|1\rangle|\varphi_{1}\rangle$, where \mathopen{}\mathclose{{}\left|\psi_{0}}\right\rangle, \mathopen{}\mathclose{{}\left|\psi_{1}}\right\rangle, \mathopen{}\mathclose{{}\left|\varphi_{0}}\right\rangle, \mathopen{}\mathclose{{}\left|\varphi_{1}}\right\rangle are unit vectors, and $\alpha,\beta\in[0,1]$. Suppose \mathopen{}\mathclose{{}\left\||\psi\rangle-|\varphi\rangle}\right\|\leq\delta<\alpha. Then $\beta\geq\alpha-\delta$.

Proof.

Note that \mathopen{}\mathclose{{}\left\||\psi\rangle-|\varphi\rangle}\right\|\leq\delta implies \mathopen{}\mathclose{{}\left\|\mathopen{}\mathclose{{}\left\langle 0|\psi}\right\rangle-\mathopen{}\mathclose{{}\left\langle 0|\varphi}\right\rangle)}\right\|\leq\delta. Then by the triangle inequality, we get

[TABLE]

as claimed. ∎