High-precision quantum algorithms for partial differential equations

Quantum computers can produce a quantum encoding of the solution of a system of differential equations exponentially faster than a classical algorithm can produce an explicit description. However, while high-precision quantum algorithms for linear ordinary differential equations are well established, the best previous quantum algorithms for linear partial differential equations (PDEs) have complexity $\mathrm{poly}(1/\epsilon)$, where $\epsilon$ is the error tolerance. By developing quantum algorithms based on adaptive-order finite difference methods and spectral methods, we improve the complexity of quantum algorithms for linear PDEs to be $\mathrm{poly}(d, \log(1/\epsilon))$, where $d$ is the spatial dimension. Our algorithms apply high-precision quantum linear system algorithms to systems whose condition numbers and approximation errors we bound. We develop a finite difference algorithm for the Poisson equation and a spectral algorithm for more general second-order elliptic equations.


Introduction
Many scientific problems involve partial differential equations (PDEs). Prominent examples include Maxwell's equations for electromagnetism, Boltzmann's equation and the Fokker-Planck equation in thermodynamics, and Schrödinger's equation in continuum quantum mechanics. While models of physics are often studied in a constant number of spatial dimensions, it is also natural to study high-dimensional PDEs, such as to model systems with many interacting particles. Classical numerical methods have complexity that grows exponentially in the dimension, a phenomenon sometimes called the curse of dimensionality [2]. This is a major challenge for attempts to solve PDEs on classical computers.
A common approach to solving PDEs on a digital computer is the finite difference method (FDM). In this approach, we discretize space into a rectangular lattice, solve a system of linear equations that approximates the PDE on the lattice, and output the solution on those grid points. If each spatial coordinate has n discrete values, then n d points are needed to discretize a d-dimensional problem. Simply outputting the solution on these grid points takes time Ω(n d ).
Beyond uniform grids, the sparse grid technique [32] has been applied to reduce the time and space complexity of outputting a sparse encoding of the solution to O(n log d n) [6,39]. While this is a significant improvement, it still scales exponentially in d. It can be shown that for a grid-based approach this complexity is optimal with respect to certain norms [6]. Reference [6] proposes alternative sparse grid algorithms whose complexities scale linearly with n but exponentially with d. Another grid-based method is the finite element method (FEM), where the differential equation is multiplied by functions with local support (restricted by the grid) and then integrated. This produces a set of equations that the solution must satisfy, which are then used to approximate the solution. In yet another grid-based approach, the finite volume method (FVM) considers a grid dividing space into volumes/cells. The field is integrated over these volumes to create auxiliary variables, and relations between these variables are derived from the differential equation.
An alternative to grid methods is the concept of spectral methods [15,28]. Spectral methods use linear combinations of basis functions (such as Fourier basis states or Chebyshev polynomials) to globally approximate the solution. These basis functions allow the construction of a linear system whose solution approximates the solution of the PDE.
These classical algorithms often consider the problem of outputting the solution at N points in space, which clearly requires Ω(N ) space and time. Quantum algorithms often (though not always) consider the alternative problem of outputting a quantum state proportional to such a vector, which requires only Ω(log N ) space-and correspondingly provides more limited access to the solution-but can potentially be done in only poly(log N ) time.
The fact that quantum states can efficiently encode exponentially long vectors has also been leveraged for the development of quantum linear system algorithms (QLSAs) [1,8,16]. For a linear system A x = b, a QLSA outputs a quantum state proportional to the solution x. To learn information about the solution x, the output of the QLSA must be post-processed. For example, to output all the entries of an N -dimensional vector x given a quantum state |x proportional to it, even a quantum computer needs time and space Ω(N ).
Because linear systems are often used in classical algorithms for PDEs such as those described above, it is natural to consider their quantum counterparts. Clader, Jacobs, and Sprouse [10] give a heuristic algorithm for using sparse preconditioners and QLSAs to solve a linear system constructed using the FEM for Maxwell's equations. The state output by the QLSA is then post-processed to compute electromagnetic scattering cross-sections.
In subsequent work, Montanaro and Pallister [22] use QLSAs to implement the FEM for d-dimensional boundary value problems and evaluate the quantum speedup that can be achieved when estimating a function of the solution within precision . This involves a careful analysis of how different algorithmic parameters (such as the dimension and condition number of the FEM linear system and the number of post-processing measurements) scale with respect to input variables (such as the spatial dimension d and desired precision ), since all of these affect the complexity. Their algorithms have complexity poly(d, 1/ ), compared to O((1/ ) d ) for the classical FEM. This exponential improvement with respect to d suggests that quantum algorithms may be notably faster when d is large. However, they also argue that for fixed d, at most a polynomial speed-up can be expected due to lower bounds on the cost of post-processing the state to estimate a function of the solution. The FDM has also been used in quantum algorithms for PDEs. References [7,35] apply the FDM to solve Poisson's equation in rectangular volumes under Dirichlet boundary conditions. Although the circuits they construct have poly(log(1/ )) gates, these circuits have success probability poly(1/ ), leading to poly(1/ ) time complexity. Additionally, they do not quantify errors resulting from the finite-difference approximation. Reference [12] applies the FDM to the problem of outputting states proportional to solutions of the wave equation, giving complexity d is poly(n) for a fixed-order FDM). The FVM is combined with the reservoir method in Reference [14] to simulate hyperbolic equations; although they achieve linear scaling with respect to the spatial dimension, they use fixed order differences, leading to poly(1/ ) scaling. These FDM, FEM, and FVM approaches can only give a total complexity poly(1/ ), even using high-precision methods for the QLSA or Hamiltonian simulation, because of the additional approximation errors in the FDM, FEM, and FVM.
The FDM is also applied in Reference [18] to simulate how a fixed number of particles evolve under the Schrödinger equation with access to an oracle for the potential term. This can be seen as a special case of quantum algorithms for PDEs. Other examples include quantum algorithms for many-body quantum dynamics [37,38] and for electronic structure problems, including for quantum chemistry (see for example References [20,25]). However, here we focus on PDEs whose dynamics are not necessarily unitary.
In this paper, we propose new quantum algorithms for linear PDEs where the boundary is the unit hypercube. In the spirit of Reference [22], we state our results in terms of the approximation error and the spatial dimension; however, we do not consider the problem of estimating a function of the PDE solution and instead focus on outputting states encoding the solution, allowing us to give algorithms with complexity poly(log(1/ )). Just as for the QLSA, this improvement is potentially significant if the given equations must be solved as a subroutine within some larger computation. The problem we address can be informally stated as follows: Given a linear PDE with boundary conditions and an error parameter , output a quantum state that is -close to one whose amplitudes are proportional to the solution of the PDE at a set of grid points in the domain of the PDE. We focus on elliptic PDEs, and we assume a technical condition that we call global strict diagonal dominance (defined in (2.8)). Our first algorithm is based on a quantum version of the FDM approach: we use a finite-difference approximation to produce a system of linear equations and then solve that system using the QLSA. We analyze our FDM algorithm as applied to Poisson's equation (which automatically satisfies global strict diagonal dominance) under periodic, Dirichlet, and Neumann boundary conditions. Whereas previous FDM approaches [7,12] considered fixed orders of truncation, we adapt the order of truncation depending on , inspired by the classical adaptive FDM [3]. As the order increases, the eigenvalues of the FDM matrix approach the eigenvalues of the continuous Laplacian, allowing for more precise approximations. The main algorithm we present uses the quantum Fourier transform (QFT) and takes advantage of the high-precision LCU-based QLSA [8]. We first consider periodic boundary conditions, but by restricting to appropriate subspaces, this approach can also be applied to homogeneous Dirichlet and Neumann boundary conditions. We state our result in Theorem 1, which (informally) says that this quantum adaptive FDM approach produces a quantum state approximating the solution of Poisson's equation with complexity d 6.5 poly(log d, log(1/ )).
We also propose a quantum algorithm for more general second-order elliptic PDEs under periodic or non-periodic Dirichlet boundary conditions. This algorithm is based on quantum spectral methods [9]. The spectral method globally approximates the solution of a PDE by a truncated Fourier or Chebyshev series (which converges exponentially for smooth functions) with undetermined coefficients, and then finds the coefficients by solving a linear system. This system is exponentially large in d, so solving it is infeasible for classical algorithms but feasible in a quantum context. To be able to apply the QLSA efficiently, we show how to make the system sparse using variants of the quantum Fourier transform. Our bound on the condition number of the linear system uses global strict diagonal dominance, and introduces a factor in the complexity that measures the extent to which this condition holds. We state our result in Theorem 2, which (informally) gives a complexity of d 2 poly(log(1/ )) for producing a quantum state approximating the solution of general second-order elliptic PDEs with Dirichlet boundary conditions.
Both of these approaches have complexity poly(d, log(1/ )), providing optimal dependence on and an exponential improvement over classical methods as a function of the spatial dimension d. Bounding the complexities of these algorithms requires analyzing how d and affect the condition numbers of the relevant linear systems (finite difference matrices and matrices relating the spectral coefficients) and accounting for errors in the approximate solution provided by the QLSA. Furthermore, the complexities of both approaches scale logarithmically with high-order derivatives of the solution and the inhomogeneity. The detailed complexity dependence is presented in Theorem 1 and Theorem 2, and is further discussed in Section 5. Table 1 compares the performance of our approaches to other classical and quantum algorithms for PDEs. Compared to classical algorithms, quantum algorithms improve the dependence on spatial dimension from exponential to polynomial (with the significant caveat that they produce a different representation of the solution). Compared to previous quantum FDM/FEM/FVM algorithms [7,12,14,22], the quantum adaptive FDM and quantum spectral method improve the error dependence from poly(1/ ) to poly(log(1/ )). Our approaches achieve the best known dependence on the parameter for the Poisson equation with homogeneous boundary conditions. Furthermore, our quantum spectral method approach not only achieves the best known dependence on d and for elliptic PDEs with inhomogeneous Dirichlet boundary conditions, but also improves the dependence on d for the Poisson equation with inhomogeneous Dirichlet boundary conditions, as compared to previous quantum algorithms.
The remainder of the paper is structured as follows. Section 2 introduces technical details about linear PDEs and formally states the problem we solve. Section 3 covers our FDM algorithm for Poisson's equation. Section 4 details the spectral algorithm for elliptic PDEs. Finally, Section 5 concludes with a brief discussion of the results, their possible applications, and some open problems.

Linear PDEs
In this paper, we focus on systems of linear PDEs. Such equations can be written in the form and the inhomogeneity f (x) ∈ C are scalar functions, and L is a linear differential operator acting on u(x). In general, L can be written in a linear combination of u(x) and its derivatives. A linear differential operator L of order h has the form 3) The problem reduces to a system of linear ordinary differential equations (ODEs) when For example, systems of first-order linear PDEs can be written in the form A well-known second-order linear PDEs is the Poisson equation A linear PDE of order h is called elliptic if its differential operator (2.2) satisfies for all nonzero ξ j = ξ j 1 1 . . . ξ j d d with ξ 1 , . . . , ξ d ∈ R m and all x. Note that ellipticity only depends on the highest-order terms. When h = 2, the linear PDE (2.5) is called a secondorder elliptic PDE if and only if A j 1 j 2 (x) is positive-definite or negative-definite for any x. In particular, the Poisson equation (2.6) is a second-order elliptic PDE.
We consider a class of elliptic PDEs that also satisfy the condition for all x. We call this condition global strict diagonal dominance, since it is a strengthening of the standard (strict) diagonal dominance condition Observe that (2.8) holds for the Poisson equation (2.6) with C = 1.
In this paper, we focus on the following boundary value problem: In the quantum PDE problem, we are given a system of second-order elliptic equations satisfying the global strict diagonal dominance condition (2.8), where the variable and the linear coefficients A j ∈ C. We are also given boundary conditions We assume there exists a weak solutionû(x) ∈ C for the boundary value problem (see Reference [13,Section 6.1.2]). Given oracles that compute the coefficients A j , and that prepare normalized states |γ(x) and |f (x) whose amplitudes are proportional to γ(x) and f (x) on a set of interpolation nodes x, the goal is to output a quantum state |u(x) whose amplitudes are proportional to u(x) on a set of interpolation nodes x.

Finite difference method
We now describe our first approach to quantum algorithms for linear PDEs, based on the finite difference method (FDM). Using this approach, we show the following.

Theorem 1.
There exists a quantum algorithm that outputs a state -close to |u that runs in timeÕ and makesÕ queries to the oracle for f .
To show this, we first construct a linear system corresponding to the finite difference approximation of Poisson's equation with periodic boundary conditions and bound the error of this high-order FDM in Section 3.1 (Lemma 1). Then we bound the condition number of this system in Section 3.2 (Lemma 2 and Lemma 3) and bound the error of approximation in Section 3.3 (Lemma 4). We use these results to give an efficient quantum algorithm in Section 3.4, establishing Theorem 1. We conclude by discussing how to use the method of images to apply this algorithm for Neumann and Dirichlet boundary conditions in Section 3.5.
The FDM approximates the derivative of a function f at a point x in terms of the values of f on a finite set of points near x. Generally there are no restrictions on where these points are located relative to x, but they are typically taken to be uniformly spaced points with respect to a certain coordinate. This corresponds to discretizing [−1, 1] d (or [0, 2π) d ) to a d-dimensional rectangular lattice (where we use periodic boundary conditions).

Linear system
To approximate the second derivatives appearing in Poisson's equation, we apply the central finite difference formula of order 2k. Taking x j = jh for a lattice with spacing h, this formula gives the approximation where the coefficients are [18,21] r j :=

(3.4)
We leave the dependence on k implicit in this notation. The following lemma characterizes the error of this formula.

Lemma 1 ([18, Theorem 7])
. Let k ≥ 1 and suppose f (x) ∈ C 2k+1 for x ∈ R. Define the coefficients r j as in (3.4). Then Since we assume periodic boundary conditions and apply the same FDM formula at each lattice site, the matrices we consider are circulant. Define the 2n × 2n matrix S to have entries S i,j = δ i,j+1 mod 2n . If we represent the solution u(x) as a vector u = 2n j=1 u(πj/n) e j , then we can approximate Poisson's equation using a central difference formula as where f = 2n j=1 f (πj/n) e j . The solution u corresponds exactly with the quantum state we want to produce, so we do not have to perform any post-processing such as in Reference [12] and other quantum differential equation algorithms. The matrix in this linear system is just the finite difference matrix, so it suffices to bound its condition number and approximation error (whereas previous quantum algorithms involved more complicated linear systems).

Condition number
The following lemma characterizes the condition number of a circulant Laplacian on 2n points.
Proof. We first upper bound L using Gershgorin's circle theorem [17] (a similar argument appears in Reference [18]). Note that The radii of the Gershgorin discs are The discs are centered at r 0 , and so L ≤ 4π 2 3 . To lower bound L −1 we lower bound the (absolute value of the) smallest non-zero eigenvalue of L (since by construction the all-ones vector is a zero eigenvector). Let ω := exp(πi/n). Since L is circulant, its eigenvalues are where the c j ∈ [0, lj] arise from the Taylor remainder theorem. Using (3.8), we have We now compute the sum Therefore, we have Finally, we see that In d dimensions, a similar analysis holds.
Proof. By the triangle inequality for spectral norms, L ≤ d L . Since L has zero-sum rows by construction, the all-ones vector lies in its kernel, and thus the smallest non-zero eigenvalue of L is the same as that of L . Therefore we have which is O(dn 2 ) provided k < (6/π 2 ) 1/3 n 2/3 .

Error analysis
There are two types of error relevant to our analysis: the FDM error and the QLSA error. We assume that we are able to perfectly generate states proportional to f . The FDM errors arise from the remainder terms in the finite difference formulas and from inexact approximations of the eigenvalues. We introduce several states for the purpose of error analysis. Let |u be the quantum state that is proportional to u = j∈Z d 2n u(πj/n) d i=1 e j i for the exact solution of the differential equation. Let |ū be the state output by a QLSA that exactly solves the linear system. Let |ũ be the state output by a QLSA with error. Then the total error of approximating |u by |ũ is bounded by and without loss of generality we can take FDM and QLSA to be of the same order of magnitude.

Lemma 4. Let u( x) be the exact solution of
where is a (2n) d dimensional vector whose entries are O(1). This implies and therefore By Lemma 2 we have λ 1 = Θ(1/n 2 ), and since h = Θ(1/n), we have as claimed.

FDM algorithm
To apply QLSAs, we must consider the complexity of simulating Hamiltonians that correspond to Laplacian FDM operators. For periodic boundary conditions, the Laplacians are circulant, so they can be diagonalized by the QFT F (or a tensor product of QFTs for the multi-dimensional Laplacian L ), i.e., D = F † LF is diagonal. In this case the simplest way to simulate exp(iLt) is to perform the inverse QFT, apply controlled phase rotations to implement exp(iDt), and perform the QFT. Reference [27] shows how to exactly implement arbitrary diagonal unitaries on m qubits using O(2 m ) gates. Since we consider Laplacians on n lattice sites, simulating exp(iLt) takes O(n) gates with the dominant contribution coming from the phase rotations (alternatively, the methods of Reference [36] or Reference [4] could also be used). Using this Hamiltonian simulation algorithm in a QLSA for the FDM linear system gives us the following theorem. We restate Theorem 1 as follows.

Theorem 1.
There exists a quantum algorithm that outputs a state -close to |u that runs in timeÕ and makesÕ queries to the oracle for f .
Proof. We use the Fourier series based QLSA from Reference [8]. By Theorem 3 of that work, the QLSA makes O(κ log(κ/ QLSA )) uses of a Hamiltonian simulation algorithm and uses of the oracle for the inhomogeneity. For Hamiltonian simulation we use d parallel QFTs and phase rotations as described in Reference [27], for a total of O(dnκ log(κ/ QLSA )) gates. The condition number for the d-dimensional Laplacian scales as κ = O(dn 2 ). We take FDM and QLSA to be of the same order and just write . Then the QLSA has time complexity O(d 2 n 3 log(dn 2 / )) and query complexity O(dn 2 log(dn 2 / )). The adjustable parameters are the number of lattice sites n and the order 2k of the finite difference formula. To keep the error below the target error of we require (3.36) or equivalently, Now we focus on the choice of adjustable n and k relying on . This procedure is inspired by the classical adaptive FDM [3], so we call it the adaptive FDM approach. We must have 2k − 1 > d/2 for the left-hand side of (3.37) to be positive for large n. Indeed, we find the best performance by taking k as large as possible subject to the assumption of Lemma 2, i.e., k = cn 2/3 where c := (6/π 2 ) 1/3 . For this choice of k and for n sufficiently large, (3.37) is equivalent to To satisfy the condition 2cn 2/3 − 1 > d/2, we must have n = Ω(d 3/2 ). Combining this observation with (3.38), we choose The QLSA then has the stated time complexitỹ and makes queries to the oracle for f .

Boundary conditions via the method of images
We can apply the method of images to deal with homogeneous Neumann and Dirichlet boundary conditions using the algorithm for periodic boundary conditions described above.
In the method of images, the domain [−1, 1] is extended to include all of R, and the boundary conditions are related to symmetries of the solutions. For a pair of Dirichlet boundary conditions there are two symmetries: the solutions are anti-symmetric about Continuity and anti-symmetry about −1 and 1 imply f (−1) = f (1) = 0, and furthermore that f (x) = 0 for all odd x ∈ Z and that f (x + 4) = f (x) for all x ∈ R. For Neumann boundary conditions, the solutions are instead symmetric about −1 and 1, which also We would like to combine the method of images with the FDM to arrive at finite difference formulas for this special case. In both cases, the method of images implies that the solutions are periodic, so without loss of generality we can consider a lattice on [0, 2π) instead of a lattice on R. It is useful to think of this lattice in terms of the cycle graph on 2n vertices, i.e., (V, E) = (Z 2n , {(i, i + 1) | i ∈ Z 2n }), which means that the vectors encoding the solution u(x) will lie in R 2n . Let each vector e j correspond to the vertex j. Then we divide R 2n into a symmetric and an anti-symmetric subspace, namely span{e j + e 2n+1−j } n j=1 and span{e j − e 2n+1−j } n j=1 , respectively. Vectors lying in the symmetric subspace correspond to solutions that are symmetric about 0 and π, so they obey Neumann boundary conditions at 0 and π; similarly, vectors in the anti-symmetric space correspond to solutions obeying Dirichlet boundary conditions at 0 and π.
Restricting to a subspace of vectors reduces the size of the FDM vectors and matrices we consider, and the symmetry of that subspace indicates how to adjust the coefficients.
If the FDM linear system is L u = f then L has entries where + (−) is chosen for Neumann (Dirichlet) boundary conditions and due to the truncation order k, r j = 0 for any j > k. This is similar to how Laplacian coefficients are modified when imposing boundary conditions in discrete variable representations [11].
For the purpose of solving the new linear systems using quantum algorithms, we still treat these cases as obeying periodic boundary conditions. We assume access to an oracle that produces states |f proportional to the inhomogeneity f (x). Then we apply the QLSA for periodic boundary conditions using |f |± to encode the inhomogeneity, which will output solutions of the form |u |± . Here the ancillary state is chosen to be |+ (|− ) for Neumann (Dirichlet) boundary conditions.
Typically, the (second-order) graph Laplacian for the path graph with Dirichlet boundary conditions has diagonal entries that are all equal to 2; however, using the above specification for the entries of L leads to the (1, 1) and (n, n) entries being 3 while the rest of the diagonal entries are 2.
To reproduce this case, we consider an alternative subspace restriction used in Reference [33] to diagonalize the Dirichlet graph Laplacian. In this case it is easiest to consider the lattice of a cycle graph on 2n + 2 vertices, where the vertices 0 and n + 1 are selected as boundary points where the field takes the value 0. The relevant antisymmetric subspace is now span({e j − e 2n+2−j } n j=1 ) (which has no support on e 0 and e n+1 ). If we again write the linear system as L u = f , then the Laplacian has entries We again assume access to an oracle producing states proportional to f (x); however, we assume that this oracle operates in a Hilbert space with one additional dimension compared to the previous approaches (i.e., whereas previously we considered implementing U , here we consider implementing U 0 0 T 1 ). With this oracle we again prepare the state |f |− and solve Poisson's equation for periodic boundary conditions to output a state |u |− (where |u lies in an (n + 1)-dimensional Hilbert space but has no support on the (n + 1)st basis state).
algorithm based on the pseudo-spectral method [15,28,34] for second-order elliptic equations with global strict diagonal dominance, under various boundary conditions. Using this approach, we show the following.

Theorem 2. Consider an instance of the quantum PDE problem as defined in Problem 1
with Dirichlet boundary conditions (4.28). Then there exists a quantum algorithm that produces a state in the form of (4.29) whose amplitudes are proportional to u(x) on a set of interpolation nodes x (with respect to the uniform grid nodes for periodic boundary conditions or the Chebyshev-Gauss-Lobatto quadrature nodes for non-periodic boundary conditions, as defined in in (4 queries to oracles as defined in Section 4.4. Here The gate complexity is larger than the query complexity by a factor of poly(log(d A Σ / )).
After introducing the method, we discuss the complexity of the quantum shifted Fourier transform (Lemma 5) and the quantum cosine transform (Lemma 6) in Section 4.1. These transforms are used as subroutines in our algorithm. Then we construct a linear system whose solution encodes the solution of the PDE in Section 4.2 (with a simple illustrative example presented in Appendix A), analyze its condition number in Section 4.3 (Lemma 10, established using Lemma 7, Lemma 8, and Lemma 9), and consider the complexity of state preparation in Section 4.4 (Lemma 11). Finally, we prove our main result (Theorem 2) in Section 4.5.
In the spectral approach, we approximate the exact solutionû(x) by a linear combination of basis functions We choose different basis functions for the case of periodic boundary conditions and for the more general case of non-periodic boundary conditions. When the boundary conditions are periodic, the algorithm implementation is more straightforward, and in some cases (e.g., for the Poisson equation), can be faster. Specifically, for any k j ∈ [n + 1] 0 and x j ∈ [−1, 1], we take φ k j (x j ) = e i(k j − n/2 )πx j , periodic conditions, T k j (x j ) := cos(k j arccos x j ), non-periodic conditions. (4.6) Here T k is the degree-k Chebyshev polynomial of the first kind. The coefficients c k are determined by demanding that u(x) satisfies the ODE and boundary conditions at a set of interpolation nodes {χ l = (χ l 1 , . . . , χ l d )} l ∞≤n with l j ∈ [n + 1] 0 , where non-periodic conditions. (4.7) Here { 2l n+1 − 1 : l ∈ [n + 1] 0 } are called the uniform grid nodes, and {cos πl n : l ∈ [n + 1] 0 } are called the Chebyshev-Gauss-Lobatto quadrature nodes.
We require the numerical solution u(x) to satisfy We would like to be able to increase the accuracy of the approximation by increasing n, so 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 = poly(1/ ). Furthermore, if the solution is in C ∞ , the spectral method approximates the solution to within using only n = poly(log(1/ )) [28]. Since we require k j ∈ [n + 1] 0 for all j ∈ [d], we have (n + 1) d terms in total. Consequently, a classical pseudo-spectral method solves multi-dimensional PDEs with complexity poly(log d (1/ )). Such classical spectral methods rapidly become infeasible since the number of coefficients (n + 1) d grows exponentially with d.
Here we develop a quantum algorithm for multi-dimensional PDEs. The algorithm applies techniques from the quantum spectral method for ODEs [9]. However, in the case of PDEs, the linear system to be solved is non-sparse. We address this difficulty using a quantum transform that restores sparsity.

Quantum shifted Fourier transform and quantum cosine transform
The well-known quantum Fourier transform (QFT) can be regarded as an analogue of the discrete Fourier transform (DFT) acting on the amplitudes of a quantum state. The QFT maps the (n + 1)-dimensional quantum state In other words, the QFT is the unitary transform  Here we also consider the quantum shifted Fourier transform (QSFT), an analogue of the classical shifted discrete Fourier transform, which maps v ∈ C n+1 tov ∈ C n+1 witĥ v l = 1 √ n + 1 We define the multi-dimensional QSFT by the tensor product, namely It is well known that F n can be implemented with gate complexity O(log n log log n), and it is straightforward to implement R n and S n with gate complexity O(log n). Thus the total complexity is O(log n log log n).
where v k ∈ C with k = (k 1 , . . . , k d ), and each k j ∈ [n] 0 for j ∈ [d]. The unitary matrix F s n can be written as the tensor product Performing the multi-dimensional QSFT is equivalent to performing the one-dimensional QSFT on each register. Thus, the gate complexity of performing F s n is O(d log n log log n).
Another efficient quantum transformation is the quantum cosine transform (QCT) [19,26]. The QCT can be regarded as an analogue of the discrete cosine transform (DCT). where In other words, the QCT is the orthogonal transform Again we define the multi-dimensional QCT by the tensor product, namely where k = (k 1 , . . . , k d ) and l = (l 1 , . . . , l d ) are d-dimensional vectors with k j , l j ∈ [n + 1] 0 . The classical DCT on (n+1)-dimensional vectors takes Θ(n log n) gates, while the QCT on (n + 1)-dimensional quantum states can be implemented with complexity poly(log n). According to Theorem 1 of Reference [19], the gate complexity of performing C n is O(log 2 n). We observe that this can be improved as follows.

Lemma 6. The quantum cosine transform C n defined by (4.22) can be performed with gate complexity O(log n log log n). More generally, the multi-dimensional QCT C n defined by (4.23) can be performed with gate complexity O(d log n log log n).
Proof. According to the quantum circuit in Figure 2 of Reference [19], C n can be decomposed into a QFT F n+1 , a permutation and additional operations with O(1) cost. The QFT F n+1 has gate complexity O(log n log log n). We then consider an alternative way to implement P n that improves over the approach in [24]. The permutation P n can be decomposed as where F n is the Fourier transform (4.11) and T n = n k=0 e − 2πik n+1 |k k| is diagonal. The gate complexities of performing F n and T n are O(log n log log n) and O(log n), respectively. It follows that C n can be implemented with circuit complexity O(log n log log n).
The matrix C n can be written as the tensor product As in Lemma 5, performing the multi-dimensional QCT is equivalent to performing a QCT on each register. Thus, the gate complexity of performing C n is O(d log n log log n).

Linear system
In this section we introduce the quantum PDE solver for the problem (2.1). We construct a linear system that encodes the solution of (2.1) according to the pseudo-spectral method introduced above, using the QSFT/QCT introduced in Section 4.1 to ensure sparsity. We consider a linear PDE problem (Problem 1) with periodic boundary conditions or non-periodic Dirichlet boundary conditions According to the elliptic regularity theorem (Theorem 6 in Section 6.3 of Reference [13]), there exists a unique solutionû(x) in C ∞ for Problem 1.
We now show how to apply the Fourier and Chebyshev pseudo-spectral methods to this problem. Our goal is to obtain the quantum state where φ k (χ l ) is defined by (4.5) using (4.6) for the appropriate boundary conditions (periodic or non-periodic). This state corresponds to a truncated Fourier/Chebyshev approximation and is -close to the exact solutionû(χ l ) with n = poly(log(1/ )) [28]. Note that this state encodes the values of the solution at the interpolation nodes (4.7) appropriate to the boundary conditions (the uniform grid nodes in the Fourier approach, for periodic boundary conditions, and the Chebyshev-Gauss-Lobatto quadrature nodes in the Chebyshev approach, for non-periodic boundary conditions). Instead of developing our algorithm for the standard basis, we aim to produce a state |c ∝ k ∞≤n c k |k 1 . . . |k d (4.30) that is the inverse QSFT/QCT of |u . We then apply the QSFT/QCT to transform back into the interpolation node basis. The truncated spectral series of the inhomogeneity f (x) and the boundary conditions γ(x) can be expressed as and respectively. We define quantum states |f and |γ by interpolating the nodes {χ l } defined by (4.7) as |f ∝ and respectively. These are the states that we assume we can produce using oracles. We perform the multi-dimensional inverse QSFT/QCT to obtain the states Having defined these states, we now detail the construction of the linear system. At a high level, we construct two linear systems: one system A x = f (where x corresponds to (4.30)) describes the differential equation, and another system B x = g describes the boundary conditions. We combine these into a linear system with the form (4.37) Even though we do not impose the two linear systems separately, we show that there exists a unique solution of (4.37) (which is therefore the solution of the simultaneous equations A x = f and B x = g), since we show that L has full rank, and indeed we upper bound its condition number in Section 4.3. Part of this linear system will correspond to just the differential equation while another part will come from imposing the boundary conditions on ∂D can be expressed as conditions on each boundary: (4.40)

Linear system from the differential equation
To evaluate the matrix corresponding to the differential operator from (4.38), it is convenient to define coefficients c (j) k and k ∞ ≤ n such that for some fixed j ∈ N d (as we explain below, such a decomposition exists for the choices of basis functions in (4.6)). Using this expression, we obtain the following linear equations for c (j) k : To determine the transformation between c (j) k and c k , we can make use of the differential properties of Fourier and Chebyshev series, namely and where D (j) n can be expressed as the tensor product with j = (j 1 , . . . , j d ). The matrix D n for the Fourier basis functions in (4.6) can be written as the (n + 1) × (n + 1) diagonal matrix with entries As detailed in Appendix A of Reference [9], the matrix D n for the Chebyshev polynomials in (4.6) can be expressed as the (n + 1) × (n + 1) upper triangular matrix with nonzero entries Substituting (4.46) into (4.42), with D n defined by (4.47) in the periodic case or (4.48) in the non-periodic case, and performing the multi-dimensional inverse QSFT/QCT (for a reason that will be explained in the next section), we obtain the following linear equations for c r : Notice that the matrices (4.47) and (4.48) are not full rank. More specifically, there exists at least one zero row in the matrix of (4.50) when using either (4.47) (k = n/2 ) or (4.48) (k = n). To obtain an invertible linear system, we next introduce the boundary conditions.

Adding the linear system from the boundary conditions
When we use the form (4.29) of u(x) to write linear equations describing the boundary conditions (4.40), we obtain a non-sparse linear system. Thus, for each x ∈ ∂D j in (4.40), we perform the (d − 1)-dimensional inverse QSFT/QCT on the d − 1 registers except the jth register to obtain the linear equations where the values of k j indicate that we place these constraints in the last two rows with respect to the jth coordinate. We combine these equations with (4.50) to obtain the linear system n defined below. In other words, D n for each j that has exactly one entry equal to 2 and all other entries 0, whereas D Alternatively, for the Chebyshev case in (4.6) used for non-periodic boundary conditions, D n comes from (4.48), and the nonzero entries of G n are The system (4.52) has the form of (4.37). For instance, the matrix in n ⊗ I ⊗d−1 + I ⊗ D (2) n ⊗ I ⊗d−2 + · · · + I ⊗d−1 ⊗ D (2) n . (4.58) For periodic boundary conditions, using (4.45), (4.47), and (4.55), the second-order differential matrix D (2) n has nonzero entries [D (2) n ] k,k = −((k − n/2 )π) 2 , k ∈ [n + 1] 0 \{ n/2 }, [D (2) n ] n/2 ,k = 1, k ∈ [n + 1] 0 .
n has nonzero entries [D (2) (4.60) To explicitly illustrate this linear system, we present a simple example in Appendix A. We discuss the invertible linear system (4.52) and upper bound its condition number in the following section.

Condition number
We now analyze the condition number of the linear system. We begin with two lemmas bounding the singular values of the matrices (4.59) and (4.60) that appear in the linear system.

Lemma 7.
Consider the case of periodic boundary conditions. Then for n ≥ 4, the largest and smallest singular values of D (2) n defined in (4.59) satisfy   4 (4.74) as claimed.
We now consider the condition number of the linear system for general elliptic PDEs.

Lemma 10. Consider an instance of the quantum PDE problem as defined in Problem 1
with Dirichlet boundary conditions (4.28). Then for n ≥ 4, the condition number of L in the linear system (4.37) satisfies |A j,j |, and C > 0 is defined in (2.8).
Recall that C quantifies the extent to which the global strict diagonal dominance condition holds.
Proof. According to (4.52), the matrix in (4.37) is (4.76) We upper bound the spectral norm of the matrix L by For the matrix D Next we lower bound Lξ for any ξ = 1. It is non-trivial to directly compute the singular values of a sum of non-normal matrices. Instead, we write L as a sum of terms L 1 and L 2 , where L 1 is a tensor sum similar to (4.58) that can be bounded by Lemma 9, and L 2 is a sum of tensor products that are easily bounded. Specifically, we have The ellipticity condition (2.7), ∀ξ j 1 =h A j (x)ξ j = 0, can only hold if the A j,j for j ∈ [d] are either all positive or all negative; we consider A j,j > 0 without loss of generality, so (4.81) Also, the global strict diagonal dominance condition (2.8) simplifies to where 0 < C ≤ 1. for each j = (j 1 , . . . , j d ) that has exactly two entries equal to 1 and all other entries 0. Specifically, consider j r 1 = j r 2 = 1 for r 1 , r 2 ∈ [d], r 1 = r 2 , and j r = 0 for r ∈ [d]\{r 1 , r 2 }. We denote We first upper bound D n and L (j) share the same singular vectors. For k ∈ [n + 1] 0 , we let v k and λ k denote the right singular vectors and corresponding singular values of D n , respectively. Then the right singular vectors of D (j) n v ≤ 1 2 L (j) v by the inequality of arithmetic and geometric means (also known as the AM-GM inequality). Since this holds for any vector v, we have Next we upper bound D 2 n by D Notice that w n/2 = 0 and w k = w k for k ∈ [n + 1] 0 \{ n/2 } for periodic conditions, and w n−1 = w n = 0 and w k = w k for k ∈ [n + 1] 0 \{n − 1, n} for non-periodic conditions. Thus, for any vector v, Therefore, We also have We can rewrite I ⊗rs−1 ⊗ D (2) n ⊗ I ⊗d−rs L −1 1 in the form (4.92) The matrices I ⊗r h −1 ⊗ D (2) n ⊗ I ⊗d−r h share the same singular values and singular vectors, so . This implies ). (4.94) Using (4.82), considering each instance of D (j) n in L 2 , we have Since L and L 1 are invertible, L −1 1 L 2 ≤ 1 − C < 1, and by Lemma 9 applied to L −1 1 , we have (4.96) Thus we have as claimed.

State preparation
We now describe a state preparation procedure for the vector f + g in the linear system (4.37). (4.98) Then the normalized quantum state with coefficients defined as in (4.35) and (4.36), can be prepared with gate and query complexity O(qd 2 log n log log n).
Proof. Starting from the initial state |0 |0 , we first perform a unitary transformation U satisfying on the first register to obtain This can be done in time O(2d + 1) by standard techniques [31]. Then we apply O x and O f to obtain Finally, observe that if we measure the first register in a basis containing the uniform superposition |0 + |1 + · · · + |2d (say, the Fourier basis) and obtain the outcome corresponding to the uniform superposition, we produce the state Since this outcome occurs with probability 1/q 2 , we can prepare this state with probability close to 1 using O(q) steps of amplitude amplification. According to Lemma 5 and Lemma 6, the d-dimensional (inverse) QSFT or QCT can be performed with gate complexity O(d log n log log n). Thus the total gate and query complexity is O(qd 2 log n log log n).
Alternatively, if it is possible to directly prepare the quantum state |B , then we may be able to avoid the factor of q in the complexity of the overall algorithm.

Main result
Having analyzed the condition number and state preparation procedure for our approach, we are now ready to establish the main result. Theorem 2 as follows. (4.28). Then there exists a quantum algorithm that produces a state in the form of (4.29) whose amplitudes are proportional to u(x) on a set of interpolation nodes x (with respect to the uniform grid nodes for periodic boundary conditions or the Chebyshev-Gauss-Lobatto quadrature nodes for non-periodic boundary conditions, as defined in in (4.7)), where u(x)/ u(x) is -close toû(x)/ û(x) in l 2 norm for all nodes x, succeeding with probability Ω(1), with a flag indicating success, using d A Σ C A * + qd 2 poly(log(g /g )) (4.1)

Theorem 2. Consider an instance of the quantum PDE problem as defined in Problem 1 with Dirichlet boundary conditions
queries to oracles as defined in Section 4.4. Here (2.8), and The gate complexity is larger than the query complexity by a factor of poly(log(d A Σ / )).
Proof. We analyze the complexity of the algorithm presented in Section 4.2.
First we choose n := log(Ω) log(log(Ω)) , (4.105) where Ω = g (1 + ) g . (4.106) By Eq. (1.8.28) of Reference [34], this choice guarantees so we can choose n to ensure that the normalized output state is -close toû(x)/ û(x) . As described in Section 4.2, the algorithm uses the high-precision QLSA from Reference [8] and the multi-dimensional QSFT/QCT (and its inverse). According to Lemma 5 and Lemma 6, the d-dimensional (inverse) QSFT or QCT can be performed with gate complexity O(d log n log log n). According to Lemma 11, the query and gate complexity for state preparation is O(qd 2 log n log log n).
For the linear system L x = f + g in (4.37), the matrix L is an (n + 1) d × (n + 1) d matrix with (n + 1) or (n + 1)d nonzero entries in any row or column for periodic or non-periodic conditions, respectively. According to Lemma 10, the condition number of L is upper bounded by A Σ C A * (2n) 4 . Consequently, by Theorem 5 of Reference [8], the QLSA produces a state proportional to x with O( d A Σ C A * (2n) 5 ) queries to the oracles, and its gate complexity is larger by a factor of poly(log(d A Σ n)). Using the value of n specified in (4.105), the overall query complexity of our algorithm is d A Σ C A * + qd 2 poly(log(g /g )), (4.109) and the gate complexity is which is larger by a factor of poly(log(d A Σ / )), as claimed.
Note that we can establish a more efficient algorithm in the special case of the Poisson equation with homogenous boundary conditions. In this case, A Σ = A * = d and C = 1. Under homogenous boundary conditions, the complexity of state preparation can be reduced to d poly(log(g /g )), since we can remove 2d applications of the QSFT or QCT for preparing a state depending on the boundary conditions, and since γ = 0 there is no need to postselect on the uniform superposition to incorporate the boundary conditions. In summary, the query complexity of the Poisson equation with homogenous boundary conditions is d poly(log(g /g )); (4.111) again the gate complexity is larger by a factor of poly(log(d A Σ / )).

Discussion and open problems
We have presented high-precision quantum algorithms for d-dimensional PDEs using the FDM and spectral methods. These algorithms use high-precision QLSAs to solve Poisson's equation and other second-order elliptic equations. Whereas previous algorithms scaled as poly(d, 1/ ), our algorithms scale as poly(d, log(1/ )). This work raises several natural open problems. First, for the quantum adaptive FDM, we only deal with Poisson's equation with homogeneous boundary conditions. Can we apply the adaptive FDM to other linear equations or to inhomogeneous boundary conditions? The quantum spectral algorithm applies to second-order elliptic PDEs with Dirichlet boundary conditions. Can we generalize it to other linear PDEs with Neumann or mixed boundary conditions? Also, can we develop algorithms for space-and timedependent PDEs? These cases are more challenging since the quantum Fourier transform cannot be directly applied to ensure sparsity. Finally, can we improve the dependence on d?
Second, the complexity scales logarithmically with high-order derivatives (of the inhomogeneity or solution) for both the adaptive FDM and the spectral method. In particular, Theorem 1 shows that the complexity of the quantum adaptive FDM scales logarithmically with d 2k+1 u dx 2k+1 , and Theorem 2 shows that the complexity of the quantum spectral method is poly(log g ), where g upper bounds û (n+1) (x) (see (4.2)). Such a logarithmic dependence on high-order derivatives of the solution is typical for classical algorithms, including the classical adaptive FDM (see for example Theorem 7 of Reference [18]) and spectral methods (see for example Eq. (1.8.28) of Reference [34]), both of which have the same logarithmic dependence on d 2k+1 u dx 2k+1 and g . This logarithmic dependence means that the algorithm is efficient even when faced with a highly oscillatory solution with an exponentially large derivative.
However, the query complexity of time-dependent Hamiltonian simulation only depends on the first-order derivatives of the Hamiltonian [5,23]. Can we develop quantum algorithms for PDEs with query complexity independent of high-order derivatives, and henceforth develop an unexpected advantage of quantum algorithms for PDEs?
Third, can we develop quantum algorithms for other types of PDEs, such as stochastic or nonlinear PDEs?
Fourth, can we use quantum algorithms for PDEs as a subroutine of other quantum algorithms? For example, some PDE algorithms have state preparation steps that require inverting finite difference matrices (such as Reference [12] using certain oracles for the initial conditions); are there other scenarios in which state preparation can be done using the solution of another system of PDEs? Can quantum algorithms for PDEs be applied to other algorithmic tasks, such as optimization?
Finally, how should these algorithms be applied? While PDEs have broad applications, much more work remains to understand the extent to which quantum algorithms can be of practical value. Answering this question will require careful consideration of various technical aspects of the algorithms. In particular: What measurements give useful information about the solutions, and how can those measurements be efficiently implemented? How should the oracles encoding the equations and boundary conditions be implemented in practice? And with these aspects taken into account, what are the resource requirements for quantum computers to solve classically intractable problems related to PDEs?

A An example for solving Poisson's equation
In this appendix, we present an example of solving Poisson's equation in two dimensions using our algorithm. The Poisson equation is We consider two kinds of boundary value problems, as follows.
• Periodic boundary conditions: • Non-periodic boundary conditions: We first present the quantum Fourier spectral method to solve (A.1) with the periodic conditions (A.2). In particular, we choose n = 2 in the specification of the linear system. The truncated Fourier series can be written as We are given an oracle for preparing the state that interpolates the uniform grid nodes (4.7). We first perform a multi-dimensional inverse QSFT on (A.5) to obtain where b k 1 ,k 2 are defined by (4.35). Then we apply the quantum linear system algorithm of Reference [8] to solve the linear system where the solution is According to (4.52),the discretized linear system from (A.1) is where the Fourier difference matrix D n is defined by (4.47) with n = 2, namely so that Therefore, the matrix (A.9) is The rank of this matrix is (n + 1) d − 1 with d = 2, n = 2. We use the boundary condition to complete the linear system: where the additional linear equation comes from u(0, 0) = 2 In some problems, we might be directly given the value of c 1,1 , say, c 1,1 = γ. Then the linear system would be We now present the quantum Chebyshev spectral method to solve (A.1) with nonperiodic conditions (A.3). Similarly, we choose n = 3 in the specification of the linear system. The truncated Chebyshev series of the solution can be written as We are given an oracle for preparing the state that interpolates the Chebyshev-Gauss-Lobatto quadrature nodes specified in (4.7). We first perform the multi-dimensional inverse QCT on (A.5) to obtain (A.6), where f k 1 ,k 2 are defined by (4.35). Then we apply the quantum linear system algorithm of Reference [8] We now use the boundary conditions to complete the linear system. The truncated Chebyshev series of the solution can be written as We are given an oracle for preparing the state by interpolating the Chebyshev-Gauss-Lobatto quadrature nodes specified in (4.7) 3 (A.23) We perform the multi-dimensional inverse QCT on (A.23) to obtain where a k 1 ,k 2 are defined by (4.36). The linear system from the boundary conditions is (A. 25) Adding the two linear systems (A.21) and (A.25) together, we obtain a full-rank linear system D 3 ⊗ I + I ⊗ D In summary, the linear system including the differential equations and the boundary conditions is f0,0 f0,1 f0,2 + gS 0 f0,3 + gN 0 f1,0 f1,1 f1,2 + gS 1 f1,3 + gN 1 f2,0 + gW 0 f2,1 + gW 1 gW 2 + gS 2 gW 3 + gN 2 f3,0 + gE 0 f3,1 + gE 1 gE 2 + gS 3 (A. 28) B Singular values of second-order differential matrices Here we present a detailed proof of the singular value estimation in Lemma 7 and Lemma 8.