# High-precision quantum algorithms for partial differential equations

Andrew M. Childs1,2,3 Jin-Peng Liu1,3,4 Aaron Ostrander1,3,5
1Joint Center for Quantum Information and Computer Science,
2Department of Computer Science,
3Institute for Advanced Computer Studies,
4Applied Mathematics, Statistics, and Scientific Computation, and
5Department of Physics,
University of Maryland

###### Abstract

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 $\operatorname{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
$\operatorname{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.

##  1 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. 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 $\Omega(n^{d})$.

Beyond uniform grids, the sparse grid technique [31] has been applied to
reduce the time and space complexity of outputting a sparse encoding of the
solution to $O(n\log^{d}n)$ [5, 37]. 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 [5].
Reference [5] 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_ [14, 27].
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 $\Omega(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 $\Omega(\log N)$ space—and correspondingly provides more limited
access to the solution—but can potentially be done in only
$\operatorname{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) [15, 1, 7]. For a linear system $A\vec{x}=\vec{b}$, a QLSA
outputs a quantum state proportional to the solution $\vec{x}$. To learn
information about the solution $\vec{x}$, the output of the QLSA must be post-
processed. For example, to output all the entries of an $N$-dimensional vector
$\vec{x}$ given a quantum state $|{x}\rangle$ proportional to it, even a
quantum computer needs time and space $\Omega(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 [9] 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 [21] 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 $\epsilon$. 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 $\epsilon$), since all of these affect the complexity. Their
algorithms have complexity $\operatorname{poly}(d,1/\epsilon)$, compared to
$O((1/\epsilon)^{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 [6, 33]
apply the FDM to solve Poisson’s equation in rectangular volumes under
Dirichlet boundary conditions. Although the circuits they construct have
$\operatorname{poly}(\log(1/\epsilon))$ gates, these circuits have success
probability $\operatorname{poly}(1/\epsilon)$, leading to
$\operatorname{poly}(1/\epsilon)$ time complexity. Additionally, they do not
quantify errors resulting from the finite-difference approximation. Reference
[11] applies the FDM to the problem of outputting states proportional to
solutions of the wave equation, giving complexity
$d^{\frac{5}{2}}\operatorname{poly}(1/\epsilon)$, a polynomial dependence on
$d$ and $1/\epsilon$ (which is $\operatorname{poly}(n)$ for a fixed-order
FDM). The FVM is combined with the reservoir method in Ref. [13] to simulate
hyperbolic equations; although they achieve linear scaling with respect to the
spatial dimension, they use fixed order differences, leading to
$\operatorname{poly}(1/\epsilon)$ scaling. These FDM, FEM, and FVM approaches
can only give a total complexity $\operatorname{poly}(1/\epsilon)$, 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 Ref. [17] 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 [36, 35] and for electronic structure problems, including for quantum
chemistry (see for example [24, 19]). 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 Ref. [21], 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
$\operatorname{poly}(\log(1/\epsilon))$. 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 $\epsilon$, output a quantum state that is $\epsilon$-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.

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 under periodic, Dirichlet, and Neumann boundary
conditions. Whereas previous FDM approaches [6, 11] considered fixed orders of
truncation, we adjust the order of truncation depending on $\epsilon$. 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 [7]. 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 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 [8]. 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.

Both of these approaches have complexity
$\operatorname{poly}(d,\log(1/\epsilon))$, providing optimal dependence on
$\epsilon$ 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 $\epsilon$ 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.

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 [6, 13, 11, 21], the quantum adaptive FDM and quantum spectral
method improve the error dependence from $\operatorname{poly}(1/\epsilon)$ to
$\operatorname{poly}(\log(1/\epsilon))$. Our approaches achieve the best known
dependence on the parameters $d$ and $\epsilon$ for the Poisson equation with
homogeneous boundary conditions. Furthermore, our quantum spectral method
approach not only achieves the best known dependence on $d$ and $\epsilon$ 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.

| Algorithm | Equation | Boundary conditions | Complexity
---|---|---|---|---
Classical | FDM/FEM/FVM | general | general | $\operatorname{poly}((1/\epsilon)^{d})$
Adaptive FDM/FEM [2] | general | general | $\operatorname{poly}(({\color[rgb]{0.00,0.07,1.00}\log(1/\epsilon)})^{d})$
Spectral method [14, 27] | general | general | $\operatorname{poly}(({\color[rgb]{0.00,0.07,1.00}\log(1/\epsilon)})^{d})$
Sparse grid FDM/FEM [5, 37] | general | general | $\operatorname{poly}((1/\epsilon)(\log(1/\epsilon))^{d})$
Sparse grid spectral method [28, 29] | elliptic | general | $\operatorname{poly}({\color[rgb]{0.00,0.07,1.00}\log(1/\epsilon)(\log\log(1/\epsilon))}^{d})$
Quantum |  FEM [21] | Poisson | homogeneous | $\operatorname{poly}(d,1/\epsilon)$
FDM [6] | Poisson | homogeneous Dirichlet | ${\color[rgb]{0.00,0.07,1.00}d}\operatorname{poly}({\color[rgb]{0.00,0.07,1.00}\log d},1/\epsilon)$
FDM [11] | wave | homogeneous | ${\color[rgb]{0.00,0.07,1.00}d^{5/2}}\operatorname{poly}(1/\epsilon)$
FVM [13] | hyperbolic | periodic | ${\color[rgb]{0.00,0.07,1.00}d}\operatorname{poly}(1/\epsilon)$
Adaptive FDM [this paper] | Poisson | periodic, homogeneous | ${\color[rgb]{0.00,0.07,1.00}d^{2}}\operatorname{poly}({\color[rgb]{0.00,0.07,1.00}\log d},{\color[rgb]{0.00,0.07,1.00}\log(1/\epsilon)})$
Spectral method [this paper] | Poisson | homogeneous Dirichlet | ${\color[rgb]{0.00,0.07,1.00}d}\operatorname{poly}({\color[rgb]{0.00,0.07,1.00}\log d},{\color[rgb]{0.00,0.07,1.00}\log(1/\epsilon)})$
Spectral method [this paper] | elliptic | inhomogeneous Dirichlet | ${\color[rgb]{0.00,0.07,1.00}d^{2}}\operatorname{poly}({\color[rgb]{0.00,0.07,1.00}\log(1/\epsilon)})$
Table 1:  Summary of the time complexities of classical and quantum algorithms
for $d$-dimensional PDEs with error tolerance $\epsilon$. Portions of the
complexity in blue represent best known dependence on that parameter.

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 and some open problems.

##  2 Linear PDEs

In this paper, we focus on systems of linear PDEs. Such equations can be
written in the form

$\mathscr{L}(u(\bm{x}))=f(\bm{x}),$ (2.1)

where the variable $\bm{x}=(x_{1},\ldots,x_{d})\in{\mathbb{C}}^{d}$ is a
$d$-dimensional vector, the solution $u(\bm{x})\in{\mathbb{C}}$ and the
inhomogeneity $f(\bm{x})\in{\mathbb{C}}$ are scalar functions, and
$\mathscr{L}$ is a linear differential operator acting on $u(\bm{x})$. In
general, $\mathscr{L}$ can be written in a linear combination of $u(\bm{x})$
and its derivatives. A linear differential operator $\mathscr{L}$ of order $h$
has the form

$\mathscr{L}(u(\bm{x}))=\sum_{\|\bm{j}\|_{1}\leq
h}A_{\bm{j}}(\bm{x})\frac{\partial^{\bm{j}}}{\partial\bm{x}^{\bm{j}}}u(\bm{x}),$
(2.2)

where $\bm{j}=(j_{1},\ldots,j_{d})$ is a $d$-dimensional non-negative vector
with $\|{\bm{j}}\|_{1}=j_{1}+\dots+j_{d}\leq h$,
$A_{j}(\bm{x})\in{\mathbb{C}}$, and

$\frac{\partial^{\bm{j}}}{\partial\bm{x}^{\bm{j}}}u(\bm{x})=\frac{\partial^{j_{1}}}{\partial
x_{1}^{j_{1}}}\cdots\frac{\partial^{j_{d}}}{\partial x_{d}^{j_{d}}}u(\bm{x}).$
(2.3)

The problem reduces to a system of linear ordinary differential equations
(ODEs) when $d=1$. For $d\geq 2$, we call (2.1) a (multi-dimensional) PDE.

For example, systems of first-order linear PDEs can be written in the form

$\sum_{j=1}^{d}A_{j}(\bm{x})\frac{\partial u(\bm{x})}{\partial
x_{j}}+A_{0}(\bm{x})u(\bm{x})=f(\bm{x}),$ (2.4)

where $A_{j}(\bm{x}),A_{0}(\bm{x}),f(\bm{x})\in{\mathbb{C}}$ for
$j\in[{d}]:=\\{1,\ldots,d\\}$. Similarly, systems of second-order linear PDEs
can be expressed in the form

$\sum_{j_{1},j_{2}=1}^{d}A_{j_{1}j_{2}}(\bm{x})\frac{\partial^{2}u(\bm{x})}{\partial
x_{j_{1}}\partial x_{j_{2}}}+\sum_{j=1}^{d}A_{j}(\bm{x})\frac{\partial
u(\bm{x})}{\partial x_{j}}+A_{0}(\bm{x})u(\bm{x})=f(\bm{x}),$ (2.5)

where
$A_{j_{1},j_{2}}(\bm{x}),A_{j}(\bm{x}),A_{0}(\bm{x}),f(\bm{x})\in{\mathbb{C}}$
for $j_{1},j_{2},j\in[{d}]$. A well-known second-order linear PDEs is the
_Poisson equation_

$\Delta u(\bm{x}):=\sum_{j=1}^{d}\frac{\partial^{2}}{\partial
x_{j}^{2}}u(\bm{x})=f(\bm{x}).$ (2.6)

A linear PDE of order $h$ is called _elliptic_ if its differential operator
(2.2) satisfies

$\sum_{\|\bm{j}\|_{1}=h}A_{\bm{j}}(\bm{x})\bm{\xi}^{\bm{j}}\neq 0,$ (2.7)

for all nonzero $\bm{\xi}^{\bm{j}}=\xi_{1}^{j_{1}}\ldots\xi_{d}^{j_{d}}$ with
$\xi_{1},\ldots,\xi_{d}\in{\mathbb{R}}^{m}$. Note that ellipticity only
depends on the highest-order terms. When $h=2$, the linear PDE (2.5) is called
a second-order elliptic PDE if and only if $A_{j_{1}j_{2}}(\bm{x})$ is
positive-definite or negative-definite for any $\bm{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 following
condition:

$C:=1-\sum_{j_{1}=1}^{d}\frac{1}{|A_{j_{1},j_{1}}(\bm{x})|}\sum_{j_{2}\in[{d}]\backslash\\{j_{1}\\}}|A_{j_{1},j_{2}}(\bm{x})|>0.$
(2.8)

We call this condition _global strict diagonal dominance_ , since it is a
strengthening of standard (strict) diagonal dominance. 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:

######  Problem 1.

In the _quantum PDE problem_ , we are given a system of second-order elliptic
equations

$\mathscr{L}(u(\bm{x}))=\sum_{\|\bm{j}\|_{1}=2}A_{\bm{j}}\frac{\partial^{\bm{j}}}{\partial\bm{x}^{\bm{j}}}u(\bm{x})=\sum_{j_{1},j_{2}=1}^{d}A_{j_{1}j_{2}}\frac{\partial^{2}u(\bm{x})}{\partial
x_{j_{1}}\partial x_{j_{2}}}=f(\bm{x})$ (2.9)

satisfying the global strict diagonal dominance condition (2.8), where the
variable $\bm{x}=(x_{1},\ldots,x_{d})\in{\mathscr{D}}=[-1,1]^{d}$ is a
$d$-dimensional vector, the inhomogeneity $f(\bm{x})\in{\mathbb{C}}$ is a
scalar function of $\bm{x}$ satisfying $f(\bm{x})\in C^{\infty}$, and the
linear coefficients $A_{\bm{j}}\in{\mathbb{C}}$. We are also given boundary
conditions $u(\bm{x})=\gamma(\bm{x})\in\partial{\mathscr{D}}$ or
$\frac{\partial u(\bm{x})}{\partial x_{j}}{\big{|}}_{x_{j}=\pm
1}=\gamma(\bm{x})|_{x_{j}=\pm 1}\in\partial{\mathscr{D}}$ where
$\gamma(\bm{x})\in C^{\infty}$. We assume there exists a weak solution
$\hat{u}(\bm{x})\in{\mathbb{C}}$ for the boundary value problem (see Ref. [12,
Section 6.1.2]). Given oracles that compute the coefficients $A_{\bm{j}}$, and
that prepare normalized states $|\gamma(\bm{x})\rangle$ and
$|f(\bm{x})\rangle$ whose amplitudes are proportional to $\gamma(\bm{x})$ and
$f(\bm{x})$ on a set of interpolation nodes $\bm{x}$, the goal is to output a
quantum state $|u(\bm{x})\rangle$ whose amplitudes are proportional to
$u(\bm{x})$ on a set of interpolation nodes $\bm{x}$.

##  3 Finite difference method

We now describe our first approach to quantum algorithms for linear PDEs,
based on the finite difference method (FDM). We first construct the linear
system corresponding to the finite difference approximation of Poisson’s
equation with periodic boundary conditions (Section 3.1). Then we bound its
condition number (Section 3.2) and the error of approximation (Section 3.3),
and we use these results to give an efficient quantum algorithm (Section 3.4).
We conclude by discussing how to use the method of images to apply this
algorithm for Neumann and Dirichlet boundary conditions (Section 3.5).

The FDM approximates the derivative of a function $f$ at a point $\bm{x}$ in
terms of the values of $f$ on a finite set of points near $\bm{x}$. Generally
there are no restrictions on where these points are located relative to
$\bm{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\pi)^{d}$) to a $d$-dimensional rectangular lattice (where we use
periodic boundary conditions).

For a scalar field, in which $u(\bm{x})\in{\mathbb{C}}$, the canonical
elliptic PDE is Poisson’s equation (2.6), which we consider solving on
$[0,2\pi)^{d}$ with periodic boundary conditions. This also implies results
for the domain $\Omega=[-1,1]^{d}$ under Dirichlet ($u(\partial\Omega)=0$) and
Neumann ($\hat{n}\cdot\nabla u(\partial\Omega)=0$ where $\hat{n}$ denotes the
normal direction to $\partial\Omega$, which for domain $\Omega=[-1,1]^{d}$ is
equivalent to $\frac{\partial u}{\partial x_{j}}\big{|}_{x_{j}=\pm 1}=0$ for
$j\in[d]$) boundary conditions.

###  3.1 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

$\displaystyle
f^{\prime\prime}(0)\approx\frac{1}{h^{2}}\sum_{j=-k}^{k}r_{j}f(jh)$ (3.1)

where the coefficients are [17, 20]

$\displaystyle
r_{j}:=\begin{cases}\frac{2(-1)^{j+1}(k!)^{2}}{j^{2}(k-j)!(k+j)!}&j\in[k]\\\
-2\sum_{j=1}^{k}r_{j}&j=0\\\ r_{-j}&j\in-[k].\end{cases}$ (3.2)

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

######  Lemma 1 ([17, Theorem 7]).

Let $k\geq 1$ and suppose $f(x)\in C^{2k+1}$ for $x\in\mathbb{R}$. Define the
coefficients $r_{j}$ as in (3.2). Then

$\displaystyle\frac{\mathrm{d}^{2}u(x_{0})}{\mathrm{d}x^{2}}=\frac{1}{h^{2}}\sum_{j=-k}^{k}r_{j}f(x_{0}+jh)+O\Bigl{(}\Bigl{|}\frac{\mathrm{d}^{2k+1}u}{\mathrm{d}x^{2k+1}}\Bigr{|}\Bigl{(}\frac{eh}{2}\Bigr{)}^{2k-1}\Bigr{)}$
(3.3)

where

$\displaystyle\Bigl{|}\frac{\mathrm{d}^{2k+1}u}{\mathrm{d}x^{2k+1}}\Bigr{|}:=\max_{y\in[x_{0}-kh,x_{0}+kh]}\Bigl{|}\frac{\mathrm{d}^{2k+1}u}{\mathrm{d}x^{2k+1}}(y)\Bigr{|}.$
(3.4)

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\times 2n$ matrix $S$ to have entries $S_{i,j}=\delta_{i,j+1\bmod 2n}$. If
we represent the solution $u(x)$ as a vector $\vec{u}=\sum_{j=1}^{2n}u(\pi
j/n)\vec{e}_{j}$, then we can approximate Poisson’s equation using a central
difference formula as

$\displaystyle\frac{1}{h^{2}}L\vec{u}=\frac{1}{h^{2}}\Bigl{(}r_{0}I+\sum_{j=1}^{k}r_{j}(S^{j}+S^{-j})\Bigr{)}\vec{u}=\vec{f}$
(3.5)

where $\vec{f}=\sum_{j=1}^{2n}f(\pi j/n)\vec{e}_{j}$. The solution $\vec{u}$
corresponds exactly with the quantum state we want to produce, so we do not
have to perform any post-processing such as in Ref. [11] 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).

###  3.2 Condition number

The following lemma characterizes the condition number of a circulant
Laplacian on $2n$ points.

######  Lemma 2.

For $k=o(n^{2/3})$, the matrix $L=r_{0}I+\sum_{j=1}^{k}r_{j}(S^{j}+S^{-j})$
with $r_{j}$ as in (3.2) has condition number $\kappa(L)=\Theta(n^{2})$.

###### Proof.

We first upper bound $\|{L}\|$ using Gershgorin’s circle theorem [16] (a
similar argument appears in Ref. [17]). Note that

$\displaystyle|r_{j}|=\frac{2(k!)^{2}}{j^{2}(k-j)!(k+j)!}\leq\frac{2}{j^{2}}$
(3.6)

since

$\displaystyle\frac{(k!)^{2}}{(k-j)!(k+j)!}=\frac{k(k-1)\cdots(k-j+1)}{(k+j)(k+j-1)\cdots(k+1)}<1.$
(3.7)

The radii of the Gershgorin discs are

$\displaystyle 2\sum_{j=1}^{k}|r_{j}|$ $\displaystyle\leq
2\sum_{j=1}^{k}\frac{2}{j^{2}}\leq\frac{2\pi^{2}}{3}.$ (3.8)

The discs are centered at $r_{0}$, and

$\displaystyle|r_{0}|$ $\displaystyle\leq
2\sum_{j=1}^{k}|r_{j}|\leq\frac{2\pi^{2}}{3},$ (3.9)

so $\|{L}\|\leq\frac{4\pi^{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 $\omega:=\exp(\pi i/n)$. Since $L$ is circulant,
its eigenvalues are

$\displaystyle\lambda_{l}$
$\displaystyle=r_{0}+\sum_{j=1}^{k}r_{j}(\omega^{lj}+\omega^{-lj})$ (3.10)
$\displaystyle=r_{0}+\sum_{j=1}^{k}2r_{j}\cos\Bigl{(}\frac{\pi lj}{n}\Bigr{)}$
(3.11)
$\displaystyle=r_{0}+\sum_{j=1}^{k}2r_{j}\mathopen{}\mathclose{{}\left(1-\frac{\pi^{2}l^{2}j^{2}}{2n^{2}}+\frac{(\pi
c_{j})^{4}}{4!n^{4}}\cos\mathopen{}\mathclose{{}\left(\frac{\pi
c_{j}}{n}}\right)}\right)$ (3.12)
$\displaystyle=-\frac{\pi^{2}l^{2}}{n^{2}}\sum_{j=1}^{k}r_{j}j^{2}+O\mathopen{}\mathclose{{}\left(\frac{l^{4}k^{3}}{n^{4}}}\right)$
(3.13)

where the $c_{j}\in[0,lj]$ arise from the Taylor remainder theorem. The last
line follows from the fact that $|r_{j}|=O(1/j^{2})$.

We now compute the sum

$\displaystyle-\sum_{j=1}^{k}r_{j}j^{2}$
$\displaystyle=\sum_{j=1}^{k}j^{2}\frac{2(-1)^{j}(k!)^{2}}{j^{2}(k+j)!(k-j)!}$
(3.14) $\displaystyle=2(k!)^{2}\sum_{j=1}^{k}\frac{(-1)^{j}}{(k+j)!(k-j)!}$
(3.15)
$\displaystyle=\frac{2(k!)^{2}}{(2k)!}\sum_{j=1}^{k}(-1)^{j}\binom{2k}{k+j}$
(3.16)
$\displaystyle=\frac{2(k!)^{2}}{(2k)!}\sum_{j=k+1}^{2k}(-1)^{j+k}\binom{2k}{j}$
(3.17) $\displaystyle=(-1)^{k}\frac{(k!)^{2}}{(2k)!}\sum_{j=0,\,j\neq
k}^{2k}(-1)^{j}\binom{2k}{j}$ (3.18)
$\displaystyle=(-1)^{k}\frac{(k!)^{2}}{(2k)!}\biggl{(}(1-1)^{2k}-(-1)^{k}\binom{2k}{k}\biggr{)}$
(3.19) $\displaystyle=-1.$ (3.20)

In particular, we have

$\displaystyle\lambda_{1}$
$\displaystyle=-\frac{\pi^{2}}{n^{2}}+O\biggl{(}\frac{k^{3}}{n^{4}}\biggr{)}$
(3.21)

and the error term is vanishingly small for $k=o(n^{2/3})$. ∎

In $d$ dimensions, a similar analysis holds.

######  Lemma 3.

For $k=o(n^{2/3})$, let $L:=r_{0}I+\sum_{j=1}^{k}r_{j}(S^{j}+S^{-j})$ with
$r_{j}$ as in (3.2). The matrix $L^{\prime}:=L\otimes I^{\otimes d-1}+I\otimes
L\otimes I^{\otimes d-2}+\dots+I^{\otimes d-1}\otimes L$ has condition number
$\kappa(L^{\prime})=\Theta(dn^{2})$.

###### Proof.

By the triangle inequality for spectral norms, $\|{L^{\prime}}\|\leq
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^{\prime}$. ∎

###  3.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 $\vec{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\rangle$
be the quantum state that is proportional to
$\vec{u}=\sum_{j\in\mathbb{Z}_{2n}^{d}}u(\pi
j/n)\bigotimes_{i=1}^{d}e_{j_{i}}$ for the exact solution of the differential
equation. Let $|\bar{u}\rangle$ be the state output by a QLSA that exactly
solves the linear system. Let $|\tilde{u}\rangle$ be the state output by a
QLSA with error. Then the total error of approximating $|u\rangle$ by
$|\tilde{u}\rangle$ is bounded by

$\displaystyle\|{|u\rangle-|\tilde{u}\rangle}\|$
$\displaystyle\leq\|{|u\rangle-|\bar{u}\rangle}\|+\|{|\bar{u}\rangle-|\tilde{u}\rangle}\|$
(3.22) $\displaystyle=\epsilon_{\mathrm{FDM}}+\epsilon_{\mathrm{QLSA}}$ (3.23)

and without loss of generality we can take $\epsilon_{\mathrm{FDM}}$ and
$\epsilon_{\mathrm{QLSA}}$ to be of the same order of magnitude.

######  Lemma 4.

Let $u(\vec{x})$ be the exact solution of
$(\sum_{i=1}^{d}\frac{\mathrm{d}^{2}}{\mathrm{d}x_{i}^{2}})u(\vec{x})=f(\vec{x})$.
Let $\vec{u}\in\mathbb{R}^{(2n)^{d}}$ encode the exact solution in the sense
that $\vec{u}=\sum_{j\in\mathbb{Z}_{2n}^{d}}u(\pi
j/n)\bigotimes_{i=1}^{d}e_{j_{i}}$. Let $\bar{u}\in\mathbb{R}^{(2n)^{d}}$ be
the exact solution of the FDM linear system
$\frac{1}{h^{2}}L^{\prime}\bar{u}=\vec{f}$, where $L^{\prime}$ is a
$d$-dimensional $(2k)$th-order Laplacian as above where $k=o(n^{2/3})$ and
$\vec{f}=\sum_{j=1}^{2n}f(\pi j/n)\vec{e}_{j}$. Then
$\|{\vec{u}-\bar{u}}\|\leq
O(2^{d/2}n^{(d/2)-2k+1}\bigl{|}\frac{\mathrm{d}^{2k+1}u}{\mathrm{d}x^{2k+1}}\bigr{|}(e^{2}/4)^{k})$.

###### Proof.

The remainder term of the central difference formula is
$O(\bigl{|}\frac{\mathrm{d}^{2k+1}u}{\mathrm{d}x^{2k+1}}\bigr{|}h^{2k-1}(e/2)^{2k})$,
so

$\displaystyle\frac{1}{h^{2}}L^{\prime}\vec{u}$
$\displaystyle=\vec{f}+O\Bigl{(}\Bigl{|}\frac{\mathrm{d}^{2k+1}u}{\mathrm{d}x^{2k+1}}\Bigr{|}(eh/2)^{2k-1}\Bigr{)}\vec{\epsilon}$
(3.24)

where $\vec{\epsilon}$ is a $(2n)^{d}$ dimensional vector whose entries are
$O(1)$. This implies

$\displaystyle\frac{1}{h^{2}}L^{\prime}(\vec{u}-\bar{u})$
$\displaystyle=O\Bigl{(}\Bigl{|}\frac{\mathrm{d}^{2k+1}u}{\mathrm{d}x^{2k+1}}\Bigr{|}(eh/2)^{2k-1}\Bigr{)}\vec{\epsilon}$
(3.25)

and therefore

$\displaystyle\|{\vec{u}-\bar{u}}\|$
$\displaystyle=O\Bigl{(}\Bigl{|}\frac{\mathrm{d}^{2k+1}u}{\mathrm{d}x^{2k+1}}\Bigr{|}(eh/2)^{2k+1}\Bigr{)}\|{(L^{\prime})^{-1}\vec{\epsilon}}\|$
(3.26)
$\displaystyle=O\Bigl{(}(2n)^{d/2}\Bigl{|}\frac{\mathrm{d}^{2k+1}u}{\mathrm{d}x^{2k+1}}\Bigr{|}(eh/2)^{2k+1}/\lambda_{1}\Bigr{)}.$
(3.27)

By Lemma 2 we have $\lambda_{1}=\Theta(1/n^{2})$, and since $h=\Theta(1/n)$,
we have

$\displaystyle\|{\vec{u}-\bar{u}}\|$
$\displaystyle=O\Bigl{(}2^{d/2}n^{(d/2)-2k+1}\Bigl{|}\frac{\mathrm{d}^{2k+1}u}{\mathrm{d}x^{2k+1}}\Bigr{|}(e/2)^{2k}\Bigr{)}$
(3.28)

as claimed. ∎

###  3.4 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 multidimensional Laplacian $L^{\prime}$), i.e.,
$D=F^{\dagger}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 [26] 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 Ref. [34] or Ref. [3] could
also be used). Using this Hamiltonian simulation algorithm in a QLSA for the
FDM linear system gives us the following theorem.

######  Theorem 1.

(Quantum FDM for Poisson equation with periodic boundary conditions) There
exists a quantum algorithm that outputs a state $\epsilon$-close to
$|u\rangle$ that runs in time

$\displaystyle\tilde{O}(d^{2}\log^{9/2+\gamma}(\sqrt{d}\Bigl{|}\frac{\mathrm{d}^{2k+1}u}{\mathrm{d}x^{2k+1}}\Bigr{|}/\epsilon)\sqrt{\log(1/\epsilon)})$
(3.29)

and makes

$\displaystyle\tilde{O}(d\log^{3+\beta}(\sqrt{d}\Bigl{|}\frac{\mathrm{d}^{2k+1}u}{\mathrm{d}x^{2k+1}}\Bigr{|}/\epsilon)\sqrt{\log(1/\epsilon)})$
(3.30)

queries to the oracle for $\vec{f}$ for arbitrarily small $\gamma,\beta>0$.

###### Proof.

We use the Fourier series based QLSA from Ref. [7]. By Theorem 3 of that work,
the QLSA makes $O(\kappa\sqrt{\log(\kappa/\epsilon_{\mathrm{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 Ref. [26], for a total of
$O(dn\kappa\sqrt{\log(\kappa/\epsilon_{\mathrm{QLSA}})})$ gates. The condition
number for the $d$-dimensional Laplacian scales as $\kappa=O(dn^{2})$.

We take $\epsilon_{\mathrm{FDM}}$ and $\epsilon_{\mathrm{QLSA}}$ to be of the
same order and just write $\epsilon$. Then the QLSA has time complexity
$O(n^{3}\sqrt{\log(n^{2}/\epsilon)})$ and query complexity
$O(n^{2}\log(n^{2}/\epsilon))$. The adjustable parameters are the number of
lattice sites $n$ and the order of the finite difference formula. To keep the
error below the target error of $\epsilon$ we require

$\displaystyle
2^{d/2}n^{(d/2)-2k+1}\Bigl{|}\frac{\mathrm{d}^{2k+1}u}{\mathrm{d}x^{2k+1}}\Bigr{|}(e/2)^{2k}$
$\displaystyle=O(\epsilon),$ (3.31)

or equivalently,

$\displaystyle(-d/2)+(-(d/2)+2k-1)\log(n)-2k\log(e/2)$
$\displaystyle=\Omega\Bigl{(}\log\Bigl{(}\Bigl{|}\frac{\mathrm{d}^{2k+1}u}{\mathrm{d}x^{2k+1}}\Bigr{|}/\epsilon\Bigr{)}\Bigr{)}.$
(3.32)

We can satisfy this condition by setting $k=dn^{b}$ for some constant $b<2/3$
(so that our condition number lemmas hold) and taking $n$ sufficiently large,
giving

$\displaystyle n^{b}\log(n)$
$\displaystyle=\Omega\Bigl{(}\frac{1}{d}\log\Bigl{(}\Bigl{|}\frac{\mathrm{d}^{2k+1}u}{\mathrm{d}x^{2k+1}}\Bigr{|}/\epsilon\Bigr{)}\Bigr{)}.$
(3.33)

We can set
$n=\Theta(\log^{3/2+\delta}(\bigl{|}\frac{\mathrm{d}^{2k+1}u}{\mathrm{d}x^{2k+1}}\bigr{|}/\epsilon))$
to satisfy this criterion for arbitrary ($b$-dependent) $\delta>0$. The QLSA
then has the stated running time and query complexity. ∎

This can be compared to the cost of using the conjugate gradient method to
solve the same linear system classically. The sparse conjugate gradient
algorithm for an $N\times N$ matrix has time complexity
$O(Ns\sqrt{\kappa}\log(1/\epsilon))$. For arbitrary dimension
$N=\Theta(n^{d})$, we have $s=dk=d^{2}n^{b}$ and $\kappa=O(dn^{2})$, so that
the time complexity is
$O(d^{2.5}\log(1/\epsilon)\log^{3+1.5d+1.5b+\gamma^{\prime}}(\sqrt{d}\bigl{|}\frac{\mathrm{d}^{2k+1}u}{\mathrm{d}x^{2k+1}}\bigr{|}/\epsilon))$
for arbitrary $\gamma^{\prime}>0$. Alternatively, $d$ fast Fourier transforms
could be used, although this will generally take
$\Omega(n^{d})=\Omega(\log^{3d/2}(\sqrt{d}\bigl{|}\frac{\mathrm{d}^{2k+1}u}{\mathrm{d}x^{2k+1}}\bigr{|}/\epsilon))$
time.

###  3.5 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 $\mathbb{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
$-1$ (i.e., $f(-x-1)=-f(x-1)$) and anti-symmetric about 1 (i.e.,
$f(1+x)=-f(1-x)$). Continuity and anti-symmetry about $-1$ and $1$ imply
$f(-1)=f(1)=0$, and furthermore that $f(x)=0$ for all odd $x\in\mathbb{Z}$ and
that $f(x+4)=f(x)$ for all $x\in{\mathbb{R}}$. For Neumann boundary
conditions, the solutions are instead symmetric about $-1$ and $1$, which also
implies $f(x+2)=f(x)$ for all $x\in{\mathbb{R}}$.

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\pi)$ instead of a lattice on $\mathbb{R}$. It is
useful to think of this lattice in terms of the cycle graph on $2n$ vertices,
i.e., $(V,E)=(\mathbb{Z}_{2n},\\{(i,i+1)\mid i\in\mathbb{Z}_{2n}\\})$, which
means that the vectors encoding the solution $u(x)$ will lie in
$\mathbb{R}^{2n}$. Let each vector $\vec{e}_{j}$ correspond to the vertex $j$.
Then we divide $\mathbb{R}^{2n}$ into a symmetric and an anti-symmetric
subspace, namely $\mathrm{span}\\{e_{j}+e_{2n+1-j}\\}_{j=1}^{n}$ and
$\mathrm{span}\\{e_{j}-e_{2n+1-j}\\}_{j=1}^{n}$, respectively. Vectors lying
in the symmetric subspace correspond to solutions that are symmetric about $0$
and $\pi$, so they obey Neumann boundary conditions at $0$ and $\pi$;
similarly, vectors in the anti-symmetric space correspond to solutions obeying
Dirichlet boundary conditions at $0$ and $\pi$.

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^{\prime\prime}\vec{u}^{\prime\prime}=\vec{f}^{\prime\prime}$ then
$L^{\prime\prime}$ has entries

$\displaystyle L^{\prime\prime}_{i,j}=\begin{cases}r_{|i-j|}\pm
r_{i+j-1}&i\leq k\\\ r_{|i-j|}&k<i\leq n-k\\\ r_{|i-j|}\pm
r_{2n-i-j+1}&n-k\leq i\end{cases}$ (3.34)

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 [10].

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^{\prime\prime}\rangle$
proportional to the inhomogeneity $f^{\prime\prime}(x)$. Then we apply the
QLSA for periodic boundary conditions using
$|f^{\prime\prime}\rangle|\pm\rangle$ to encode the inhomogeneity, which will
output solutions of the form $|u^{\prime\prime}\rangle|\pm\rangle$. Here the
ancillary state is chosen to be $|+\rangle$ ($|-\rangle$) 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 Ref. [32] 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
$\mathrm{span}(\\{e_{j}-e_{2n+2-j}\\}_{j=1}^{n})$ (which has no support on
$e_{0}$ and $e_{n+1}$).

If we again write the linear system as
$L^{\prime\prime}\vec{u}^{\prime\prime}=\vec{f}^{\prime\prime}$, then the
Laplacian has entries

$L^{\prime\prime}_{i,j}=\begin{cases}r_{|i-j|}-r_{i+j}&i\leq k\\\
r_{|i-j|}&k<i\leq n-k\\\ r_{|i-j|}-r_{2n-i-j+2}&n-k\leq i.\end{cases}$

We again assume access to an oracle producing states proportional to
$f^{\prime\prime}(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 $\Bigl{(}\begin{smallmatrix}U&\vec{0}\\\
\vec{0}^{T}&1\end{smallmatrix}\Bigr{)}$). With this oracle we again prepare
the state $|f^{\prime\prime}\rangle|-\rangle$ and solve Poisson’s equation for
periodic boundary conditions to output a state
$|u^{\prime\prime}\rangle|-\rangle$ (where $|u^{\prime\prime}\rangle$ lies in
an ($n+1$)-dimensional Hilbert space but has no support on the $(n+1)$st basis
state).

##  4 Multi-dimensional spectral method

We now turn our attention to the spectral method for multi-dimensional PDEs.
Since interpolation facilitates constructing a straightforward linear system,
we develop a quantum algorithm based on the _pseudo-spectral method_ [14, 27]
with various boundary conditions. After introducing the method, we discuss the
quantum shifted Fourier transform and the quantum cosine transform (Section
4.1), which are used as subroutines in our algorithm. Then we construct a
linear system whose solution encodes the solution of the PDE (Section 4.2),
analyze its condition number (Section 4.3), and consider the complexity of
state preparation (Section 4.4). Finally, we present our main result in
Section 4.5.

In the spectral approach, we approximate the exact solution $\hat{u}(\bm{x})$
by a linear combination of basis functions

$u(\bm{x})=\sum_{\|\bm{k}\|_{\infty}\leq n}c_{\bm{k}}\phi_{\bm{k}}(\bm{x})$
(4.1)

for some $n\in{\mathbb{Z}}^{+}$. Here $\bm{k}=(k_{1},\ldots,k_{d})$ with
$k_{j}\in[{n+1}]_{0}:=\\{0,1,\ldots,n\\}$, $c_{\bm{k}}\in{\mathbb{C}}$, and

$\phi_{\bm{k}}(\bm{x})=\prod_{j=1}^{d}\phi_{k_{j}}(x_{j}),\quad j\in[{d}].$
(4.2)

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}\in[{n+1}]_{0}$ and $x_{j}\in[-1,1]$,
we take

$\phi_{k_{j}}(x_{j})=\begin{cases}e^{i(k_{j}-\mathopen{}\mathclose{{}\left\lfloor
n/2}\right\rfloor)\pi x_{j}},&\text{periodic conditions},\\\
T_{k_{j}}(x_{j}):=\cos(k_{j}\arccos x_{j}),&\text{non-periodic
conditions}.\end{cases}$ (4.3)

Here $T_{k}$ is the degree-$k$ Chebyshev polynomial of the first kind.

The coefficients $c_{\bm{k}}$ are determined by demanding that $u(\bm{x})$
satisfies the ODE and boundary conditions at a set of _interpolation nodes_
$\\{\bm{\chi}_{\bm{l}}=(\chi_{l_{1}},\ldots,\chi_{l_{d}})\\}_{\|\bm{l}\|_{\infty\leq
n}}$ with $l_{j}\in[{n+1}]_{0}$, where

$\chi_{l_{j}}=\begin{cases}\frac{2l_{j}}{n+1}-1,&\text{periodic
conditions},\\\ \cos\frac{\pi l_{j}}{n},&\text{non-periodic
conditions}.\end{cases}$ (4.4)

Here $\\{\frac{2l}{n+1}-1:l\in[{n+1}]_{0}\\}$ are called the _uniform grid
nodes_ , and $\\{\cos\frac{\pi l}{n}:l\in[{n+1}]_{0}\\}$ are called the
_Chebyshev-Gauss-Lobatto quadrature nodes_.

We require the numerical solution $u(\bm{x})$ to satisfy

$\mathscr{L}(u(\bm{\chi}_{\bm{l}}))=f(\bm{\chi}_{\bm{l}}),\quad\forall\,l_{j}\in[{n+1}]_{0},~{}j\in[{d}].$
(4.5)

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

$\|\hat{u}(\bm{x})-u(\bm{x})\|\rightarrow 0\quad\text{as}\quad
n\rightarrow\infty.$ (4.6)

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 in $C^{\infty}$, the spectral method approximates the solution to
within $\epsilon$ using only $n=\operatorname{poly}(\log(1/\epsilon))$ [27].
Since we require $k_{j}\in[{n+1}]_{0}$ for all $j\in[{d}]$, we have
$(n+1)^{d}$ terms in total. Consequently, a classical pseudo-spectral method
solves multi-dimensional PDEs with complexity
$\operatorname{poly}(\log^{d}(1/\epsilon))$. 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 [8]. 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.

###  4.1 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
$v=(v_{0},v_{1},\ldots,v_{n})\in{\mathbb{C}}^{n+1}$ to the state
$\hat{v}=(\hat{v}_{0},\hat{v}_{1},\ldots,\hat{v}_{n})\in{\mathbb{C}}^{n+1}$
with

$\hat{v}_{l}=\frac{1}{\sqrt{n+1}}\sum_{k=0}^{n}\exp\Bigl{(}\frac{2\pi
ikl}{n+1}\Bigr{)}v_{k},\quad l\in[{n+1}]_{0}.$ (4.7)

In other words, the QFT is the unitary transform

$F_{n}:=\frac{1}{\sqrt{n+1}}\sum_{k,l=0}^{n}\exp\Bigl{(}\frac{2\pi
ikl}{n+1}\Bigr{)}|l\rangle\langle k|.$ (4.8)

Here we also consider the _quantum shifted Fourier transform_ (QSFT), an
analogue of the classical shifted discrete Fourier transform, which maps
$v\in{\mathbb{C}}^{n+1}$ to $\hat{v}\in{\mathbb{C}}^{n+1}$ with

$\hat{v}_{l}=\frac{1}{\sqrt{n+1}}\sum_{k=0}^{n}\exp\Bigl{(}\frac{2\pi
i(k-\mathopen{}\mathclose{{}\left\lfloor
n/2}\right\rfloor)(l-(n+1)/2)}{n+1}\Bigr{)}v_{k},\quad l\in[{n+1}]_{0}.$ (4.9)

In other words, the QSFT is the unitary transform

$F^{s}_{n}:=\frac{1}{\sqrt{n+1}}\sum_{k,l=0}^{n}\exp\Bigl{(}\frac{2\pi
i(k-\mathopen{}\mathclose{{}\left\lfloor
n/2}\right\rfloor)(l-(n+1)/2)}{n+1}\Bigr{)}|l\rangle\langle k|.$ (4.10)

We define the multi-dimensional QSFT by the tensor product, namely

${\bm{F}}^{s}_{n}:=\frac{1}{\sqrt{(n+1)^{d}}}\sum_{\|\bm{k}\|_{\infty},\|\bm{l}\|_{\infty}\leq
n}\prod_{j=1}^{d}\exp\Bigl{(}\frac{2\pi
i(k_{j}-\mathopen{}\mathclose{{}\left\lfloor
n/2}\right\rfloor)(l_{j}-(n+1)/2)}{n+1}\Bigr{)}|l_{1}\rangle\ldots|l_{d}\rangle\langle
k_{1}|\ldots\langle k_{d}|,$ (4.11)

where $\bm{k}=(k_{1},\ldots,k_{d})$ and $\bm{l}=(l_{1},\ldots,l_{d})$ are
$d$-dimensional vectors with $k_{j},l_{j}\in[{n}]_{0}$.

The QSFT can be efficiently implemented as follows:

######  Lemma 5.

The QSFT $F_{n}^{s}$ defined by (4.10) can be performed with gate complexity
$O(\log n\log\log n)$. More generally, the $d$-dimensional QSFT
${\bm{F}}_{n}^{s}$ defined by (4.11) can be performed with gate complexity
$O(d\log n\log\log n)$.

###### Proof.

The unitary matrix $F^{s}_{n}$ can be written as the product of three unitary
matrices

$F^{s}_{n}=S_{n}F_{n}R_{n},$ (4.12)

where

$R_{n}=\sum_{k=0}^{n}\exp\Bigl{(}-\frac{2\pi
ik(n+1)/2}{n+1}\Bigr{)}|k\rangle\langle k|$ (4.13)

and

$S_{n}=\sum_{l=0}^{n}\exp\Bigl{(}-\frac{2\pi
i\mathopen{}\mathclose{{}\left\lfloor
n/2}\right\rfloor(l-(n+1)/2)}{n+1}\Bigr{)}|l\rangle\langle l|.$ (4.14)

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)$.

We rewrite $\bm{v}$ in the form

$\bm{v}=\sum_{\|\bm{k}\|_{\infty}\leq
n}\bm{v}_{\bm{k}}|k_{1}\rangle\ldots|k_{d}\rangle,$ (4.15)

where $\bm{v}_{\bm{k}}\in{\mathbb{C}}$ with $\bm{k}=(k_{1},\ldots,k_{d})$, and
each $k_{j}\in[{n}]_{0}$ for $j\in[{d}]$. The unitary matrix
${\bm{F}}^{s}_{n}$ can be written as the tensor product

${\bm{F}}^{s}_{n}=\bigotimes_{j=1}^{d}F^{s}_{n}.$ (4.16)

Performing the multi-dimensional QSFT is equivalent to performing the one-
dimensional QSFT on each register. Thus, the gate complexity of performing
${\bm{F}}^{s}_{n}$ is $O(d\log n\log\log n)$. ∎

Another efficient quantum transformation is the _quantum cosine transform_
(QCT) [18, 25]. The QCT can be regarded as an analogue of the discrete cosine
transform (DCT). The QCT maps $v\in{\mathbb{C}}^{n+1}$ to
$\hat{v}\in{\mathbb{C}}^{n+1}$ with

$\hat{v}_{l}=\sqrt{\frac{2}{n}}\sum_{k=0}^{n}\delta_{k}\delta_{l}\cos\frac{kl\pi}{n}v_{k},\quad
l\in[{n+1}]_{0},$ (4.17)

where

$\delta_{l}:=\begin{cases}\frac{1}{\sqrt{2}}&l=0,n\\\
1&l\in[{n-1}].\end{cases}$ (4.18)

In other words, the QCT is the orthogonal transform

$C_{n}:=\sqrt{\frac{2}{n}}\sum_{k,l=0}^{n}\delta_{l}\delta_{k}\cos\frac{kl\pi}{n}|l\rangle\langle
k|.$ (4.19)

Again we define the multi-dimensional QCT by the tensor product, namely

${\bm{C}}_{n}:=\sqrt{\Bigl{(}\frac{2}{n}\Bigr{)}^{d}}\sum_{\|\bm{k}\|_{\infty},\|\bm{l}\|_{\infty}\leq
n}\prod_{j=1}^{d}\delta_{k_{j}}\delta_{l_{j}}\cos\frac{k_{j}l_{j}\pi}{n}|l_{1}\rangle\ldots|l_{d}\rangle\langle
k_{1}|\ldots\langle k_{d}|,$ (4.20)

where $\bm{k}=(k_{1},\ldots,k_{d})$ and $\bm{l}=(l_{1},\ldots,l_{d})$ are
$d$-dimensional vectors with $k_{j},l_{j}\in[{n+1}]_{0}$.

The classical DCT on $(n+1)$-dimensional vectors takes $\Theta(n\log n)$
gates, while the QCT on $(n+1)$-dimensional quantum states can be implemented
with complexity $\operatorname{poly}(\log n)$. According to [18, Theorem 1],
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.19) can be performed with
gate complexity $O(\log n\log\log n)$. More generally, the multi-dimensional
QCT $\bm{C}_{n}$ defined by (4.20) can be performed with gate complexity
$O(d\log n\log\log n)$.

###### Proof.

According to the quantum circuit in Figure 2 of [18], $C_{n}$ can be
decomposed into a QFT $F_{n+1}$, a permutation

$P_{n}=\begin{pmatrix}&&&&1\\\ 1&&&&\\\ &1&&&\\\ &&\ddots&&\\\ &&&1&\\\
\end{pmatrix},$ (4.21)

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 [23].

The permutation $P_{n}$ can be decomposed as

$P_{n}=F_{n}T_{n}F_{n}^{-1},$ (4.22)

where $F_{n}$ is the Fourier transform (4.8) and
$T_{n}=\sum_{k=0}^{n}e^{-\frac{2\pi ik}{n+1}}|k\rangle\langle 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 ${\bm{C}}_{n}$ can be written as the tensor product

${\bm{C}}_{n}=\bigotimes_{j=1}^{d}C_{n}.$ (4.23)

As in Lemma 5, performing the multi-dimensional QCT is equivalent to
performing a QCT on each register. Thus, the gate complexity of performing
${\bm{C}}_{n}$ is $O(d\log n\log\log n)$. ∎

###  4.2 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

$u(\bm{x}+2\bm{v})=u(\bm{x})\quad\forall\,\bm{x}\in{\mathscr{D}},~{}\forall\,\bm{v}\in{\mathbb{Z}}^{d}$
(4.24)

or non-periodic Dirichlet boundary conditions

$u(\bm{x})=\gamma(\bm{x})\quad\forall\,\bm{x}\in\partial{\mathscr{D}}.$ (4.25)

According to the elliptic regularity theorem (Theorem 6 in Section 6.3 of
[12]), there exists a unique solution $\hat{u}(\bm{x})$ in $C^{\infty}$ 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

$|u\rangle\propto\sum_{\|\bm{k}\|_{\infty},\|\bm{l}\|_{\infty}\leq
n}c_{\bm{k}}\phi_{\bm{k}}(\bm{\chi}_{\bm{l}})|l_{1}\rangle\ldots|l_{d}\rangle,$
(4.26)

where $\phi_{\bm{k}}(\bm{\chi}_{\bm{l}})$ is defined by (4.2) using (4.3) for
the appropriate boundary conditions (periodic or non-periodic). This state
corresponds to a truncated Fourier/Chebyshev approximation and is
$\epsilon$-close to the exact solution $\hat{u}(\bm{\chi}_{\bm{l}})$ with
$n=\operatorname{poly}(\log(1/\epsilon))$ [27]. Note that this state encodes
the values of the solution at the interpolation nodes (4.4) 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\rangle\propto\sum_{\|\bm{k}\|_{\infty}\leq
n}c_{\bm{k}}|k_{1}\rangle\ldots|k_{d}\rangle$ (4.27)

that is the inverse QSFT/QCT of $|u\rangle$. We then apply the QSFT/QCT to
transform back into the interpolation node basis.

The truncated spectral series of the inhomogeneity $f(\bm{x})$ and the
boundary conditions $\gamma(\bm{x})$ can be expressed as

$f(\bm{x})=\sum_{\|\bm{k}\|_{\infty}\leq
n}\hat{f}_{\bm{k}}\phi_{\bm{k}}(\bm{x})$ (4.28)

and

$\gamma(\bm{x})=\sum_{\|\bm{k}\|_{\infty}\leq
n}\hat{\gamma}_{\bm{k}}\phi_{\bm{k}}(\bm{x}),$ (4.29)

respectively. We define quantum states $|f\rangle$ and $|\gamma\rangle$ by
interpolating the nodes $\\{\bm{\chi}_{\bm{l}}\\}$ defined by (4.4) as

$|f\rangle\propto\sum_{\|\bm{k}\|_{\infty},\|\bm{l}\|_{\infty}\leq
n}\phi_{k_{j}}(\bm{\chi}_{\bm{l}})\hat{f}_{\bm{k}}|l_{1}\rangle\ldots|l_{d}\rangle,$
(4.30)

and

$|\gamma\rangle\propto\sum_{\|\bm{k}\|_{\infty},\|\bm{l}\|_{\infty}\leq
n}\phi_{k_{j}}(\bm{\chi}_{\bm{l}})\hat{\gamma}_{\bm{k}}|l_{1}\rangle\ldots|l_{d}\rangle,$
(4.31)

respectively. These are the states that we assume we can produce using
oracles. We perform the multidimensional inverse QSFT/QCT to obtain the states

$|\hat{f}\rangle\propto\sum_{\|\bm{k}\|_{\infty}\leq
n}\hat{f}_{\bm{k}}|k_{1}\rangle\ldots|k_{d}\rangle,$ (4.32)

and

$|\hat{\gamma}\rangle\propto\sum_{\|\bm{k}\|_{\infty}\leq
n}\hat{\gamma}_{\bm{k}}|k_{1}\rangle\ldots|k_{d}\rangle.$ (4.33)

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\vec{x}=\vec{f}$ (where $\vec{x}$ corresponds to (4.27)) describes the
differential equation, and another system $B\vec{x}=\vec{g}$ describes the
boundary conditions. We combine these into a linear system with the form

$L\vec{x}=(A+B)\vec{x}=\vec{f}+\vec{g}.$ (4.34)

Even though we do not impose the two linear systems separately, we show that
there exists a unique solution of (4.34) (which is therefore the solution of
the simultaneous equations $A\vec{x}=\vec{f}$ and $B\vec{x}=\vec{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

$\mathscr{L}(u(\bm{\chi}_{\bm{l}}))=\sum_{\|\bm{j}\|_{1}=2}A_{\bm{j}}\frac{\partial^{\bm{j}}}{\partial\bm{x}^{\bm{j}}}u(\bm{\chi}_{\bm{l}})=f(\bm{\chi}_{\bm{l}}),$
(4.35)

while another part will come from imposing the boundary conditions on
$\partial{\mathscr{D}}=\bigcup_{j\in[d]}\partial{\mathscr{D}}_{j}$, where
$\partial{\mathscr{D}}_{j}:=\\{\bm{x}\in{\mathscr{D}}\mid x_{j}=\pm 1\\}$ is a
$(d-1)$-dimensional subspace. More specifically, the boundary conditions

$u(\bm{\chi}_{\bm{l}})=\gamma(\bm{\chi}_{\bm{l}})\quad\forall\,\bm{\chi}_{\bm{l}}\in\partial{\mathscr{D}}$
(4.36)

can be expressed as conditions on each boundary:

$\displaystyle
u(x_{1},\ldots,x_{j-1},1,x_{j+1},\ldots,x_{d})=\gamma^{j+},\quad\bm{x}\in\partial{\mathscr{D}}_{j},\quad
j\in[d]$ (4.37) $\displaystyle
u(x_{1},\ldots,x_{j-1},-1,x_{j+1},\ldots,x_{d})=\gamma^{j-},\quad\bm{x}\in\partial{\mathscr{D}}_{j},\quad
j\in[d].$

####  4.2.1 Linear system from the differential equation

To evaluate the matrix corresponding to the differential operator from (4.35),
it is convenient to define coefficients $c^{(\bm{j})}_{\bm{k}}$ and
$\|\bm{k}\|_{\infty}\leq n$ such that

$\frac{\partial^{\bm{j}}}{\partial\bm{x}^{\bm{j}}}u(\bm{x})=\sum_{\|\bm{k}\|_{\infty}\leq
n}c^{(\bm{j})}_{\bm{k}}\phi_{\bm{k}}(\bm{x})$ (4.38)

for some fixed $\bm{j}\in{\mathbb{N}}^{d}$ (as we explain below, such a
decomposition exists for the choices of basis functions in (4.3)). Using this
expression, we obtain the following linear equations for
$c^{(\bm{j})}_{\bm{k}}$:

$\sum_{\|\bm{j}\|_{1}=2}A_{\bm{j}}\sum_{\|\bm{k}\|_{\infty},\|\bm{l}\|_{\infty}\leq
n}\phi_{\bm{k}}(\bm{\chi}_{\bm{l}})c^{(\bm{j})}_{\bm{k}}|l_{1}\rangle\ldots|l_{d}\rangle=\sum_{\|\bm{k}\|_{\infty},\|\bm{l}\|_{\infty}\leq
n}\phi_{\bm{k}}(\bm{\chi}_{\bm{l}})\hat{f}_{\bm{k}}|l_{1}\rangle\ldots|l_{d}\rangle.$
(4.39)

To determine the transformation between $c^{(\bm{j})}_{\bm{k}}$ and
$c_{\bm{k}}$, we can make use of the differential properties of Fourier and
Chebyshev series, namely

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}x}e^{ik\pi x}=ik\pi e^{ik\pi x}$
(4.40)

and

$\displaystyle
2T_{k}(t)=\frac{T^{\prime}_{k+1}(t)}{k+1}-\frac{T^{\prime}_{k-1}(t)}{k-1},$
(4.41)

respectively. We have

$c^{(\bm{j})}_{\bm{k}}=\sum_{\|\bm{r}\|_{\infty}\leq
n}[\bm{D}^{(\bm{j})}_{n}]_{\bm{kr}}c_{\bm{r}.},\quad\|\bm{k}\|_{\infty}\leq
n,$ (4.42)

where $\bm{D}^{(\bm{j})}_{n}$ can be expressed as the tensor product

$\bm{D}^{(\bm{j})}_{n}=D^{j_{1}}_{n}\otimes D^{j_{2}}_{n}\otimes\cdots\otimes
D^{j_{d}}_{n},$ (4.43)

with $\bm{j}=(j_{1},\ldots,j_{d})$. The matrix $D_{n}$ for the Fourier basis
functions in (4.3) can be written as the $(n+1)\times(n+1)$ diagonal matrix
with entries

$[D_{n}]_{kk}=i(k-\mathopen{}\mathclose{{}\left\lfloor n/2}\right\rfloor)\pi.$
(4.44)

As detailed in Appendix A of [8], the matrix $D_{n}$ for the Chebyshev
polynomials in (4.3) can be expressed as the $(n+1)\times(n+1)$ upper
triangular matrix with nonzero entries

$[D_{n}]_{kr}=\frac{2r}{\sigma_{k}},\qquad\text{$k+r$ odd},~{}r>k,$ (4.45)

where

$\sigma_{k}:=\begin{cases}2&k=0\\\ 1&k\in[{n}].\end{cases}$ (4.46)

Substituting (4.43) into (4.39), with $D_{n}$ defined by (4.44) in the
periodic case or (4.45) in the non-periodic case, and performing the
multidimensional inverse QSFT/QCT (for a reason that will be explained in the
next section), we obtain the following linear equations for $c_{\bm{r}}$:

$\sum_{\|\bm{j}\|_{1}=2}A_{\bm{j}}\sum_{\|\bm{k}\|_{\infty},\|\bm{l}\|_{\infty},\|\bm{r}\|_{\infty}\leq
n}[\bm{D}^{(\bm{j})}_{n}]_{\bm{kr}}c_{\bm{r}}|l_{1}\rangle\ldots|l_{d}\rangle=\sum_{\|\bm{k}\|_{\infty},\|\bm{l}\|_{\infty}\leq
n}\hat{f}_{\bm{k}}|l_{1}\rangle\ldots|l_{d}\rangle.$ (4.47)

Notice that the matrices (4.44) and (4.45) are not full rank. More
specifically, there exists at least one zero row in the matrix of (4.47) when
using either (4.44) ($k=\mathopen{}\mathclose{{}\left\lfloor
n/2}\right\rfloor$) or (4.45) ($k=n$). To obtain an invertible linear system,
we next introduce the boundary conditions.

####  4.2.2 Adding the linear system from the boundary conditions

When we use the form (4.26) of $u(\bm{x})$ to write linear equations
describing the boundary conditions (4.37), we obtain a non-sparse linear
system. Thus, for each $\bm{x}\in\partial{\mathscr{D}}_{j}$ in (4.37), we
perform the $(d-1)$-dimensional inverse QSFT/QCT on the $d-1$ registers except
the $j$th register to obtain the linear equations

$\displaystyle\sum_{\begin{subarray}{c}\|\bm{k}\|_{\infty}\leq n\\\
k_{j}=n\end{subarray}}c_{\bm{k}}|k_{1}\rangle\ldots|k_{d}\rangle$
$\displaystyle=\sum_{\begin{subarray}{c}\|\bm{k}\|_{\infty}\leq n\\\
k_{j}=n\end{subarray}}\hat{\gamma}^{1+}_{\bm{k}}|k_{1}\rangle\ldots|k_{d}\rangle,\quad\hat{\gamma}^{j+}_{\bm{k}}\in\partial{\mathscr{D}}_{j},$
(4.48) $\displaystyle\sum_{\begin{subarray}{c}\|\bm{k}\|_{\infty}\leq n\\\
k_{j}=n-1\end{subarray}}(-1)^{k_{j}}c_{\bm{k}}|k_{1}\rangle\ldots|k_{d}\rangle$
$\displaystyle=\sum_{\begin{subarray}{c}\|\bm{k}\|_{\infty}\leq n\\\
k_{j}=n-1\end{subarray}}\hat{\gamma}^{1-}_{\bm{k}}|k_{1}\rangle\ldots|k_{d}\rangle,\quad\hat{\gamma}^{j-}_{\bm{k}}\in\partial{\mathscr{D}}_{j}$

for all $j\in[d]$, where the values of $k_{j}$ indicate that we place these
constraints in the last two rows with respect to the $j$th coordinate. We
combine these equations with (4.47) to obtain the linear system

$\sum_{\|\bm{j}\|_{1}=2}A_{\bm{j}}\sum_{\|\bm{k}\|_{\infty},\|\bm{r}\|_{\infty}\leq
n}[\overline{\bm{D}}^{(\bm{j})}_{n}]_{\bm{kr}}c_{\bm{r}}|k_{1}\rangle\ldots|k_{d}\rangle=\sum_{\|\bm{k}\|_{\infty}\leq
n}\sum_{j=1}^{d}(A_{j,j}\hat{\gamma}^{j+}_{\bm{k}}+A_{j,j}\hat{\gamma}^{j-}_{\bm{k}}+\hat{f}_{\bm{k}})|k_{1}\rangle\ldots|k_{d}\rangle,$
(4.49)

where

$\overline{\bm{D}}^{(\bm{j})}_{n}=\begin{cases}\bm{D}^{(\bm{j})}_{n}+\bm{G}^{(\bm{j})}_{n},\quad&\|\bm{j}\|_{1}=2,~{}\|\bm{j}\|_{\infty}=2;\\\
\bm{D}^{(\bm{j})}_{n},&\|\bm{j}\|_{1}=2,~{}\|\bm{j}\|_{\infty}=1\end{cases}$
(4.50)

with $\bm{G}^{(\bm{j})}_{n}$ defined below. In other words,
$\overline{\bm{D}}^{(\bm{j})}_{n}=\bm{D}^{(\bm{j})}_{n}+\bm{G}^{(\bm{j})}_{n}$
for each $\bm{j}$ that has exactly one entry equal to $2$ and all other
entries $0$, whereas $\overline{\bm{D}}^{(\bm{j})}_{n}=\bm{D}^{(\bm{j})}_{n}$
for each $\bm{j}$ that has exactly two entries equal to $1$ and all other
entries $0$. Here $\bm{G}^{(\bm{j})}_{n}$ can be expressed as the tensor
product

$\bm{G}^{(\bm{j})}_{n}=I^{\otimes r-1}\otimes G_{n}\otimes I^{\otimes d-r}$
(4.51)

where the $r$th entry of $\bm{j}$ is 2 and all other entries are 0. For the
Fourier case in (4.3) used for periodic boundary conditions, $D_{n}$ comes
from (4.44), and the nonzero entries of $G_{n}$ are

$[G_{n}]_{\mathopen{}\mathclose{{}\left\lfloor n/2}\right\rfloor,k}=1,\quad
k\in[{n+1}]_{0}.$ (4.52)

Alternatively, for the Chebyshev case in (4.3) used for non-periodic boundary
conditions, $D_{n}$ comes from (4.45), and the nonzero entries of $G_{n}$ are

$\displaystyle[G_{n}]_{n,k}=1,\quad$ $\displaystyle k\in[{n+1}]_{0},$ (4.53)
$\displaystyle[G_{n}]_{n-1,k}=(-1)^{k},\quad$ $\displaystyle k\in[{n+1}]_{0}.$

The system (4.49) has the form of (4.34). For instance, the matrix in (4.34)
for Poisson’s equation (2.6) is

$\displaystyle L_{\mathrm{Poisson}}$
$\displaystyle:=\overline{\bm{D}}_{n}^{(2,0,\ldots,0)}+\overline{\bm{D}}_{n}^{(0,2,\ldots,0)}+\cdots+\overline{\bm{D}}_{n}^{(0,0,\ldots,2)}$
(4.54)
$\displaystyle=\bigoplus_{j=1}^{d}\overline{D}^{(2)}_{n}=\overline{D}^{(2)}_{n}\otimes
I^{\otimes d-1}+I\otimes\overline{D}^{(2)}_{n}\otimes I^{\otimes
d-2}+\cdots+I^{\otimes d-1}\otimes\overline{D}^{(2)}_{n}.$ (4.55)

For periodic boundary conditions, using (4.42), (4.44), and (4.52), the
second-order differential matrix $\overline{D}^{(2)}_{n}$ has nonzero entries

k,k $\displaystyle=-((k-\mathopen{}\mathclose{{}\left\lfloor
n/2}\right\rfloor)\pi)^{2},$ $\displaystyle k$
$\displaystyle\in[{n+1}]_{0}\backslash\\{\mathopen{}\mathclose{{}\left\lfloor
n/2}\right\rfloor\\},$ (4.56)
$\displaystyle[\overline{D}^{(2)}_{n}]_{\mathopen{}\mathclose{{}\left\lfloor
n/2}\right\rfloor,k}$ $\displaystyle=1,$ $\displaystyle k$
$\displaystyle\in[{n+1}]_{0}.$

For non-periodic boundary conditions, using (4.42), (4.45), and (4.53),
$\overline{D}^{(2)}_{n}$ has nonzero entries

kr $\displaystyle=\sum_{\begin{subarray}{c}l=k+1\\\ \text{$k+l$ odd}\\\
\text{$l+r$
odd}\end{subarray}}^{r-1}[D_{n}]_{kl}[D_{n}]_{lr}=\frac{2r}{\sigma_{k}}\sum_{\begin{subarray}{c}l=k+1\\\
\text{$k+l$ odd}\\\ \text{$l+r$
odd}\end{subarray}}^{r-1}\frac{2l}{\sigma_{l}}=\frac{r(r^{2}-k^{2})}{\sigma_{k}},$
$\displaystyle\text{$k+r$ even},\,r>k+1,$ (4.57)
$\displaystyle[\overline{D}^{(2)}_{n}]_{n,k}$ $\displaystyle=1,$
$\displaystyle k\in[{n+1}]_{0},$
$\displaystyle[\overline{D}^{(2)}_{n}]_{n-1,k}$ $\displaystyle=(-1)^{k},$
$\displaystyle k\in[{n+1}]_{0}.$

To explicitly illustrate this linear system, we present a simple example in
Appendix A. We discuss the invertible linear system (4.49) and upper bound its
condition number in the following section.

###  4.3 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.56) and (4.57) that
appear in the linear system.

######  Lemma 7.

Consider the case of periodic boundary conditions. Then the largest and
smallest singular values of $\overline{D}^{(2)}_{n}$ defined in (4.56) satisfy

$\displaystyle\sigma_{\max}(\overline{D}^{(2)}_{n})\leq(2n)^{2.5},$ (4.58)
$\displaystyle\sigma_{\min}(\overline{D}^{(2)}_{n})\geq\frac{1}{\sqrt{2}}.$

######  Lemma 8.

Consider the case of non-periodic boundary conditions. Then the largest and
smallest singular values of $\overline{D}^{(2)}_{n}$ defined in (4.57) satisfy

$\displaystyle\sigma_{\max}(\overline{D}^{(2)}_{n})\leq n^{4},$ (4.59)
$\displaystyle\sigma_{\min}(\overline{D}^{(2)}_{n})\geq\frac{1}{16}.$

The proofs of Lemma 7 and Lemma 8 appear in Appendix B. Using these two
lemmas, we first upper bound the condition number of the linear system for
Poisson’s equation, and then extend the result to general elliptic PDEs.

For the case of the Poisson equation, we use the following simple bounds on
the extreme singular values of a Kronecker sum.

######  Lemma 9.

Let

$L=\bigoplus_{j=1}^{d}M_{j}=M_{1}\otimes I^{\otimes d-1}+I\otimes M_{2}\otimes
I^{\otimes d-2}+\cdots+I^{\otimes d-1}\otimes M_{d},$ (4.60)

where $\\{M_{j}\\}_{j=1}^{d}$ are square matrices. If the largest and smallest
singular values of $M_{j}$ satisfy

$\displaystyle\sigma_{\max}(M_{j})\leq s_{j}^{\max},$ (4.61)
$\displaystyle\sigma_{\min}(M_{j})\geq s_{j}^{\min},$

respectively, then the condition number of $L$ satisfies

$\kappa_{L}\leq\frac{\sum_{j=1}^{d}s_{j}^{\max}}{\sum_{j=1}^{d}s_{j}^{\min}}.$
(4.62)

###### Proof.

We bound the singular values of the matrix exponential $\exp(M_{j})$ by

$\displaystyle\sigma_{\max}(\exp(M_{j}))\leq e^{s_{j}^{\max}},$ (4.63)
$\displaystyle\sigma_{\min}(\exp(M_{j}))\geq e^{s_{j}^{\min}}$

using (4.61). The singular values of the Kronecker product
$\bigotimes_{j=1}^{d}\exp(M_{j})$ are

$\sigma_{k_{1},\ldots,k_{d}}\biggl{(}\bigotimes_{j=1}^{d}\exp(M_{j})\biggr{)}=\prod_{j=1}^{d}\sigma_{k_{j}}(\exp(M_{j}))$
(4.64)

where $\sigma_{k_{j}}(\exp(M_{j}))$ are the singular values of the matrix
$\exp(M_{j})$ for each $j\in[{d}]$, where $k_{j}$ runs from $1$ to the
dimension of $M_{j}$. Using the property of the Kronecker sum that

$\exp(L)=\exp\biggl{(}\bigoplus_{j=1}^{d}M_{j}\biggr{)}=\bigotimes_{j=1}^{d}\exp(M_{j}),$
(4.65)

we bound the singular values of the matrix exponential of (4.55) by

$\displaystyle\sigma_{\max}(\exp(L))\leq e^{\sum_{j=1}^{d}s_{j}^{\max}},$
(4.66) $\displaystyle\sigma_{\min}(\exp(L))\geq
e^{\sum_{j=1}^{d}s_{j}^{\min}}.$

Finally, we bound the singular values of the matrix logarithm of (4.66) by

$\displaystyle\sigma_{\max}(L)\leq\sum_{j=1}^{d}s_{j}^{\max},$ (4.67)
$\displaystyle\sigma_{\min}(L)\geq\sum_{j=1}^{d}s_{j}^{\min}.$

Thus the condition number of $L$ satisfies

$\kappa_{L}\leq\frac{\sum_{j=1}^{d}s_{j}^{\max}}{\sum_{j=1}^{d}s_{j}^{\min}}$
(4.68)

as claimed. ∎

This lemma easily implies a bound on the condition number of the linear system
for Poisson’s equation:

######  Corollary 1.

Consider an instance of the quantum PDE problem as defined in Problem 1 for
Poisson’s equation (2.6) with Dirichlet boundary conditions (4.25). The
condition number of $L_{\mathrm{Poisson}}$ in the linear system (4.34)
satisfies

$\kappa_{L_{\mathrm{Poisson}}}\leq(2n)^{4}.$ (4.69)

###### Proof.

The matrix in (4.34) for Poisson’s equation (2.6) is $L_{\mathrm{Poisson}}$
defined in (4.55). For both the periodic and the non-periodic case, we have

$\displaystyle\sigma_{\max}(\overline{D}^{(2)}_{n})\leq n^{4},$ (4.70)
$\displaystyle\sigma_{\min}(\overline{D}^{(2)}_{n})\geq\frac{1}{16}$

by Lemma 7 and Lemma 8. Let $M_{j}=\overline{D}^{(2)}_{n}$ for $j\in[{d}]$ in
(4.60), and apply Lemma 9 with $s_{j}^{\max}=n^{4}$ and
$s_{j}^{\min}={1}/{16}$ in (4.62). Then the condition number of
$L_{\mathrm{Poisson}}$ is bounded by

$\kappa_{L_{\mathrm{Poisson}}}\leq\frac{\sigma_{\max}(\overline{D}^{(2)}_{n})}{\sigma_{\min}(\overline{D}^{(2)}_{n})}\leq(2n)^{4}$
(4.71)

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.25). The condition number of $L$ in the
linear system (4.34) satisfies

$\kappa_{L}\leq\frac{\|A\|_{\Sigma}}{C\|A\|_{\ast}}(2n)^{4},$ (4.72)

where $\|A\|_{\Sigma}:=\sum_{\|\bm{j}\|_{1}\leq
2}|A_{\bm{j}}|=\sum_{j_{1},j_{2}=1}^{d}|A_{j_{1},j_{2}}|$,
$\|A\|_{\ast}:=\sum_{j=1}^{d}|A_{j,j}|$, and the constant $C$ comes from
(2.8).

###### Proof.

According to (4.49), the matrix in (4.34) is

$L=\sum_{\|\bm{j}\|_{1}=2}A_{\bm{j}}\overline{\bm{D}}^{(\bm{j})}_{\bm{n}}.$
(4.73)

We upper bound the spectral norm of the matrix $L$ by

$\|L\|\leq\sum_{\|\bm{j}\|_{1}=2}|A_{\bm{j}}|\|\overline{\bm{D}}^{(\bm{j})}_{\bm{n}}\|.$
(4.74)

For the matrix $\overline{\bm{D}}^{(\bm{j})}_{\bm{n}}$ defined by (4.50),
Lemma 7 (in the periodic case) and Lemma 8 (in the non-periodic case) give the
inequality

$\|\overline{\bm{D}}^{(\bm{j})}_{\bm{n}}\|\leq n^{4},$ (4.75)

so we have

$\|L\|\leq\sum_{\|\bm{j}\|_{1}=2}|A_{\bm{j}}|n^{4}=\|A\|_{\Sigma}\,n^{4}.$
(4.76)

Next we lower bound $\|L\xi\|$ for any $\|\xi\|=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.55) that can be bounded by Lemma
9, and $L_{2}$ is a sum of tensor products that are easily bounded.
Specifically, we have

$\displaystyle L_{1}$ $\displaystyle=A_{1,1}\overline{D}^{(2)}_{n}\otimes
I^{\otimes d-1}+\cdots+A_{d,d}I^{\otimes d-1}\otimes\overline{D}^{(2)}_{n}$
(4.77) $\displaystyle L_{2}$ $\displaystyle=L-L_{1}.$

The ellipticity condition (2.7) can only hold if the $A_{j,j}$ for $j\in[{d}]$
are either all positive or all negative; we consider $A_{j,j}>0$ without loss
of generality, so

$\|A\|_{\ast}=\sum_{j=1}^{d}|A_{j,j}|=\sum_{j=1}^{d}A_{j,j}.$ (4.78)

then (2.8) can be rewritten as

$1-\sum_{j_{1}=1}^{d}\frac{1}{A_{j_{1},j_{1}}}\sum_{j_{2}\in[{d}]\backslash\\{j_{1}\\}}|A_{j_{1},j_{2}}|=C>0,$
(4.79)

where $0<C\leq 1$.

We now upper bound $\|L_{2}L_{1}^{-1}\|$ by bounding
$\|\bm{D}^{(\bm{j})}_{\bm{n}}L_{1}^{-1}\|$ for each
$\bm{j}=(j_{1},\ldots,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}\in[{d}]$, $r_{1}\neq r_{2}$, and $j_{r}=0$ for
$r\in[{d}]\backslash\\{r_{1},r_{2}\\}$. We denote

$L^{(\bm{j})}:=I^{\otimes r_{1}-1}\otimes D^{2}_{n}\otimes I^{\otimes
d-r_{1}}+I^{\otimes r_{2}-1}\otimes D^{2}_{n}\otimes I^{\otimes d-r_{2}}.$
(4.80)

We first upper bound $\|\bm{D}^{(\bm{j})}_{\bm{n}}\|$ by
$\frac{1}{2}\|L^{(\bm{j})}\|$. Notice the matrices
$\bm{D}^{(\bm{j})}_{\bm{n}}$ and $L^{(\bm{j})}$ share the same singular
vectors. For $k\in[{n+1}]_{0}$, we let $v_{k}$ and $\lambda_{k}$ denote the
right singular vectors and corresponding singular values of $D_{n}$,
respectively. Then the right singular vectors of $\bm{D}^{(\bm{j})}_{\bm{n}}$
and $L^{(\bm{j})}$ are $\bm{v}_{\bm{k}}:=\bigotimes_{j=1}^{d}v_{k_{j}}$, where
$\bm{k}=(k_{1},\ldots,k_{d})$ with $k_{j}\in[{n+1}]_{0}$ for $j\in[{d}]$. For
any vector $v=\sum_{\|\bm{k}\|_{\infty}\leq n}\alpha_{\bm{k}}\bm{v}_{\bm{k}}$,
we have

$\|\bm{D}^{(\bm{j})}_{\bm{n}}v\|^{2}=\sum_{\|\bm{k}\|_{\infty}\leq
n}|\alpha_{\bm{k}}|^{2}\|\bm{D}^{(\bm{j})}_{\bm{n}}\bm{v}_{\bm{k}}\|^{2}=\sum_{\|\bm{k}\|_{\infty}\leq
n}|\alpha_{\bm{k}}|^{2}(\lambda_{k_{j_{r_{1}}}}\lambda_{k_{j_{r_{2}}}})^{2},$
(4.81) $\|L^{(\bm{j})}v\|^{2}=\sum_{\|\bm{k}\|_{\infty}\leq
n}|\alpha_{\bm{k}}|^{2}\|L^{(\bm{j})}\bm{v}_{\bm{k}}\|^{2}=\sum_{\|\bm{k}\|_{\infty}\leq
n}|\alpha_{\bm{k}}|^{2}(\lambda_{k_{j_{r_{1}}}}^{2}+\lambda_{k_{j_{r_{2}}}}^{2})^{2},$
(4.82)

which implies
$\|\bm{D}^{(\bm{j})}_{\bm{n}}v\|\leq\frac{1}{2}\|L^{(\bm{j})}v\|$ by the AM-GM
inequality. Since this holds for any vector $v$, we have

$\|\bm{D}^{(\bm{j})}_{\bm{n}}L_{1}^{-1}\|\leq\frac{1}{2}\|L^{(\bm{j})}L_{1}^{-1}\|.$
(4.83)

Next we upper bound $\|D^{2}_{n}\|$ by $\|\overline{D}^{(2)}_{n}\|$. For any
vector $u=[u_{0},\ldots,u_{n}]^{T}$, define two vectors
$w=[w_{0},\ldots,w_{n}]^{T}$ and
$\overline{w}=[\overline{w}_{0},\ldots,\overline{w}_{n}]^{T}$ such that

$D^{2}_{n}[u_{0},\ldots,u_{n}]^{T}=[w_{0},\ldots,w_{n}]^{T}$ (4.84)

and

$\overline{D}^{2}_{n}[u_{0},\ldots,u_{n}]^{T}=[\overline{w}_{0},\ldots,\overline{w}_{n}]^{T}.$
(4.85)

Notice that $w_{\mathopen{}\mathclose{{}\left\lfloor n/2}\right\rfloor}=0$ and
$w_{k}=\overline{w}_{k}$ for
$k\in[{n+1}]_{0}\backslash\\{\mathopen{}\mathclose{{}\left\lfloor
n/2}\right\rfloor\\}$ for periodic conditions, and $w_{n-1}=w_{n}=0$ and
$w_{k}=\overline{w}_{k}$ for $k\in[{n+1}]_{0}\backslash\\{n-1,n\\}$ for non-
periodic conditions. Thus, for any vector $v$,

$\|D^{2}_{n}v\|^{2}=\|w\|^{2}=\sum_{k=0}^{n}w_{k}^{2}\leq\sum_{k=0}^{n}\overline{w}_{k}^{2}=\|\overline{w}\|^{2}=\|\overline{D}^{(2)}_{n}v\|^{2}.$
(4.86)

Therefore,

$\|L^{(\bm{j})}L_{1}^{-1}\|\leq\sum_{s=1}^{2}\|I^{\otimes r_{s}-1}\otimes
D^{2}_{n}\otimes I^{\otimes d-r_{s}}L_{1}^{-1}\|\leq\sum_{s=1}^{2}\|I^{\otimes
r_{s}-1}\otimes\overline{D}^{(2)}_{n}\otimes I^{\otimes d-r_{s}}L_{1}^{-1}\|.$
(4.87)

We also have

$\|\bm{D}^{(\bm{j})}_{\bm{n}}L_{1}^{-1}\|\leq\frac{1}{2}\sum_{s=1}^{2}\|I^{\otimes
r_{s}-1}\otimes\overline{D}^{(2)}_{n}\otimes I^{\otimes d-r_{s}}L_{1}^{-1}\|.$
(4.88)

We can rewrite $I^{\otimes r_{s}-1}\otimes\overline{D}^{(2)}_{n}\otimes
I^{\otimes d-r_{s}}L_{1}^{-1}$ in the form

$I^{\otimes r_{s}-1}\otimes\overline{D}^{(2)}_{n}\otimes I^{\otimes
d-r_{s}}\mathopen{}\mathclose{{}\left(\sum_{h=1}^{d}A_{h,h}I^{\otimes
r_{h}-1}\otimes\overline{D}^{(2)}_{n}\otimes I^{\otimes
d-r_{h}}}\right)^{-1}.$ (4.89)

The matrices $I^{\otimes r_{h}-1}\otimes\overline{D}^{(2)}_{n}\otimes
I^{\otimes d-r_{h}}$ share the same singular values and singular vectors, so

$\displaystyle\|I^{\otimes r_{s}-1}\otimes\overline{D}^{(2)}_{n}\otimes
I^{\otimes
d-r_{s}}L_{1}^{-1}\|=\max_{\overline{\lambda}_{k_{r}}}\frac{\overline{\lambda}_{k_{r_{s}}}}{\sum_{h=1}^{d}A_{h,h}\overline{\lambda}_{k_{h}}}<\frac{1}{A_{{r_{s}},{r_{s}}}},$
(4.90)

where $\overline{\lambda}_{k_{h}}$ are singular values of $I^{\otimes
r_{h}-1}\otimes\overline{D}^{(2)}_{n}\otimes I^{\otimes d-r_{h}}$ for
$k_{h}\in[{n}]_{0}$, $h\in[{d}]$. This implies

$\|\bm{D}^{(\bm{j})}_{\bm{n}}L_{1}^{-1}\|\leq\frac{1}{2}(\frac{1}{A_{{r_{1}},{r_{1}}}}+\frac{1}{A_{{r_{2}},{r_{2}}}}).$
(4.91)

Using (4.79), considering each instance of $\bm{D}^{(\bm{j})}_{\bm{n}}$ in
$L_{2}$, we have

$\|L_{2}L_{1}^{-1}\|\leq\sum_{j_{1}\neq
j_{2}}|A_{j_{1},j_{2}}|\|\bm{D}^{(\bm{j})}_{\bm{n}}L_{1}^{-1}\|\leq\sum_{j_{1}=1}^{d}\frac{1}{A_{j_{1},j_{1}}}\sum_{j_{2}\in[{d}]\backslash\\{j_{1}\\}}|A_{j_{1},j_{2}}|\leq
1-C.$ (4.92)

Since $L$ and $L_{1}$ are invertible, $\|L_{1}^{-1}L_{2}\|\leq 1-C<1$, and by
Lemma 9 applied to $\|L_{1}^{-1}\|$, we have

$\|L^{-1}\|=\|(L_{1}+L_{2})^{-1}\|\leq\|(I+L_{2}L_{1}^{-1})^{-1}\|\|L_{1}^{-1}\|\leq\frac{\|L_{1}^{-1}\|}{1-\|L_{2}L_{1}^{-1}\|}\leq\frac{1/\frac{1}{16}\|A\|_{\ast}}{C}=\frac{16}{C\|A\|_{\ast}}.$
(4.93)

Thus we have

$\kappa_{L}=\|L\|\|L^{-1}\|\leq\frac{\|A\|_{\Sigma}}{C\|A\|_{\ast}}(2n)^{4}$
(4.94)

as claimed. ∎

###  4.4 State preparation

We now describe a state preparation procedure for the vector $\vec{f}+\vec{g}$
in the linear system (4.34).

######  Lemma 11.

Let $O_{f}$ be a unitary oracle that maps $|0\rangle|0\rangle$ to a state
proportional to $|0\rangle|f\rangle$, and $|\phi\rangle|0\rangle$ to
$|\phi\rangle|0\rangle$ for any $|\phi\rangle$ orthogonal to $|0\rangle$; let
$O_{x}$ be a unitary oracle that maps $|0\rangle|0\rangle$ to
$|0\rangle|0\rangle$, $|j\rangle|0\rangle$ to a state proportional to
$|j\rangle|\gamma^{j+}\rangle$ for $j\in[{d}]$, and $|j+d\rangle|0\rangle$ to
a state proportional to $|j+d\rangle|\gamma^{j-}\rangle$ for $j\in[{d}]$.
Suppose $\|{|f\rangle}\|,\|{|\gamma^{j+}\rangle}\|,\|{|\gamma^{j-}\rangle}\|$
and $A_{j,j}$ for $j\in[{d}]$ are known. Define the parameter

$q:=\sqrt{\frac{\sum_{\|\bm{k}\|_{\infty}\leq
n}\sum_{j=1}^{d}[\hat{f}^{2}_{\bm{k}}+(A_{j,j}\hat{\gamma}^{j+}_{\bm{k}})^{2}+(A_{j,j}\hat{\gamma}^{j-}_{\bm{k}})^{2}]}{\sum_{\|\bm{k}\|_{\infty}\leq
n}\sum_{j=1}^{d}|\hat{f}_{\bm{k}}+A_{j,j}\hat{\gamma}^{j+}_{\bm{k}}+A_{j,j}\hat{\gamma}^{j-}_{\bm{k}}|^{2}}}.$
(4.95)

Then the normalized quantum state

$|B\rangle\propto\sum_{\|\bm{k}\|_{\infty}\leq
n}\sum_{j=1}^{d}(\hat{f}_{\bm{k}}+A_{j,j}\hat{\gamma}^{j+}_{\bm{k}}+A_{j,j}\hat{\gamma}^{j-}_{\bm{k}})|k_{1}\rangle\ldots|k_{d}\rangle,$
(4.96)

with coefficients defined as in (4.32) and (4.33), can be prepared with gate
and query complexity $O(qd^{2}\log n\log\log n)$.

###### Proof.

Starting from the initial state $|0\rangle|0\rangle$, we first perform a
unitary transformation $U$ satisfying

$\displaystyle U|0\rangle=$
$\displaystyle\frac{\||f\rangle\|}{\sqrt{\||f\rangle\|^{2}+\sum_{j=1}^{d}\bigl{(}A_{j,j}^{2}\||\gamma^{j+}\rangle\|^{2}+A_{j,j}^{2}\||\gamma^{j-}\rangle\|^{2}\bigr{)}}}|0\rangle$
(4.97)
$\displaystyle+\sum_{j=1}^{d}\frac{A_{j,j}\||\gamma^{j+}\rangle\|}{\sqrt{\||f\rangle\|^{2}+\sum_{j=1}^{d}\bigl{(}A_{j,j}^{2}\||\gamma^{j+}\rangle\|^{2}+A_{j,j}^{2}\||\gamma^{j-}\rangle\|^{2}\bigr{)}}}|j\rangle$
$\displaystyle+\sum_{j=1}^{d}\frac{A_{j,j}\||\gamma^{j-}\rangle\|}{\sqrt{\||f\rangle\|^{2}+\sum_{j=1}^{d}\bigl{(}A_{j,j}^{2}\||\gamma^{j+}\rangle\|^{2}+A_{j,j}^{2}\||\gamma^{j-}\rangle\|^{2}\bigr{)}}}|j+d\rangle$

on the first register to obtain

$\frac{\||f\rangle\||0\rangle+A_{1,1}\||\gamma^{1+}\rangle\||1\rangle+\cdots+A_{d,d}\||\gamma^{d-}\rangle\||2d\rangle}{\sqrt{\||f\rangle\|^{2}+\sum_{j=1}^{d}\bigl{(}A_{j,j}^{2}\||\gamma^{j+}\rangle\|^{2}+A_{j,j}^{2}\||\gamma^{j-}\rangle\|^{2}\bigr{)}}}|0\rangle.$
(4.98)

This can be done in time $O(2d+1)$ by standard techniques [30]. Then we apply
$O_{x}$ and $O_{f}$ to obtain

$\displaystyle|0\rangle|f\rangle+A_{1,1}|1\rangle|\gamma^{1+}\rangle+\cdots+A_{d,d}|2d\rangle|\gamma^{d-}\rangle,$
(4.99)
$\displaystyle\quad\propto\sum_{\|\bm{k}\|_{\infty},\|\bm{l}\|_{\infty}\leq
n}\phi_{\bm{k}}(\bm{\chi}_{\bm{l}})(\hat{f}_{\bm{k}}|0\rangle+A_{1,1}\hat{\gamma}^{1+}_{\bm{k}}|1\rangle+\cdots+A_{d,d}\hat{\gamma}^{d-}_{\bm{k}}|2d\rangle)|l_{1}\rangle\ldots|l_{d}\rangle,$

according to (4.30) and (4.31). We then perform the $d$-dimensional inverse
QSFT (for periodic boundary conditions) or inverse QCT (for non-periodic
boundary conditions) on the last $d$ registers, obtaining

$\sum_{\|\bm{k}\|_{\infty}\leq
n}(\hat{f}_{\bm{k}}|0\rangle+A_{1,1}\hat{\gamma}^{1+}_{\bm{k}}|1\rangle+\cdots+A_{d,d}\hat{\gamma}^{d-}_{\bm{k}}|2d\rangle)|k_{1}\rangle\ldots|k_{d}\rangle.$
(4.100)

Finally, observe that if we measure the first register in a basis containing
the uniform superposition $|0\rangle+|1\rangle+\cdots+|2d\rangle$ (say, the
Fourier basis) and obtain the outcome corresponding to the uniform
superposition, we produce the state

$\sum_{\|\bm{k}\|_{\infty}\leq
n}\sum_{j=1}^{d}(\hat{f}_{\bm{k}}+A_{j,j}\hat{\gamma}^{j+}_{\bm{k}}+A_{j,j}\hat{\gamma}^{j-}_{\bm{k}})|k_{1}\rangle\ldots|k_{d}\rangle.$
(4.101)

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\rangle$, then we may be able to avoid the factor of $q$ in the complexity
of the overall algorithm.

###  4.5 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.

Consider an instance of the quantum PDE problem as defined in Problem 1 with
Dirichlet boundary conditions (4.25). Then there exists a quantum algorithm
that produces a state in the form of (4.26) whose amplitudes are proportional
to $u(\bm{x})$ on a set of interpolation nodes $\bm{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.4)), where $u(\bm{x})/\|u(\bm{x})\|$ is $\epsilon$-close to
$\hat{u}(\bm{x})/\|\hat{u}(\bm{x})\|$ in $l_{2}$ norm for all nodes $\bm{x}$,
succeeding with probability $\Omega(1)$, with a flag indicating success, using

$\biggl{(}\frac{d\|A\|_{\Sigma}}{C\|A\|_{\ast}}+qd^{2}\biggr{)}\operatorname{poly}(\log(g^{\prime}/g\epsilon))$
(4.102)

queries to oracles as defined in Section 4.4. Here
$\|A\|_{\Sigma}:=\sum_{\|\bm{j}\|_{1}\leq h}\|A_{\bm{j}}\|$,
$\|A\|_{\ast}:=\sum_{j=1}^{d}|A_{j,j}|$, $C$ is defined in (2.8), and

$\displaystyle g=\min_{\bm{x}}\|\hat{u}(\bm{x})\|,\qquad
g^{\prime}:=\max_{\bm{x}}\max_{n\in{\mathbb{N}}}\|\hat{u}^{(n+1)}(\bm{x})\|,\qquad$
(4.103) $\displaystyle q=\sqrt{\frac{\sum_{\|\bm{k}\|_{\infty}\leq
n}\sum_{j=1}^{d}\hat{f}^{2}_{\bm{k}}+(A_{j,j}\hat{\gamma}^{j+}_{\bm{k}})^{2}+(A_{j,j}\hat{\gamma}^{j-}_{\bm{k}})^{2}}{\sum_{\|\bm{k}\|_{\infty}\leq
n}\sum_{j=1}^{d}(\hat{f}_{\bm{k}}+A_{j,j}\hat{\gamma}^{j+}_{\bm{k}}+A_{j,j}\hat{\gamma}^{j-}_{\bm{k}})^{2}}}.$
(4.104)

The gate complexity is larger than the query complexity by a factor of
$\operatorname{poly}(\log(d\|A\|_{\Sigma}/\epsilon))$.

###### Proof.

We analyze the complexity of the algorithm presented in Section 4.2.

First we choose

$n:=\biggl{\lfloor}\frac{\log(\Omega)}{\log(\log(\Omega))}\biggr{\rfloor},$
(4.105)

where

$\Omega=\frac{g^{\prime}(1+\epsilon)}{g\epsilon}.$ (4.106)

By Eq. (1.8.28) of [27], this choice guarantees

$\|\hat{u}(\bm{x})-u(\bm{x})\|\leq\max_{t\in[0,T]}\|\hat{u}^{(n+1)}(\bm{x})\|\frac{e^{n}}{(2n)^{n}}\leq\frac{g^{\prime}}{\Omega}=\frac{g\epsilon}{1+\epsilon}=:\delta.$
(4.107)

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

$\biggl{\|}\frac{\hat{u}(\bm{x})}{\|\hat{u}(\bm{x})\|}-\frac{u(\bm{x})}{\|u(\bm{x})\|}\biggr{\|}\leq\frac{\delta}{\min\\{\|\hat{u}(\bm{x})\|,\|u(\bm{x})\|\\}}\leq\frac{\delta}{g-\delta}=\epsilon,$
(4.108)

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

As described in Section 4.2, the algorithm uses the high-precision QLSA from
reference [7] 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\vec{x}=\vec{f}+\vec{g}$ in (4.34), the matrix $L$ is
an $(n+1)^{d}\times(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
$\frac{\|A\|_{\Sigma}}{C\|A\|_{\ast}}(2n)^{4}$. Consequently, by Theorem 5 of
[7], the QLSA produces a state proportional to $\vec{x}$ with
$O(\frac{d\|A\|_{\Sigma}}{C\|A\|_{\ast}}(2n)^{5})$ queries to the oracles, and
its gate complexity is larger by a factor of
$\operatorname{poly}(\log(d\|A\|_{\Sigma}\,n))$. Using the value of $n$
specified in (4.105), the overall query complexity of our algorithm is

$\biggl{(}\frac{d\|A\|_{\Sigma}}{C\|A\|_{\ast}}+qd^{2}\biggr{)}\operatorname{poly}(\log(g^{\prime}/g\epsilon)),$
(4.109)

and the gate complexity is

$\biggl{(}\frac{d\|A\|_{\Sigma}}{C\|A\|_{\ast}}\operatorname{poly}(\log(d\|A\|_{\Sigma}/\epsilon))+qd^{2}\biggr{)}\operatorname{poly}(\log(g^{\prime}/g\epsilon))$
(4.110)

which is larger by a factor of
$\operatorname{poly}(\log(d\|A\|_{\Sigma}/\epsilon))$, 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\|_{\Sigma}=\|A\|_{\ast}=d$ and $C=1$. Under homogenous boundary
conditions, the complexity of state preparation can be reduced to
$d\operatorname{poly}(\log(g^{\prime}/g\epsilon))$, since we can remove $2d$
applications of the QSFT or QCT for preparing a state depending on the
boundary conditions, and since $\gamma=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\operatorname{poly}(\log(g^{\prime}/g\epsilon));$ (4.111)

again the gate complexity is larger by a factor of
$\operatorname{poly}(\log(d\|A\|_{\Sigma}/\epsilon))$.

##  5 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 second-order elliptic equations. Whereas
previous algorithms scaled as $\operatorname{poly}(d,1/\epsilon)$, our
algorithms scale as $\operatorname{poly}(d,\log(1/\epsilon))$.

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 time-dependent PDEs? These
cases are more challenging since the quantum Fourier transform cannot be
directly applied to ensure sparsity.

Second, the complexity scales logarithmically with high-order derivatives (of
the inhomogeneity or solution) for both the adative FDM and spectral method.
However, for time-dependent Hamiltonian simulation, the query complexity only
depends on the first-order derivatives of the Hamiltonian [22, 4]. Can we
develop quantum algorithms for PDEs with query complexity independent of high-
order derivatives?

Third, can we develop quantum algorithms for stochastic PDEs or for 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 [11] 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?

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?

## Acknowledgments

JPL thanks Jacob Bedrossian, Dominic Berry, Stephen Jordan, Ashley Montanaro,
and Konstantina Trivisa for valuable discussions. We thank anonymous referees
for their comments on an earlier draft.

The authors acknowledge support from National Science Foundation grant
CCF-1813814 and from the U.S. Department of Energy, Office of Science, Office
of Advanced Scientific Computing Research, Quantum Algorithms Teams and
Accelerated Research in Quantum Computing programs.

## References

  * [1] Andris Ambainis, _Variable time amplitude amplification and quantum algorithms for linear algebra problems_ , 29th Symposium on Theoretical Aspects of Computer Science, vol. 14, pp. 636–647, LIPIcs, 2012, arXiv:1010.4458.
  * [2] Ivo Babuška and Manil Suri, _The $h-p$ version of the finite element method with quasiuniform meshes_, Mathematical Modelling and Numerical Analysis 21 (1987), no. 2, 199–238.
  * [3] Dominic W. Berry, Andrew M. Childs, Richard Cleve, Robin Kothari, and Rolando D. Somma, _Simulating Hamiltonian dynamics with a truncated Taylor series_ , Physical Review Letters 114 (2015), no. 9, 090502, arXiv:1412.4687.
  * [4] , _Exponential improvement in precision for simulating sparse Hamiltonians_ , Forum of Mathematics, Sigma 5 (2017), e8, arXiv:1312.1414.
  * [5] Hans-Joachim Bungartz and Michael Griebel, _Sparse grids_ , Acta Numerica 13 (2004), 147–269.
  * [6] Yudong Cao, Anargyros Papageorgiou, Iasonas Petras, Joseph Traub, and Sabre Kais, _Quantum algorithm and circuit design solving the Poisson equation_ , New Journal of Physics 15 (2013), no. 1, 013021, arXiv:1207.2485.
  * [7] Andrew M. Childs, Robin Kothari, and Rolando D. Somma, _Quantum algorithm for systems of linear equations with exponentially improved dependence on precision_ , SIAM Journal on Computing 46 (2017), no. 6, 1920–1950, arXiv:1511.02306.
  * [8] Andrew M. Childs and Jin-Peng Liu, _Quantum spectral methods for differential equations_ , arXiv:1901.00961 (2019).
  * [9] B. David Clader, Bryan C. Jacobs, and Chad R. Sprouse, _Preconditioned quantum linear system algorithm_ , Physical Review Letters 110 (2013), no. 25, 250504, arXiv:1301.2340.
  * [10] Daniel T. Colbert and William H. Miller, _A novel discrete variable representation for quantum mechanical reactive scattering via the S-matrix Kohn method_ , The Journal of Chemical Physics 96 (1992), no. 3, 1982–1991.
  * [11] Pedro C. S. Costa, Stephen Jordan, and Aaron Ostrander, _Quantum algorithm for simulating the wave equation_ , Physical Review A 99 (2019), no. 1, 012323, arXiv:1711.05394.
  * [12] Lawrence C. Evans, _Partial differential equations (2nd ed.)_ , American Mathematical Society, 2010.
  * [13] François Fillion-Gourdeau and Emmanuel Lorin, _Simple digital quantum algorithm for symmetric first-order linear hyperbolic systems_ , Numerical Algorithms 82 (2019), 1009–1045, arXiv:1705.09361.
  * [14] Călin Ioan Gheorghiu, _Spectral methods for differential problems_ , Casa Cărţii de Ştiinţă Cluj-Napoca, 2007.
  * [15] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd, _Quantum algorithm for linear systems of equations_ , Physical Review Letters 103 (2009), no. 15, 150502, arXiv:0811.3171.
  * [16] Roger A. Horn and Charles R. Johnson, _Matrix analysis_ , Cambridge University Press, 2012.
  * [17] Ian D. Kivlichan, Nathan Wiebe, Ryan Babbush, and Alán Aspuru-Guzik, _Bounding the costs of quantum simulation of many-body physics in real space_ , Journal of Physics A: Mathematical and Theoretical 50 (2017), no. 30, 305301, arXiv:1608.05696.
  * [18] Andreas Klappenecker and Martin Rotteler, _Discrete cosine transforms on quantum computers_ , Proceedings of the 2nd International Symposium on Image and Signal Processing and Analysis, pp. 464–468, 2001, arXiv:quant-ph/0111038.
  * [19] Benjamin P. Lanyon, James D. Whitfield, Geoff G. Gillett, Michael E. Goggin, Marcelo P. Almeida, Ivan Kassal, Jacob D. Biamonte, Masoud Mohseni, Ben J. Powell, Marco Barbieri, Alán Aspuru-Guzik, and Andrew G. White, _Towards quantum chemistry on a quantum computer_ , Nature Chemistry 2 (2010), no. 2, 106, arXiv:0905.0887.
  * [20] Jianping Li, _General explicit difference formulas for numerical differentiation_ , Journal of Computational and Applied Mathematics 183 (2005), no. 1, 29–52.
  * [21] Ashley Montanaro and Sam Pallister, _Quantum algorithms and the finite element method_ , Physical Review A 93 (2016), no. 3, 032324, arXiv:1512.05903.
  * [22] David Poulin, Angie Qarry, Rolando D. Somma, and Frank Verstraete, _Quantum simulation of time-dependent Hamiltonians and the convenient illusion of Hilbert space_ , Physical Review Letters 106 (2011), no. 17, 170501, arXiv:1102.1360.
  * [23] Markus Püschel, Martin Rötteler, and Thomas Beth, _Fast quantum Fourier transforms for a class of non-abelian groups_ , International Symposium on Applied Algebra, Algebraic Algorithms, and Error-Correcting Codes, pp. 148–159, 1999, arXiv:quant-ph/9807064.
  * [24] Markus Reiher, Nathan Wiebe, Krysta M. Svore, Dave Wecker, and Matthias Troyer, _Elucidating reaction mechanisms on quantum computers_ , Proceedings of the National Academy of Sciences 114 (2017), no. 29, 7555–7560, arXiv:1605.03590.
  * [25] Martin Rötteler, Markus Püschel, and Thomas Beth, _Fast signal transforms for quantum computers_ , Proceedings of the Workshop on Physics and Computer Science, pp. 31–43, 1999.
  * [26] Norbert Schuch and Jens Siewert, _Programmable networks for quantum algorithms_ , Physical Review Letters 91 (2003), no. 2, 027902, arXiv:quant-ph/0303063.
  * [27] Jie Shen, Tao Tang, and Li-Lian Wang, _Spectral methods: algorithms, analysis and applications_ , vol. 41, Springer Science & Business Media, 2011\.
  * [28] Jie Shen and Haijun Yu, _Efficient spectral sparse grid methods and applications to high-dimensional elliptic problems_ , SIAM Journal on Scientific Computing 32 (2010), no. 6, 3228–3250.
  * [29] Jie Shen and Haijun Yu, _Efficient spectral sparse grid methods and applications to high-dimensional elliptic equations II. Unbounded domains_ , SIAM Journal on Scientific Computing 34 (2012), no. 2, A1141–A1164.
  * [30] Vivek V. Shende, Stephen S. Bullock, and Igor L. Markov, _Synthesis of quantum-logic circuits_ , IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems 25 (2006), no. 6, 1000–1010, arXiv:quant-ph/0406176.
  * [31] Sergei Abramovich Smolyak, _Quadrature and interpolation formulas for tensor products of certain classes of functions_ , Doklady Akademii Nauk, vol. 148, pp. 1042–1045, 1963.
  * [32] Daniel Spielman, _Rings, paths, and Cayley graphs (course notes)_ , 2014, http://www.cs.yale.edu/homes/spielman/561/lect05-15.pdf.
  * [33] Shengbin Wang, Zhimin Wang, Wendong Li, Lixin Fan, Zhiqiang Wei, and Yongjian Gu, _Quantum fast Poisson solver: the algorithm and modular circuit design_ , arXiv preprint arXiv:1910.09756 (2019).
  * [34] Jonathan Welch, Daniel Greenbaum, Sarah Mostame, and Alán Aspuru-Guzik, _Efficient quantum circuits for diagonal unitaries without ancillas_ , New Journal of Physics 16 (2014), no. 3, 033040, arXiv:1306.3991.
  * [35] Stephen Wiesner, _Simulations of many-body quantum systems by a quantum computer_ , arXiv:quant-ph/9603028 (1996).
  * [36] Christof Zalka, _Efficient simulation of quantum systems by quantum computers_ , Fortschritte der Physik 46 (1998), no. 6-8, 877–879, arXiv:quant-ph/9603026.
  * [37] Christoph Zenger, _Sparse grids_ , (1991), https://www5.in.tum.de/pub/zenger91sg.pdf.

##  Appendix 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

$\displaystyle\Delta u(x_{1},x_{2})$
$\displaystyle=\mathopen{}\mathclose{{}\left(\frac{\partial^{2}}{\partial
x_{1}^{2}}+\frac{\partial^{2}}{\partial
x_{2}^{2}}}\right)u(x_{1},x_{2})=f(x_{1},x_{2}),$ $\displaystyle(x_{1},x_{2})$
$\displaystyle\in\Omega=[-1,1]^{2}.$ (A.1)

We consider two kinds of boundary value problems, as follows.

  * •

_Periodic boundary conditions:_

$\displaystyle u(x_{1},x_{2})$
$\displaystyle=u(x_{1}+2v,x_{2})=u(x_{1},x_{2}+2v),$
$\displaystyle(x_{1},x_{2})$
$\displaystyle\in\Omega=[-1,1]^{2},v\in{\mathbb{Z}}$ (A.2) $\displaystyle
u(0,0)$ $\displaystyle=\gamma.$

  * •

_Non-periodic boundary conditions:_

$\displaystyle u(x_{1},1)$ $\displaystyle=\gamma_{N}(x_{1}),$ $\displaystyle
u(1,x_{2})$ $\displaystyle=\gamma_{E}(x_{2}),$ (A.3) $\displaystyle
u(x_{1},-1)$ $\displaystyle=\gamma_{S}(x_{1}),$ $\displaystyle u(-1,x_{2})$
$\displaystyle=\gamma_{W}(x_{2}).$

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

$u(x_{1},x_{2})=\sum_{k_{1}=0}^{2}\sum_{k_{2}=0}^{2}c_{k_{1},k_{2}}e^{i(k_{1}-1)\pi
x_{1}}e^{i(k_{2}-1)\pi x_{2}}.$ (A.4)

We are given an oracle for preparing the state

$\sum_{l_{1}=0}^{2}\sum_{l_{2}=0}^{2}f\mathopen{}\mathclose{{}\left(\frac{2l_{1}}{3}-1,\frac{2l_{2}}{3}-1}\right)|l_{1}\rangle|l_{2}\rangle$
(A.5)

that interpolates the uniform grid nodes (4.4). We first perform a multi-
dimensional inverse QSFT on (A.5) to obtain

$\sum_{k_{1}=0}^{2}\sum_{k_{2}=0}^{2}f_{k_{1},k_{2}}|k_{1}\rangle|k_{2}\rangle,$
(A.6)

where $b_{k_{1},k_{2}}$ are defined by (4.32). Then we apply the quantum
linear system algorithm of [7] to solve the linear system

$L_{p}|X\rangle=|B\rangle,$ (A.7)

where the solution is

$|X\rangle=\sum_{k_{1}=0}^{2}\sum_{k_{2}=0}^{2}c_{k_{1},k_{2}}|k_{1}\rangle|k_{2}\rangle.$
(A.8)

According to (4.49),the discretized linear system from (A.1) is

$D^{2}_{n}\otimes I+I\otimes D^{2}_{n},$ (A.9)

where the Fourier difference matrix $D_{n}$ is defined by (4.44) with $n=2$,
namely

$D_{2}=\begin{pmatrix}-i\pi&0&0\\\ 0&0&0\\\ 0&0&i\pi\\\ \end{pmatrix},$ (A.10)

so that

$D_{2}^{2}=\begin{pmatrix}-\pi^{2}&0&0\\\ 0&0&0\\\ 0&0&-\pi^{2}\\\
\end{pmatrix}.$ (A.11)

Therefore, the matrix (A.9) is

$D^{2}_{2}\otimes I+I\otimes D^{2}_{2}=\begin{pmatrix}-2\pi^{2}&&&&&&&&\\\
&-\pi^{2}&&&&&&&\\\ &&-2\pi^{2}&&&&&&\\\ &&&-\pi^{2}&&&&&\\\ &&&&0&&&&\\\
&&&&&-\pi^{2}&&&\\\ &&&&&&-2\pi^{2}&&\\\ &&&&&&&-\pi^{2}&\\\
&&&&&&&&-2\pi^{2}\\\ \end{pmatrix}.$ (A.12)

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:

$\begin{pmatrix}-2\pi^{2}&&&&&&&&\\\ &-\pi^{2}&&&&&&&\\\ &&-2\pi^{2}&&&&&&\\\
&&&-\pi^{2}&&&&&\\\ &&&1&1&1&&&\\\ &&&&&-\pi^{2}&&&\\\ &&&&&&-2\pi^{2}&&\\\
&&&&&&&-\pi^{2}&\\\ &&&&&&&&-2\pi^{2}\\\
\end{pmatrix}\begin{pmatrix}c_{0,0}\\\ c_{0,1}\\\ c_{0,2}\\\ c_{1,0}\\\
c_{1,1}\\\ c_{1,2}\\\ c_{2,0}\\\ c_{2,1}\\\ c_{2,2}\\\
\end{pmatrix}=\begin{pmatrix}f_{0,0}\\\ f_{0,1}\\\ f_{0,2}\\\ f_{1,0}\\\
\gamma\\\ f_{1,2}\\\ f_{2,0}\\\ f_{2,1}\\\ f_{2,2}\\\ \end{pmatrix},$ (A.13)

where the additional linear equation comes from
$u(0,0)=\sum_{k_{1}=0}^{2}\sum_{k_{2}=0}^{2}c_{k_{1},k_{2}}=\gamma$. In some
problems, we might be directly given the value of $c_{1,1}$, say,
$c_{1,1}=\gamma$. Then the linear system would be

$\begin{pmatrix}-2\pi^{2}&&&&&&&&\\\ &-\pi^{2}&&&&&&&\\\ &&-2\pi^{2}&&&&&&\\\
&&&-\pi^{2}&&&&&\\\ &&&&1&&&&\\\ &&&&&-\pi^{2}&&&\\\ &&&&&&-2\pi^{2}&&\\\
&&&&&&&-\pi^{2}&\\\ &&&&&&&&-2\pi^{2}\\\
\end{pmatrix}\begin{pmatrix}c_{0,0}\\\ c_{0,1}\\\ c_{0,2}\\\ c_{1,0}\\\
c_{1,1}\\\ c_{1,2}\\\ c_{2,0}\\\ c_{2,1}\\\ c_{2,2}\\\
\end{pmatrix}=\begin{pmatrix}f_{0,0}\\\ f_{0,1}\\\ f_{0,2}\\\ f_{1,0}\\\
\gamma\\\ f_{1,2}\\\ f_{2,0}\\\ f_{2,1}\\\ f_{2,2}\\\ \end{pmatrix}.$ (A.14)

We now present the quantum Chebyshev spectral method to solve (A.1) with non-
periodic 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

$u(x_{1},x_{2})=\sum_{k_{1}=0}^{3}\sum_{k_{2}=0}^{3}c_{k_{1},k_{2}}T_{k_{1}}(x_{1})T_{k_{2}}(x_{2}).$
(A.15)

We are given an oracle for preparing the state

$\sum_{l_{1}=0}^{3}\sum_{l_{2}=0}^{3}f\mathopen{}\mathclose{{}\left(\cos\frac{\pi
l_{1}}{3},\cos\frac{\pi l_{2}}{3}}\right)|l_{1}\rangle|l_{2}\rangle$ (A.16)

that interpolates the Chebyshev-Gauss-Lobatto quadrature nodes specified in
(4.4). 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.32). Then we apply the
quantum linear system algorithm of [7] to solve the linear system (A.7) with
the solution (A.8). The discretized linear system from (A.1) is (A.9), where
the Chebyshev difference matrix $D_{n}$ is defined by (4.45) with $n=3$,
namely

$D_{3}=\begin{pmatrix}0&1&0&3\\\ 0&0&4&0\\\ 0&0&0&6\\\ 0&0&0&0\\\
\end{pmatrix},$ (A.17)

and

$D_{3}^{2}=\begin{pmatrix}0&0&4&0\\\ 0&0&0&24\\\ 0&0&0&0\\\ 0&0&0&0\\\
\end{pmatrix}.$ (A.18)

Notice that the rank of $D_{n}^{2}$ is $n-1$, which implies the second
derivative for $d=1$ can be represented as

$u^{\prime\prime}(x)=c^{\prime\prime}_{0}T_{0}(x)+c^{\prime\prime}_{1}T_{1}(x)=4c_{2}T_{0}(x)+24c_{3}T_{1}(x),$
(A.19)

where the truncation order of $u^{\prime\prime}(x)$ is $n-2$, and the
coefficients $c^{\prime\prime}_{0},\ldots,c^{\prime\prime}_{n-2}$ are
determined by $c_{2},\ldots,c_{n}$. Similarly for the case $d=2$, the
coefficients of $\Delta u(x)$ are determined by

$\displaystyle c^{\prime\prime}_{00}$ $\displaystyle=4c_{02}+4c_{20},$
$\displaystyle c^{\prime\prime}_{01}$ $\displaystyle=24c_{03}+4c_{21},$
$\displaystyle c^{\prime\prime}_{02}$ $\displaystyle=4c_{22},$ $\displaystyle
c^{\prime\prime}_{03}$ $\displaystyle=4c_{23},$ (A.20) $\displaystyle
c^{\prime\prime}_{10}$ $\displaystyle=4c_{12}+24c_{30},$ $\displaystyle
c^{\prime\prime}_{11}$ $\displaystyle=24c_{13}+24c_{31},$ $\displaystyle
c^{\prime\prime}_{12}$ $\displaystyle=24c_{32},$ $\displaystyle
c^{\prime\prime}_{13}$ $\displaystyle=24c_{33},$ $\displaystyle
c^{\prime\prime}_{20}$ $\displaystyle=4c_{22},$ $\displaystyle
c^{\prime\prime}_{21}$ $\displaystyle=24c_{23},$ $\displaystyle
c^{\prime\prime}_{30}$ $\displaystyle=4c_{32},$ $\displaystyle
c^{\prime\prime}_{33}$ $\displaystyle=24c_{33},$

so the matrix $D^{2}_{3}\otimes I+I\otimes D^{2}_{3}$ gives the linear system

$\setcounter{MaxMatrixCols}{16}\begin{pmatrix}&&4&&&&&&4&&&&&&&\\\
&&&24&&&&&&4&&&&&&\\\ &&&&&&&&&&4&&&&&\\\ &&&&&&&&&&&4&&&&\\\
&&&&&&4&&&&&&24&&&\\\ &&&&&&&24&&&&&&24&&\\\ &&&&&&&&&&&&&&24&\\\
&&&&&&&&&&&&&&&24\\\ &&&&&&&&&&4&&&&&\\\ &&&&&&&&&&&24&&&&\\\
&&&&&&&&&&&&&&4&\\\ &&&&&&&&&&&&&&&24\\\
\end{pmatrix}\begin{pmatrix}c_{0,0}\\\ c_{0,1}\\\ c_{0,2}\\\ c_{0,3}\\\
c_{1,0}\\\ c_{1,1}\\\ c_{1,2}\\\ c_{1,3}\\\ c_{2,0}\\\ c_{2,1}\\\ c_{2,2}\\\
c_{2,3}\\\ c_{3,0}\\\ c_{3,1}\\\ c_{3,2}\\\ c_{3,3}\\\
\end{pmatrix}=\begin{pmatrix}f_{0,0}\\\ f_{0,1}\\\ f_{0,2}\\\ f_{0,3}\\\
f_{1,0}\\\ f_{1,1}\\\ f_{1,2}\\\ f_{1,3}\\\ f_{2,0}\\\ f_{2,1}\\\ f_{3,0}\\\
f_{3,1}\\\ \end{pmatrix}.$ (A.21)

We now use the boundary conditions to complete the linear system. The
truncated Chebyshev series of the solution can be written as

$\displaystyle\gamma_{N}(x_{1})$
$\displaystyle=\sum_{k_{1}=0}^{2}{g_{N}}_{k_{1}}T_{k_{1}}(x_{1}),$ (A.22)
$\displaystyle\gamma_{S}(x_{1})$
$\displaystyle=\sum_{k_{1}=0}^{2}{g_{S}}_{k_{1}}T_{k_{1}}(x_{1}),$
$\displaystyle\gamma_{E}(x_{2})$
$\displaystyle=\sum_{k_{2}=0}^{2}{g_{E}}_{k_{1}}T_{k_{1}}(x_{2}),$
$\displaystyle\gamma_{W}(x_{2})$
$\displaystyle=\sum_{k_{2}=0}^{2}{g_{W}}_{k_{1}}T_{k_{1}}(x_{2}).$

We are given an oracle for preparing the state by interpolating the Chebyshev-
Gauss-Lobatto quadrature nodes specified in (4.4)

$\displaystyle\sum_{l_{1}=0}^{3}\gamma_{N}\mathopen{}\mathclose{{}\left(\cos\frac{\pi
l_{1}}{3}}\right)|l_{1}\rangle,$
$\displaystyle\sum_{l_{2}=0}^{3}\gamma_{E}\mathopen{}\mathclose{{}\left(\cos\frac{\pi
l_{2}}{3}}\right)|l_{2}\rangle,$ (A.23)
$\displaystyle\sum_{l_{1}=0}^{3}\gamma_{S}\mathopen{}\mathclose{{}\left(\cos\frac{\pi
l_{1}}{3}}\right)|l_{1}\rangle,$
$\displaystyle\sum_{l_{2}=0}^{3}\gamma_{W}\mathopen{}\mathclose{{}\left(\cos\frac{\pi
l_{2}}{3}}\right)|l_{2}\rangle.$

We perform the multi-dimensional inverse QCT on (A.23) to obtain

$\displaystyle\sum_{k_{1}=0}^{3}{g_{N}}_{k_{1}}|k_{1}\rangle,$
$\displaystyle\sum_{k_{2}=0}^{3}{g_{E}}_{k_{2}}|k_{2}\rangle,$ (A.24)
$\displaystyle\sum_{k_{1}=0}^{3}{g_{S}}_{k_{1}}|k_{1}\rangle,$
$\displaystyle\sum_{k_{2}=0}^{3}{g_{W}}_{k_{2}}|k_{2}\rangle,$

where $a_{k_{1},k_{2}}$ are defined by (4.33). The linear system from the
boundary conditions is

$\setcounter{MaxMatrixCols}{16}\begin{pmatrix}1&1&1&1&&&&&&&&&&&&\\\
&&&&1&1&1&1&&&&&&&&\\\ &&&&&&&&1&1&1&1&&&&\\\ &&&&&&&&&&&&1&1&1&1\\\
1&-1&1&-1&&&&&&&&&&&&\\\ &&&&1&-1&1&-1&&&&&&&&\\\ &&&&&&&&1&-1&1&-1&&&&\\\
&&&&&&&&&&&&1&-1&1&-1\\\ 1&&&&1&&&&1&&&&1&&&\\\ &1&&&&1&&&&1&&&&1&&\\\
&&1&&&&1&&&&1&&&&1&\\\ &&&1&&&&1&&&&1&&&&1\\\ 1&&&&-1&&&&1&&&&-1&&&\\\
&1&&&&-1&&&&1&&&&-1&&\\\ &&1&&&&-1&&&&1&&&&-1&\\\ &&&1&&&&-1&&&&1&&&&-1\\\
\end{pmatrix}\begin{pmatrix}c_{0,0}\\\ c_{0,1}\\\ c_{0,2}\\\ c_{0,3}\\\
c_{1,0}\\\ c_{1,1}\\\ c_{1,2}\\\ c_{1,3}\\\ c_{2,0}\\\ c_{2,1}\\\ c_{2,2}\\\
c_{2,3}\\\ c_{3,0}\\\ c_{3,1}\\\ c_{3,2}\\\ c_{3,3}\\\
\end{pmatrix}=\begin{pmatrix}{g_{N}}_{0}\\\ {g_{N}}_{1}\\\ {g_{N}}_{2}\\\
{g_{N}}_{3}\\\ {g_{S}}_{0}\\\ {g_{S}}_{1}\\\ {g_{S}}_{2}\\\ {g_{S}}_{3}\\\
{g_{E}}_{0}\\\ {g_{E}}_{1}\\\ {g_{E}}_{2}\\\ {g_{E}}_{3}\\\ {g_{W}}_{0}\\\
{g_{W}}_{1}\\\ {g_{W}}_{2}\\\ {g_{W}}_{3}\\\ \end{pmatrix}.$ (A.25)

Adding the two linear systems (A.21) and (A.25) together, we obtain a full-
rank linear system

$\overline{D}^{(2)}_{3}\otimes I+I\otimes\overline{D}^{(2)}_{3},$ (A.26)

where

$\overline{D}_{3}^{(2)}=\begin{pmatrix}0&0&4&0\\\ 0&0&0&24\\\ 1&-1&1&-1\\\
1&1&1&1\\\ \end{pmatrix}.$ (A.27)

In summary, the linear system including the differential equations and the
boundary conditions is

$\small\setcounter{MaxMatrixCols}{16}\begin{pmatrix}&&4&&&&&&4&&&&&&&\\\
&&&24&&&&&&4&&&&&&\\\ 1&-1&1&-1&&&&&&&4&&&&&\\\ 1&1&1&1&&&&&&&&4&&&&\\\
&&&&&&4&&&&&&24&&&\\\ &&&&&&&24&&&&&&24&&\\\ &&&&1&-1&1&-1&&&&&&&24&\\\
&&&&1&1&1&1&&&&&&&&24\\\ 1&&&&-1&&&&1&&4&&-1&&&\\\ &1&&&&-1&&&&1&&24&&-1&&\\\
&&1&&&&-1&&1&-1&2&-1&&&-1&\\\ &&&1&&&&-1&1&1&1&2&&&&-1\\\
1&&&&1&&&&1&&&&1&&4&\\\ &1&&&&1&&&&1&&&&1&&24\\\ &&1&&&&1&&&&1&&1&-1&2&-1\\\
&&&1&&&&1&&&&1&1&1&1&2\\\ \end{pmatrix}\begin{pmatrix}c_{0,0}\\\ c_{0,1}\\\
c_{0,2}\\\ c_{0,3}\\\ c_{1,0}\\\ c_{1,1}\\\ c_{1,2}\\\ c_{1,3}\\\ c_{2,0}\\\
c_{2,1}\\\ c_{2,2}\\\ c_{2,3}\\\ c_{3,0}\\\ c_{3,1}\\\ c_{3,2}\\\ c_{3,3}\\\
\end{pmatrix}=\begin{pmatrix}f_{0,0}\\\ f_{0,1}\\\ f_{0,2}+{g_{S}}_{0}\\\
f_{0,3}+{g_{N}}_{0}\\\ f_{1,0}\\\ f_{1,1}\\\ f_{1,2}+{g_{S}}_{1}\\\
f_{1,3}+{g_{N}}_{1}\\\ f_{2,0}+{g_{W}}_{0}\\\ f_{2,1}+{g_{W}}_{1}\\\
{g_{W}}_{2}+{g_{S}}_{2}\\\ {g_{W}}_{3}+{g_{N}}_{2}\\\ f_{3,0}+{g_{E}}_{0}\\\
f_{3,1}+{g_{E}}_{1}\\\ {g_{E}}_{2}+{g_{S}}_{3}\\\ {g_{E}}_{3}+{g_{N}}_{3}\\\
\end{pmatrix}.$ (A.28)

##  Appendix 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.

See 7

###### Proof.

By direct calculation of the $l_{\infty}$ norm (i.e., the maximum absolute
column sum) of (4.56), we have

$\|\overline{D}^{(2)}_{n}\|_{\infty}\leq\mathopen{}\mathclose{{}\left(\frac{(n+1)\pi}{2}}\right)^{2}\leq(2n)^{2}.$
(B.1)

Then the inverse of the matrix (4.56) is

$\displaystyle[(\overline{D}^{(2)}_{n})^{-1}]_{k,k}=-\frac{1}{((k-\mathopen{}\mathclose{{}\left\lfloor
n/2}\right\rfloor)\pi)^{2}},\quad$ $\displaystyle
k\in[{n+1}]_{0}\backslash\\{\mathopen{}\mathclose{{}\left\lfloor
n/2}\right\rfloor\\},$ (B.2)
$\displaystyle[(\overline{D}^{(2)}_{n})^{-1}]_{\mathopen{}\mathclose{{}\left\lfloor
n/2}\right\rfloor,k}=\frac{1}{((k-\mathopen{}\mathclose{{}\left\lfloor
n/2}\right\rfloor)\pi)^{2}},\quad$ $\displaystyle k\in[{n+1}]_{0}$

as can easily be verified by a direct calculation.

By direct calculation of the Frobenius norm of (4.56), we have

$\|(\overline{D}^{2}_{n})^{-1}\|_{F}^{2}\leq
1+2\sum_{k=1}^{\infty}\frac{1}{k^{4}\pi^{4}}=1+\frac{2}{\pi^{4}}\frac{\pi^{4}}{90}\leq
2.$ (B.3)

Thus we have the result in (4.58):

$\displaystyle\sigma_{\max}(\overline{D}^{(2)}_{n})\leq\sqrt{n+1}\|\overline{D}^{2}_{n}\|_{\infty}\leq(2n)^{2.5},$
(B.4)
$\displaystyle\sigma_{\min}(\overline{D}^{(2)}_{n})\geq\frac{1}{\|(\overline{D}^{2}_{n})^{-1}\|_{F}}\geq\frac{1}{\sqrt{2}}$

as claimed. ∎

See 8

###### Proof.

By direct calculation of the Frobenius norm of (4.57), we have

$\|\overline{D}^{(2)}_{n}\|_{F}^{2}\leq
n^{2}\max_{k,r}\mathopen{}\mathclose{{}\left(\frac{r(r^{2}-k^{2})}{\sigma_{k}}}\right)^{2}\leq
n^{2}\cdot n^{6}=n^{8}.$ (B.5)

Next we upper bound $\|(\overline{D}^{(2)}_{n})^{-1}\|$. By definition,

$\|(\overline{D}^{(2)}_{n})^{-1}\|=\sup_{\|b\|\leq
1}\|(\overline{D}^{2}_{n})^{-1}b\|.$ (B.6)

Given any vector $b$ satisfying $\|b\|\leq 1$, we estimate $\|x\|$ defined by
the full-rank linear system

$\overline{D}^{(2)}_{n}x=b.$ (B.7)

Notice that $\overline{D}^{(2)}_{n}$ is the sum of the upper triangular matrix
${D}^{2}_{n}$ and (4.53), the coordinates $x_{2},\ldots,x_{n}$ are only
defined by coordinates $b_{0},\ldots,b_{n-2}$. So we only focus on the partial
system

${D}^{(2)}_{n}[0,0,x_{2},\ldots,x_{n}]^{T}=[b_{0},\ldots,b_{n-2},0,0]^{T}.$
(B.8)

Given the same $b$, we also define the vector $y$ by

${D}_{n}[0,y_{1},\ldots,y_{n-1},0]^{T}=[b_{0},\ldots,b_{n-2},0,0]^{T},$ (B.9)

where each coordinate of $y$ can be expressed by

$b_{k}=\sum_{l=1}^{n-1}[{D}_{n}]_{kl}y_{l}=\sum_{l=1}^{n-1}\frac{2l}{\sigma_{k}}y_{l},\quad\text{$k+l$
odd},~{}l>k,~{}k\in[{n-1}]_{0}.$ (B.10)

Using this equation with $k=l-1$ and $k=l+1$, we can express $y_{l}$ in terms
of $b_{l-1}$ and $b_{l+1}$:

$\frac{2l}{\sigma_{l-1}}y_{l}=b_{l-1}-\frac{1}{\sigma_{l-1}}b_{l+1},\quad
l\in[{n-1}],$ (B.11)

where we let $b_{n-1}=b_{n}=0$. Thus we have

$\displaystyle\sum_{l=1}^{n-1}y_{l}^{2}$
$\displaystyle=\sum_{l=1}^{n-1}\mathopen{}\mathclose{{}\left(\frac{\sigma_{l-1}}{2l}\mathopen{}\mathclose{{}\left(b_{l-1}-\frac{1}{\sigma_{l-1}}b_{l+1}}\right)}\right)^{2}$
(B.12)
$\displaystyle\leq\sum_{l=1}^{n-1}\frac{\sigma_{l-1}^{2}}{4l^{2}}\mathopen{}\mathclose{{}\left(1+\frac{1}{\sigma_{l-1}^{2}}}\right)(b_{l-1}^{2}+b_{l+1}^{2})$
$\displaystyle\leq\frac{5}{4}(b_{0}^{2}+b_{2}^{2})+\frac{1}{16}\sum_{l=2}^{n-2}(b_{l-1}^{2}+b_{l+1}^{2})$
$\displaystyle\leq 2\sum_{l=0}^{n-2}b_{l}^{2}.$

We notice that $y$ also satisfies

$[0,y_{1},\ldots,y_{n-1},0]^{T}={D}_{n}[0,0,x_{2},\ldots,x_{n}]^{T},$ (B.13)

where each coordinate of $y$ can be expressed by

$y_{l}=\sum_{r=1}^{n}[{D}_{n}]_{lr}x_{r}=\sum_{r=1}^{n}\frac{2r}{\sigma_{l}}x_{r},\quad\text{$l+r$
odd},~{}r>l,~{}l\in[{n-1}].$ (B.14)

Substituting the $(r-1)$st and the $(r+1)$st equations of (B.14), we can
express $x$ in terms of $y$:

$\frac{2r}{\sigma_{r-1}}x_{r}=y_{r-1}-\frac{1}{\sigma_{r-1}}y_{r+1},\quad
r\in[{n}]\backslash\\{1\\},$ (B.15)

where we let $y_{n}=y_{n+1}=0$. Similarly, according to (B.15), we also have

$\sum_{l=2}^{n}x_{l}^{2}\leq 2\sum_{l=1}^{n-1}y_{l}^{2}.$ (B.16)

Then we calculate $x_{0}^{2}+x_{1}^{2}$ based on the last two equations of
(B.7), (B.12), and (B.15), giving

$\displaystyle x_{0}^{2}+x_{1}^{2}$
$\displaystyle=\frac{1}{2}[(x_{0}+x_{1})^{2}+(x_{0}-x_{1})^{2}]$ (B.17)
$\displaystyle=\frac{1}{2}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(b_{n}-\sum_{l=2}^{n}x_{l}}\right)^{2}+\mathopen{}\mathclose{{}\left(b_{n-1}-\sum_{l=2}^{n}(-1)^{l}x_{l}}\right)^{2}}\right]$
$\displaystyle=\frac{1}{2}\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left(b_{n}-\sum_{l=2}^{n}\frac{\sigma_{l-1}}{2l}\mathopen{}\mathclose{{}\left(y_{l-1}-\frac{1}{\sigma_{l-1}}y_{l+1}}\right)}\right)^{2}+\mathopen{}\mathclose{{}\left(b_{n-1}-\sum_{l=2}^{n}(-1)^{l}\frac{\sigma_{l-1}}{2l}\mathopen{}\mathclose{{}\left(y_{l-1}-\frac{1}{\sigma_{l-1}}y_{l+1}}\right)}\right)^{2}}\right]$
$\displaystyle\leq\frac{1}{2}\mathopen{}\mathclose{{}\left(1+\sum_{l=2}^{n}\frac{\sigma_{l-1}^{2}}{4l^{2}}}\right)\mathopen{}\mathclose{{}\left[b_{n}^{2}+\sum_{l=2}^{n}\mathopen{}\mathclose{{}\left(y_{l-1}-\frac{1}{\sigma_{l-1}}y_{l+1}}\right)^{2}+b_{n-1}^{2}+\sum_{l=2}^{n}\mathopen{}\mathclose{{}\left(y_{l-1}-\frac{1}{\sigma_{l-1}}y_{l+1}}\right)^{2}}\right]$
$\displaystyle\leq\frac{1}{2}\mathopen{}\mathclose{{}\left(1+\frac{1}{4}\sum_{l=2}^{\infty}\frac{1}{l^{2}}}\right)\mathopen{}\mathclose{{}\left[b_{n}^{2}+b_{n-1}^{2}+\sum_{l=2}^{n}\mathopen{}\mathclose{{}\left(1+\frac{1}{\sigma_{l-1}^{2}}}\right)(y_{l-1}^{2}+y_{l+1}^{2})}\right]$
$\displaystyle\leq\frac{1}{2}\mathopen{}\mathclose{{}\left(1+\frac{\pi^{2}}{24}}\right)\mathopen{}\mathclose{{}\left[b_{n}^{2}+b_{n-1}^{2}+4\sum_{l=1}^{n-1}y_{l}^{2}}\right]$
$\displaystyle\leq b_{n}^{2}+b_{n-1}^{2}+8\sum_{l=0}^{n-2}b_{l}^{2}.$

Thus, based on (B.12), (B.16), and (B.17), the inequality

$\displaystyle\sum_{l=0}^{n}x_{l}^{2}$
$\displaystyle=x_{0}^{2}+x_{1}^{2}+\sum_{l=2}^{n}x_{l}^{2}$ (B.18)
$\displaystyle\leq
b_{n}^{2}+b_{n-1}^{2}+8\sum_{l=0}^{n-2}b_{l}^{2}+4\sum_{l=0}^{n-2}b_{l}^{2}$
$\displaystyle\leq b_{n}^{2}+b_{n-1}^{2}+12\sum_{l=0}^{n-2}b_{l}^{2}\leq 12$

holds for any vectors $b$ satisfying $\|b\|\leq 1$. Thus

$\|(\overline{D}^{(2)}_{n})^{-1}\|=\sup_{\|b\|\leq 1}\|x\|\leq 12<16.$ (B.19)

Altogether, we have

$\displaystyle\sigma_{\max}(\overline{D}^{(2)}_{n})\leq\|\overline{D}^{2}_{n}\|_{F}\leq
n^{4},$ (B.20)
$\displaystyle\sigma_{\min}(\overline{D}^{(2)}_{n})\geq\frac{1}{\|(\overline{D}^{2}_{n})^{-1}\|}\geq\frac{1}{16}$

as claimed in (4.59). ∎