Quantum spectral methods for differential equations

TL;DR

This paper introduces a quantum spectral method algorithm for solving linear ordinary differential equations with time-dependent coefficients, achieving exponential speedup in problem size and accuracy compared to classical methods.

Contribution

It develops the first quantum algorithm for time-dependent linear differential equations using spectral methods, expanding quantum numerical analysis capabilities.

Findings

Quantum algorithm solves time-dependent ODEs with poly(log d, log 1/ε) complexity.

Spectral methods provide a global approximation approach in quantum algorithms.

Algorithm extends quantum linear algebra techniques to differential equations with time-dependent coefficients.

Abstract

Recently developed quantum algorithms address computational challenges in numerical analysis by performing linear algebra in Hilbert space. Such algorithms can produce a quantum state proportional to the solution of a $d$-dimensional system of linear equations or linear differential equations with complexity $\mathrm{poly}(\log d)$. While several of these algorithms approximate the solution to within $\epsilon$ with complexity $\mathrm{poly}(\log(1/\epsilon))$, no such algorithm was previously known for differential equations with time-dependent coefficients. Here we develop a quantum algorithm for linear ordinary differential equations based on so-called spectral methods, an alternative to finite difference methods that approximates the solution globally. Using this approach, we give a quantum algorithm for time-dependent initial and boundary value problems with complexity $\mathrm{poly}(\log d, \log(1/\epsilon))$.

Equations

$$\frac{\mathrm{d}x(t)}{\mathrm{d}t}=A(t)x(t)+f(t)$$

$$\frac{\mathrm{d}x(t)}{\mathrm{d}t}=\frac{x(t+h)-x(t-h)}{2h}+O(h^2).$$

$$x(0)=\gamma$$

$$x_i(t)=\sum _{k=0}^n{c}_{i,k}{T}_k(t),i\in [d{]}_0:=\{0,1,\dots ,d-1\}$$

$$\frac{\mathrm{d}x(t_l)}{\mathrm{d}t}=A(t_l)x(t_l)+f(t_l),\forall l\in [n+1],t\in [-1,1],$$

$$x_i(t_0)={\gamma }_i,i\in [d{]}_0.$$

$$\Vert \hat{x}(t)-x(t)\Vert \to 0\text{as}n\to \infty .$$

$$\frac{\mathrm{d}{x}_i(t)}{\mathrm{d}t}=\sum _{k=0}^n{c}_{i,k}^{\prime }{T}_k(t).$$

$$2T_k(t)=\frac{T_{k+1}^{\prime }(t)}{k+1}-\frac{T_{k-1}^{\prime }(t)}{k-1},$$

$$c_{i,k}^{\prime }=\sum _{j=0}^n[D_n{]}_{kj}{c}_{i,j},i\in [d{]}_0,k\in [n+1{]}_0,$$

$$[D_n{]}_{kj}=\frac{2j}{{\sigma }_k},\text{k+j odd},j>k,$$

$${\sigma }_k:=\{\begin{array}{ll}2& k=0\\ 1& k\in [n]:=\{1,2,\dots ,n\}.\end{array}$$

$$\sum _{k=0}^n{T}_k(t_l)c_{i,k}^{\prime }=\sum _{j=0}^{d-1}{A}_{ij}(t_l)\sum _{k=0}^n{T}_k(t_l)c_{j,k}+f(t_l{)}_i,i\in [d{]}_0,l\in [n+1{]}_0.$$

$$\underset{t\in [-1,1]}{\max }\Vert \hat{x}(t)-x(t)\Vert \le C\underset{t\in [-1,1]}{\max }\frac{\Vert {\hat{x}}^{(n+1)}(t)\Vert }{n^{r-2}}.$$

$$\underset{t\in [-1,1]}{\max }\Vert \hat{x}(t)-x(t)\Vert \le \sqrt{\frac{2}{\pi }}\underset{t\in [-1,1]}{\max }\Vert {\hat{x}}^{(n+1)}(t)\Vert {\left(\frac{e}{2n}\right)}^n.$$

$$\frac{\mathrm{d}x(t)}{\mathrm{d}t}=A(t)x(t)+f(t)$$

$$K_h:t\mapsto 1-\frac{2(t-{\Gamma }_h)}{{\Gamma }_{h+1}-{\Gamma }_h},$$

$$A_h(t)$$

$$x_h(t)$$

$$f_h(t)$$

$$\frac{\mathrm{d}{x}_h}{\mathrm{d}t}=A_h(t)x_h(t)+f_h(t)$$

$$x_h(1)=\{\begin{array}{ll}\gamma & h=0\\ x_{h-1}(-1)& h\in [m].\end{array}$$

$${\tau }_h\le \frac{2}{{\max }_{t\in [{\Gamma }_h,{\Gamma }_{h+1}]}\Vert A(t)\Vert }$$

$$\tau :=\underset{h\in \{0,1,\dots ,m-1\}}{\max }{\tau }_h\le \frac{2}{{\max }_{t\in [0,T]}\Vert A(t)\Vert }.$$

$$\frac{\mathrm{d}x(t_l)}{\mathrm{d}t}=A_h(t_l)x(t_l)+f_h(t_l),h\in [m{]}_0,l\in [n+1]$$

$$x(t_0)=\gamma .$$

$$L|X⟩=|B⟩$$

$$|X⟩=\sum _{h=0}^{m-1}\sum _{i=0}^{d-1}\sum _{l=0}^n{c}_{i,l}({\Gamma }_{h+1})|hil⟩+\sum _{h=m}^{m+p}\sum _{i=0}^{d-1}\sum _{l=0}^n{x}_i|hil⟩$$

$$|B⟩=\sum _{h=0}^{m-1}|h⟩|B(f_h)⟩$$

$$|B(f_h)⟩=\sum _{i=0}^{d-1}{\gamma }_i|i0⟩+\sum _{i=0}^{d-1}\sum _{l=1}^n{f}_h(\cos \frac{l\pi }{n}{)}_i|il⟩,h\in [m-1].$$

$$L=\sum _{h=0}^{m-1}|h⟩⟨h|\otimes (L_1+L_2(A_h))+\sum _{h=1}^m|h⟩⟨h-1|\otimes L_3+\sum _{h=m}^{m+p}|h⟩⟨h|\otimes L_4+\sum _{h=m+1}^{m+p}|h⟩⟨h-1|\otimes L_5.$$

Full text

1 Introduction

Differential equations have extensive applications throughout mathematics, science, and engineering. Numerical methods for differential equations have been widely studied (see for example [26]), giving fast algorithms for solving them using classical computers.

Recent work has developed quantum algorithms with the potential to extract information about solutions of systems of differential equations even faster than is possible classically. This body of work grew from the quantum linear systems algorithm (QLSA) [21], which produces a quantum state proportional to the solution of a sparse system of $d$ linear equations in time $\operatorname{poly}(\log d)$. Subsequent work improved the performance of that algorithm [1, 15] and applied it to develop similar quantum algorithms for differential equations.

To achieve this improvement, quantum algorithms must operate under different assumptions than those made for algorithms in classical numerical analysis. To represent the output using $\operatorname{poly}(\log d)$ qubits, the output is produced as a quantum state, not as an explicit description of a vector. Furthermore, to facilitate applying Hamiltonian simulation, these quantum algorithms use implicit access to the system of equations (say, through a matrix specified by a sparse Hamiltonian oracle and the ability to prepare quantum states encoding certain vectors). While these assumptions restrict the types of equations that can be addressed and the types of information that can be extracted from the final state, they nevertheless appear capable of producing useful information that cannot be efficiently computed classically.

In this paper, we focus on systems of first-order linear ordinary differential equations (ODEs). Such equations can be written in the form

[TABLE]

where $t\in[0,T]$ for some $T>0$, the solution $x(t)\in{\mathbb{C}}^{d}$ is a $d$-dimensional vector, and the system is determined by a time-dependent matrix $A(t)\in{\mathbb{C}}^{d\times d}$ and a time-dependent inhomogeneity $f(t)\in{\mathbb{C}}^{d}$. Provided $A(t)$ and $f(t)$ are continuous functions of $t$, the initial value problem (i.e., the problem of determining $x(t)$ for a given initial condition $x(0)$) has a unique solution [2].

The Hamiltonian simulation problem is simply the special case of the quantum ODE problem where $A$ is antihermitian and $f$ is zero. A substantial body of work has developed fast quantum algorithms for that special case [6, 31, 7, 8, 9, 13, 11, 28, 29, 30]. Hamiltonian simulation underlies the QLSA [21, 15] which in turn gives algorithms for more general differential equations.

Berry presented the first efficient quantum algorithm for general linear ODEs [5]. His algorithm represents the system of differential equations as a system of linear equations using a linear multistep method and solves that system using the QLSA. This approach achieves complexity logarithmic in the dimension $d$ and, by using a high-order integrator, close to quadratic in the evolution time $T$. While this method could in principle be applied to handle time-dependent equations, the analysis of [5] only explicitly considers the time-independent case for simplicity.

Since it uses a finite difference approximation, the complexity of Berry’s algorithm as a function of the solution error $\epsilon$ is $\operatorname{poly}(1/\epsilon)$ [5]. Reference [10] improved this to $\operatorname{poly}(\log(1/\epsilon))$ by using a high-precision QLSA based on linear combinations of unitaries [15] to solve a linear system that encodes a truncated Taylor series. However, this approach assumes that $A(t)$ and $f(t)$ are time-independent so that the solution of the ODE can be written as an explicit series, and it is unclear how to generalize the algorithm to time-dependent ODEs.

While we focus here on extending the above line of work, several other approaches have been proposed for addressing differential equations with quantum computers. Reference [27] used a quantum version of the Euler method to handle nonlinear ODEs with polynomial nonlinearities. This algorithm has complexity logarithmic in the dimension but exponential in the evolution time (as is inevitable for general nonlinear ODEs). Other work has developed quantum algorithms for partial differential equations (PDEs). Reference [16] described a quantum algorithm that applies the QLSA to implement a finite element method for Maxwell’s equations. Reference [17] applied Hamiltonian simulation to a finite difference approximation of the wave equation. Most recently, reference [3] presented a continuous-variable quantum algorithm for initial value problems with non-homogeneous linear PDEs.

Most of the aforementioned algorithms use a local approximation: they discretize the differential equations into small time intervals to obtain a system of linear equations or linear differential equations that can be solved by the QLSA or Hamiltonian simulation. For example, the central difference scheme approximates the time derivative at the point $x(t)$ as

[TABLE]

High-order finite difference or finite element methods can reduce the error to $O(h^{k})$, where $k-1$ is the order of the approximation. However, when solving an equation over the interval $[0,T]$, the number of iterations is $T/h=\Theta(\epsilon^{-1/k})$ for fixed $k$, giving a total complexity that is $\operatorname{poly}(1/\epsilon)$ even using high-precision methods for the QLSA or Hamiltonian simulation.

For ODEs with special structure, some prior results already show how to avoid a local approximation and thereby achieve complexity $\operatorname{poly}(\log(1/\epsilon))$. When $A(t)$ is anti-Hermitian and $f(t)=0$, we can directly apply Hamiltonian simulation [9]; if $A$ and $f$ are time-independent, then [10] uses a Taylor series to achieve complexity $\operatorname{poly}(\log({1}/{\epsilon}))$. However, the case of general time-dependent linear ODEs had remained elusive.

In this paper, we use a nonlocal representation of the solution of a system of differential equations to give a new quantum algorithm with complexity $\operatorname{poly}(\log(1/\epsilon))$ even for time-dependent equations. Time-dependent linear differential equations describe a wide variety of systems in science and engineering. Examples include the wave equation and the Stokes equation (i.e., creeping flow) in fluid dynamics [25], the heat equation and the Boltzmann equation in thermodynamics [34, 20], the Poisson equation and Maxwell’s equations in electromagnetism [23, 35], and of course Schrödinger’s equation in quantum mechanics. Moreover, some nonlinear differential equations can be studied by linearizing them to produce time-dependent linear equations (e.g., the linearized advection equation in fluid dynamics [12]).

We focus our discussion on first-order linear ODEs. Higher-order ODEs can be transformed into first-order ODEs by standard methods. Also, by discretizing space, PDEs with both time and space dependence can be regarded as sparse linear systems of time-dependent ODEs. Thus we focus on an equation of the form (1.1) with initial condition

[TABLE]

for some specified $\gamma\in{\mathbb{C}}^{d}$. We assume that $A(t)$ is $s$-sparse (i.e., has at most $s$ nonzero entries in any row or column) for any $t\in[0,T]$. Furthermore, we assume that $A(t)$, $f(t)$, and $\gamma$ are provided by black-box subroutines (which serve as abstractions of efficient computations). In particular, following essentially the same model as in [10] (see also Section 1.1 of [15]), suppose we have an oracle $O_{A}(t)$ that, for any $t\in[0,T]$ and any given row or column specified as input, computes the locations and values of the nonzero entries of $A(t)$ in that row or column. We also assume oracles $O_{x}$ and $O_{f}(t)$ that, for any $t\in[0,T]$, prepare normalized states $|\gamma\rangle$ and $|f(t)\rangle$ proportional to $\gamma$ and $f(t)$, and that also compute $\|\gamma\|$ and $\|f(t)\|$, respectively. Given such a description of the instance, the goal is to produce a quantum state $\epsilon$-close to $|x(T)\rangle$ (a normalized quantum state proportional to $x(T)$).

As mentioned above, our main contribution is to implement a method that uses a global approximation of the solution. We do this by developing a quantum version of so-called spectral methods, a technique from classical numerical analysis that (approximately) represents the components of the solution $x(t)_{i}\approx\sum_{j}c_{ij}\phi_{j}(t)$ as linear combinations of basis functions $\phi_{j}(t)$ expressing the time dependence. Specifically, we implement a Chebyshev pseudospectral method [4, 22] using the QLSA. This approach approximates the solution by a truncated Chebyshev series with undetermined coefficients and solves for those coefficients using a linear system that interpolates the differential equations. According to the convergence theory of spectral methods, the solution error decreases exponentially provided the solution is sufficiently smooth [19, 32]. We use the LCU-based QLSA to solve this linear system with high precision [15]. To analyze the algorithm, we upper bound the solution error and condition number of the linear system and lower bound the success probability of the final measurement. Overall, we show that the total complexity of this approach is $\operatorname{poly}(\log(1/\epsilon))$ for general time-dependent ODEs.

In addition to initial value problems, our approach can also address boundary value problems (BVPs). Given an oracle for preparing a state $\alpha|x(0)\rangle+\beta|x(T)\rangle$ expressing a general boundary condition, the goal of the quantum BVP is to produce a quantum state $\epsilon$-close to $|x(t)\rangle$ (a normalized state proportional to $x(t)$) for any desired $t\in[0,T]$. We also give a quantum algorithm for this problem with complexity $\operatorname{poly}(\log(1/\epsilon))$.

The remainder of this paper is organized as follows. Section 2 introduces the spectral method and Section 3 shows how to encode it into a quantum linear system. Then Section 4 analyzes the exponential decrease of the solution error, Section 5 bounds the condition number of the linear system, Section 6 lower bounds the success probability of the final measurement, and Section 7 describes how to prepare the initial quantum state. We combine these bounds in Section 8 to establish the main result. We then extend the analysis for initial value problems to boundary value problems in Section 9. Finally, we conclude in Section 10 with a discussion of the results and some open problems.

2 Spectral method

Spectral methods provide a way of solving differential equations using global approximations [19, 32]. The main idea of the approach is as follows. First, express an approximate solution as a linear combination of certain basis functions with undetermined coefficients. Second, construct a system of linear equations that such an approximate solution should satisfy. Finally, solve the linear system to determine the coefficients of the linear combination.

Spectral methods offer a flexible approach that can be adapted to different settings by careful choice of the basis functions and the linear system. A Fourier series provides an appropriate basis for periodic problems, whereas Chebyshev polynomials can be applied more generally. The linear system can be specified using Gaussian quadrature (giving a spectral element method or Tau method), or one can simply interpolate the differential equations using quadrature nodes (giving a pseudo-spectral method) [32]. Since general linear ODEs are non-periodic, and interpolation facilitates constructing a straightforward linear system, we develop a quantum algorithm based on the Chebyshev pseudo-spectral method [4, 22].

In this approach, we consider a truncated Chebyshev approximation $x(t)$ of the exact solution $\hat{x}(t)$, namely

[TABLE]

for any $n\in{\mathbb{Z}}^{+}$. (See Appendix A for the definition of $T_{k}(t)$ and a discussion of its properties.) The coefficients $c_{i,k}\in{\mathbb{C}}$ for all $i\in[{d}]_{0}$ and $k\in[{n+1}]_{0}$ are determined by demanding that $x(t)$ satisfies the ODE and initial conditions at a set of interpolation nodes $\{t_{l}\}_{l=0}^{n}$ (with $1=t_{0}>t_{1}>\cdots>t_{n}=-1$), where $x(t_{0})$ and $x(t_{n})$ are the initial and final states, respectively. In other words, we require

[TABLE]

and

[TABLE]

We would like to be able to increase the accuracy of the approximation by increasing $n$, so that

[TABLE]

There are many possible choices for the interpolation nodes. Here we use the Chebyshev-Gauss-Lobatto quadrature nodes, $t_{l}=\cos\frac{l\pi}{n}$ for $l\in[{n+1}]_{0}$, since these nodes achieve the highest convergence rate among all schemes with the same number of nodes [24, 26]. These nodes also have the convenient property that $T_{k}(t_{l})=\cos\frac{kl\pi}{n}$, making it easy to compute the values $x_{i}(t_{l})$.

To evaluate the condition (2.2), it is convenient to define coefficients $c^{\prime}_{i,k}$ for $i\in[{d}]_{0}$ and $k\in[{n+1}]_{0}$ such that

[TABLE]

We can use the differential property of Chebyshev polynomials,

[TABLE]

to determine the transformation between $c_{i,k}$ and $c^{\prime}_{i,k}$. As detailed in Appendix A, we have

[TABLE]

where $D_{n}$ is the $(n+1)\times(n+1)$ upper triangular matrix with nonzero entries

[TABLE]

where

[TABLE]

Using this expression in (2.2), (2.7), and (2.8), we obtain the following linear equations:

[TABLE]

We also demand that the Chebyshev series satisfies the initial condition $x_{i}(1)=\gamma_{i}$ for all $i\in[{d}]_{0}$. This system of linear equations gives a global approximation of the underlying system of differential equations. Instead of locally approximating the ODE at discretized times, these linear equations use the behavior of the differential equations at the $n+1$ times $\{t_{l}\}_{l=0}^{n}$ to capture their behavior over the entire interval $[-1,1]$.

Our algorithm solves this linear system using the high-precision QLSA [15]. Given an encoding of the Chebyshev coefficients $c_{ik}$, we can obtain the approximate solution $x(t)$ as a suitable linear combination of the $c_{ik}$, a computation that can also be captured within a linear system. The resulting approximate solution $x(t)$ is close to the exact solution $\hat{x}(t)$:

Lemma 1** (Lemma 19 of [19]).**

Let $\hat{x}(t)\in C^{r+1}(-1,1)$ be the solution of the differential equations (1.1) and let $x(t)$ satisfy (2.2) and (2.3) for $\{t_{l}=\cos\frac{l\pi}{n}\}_{l=0}^{n}$. Then there is a constant $C$, independent of $n$, such that

[TABLE]

This shows that the convergence behavior of the spectral method is related to the smoothness of the solution. For a solution in $C^{r+1}$, the spectral method approximates the solution with $n=\operatorname{poly}({1}/{\epsilon})$. Furthermore, if the solution is smoother, we have an even tighter bound:

Lemma 2** (Eq. (1.8.28) of [32]).**

Let $\hat{x}(t)\in C^{\infty}(-1,1)$ be the solution of the differential equations (1.1) and let $x(t)$ satisfy (2.2) and (2.3) for $\{t_{l}=\cos\frac{l\pi}{n}\}_{l=0}^{n}$. Then

[TABLE]

For simplicity, we replace the value $\sqrt{{2}/{\pi}}$ by the upper bound of $1$ in the following analysis.

This result implies that if the solution is in $C^{\infty}$, the spectral method approximates the solution to within $\epsilon$ using only $n=\operatorname{poly}(\log(1/\epsilon))$ terms in the Chebyshev series. Consequently, this approach gives a quantum algorithm with complexity $\operatorname{poly}(\log(1/\epsilon))$.

3 Linear system

In this section we construct a linear system that encodes the solution of a system of differential equations via the Chebyshev pseudospectral method introduced in Section 2. We consider a system of linear, first-order, time-dependent ordinary differential equations, and focus on the following initial value problem:

Problem 1**.**

In the quantum ODE problem, we are given a system of equations

[TABLE]

where $x(t)\in{\mathbb{C}}^{d}$, $A(t)\in{\mathbb{C}}^{d\times d}$ is $s$-sparse, and $f(t)\in{\mathbb{C}}^{d}$ for all $t\in[0,T]$. We assume that $A_{ij},f_{i}\in C^{\infty}(0,T)$ for all $i,j\in[{d}]$. We are also given an initial condition $x(0)=\gamma\in{\mathbb{C}}^{d}$. Given oracles that compute the locations and values of nonzero entries of $A(t)$ for any $t$, and that prepare normalized states $|\gamma\rangle$ proportional to $\gamma$ and $|f(t)\rangle$ proportional to $f(t)$ for any $t\in[0,T]$, the goal is to output a quantum state $|x(T)\rangle$ that is proportional to $x(T)$.

Without loss of generality, we rescale the interval $[0,T]$ onto $[-1,1]$ by the linear map $t\mapsto 1-2t/T$. Under this rescaling, we have $\frac{\mathrm{d}}{\mathrm{d}{t}}\mapsto-\frac{T}{2}\frac{\mathrm{d}}{\mathrm{d}{t}}$, so $A\mapsto-\frac{T}{2}A$, which can dramatically increase the spectral norm. To avoid this phenomenon, we divide the interval $[0,T]$ into subintervals $[0,\Gamma_{1}],[\Gamma_{1},\Gamma_{2}],\ldots,[\Gamma_{m-1},T]$ with $\Gamma_{0}:=0,\Gamma_{m}:=T$. Each subinterval $[\Gamma_{h},\Gamma_{h+1}]$ for $h\in[{m}]_{0}$ is then rescaled onto $[-1,1]$ with the linear map $K_{h}\colon[\Gamma_{h},\Gamma_{h+1}]\to[-1,1]$ defined by

[TABLE]

which satisfies $K_{h}(\Gamma_{h})=1$ and $K_{h}(\Gamma_{h+1})=-1$. To solve the overall initial value problem, we simply solve the differential equations for each successive interval (as encoded into a single system of linear equations).

Now let $\tau_{h}:=|\Gamma_{h+1}-\Gamma_{h}|$ and define

[TABLE]

Then, for each $h\in[{m}]_{0}$, we have the rescaled differential equations

[TABLE]

for $t\in[-1,1]$ with the initial conditions

[TABLE]

By taking

[TABLE]

where $\|{\cdot}\|$ denotes the spectral norm, we can ensure that $\|{A_{h}(t)}\|\leq 1$ for all $t\in[-1,1]$. In particular, it suffices to take

[TABLE]

Having rescaled the equations to use the domain $[-1,1]$, we now apply the Chebyshev pseudospectral method. Following Section 2, we substitute the truncated Chebyshev series of $x(t)$ into the differential equations with interpolating nodes $\{t_{l}=\cos\frac{l\pi}{n}:l\in[{n}]\}$, giving the linear system

[TABLE]

with initial condition

[TABLE]

Note that in the following, terms with $l=0$ refer to this initial condition.

We now describe a linear system

[TABLE]

that encodes the Chebyshev pseudospectral approximation and uses it to produce an approximation of the solution at time $T$.

The vector $|X\rangle\in{\mathbb{C}}^{m+p}\otimes{\mathbb{C}}^{d}\otimes{\mathbb{C}}^{n+1}$ represents the solution in the form

[TABLE]

where $c_{i,l}(\Gamma_{h+1})$ are the Chebyshev series coefficients of $x(\Gamma_{h+1})$ and $x_{i}:=x(\Gamma_{m})_{i}$ is the $i$th component of the final state $x(\Gamma_{m})$.

The right-hand-side vector $|B\rangle$ represents the input terms in the form

[TABLE]

where

[TABLE]

Here $\gamma$ is the initial condition and $f_{h}(\cos\frac{l\pi}{n})_{i}$ is $i$th component of $f_{h}$ at the interpolation point $t_{l}=\cos\frac{l\pi}{n}$.

We decompose the matrix $L$ in the form

[TABLE]

We now describe each of the matrices $L_{i}$ for $i\in[{5}]$ in turn.

The matrix $L_{1}$ is a discrete representation of $\frac{\mathrm{d}{x}}{\mathrm{d}{t}}$, satisfying

[TABLE]

(recall from (2.5) and (2.7) that $D_{n}$ encodes the action of the time derivative on a Chebyshev expansion). Thus $L_{1}$ has the form

[TABLE]

where the interpolation matrix is a discrete cosine transform matrix:

[TABLE]

The matrix $L_{2}(A_{h})$ discretizes $A_{h}(t)$, i.e.,

[TABLE]

Thus

[TABLE]

Note that if $A_{h}$ is time-independent, then

[TABLE]

The matrix $L_{3}$ combines the Chebyshev series coefficients $c_{i,l}$ to produce $x_{i}$ for each $i\in[{d}]_{0}$. To express the final state $x(-1)$, $L_{3}$ represents the linear combination $x_{i}(-1)=\sum_{k=0}^{n}c_{i,k}T_{k}(-1)=\sum_{k=0}^{n}(-1)^{k}c_{i,k}$. Thus we take

[TABLE]

Notice that $L_{3}$ has zero rows for $l\in[{n}]$.

When $h=m$, $L_{4}$ is used to construct $x_{i}$ from the output of $L_{3}$ for $l=0$, and to repeat $x_{i}$ $n$ times for $l\in[{n}]$. When $m+1\leq h\leq m+p$, both $L_{4}$ and $L_{5}$ are used to repeat $x_{i}$ $(n+1)p$ times for $l\in[{n}]$. This repetition serves to increase the success probability of the final measurement. In particular, we take

[TABLE]

and

[TABLE]

In summary, the linear system is as follows. For each $h\in[{m}]_{0}$, $(L_{1}+L_{2}(A_{h}))|X\rangle=|B_{h}\rangle$ solves the differential equations over $[\Gamma_{h},\Gamma_{h+1}]$, and the coefficients $c_{i,l}(\Gamma_{h+1})$ are combined by $L_{3}$ into the $(h+1)$st block as initial conditions. When $h=m$, the final coefficients $c_{i,l}(\Gamma_{m})$ are combined by $L_{3}$ and $L_{4}$ into the final state with coefficients $x_{i}$, and this solution is repeated $(p+1)(n+1)$ times by $L_{4}$ and $L_{5}$.

To explicitly illustrate the structure of this system, we present a simple example in Appendix B.

4 Solution error

In this section, we bound how well the solution of the linear system defined above approximates the actual solution of the system of differential equations.

Lemma 3**.**

For the linear system $L|X\rangle=|B\rangle$ defined in (3.12), let $x$ be the approximate ODE solution specified by the linear system and let $\hat{x}$ be the exact ODE solution. For all $t\in[0,T]$, assume $A(t)$ can be diagonalized as $A(t)=V(t)\Lambda(t)V^{-1}(t)$ for some $\Lambda(t)=\operatorname{diag}(\lambda_{0}(t),\ldots,\lambda_{d}(t))$ with $\mathop{\mathrm{Re}}(\lambda_{i}(t))\leq 0$ for all $i\in[{d}]_{0}$. Let $\kappa_{V}:=\max_{t}\kappa_{V}(t)$ be an upper bound of the condition number of $V(t)$. Then for $n$ sufficiently large, the error in the solution at time $T$ satisfies

[TABLE]

Proof.

First we carefully choose $n$ satisfying

[TABLE]

where

[TABLE]

to ensure that

[TABLE]

According to the quantum spectral method defined in Section 3, we solve

[TABLE]

We denote the exact solution by $\hat{x}(\Gamma_{h+1})$, and we let $x(\Gamma_{h+1})=\sum_{i=0}^{d}\sum_{l=0}^{n}(-1)^{n}c_{i,l}(\Gamma_{h+1})$, where $c_{i,l}(\Gamma_{h+1})$ is defined in (3.13). Define

[TABLE]

For $h=0$, Lemma 2 implies

[TABLE]

For $h\in[{m}]$, the error in the approximate solution of $\frac{\mathrm{d}{x}}{\mathrm{d}{t}}=A_{h}(t)x(t)+f_{h}(t)$ has two contributions: the error from the linear system and the error in the initial condition. We let $\widetilde{x}(\Gamma_{h+1})$ denote the solution of the linear system \bigl{(}L_{1}+L_{2}(A_{h})\bigr{)}|\widetilde{x}(\Gamma_{h+1})\rangle=|B(f_{h})\rangle under the initial condition $\hat{x}(\Gamma_{h})$. Then

[TABLE]

The first term can be bounded using Lemma 2, giving

[TABLE]

The second term comes from the initial error $\Delta_{h}$, which is transported through the linear system. Let

[TABLE]

where $E_{h+1}$ is the solution of the linear system with input $\Delta_{h}$ and $\hat{E}_{h+1}$ is the exact solution of $\frac{\mathrm{d}{x}}{\mathrm{d}{t}}=A_{h+1}(t)x(t)+f_{h+1}(t)$ with initial condition $x(\Gamma_{h})=\Delta_{h}$. Then by Lemma 2,

[TABLE]

so

[TABLE]

Thus, we have an inequality recurrence for bounding $\Delta_{h}$:

[TABLE]

Now we iterate $h$ from $1$ to $m$. Equation (4.4) implies

[TABLE]

so

[TABLE]

Therefore

[TABLE]

which shows that the solution error decreases exponentially with $n$. In other words, the quantum spectral method approximates the solution with error $\epsilon$ using $n=\operatorname{poly}(\log({1}/{\epsilon}))$. ∎

Note that for time-independent differential equations, we can directly estimate $\|\hat{x}^{(n+1)}(t)\|$ using

[TABLE]

Writing $A_{h}=V_{h}\Lambda_{h}V_{h}^{-1}$ where $\Lambda_{h}=\operatorname{diag}(\lambda_{0},\ldots,\lambda_{d-1})$, we have $e^{A_{h}}=V_{h}e^{\Lambda_{h}}V_{h}^{-1}$. Thus the exact solution of time-independent equation with initial condition $\hat{x}(1)=\gamma$ is

[TABLE]

Since $\mathop{\mathrm{Re}}(\lambda_{i})\leq 0$ for all eigenvalues $\lambda_{i}$ of $A_{h}$ for $i\in[{d}]_{0}$, we have $\|e^{\Lambda_{h}}\|\leq 1$. Therefore

[TABLE]

Furthermore, since $\max_{h,t}\|A_{h}(t)\|\leq 1$, we have

[TABLE]

Thus the solution error satisfies

[TABLE]

Note that, although we represent the solution differently, this bound is similar to the corresponding bound in [10, Theorem 6].

5 Condition number

We now analyze the condition number of the linear system.

Lemma 4**.**

Consider an instance of the quantum ODE problem as defined in Problem 1. For all $t\in[0,T]$, assume $A(t)$ can be diagonalized as $A(t)=V(t)\Lambda(t)V^{-1}(t)$ for some $\Lambda(t)=\operatorname{diag}(\lambda_{0}(t),\ldots,\lambda_{d}(t))$, with $\mathop{\mathrm{Re}}(\lambda_{i}(t))\leq 0$ for all $i\in[{d}]_{0}$. Let $\kappa_{V}:=\max_{t\in[0,T]}\kappa_{V}(t)$ be an upper bound on the condition number of $V(t)$. Then for $m,p\in{\mathbb{Z}}^{+}$ and $n$ sufficiently large, the condition number of $L$ in the linear system (3.12) satisfies

[TABLE]

Proof.

We begin by bounding the norms of some operators that appear in the definition of $L$. First we consider the $l_{\infty}$ norm of $D_{n}$ since this is straightforward to calculate:

[TABLE]

Thus we have the upper bound

[TABLE]

Next we upper bound the spectral norm of the discrete cosine transform matrix $P_{n}$:

[TABLE]

Therefore

[TABLE]

Thus we can upper bound the norm of $L_{1}$ as

[TABLE]

Next we consider the spectral norm of $L_{2}(A_{h})$ for any $h\in[{m}]_{0}$. We have

[TABLE]

Since the eigenvalues of each $A_{h}(t_{l})$ for $l\in[{n+1}]_{0}$ are all eigenvalues of

[TABLE]

we have

[TABLE]

by (3.8). Therefore

[TABLE]

By direct calculation, we have

[TABLE]

Thus, for $n\geq 5$, we find

[TABLE]

Next we upper bound $\|L^{-1}\|$. By definition,

[TABLE]

We express $|B\rangle$ as

[TABLE]

where $|b_{hl}\rangle:=\sum_{i=0}^{d-1}\beta_{hil}|hil\rangle$ satisfies $\||b_{hl}\rangle\|^{2}=\sum_{i=0}^{d-1}|\beta_{hil}|^{2}\leq 1$. For any fixed $h\in[{m+p+1}]_{0}$ and $l\in[{n+1}]_{0}$, we first upper bound $\|L^{-1}|b_{hl}\rangle\|$ and use this to upper bound the norm of $L^{-1}$ applied to linear combinations of such vectors.

Recall that the linear system comes from (2.10), which is equivalent to

[TABLE]

For fixed $h\in[{m+p+1}]_{0}$ and $r\in[{n+1}]_{0}$, define vectors $x_{hr},x_{hr}^{\prime}\in{\mathbb{C}}^{d}$ with

[TABLE]

for $i\in[{d}]_{0}$. We claim that $x_{hr}=x_{hr}^{\prime}=0$ for any $r\neq l$. Combining only the equations from (5.17) with $r\neq l$ gives the system

[TABLE]

Consider a corresponding system of differential equations

[TABLE]

where $\hat{x}_{hr}(t)\in{\mathbb{C}}^{d}$ for all $t\in[-1,1]$. The solution of this system with $b=0$ and initial condition $\hat{x}_{hr}(1)=0$ is clearly $\hat{x}_{hr}(t)=0$ for all $t\in[-1,1]$. Then the $n$th-order truncated Chebyshev approximation of (5.20), which should satisfy the linear system (5.19) by (2.1) and (2.2), is exactly $x_{hr}$. Using Lemma 3 and observing that $\hat{x}^{(n+1)}(t)=0$, we have

[TABLE]

When $t=t_{l}$, we let $|B\rangle=|b_{hl}\rangle$ denote the first nonzero vector. Combining only the equations from (5.17) with $r=l$ gives the system

[TABLE]

Consider a corresponding system of differential equations

[TABLE]

with $\gamma=b_{h0},b=0$ for $l=0$; or $\gamma=0,b=b_{hl}$ for $l\in[{n}]$.

Using the diagonalization $A_{h}(t_{l})=V(t_{l})\Lambda_{h}(t_{l})V^{-1}(t_{l})$, we have $e^{A}=V(t_{l})e^{\Lambda_{h}(t_{l})}V^{-1}(t_{l})$. Thus the exact solution of the differential equations (5.20) with $r=l$ and initial condition $\hat{x}_{hr}(1)=\gamma$ is

[TABLE]

According to equation (4.4) in the proof of Lemma 3, we have

[TABLE]

where

[TABLE]

Now for $h\in[{m+1}]_{0}$, we take $x_{hl}$ to be the initial condition $\gamma$ for the next subinterval to obtain $x_{(h+1)l}$. Using (5.24) and (5.25), starting from $\gamma=b_{h0},b=0$ for $l=0$, we find

[TABLE]

Since $\|\Lambda_{h}(t_{l})\|\leq\|\Lambda\|\leq 1$ and $\Lambda_{h}(t_{l})=\operatorname{diag}(\lambda_{0},\ldots,\lambda_{d-1})$ with $\mathop{\mathrm{Re}}(\lambda_{i})\leq 0$ for $i\in[{d}]_{0}$, we have $\|e^{2\Lambda_{h}(t_{l})}\|\leq 1$. Therefore

[TABLE]

On the other hand, with $\gamma=0,b=b_{hl}$ for $l\in[{n}]$, we have

[TABLE]

so

[TABLE]

For $h\in\{m,m+1,\ldots,m+p\}$, according to the definition of $L_{4}$ and $L_{5}$, we similarly have

[TABLE]

According to (5.24), $\hat{x}_{hl}(t)$ is a monotonic function of $t\in[-1,1]$, which implies

[TABLE]

Using the identity

[TABLE]

we have

[TABLE]

Consider the Chebyshev expansion of $\hat{x}_{hl}(t)$ as in (2.1):

[TABLE]

By the orthogonality of Chebyshev polynomials (as specified in (A.7)), we have

[TABLE]

Using (5.34), this gives

[TABLE]

Now we compute $\||X\rangle\|$, summing the contributions from all $c_{i,r}(\Gamma_{h})$ and $x_{mr}$, and notice that $c_{i,r}=0$ and $x_{mr}=0$ for all $r\neq l$, giving

[TABLE]

Finally, considering all $h\in[{m+p+1}]_{0}$ and $l\in[{n+1}]_{0}$, from (5.16) we have

[TABLE]

so

[TABLE]

and therefore

[TABLE]

Finally, combining (5.14) and (5.41) gives

[TABLE]

as claimed. ∎

6 Success probability

We now evaluate the success probability of our approach to the quantum ODE problem.

Lemma 5**.**

Consider an instance of the quantum ODE problem as defined in Problem 1 with the exact solution $\hat{x}(t)$ for $t\in[0,T]$, and its corresponding linear system (3.12) with $m,p\in{\mathbb{Z}}^{+}$ and $n$ sufficiently large. When applying the QLSA to this system, the probability of measuring a state proportional to $|x(T)\rangle=\sum_{i=0}^{d-1}x_{i}|i\rangle$ is

[TABLE]

where $x_{i}$ is defined in (3.13), $\tau$ is defined in (3.9), and

[TABLE]

Proof.

After solving the linear system (3.12) using the QLSA, we measure the first and third registers of $|X\rangle$ (as defined in (3.13)). We decompose this state as

[TABLE]

where

[TABLE]

When the first register is observed in some $h\in\{m,m+1,\ldots,m+p\}$ (no matter what outcome is seen for the third register), we output the second register, which is then in a normalized state proportional to the final state:

[TABLE]

with

[TABLE]

Notice that

[TABLE]

and

[TABLE]

Considering the definition of $q$, the contribution from time interval $h$ under the rescaling (3.2), and the identity (5.33), we have

[TABLE]

where $\hat{x}_{h}(\overline{t})$ is the solution of (4.5) with the rescaling in (3.4). By the orthogonality of Chebyshev polynomials (as specified in (A.7)),

[TABLE]

For all $h\in[{m}]_{0}$, we have

[TABLE]

and therefore

[TABLE]

Thus we see that the success probability of the measurement satisfies

[TABLE]

as claimed. ∎

7 State preparation

We now describe a procedure for preparing the vector $|B\rangle$ in the linear system (3.12) (defined in (3.14) and (3.15)) using the given ability to prepare the initial state of the system of differential equations. We also evaluate the complexity of this procedure.

Lemma 6**.**

Consider state preparation oracles acting on a state space with basis vectors $|h\rangle|i\rangle|l\rangle$ for $h\in[{m}]_{0},i\in[{d}]_{0},l\in[{n}]_{0}$, where $m,d,n\in{\mathbb{N}}$, encoding an initial condition $\gamma\in{\mathbb{C}}^{d}$ and function $f_{h}(\cos\frac{l\pi}{n})\in{\mathbb{C}}^{d}$ as in (3.15). Specifically, for any $h\in[{m}]_{0}$ and $l\in[{n}]$, let $O_{x}$ be a unitary oracle that maps $|0\rangle|0\rangle|0\rangle$ to a state proportional to $|0\rangle|\gamma\rangle|0\rangle$ and $|h\rangle|\phi\rangle|l\rangle$ to $|h\rangle|\phi\rangle|l\rangle$ for any $|\phi\rangle$ orthogonal to $|0\rangle$; let $O_{f}(h,l)$ be a unitary that maps $|h\rangle|0\rangle|l\rangle$ to a state proportional to $|h\rangle|f_{h}(\cos\frac{l\pi}{n})\rangle|l\rangle$ and maps $|0\rangle|\phi\rangle|0\rangle$ to $|0\rangle|\phi\rangle|0\rangle$ for and $|\phi\rangle$ orthogonal to $|0\rangle$. Suppose $\|\gamma\|$ and $\|f_{h}(\cos\frac{l\pi}{n})\|$ are known. Then the normalized quantum state

[TABLE]

can be prepared with gate and query complexity $O(mn)$.

Proof.

We normalize the components of the state using the coefficients

[TABLE]

so that

[TABLE]

First we perform a unitary transformation mapping

[TABLE]

This can be done in time complexity $O(mn)$ by standard techniques [33]. Then we perform $O_{x}$ and $O_{f}(h,l)$ for all $h\in[{m}]_{0},l\in[{n}]$, giving

[TABLE]

using $O(mn)$ queries. ∎

8 Main result

Having analyzed the solution error, condition number, success probability, and state preparation procedure for our approach, we are now ready to establish the main result.

Theorem 1**.**

Consider an instance of the quantum ODE problem as defined in Problem 1. Assume $A(t)$ can be diagonalized as $A(t)=V(t)\Lambda(t)V^{-1}(t)$ where $\Lambda(t)=\operatorname{diag}(\lambda_{1}(t),\ldots,\lambda_{d}(t))$ with $\mathop{\mathrm{Re}}(\lambda_{i}(t))\leq 0$ for each $i\in[{d}]_{0}$ and $t\in[0,T]$. Then there exists a quantum algorithm that produces a state $x(T)/\|x(T)\|$ $\epsilon$-close to $\hat{x}(T)/\|\hat{x}(T)\|$ in $l^{2}$ norm, succeeding with probability $\Omega(1)$, with a flag indicating success, using

[TABLE]

queries to oracles $O_{A}(h,l)$ (a sparse matrix oracle for $A_{h}(t_{l})$ as defined in (3.3)) and $O_{x}$ and $O_{f}(h,l)$ (as defined in Lemma 6). Here $\|A\|:=\max_{t\in[0,T]}\|A(t)\|$; $\kappa_{V}:=\max_{t}\kappa_{V}(t)$, where $\kappa_{V}(t)$ is the condition number of $V(t)$; and

[TABLE]

The gate complexity is larger than the query complexity by a factor of $\operatorname{poly}(\log(\kappa_{V}ds\|A\|g^{\prime}T/\epsilon))$.

Proof.

We first present the algorithm and then analyze its complexity.

Statement of the algorithm. First, we choose $m$ to guarantee

[TABLE]

Then, as in Section 3, we divide the interval $[0,T]$ into small subintervals $[0,\Gamma_{1}],[\Gamma_{1},\Gamma_{2}],\ldots,[\Gamma_{m-1},T]$ with $\Gamma_{0}=0,\Gamma_{m}=T$, and define

[TABLE]

Each subinterval $[\Gamma_{h},\Gamma_{h+1}]$ for $h\in[{m-1}]$ is mapped onto $[-1,1]$ with a linear mapping $K_{h}$ satisfying $K_{h}(\Gamma_{h})=1,K_{h}(\Gamma_{h+1})=-1$:

[TABLE]

We choose

[TABLE]

where

[TABLE]

and

[TABLE]

Since $\max_{t\in[0,T]}\|\hat{x}^{(n+1)}(t)\|\leq g^{\prime}$ and Lemma 3, this choice guarantees

[TABLE]

and

[TABLE]

Now $\|\hat{x}(T)-x(T)\|\leq\delta$ implies

[TABLE]

so we can choose such $n$ to ensure that the normalized output state is $\epsilon$-close to $\hat{x}(T)/\|\hat{x}(T)\|$.

Following Section 3, we build the linear system $L|X\rangle=|B\rangle$ (see (3.12)) that encodes the quantum spectral method. By Lemma 4, the condition number of this linear system is at most $(\pi m+p+1)(n+1)^{3.5}(2\kappa_{V}+e\|\gamma\|_{\infty})$. Then we use the QLSA from reference [15] to obtain a normalized state $|X\rangle$ and measure the first and third register of $|X\rangle$ in the standard basis. If the measurement outcome for the first register belongs to

[TABLE]

we output the state of the second register, which is a normalized state $|x(T)\rangle/\||x(T)\rangle\|$ satisfying (8.11). By Lemma 5, the probability of this event happening is at least $\frac{(p+1)(n+1)}{\pi mq^{2}+(p+1)(n+1)}$. To ensure $m+p=O(\|A\|T)$, we can choose

[TABLE]

so we can achieve success probability $\Omega(1)$ with $O(q/\sqrt{n})$ repetitions of the above procedure.

Analysis of the complexity. The matrix $L$ is an $(m+p+1)d(n+1)\times(m+p+1)d(n+1)$ matrix with $O(ns)$ nonzero entries in any row or column. By Lemma 4 and our choice of parameters, the condition number of $L$ is O\bigl{(}\kappa_{V}(m+p)n^{3.5}\bigr{)}. Consequently, by Theorem 5 of [15], the QLSA produces the state $|x(T)\rangle$ with

[TABLE]

queries to the oracles $O_{A}(h,l)$, $O_{x}$, and $O_{f}(h,l)$, and its gate complexity is larger by a factor of $\operatorname{poly}(\log(\kappa_{V}mnds/\delta))$. Using $O(q/\sqrt{n})$ steps of amplitude amplification to achieve success probability $\Omega(1)$, the overall query complexity of our algorithm is

[TABLE]

and the gate complexity is larger by a factor of

[TABLE]

as claimed. ∎

In general, $g^{\prime}$ could be unbounded above as $n\rightarrow\infty$. However, we could obtain a useful bound in such a case by solving the implicit equations (8.9) and (8.10).

Note that for time-independent differential equations, we can replace $g^{\prime}$ by $\|\gamma\|+2\tau\|f\|$ as shown in (4.20). In place of (8.7) and (8.8), we choose

[TABLE]

and

[TABLE]

By Lemma 3, this choice guarantees

[TABLE]

and

[TABLE]

Thus we have the following:

Corollary 1**.**

For time-independent differential equations, under the same assumptions of Theorem 1, there exists a quantum algorithm using

[TABLE]

queries to $O_{A}(h,l)$, $O_{x}$, and $O_{f}(h,l)$. The gate complexity of this algorithm is larger than its query complexity by a factor of $\operatorname{poly}(\log(\kappa_{V}ds\gamma\|A\|\|f\|T/\epsilon))$.

The complexity of our algorithm depends on the parameter $q$ in defined in (8.2), which characterizes the decay of the final state relative to the initial state. As discussed in Section 8 of [10], it is unlikely that the dependence on $q$ can be significantly improved, since renormalization of the state effectively implements postselection and an efficient procedure for performing this would have the unlikely consequence $\mathrm{BQP}=\mathrm{PP}$.

We also require the real parts of the eigenvalues of $A(t)$ to be non-positive for all $t\in[0,T]$ so that the solution cannot grow exponentially. This requirement is essentially the same as in the time-independent case considered in [10] and improves upon the analogous condition in [5] (which requires an additional stability condition). Also as in [10], our algorithm can produce approximate solutions for non-diagonalizable $A(t)$, although the dependence on $\epsilon$ degrades to $\operatorname{poly}(1/\epsilon)$. For further discussion of these considerations, see Sections 1 and 8 of [10].

9 Boundary value problems

So far we have focused on initial value problems (IVPs). Boundary value problems (BVPs) are another widely studied class of differential equations that appear in many applications, but that can be harder to solve than IVPs.

Consider a sparse, linear, time-dependent system of differential equations as in Problem 1 but with a constraint on some linear combination of the initial and final states:

Problem 2**.**

In the quantum BVP, we are given a system of equations

[TABLE]

where $x(t)\in{\mathbb{C}}^{d}$, $A(t)\in{\mathbb{C}}^{d\times d}$ is $s$-sparse, and $f(t)\in{\mathbb{C}}^{d}$ for all $t\in[0,T]$, and a boundary condition $\alpha x(0)+\beta x(T)=\gamma$ with $\alpha,\beta,\gamma\in{\mathbb{C}}^{d}$. Suppose there exists a unique solution $\hat{x}\in C^{\infty}(0,T)$ of this boundary value problem. Given oracles that compute the locations and values of nonzero entries of $A(t)$ for any $t$, and that prepare quantum states $\alpha|x(0)\rangle+\beta|x(T)\rangle=|\gamma\rangle$ and $|f(t)\rangle$ for any $t$, the goal is to output a quantum state $|x(t^{*})\rangle$ that is proportional to $x(t^{*})$ for some specified $t^{*}\in[0,T]$.

As before, we can rescale $[0,T]$ onto $[-1,1]$ by a linear mapping. However, since we have boundary conditions at $t=0$ and $t=T$, we cannot divide $[0,T]$ into small subintervals. Instead, we directly map $[0,T]$ onto $[-1,1]$ with a linear map $K$ satisfying $K(0)=1$ and $K(T)=-1$:

[TABLE]

Now the new differential equations are

[TABLE]

If we define $A_{K}(\overline{t}):=-\frac{T}{2}A(t)$ and $f_{K}(\overline{t})=-\frac{T}{2}f(t)$, we have

[TABLE]

for $\overline{t}\in[-1,1]$. Now the boundary condition takes the form

[TABLE]

Since we only have one solution interval, we need to choose a larger order $n$ of the Chebyshev series to reduce the solution error. In particular, we take

[TABLE]

where $\Omega$ and $\omega$ are the same as Theorem 1.

As in Section 3, we approximate $x(t)$ by a finite Chebyshev series with interpolating nodes $\{t_{l}=\cos\frac{l\pi}{n}:l\in[{n}]\}$ and thereby obtain a linear system

[TABLE]

with the boundary condition

[TABLE]

Observe that the linear equations have the same form as in (3.10). Instead of (3.11), the term with $l=0$ encodes the condition (9.8) expanded in a Chebyshev series, namely

[TABLE]

for each $i\in[{d}]_{0}$. Since $T_{k}(t_{0})=1$ and $T_{k}(t_{n})=(-1)^{k}$, this can be simplified as

[TABLE]

If $\alpha_{i}+(-1)^{k}\beta_{i}=0$, the element of $|il\rangle\langle ik|$ of $L_{2}(A_{K})$ is zero; if $\alpha_{i}+(-1)^{k}\beta_{i}\neq 0$, without loss of generality, the two sides of this equality can be divided by $\alpha_{i}+(-1)^{k}\beta_{i}$ to guarantee that the terms with $l=0$ can be encoded as in (3.11).

Now this system can be written in the form of equation (3.12) with $m=1$. Here $L$, $|X\rangle$, and $|B\rangle$ are the same as in (3.16), (3.13), and (3.14), respectively, with $m=1$, except for adjustments to $L_{3}$ that we now describe.

The matrix $L_{3}$ represents the linear combination $x_{i}(t^{*})=\sum_{k=0}^{n}c_{i,k}T_{k}(t^{*})$. Thus we take

[TABLE]

Since $|T_{k}(t^{*})|\leq 1$, we have

[TABLE]

and it follows that Lemma 3 also holds for boundary value problems. Similarly, Lemma 4 still holds with $m=1$.

We are now ready to analyze the complexity of the quantum BVP algorithm. The matrix $L$ defined above is a $(p+2)d(n+1)\times(p+2)d(n+1)$ matrix with $O(ns)$ nonzero entries in any row or column, with condition number $O(\kappa_{V}pn^{3.5})$. By Lemma 5 with $p=O(1)$, $O(q/\sqrt{n})$ repetitions suffice to ensure success probability $\Omega(1)$. By (9.6), $n$ is linear in $\|A\|T$ and poly-logarithmic in $\Omega$ and $\omega$. Therefore, we have the following:

Theorem 2**.**

Consider an instance of the quantum BVP as defined in Problem 2. Assume $A(t)$ can be diagonalized as $A(t)=V(t)\Lambda(t)V^{-1}(t)$ where $\Lambda(t)=\operatorname{diag}(\lambda_{1}(t),\ldots,\lambda_{d}(t))$ with $\mathop{\mathrm{Re}}(\lambda_{i}(t))\leq 0$ for each $i\in[{d}]_{0}$ and $t\in[0,T]$. Then there exists a quantum algorithm that produces a state $x(t^{*})/\|x(t^{*})\|$ $\epsilon$-close to $\hat{x}(t^{*})/\|\hat{x}(t^{*})\|$ in $l^{2}$ norm, succeeding with probability $\Omega(1)$, with a flag indicating success, using

[TABLE]

queries to $O_{A}(h,l)$, $O_{x}$, and $O_{f}(h,l)$. Here $\|A\|$, $\kappa_{V}$, $g$, $g^{\prime}$ and $q$ are defined as in Theorem 1. The gate complexity is larger than the query complexity by a factor of $\operatorname{poly}(\log(\kappa_{V}ds\|A\|g^{\prime}T/\epsilon))$.

As for initial value problems, we can simplify this result in the time-independent case.

Corollary 2**.**

For a time-independent boundary value problem, under the same assumptions of Theorem 2, there exists a quantum algorithm using

[TABLE]

queries to $O_{A}(h,l)$, $O_{x}$, and $O_{f}(h,l)$. The gate complexity of this algorithm is larger than its query complexity by a factor of $\operatorname{poly}(\log(\kappa_{V}ds\gamma\|A\|\|f\|T/\epsilon))$.

10 Discussion

In this paper, we presented a quantum algorithm to solve linear, time-dependent ordinary differential equations. Specifically, we showed how to employ a global approximation based on the spectral method as an alternative to the more straightforward finite difference method. Our algorithm handles time-independent differential equations with almost the same complexity as [10], but unlike that approach, can also handle time-dependent differential equations. Compared to [5], our algorithm improves the complexity of solving time-dependent linear differential equations from $\operatorname{poly}(1/\epsilon)$ to $\operatorname{poly}(\log(1/\epsilon))$.

This work raises several natural open problems. First, our algorithm must assume that the solution is smooth. If the solution is in $C^{r}$, the solution error is $O(\frac{1}{n^{r-2}})$ by Lemma 1. Can we improve the complexity to $\operatorname{poly}(\log(1/\epsilon))$ under such weaker smoothness assumptions?

Second, the complexity of our algorithm is logarithmic in the parameter $g^{\prime}$ defined in (8.2), which characterizes the amount of fluctuation in the solution. However, the query complexity of Hamiltonian simulation is independent of that parameter [31, 7]. Can we develop quantum algorithms for general differential equations with query complexity independent of $g^{\prime}$?

Third, our algorithm has nearly optimal dependence on $T$, scaling as $O(T\operatorname{poly}(\log T))$. According to the no-fast-forwarding theorem [6], the complexity must be at least linear in $T$, and indeed linear complexity is achievable for the case of Hamiltonian simulation [14]. Can we handle general differential equations with complexity linear in $T$? Furthermore, can we achieve an optimal tradeoff between $T$ and $\epsilon$ as shown for Hamiltonian simulation in [29]?

Finally, can the techniques developed here be applied to give improved quantum algorithms for linear partial differential equations, or even for nonlinear ODEs or PDEs?

Acknowledgements

We thank Stephen Jordan and Aaron Ostrander for valuable discussions of quantum algorithms for time-dependent linear differential equations.

This work was supported in part by the Army Research Office (MURI award W911NF-16-1-0349), the Canadian Institute for Advanced Research, the National Science Foundation (grants 1526380 and 1813814), and the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Quantum Algorithms Teams and Quantum Testbed Pathfinder programs.

Appendix A Chebyshev polynomials

This appendix defines the Chebyshev polynomials and presents some of their properties that are useful for our analysis.

For any $k\in{\mathbb{N}}$, the Chebyshev polynomial of the first kind can be defined as the function

[TABLE]

It can be shown that this is a polynomial of degree $k$ in $x$. For example, we have

[TABLE]

Using the trigonometric addition formula $\cos(k+1)\theta+\cos(k-1)\theta=2\cos\theta\cos k\theta$, we have the recurrence

[TABLE]

(which also provides an alternative definition of the Chebyshev polynomials, starting from the initial conditions $T_{0}(x)=1$ and $T_{1}(x)=x$). We also have the bounds

[TABLE]

Chebyshev polynomials are orthogonal polynomials on $[-1,1]$ with the weight function $w(x):=(1-x^{2})^{-1/2}$. More concretely, defining an inner product on $L_{w}^{2}(-1,1)$ by

[TABLE]

we have

[TABLE]

where

[TABLE]

It is well known from the approximation theorem of Weierstrass that $\{T_{k}(x):k\in{\mathbb{N}}\}$ is complete on the space $L_{w}^{2}(-1,1)$. In other words, we have the following:

Lemma 7**.**

Any function $u\in L_{w}^{2}(-1,1)$ can be expanded by a unique Chebyshev series as

[TABLE]

where the coefficients are

[TABLE]

For any $N\in{\mathbb{N}}$, we introduce the orthogonal projection $P_{N}\colon L_{w}^{2}(-1,1)\rightarrow\mathbb{P}_{N}$ (where $\mathbb{P}_{N}$ denotes the set of polynomials of degree at most $N$) by

[TABLE]

By the completeness of the Chebyshev polynomials, we have

[TABLE]

and

[TABLE]

Finally, we compute the Chebyshev series of $u^{\prime}(x)$ in terms of the Chebyshev series of $u(x)$. Since $T_{k}(x)=\cos k\theta$ where $\theta=\arccos x$, we have

[TABLE]

Since

[TABLE]

we obtain

[TABLE]

and

[TABLE]

Since $P_{N}u(x)\in\mathbb{P}_{N}$, the derivative of this projection should be in $\mathbb{P}_{N-1}$. Indeed, we have

[TABLE]

Comparing the coefficients of both sides, we find

[TABLE]

where $\sigma_{k}$ is defined in (A.8).

Since $\hat{c}^{\prime}_{k}=0$ for $k\geq N$, we can calculate $\hat{c}^{\prime}_{N-1}$ from $\hat{c}_{N}$ and then successively calculate $\hat{c}^{\prime}_{N-2},\ldots,\hat{c}^{\prime}_{1},\hat{c}^{\prime}_{0}$. This recurrence gives

[TABLE]

Since $\hat{c}^{\prime}_{k}$ only depends on $\hat{c}_{j}$ for $j>k$, the transformation matrix $D_{N}$ between the values $\hat{c}^{\prime}_{k}$ and $\hat{c}_{k}$ for $k\in[{n+1}]_{0}$ is an upper triangular matrix with all zero diagonal elements, namely

[TABLE]

Appendix B An example of the quantum spectral method

Section 3 defines a linear system that implements the quantum spectral method for solving a system of $d$ time-dependent differential equations. Here we present a simple example of this system for the case $d=1$, namely

[TABLE]

where $x(t),A(t),f(t)\in{\mathbb{C}}$, $t\in[0,T]$, and we have the initial condition

[TABLE]

In particular, we choose $m=3$, $n=2$, and $p=1$ in the specification of the linear system. We divide $[0,T]$ into $m=3$ intervals $[0,\Gamma_{1}],[\Gamma_{1},\Gamma_{2}],[\Gamma_{2},T]$ with $\Gamma_{0}=0,\Gamma_{m}=T$, and map each one onto $[-1,1]$ with the linear mapping $K_{h}$ satisfying $K_{h}(\Gamma_{h})=1$ and $K_{h}(\Gamma_{h+1})=-1$. Then we take the finite Chebyshev series of $x(t)$ with $n=2$ into the differential equation with interpolating nodes $\{t_{l}=\cos\frac{l\pi}{n}:l\in[{2}]\}=\{0,-1\}$ to obtain a linear system. Finally, we repeat the final state $p=1$ time to increase the success probability.

With these choices, the linear system has the form

[TABLE]

with

[TABLE]

The vector $|X\rangle$ has the form

[TABLE]

where $c_{l}(\Gamma_{h+1})$ are the Chebyshev series coefficients of $x(\Gamma_{h+1})$ and $x$ is the final state $x(\Gamma_{m})=x(-1)$.

Finally, the vector $|B\rangle$ has the form

[TABLE]

where $\gamma$ comes from the initial condition and $f_{h}(\cos\frac{l\pi}{n})$ is the value of $f_{h}$ at the interpolation point $t_{l}=\cos\frac{l\pi}{n}\in\{0,-1\}$.