widerequation

$\BODY$ (1)

# Quantum algorithms for computing observables of nonlinear partial
differential equations

Shi Jin  Institute of Natural Sciences, Shanghai Jiao Tong University,
Shanghai 200240, China  School of Mathematical Sciences, Shanghai Jiao Tong
University, Shanghai 200240, China  Ministry of Education, Key Laboratory in
Scientific and Engineering Computing, Shanghai Jiao Tong University, Shanghai
200240, China  Nana Liu  Corresponding author<EMAIL_ADDRESS>Institute of Natural Sciences, Shanghai Jiao Tong University, Shanghai 200240,
China  Ministry of Education, Key Laboratory in Scientific and Engineering
Computing, Shanghai Jiao Tong University, Shanghai 200240, China  University
of Michigan-Shanghai Jiao Tong University Joint Institute, Shanghai 200240,
China

###### Abstract

We construct quantum algorithms to compute physical observables of nonlinear
PDEs with $M$ initial data. Based on an exact mapping between nonlinear and
linear PDEs using the level set method, these new quantum algorithms for
nonlinear Hamilton-Jacobi and scalar hyperbolic PDEs can be performed with a
computational cost that is independent of $M$, for arbitrary nonlinearity.
Depending on the details of the initial data, it can also display up to
exponential advantage in both the dimension of the PDE and the error in
computing its observables. For general nonlinear PDEs, quantum advantage with
respect to $M$ is possible in the large $M$ limit.

##  I Introduction

Nonlinear ordinary and partial differential equations (ODEs and PDEs) have
been central to modelling of some of the most significant problems in physics,
chemistry, engineering, biology and finance, including climate modelling,
aircraft design, molecular dynamics and drug design, deep learning neural
networks and financial markets. In physics, the most important mathematical
equations–from quantum mechanics, classical mechanics to kinetic theory and
hydrodynamics–are all modelled by linear or nonlinear (integro)-differential
equations. Although quantum algorithms can be potentially advantageous for
certain linear problems like linear PDEs (e.g. [1, 2, 3, 4, 5, 6]) it is still
unclear to what extent quantum algorithms can be leveraged for nonlinear
problems. Although quantum mechanics itself is fundamentally linear (as far as
we know), most natural phenomena–and their associated mathematical equations
or models–are nonlinear, hence the ability to simulate nonlinear
problems–including the nonlinear PDEs–will significantly extend the horizon of
quantum computing.

The most natural way to approach a nonlinear problem using quantum algorithms
is to find a way to represent the nonlinear problem in a linear way, where
quantum computational advantage in the former problem can still be maintained.
Here we distinguish between two types of approaches that converts a nonlinear
PDE into a linear PDE. One approach involves approximations (e.g. either
through linearisation of the nonlinearity or through discretisation). They
include the linear approximation and the linear representation of nonlinear
ODEs methods. The second approach is the linear representation for nonlinear
PDEs, where no approximations are required to map between the nonlinear and
linear PDEs.

In the linear approximation approach, in which the nonlinear term is
linearized, errors are introduced, so the approach may only be valid for a
short time, for weak nonlinearities and consequently may lose significant
nonlinear features of the problem after a long time. On the other hand,
approaches like Carlemann linearisation [7] or in [8, 9], require that the
size of the corresponding linear problem to increase in a way that is
dependent on the degree of nonlinearity and is also restricted to polynomial
nonlinearities. This requirement of expanding into ever higher dimensions to
deal with strong nonlinearities comes at a sacrifice of the (sometimes
significant) resource cost. For general nonlinear functions that do not have
low-order polynomial expansions, this can be infeasible and important
nonlinear features can be lost in the long time limit due to the truncation in
the Carlemann linearization. This is similarly true for methods in [8, 9].

The Koopman-von Neumann approach [10] (a similar approach taken in [11]), is a
type of linear representation method for nonlinear ODEs, though not for
nonlinear PDEs. This means it allows one to convert nonlinear ODEs into a
linear transport equation, without any loss of information, thus is valid
globally in time, for any nonlinearity. However, to apply this method to solve
nonlinear PDEs, one needs to first discretise its spatial derivatives, giving
rise to a system of nonlinear ODEs which can then be implemented via the
Koopman-von Neumann approach. This, unfortunately, would give rise to ODE
systems with exceedingly large dimensions, where the dimension depends on the
discretisation error, hence the quantum advantages can be lost. Since one
needs to discretize the nonlinear PDE first, in this regard this method is not
a linear representation method for PDEs, as we defined above.

In this paper, for several important classes of nonlinear PDEs, including the
Hamilton-Jacobi equations and scalar nonlinear hyperbolic equations, we map
the underlying nonlinear PDEs of $(d+1)$-dimension to a (not more than)
$(2d+1)$-dimensional linear PDEs, by utilizing the level set formalism [12].
One can then solve the linear PDEs – whose dimension is at most twice that of
the original nonlinear PDE – on a quantum computer, without losing any
physical information. This procedure is exact, so no approximations are made.
This is an example of linear representation method for nonlinear PDEs.

We also find another pleasing bonus of this approach, thanks to the linearity
being valid globally in time, in that one can run problems with many ($M>>1$)
different initial data, and obtain ensemble averages of the physical
observables at a later time, by solving the PDE just once!

Figure 1: Computing observables from $(d+1)$-dimensional nonlinear PDE with
$M$ initial data (a) Linear representation for nonlinear PDE approach (Green):
Using the level set formalism, $(d+1)$-dimensional Hamilton-Jacobi PDEs
(Section III) with $M$ initial data are mapped onto a $(2d+1)$-dimensional
linear transport PDE with a single initial datum. Similarly,
$(d+1)$-dimensional scalar hyperbolic PDEs (Section IV) with $M$ initial data
are mapped onto a $(d+2)$-dimensional linear transport PDE with a single
initial datum; (b) Linear representation for nonlinear ODE approach (Blue):
$(d+1)$-dimensional general nonlinear PDEs can be mapped onto $D$ nonlinear
ODEs with $M$ initial data using either the Lagrangian or Eulerian
discretisation methods (Section VI). The system of $D$ nonlinear ODEs with $M$
initial data can be mapped onto a $(D+1)$-dimensional linear transport PDE
with a single initial datum and ensemble averages can be computed with a
quantum algorithm (Section V).

Using this linear representation for nonlinear PDEs, we develop quantum
algorithms to compute physical observables of general $(d+1)$-dimensional
nonlinear Hamilton-Jacobi and scalar hyperbolic PDEs with multiple ($M$)
initial data. In addition to being of practical interest, computing
observables also allows one to compare quantum computational costs with
classical computational costs on an equal footing, since they are solving the
same problem. This is in contrast to many previous quantum algorithms for ODEs
and PDEs (e.g. [13, 14, 2]) that focus only on the cost of quantum state
preparation of the PDE solutions, without considering (potentially large)
costs in extracting those solutions and interpreting them. Furthermore, this
level set linear representation allows us to introduce a new embedding of
classical data into a quantum state which we call the level set embedding.
This embedding is shown to be advantageous in computing physical observables
compared to the amplitude embedding commonly used.

When $M=1$, we show that, for certain classes of initial data, there can be up
to exponential quantum advantage in computing observables to precision
$\epsilon$, with respect to both $d$ and $\epsilon$, with no large overheads
in any other parameter, including time $T$. The quantum cost is also
independent of the nature of the nonlinearity. For $M>1$, our algorithm
computes ensemble averages over $M$ different initial data in one computation,
instead of computing the problem $M$-times. Hence the quantum resource cost is
independent of $M$, while classical costs are linear in $M$. The same
advantages in $d$ and $\epsilon$ still apply.

There are important applications where computing a PDE or ODE with multiple
initial data ($M>1$) is of interest. For example, one needs to run numerical
simulations with many different initial data to obtain ensemble averaged
solutions in uncertainty quantification with random initial data, using Monte-
Carlo sampling techniques or stochastic collocation [15]. In Rayleigh-Taylor
instability one needs to run many experiments with different initial data to
obtain an ensemble averaged numerical solution that can converge (to the
Young-measure solution [16]). In geometric optics one needs to solve multiple
rays that follows a Hamiltonian system (eikonal equation) with different
initial data [17]. In quantum wave packet methods for quantum dynamics
simulations, one needs to solve for multiple Gaussian wave packets, each
evolving by a system of nonlinear ODEs [18], and later seek their linear
superposition. We show that our quantum algorithms always have quantum
advantage in $M$, under various conditions.

For more general nonlinear PDEs, the idea of first discretising the spatial
derivatives to convert them into system of nonlinear ODEs is also discussed.
Two approximation techniques, the Lagrangian methods, such as particle, vortex
or mesh-free methods that do not use grids, and Eulerian methods which use
grid-based discretisations of spatial derivatives, are studied. We show that
using this approach, there can only be quantum advantage in the large $M$
limit, while there are no quantum advantages in $d$ and $\epsilon$. Thus
computing ensemble averages for general nonlinear PDEs using quantum
algorithms still requires further work to demonstrate useful quantum
advantages.

The outline of the paper is as follows. We begin in section II with the
background to nonlinear Hamilton-Jacobi and scalar hyperbolic PDEs and a
system of ODEs with multiple initial conditions. We finish with a basic
summary of quantum protocols for solving linear systems of equations. Readers
familiar with these areas can skip this section. In Section III we present the
linear representation method using level sets and the quantum algorithm used
to compute the physical observables. We compare the classical and quantum
computational costs and derive the conditions under which there are quantum
advantages. In Section IV we repeat the analysis for scalar nonlinear
hyperbolic equations. The quantum query and gate complexities for computing
ensemble averages of a system of nonlinear ODEs is presented in Section V, and
the application of this method to more general nonlinear PDEs is in Section
VI.

In Figure 1 we give a basic outline of the methods considered in this paper.
Table 1 gives a summary of the comparison between classical and quantum
resource costs in computing ensemble averages for nonlinear PDEs and ODEs
studied in this paper. Throughout the paper, we always use first order
approximations to space and time derivatives for all systems under
consideration. Extensions to higher order approximations can be done in a
straightforward way and won’t be pursued here.

Table 1: Quantum ($\mathcal{Q}$) and classical ($\mathcal{C}$) cost comparison
in computing observables at time $T$ to precision $\epsilon$

$\displaystyle\text{PDE:
}\mathcal{O}\left(\frac{\mathcal{C}}{\mathcal{Q}}\right)=\tilde{O}\left(\frac{M}{T^{2}}d^{r_{1}}\left(\frac{1}{\epsilon}\right)^{r_{2}}\right)$
$\displaystyle\text{ODE:
}\mathcal{O}\left(\frac{\mathcal{C}}{\mathcal{Q}}\right)=\tilde{O}\left(\frac{M}{T^{2}}D^{r_{1}}\left(\frac{1}{\epsilon}\right)^{r_{2}}\right)$

Nonlinear equations | $r_{1}$ |  | $r_{2}$ |  $b$ range |  | Quantum
---|---|---|---|---|---|---
($M$ initial data) |  |  |  | (initial data- |  | advantage
|  |  |  | dependent) |  | (possible)
$(d+1)$-dimensional | $d-4-b$ |  | $d-9-3b$ | $b\in[0,\frac{d}{3}-3)$ |  |  $M$, $d$, $\epsilon$
Hamilton-Jacobi PDE |  |  |  |  |  |
$(d+1)$-dimensional | $d-5-b$ |  | $d-9-3b$ | $b\in[0,\frac{d}{3}-3)$ |  |  $M$, $d$, $\epsilon$
hyperbolic PDE |  |  |  |  |  |
System of $D$ ODEs | $-5$ |  | $-9$ | $b=0$ |  | $M$
$(d+1)$-dimensional general PDE | $-7$ |  | $-13$ | $b=0$ |  | $M$
(Lagrangian discretisation) |  |  |  |  |  |
$(d+1)$-dimensional general PDE | $-4d$ |  | $-9-4d$ | $b=0$ |  | Large $M$
(Eulerian discretisation) |  |  |  |  |  |

##  II Background

We introduce nonlinear PDEs, in particular Hamilton-Jacobi and hyperbolic
PDEs, and a system of nonlinear ODEs. We also provide the background of the
relevant quantum algorithms.

###  II.1 Nonlinear PDEs

Nonlinear partial difference equations can be written in the following general
form

$\frac{\partial u}{\partial t}+F(u,\nabla u,\nabla^{2}u,\cdots)=0,\quad
t\in\mathbb{R}^{+},x\in\mathbb{R}^{d},u\in\mathbb{R}^{d}$ (2)

Here $t\geq 0$ is time, $x$ is the spatial variable, while $F$ is a nonlinear
function or functional.

For most differential equations, linear or nonlinear, numerical computations
have been the most important tools to solve them, since analytical solutions
seldom exist. In fact, classical algorithms–referring to numerical algorithms
using classical computers–for differential equations have been among the most
important achievements in scientific computing in more than half a decade.

Different partial differential equations have drastically different behavior
in their solutions, from different regularities (smoothness) to different
physical behaviors (e.g. conservation, invariances, entropy conditions), thus
they call for drastically different numerical strategies. In this paper, we
will mainly focus on first order quasi-linear PDEs (such as nonlinear
Hamilton-Jacobi and hyperbolic equations).

For more general PDEs, one can write it as a system of ODEs by first
discretising in space. There are two classes of discretisation methods. One is
the so-called Eulerian framework, in which one discretises the spatial
variables by a finite difference, finite element, finite volume, or spectral
method, on a fixed mesh. The advantage of Eulerian methods is their high order
accuracy but they suffer from the curse-of-dimensionality in high spatial
dimensions. The other is the Lagrangian method, including the so-called
particle or mesh-free method. This is popular for high dimensional problems,
for example particle or Monte-Carlo methods for kinetic equations (the
Boltzmann equation and Vlasov type equations), and vortex methods for
incompressible Euler and Navier-Stokes equations in fluid dynamics. The
Lagrangian methods do not suffer from the curse-of-dimensionality but they are
of lower order (typically only half to first order) methods. We briefly
discuss the two methods in Section VI.

In our notation throughout the paper, we use $O(\cdot)$ to denote the case
where all constant factors are suppressed unless otherwise stated.
$\mathcal{O}(\cdot)$ is the standard big-O notation denoting a tight
asymptotic upper bound and ${\scriptstyle\mathcal{O}}(\cdot)$ is the standard
little-O notation, which is a looser asymptotic upper bound.
$\tilde{O}(\cdot)$ indicates that in addition, all logarithmic factors are
suppressed. The notation $``(\cdot)\sim(\cdot)"$ denotes equivalence with
constant factors ignored.

####  II.1.1 Hamilton-Jacobi PDEs

A distinct feature of first-order quasi-linear PDEs, including Hamilton-Jacobi
or nonlinear hyperbolic equations, is that solutions can become singular even
if the initial data is smooth. For nonlinear hyperbolic equations, shocks
(solution becomes discontinuous) may develop, while for Hamilton-Jacobi
equations the solution may form cusps (at which points the derivatives of the
solution becomes discontinuous) corresponding to caustics in geometric optics.

When solutions become singular, one needs to make sense of the equations since
the derivative terms are no longer well-defined. In this regard one either
uses the notion of viscosity solutions, which are physically relevant in
applications such in gas dynamics for compressible Euler equations [19], or
optimal control using the Hamilton-Jacobi-Bellman equation [20, 21, 22].
Another notion is the multi-valued solutions, which are relevant to
applications such as semi-classical quantum dynamics, and geometric optics, in
which the dynamics is time-irreversible and the solution satisfies the linear
superposition principle [23, 24, 25]. When the solutions are smooth, both
notions define the same solution, but not while singularities emerge. It is
the latter case that will be considered in this paper.

Hamilton-Jacobi equations arise for instance in geometric optics, the
semiclassical limit of the Schrödinger equation, the level set formulation of
front propagation, optimal control, mean-field games, sticky particles or
pressureless gases and KPZ equations. It has the following general form

$\displaystyle\partial_{t}S^{[k]}+H(\nabla S^{[k]},x)=0,\quad
t\in\mathbb{R}^{+},x\in\mathbb{R}^{d},S^{[k]}(t,x)\in\mathbb{R},$ (3)
$\displaystyle S^{[k]}(0,x)=S_{0}^{[k]}(x),\quad k=1,...,M$ (4)

subject to $M$ different initial data, where the same Hamilton-Jacobi PDE is
satisfied for each $k$.

A direct numerical approximation to Eq. (3) usually gives rise to the so-
called viscosity solution [20, 21]. This notion of the solution is not valid
in geometric optics (multiple arrivals in seismic waves [26], for example), in
the semiclassical limit of quantum dynamics, or in the high frequency limit of
linear wave equations (elastic waves, electromagnetic waves, etc.) [27, 24,
28], since it violates the linear superposition principle. In these
applications, one is interested in computing the multi-valued solution, [24,
12, 25, 29]. There are several classes of algorithms that were developed to
capture such solutions:

  * •

Ray tracing [30, 31]. This is based on solving the characteristics of the
system, which is a Hamiltonian system

$\partial_{t}x(t)=\nabla_{p}H(x,p)\,,\quad\partial_{t}p(t)=-\nabla_{x}H(x,p)\,.$
(5)

The advantage of this method is its simplicity since it just solves a system
of ODEs. The disadvantage of such Lagrangian type methods is that particles
are not uniformly distributed, hence at later time there may be regions where
there are not enough particles to guarantee numerical accuracy and one needs
to add more particles and then use interpolations to define these particle,
which are quite delicate.

  * •

Moment methods [32, 12]. For multivalued solutions one can use moment systems,
which are superposition of Hamilton-Jacobi equations. The advantage of this
method is that one stays in the physical space, however, the moment systems
are difficult to derive in higher-dimensions and, in particular, one needs to
know, a priori, the number of branches in order to have the right number of
moments. This is pretty much an impossible task.

  * •

Level set formulation. Here one builds the gradient of $S$ of Eq. (3) into the
zero level sets (defined later) of functions which solves a system of
Liouville equations [25, 33, 29]. This method is globally valid and one solves
a linear system of PDEs. Its disadvantage, for classical computers, is the
curse-of-dimensionality since the equations are defined in the phase space,
hence the dimension is doubled. Since the curse-of-dimensionality can in cases
be resolved with quantum computers, and thanks to the linearity of the system,
this is the approach we advocate in this paper.

Define $u^{[k]}=\nabla S^{[k]}\in\mathbb{R}^{d}$. Then $u^{[k]}$ solves a
hyperbolic system of conservation laws in gradient form:

$\displaystyle\partial_{t}u^{[k]}+\nabla H(u^{[k]},x)=0,$ (6) $\displaystyle
u^{[k]}(0,x)=\nabla S^{[k]}_{0}(x).$ (7)

Two classical examples of Hamiltonians are

$H(x,p)=\frac{p^{2}}{2}+V(x),$ (8)

corresponding to classical Newtonian particles, and

$H(x,p)=c(x)|p|$ (9)

that arises in geometric optics and the level set formulation of front
propagation [27, 21], in which $c(x)$ is the reciprocal of the index of
reflection or the speed of a propagating front in the normal direction.

When we encounter multi-valued solutions of Eq. (6), we define an ensemble
average in Section III, which can be identified as physical observables of the
system. Then we show how it can be computed with a quantum algorithm. The
classical cost for computing these observables is in the following lemma.

######  Lemma 1.

If a finite difference or finite volume method is used based on a regular mesh
using $N_{\text{HJ}}$ spatial points in each dimension, with mesh size
$h_{\text{HJ}}=1/N_{\text{HJ}}$, and time step $\Delta
t_{\text{HJ}}=O(h_{\text{HJ}}/d)$ (due to the CFL stability condition), then
the computational cost of solving the Hamilton-Jacobi equation Eq. (3) is
$O(dN_{t,\text{HJ}}N_{\text{HJ}}^{d})=O(d^{2}TN_{\text{HJ}}^{d+1})$. To reach
an error tolerance of $\epsilon_{\text{HJ}}$ one needs
$N_{\text{HJ}}=O(d/\epsilon_{\text{HJ}})$, hence the total cost is
$O(Td^{d+3}(1/\epsilon_{\text{HJ}})^{d+1})$. If one is interested in computing
the ensemble average of $S^{[k]},(k=1,\cdots,M)$ in the presence of $M$
different initial data, then the cost will be
$O(MTd^{d+3}(1/\epsilon_{\text{HJ}})^{d+1})$. The cost of computing the
ensemble average of $u^{[k]}$ is $O(MTd^{d+4}(1/\epsilon_{\text{HJ}})^{d+1})$.

###### Proof.

The cost of $O(dN_{t,\text{HJ}}N_{\text{HJ}}^{d})$ is obtained by directly
counting. If one computes to time $T=N_{t,{\text{HJ}}}\Delta t_{\text{HJ}}$,
since $\Delta t=O(h_{\text{HJ}}/d)$ (due to the CFL condition) and
$h_{\text{HJ}}=1/N_{\text{HJ}}$ one also obtains
$O(d^{2}TN_{\text{HJ}}^{d+1})$. Since the truncation error of a first order
method to approximate the derivative in each dimension is of
$O(1/N_{\text{HJ}})$, and there are $d$ differential operators to be
discretised, the total truncation error is of $O(d/N_{\text{HJ}})$. Hence to
reach the error $\epsilon_{\text{HJ}}$ one needs
$N_{\text{HJ}}=O(d/\epsilon_{\text{HJ}})$. This means the cost is
$O(d^{2}TN_{\text{HJ}}^{d+1})=O(Td^{d+3}(1/\epsilon_{\text{HJ}})^{d+1})$. If
one is interested in computing the ensemble average of
$S^{[k]}.(k=1,\cdots,M)$ in the presence of $M$ different initial data, the
total cost will then be $O(MTd^{d+3}(1/\epsilon_{\text{HJ}})^{d+1})$. The cost
of computing the ansemble average of $u^{[k]}$ is then
$O(MTd^{d+4}(1/\epsilon_{\text{HJ}})^{d+1})$ since $u^{[k]}$ is
$d$-dimensional. ∎

Remark: Spectral methods are usually not used for Hamilton-Jacobi solutions
since the solutions develop singularities (caustics) and spectral methods
introduce numerical oscillations.

####  II.1.2 Nonlinear scalar hyperbolic PDEs

Nonlinear hyperbolic PDEs arise for instance in gas dynamics, combustion,
magnetohydrodynamics, shallow water and traffic flows. Here we focus on the
scalar equation where $u^{[k]}(t,x)\in\mathbb{R}$ is a scalar solving an
initial value problem of an $(d+1)$-dimensional first-order hyperbolic PDE
with a non-zero source term

$\displaystyle\partial_{t}u^{[k]}+F(u^{[k]})\cdot\nabla_{x}u^{[k]}+Q(x,u^{[k]})=0,\quad
t\in\mathbb{R}^{+},\quad x\in\mathbb{R}^{d},$ (10) $\displaystyle
u^{[k]}(0,x)=u^{[k]}_{0}(x),\quad k=1,...,M.$ (11)

Here $F(u^{[k]}):\mathbb{R}\to\mathbb{R}^{d}$ is a vector and
$Q:\mathbb{R}^{d+1}\to\mathbb{R}$ is the source term. This equation includes
any such hyperbolic PDE in conservative or non-conservative form. The cost in
solving this equation is essentially the same as those in Lemma 1 for
$S^{[k]}$.

######  Lemma 2.

If a finite difference or finite volume method is used with a regular mesh
using $N_{hyp}$ points per space dimension, the computational cost of solving
Eq. (10) with $M$ initial data in Eq. (63) and error tolerance of
$\epsilon_{\text{hyp}}$ is $O(MTd^{d+3}(1/\epsilon_{\text{hyp}})^{d+1})$.

Finally we point out that for Hamilton-Jacobi equations and general nonlinear
hyperbolic PDEs, not all initial data lead to caustics or shocks, hence the
solutions may remain smooth for all time, thus multivalued-solutions will not
appear. In such cases, the solutions computed are exactly the smooth solutions
to the original nonlinear PDEs without needing to use the notion of viscosity
or multi-valued solutions.

###  II.2 System of nonlinear ODEs

A system of $D$ nonlinear ODEs subject to $M$ different initial data can be
written as

$\displaystyle\frac{dX^{[k]}(t)}{dt}=F(X^{[k]}(t)),\qquad
X^{[k]}\in\mathbb{R}^{D}\,,$ (12) $\displaystyle X^{[k]}(0)=X_{0}^{[k]},\qquad
k=1,\cdots,M$ (13)

where the same ODE is satisfied for each $k$. This can be interpreted as a
system of $M$ non-interacting particles in $D$ dimensions, each with a
trajectory described by $X^{[k]}(t)$. For nonlinear ODEs,
$F(X)\in\mathbb{R}^{D}$ are nonlinear functions of its argument $X$. To
compute ensemble averages

$\displaystyle\langle A(t)\rangle=\frac{1}{M}\sum_{k=1}^{M}A(X^{[k]}(t))$ (14)

one can directly solve for the system of ODEs starting from $M$ different
initial data $X_{0}^{[k]}$, and then carry out the summation in Eq. (14).

######  Lemma 3.

The classical computational cost in computing the ensemble observable $\langle
A(T)\rangle$ for $M$ initial data, to precision $\epsilon_{\text{ODE}}$ at
time $T=N_{t,\text{ODE}}\Delta t_{\text{ODE}}$ for time-step size $\Delta
t_{ODE}$ is $O({D^{3}}MT/\epsilon_{\text{ODE}})$.

###### Proof.

Assume the evaluation of each component of $F(X)$ costs, at most, $O(D)$
operations. Then for the first order method (say the forward Euler method),
the error is of $O(D\Delta t_{\text{ODE}})$. To reach an error
$\epsilon_{\text{ODE}}$ one needs $\Delta
t_{\text{ODE}}=O(\epsilon_{\text{ODE}}/D)$. The ensemble average step in Eq.
(14) costs $O(M)$ since it is just a sum of $M$ terms. The computational cost
of solving Eq. (12) and computing the ensemble average in Eq. (14) is
$O(D^{2}MN_{t,{\text{ODE}}})=O({D^{3}}MT/\epsilon_{\text{ODE}})$. ∎

Remark: The ODE solver, by the spectral deferred correction method [34], can
achieve the complexity of
$O(DMm\log(m)/\epsilon_{\text{ODE}}^{1/m})$ for any pre-chosen positive
integer $m$. So if one wants to choose a large $m$, the cost is about
$O(DMm\log(m))$. For simplicity, we won’t consider spectral methods in this
paper, but they can be the investigation of future work.

###  II.3 Quantum subroutines

We now briefly review the quantum linear systems problem (QLSP), which are
useful for solving ODEs and PDEs. It should be emphasized that the output of
QLSP for ODEs and PDEs are quantum states and not the classical solutions of
the ODEs and PDEs, thus making these quantum subroutines rather than full
quantum algorithms. We then describe the system of linear equations problem
(SLEP) which allows one to compute observables at the output of QLSP, thus
allowing the quantum algorithm to solve the same problem as the classical
algorithm.

####  II.3.1 Quantum linear systems problem: QLSP

The quantum linear systems problem (QLSP) [35, 36] can be stated informally in
the following way.

######  Problem 1.

(QLSP) Let $\mathcal{M}$ be a $2^{m}\times 2^{m}$ Hermitian matrix such that
$\|\mathcal{M}\|\leq 1$. Assume vectors $x$ and $y$ with elements
$\\{x_{i}\\}$, $\\{y_{i}\\}$ that satisfy $\mathcal{M}x=y$. One can then
define the following $m$-qubit quantum states
$|x\rangle\equiv\sum_{i}x_{i}/N_{x}|i\rangle$,
$|y\rangle\equiv\sum_{i}y_{i}/N_{y}|i\rangle$ where
$N_{x}=\sqrt{\sum_{i}|x_{i}|^{2}}$, $N_{y}=\sqrt{\sum_{i}|y_{i}|^{2}}$ are
normalisation constants. The aim of any QLSP algorithm is, when given access
to $\mathcal{M}$ and unitary $U_{initial}$ (where
$U_{initial}|0\rangle=|y\rangle$), to prepare the quantum state
$|x^{\prime}\rangle$ that is $\eta$-close to $|x\rangle$, i.e.,
$\||x^{\prime}\rangle-|x\rangle\|\leq\eta$.

The most notable algorithms to solve QLSP are the HHL algorithm [37] built on
quantum phase estimation and the alternative CKS algorithm [35] that bypasses
quantum phase estimation, where the latter can provide an exponential
improvement in precision. Their output is an approximation to the quantum
state $|x\rangle$ instead of an approximation to the solution of the original
vector problem $x=\mathcal{M}^{-1}y$. Thus these are often termed quantum
subroutines instead of full quantum algorithms, since they provide a stepping
stone but do not solve the same problem as the corresponding classical
algorithm.

In Problem 1, the informal phrase ‘when given access to $\mathcal{M}$’ in
order to solve QLSP, refers to a description of how the entries of
$\mathcal{M}$ can be accessed during the protocol. The total cost of the
protocol would also be computed with respect to the type of access one has.
For both algorithms [37, 35] the assumption is of sparse access to
$\mathcal{M}$, defined in the following way [38, 36].

######  Definition 4.

Sparse access to a Hermitian matrix $\mathcal{M}$ is a $4$-tuple
$(s,\|\mathcal{M}\|_{max},O_{\mathcal{M}},O_{F})$ and the $(i,j)^{\text{th}}$
entry of $\mathcal{M}$ is denoted $\mathcal{M}_{ij}$. Here $s$ is the sparsity
of $\mathcal{M}$ and $\|\mathcal{M}\|_{max}=\max_{i,j}(|\mathcal{M}_{ij}|)$ is
the max-norm of $\mathcal{M}$. $O_{M}$ and $O_{F}$ are unitary black boxes
which can access the matrix elements $\mathcal{M}_{ij}$ such that

$\displaystyle
O_{M}|j\rangle|k\rangle|z\rangle=|j\rangle|k\rangle|z\oplus\mathcal{M}_{jk}\rangle$
$\displaystyle O_{F}|j\rangle|l\rangle=|j\rangle|F(j,l)\rangle$ (15)

where the function $F$ takes the row index $j$ and a number $l=1,2,...,s$ and
outputs the column index of the $l^{\text{th}}$ non-zero elements in row $j$.

Then the HHL algorithm [37] has the following query complexity, which denotes
the number of times oracles $O_{M}$, $O_{F}$ and $U_{initial}$ are used
throughout the protocol. The gate complexity refers to the number of $2$-qubit
gates required in the algorithm.

######  Lemma 5.

[37, 35] Let $\mathcal{M}$ be Hermitian and each copy of $|y\rangle$ is
provided by the oracle $U_{initial}$. To create a state that is $\eta$-close
to $|x\rangle$, it is sufficient that the oracles $O_{M},O_{F}$ are queried
$O((s\kappa^{2}\|\mathcal{M}\|_{max}/\eta)\text{poly}\log(s\kappa\|\mathcal{M}\|_{max}/\eta))$
times, where $\kappa$ is the condition number of $\mathcal{M}$. The number of
2-qubit gates required in this algorithm is at most logarithmically larger
than the query complexities of $O_{M},O_{F}$. The query complexity for the
oracle $U_{initial}$ is
$\mathcal{O}(s\kappa\|\mathcal{M}\|_{max}\text{poly}\log(s\kappa\|\mathcal{M}\|_{max}/\eta))$.

The $1/\eta$ dependence in the query complexities of $O_{M},O_{F}$ comes from
using quantum phase estimation. To improve upon this factor, the CKS algorithm
in [35] bypasses phase estimation and instead uses a sequence of unitaries
whose sum approximates $\mathcal{M}^{-1}$.

######  Lemma 6.

[35] A state $\eta$-close to $|x\rangle$ can be created by querying the
$O_{M},O_{F}$ and $U_{initial}$ oracles
$\mathcal{O}(s\kappa\|\mathcal{M}\|_{max}\text{poly}\log(s\kappa\|\mathcal{M}\|_{max}/\eta))$
times, where the number of 2-qubit gates are at most logarithmally larger than
the query complexity.

An important application of QLSP is preparing quantum states whose amplitudes
are proportional to the solutions for linear ODEs and PDEs, known as
amplitude-encoding of the solutions. The first step is to discretise the ODEs
and PDEs and transform the equations into a linear algebra problem of the form
$x=\mathcal{M}^{-1}y$, where the size of $\mathcal{M}$ can be very large. The
quantum subroutines for QLSP is then in the matrix inversion process to
prepare the corresponding $|x\rangle$ states. The exact query and gate
complexities with respect to error $\eta$, time $T$ and the dimension $d$ of
the problem would also depend on details of the discretisation procedure.

We give full details of our discretisation procedures and the corresponding
matrix inversion problem in the main body of the paper.

####  II.3.2 System of linear equations problem: SLEP

When solving ODEs and PDEs, the actual desired outcomes of the problem are the
observables associated with the solutions of the ODEs and PDEs. Solving QLSP
only prepares $|x\rangle$, whereas the system of linear equations problem
(SLEP) [36] aims to compute observables from $x$. Then given the same
definitions as Problem 1, one can state SLEP in the following way.

######  Problem 2.

(SLEP) Given a Hermitian matrix $\mathcal{G}$, which is of the same size as
$\mathcal{M}$, access to $\mathcal{M}$ and $U_{initial}$, the aim of SLEP is
to compute the expectation value
$(\mathcal{M}^{-1}y)^{T}\mathcal{G}(\mathcal{M}^{-1}y)=x^{T}\mathcal{G}x$ to
precision $\epsilon^{\prime}$.

There are various different methods of measuring the outcome
$x^{T}\mathcal{G}x$ directly after obtaining
$|x\rangle\propto\mathcal{M}^{-1}|y\rangle$ from the output of QLSP. Many of
these methods have a query complexity with error scaling as
$1/(\epsilon^{\prime})^{2}$ [37], for instance the quantum swap test [39] when
$\mathcal{G}$ is a density matrix or applying the Hadamard test when
$\mathcal{G}$ is decomposed as a sum of two unitary operators [40]. An elegant
formalism that generalises the quantum methods for matrix inversion (that uses
sparse access to $\mathcal{M}$), as well as neatly achieving the improved
optimal $1/\epsilon^{\prime}$ scaling via amplitude estimation [41], without
too many extra assumptions, is the formalism of block access [42, 43, 36],
defined below. It is possible to create block access to $\mathcal{M}$ from
sparse access to $\mathcal{M}$ (see Lemma 32), a fact which we later exploit
to approximate $x^{T}\mathcal{G}x$ while still beginning from sparse access to
$\mathcal{M}$.

######  Definition 7.

Let $\mathcal{M}$ be a $m$-qubit Hermitian matrix, $\delta_{\mathcal{M}}>0$
and $n_{\mathcal{M}}$ is a positive integer. A $(m+n_{\mathcal{M}})$-qubit
unitary matrix $U_{\mathcal{M}}$ is a
$(\alpha_{\mathcal{M}},n_{\mathcal{M}},\delta_{\mathcal{M}})$-block encoding
of $\mathcal{M}$ if

$\displaystyle\|\mathcal{M}-\alpha_{\mathcal{M}}\langle
0^{n_{\mathcal{M}}}|U_{\mathcal{M}}|0^{n_{\mathcal{M}}}\rangle\|\leq\delta_{\mathcal{M}}.$
(16)

Block access to $\mathcal{M}$ is then the 4-tuple
$(\alpha_{\mathcal{M}},n_{\mathcal{M}},\delta_{\mathcal{M}},U_{\mathcal{M}})$
where $U_{\mathcal{M}}$ is the unitary black-box block-encoding of
$\mathcal{M}$.

In the rest of the paper, we assume that if the block access
$(\alpha_{\mathcal{M}},n_{\mathcal{M}},\delta_{\mathcal{M}},U_{\mathcal{M}})$
to $\mathcal{M}$ is given, then $U^{\dagger}_{\mathcal{M}}$,
controlled-$U_{\mathcal{M}}$ and controlled-$U^{\dagger}_{\mathcal{M}}$ are
also given.

Then in Problem 2, if one assumes that access to $\mathcal{M}$ refers to block
access to $\mathcal{M}$ instead of sparse access, there exists an algorithm
[36] that solves SLEP with the following query and gate complexities.

######  Lemma 8.

[36] A quantum algorithm can be constructed that takes $m$ qubits, block
access $(\alpha_{\mathcal{M}},n_{\mathcal{M}},0,U_{\mathcal{M}})$ to a
$2^{m}\times 2^{m}$ invertible Hermitian matrix $\mathcal{M}$ with condition
number $\kappa$, block access
$(\alpha_{\mathcal{G}},n_{\mathcal{G}},0,U_{\mathcal{G}})$ to a $2^{m}\times
2^{m}$ Hermitian matrix $\mathcal{G}$, an accuracy
$\epsilon\in[\alpha_{\mathcal{G}}/2^{m},\alpha_{\mathcal{G}}]$, and a
$m$-qubit unitary black box $U_{initial}$, and returns with probability at
least $2/3$ an $\epsilon^{\prime}$-additive approximation to
$x^{T}\mathcal{G}x$, where $x=\mathcal{M}^{-1}y$ and $y$ is the
$2^{m}$-dimensional vector with entries $y_{i}=\langle
i|U_{initial}|0\rangle=\langle i|y\rangle$, by making
$\mathcal{O}(\alpha_{\mathcal{G}}\kappa^{2}/\epsilon^{\prime})$ queries to
$U_{\mathcal{G}}$,
$\mathcal{O}(\alpha_{\mathcal{G}}\kappa^{3}\log(\alpha_{\mathcal{G}}\kappa^{2}/\epsilon)/\epsilon^{\prime})$
queries to $U_{\mathcal{M}}$,
$\mathcal{O}(\alpha_{\mathcal{G}}\kappa^{2}/\epsilon^{\prime})$ queries to
$U_{initial}$ and
$\mathcal{O}(\alpha_{\mathcal{G}}\kappa^{2}(m+n_{\mathcal{G}}+n_{\mathcal{M}}+n_{\mathcal{M}}\kappa\log(\alpha_{\mathcal{G}}\kappa^{2}/\epsilon^{\prime}))/\epsilon^{\prime})$
additional $2$-qubit gates.

We later use a modified version of this result, when given sparse access to
$\mathcal{M}$, for computing ensemble averages from nonlinear ODEs and PDEs.

##  III Solving Hamilton-Jacobi equations

In this section, we develop a quantum algorithm for computing observables from
Hamilton-Jacobi equations. The first step is to transform the nonlinear
equation into a linear equation without making any approximations. One
convenient method is to use the level set formulation [25]. We demonstrate
below (Section III.1) how to use this technique to translate a
$(d+1)$-dimensional nonlinear Hamilton-Jacobi equations into
$(2d+1)$-dimensional (linear) Liouville equations. By discretising this new
PDE, we can convert this into a linear algebra problem. We also present
physical observables that can be computed in our formalism in Section III.2
and present explicit error bounds. Finally in Section III.3 we show how to
compute the observable with a quantum algorithm and present its query and gate
complexities.

###  III.1 Linear represenation of the nonlinear Hamilton-Jacobi equation

Here we show how to transform a nonlinear Hamilton-Jacobi equation into a
linear algebra problem by a two-step procedure: using a level set function to
convert the nonlinear equation into a linear equation and then to discretise
that linear equation.

The level set function $\phi^{[k]}_{i}(t,x,p)$ can be defined by

$\phi^{[k]}_{i}(t,x,p=u^{[k]}(t,x))=0$ (17)

where $i=1,\cdots,d$ and $\,x,p\in\mathbb{R}^{d}$, $k=1,...,M$ and
$\\{u^{[k]}(t,x)\\}$ are the solutions of the nonlinear Hamilton-Jacobi
equation in gradient form with $M$ initial conditions in Eq. (6). The zero
level set of $\phi^{[k]}_{i}$ is the set
$\\{(t,x,p)|\phi^{[k]}_{i}(t,x,p)=0\\}$.

Since $u^{[k]}(t,x)$ solves Eq. (6), then one can show that
$\phi^{[k]}=(\phi^{[k]}_{1},\cdots,\phi^{[k]}_{d})\in\mathbb{R}^{d}$ solves a
(linear!) Liouville equation [25]

$\partial_{t}\phi^{[k]}+\nabla_{p}H\cdot\nabla_{x}\phi^{[k]}-\nabla_{x}H\cdot\nabla_{p}\phi^{[k]}=0.$
(18)

Note that the (bi)-characteristics of the Liouville equation in Eq. (18) is
the Hamiltonian system in Eq. (5). The initial data can be chosen as

$\phi^{[k]}_{i}(0,x,p)=p_{i}-u^{[k]}_{i}(0,x),\quad i=1,\cdots,d.$ (19)

Then $u^{[k]}$ can be recovered from the intersection of the zero level sets
of $\phi_{i}\,(i=1,\cdots,d)$, as in (21), namely

$u^{[k]}(t,x)=\\{p(t,x)|\,\phi_{i}^{[k]}(t,x,p)=0,\,i=1,\cdots,d\\}.$ (20)

Note that $\phi_{i}^{[k]}(t,x,p)=0$ may have multiple (say $J_{k}$) roots,
denoted by $p_{\gamma(t,x)}(\gamma=1,\cdots,J_{k})$, hence the so-called
multi-valued solutions will arise, which is denoted by

$u^{[k]}_{\gamma}(t,x)=\cap_{i=1}^{d}\\{p_{\gamma}(t,x)|\phi^{[k]}_{i}(t,x,p_{\gamma})=0\\},\quad\gamma=1,\cdots,J_{k}.$
(21)

We include this possibility when defining our observable in Section III.2.

It is crucial to observe here that we have now transformed a
$(d+1)$-dimensional nonlinear PDE– the Hamilton-Jacobi equation– to a
$(2d+1)$-dimensional linear PDE –the Liouville equation–without any
assumptions on either the form or extent of the original nonlinearity. No
linear approximation is made. The mapping is exact. The cost in going from a
nonlinear to a linear system is only at the expense of doubling the dimension.
Doubling the dimension of the problem may seem too costly for a classical
device because the cost increases exponentially with dimension. However, we
will see that for quantum algorithms, the relative overhead in doubling the
dimension can be up to exponentially smaller.

While in principle one could apply quantum subroutines for QLSP to Eq. (18)
and generate the quantum states $|\phi\rangle$ whose amplitudes are
proportional to $\phi_{i}$, it would be too costly to recover our desired
solutions $u$ or observables from $|\phi\rangle$ due to the extra measurement
costs in finding the zero level set of $\phi_{i}^{[k]}\,(i=1,\cdots d)$. An
alternative method is to solve for $\psi$, defined by the following problem

$\displaystyle\partial_{t}\psi+\nabla_{p}H\cdot\nabla_{x}\psi-\nabla_{x}H\cdot\nabla_{p}\psi=0$
(22)

with the initial condition

$\displaystyle\psi(0,x,p)=\frac{1}{M}\sum_{k=1}^{M}\prod_{i=1}^{d}\delta(p_{i}-u^{[k]}_{i}(0,x)).$
(23)

Then we have the following result.

######  Lemma 9.

The analytical solution of (22) with initial data (23) is

$\displaystyle\psi(t,x,p)=\frac{1}{M}\sum_{k=1}^{M}\delta(\phi^{[k]}(t,x,p))\,.$
(24)

###### Proof.

See Appendix A ∎

An important observation here is that all $M$ distinct initial conditions of
$u^{[k]}$ in the original problem have now been converted into a single
initial condition in $\psi$. We will exploit this property later to show that
the resource cost for the quantum algorithm is independent of $M$.

Remark: Instead of (23), one can define a more general average

$\displaystyle\psi(0,x,p)=\frac{1}{M}\sum_{k=1}^{M}c_{k}\prod_{i=1}^{d}\delta(p_{i}-u^{[k]}_{i}(0,x)).$
(25)

where $0<c_{k}<1$ and $\sum_{k}c_{k}=1$. These weights will persist to all
later times, and one has

$\displaystyle\psi(t,x,p)=\frac{1}{M}\sum_{k=1}^{M}c_{k}\delta(\phi^{[k]}(t,x,p))$
(26)

and all physical observables, which are moments of $\psi$, are weighted
average –with the same weights– of individual corresponding observables.

One can convert the linear PDE in Eq. (22) with initial conditions in Eq. (23)
into a linear algebra problem by discretising the function

$\displaystyle\psi(t,x,p)\rightarrow\psi^{\omega}_{n,\boldsymbol{j},\boldsymbol{l}}$
(27)

where $n$ denotes the time step from $t\rightarrow t_{n}\equiv n\Delta t$ and
$n=0,...,N_{t}$. The vectors $\boldsymbol{j}=(j_{1},...,j_{d})$,
$\mathbf{l}=(l_{1},...,l_{d})$ denote spatial grid indices with grid size
$h=1/N$ where $x\rightarrow h\boldsymbol{j}$, $p\rightarrow h\boldsymbol{l}$
are grid points with $j_{i},l_{i}=1,..,N$ for $i=1,..,d$.

Throughout the paper we always assume, without loss of generality, the
computational domain to be in a box of $[0,1]^{D}$, where $D$ is the spatial
dimension of the problem. Since all the linear PDEs considered in this paper
are transport equations, due to their finite propagation speeds, by time
$T<\infty$, as long as the initial data have compact support in $x$–which we
assume here for $u_{0}(x)$– solution will still have compact support for
$t\leq T$ in $x$ domain. As far as the domain in $p$ is concerned, notice that
we always start with delta functions in $p$, which has a compact support as
long as $u(t,x)$ is bounded, so the computational domain in $p$ can also be
restricted to a finite domain as long as $T<\infty$. Upon a suitable scaling
we confine our computational domain, in both $x$ and $p$, within the box
$[0,1]^{D}$.

Since the initial condition in Eq. (23) involves a delta function, one also
needs a discretised delta function $\delta_{\omega}$ where $\omega$ is a
smoothing parameter of the delta function. Hence for the initial state at
$n=0$ we have

$\displaystyle\psi^{\omega}_{0,\boldsymbol{j},\boldsymbol{l}}=\frac{1}{M}\sum_{k=1}^{M}\prod_{i=1}^{d}\delta_{\omega}(l_{i}h-u^{[k]}_{i}(n=0,h\boldsymbol{j})).$
(28)

As conventionally done, we choose $\delta_{\omega}$ to be smooth and to
satisfy, for $x\in\mathbb{R}^{1}$,

$\delta_{\omega}(x)=0\quad{\text{if}}\quad|x|>\omega;\quad\int_{|x|\leq\omega}\delta_{\omega}(x)\,dx=1.$
(29)

One usually approximates $\delta_{\omega}$ by the form

$\delta_{\omega}(x)=\begin{cases}\frac{1}{\omega}\beta(x/\omega)\quad&|x|\leq\omega;\\\
0\quad&|x|>\omega\end{cases}$ (30)

where typical choices of $\beta(x)$ include $\beta(x)=1-|\beta|$ and
$\beta(x)=\frac{1}{2}(1+\cos(\pi x))$ [44]. Here one can choose $\omega=mh$
where $m$ is the number of mesh points within the support of
$\delta_{\omega}$. For $x=(x_{1},\cdots,x_{d})\in\mathbb{R}^{d}$, one defines
$\delta_{\omega}(x)\equiv\Pi_{i=1}^{d}\delta_{\omega}(x_{i})$.

The solution of the discretised PDE can then be written as the following
matrix equation

$\displaystyle\begin{pmatrix}\psi_{1,\boldsymbol{j},\boldsymbol{l}}\\\
\psi_{2,\boldsymbol{j},\boldsymbol{l}}\\\ \vdots\\\
\psi_{N_{t}-1,\boldsymbol{j},\boldsymbol{l}}\\\
\psi_{N_{t},\boldsymbol{j},\boldsymbol{l}}\end{pmatrix}=\mathcal{K}^{-1}\begin{pmatrix}\psi_{0,\boldsymbol{j},\boldsymbol{l}}\\\
0\\\ \vdots\\\ 0\\\ 0\end{pmatrix}$ (31)

where $\mathcal{K}$ is a $N_{t}N^{2d}\times N_{t}N^{2d}$ Toeplitz matrix. See
Appendix B for the form of $\mathcal{K}$ and the details of the discretisation
procedure.

Both QLSP and SLEP involve the inversion of an Hermitian matrix $\mathcal{M}$.
Since $\mathcal{K}$ is not Hermitian, one can define a new Hermitian matrix

$\displaystyle\mathcal{M}=\begin{pmatrix}0&\mathcal{K}\\\
\mathcal{K}^{\dagger}&0\end{pmatrix}$ (32)

which has the same sparsity and condition number as $\mathcal{K}$. Using
$\mathcal{M}$, the matrix inversion problem to solve
$\psi_{n,\boldsymbol{j},\boldsymbol{l}}$ becomes

$\displaystyle\begin{pmatrix}\boldsymbol{0}\\\
\psi_{n,\boldsymbol{j},\boldsymbol{l}}\end{pmatrix}=\mathcal{M}^{-1}\begin{pmatrix}\psi_{0,\boldsymbol{j},\boldsymbol{l}}\\\
\boldsymbol{0}\end{pmatrix}$ (33)

where $\boldsymbol{0}$ is a zero-vector of the same dimension as
$\mathcal{K}$.

######  Lemma 10.

The condition number of $\mathcal{M}$ is $\kappa\leq O(dNT)$ where
$T=N_{t}\Delta t$ is the stopping time, and sparsity is $s=O(d)$.

###### Proof.

See Appendix C. ∎

###  III.2 The observables

Given any function $G:\mathbb{R}^{d}\rightarrow\mathbb{R}$ one can define the
following ensemble average which we call the observable.

######  Definition 11.

The following ensemble average we define as the observable

$\langle
G(t,x)\rangle\equiv\int_{\mathbb{R}^{d}}G(p)\psi(t,x,p)dp=\frac{1}{M}\sum_{k=1}^{M}\int_{\mathbb{R}^{d}}G(p)\delta(\phi^{[k]}(t,x,p))dp=\frac{1}{M}\sum_{k=1}^{M}\sum_{\gamma=1}^{J_{k}}\frac{G(u^{[k]}_{\gamma}(t,x))}{\mathcal{J}^{[k]}_{\gamma}},$
(34)

where the Jacobian
$\mathcal{J}^{[k]}_{\gamma}\equiv|\det(\partial{\phi^{[k]}}/\partial
p)|_{p=u^{[k]}_{\gamma}(t,x)}$.

This is an ensemble average of solution over the $M$ different initial data,
each with the multi-valued solutions of $J_{k}$ branches. The multi-valued
solution is a weighted average of each branch with weights
$1/\mathcal{J}^{[k]}_{\gamma}$ that depends on the level set function
$\phi^{[k]}$.

Remark: This is easily extended to the case of more general weighting of the
initial conditions in Eq. (26).

Assume $\psi$ is the solution to Eq. (22), and $u^{[k]}$ is the solution to
Eq. (6). One can then compute $\langle G(x,t)\rangle$ by using the numerical
quadrature rule

$\displaystyle\langle
G(t_{n},x=\boldsymbol{j}/N)\rangle=\int_{\mathbb{R}^{d}}G(p)\psi(t_{n},\boldsymbol{j}/N,p)dp\approx\frac{1}{N^{d}}\sum_{\boldsymbol{l}}^{N}G_{\boldsymbol{l}}\psi^{\omega}_{n,\boldsymbol{j},\boldsymbol{l}}\equiv\langle
G^{\omega}_{n,\boldsymbol{j}}\rangle$ (35)

where
$\psi^{\omega}_{n,\boldsymbol{j},\boldsymbol{l}}\equiv\frac{1}{M}\sum_{k=1}^{M}\delta_{\omega}(\phi^{[k]}(t_{n},\boldsymbol{j}/N,\boldsymbol{l}/N))$
and $G_{\boldsymbol{l}}\equiv G(\boldsymbol{l}/N)$. We also use the notation
$\sum_{\boldsymbol{l}}^{N}\equiv\sum_{l_{1}=1}^{N}...\sum_{l_{d}=1}^{N}$.

We can also define an ensemble average normalised by its zeroth moment:

$\displaystyle G_{O}(t,x)\equiv\frac{\langle
G(t,x)\rangle}{\langle\mathbf{1}\rangle}=\frac{\int_{\mathbb{R}^{n}}G(p)\psi(t,x,p)dp}{\int_{\mathbb{R}^{n}}\psi(t,x,p)dp}$
(36)

where $\mathbf{1}$ is the identity function, corresponding to the zeroth
moment of $\psi$. In the cases of $M=1$ and when there are no multi-valued
solutions, one has

$\displaystyle G_{O}(t,x)=G(u(t,x)).$ (37)

Our goal is to devise a quantum algorithm to compute the observable $\langle
G(t,x)\rangle$ (also $G_{O}(t,x)$) to precision $\epsilon$ and demonstrate
that it can be more efficient on a quantum device with respect to parameters
$M$, $d$ and $\epsilon$ compared to a purely classical algorithm. In the next
subsection we show how different $G$ leads to different physical observables
captured by $\langle G(t,x)\rangle$ and $G_{O}$, such as density, momentum and
energy.

####  III.2.1 Physical interpretation

We now give physical interpretations for the observables defined in Eq. (34)
through several physically important examples. We show in these cases that for
nonlinear Hamilton-Jacobi where the WKB approximation is applicable, the
corresponding Wigner function, in the semiclassical limit of the Schrödinger
equation, obeys exactly the same Liouville equation satisfied by
$\psi(t,x,p)$. Given different choices of $G(p)$, the observable $\langle
G(t,x)\rangle$ correspond to moments of this Wigner function and they are
exactly the physical observables like density, momentum and energy.

We first consider the classical limit (WKB approximation) of the Schrödinger
equation with $M=1$:

$\displaystyle
i\hbar\partial_{t}\Psi=-\frac{\hbar^{2}}{2}\Delta\Psi+V(x)\Psi\,,\qquad\Psi(0,x)=A_{0}(x)e^{i\frac{S_{0}(x)}{\hbar}},$
(38)

with wavefunction $\Psi(t,x)$. The WKB analysis uses the ansatz
$\Psi(t,x)=A(t,x)\exp(iS(t,x)/\hbar)$, where $A(t,x)$ and $S(t,x)$ are the
amplitude and phase respectively. Ignoring $O(\hbar^{2})$ terms, this ansatz
results in two independent PDEs: the eikonal equation for $S$, which is a
Hamilton-Jacobi PDE in Eq. (3) with Hamiltonian $H(x,\nabla S)=(1/2)|\nabla
S|^{2}+V(x)$, and the transport equation for $|A|^{2}$ . These two equations
can be deduced from the moment-closure of the Liouville equation

$\displaystyle\partial_{t}w+\nabla_{p}H\cdot\nabla_{x}w-\nabla_{x}H\cdot\nabla_{p}w=0$
(39)

with initial data

$w(0,x,p)=|A_{0}(x)|^{2}\delta(p-\nabla S_{0}(x)),$ (40)

with mono-kinetic ansatz $w(t,x,p)=|A(t,x)|^{2}\delta(p-\nabla_{x}S(t,x))$,
but are not valid beyond caustics since $\delta(p-\nabla_{x}S(t,x))$ is not
well-defined when $\nabla_{x}S(t,x)=u(t,x)$ becomes discontinuous. But EQ.
(39) is valid globally in time [45, 46], since it unfolds the caustics in the
phase space. The problem defined by Eqs. (39)-(40) can be solved by first
solving

$\mathcal{L}\tilde{\phi}=0\,,\qquad\tilde{\phi}(0,x,p)=|A_{0}(x)|^{2}\,,$ (41)

and

$\mathcal{L}\phi=0\,,\quad\phi(0,x,p)=p-\nabla S_{0}(x)\,,$ (42)

independently and then $w$ can be obtained by $w=\tilde{\phi}\delta(\phi)$
[29]. Here $\phi$ is exactly the level set function defined by (18) and (19),
while $\psi$ defined (22) and (23) is exactly $w$ in the case of $A_{0}(x)=1$
for $M=1$. Therefore different choices of $G$ can be used to recover moments
of the Wigner function. These moments provide the classical limits of the
observables of the original Schrödinger equation. For instance, the zeroth,
first and second moments of the Wigner function $w$ are

$\rho(t,x)=\int w\,dp\,,\quad\rho(t,x)u(t,x)=\int
pw\,dp\,,\quad\frac{1}{2}\rho(t,x)u^{2}(t,x)=\int\frac{|p|^{2}}{2}w\,dp\quad$
(43)

which correspond to using $G(p)=1,p,|p|^{2}$ respectively. Here $\rho(t,x)$ is
the classical limit to the position density $|\Psi|^{2}$, $\rho u$ is the
momentum, or the classical limit to the current density
$\hbar\,\text{Im}(\overline{\Psi}\nabla{\Psi})$, while the second moment is
the classical limit of the kinetic energy $(\hbar^{2}/2)|\nabla\Psi|^{2}$. The
total energy $E(t,x)$ can be recovered by combining the second and zeroth
moments to choose $G(p)=(1/2)|p|^{2}+V(x)$ to obtain
$E(t,x)=\int(|p|^{2}/2+V(x))wdp$. Likewise, the moments of $\psi$, if
similarly defined as in (43), give the same physical observables in the
special case of $A_{0}(x)=1.$

If one begins with more general $M>1$ initial data

$w(0,x,p)=\frac{1}{M}\sum_{k=1}^{M}|A^{[k]}_{0}(x)|^{2}\delta(p-\nabla
S^{[k]}_{0}(x))$ (44)

then $w(t,x,p)=\frac{1}{M}\sum_{k=1}^{M}w^{[k]}(t,x,p)$ where $w^{[k]}(t,x,p)$
solves equation (39) with initial condition
$w^{k]}(0,x,p)=|A^{[k]}_{0}(x)|^{2}\delta(p-\nabla S^{[k]}_{0}(x))$, thus
consequently, due to the linear superposition principle, the observables

$\int w\,dp=\frac{1}{M}\sum_{k=1}^{M}\rho^{[k]}(t,x)\,,\quad\int
pw\,dp=\frac{1}{M}\sum_{k=1}^{M}\rho^{[k]}(t,x)u^{[k]}(t,x)\,,\quad\int\left(\frac{p^{2}}{2}+V(x)\right)w\,dp=\frac{1}{M}\sum_{k=1}^{M}E^{[k]}(t,x)\quad$
(45)

which are the ensemble average of the physical observables (position density,
momentum and kinetic energy) of each individual observables.

We can also consider examples for general symmetric hyperbolic systems. The
high frequency limit of general symmetric hyperbolic systems–for which
geometric optics is one example– possesses a strong similarity with the semi-
classical limit of the Schrödinger equation via the Wigner transform [47].

In the example of acoustic waves, the WKB approximation of such systems gives
a Hamilton-Jacobi equation in Eq. (3) for the phase $S$ with Hamiltonian

$H(x,p)=\frac{1}{\sigma(x)\tau(x)}|p|$ (46)

where $\sigma(x)$ and $\tau(x)$ are the density and compressibility of the
wave. The WKB approximation also provides a transport equation for the
amplitude $|A|$. The corresponding Wigner function for this approximation also
satisfies the Liouville equation in Eq. (39) with the Hamiltonian in Eq. (46)
[46]. Then its zeroth and first moments are

$\rho(t,x)=\int w\,dp\,,\quad\rho(t,x)u(t,x)=\int pw\,dp,$ (47)

which can be used to obtain the intensity $|A|=\sqrt{\rho}$ and the velocity
$u=\rho(t,x)u(t,x)/\rho(t,x)=\nabla S$, where $\nabla S$ is also known as the
slowness vector.

The second example is Maxwell’s equation in an isotropic medium. One can
similarly retrieve the moments of its high frequency limit via a Wigner
analysis. The observables recovered include the energy and the Poynting vector
[47]. The third example is elastic waves, where a similar analysis can allow
one to retrieve the kinetic energy, strain energy and energy flux of the
elastic wave [47]. The fourth example is the Dirac equation, the relativistic
version of the Schrödinger equation that describes very fast electrons in an
electromagnetic field. In the semiclassical limit, it gives rise to a similar
Liouville equation, like in geometric optics, but with a Lorentz term due to
relativistic effects [46]. The Wigner function can similarly recover physical
observables such as the density of positrons and electrons. [46].

We can also explain the presence of the Jacobian weights
$\mathcal{J}^{[k]}_{\gamma}$ in Eq. (34). The Wigner function $w$ for $M=1$
has the general form $w(t,x,p)=a(t,x,p)\delta(\phi(t,x,p))$, where in general
$\phi(t,x,p)\neq p-u(t,x)$ for $t>0$. By definition
$\delta(\phi(t,x,p))=\sum_{\gamma=1}^{J_{k}}\delta(p-u_{\gamma}(t,x))/\mathcal{J}_{\gamma}$,
in the presence of multiple zeros of $\phi(t,x,p)$, corresponding to multi-
valued solutions. Thus one can rewrite the Wigner function as
$w(t,x,p)=\sum_{\gamma=1}^{J_{k}}\rho_{\gamma}(t,x,p)\delta(p-u_{\gamma}(t,x))$,
where $\rho_{\gamma}(t,x,p)=a(t,x,p)/\mathcal{J}_{\gamma}$ [46]. This means
$\rho_{\gamma}$ (with the Jacobian term contained in its definition) is the
true density associated with each multi-valued solution, since $\rho(t,x)=\int
wdp=\sum_{\gamma=1}^{J_{k}}\rho_{\gamma}$. The extension to $M>1$ and other
moments is straightforward. The fact that this non-trivial Jacobian term (that
is naturally captured by using $\psi(t,x,p)$) cannot be ignored in order to
obtain the correct physical quantities is important, an issue we will return
to at the end of Section III.3.

####  III.2.2 Classical error bounds

Suppose one is interested in computing observables coming from solutions to
Eq. (6) with $M$ different initial data. Several approximations will be used
in our computation. (1) We solve instead the corresponding linear PDE problem
in Eq. (22) with initial data in Eq. (23), where the PDE will be approximated
by some classical finite difference or finite volume method on a
$2d$-dimensional phase space mesh on $[0,1]^{2d}$ with mesh size $h=1/N$ and
forward Euler method in time with time step $\Delta t$; (2) The delta function
$\delta$ in the initial condition will be approximated by a discrete delta
function $\delta_{\omega}$ defined in Eq. (30), with smoothing parameter
$\omega=mh$, for some small integer $m>0$; (3) The observable defined by the
integral in Eq. (34) will be approximated by a quadrature rule in Eq. (35).
Taking all three sources of error into account, we obtain the following error
bound in estimating the ensemble average.

######  Lemma 12.

The ensemble average $\langle G(t_{n},\boldsymbol{j}/N)\rangle$ can be
estimated by $\langle G^{\omega}_{n,\boldsymbol{j}}\rangle$ with error

$\displaystyle\epsilon_{CL}\equiv|\langle
G(t_{n},\boldsymbol{j}/N)\rangle-\langle
G^{\omega}_{n,\boldsymbol{j}}\rangle|\leq C(\omega+dh/\omega^{2})$ (48)

where $C>0$ is independent of $h,\omega$. By choosing $\omega=(dh)^{1/3}$ one
gets

$\displaystyle\epsilon_{CL}\equiv|\langle
G(t_{n},\boldsymbol{j}/N)\rangle-\langle
G^{\omega}_{n,\boldsymbol{j}}\rangle|\leq C(dh)^{1/3}.$ (49)

###### Proof.

See Appendix D. ∎

Remark: By using a second order finite difference scheme one can improve the
above error to $O(\sqrt{dh})$ with $\omega=\sqrt{dh}$. Higher order finite
difference approximations will further reduce $\omega$.

This bound can also be used to estimate the total computational cost for a
classical computer to solve the $(2d+1)$-dimensional linear) Liouville PDE.

######  Lemma 13.

If one uses $N_{\text{CL}}$ points in each dimension of the $2d$ phase space
(hence $h_{\text{CL}}=1/N_{\text{CL}}$) to approximate Eq. (22) to time $T$ by
a first order finite difference or finite volume scheme, then
$N_{\text{CL}}=O(d/\epsilon_{\text{CL}}^{3})$ for
$\omega=(dh_{\text{CL}})^{1/3}$, while the overall computational cost will be
of $O(d^{2d+3}T(1/\epsilon_{\text{CL}})^{6d+3})$.

###### Proof.

If one uses $N_{\text{CL}}$ points in each dimension of the $2d$ phase space
(hence $h_{\text{CL}}=1/N_{\text{CL}}$) to approximate (22) to time $T$ by a
first order finite difference or finite volume scheme, the CFL condition will
require $\Delta t_{\text{CL}}=O(h_{\text{CL}}/(2d))$. To reach an error of
$O(\epsilon_{\text{CL}})$ one needs
$N_{\text{CL}}=O(d/(\omega^{2}(O(\epsilon_{\text{CL}})-O(\omega))))$. (If one
chooses $\omega=(dh_{\text{CL}})^{1/3}$ then
$h_{\text{CL}}=O(\epsilon_{\text{CL}}^{3}/d)$,
$N_{\text{CL}}=O(d^{1/3}/(h_{\text{CL}}^{2/3}\epsilon_{\text{CL}}))=O(d/\epsilon_{\text{CL}}^{3})$.)
The overall computational cost will be of
$O(dN_{t,\text{CL}}N_{\text{CL}}^{2d})=O(d^{2}TN_{\text{CL}}^{2d+1})=O(d^{2d+3}T(1/\epsilon_{\text{CL}})^{6d+3})$.
∎

Remark: The above computational cost is independent of $M$. However, the cost
now is much greater than that in Lemma 1, since one now solves a higher
dimensional PDE, and with much smaller ($\omega$-dependent) mesh size to
numerically resolve the discrete delta function $\delta_{\omega}$.

###  III.3 The quantum algorithm to approximate physical observables

Our aim is to approximate the observable $\langle
G(t_{n},\boldsymbol{j}/N)\rangle$ by devising the corresponding SLEP quantum
algorithm. Consider $\langle G^{\omega}_{n,\boldsymbol{j}}\rangle$ defined by

$\displaystyle\langle
G^{\omega}_{n,\boldsymbol{j}}\rangle\equiv\frac{1}{N^{d}}\sum_{\boldsymbol{l}}^{N}G_{\boldsymbol{l}}\psi_{n,\boldsymbol{j},\boldsymbol{l}}.$
(50)

The initial condition in its discretised form in Eq. (28) can be represented
by a quantum state

$\displaystyle|\psi_{0}\rangle=\frac{1}{N_{\psi_{0}}}\sum_{\boldsymbol{j}^{N}}\sum_{\boldsymbol{l}}^{N}\psi_{0,\boldsymbol{j},\boldsymbol{l}}|\boldsymbol{j}\rangle|\boldsymbol{l}\rangle|n=0\rangle$
(51)

and the normalisation is given by
$N_{\psi_{0}}=\sqrt{\sum_{\boldsymbol{j}}^{N}\sum_{\boldsymbol{l}}^{N}|\psi_{0,\boldsymbol{j},\boldsymbol{l}}|^{2}}$.
In this paper, we assume access to a unitary operation $U_{initial}$ that
prepares $U_{initial}|0\rangle=|\psi_{0}\rangle$. We also define the state

$\displaystyle|G_{n,\boldsymbol{j}}\rangle\equiv\frac{1}{N_{G}}\sum_{\boldsymbol{l}}^{N}G^{*}_{\boldsymbol{l}}|\boldsymbol{l}\rangle|\boldsymbol{j}\rangle|n\rangle$
(52)

where we note that here we do not sum over the time step index $n$ or the
spatial index $\boldsymbol{j}$. The normalisation is
$N_{G}=\sqrt{\sum_{\boldsymbol{l}}^{N}|G^{*}_{\boldsymbol{l}}|^{2}}$.

Using the states $|G_{n,\boldsymbol{j},\boldsymbol{l}}\rangle$ and applying
matrix inversion algorithms on $|\psi_{0}\rangle$, the observable can be
recovered. Given the density matrix
$\mathcal{G}\equiv|G_{n,\boldsymbol{j}}\rangle\langle G_{n,\boldsymbol{j}}|$,
one observes that the expectation value

$\displaystyle\Upsilon\equiv\langle\psi_{0}|(\mathcal{M}^{-1})^{\dagger}\mathcal{G}\mathcal{M}^{-1}|\psi_{0}\rangle=\frac{1}{N^{2}_{\psi_{0}}}\sum_{\boldsymbol{j}^{\prime},\boldsymbol{l}^{\prime}}(\mathcal{M}^{-1}\psi_{0,\boldsymbol{j}^{\prime},\boldsymbol{l}^{\prime}})^{T}\mathcal{G}(\mathcal{M}^{-1}\psi_{0,\boldsymbol{j}^{\prime},\boldsymbol{l}^{\prime}})$
$\displaystyle=\frac{1}{N^{2}_{\psi_{0}}N_{G}^{2}}\Big{|}\sum_{\boldsymbol{l}}G_{\boldsymbol{l}}\psi_{n,\boldsymbol{j},\boldsymbol{l}}\Big{|}^{2}=\frac{N^{2d}}{N^{2}_{\psi_{0}}N_{G}^{2}}\langle
G^{\omega}_{n,\boldsymbol{j}}\rangle^{2}$ (53)

where the last equality comes from Eq. (50). Our aim is to use a quantum
algorithm to extract $\Upsilon$, from which we can approximate the ensemble
average

$\displaystyle\langle G(t_{n},\boldsymbol{j}/N)\rangle\approx\langle
G^{\omega}_{n,\boldsymbol{j}}\rangle\equiv\frac{1}{N^{d}}\sum_{\boldsymbol{l}=1}^{N}G_{\boldsymbol{l}}\psi_{n,\boldsymbol{j},\boldsymbol{l}}=\frac{1}{N^{d}}N_{\psi_{0}}N_{G}|\sqrt{\Upsilon}|=n_{\psi_{0}}n_{G}|\sqrt{\Upsilon}|,$
(54)

where we define $n_{G}\equiv N_{G}/N^{d/2}$ and $n_{\psi_{0}}\equiv
N_{\psi_{0}}/N^{d/2}$. Then we have the following lemma.

######  Lemma 14.

The constant $n_{G}\equiv N_{G}/N^{d/2}=O(1)$ and depending on the
application, the range of $n_{\psi_{0}}$ lies in $O(1)\leq n_{\psi_{0}}\equiv
N_{\psi_{0}}/N^{d/2}\leq O(N^{d/2})$. Different $n_{\psi_{0}}$ corresponds to
different initial data. If we assume the initial data has support in a box of
size $\beta$, then $n_{\psi_{0}}=O((\beta\sqrt{N})^{d})$.

###### Proof.

See Appendix E. ∎

Remark: The definition $n_{\psi_{0}}=N_{\psi_{0}}/N^{d/2}$ can be an
underestimate in situations where the $N^{d}$ term in Eq. (54) overestimates
the number of non-zero factors in the summation. For instance, for a point
source where $\beta\sim 1/N$, $n_{\psi_{0}}<1$. To deal with similar
scenarios, more information about the problem is required to improve the
quadrature rule approximation in Eq. (54).

However, since the quantum device can only output an approximation
$\tilde{\Upsilon}$ of $\Upsilon$ to finite precision, we can only access the
estimate

$\displaystyle\langle\tilde{G}^{\omega}_{n,\boldsymbol{j}}\rangle\equiv
n_{\psi_{0}}n_{G}|\sqrt{\tilde{\Upsilon}}|$ (55)

where $|\sqrt{\tilde{\Upsilon}}-\sqrt{\Upsilon}|\leq\epsilon_{G}$. Then we
have the following contributions to the total error in estimating $\langle
G(t_{n},\boldsymbol{j}/N))\rangle$.

######  Lemma 15.

The error $|\langle
G(t_{n},\boldsymbol{j}/N)\rangle-\langle\tilde{G}^{\omega}_{n,\boldsymbol{j}}\rangle|$
can be broken into two independent sources of error

$\displaystyle|\langle
G(t_{n},\boldsymbol{j}/N)\rangle-\langle\tilde{G}^{\omega}_{n,\boldsymbol{j}}\rangle|\leq|\langle
G(t_{n},\boldsymbol{j}/N)\rangle-\langle
G^{\omega}_{n,\boldsymbol{j}}\rangle|+|\langle
G^{\omega}_{n,\boldsymbol{j}}\rangle-\langle\tilde{G}^{\omega}_{n,\boldsymbol{j}}\rangle|$
$\displaystyle=\epsilon_{\text{CL}}+\epsilon_{Q}\leq\epsilon.$ (56)

Here $\epsilon_{\text{CL}}$ comes from approximating the solutions of the
discretised linear PDE (Liouville equation) corresponding to the original
Liouville equation, as estimated in Lemma 12. The contributions to the quantum
sources of error is denoted $\epsilon_{Q}$. If
$|\sqrt{\tilde{\Upsilon}}-\sqrt{\Upsilon}|\leq\epsilon_{G}$ and we impose
$\epsilon_{\text{CL}}\sim\epsilon_{Q}$, then it is sufficient to choose
$\epsilon_{G}\sim\epsilon/(n_{G}n_{\psi_{0}})\sim\epsilon_{\text{CL}}/(n_{G}n_{\psi_{0}})$
where $n_{G}=O(1)$, $n_{\psi_{0}}\geq O(1)$.

###### Proof.

The quantum error can be written as

$\displaystyle\epsilon_{Q}=|\langle
G^{\omega}_{n,\boldsymbol{j}}\rangle-\langle\tilde{G}^{\omega}_{n,\boldsymbol{j}}\rangle|=n_{G}n_{\psi_{0}}(|\sqrt{\Upsilon}|-|\sqrt{\tilde{\Upsilon}}|)\leq\epsilon_{G}n_{G}n_{\psi_{0}}.$
(57)

If $\epsilon_{\text{CL}}\sim\epsilon_{Q}$, it’s sufficient to choose
$\epsilon_{G}n_{G}n_{\psi_{0}}\sim\epsilon$ and the results follow. ∎

Before stating our main theorems, we prove the following lemma.

######  Lemma 16.

A quantum algorithm can be constructed that takes the following inputs: (i)
sparse access
$(s,\|\mathcal{M}\|_{max},O_{M},O_{F})$ to a $2^{m}\times 2^{m}$ invertible
Hermitian matrix $\mathcal{M}$ such that $\|\mathcal{M}\|\leq
s\|\mathcal{M}_{max}\|$ (ii) $m$-qubit unitary $L(t_{n},\boldsymbol{j}/n)$
where $L(t_{n},\boldsymbol{j}/N)|0\rangle=|G_{n,\boldsymbol{j}}\rangle$ (iii)
an accuracy $\epsilon^{\prime}\in[1/2^{m},1]$ and (iv) $m$-qubit unitary black
box $U_{initial}$ where $U_{initial}|0\rangle=|\psi_{0}\rangle$. The algorithm
then returns with probability at least $2/3$ an $\epsilon^{\prime}$-additive
approximation to
$\Upsilon\equiv\langle\psi_{0}|(\mathcal{M}^{-1})^{\dagger}\mathcal{G}\mathcal{M}^{-1}|\psi_{0}\rangle$
where $\mathcal{G}\equiv|G_{n,\boldsymbol{j}}\rangle\langle
G_{n,\boldsymbol{j}}|$. This algorithm makes
$\mathcal{O}(\kappa^{2}/(\|\mathcal{M}\|\epsilon^{\prime}))$ queries to
$U_{\mathcal{G}}$ and $U_{initial}$,
$\mathcal{O}(s\|\mathcal{M}\|_{max}\kappa^{3}\log(\kappa^{2}/(\|\mathcal{M}\|\epsilon^{\prime}))/(\|\mathcal{M}\|\epsilon^{\prime}))$
queries to sparse oracles for $\mathcal{M}$ and
$\mathcal{O}(\kappa^{2}(m+s\|\mathcal{M}\|_{max}\kappa\log(\kappa^{2}/(\|\mathcal{M}\|\epsilon^{\prime}))))/(\|\mathcal{M}\|\epsilon^{\prime}))$
additional $2$-qubit gates.

###### Proof.

See Appendix F for details of the proof. ∎

To estimate the observable in Definition 11, a quantum algorithm is employed
to estimate $\Upsilon$ (see Appendix H for a summary of the basic steps) with
the following quantum query and gate complexities.

######  Theorem 17.

A quantum algorithm that takes sparse access
$(s=O(d),\|\mathcal{M}\|_{max}=O(1),O_{M},O_{F})$ to $\mathcal{M}$, where
$\|\mathcal{M}\|=O(1)$, and access to the unitaries
$L(t_{n},\boldsymbol{j}/N)$, where
$L(t_{n},\boldsymbol{j}/N)|0\rangle=|G_{n,\boldsymbol{j}}\rangle$ and
$U_{initial}$, where $U_{initial}|0\rangle=|\psi_{0}\rangle$, is able to
estimate the ensemble average $\langle G(T,x)\rangle$ at time $T=t_{n}$, with
$M$ initial data, to precision $\epsilon$ with an upper bound on the query
complexity $\mathcal{Q}$

$\displaystyle\mathcal{Q}=\mathcal{O}\left(\frac{n_{\psi_{0}}^{2}d^{7}T^{3}}{\epsilon^{10}}\log\left(\frac{n_{\psi_{0}}^{2}d^{4}T^{2}}{\epsilon^{7}}\right)\right)$
(58)

and the same order of additional $2$-qubit gates, where we suppress all $O(1)$
terms except $T=O(1)$ and $n_{\psi_{0}}\geq O(1)$.

###### Proof.

See Appendix G for details of the proof. ∎

The complexity in Eq. (58) has several notable features. Firstly, we see that
$\mathcal{Q}$ is independent of $M$, whereas the classical cost
$\mathcal{C}=O(MTd^{d+4}(1/\epsilon)^{d+1})$ to solve the Hamilton-Jacobi
equation from Lemma 1 is linear in $M$. Secondly, the classical cost
$\mathcal{C}$ contains exponential terms in dimension like $d^{d}$ and
$(1/\epsilon)^{d}$, which is absent in applications where $n_{\psi_{0}}$ does
not grow so quickly. Thirdly, this quantum algorithm is valid for all $T$,
whereas some previous quantum algorithms for nonlinear PDEs the linear
approximations are no longer valid for larger $T$, for any nonlinearity (e.g
[8]). This is particularly important in the case of strong nonlinearity where
it maybe interesting to study the large $T$ behavior where the nonlinearity
cannot be well-approximated by low-order polynomials. However, there is no
quantum advantage in $T$ for $T>1$ since $\mathcal{C}$ is linear in $T$ while
the quantum algorithm depends on $T^{3}$. We note that this absence of quantum
advantage in $T$ also holds for existing quantum algorithms for both linear
and nonlinear ODEs and PDEs. Thus, there is potential quantum advantage on the
three parameters $M,d,\epsilon$.

######  Corollary 18.

Let $\mathcal{C}$ be the cost to compute an observable from a
$(d+1)$-dimensional Hamilton-Jacobi equation with purely classical methods and
$\mathcal{Q}$ be the cost of our quantum algorithm. We say there is a quantum
advantage in estimating the observables when
$\mathcal{Q}={\scriptstyle\mathcal{O}}(\mathcal{C})$. To attain a quantum
advantage it is then sufficient for the following condition to hold

$\displaystyle{\scriptstyle\mathcal{O}}\left(\frac{Md^{d-4}}{n_{\psi_{0}}^{2}T^{2}}\left(\frac{1}{\epsilon}\right)^{d-9}\right)=\tilde{O}(1)$
(59)

where $\tilde{O}$ suppresses all logarithmic terms in $d,1/\epsilon,T$.

###### Proof.

The requirement $\mathcal{Q}={\scriptstyle\mathcal{O}}(\mathcal{C})$ implies
that $\mathcal{C}$ increases much more quickly than $\mathcal{Q}$, i.e., in
the asymptotic limit, the ratio $\mathcal{C}/\mathcal{Q}$ goes to infinity.
This can be equivalently expressed as
${\scriptstyle\mathcal{O}}(\mathcal{C}/\mathcal{Q})=1$. For instance, this can
include both polynomial and exponential advantages for the quantum algorithm.
From Lemma 14, the constant $n_{\psi_{0}}$ cannot be greater than $N^{d/2}$.
This means that the extra logarithmic terms in the upper bound of
$\mathcal{Q}$ is at most linear in $d$ multiplied by logarithmic factors in
$d,1/\epsilon,T$. Using the tight upper bound for $\mathcal{Q}$ from Theorem
17 and $\mathcal{C}=O(MTd^{d+4}(1/\epsilon)^{d+1})$ from Lemma 1, the result
follows. ∎

In this case, it is clear that even when $M=1$, there is potential quantum
advantage in parameters $d$ and $\epsilon$. Depending on the initial data, we
have $n_{\psi_{0}}^{2}=O(N^{b})$, where $b$ ranges from $0$ to $d$ from Lemma
14. Ignoring $T=O(1)$ factors, using Lemma 13, we see that Corollary 18
requires
$\tilde{O}(d^{d-4-b}/\epsilon^{d-9-3b})={\scriptstyle\mathcal{O}}(1)$. Thus if
$0\leq b<d/3-3$ and $d$ is large, one can obtain quantum advantage in both $d$
and $\epsilon$. For instance, in the best case scenario where $b=O(1)$, then
quantum advantage for $d$ and $\epsilon$ is always possible for high enough
$d$. However, in the worst-case scenario where $b=d$, we see that no advantage
is possible and the quantum algorithm actually performs worse than the
classical algorithm when $M=1$.

Different $b$ corresponds to different physical scenarios. From Lemma 14 we
see that if the initial data has support in a box of size $\beta$, then
$n_{\psi_{0}}=O((\beta\sqrt{N})^{d})=O(N^{b/2})$. For a Gaussian source of
width $\beta\sim 1/\sqrt{N}$ for instance, one has $n_{\psi_{0}}=O(1)$ so
$b=0$. For $b=k$ a constant, one needs $\beta=O(1/N^{1/2-k/(2d)})$. The
requirement $b<d/3-3$ for quantum advantage implies a support of the initial
data with the upper bound $\beta<O(1/N^{1/3+3/(2d)})$.

Often one needs to solve a PDE with many different initial data, such as those
in Monte-Carlo simulation, uncertainty quantification or Bayesian inference-
based inverse problems in uncertainty quantification, which demands large $M$.
Since $\mathcal{Q}$ is independent of $M$ and $\mathcal{C}$ is linear in $M$,
it is always (theoretically) possible to have large enough $M$ for a quantum
advantage. The reason for this $M$-independence lies in the fact that, in our
formulation, the initial data for $\psi$ is written as a single sum over $M$
initial data. This means that even if we consider state preparation costs for
$|\psi_{0}\rangle$, it would be still independent of $M$ since the sum over
$M$ terms lies within a single amplitude. In addition, the $M$ different
initial data can be arbitrarily distributed and there are no constraints on
neighbouring initial points to be close together to form a good mesh for high
accuracy in solving the PDE.

It is important to point out that, since the Hamilton-Jacobi equation is no
longer valid beyond the formation of caustics, the comparison in Corollary 18
only makes sense when the solution to the Hamilton-Jacobi equation is smooth.
Beyond the time of caustics, one should compare the quantum cost with the
classical cost of solving the Liouville equation (or other more sophisticated
classical algorithms laid out in at the beginning of section II.1.1). As shown
in Corollary 19, our quantum advantage is much bigger when compared with the
classical solver for the Liouville equation, except now there is no advantage
in $M$.

######  Corollary 19.

If $\mathcal{C}_{L}$ is the cost to compute an observable of a
$(2d+1)$-dimensional Liouville equation using purely classical methods and
$\mathcal{Q}$ is the cost with a quantum algorithm, then
$\mathcal{Q}={\scriptstyle\mathcal{O}}(\mathcal{C}_{L})$ when

$\displaystyle{\scriptstyle\mathcal{O}}\left(\frac{d^{2d-5}}{n^{2}_{\psi_{0}}T^{2}}\left(\frac{1}{\epsilon}\right)^{6d-7}\right)>\tilde{O}(1).$
(60)

###### Proof.

From Lemma 13, $\mathcal{C}_{L}=O(Td^{2d+3}(1/\epsilon)^{6d+3})$, which we
note is also independent of $M$.
Since $\mathcal{Q}=\tilde{O}(n^{2}_{\psi_{0}}d^{8}T^{3}(1/\epsilon)^{10})$
where $\tilde{O}$ suppresses all logarithmic terms in $d,T,1/\epsilon$, we
have our result. ∎

In this case, the classical cost is also independent of $M$, but has a large
overhead in $d$ and $1/\epsilon$. Since $b$ ranges from $0$ to $d$ when
$n^{2}_{\psi_{0}}=O(N^{b})$, we see that $\mathcal{C}_{L}>\mathcal{Q}$ if
$\tilde{O}(d^{2d-5-b}/\epsilon^{6d-7-3b})>O(1)$, ignoring $T\sim O(1)$ terms.
This means that even in the worst-case scenario where $s=b$ (when the initial
condition has the largest support $\beta=O(1)$), there is quantum advantage in
both $d$ and $\epsilon$ for large $d$.

We remark that it is also always possible to have a quantum subroutine for the
PDE problem (QLSP), where we use the HHL or CKS algorithms to prepare the
quantum state

$\displaystyle|\psi\rangle\equiv\frac{1}{N_{\psi}}\sum_{n=0}^{N_{t}-1}\sum_{\boldsymbol{j}}^{N}\sum_{\boldsymbol{l}}^{N}\psi_{n,\boldsymbol{j},\boldsymbol{l}}|\boldsymbol{j}\rangle|\boldsymbol{l}\rangle|n\rangle$
(61)

with $\psi_{n,\boldsymbol{j},\boldsymbol{l}}$ being the solutions to the
discretised version of the Liouville equation in Eq. (22) and $N_{\psi}$ is
the normalisation constant. This encoding of the solutions of the PDE, which
we can call the level-set encoding, differs from the amplitude encoded states
usually considered in the literature, where the amplitudes of the quantum
state are directly proportional to the solutions of the original PDE one wants
to solve. Although the amplitudes of this state
$\psi_{n,\boldsymbol{j},\boldsymbol{l}}$ likewise contains all the information
required for the solution, this encoding is in fact much more appropriate in
this instance to allow observables to be extracted. Given this quantum
subroutine, one can for instance apply the quantum swap-test between states
$|\psi\rangle$ and $|G\rangle$ to recover $\Upsilon$ (although with worse
query complexity compared to the current algorithm in Theorem 17). However,
any similar approach to recover physical observables from quantum states whose
amplitudes are directly proportional to the original solutions $u$ themselves,
requires one to inject the Jacobian factors in Definition 11 explicitly.
However, these factors are not generally a priori known since they depend on
the solutions of the PDE and also cannot be ignored, since they can have real
physical consequences as discussed at the end of Section III.2.1. The benefit
of a level-set encoded state is that these Jacobian factors are automatically
taken care of. This highlights the importance of choosing the right encoding
of the PDE solutions, even on the level of a quantum subroutine.

##  IV Solving scalar hyperbolic PDEs

We begin with a $(d+1)$-dimensional scalar nonlinear hyperbolic PDE

$\displaystyle\partial_{t}u^{[k]}+F(u^{[k]})\cdot\nabla_{x}u^{[k]}+Q(x,u^{[k]})=0,\quad
u\in\mathbb{R},$ (62) $\displaystyle u^{[k]}(0,x)=u^{[k]}_{0}(x),\quad
k=1,...,M$ (63)

subject to $M$ different initial data, where $x\in\mathbb{R}^{d}$. We can also
employ the level set formalism for this PDE and introduce a level set function
$\phi(t,x,p)$ in $(d+1)+1=d+2$ dimensions, where $p\in\mathbb{R}^{1}$. Its
zero level set is the solution $u$:

$\phi^{[k]}(t,x,p)=0\quad{\text{at}}\quad p=u^{[k]}(t,x).$ (64)

Then $\phi$ satisfies [25]

$\displaystyle\partial_{t}\phi^{[k]}+F(p)\cdot\nabla_{x}\phi^{[k]}-Q(x,p)\partial_{p}\phi^{[k]}=0,$
(65) $\displaystyle\phi^{[k]}(0,x,p)=p-u^{[k]}_{0}(x).$ (66)

Like for the Hamilton-Jacobi PDEs, we can similarly define a function $\psi$
by the following problem

$\displaystyle\partial_{t}\psi+F(p)\cdot\nabla_{x}\psi-Q(x,p)\partial_{p}\psi=0,$
(67)
$\displaystyle\psi(0,x,p)=\frac{1}{M}\sum_{k=1}^{M}\delta(p-u^{[k]}_{0}(x)),$
(68)

and one has

$\displaystyle\psi(t,x,p)=\frac{1}{M}\sum_{k=1}^{M}\delta(\phi^{[k]}(t,x,p)).$
(69)

The only difference compared to the Hamilton-Jacobi equation is that now
$p\in\mathbb{R}^{1}$ instead of being a $d$-dimensional vector. This means the
observable is now an integral over $\mathbb{R}$:

$\displaystyle\langle
g(t,x)\rangle=\int_{\mathbb{R}}g(p)\psi(t,x,p)dp\approx\frac{1}{N}\sum_{l=1}^{N}g_{l}\psi_{n,\boldsymbol{j},l}\equiv\langle
g^{\omega}_{n,\boldsymbol{j}}\rangle$ (70)

where after discretisation of the level set PDE
$\psi(t,x,p)\rightarrow\psi_{n,\boldsymbol{j},l}$, where $l$ is a scalar,
while $\boldsymbol{j}$ remains a vector. The initial condition can be encoded
in the quantum state

$\displaystyle|\psi^{\prime}_{0}\rangle=\frac{1}{N^{\prime}_{\psi_{0}}}\sum_{\boldsymbol{j}}^{N}\sum_{l}^{N}\psi_{0,\boldsymbol{j},l}|\boldsymbol{j}\rangle|l\rangle|n=0\rangle$
(71)

with the normalisation
$N^{\prime}_{\psi_{0}}=\sqrt{\sum_{\boldsymbol{|j|=1}}^{N}\sum_{l=1}^{N}|\psi_{0,\boldsymbol{j},l}|^{2}}$
and we also define the state

$\displaystyle|g_{n,\boldsymbol{j}}\rangle\equiv\frac{1}{N_{g}}\sum_{l=1}^{N}g^{*}_{l}|l\rangle|\boldsymbol{j}\rangle|n\rangle$
(72)

with the normalisation $N_{g}=\sqrt{\sum_{l}^{N}|g^{*}_{l}|^{2}}$. Defining
$\Upsilon$ in the same way as Eq. (III.3), with
$\mathcal{G}\equiv|g_{n,\boldsymbol{j}}\rangle\langle g_{n,\boldsymbol{j}}|$,
we can write

$\displaystyle\langle
g(t,x)\rangle\approx\frac{1}{N}N^{\prime}_{\psi_{0}}N_{g}|\sqrt{\Upsilon}|=n^{\prime}_{\psi_{0}}n_{g}|\sqrt{\Upsilon}|$
(73)

where $n^{\prime}_{\psi_{0}}\equiv N^{\prime}_{\psi_{0}}/\sqrt{N}$,
$n_{g}\equiv N_{G}/\sqrt{N}$. The normalisations are now the following.

######  Lemma 20.

The constant $n_{g}\equiv N_{G}/\sqrt{N}=O(1)$ and the range of
$n^{\prime}_{\psi_{0}}$ lies in $O(1)\leq n^{\prime}_{\psi_{0}}\equiv
N^{\prime}_{\psi_{0}}/\sqrt{N}\leq O(N^{d/2})$. If we assume the initial data
has support in a box of size $\beta$, then
$n_{\psi_{0}}=O((\beta\sqrt{N})^{d/2})$.

###### Proof.

The proof proceeds in a similar way to Lemma 14. See Appendix E. ∎

The quantum query and gate complexities are equivalent to the Hamilton-Jacobi
case up to constants and the proof is identical.

######  Theorem 21.

A quantum algorithm that takes sparse access
$(s=O(d),\|\mathcal{M}\|_{max}=O(1),O_{\mathcal{M}},O_{F})$ to $\mathcal{M}$
with $\|\mathcal{M}\|=O(1)$ and access to the unitaries
$l(t_{n},\boldsymbol{j}/N)$ where
$l(t_{n},\boldsymbol{j}/N)|0\rangle=|g_{n,\boldsymbol{j}}\rangle$, and
$U_{initial}$ where $U_{initial}|0\rangle=|\psi^{\prime}_{0}\rangle$, is able
to estimate the ensemble average $\langle g(T,\boldsymbol{j}/N)\rangle$ at
time $T$, with $M$ initial data, to precision $\epsilon$ with an upper bound
on the query complexity $\mathcal{Q}$

$\displaystyle\mathcal{Q}=\mathcal{O}\left(\frac{(n^{\prime}_{\psi_{0}})^{2}d^{7}T^{3}}{\epsilon^{10}}\log\left(\frac{(n^{\prime}_{\psi_{0}})^{2}d^{4}T^{2}}{\epsilon^{7}}\right)\right)$
(74)

and the same order of additional $2$-qubit gates, where we suppress all $O(1)$
terms except $T=O(1)$ and $n^{\prime}_{\psi_{0}}\geq O(1)$.

The classical cost for the scalar hyperbolic PDE is just one factor of $d$
smaller than for the Hamilton-Jacobi PDE we considered.

######  Corollary 22.

Let $\mathcal{C}$ be the cost to compute an observable from a
$(d+1)$-dimensional scalar hyperbolic equation with purely classical methods,
and $\mathcal{Q}$ be the cost to solve the identical problem on a quantum
algorithm. We say there is a quantum advantage in estimating physical
observables when $\mathcal{Q}={\scriptstyle\mathcal{O}}({\mathcal{C}})$. This
requires

$\displaystyle{\scriptstyle\mathcal{O}}\left(\frac{Md^{d-5}}{(n^{\prime}_{\psi_{0}})^{2}T^{2}}\left(\frac{1}{\epsilon}\right)^{d-9}\right)>\tilde{O}(1)$
(75)

where the $\tilde{O}$ notation suppresses all factors logarithmic in
$1/\epsilon,d,T$.

###### Proof.

Identical to Corollary 18 except using
$\mathcal{C}=O(MTd^{d+3}/\epsilon^{d+1})$ from Lemma 2. ∎

The consequences are very similar to that of Hamilton-Jacobi equations. When
$M=1$ and $(n^{\prime}_{\psi_{0}})^{2}=O(N^{b})$, for small enough $b$ where
$0\leq b<d/3-3$ and large $d$ one can expect quantum advantage in $d$ and
$\epsilon$. From Lemma 20,
$n^{\prime}_{\psi_{0}}=O((\beta\sqrt{N})^{d/2})=O(N^{b/2})$. A quantum
advantage thus requires $\beta<O(1/N^{1/3+3/(2d)})$. Since $\mathcal{Q}$ is
independent of $M$, the quantum algorithm will always be advantageous with
respect to $M$ in the large $M$ limit.

##  V Solving system of nonlinear ODEs

We aim to compute observables from the following system of $D$ nonlinear ODEs

$\displaystyle\frac{dX^{[k]}(t)}{dt}=F(X^{[k]}(t)),\qquad
X^{[k]}\in\mathbb{R}^{D}\,,$ (76) $\displaystyle X^{[k]}(0)=X_{0}^{[k]},\qquad
k=1,\cdots,M$ (77)

subject to $M$ different initial data, where the same ODE (i.e., same function
$F(X)$ which can be nonlinear) is satisfied for each $k$. This can be viewed
as a system of $M$ non-interacting particles in $D$ dimensions, each with a
trajectory described by $X^{[k]}(t)$. A traditional way to compute ensemble
averages

$A(t)=\frac{1}{M}\sum_{k=1}^{M}A(X^{[k]}(t))$ (78)

is to directly solve for the system of ODEs starting from $M$ different
initial data $X_{0}^{[k]}$, and then carry out the summation in Eq. (78).

We note that the proposal below has appeared in [11], but no explicit details
on the quantum costs were given, which we provide below. This is also very
similar to the Koopman-von Neumann approach in [10], except here we consider
modelling the $M$ initial conditions as trajectories of $M$ non-interacting
particles.

###  V.1 Mapping nonlinear ODEs to a linear PDE

Assume $X^{[k]}(t)$ is the solution to Eq. (76). Then one can define a
function $\Phi(t,q):\mathbb{R}^{+}\times\mathbb{R}^{D}\rightarrow\mathbb{R}$

$\displaystyle\Phi(t,q)=\frac{1}{M}\sum_{k=1}^{M}\delta(q-X^{[k]}(t))$ (79)

One can easily check (Lemma 23) that $\Phi$ satisfies, in the weak sense, the
linear $D+1$-dimensional PDE

$\displaystyle\frac{\partial\Phi(t,q)}{\partial
t}+\nabla_{q}\cdot[F(q)\Phi(t,q)]=0$
$\displaystyle\Phi(0,q)=\frac{1}{M}\sum_{k=1}^{M}\delta(q-X_{0}^{[k]}).$ (80)

Remark: If $\nabla\cdot F=0$, which is the case for the Hamiltonian system,
then $J=1$, the dynamics in Eq. (76) is volume preserving, and equation (V.1)
becomes

$\frac{\partial\Phi(t,q)}{\partial t}+F(q)\cdot\nabla_{q}\Phi(t,q)=0$

To justify the correctness of the method, we have the following lemma.

######  Lemma 23.

If $X^{[k]}(t)$ is the solution to Eq. (76), then Eq. (79) is the weak
solution, in the sense of distribution, to the problem in Eq. (V.1).

###### Proof.

See Appendix I. ∎

Suppose we are interested in solving (76) with $M=N^{D}$ initial data, equally
spaced on a $D$-dimensional uniform mesh on $[0,1]^{D}$, with a uniform mesh
size $h=1/N$. This will be the mesh to solve the PDE problem in Eq.(V.1). The
discretisation process proceeds in the same way as Section III and we can
similarly convert the linear PDE in Eq. (V.1) into a linear algebra problem.
Let the discretisation of $\Phi(t,q)$ be

$\displaystyle\Phi(t,q)\rightarrow\Phi^{\omega}_{n,\boldsymbol{j}}$ (81)

where $n$ denotes the time step from $t\rightarrow n\Delta t$ and
$n=0,...,N_{t}$. The vector $\boldsymbol{j}=(j_{1},...,j_{D})$, denotes
spatial grid positions with grid size $h=1/N$ where $q\rightarrow
h\boldsymbol{j}$. The initial state at $n=0$ can be approximated by

$\displaystyle\Phi^{\omega}_{0,\boldsymbol{j}}=\frac{1}{M}\sum_{k=1}^{M}\delta_{\omega}(\boldsymbol{j}h-X^{[k]}_{0})$
(82)

where $\delta_{\omega}$ is the smoothed delta function, defined in Section
III.1, and $\omega$ is a smoothing parameter. The linear algebra problem to
solve for $\Phi^{\omega}_{n,\boldsymbol{j}}$ is then identical to Eq. (33),
with the replacement
$\psi_{n,\boldsymbol{j},\boldsymbol{l}}\rightarrow\Phi_{n,\boldsymbol{j}}$ and
$\mathcal{M}\rightarrow\mathcal{M}_{\text{ODE}}$. Here
$\mathcal{M}_{\text{ODE}}$ is an $N_{t}N^{D}\times N_{t}N^{D}$ Hermitian
matrix. See Appendix K for the precise form of $\mathcal{M}_{\text{ODE}}$ and
details of the discretisation procedure.

######  Lemma 24.

The condition number of $\mathcal{M}_{\text{ODE}}$ is $\kappa\leq O(DNT)$ and
sparsity is $s=O(D)$.

###### Proof.

The proof is identical to the proof in Appendix C with the replacement
$d\rightarrow D/2$. ∎

###  V.2 Observables of the ODE

Any ensemble-averaged observable obtained from the solutions $X^{[k]}$ of the
system of nonlinear ODEs can be written in the following form.

######  Definition 25.

For any function $A(X)$, one can define the ensemble average via

$\displaystyle\langle A(t)\rangle\equiv\int
A(q)\Phi(t,q)\,dq=\frac{1}{M}\sum_{k=1}^{M}A(X^{[k]}(t))$ (83)

which we identify to be the observable associated with the ODE whose solutions
are $\\{X^{[k]}(t)\\}$.

The second equality in Definition 25 is justified in Lemma 26, where we prove
that, if $\Phi(x,q)$ is the solution to (V.1), and $\langle A(t)\rangle$ is
defined as in Eq. (83), then

$\langle A(t)\rangle=\frac{1}{M}\sum_{k=1}^{M}A(X^{[k]}(t)),$ (84)

where $X^{[k]}$ is the solution to Eq. (76). Since $\langle A(t)\rangle$ is an
integral of the product of $A(q)$ and $\Phi(t,q)$, we can approximate this
quantity with a very similar quantum algorithm that we used in Sections III
and IV.

######  Lemma 26.

Assume $\Phi$ is the solution to Eq. (V.1), and $X^{[k]}$ is the solution to
Eq. (12), then Eq. (84) holds.

###### Proof.

See Appendix J. ∎

We remark that here the observable $\langle A(t)\rangle$ is not a sum of terms
$A(X^{[k]}(t))$ weighted by Jacobians, unlike in the case of nonlinear PDEs.
This key difference arises from the fact that during the evolution of the
function $\Phi(t,q)$, the form of the delta function is preserved, i.e., from
$t=0$ to $T$, the delta function terms in $\Phi(t,q)$ remains the form of
$\delta(q-X^{[k]}(t))$, so there is no Jacobian term (or, in fact, the
determinant of the Jacobian is always $1$ here). This is unlike the case for
Hamilton-Jacobi equation and hyperbolic PDEs, where for general $t$ the form
of the delta function terms at $t=0$ are not preserved (from Eq. (23) to Eq.
(24)).

Similarly to Section III.2 we can also define a new observable $A_{O}(t)$ that
is the ratio of two observables where the second observable is the identity
function. This can be considered as a normalised observable

$\displaystyle A_{O}(t)\equiv\frac{\langle
A(t)\rangle}{\mathbf{1}(t)}=\frac{\int A(q)\Phi(t,q)\,dq}{\int\Phi(t,q)\,dq}.$
(85)

Remark: Our algorithm does not recover individual $X^{[k]}(t)$, only the
ensemble average. Moreover, we can recover the ensemble average in any part of
the domain, say $\mathcal{D}$, by

$A_{P}^{\mathcal{D}}(t)=\frac{1}{\text{vol}(\mathcal{D})}\int
A(q)\chi_{\mathcal{D}}(q)\Phi(t,q)\,dq$

where $\chi(\mathcal{D})$ is the characteristic function on $\mathcal{D}$ and
$\text{vol}(\mathcal{D})$ is the volume of $\mathcal{D}$. We can do the
simulation in some bounded domain in $q$, so the support of $\Phi$ remains
inside the domain throughout the computational time. This is possible since
the transport equation (V.1) has a finite propagation speed, if the support of
$\Phi(0,q)$ is bounded.

The observable defined by the integral in Eq. (83) can then be approximated by
a quadrature method:

$\langle A(t)\rangle\approx\langle
A^{\omega}_{n}\rangle=\frac{1}{N^{D}}\sum_{\boldsymbol{j}}^{N}A_{\boldsymbol{j}}\Phi^{\omega}_{n,\boldsymbol{j}},$
(86)

where one uses the discretisation $A(q)\rightarrow A_{\boldsymbol{j}}$. For
details of the discretisation see Appendix K.

The total computational cost in computing observables from Eq. (V.1) is
similar that in Lemma 13, where Eq. (V.1) has dimension $D+1$.

######  Lemma 27.

The overall computational cost will be of
$O(DN_{t,\text{CL}}N_{\text{CL}}^{D})=O(D^{2}TN_{\text{CL}}^{D+1})=O(D^{D+3}T(1/\epsilon_{\text{CL}})^{3D+3})$.

Remark: This method is related to the particle methods [48] which is a popular
classical method used to solve linear transport equations in high dimensions.
There one uses particles, which satisfy the ODEs (76), to solve the linear
transport (Liouville (V.1) here) equation. Our method here is, in a sense,
like an inverted particle method–here we use the Liouville equation (V.1) to
approximate the particle (ODE) system (76), instead of the other way around.

###  V.3 The quantum algorithm to estimate the ensemble average

We proceed in a similar way to Section III.2 and define the initial quantum
state that encodes the initial condition

$\displaystyle|\Phi_{0}\rangle=\frac{1}{N_{\Phi_{0}}}\sum_{\boldsymbol{j}}^{N}\Phi_{0,\boldsymbol{j}}|\boldsymbol{j}\rangle|n=0\rangle$
(87)

where $\Phi_{0,\boldsymbol{j}}$ is defined in Eq. (159). The normalisation is
given by
$N_{\Phi_{0}}=\sqrt{\sum_{\boldsymbol{j}}^{N}|\Phi_{0,\boldsymbol{j}}|^{2}}$.
We can also define the state

$\displaystyle|A^{n}\rangle\equiv\frac{1}{N_{A}}\sum_{\boldsymbol{j}}^{N}A^{*}_{\boldsymbol{j}}|\boldsymbol{j}\rangle|n\rangle$
(88)

where we note that here we do not sum over the time step index $n$, where
$n\Delta t$ is the time-step that we later want to compute $\langle
A^{\omega}_{n}\rangle$ in Eq. (89). The normalisation is
$N_{A}=\sqrt{\sum_{\boldsymbol{j}}^{N}|A^{*}_{\boldsymbol{j}}|^{2}}$.

Using Eq. (86) and given the states $|\Phi\rangle$ and $|A^{n}\rangle$, one
can recover the observable from

$\displaystyle\langle A(t_{n})\rangle\approx\langle
A^{\omega}_{n}\rangle\equiv\frac{1}{N^{D}}\sum_{\boldsymbol{j}=1}^{N}A_{\boldsymbol{j}}\Phi_{\omega}(t_{n},\boldsymbol{j})=n_{\Phi_{0}}n_{A}|\sqrt{\Upsilon_{\text{ODE}}}|$
(89)

where $n_{\Phi_{0}}\equiv N_{\Phi_{0}}/N^{D/2}$, $n_{A}\equiv N_{A}/N^{D/2}$
and a quantum algorithm allows one to compute

$\displaystyle\Upsilon_{\text{ODE}}\equiv\langle\Phi_{0}|(\mathcal{M}_{\text{ODE}}^{-1})^{\dagger}\mathcal{A}\mathcal{M}_{\text{ODE}}^{-1}|\Phi_{0}\rangle.$
(90)

Now the normalisations can similarly be estimated like in Section III.2

######  Lemma 28.

The constants $n_{A}\equiv N_{A}/N^{d/2}=O(1)$ and $n_{\Phi_{0}}\equiv
N_{\Phi_{0}}/N^{D/2}=O(1)$.

###### Proof.

See Appendix E. ∎

The quantum algorithm only allows one to estimate $\Upsilon_{\text{ODE}}$ to
finite precision, so we denote the quantum estimate as

$\displaystyle\langle\tilde{A}^{\omega}_{n}\rangle\equiv
n_{\Phi_{0}}n_{A}|\sqrt{\tilde{\Upsilon}_{\text{ODE}}}|$ (91)

where
$|\sqrt{\Upsilon_{\text{ODE}}}-\sqrt{\tilde{\Upsilon}_{\text{ODE}}}|\leq\epsilon_{A}$.
The error analysis then proceeds in exactly the same way as Section III.2.
Then we have the following contributions to the total error in estimating
$A_{O}(t_{n})$.

######  Lemma 29.

The error $|\langle A(t_{n})\rangle-\langle\tilde{A}^{\omega}_{n}\rangle|$ can
be broken into two independent sources of error

$\displaystyle|\langle
A(t_{n})\rangle-\langle\tilde{A}^{\omega}_{n}\rangle|\leq|\langle
A(t_{n})\rangle-\langle A^{\omega}_{n}\rangle|+|\langle
A^{\omega}_{n}\rangle-\langle\tilde{A}^{\omega}_{n}\rangle|$
$\displaystyle=\epsilon_{C}+\epsilon_{Q}=\epsilon.$ (92)

Here $\epsilon_{C}$ comes from approximating the solutions of the discretised
linear PDE (V.1), as estimated in Eq. (165). The quantum error comes from the
estimate of $\Upsilon_{\text{ODE}}$, where
$|\sqrt{\Upsilon_{\text{ODE}}}-\sqrt{\tilde{\Upsilon}_{\text{ODE}}}|\leq\epsilon_{A}$.
Imposing $\epsilon_{C}\sim\epsilon_{Q}\sim\epsilon$, then we require
$\epsilon_{A}\sim\epsilon/(n_{\Phi_{0}}n_{A})\sim\epsilon_{C}/(n_{\Phi_{0}}n_{A})$.

The query and gate complexity for estimating the observable $\langle
A(t_{n})\rangle$ is then identical to Theorem 17 with the replacement
$d\rightarrow D/2$.

######  Theorem 30.

The worst-case total query complexity $\mathcal{Q}$ to estimate the observable
$\langle A(t_{n})\rangle$ to precision $\epsilon$ on a quantum algorithm that
takes sparse access $(s=O(D),\|\mathcal{M}\|_{max}=O(1),O_{M},O_{F})$ to
$\mathcal{M}_{\text{ODE}}$, access to the unitaries $J(n)$ where
$J(n)|0\rangle=|A^{n}\rangle$ and $U_{\Phi_{0}}$ where
$U_{\Phi_{0}}|0\rangle=|\Phi_{0}\rangle$, is

$\displaystyle\mathcal{Q}=\mathcal{O}\left(\frac{n_{\Phi_{0}}^{2}D^{7}T^{3}}{\epsilon^{10}}\log\left(\frac{n_{\Phi_{0}}^{2}D^{4}T^{2}}{\epsilon^{7}}\right)\right)$
(93)

where all constant terms $O(1)$ are suppressed except $T=O(1)$ and
$n_{\Phi_{0}}=O(1)$. This complexity is independent of $M$.

###### Proof.

The proof is the same as Theorem 17 except with the replacement $D=2d$. ∎

######  Corollary 31.

We say there is a quantum advantage in estimating the observables $\langle
A(T)\rangle$ to precision $\epsilon$ when
$\mathcal{Q}={\scriptstyle\mathcal{O}}(\mathcal{C})$, which requires

$\displaystyle{\scriptstyle\mathcal{O}}\left(\frac{M}{n^{2}_{\Phi_{0}}D^{4}T^{2}(1/\epsilon)^{9}}\right)=\tilde{O}(1).$
(94)

###### Proof.

From Lemma 3, the classical cost is $\mathcal{C}=O(MD^{3}T/\epsilon)$. Since
$n_{\Phi_{0}}=O(1)$, then from Theorem 30, we see
$\mathcal{Q}=\tilde{O}(n^{2}_{\Phi_{0}}D^{7}T^{3}/\epsilon^{10})$ where
$\tilde{O}$ suppresses all logarithmic terms in $D,\epsilon$. ∎

Here the only possible quantum advantage is with respect to $M$. Ignoring
$n^{2}_{\Phi_{0}}=O(1)$ and $T=O(1)$ constant factors, quantum advantage is
only possible for $M>\tilde{O}(D^{4}(1/\epsilon)^{9})$.

We note that here, unlike in the Hamilton-Jacobi and hyperbolic PDEs where the
level set formalism is used, the $M$ initial data in the ODE case needs to be
close to each other for neighboring initial points to form a mesh good enough
for the accuracy of solving the PDE in Eq. (V.1). There is no such constraint
for the nonlinear PDE problems in Section III.

In addition, the quantum algorithm for the nonlinear PDEs are much more
efficient than the quantum ODE solver, relatively speaking, when compared to
their respective classical counterparts. In the ODE case, $D$ ODEs need to be
solved by a $(D+1)$-dimensional linear PDE, while a $(d+1)$-dimensional
Hamilton-Jacobi PDE, describes an infinite dimensional dynamical system, can
be represented by only a $(2d+1)$-dimensional PDE.

##  VI Solving more general nonlinear PDEs

In the special cases of nonlinear Hamilton-Jacobi and hyperbolic partial
differential equations we saw how it is possible to reformulate the problem to
linear partial differential equations by just (at most) doubling the dimension
using the level set formalism. This same result, however, cannot be done
analytically for general nonlinear PDEs.

To devise quantum algorithms to solve more general $(d+1)$-dimensional
nonlinear PDEs, like the Euler and Navier-Stokes equations in fluid dynamics,
a naive way is to first numerically approximating the system so they become a
system of nonlinear ODEs, and then use the quantum algorithm in Section V.
Below we present two methods to achieve this: (i) the Lagrangian
discretisation method and (ii) the Eulerian discretisation method.

###  VI.1 The Lagrangian discretisations

One of the most important Lagrangian discretisation method in incompressible
flows is the vortex method [49]. For the incompressible Euler equations the
vortex method solves a Hamiltonian particle system like in Eq. (12). Although
the diffusion term in the Navier-Stokes equation is usually modelled by
stochastic ODEs (with Brownian motion for diffusion), there are also
deterministic versions. This method of solving a $(d+1)$-dimensional nonlinear
Euler or Navier-Stokes equation involves a description of $dN$ particles each
obeying an ODE, where $N=O(d/\epsilon)$ and $\epsilon$ is the error in
computing observables from the resulting ODE solutions. This means these PDEs
can be reduced to a system of $O(d^{2}/\epsilon)$ nonlinear ODEs.

One can employ the quantum algorithm from Section V with quantum query and
gate complexity cost $\mathcal{Q}=\tilde{\mathcal{O}}(d^{14}/\epsilon^{17})$,
where we used Theorem 30 with the replacement $D\rightarrow d^{2}/\epsilon$
and $\tilde{\mathcal{O}}$ suppressing logarithmic terms in $d$, $\epsilon$ and
all $O(1)$ terms. The cost of classically solving the $(d+1)$-dimensional
nonlinear PDE with a system of $dN$ nonlinear ODEs is
$\mathcal{C}=O(MD^{3}/\epsilon)=O(Md^{6}/\epsilon^{4})$ from Lemma 3. This
quantum cost is independent of $M$ whereas the classical algorithm would
depend linearly on $M$. Thus for a large enough
$M>\tilde{O}(d^{7}/\epsilon^{13})$, one can expect an advantage in employing
the quantum algorithm.

There are other particle or mesh-free methods used for some other nonlinear
PDEs like smoothed particle hydrodynamics (SPH) in solid mechanics and fluid
flows [50]. The Boltzmann equation for rarefied gas is often solved by the
Direct Simulation Monte-Carlo method, which is a stochastic particle method
[51], and particle-in-cell (PIC) methods are often used for the Vlasov-Poisson
or Vlasov-Maxwell systems [52] in plasma physics.

###  VI.2 The Eulerian discretisations

Grid-based Eulerian discretisations solve the PDEs on a fixed grid. They offer
higher order accuracies but suffer from the curse-of-dimensionality. When
given the most general nonlinear $(d+1)$-dimensional PDE, one can first
discretise the $d$ spatial dimensions with mesh size $1/N$. This then becomes
a system of $D=N^{d}$ nonlinear ODEs. Since $N=O(d/\epsilon)$, where
$O(\epsilon)$ is the error in computing observables from the ODE system as
well as the error in the original PDE, this gives $D=(d/\epsilon)^{d}$. Using
Lemma 3 and Theorem 30 with the replacement $D\rightarrow(d/\epsilon)^{d}$,
one obtains $\mathcal{Q}=\tilde{\mathcal{O}}(d^{7d}/\epsilon^{10+7d})$ and
$\mathcal{C}=O(Md^{3d}/\epsilon^{1+3d})$. Thus one requires
$M>\tilde{O}(d^{4d}/\epsilon^{9+4d})$ for a quantum advantage. Since the size
of $M$ scales exponentially with $d$, it does not fit known realistic
problems. Better quantum methods to tackle these more general PDEs are thus
still required.

##  VII Summary and discussion

We introduced quantum algorithms for computing physical observables of
nonlinear Hamilton-Jacobi and scalar hyperbolic PDEs with arbitrary
nonlinearities. There are potential quantum speedups with respect to dimension
$d$, precision $\epsilon$, the number of initial data $M$ and for arbitrary
nonlinearity, with no large overheads in any other parameters. This would be
of potential interest for the field of high-dimensional Hamilton-Jacobi
equations, which appear in many areas including optimal control, mean-field
games and machine learning. We found tight upper asymptotic query and gate
complexities of quantum algorithms that compute ensemble averages of nonlinear
PDEs and ODEs. This allows us to find a lower bound on the ratio between the
classical and quantum cost to obtain an identical output. This can be
equivalently expressed as

$\displaystyle\mathcal{O}\left(\frac{\mathcal{C}}{\mathcal{Q}}\right)=\tilde{O}\left(\frac{M}{T^{2}}d^{r_{1}}\left(\frac{1}{\epsilon}\right)^{r_{2}}\right)$
(95)

where $(d+1)$ denotes the dimension of the PDE and $\tilde{O}$ suppresses all
logarithmic factors in $d$, $T$, $1/\epsilon$. For a system of $D$ nonlinear
ODEs, we replace $d$ with $D$. The different values of the exponents $r_{1}$,
$r_{2}$ depends on the type of PDE and its initial conditions. Here one sees
that there is always a quantum advantage with respect to $M$ and no quantum
advantage for any $T>1$. For a quantum advantage in $d$ and $\epsilon$, it is
sufficient to have $r_{1}>0$ and $r_{2}>0$ respectively in the asymptotic
limit. These results are summarised in Table 1.

The query and gate complexities presented could be improved for instance by
using higher order methods or using preconditioners for solving the Liouville
equation. The quantum query and gate complexity with respect to the condition
number can also be further optimised. It is also worthwhile to consider these
quantum algorithms for real applications with large $M$ and to extend these
methods to stochastic PDEs.

We see that even though quantum algorithms that utilise the linear
representation for nonlinear PDEs (where the level set is one example) have
advantages for Hamilton-Jacobi and scalar hyperbolic equations, this method is
not necessarily generally applicable. For more general nonlinear PDEs, we see
that even though quantum algorithms with quantum advantage are theoretically
possible in regimes of large $M$, there are not necessarily generally useful
for real applications, so other methods are required. Other linear
representations for nonlinear PDEs, rather than linear approximations or
linear representations for ODEs, are desired if it is possible to find them.
It may be that some nonlinear PDEs are more amenable to a quantum treatment
than others and different classes of PDEs can be treated more efficiently with
different quantum methods, rather than trying to find a one-method-fits-all
strategy.

## Acknowledgements

N. Liu thanks Barry Sanders for very interesting and insightful discussions.
S. Jin was partially supported by the NSFC grant No. 12031013 and the Shanghai
Municipal Science and Technology Major Project (2021SHZDZX0102). N. Liu
acknowledges funding from the NSFC International Young Scientists Project (no.
12050410230), the Science and Technology Program of Shanghai (no. 21JC1402900)
and the Natural Science Foundation of Shanghai grant 21ZR1431000.

##  Appendix A Proof of Lemma 9

Denote the bicharacteristics of Eq. (22) by

$\frac{\partial x}{\partial t}=\nabla_{p}H,\quad\frac{\partial p}{\partial
t}=-\nabla_{x}H,\quad\qquad x(0)=x^{(0)},\quad p(0)=p^{(0)}.$ (96)

Assume sufficient smoothness of $H$, then the above ODE problem has a unique
solution, denoted by

$(x(t;x^{(0)},p^{(0)}),p(t;x^{(0)},p^{(0)})),$

which can be inverted to get the inverse functions
$(x^{(0)}(t;x,p),p^{(0)}(t;x,p))$. This is possible since the Jacobian matrix
of the map from $(x^{(0)},p^{(0)})$ to $(x(t),p(t))$ has determinant $1$. Then
by the method of characteristics,

$\phi_{i}^{[k]}(t,x,p)=\phi_{i}^{[k]}(0,x^{(0)}(t;x,p),p^{(0)}(t;x,p))=p_{i}^{(0)}(t;x,p)-u_{i}^{[k]}(0,x^{(0)}(t;x,p))$
(97)

and, also by the method of characteristics and using the above solution for
$\phi^{[k]}$,

$\psi(t,x,p)=\psi(0,x^{(0)}(t;x,p),p^{(0)}(t;x,p))=\frac{1}{M}\sum_{k=1}^{M}\prod_{i=1}^{d}\delta\left(p_{i}^{(0)}(t;x,p)-u_{i}^{[k]}(0,x^{(0)}(t;x,p))\right)=\frac{1}{M}\sum_{k=1}^{M}\delta(\phi^{[k]}(t,x,p))$
(98)

##  Appendix B Discretised PDE

We proceed by first discretising the PDE in Eq. (22) by, for example, finite
difference schemes. The state with amplitudes proportional to the solutions of
this discretised PDE can be created by using the quantum linear systems solver
subroutine. For instance, we use the forward Euler method in time and upwind
approximation in $x$ and $p$-derivatives,

$\displaystyle\frac{\partial\psi(t_{n},x,p)}{\partial
t}\rightarrow\frac{\psi_{n+1,\boldsymbol{j},\boldsymbol{l}}-\psi_{n,\boldsymbol{j},\boldsymbol{l}}}{\Delta
t}$ $\displaystyle\frac{\partial H}{\partial
p_{i}}\frac{\partial\psi(t_{n},x,p)}{\partial
x_{i}}\rightarrow\frac{1}{h}\left\\{\frac{\partial H}{\partial
p_{i}}\right\\}_{-}\left[T^{+}_{i}\psi_{n}-\psi_{n}\right]_{\boldsymbol{j},\boldsymbol{l}}+\frac{1}{h}\left\\{\frac{\partial
H}{\partial
p_{i}}\right\\}_{+}\left[\psi_{n}-T^{-}_{i}\psi_{n}\right]_{\boldsymbol{j},\boldsymbol{l}}$
$\displaystyle\frac{\partial H}{\partial
x_{i}}\frac{\partial\psi(t_{n},x,p)}{\partial
q_{i}}\rightarrow\frac{1}{h}\left\\{\frac{\partial H}{\partial
x_{i}}\right\\}_{+}\left[P^{+}_{i}\psi_{n}-\psi_{n}\right]_{\boldsymbol{j},\boldsymbol{l}}+\frac{1}{h}\left\\{\frac{\partial
H}{\partial
x_{i}}\right\\}_{-}\left[\psi_{n}-P^{-}_{i}\psi_{n}\right]_{\boldsymbol{j},\boldsymbol{l}}$

Here $\alpha_{+}=\max\\{\alpha,0\\},\alpha_{-}=\min\\{\alpha,0\\}$ for a
general quantity $\alpha$, $\boldsymbol{j}\equiv(j_{1},...,j_{d})$ and
$\boldsymbol{l}\equiv(l_{1},...,l_{d})$ denote vectors corresponding to the
discretised $x_{i}\rightarrow j_{i}h$ and $p_{i}\rightarrow l_{i}h$
respectively, where $j_{i},l_{i}=1,...,N$ and $i=1,...,d$, $n=0,1,...,N_{t}$.
The time step is denoted $t_{n}=n\Delta t$. In addition,
$T^{\pm}_{i}\psi_{n,\boldsymbol{j},\boldsymbol{l}}=\psi_{n,j_{1},...,j_{i\pm
1},...,j_{d},\boldsymbol{l}}$ are translation operators with respect to $x$
while
$P^{\pm}_{i}\psi_{n,\boldsymbol{j},\boldsymbol{l}}=\psi_{n,\boldsymbol{j},l_{1},...,l_{i\pm
1},...d}$ are translation operators with respect to $p$ while operating on
$\psi_{n,\boldsymbol{j},\boldsymbol{l}}$. Note here $\frac{\partial
H}{\partial p_{i}}$ and $\frac{\partial H}{\partial x_{i}}$ can be evaluated
analytically, for given $H(x,p)$, thus do not need to be approximated.

Define $\lambda=\Delta t/h$. We require the CFL condition

$d\lambda\,\max_{i}\sup_{x,p}\left\\{\left|\frac{\partial H}{\partial
x_{i}}\right|,\left|\frac{\partial H}{\partial p_{i}}\right|\right\\}\leq 1$
(99)

for numerical stability. Then the discretised version of Eq. (22) can be
rewritten as

$\displaystyle\psi_{n+1,\boldsymbol{j},\boldsymbol{l}}-\psi^{n}_{\boldsymbol{j},\boldsymbol{l}}+\lambda\sum_{i=1}^{d}\left[\left\\{\frac{\partial
H}{\partial p_{i}}\right\\}_{-}(T^{+}_{i}-I)+\left\\{\frac{\partial
H}{\partial p_{i}}\right\\}_{+}(I-T^{-}_{i})\right.$ (100)
$\displaystyle\left.-\left\\{\frac{\partial H}{\partial
x_{i}}\right\\}_{+}(P^{+}_{i}-I)-\left\\{\frac{\partial H}{\partial
x_{i}}\right\\}_{-}(I-P^{-}_{i})\right]\psi_{n,\boldsymbol{j},\boldsymbol{l}}=0,$
(101)

where $I$ is the identity operator, with the initial condition (for $n=0$)

$\displaystyle\psi_{0,\boldsymbol{j},\boldsymbol{l}}=\frac{1}{M}\sum_{k=1}^{M}\prod_{i=1}^{d}\delta_{\omega}(l_{i}h-u^{[k]}_{i}(n=0,\boldsymbol{j})),$
(102)

where $\delta_{\omega}$ is the smoothed delta function.

Now the discretised PDEs can be written as a matrix equation

$\displaystyle\mathcal{K}\begin{pmatrix}\psi_{1,\boldsymbol{j},\boldsymbol{l}}\\\
\psi_{2,\boldsymbol{j},\boldsymbol{l}}\\\ \vdots\\\
\psi_{N_{t}-1,\boldsymbol{j},\boldsymbol{l}}\\\
\psi_{N_{t},\boldsymbol{j},\boldsymbol{l}}\end{pmatrix}=\begin{pmatrix}\psi_{0,\boldsymbol{j},\boldsymbol{l}}\\\
0\\\ \vdots\\\ 0\\\ 0\end{pmatrix}$ (103)

and $\mathcal{K}$ is a $N_{t}N^{2d}\times N_{t}N^{2d}$ Toeplitz matrix of the
form

$\displaystyle\mathcal{K}=\begin{pmatrix}\mathbf{1}&0&0&\ldots&0&0&0\\\
K&\mathbf{1}&0&\ldots&0&0&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&\ldots&K&\mathbf{1}&0\\\
0&0&0&\ldots&0&-\mathbf{1}&\mathbf{1}\end{pmatrix}$ (104)

where each $\mathbf{1}$ above is the $N^{2d}\times N^{2d}$ identity matrix and
$K$ is the $N^{2d}\times N^{2d}$ matrix:

$\displaystyle K=-I+\lambda\sum_{i=1}^{d}\left[\left\\{\frac{\partial
H}{\partial p_{i}}\right\\}_{-}(T^{+}_{i}-I)+\left\\{\frac{\partial
H}{\partial p_{i}}\right\\}_{+}(I-T^{-}_{i})-\left\\{\frac{\partial
H}{\partial x_{i}}\right\\}_{+}(P^{+}_{i}-I)-\left\\{\frac{\partial
H}{\partial x_{i}}\right\\}_{-}(I-P^{-}_{i})\right].$ (105)

We can then solve for $\psi_{n,\boldsymbol{j},\boldsymbol{l}}$ by matrix
inversion

$\displaystyle\begin{pmatrix}\psi_{1,\boldsymbol{j},\boldsymbol{l}}\\\
\psi_{2,\boldsymbol{j},\boldsymbol{l}}\\\ \vdots\\\
\psi_{N_{t}-1,\boldsymbol{j},\boldsymbol{l}}\\\
\psi_{N_{t},\boldsymbol{j},\boldsymbol{l}}\end{pmatrix}=\mathcal{K}^{-1}\begin{pmatrix}\psi_{0,\boldsymbol{j},\boldsymbol{l}}\\\
0\\\ \vdots\\\ 0\\\ 0\end{pmatrix}.$ (106)

Both QLSP and SLEP involve the inversion of an Hermitian matrix $\mathcal{M}$.
Since $\mathcal{K}$ from Eq. (104) is not Hermitian, we can define a new
Hermitian matrix

$\displaystyle\mathcal{M}=\begin{pmatrix}0&\mathcal{K}\\\
\mathcal{K}^{\dagger}&0\end{pmatrix}$ (107)

which has the same sparsity and condition number as $\mathcal{K}$.

Using $\mathcal{M}$, the matrix inversion problem in order to solve
$\psi_{n,\boldsymbol{j},\boldsymbol{l}}$ s

$\displaystyle\begin{pmatrix}\boldsymbol{0}\\\
\psi_{n,\boldsymbol{j},\boldsymbol{l}}\end{pmatrix}=\mathcal{M}^{-1}\begin{pmatrix}\psi_{0,\boldsymbol{j},\boldsymbol{l}}\\\
\boldsymbol{0}\end{pmatrix}$ (108)

where $\boldsymbol{0}$ is a zero-vector of the same dimension as
$\mathcal{K}$.

##  Appendix C Condition number of matrix $\mathcal{M}$

For clarity we will only derive this for linear constant coefficient transport
equations. We start with the simple one-dimensional equation

$\frac{\partial}{\partial t}u+\frac{\partial}{\partial x}u=0,\quad 0<x<1,\quad
u(t,0)=g(0),$ (109)

which will be discretised by the upwind scheme

$u^{n+1}_{j}+(\lambda-1)u_{j}^{n}-\lambda u^{n}_{j-1}=0,$ (110)

where $\lambda=\Delta t/h,j=1,\cdots,N,n=0,\cdots N_{t}-1$. Now the matrix $K$
in (105) is

$\displaystyle K=\begin{bmatrix}\lambda-1&&&&\\\ -\lambda&\lambda-1&&&\\\
&\ddots&\ddots&&\\\ &&\ddots&\ddots&\\\ &&&-\lambda&\lambda-1\\\
\end{bmatrix}_{(N)\times(N_{t})}.$ (111)

Since the eigenvalues of $\mathcal{M}$ is just the singular values $\sigma$ of
$\mathcal{K}$ we now estimate the latter. Assume the CFL condition
$\lambda\leq 1$. By Gershogorin’s theorem it is easy to see
$\sigma_{\text{max}}\leq 2$. To estimate the smallest singular value, let
$L=-K$. Then

$\displaystyle\mathcal{K}^{-1}=\begin{bmatrix}I&&&&\\\ L&I&&&\\\
L^{2}&L&I&&\\\ \vdots&\vdots&\ddots&\ddots&\\\
L^{N_{t}-1}&L^{N_{t}-2}&\cdots&L&I\\\
\end{bmatrix}_{(N)\times(N_{t})}=\begin{bmatrix}I&&&&\\\ &I&&&\\\ &&I&&\\\
&&&\ddots&\\\ &&&&I\\\ \end{bmatrix}+\begin{bmatrix}&&&&\\\ L&&&&\\\ &L&&&\\\
&&\ddots&&\\\ &&&L&\\\ \end{bmatrix}+\begin{bmatrix}&&&&\\\ &&&&\\\
L^{2}&&&&\\\ &\ddots&&&\\\ &&L^{2}&&\\\ \end{bmatrix}+\cdots$ (112)

Hence

$\sigma_{\text{max}}(\mathcal{K}^{-1})=\|\mathcal{K}^{-1}\|_{2}\leq\|I\|_{2}+\|L\|_{2}+\cdots+\|L^{N_{t}-1}\|_{2}.$
(113)

By Gershigorin’s theorem,

$\eta:=\|L\|_{2}=\|K\|_{2}\leq(1-\lambda)+\lambda=1,$

therefore

$\sigma_{\text{max}}(\mathcal{K}^{-1})\leq 1+\eta+\cdots\eta^{N_{t}-1}\leq
N_{t}.$ (114)

Consequently

$\kappa(\mathcal{K})=\frac{\sigma_{\text{max}}}{\sigma_{\text{min}}}\leq
N_{t}/2.$

We next consider the $d$-dimensional equation:

$\frac{\partial}{\partial t}u+\sum_{l=1}^{d}\frac{\partial}{\partial
x_{l}}u=0,\quad 0<x_{l}<1,\quad{\text{with incoming boundary conditions,}}$
(115)

which will be discretised by the upwind scheme

$u^{n+1}_{\boldsymbol{j}}-u^{n}_{\boldsymbol{j}}+\lambda\sum_{l=1}^{d}(u_{\boldsymbol{j}}^{n}-u^{n}_{\boldsymbol{j}-e_{l}})=0,$
(116)

where $e_{l}$ is the unit vector with $1$ in the $l$-th entry and $0$
elsewhere, $\lambda=\Delta t/h,j_{l}=1,\cdots,N,1\leq l\leq N,n=0,\cdots
N_{t}-1$.

Define

$\displaystyle{u}=(u_{1,1,\cdots,1},\cdots,u_{N,1,\cdots,1},u_{1,2,\cdots,1},\cdots,u_{1,N,\cdots,1},\cdots,u_{N,N,\cdots,1},\cdots,u_{N,N,\cdots,N})^{T}$

then scheme (116) can be written in vector form as

$\displaystyle{u}^{n+1}-B{u}^{n}=0,$

where

$\displaystyle B=\lambda[T_{h}\otimes I^{\otimes(n-1)}+I\otimes T_{h}\otimes
I^{\otimes(n-2)}+\cdots+I^{\otimes(n-1)}\otimes T_{h}]+(1-d\lambda)I^{\otimes
n},$

with

$\displaystyle T_{h}=\begin{bmatrix}0&&&&\\\ 1&0&&&\\\ &\ddots&\ddots&&\\\
&&&1&0\\\ \end{bmatrix}_{N\times N}.$

Let

$\tilde{T}_{h}=(1-d\lambda)I+\lambda T_{h}$

For clarify, consider the case of $d=3$. A direction calculation shows that
the first few rows of $B$ are

$B_{1}=\begin{bmatrix}\tilde{T}_{h}&&&&&&&&\\\ \lambda
T_{h}&\tilde{T}_{h}&&&\\\ &\ddots&\ddots&&&&&&\\\ &&&\lambda T_{h}&&&&&\\\
\lambda I&&&&\tilde{T}_{h}&&&&\\\ &\lambda I&&&\lambda
T_{h}&\tilde{T}_{h}&&&\\\ &&\ddots&&&\ddots&\ddots&&\\\ &&&\lambda I&&&\lambda
T_{h}&\tilde{T}_{h}&\\\ \end{bmatrix}$

after ignoring subsequent repeating blocks. Assume the stability (CFL)
condition $3\lambda\leq 1$, then Gershgorin’s theorem implies

$\|B\|_{2}\leq[(1-3\lambda)+\lambda]+\lambda+\lambda\leq 1\,.$

The rest of the proof of the condition number of the HHL matrix will be the
same as the case of $d=1$.

Similarly, for general $d$, assume the stability (CFL) conditon

$\mathrm{d}\lambda\leq 1,$ (117)

one gets the condition number of the HHL matrix

$\kappa(\mathcal{M})\lesssim N_{t}$ (118)

while the sparsity of $\mathcal{M}$ is obviously $\sim d$.

Assume one wants to compute to time $T$, namely $N_{t}\Delta t=T$. Using the
stability condition (117), one gets

$\kappa(\mathcal{M})\lesssim N_{t}\sim dNT\,.$

##  Appendix D Proof of Lemma 12

Let $\psi^{\omega}$ be the analytical solution of (22) with initial data (23)
in which $\delta$ is replaced by $\delta_{\omega}$, and $\psi^{\omega,h}$ is
the upwind discretisation of $\psi^{\omega}$. Then

$\displaystyle\epsilon_{C}\equiv|\langle
G(t_{n},x_{\boldsymbol{j}})\rangle-G^{n}_{\omega}(\boldsymbol{j})|$
$\displaystyle\leq$
$\displaystyle\left|\int_{\Omega_{p}}G(p)(\psi(t_{n},x_{\boldsymbol{j}},p)-\psi^{\omega}(t_{n},x_{\boldsymbol{j}},p))\,dp\right|+\left|\int_{\Omega_{p}}G(p)\psi^{\omega}(t_{n},x_{\boldsymbol{j}},p))\,dp-\frac{1}{N_{p}}\sum_{\boldsymbol{l}}G(p_{\boldsymbol{l}})\psi^{\omega}_{n,\boldsymbol{j},\boldsymbol{l}}\right|$
$\displaystyle\qquad+\left|\frac{1}{N_{p}}\sum_{\boldsymbol{l}}G(p_{\boldsymbol{l}})(\psi^{\omega}_{n,\boldsymbol{j},\boldsymbol{l}}-\psi^{\omega,h}_{n,\boldsymbol{j},\boldsymbol{l}})\right|$
$\displaystyle=$ $\displaystyle I+II+III$

Denote the bicharacteristics of (22) by

$\frac{\partial x}{\partial t}=\nabla_{p}H,\quad\frac{\partial p}{\partial
t}=-\nabla_{x}H,\quad\qquad x(0)=x^{(0)},\quad p(0)=p^{(0)}$ (119)

Assume sufficient smoothness of $H$, and denoted by
$(x(t;x^{(0)},p^{(0)}),p(t;x^{(0)},p^{(0)}))$ the unique solution of the above
ODE problem, which can be inverted to get the inverse functions
$(x^{(0)}(t;x,p),p^{(0)}(t;x,p))$. Note that this is possible since the
Jacobian matrix of the map from $(x^{(0)},p^{(0)})$ to $(x(t),p(t))$ has
determinant $1$. Then by method of characteristics,

$\displaystyle I=$
$\displaystyle\left|\int_{\Omega_{p}}G(p)\left[\psi(0,x^{(0)}(t_{n},x_{\boldsymbol{j}},p),p^{(0)}(t_{n},x_{\boldsymbol{j}},p))-\psi^{\omega}(0,x^{(0)}(t_{n},x_{\boldsymbol{j}},p),p^{(0)}(t_{n},x_{\boldsymbol{j}},p))\right]\,dp\right|$
$\displaystyle=$
$\displaystyle\left|\int_{\Omega_{p}}G(p)\left[\delta(p^{(0)}(t_{n},x_{\boldsymbol{j}},p)-u_{0}(x^{(0)}(t_{n},x_{\boldsymbol{j}},p)))-\delta_{\omega}(p^{(0)}(t_{n},x_{\boldsymbol{j}},p)-u_{0}(x^{(0)}(t_{n},x_{\boldsymbol{j}},p)))\right]\,dp\right|.$

Here for notation clarity we assume $M=1$. The more general case is obtained
simply by linear superposition.

Assume
$\chi(t_{n},x_{\boldsymbol{j}},p)=p^{(0)}(t_{n},x_{\boldsymbol{j}},p)-u_{0}(x^{(0)}(t_{n},x_{\boldsymbol{j}},p))$
has (just one) root $p_{*}\in\Omega_{p}$. The case of multiple roots can be
dealt with similarly. Then

$\int_{\Omega_{p}}G(p)\delta(p^{(0)}(t_{n},x_{\boldsymbol{j}},p)-u_{0}(x^{(0)}(t_{n},x_{\boldsymbol{j}},p)))dp=\frac{G(p_{*})}{\left|\frac{\partial\chi}{\partial
p}(p_{*})\right|},$ (120)

while, by the mean value theorem,

$\displaystyle\int_{\Omega_{p}}G(p)\delta_{\omega}(p^{(0)}(t_{n},x_{\boldsymbol{j}},p)-u_{0}(x^{(0)}(t_{n},x_{\boldsymbol{j}},p)))\,dp$
(121) $\displaystyle=$
$\displaystyle\int_{\Omega_{p}}G(p)\delta_{\omega}\left(\frac{\partial\chi}{\partial
p}(p_{**})(p-p_{*})\right)\,dp=\int_{\Omega_{p}}G\left(\frac{q}{\frac{\partial\chi}{\partial
p}(p_{**}(p))+p_{*}}\right)\frac{\delta_{\omega}(q)}{\left|\frac{\partial\chi}{\partial
p}(p_{**}(p))\right|}\,dp$ (122) $\displaystyle=$
$\displaystyle\int_{|q|\leq\omega}G\left(\frac{q}{\frac{\partial\chi}{\partial
p}(p_{**})}+p_{*}\right)\frac{\delta_{\omega}(q)}{\left|\frac{\partial\chi}{\partial
p}(p_{**})\right|}\,dp=G\left(\frac{\tilde{q}}{\frac{\partial\chi}{\partial
p}(\hat{q})}+p_{*}\right)\frac{1}{\left|\frac{\partial\chi}{\partial
p}(\hat{q})\right|}\int_{|q|\leq\omega}\delta_{\omega}(q)\,dq$ (123)
$\displaystyle=$ $\displaystyle
G\left(\frac{\tilde{q}}{\frac{\partial\chi}{\partial
p}(\hat{q})}+p_{*}\right)\frac{1}{\left|\frac{\partial\chi}{\partial
p}(\hat{q})\right|}=\frac{G(p_{*})}{\left|\frac{\partial\chi}{\partial
p}(p_{*})\right|}+O(\omega),$

since $|\tilde{q}|\leq\omega,|\hat{q}-q_{*}|\leq\omega$. Combining the above
two estimates gives $I=O(\omega)$.

II is just the error of the (midpoint) quadrature rule, hence
$II=O(d(h/\omega)^{2})$, where the $1/\omega^{2}$ factor comes from the second
derivative of $\psi^{\omega}$ in the truncation error.

By standard error analysis for linear hyperbolic equation [53], if the first
order upwind scheme is used, under the CFL condition, one has

$|\psi^{\omega}_{n,\boldsymbol{j},\boldsymbol{l}}-\psi^{\omega,h}_{n,\boldsymbol{j},\boldsymbol{l}}|=O(dh/\omega^{2})$

where $1/\omega^{2}$ factor again comes from the second derivative of
$\psi^{\omega}$ in the consistency error.

By combining I, II and III the lemma is proved.

##  Appendix E Estimating normalisation constants: proof of Lemmas 14, 20 and
28

Here we prove Lemma 14.

###### Proof.

To estimate $n_{G}$, one can use $O(1)\sim\int dp|G(p)|^{2}$. This is possible
since to compute any observables we only need to consider the domain of
$\phi(t,x,p)$, which has compact support in $p$ with support of size $O(1)$
domain, thus we only need to integrate $p$ in an $O(1)$ domain. Then one can
write

$\displaystyle O(1)\sim\int
dp|G(p)|^{2}\approx\frac{1}{N^{d}}\sum_{\boldsymbol{l}}|G_{\boldsymbol{l}}|^{2}=\frac{N_{G}^{2}}{N^{d}}=n^{2}_{G},$
(124)

thus $n_{G}\sim O(1)$. For a lower bound of $n^{2}_{\psi_{0}}$, we use
$\Upsilon\leq 1$ since we estimate $\Upsilon\approx\langle
0|\mathcal{U}|0\rangle$ from a large unitary matrix $\mathcal{U}$ (see
Appendix H) and can therefore view it as a quantum fidelity. Since the
ensemble average is an observable and should be also independent of $N$ and we
assume to be independent of $d$, then $O(1)\sim
G^{\omega}_{n,\boldsymbol{j}}$. Together with $\Upsilon\leq 1$ and $n_{G}\sim
O(1)$, we have $n_{\psi_{0}}\gtrsim O(1)$.

To obtain an upper bound of $n_{\psi_{0}}$, one considers the case where there
are no assumption on the initial data. One can estimate $N_{\psi_{0}}$ with
$\psi(0,x,p)$ defined by Eq. (23) with the delta function approximated by
$\delta_{\omega}$. We start with $M=1$. Without loss of generality we assume
$(p,x)\in[0,1]^{2d}$ and

$\displaystyle\int\int|\psi(0,x,p)|^{2}dp\,dx\approx\int\int\delta_{\omega}(p-u^{[k]})^{2}dp\,dx\leq\frac{1}{\omega^{d}}\int\int\delta_{\omega}(p-u^{[k]})dp\,dx$
$\displaystyle=$ $\displaystyle\frac{1}{\omega^{d}}\int
1\,dx=\frac{1}{\omega^{d}}.$ (125)

To obtain an upper bound of $n_{\psi_{0}}$, one considers the case where there
are no assumptions on the initial conditions. One can estimate $N_{\psi_{0}}$
with $\psi(0,x,p)$ defined by Eq. (23) with the delta function approximated by
$\delta_{\omega}$. We start with $M=1$. Without loss of generality we assume
$(p,x)\in[0,1]^{2d}$ and

$\displaystyle\int\int|\psi(0,x,p)|^{2}dp\,dx\approx\int\int\delta_{\omega}(p-u^{[k]})^{2}dp\,dx\leq\frac{1}{\omega^{d}}\int\int\delta_{\omega}(p-u^{[k]})dp\,dx$
$\displaystyle=$ $\displaystyle\frac{1}{\omega^{d}}\int
1\,dx=\frac{1}{\omega^{d}}$ (126)

Next we approximate the above integral by the quadrature rule. Here we use
$N^{d}$ mesh points in $x$. Note $\delta_{\omega}(p-u^{[k]})$ has support of
the size of only $O(\omega^{d})=O((mh)^{d})=O((m/N)^{d})$ in $p$, thus the
number of mesh points in $p$ is of order $m^{d}$, where $m$ is the number of
non-zero entries in the discrete delta function $\delta_{\omega}$. Then using
the quadrature rule,

$\int\int|\psi(0,x,p)|^{2}dp\,dx\approx\frac{1}{(mN)^{d}}\sum_{\boldsymbol{j}}^{N}\sum_{\boldsymbol{l}}^{N}|\psi_{0,\boldsymbol{j},\boldsymbol{l}}|^{2}=\frac{1}{(mN)^{d}}N_{\psi_{0}}^{2}.$
(127)

From Eqs. (132) and (127) one gets

$N_{\psi_{0}}=O\left(\frac{(mN)^{d/2}}{\omega^{1/2}}\right)=O\left(\frac{(mN)^{d/2}}{(m/N)^{d/2}}\right)=O\left(N^{d}\right)$
(128)

More generally, suppose the initial data has support in a box $\Omega_{\beta}$
of size $\beta$. This means one can write

$\int_{\Omega_{\beta}}\int|\psi(0,x,p)|^{2}dp\,dx\leq\frac{1}{\omega^{d}}\int_{\Omega_{\beta}}1\,dx=\frac{\beta^{d}}{\omega^{d}}$
(129)

while its numerical approximation is

$\int\int|\psi(0,x,p)|^{2}dp\,dx\approx\frac{1}{(m\beta
N)^{d}}\sum_{\boldsymbol{j}}^{N}\sum_{\boldsymbol{l}}^{N}|\psi_{0,\boldsymbol{j},\boldsymbol{l}}|^{2}=\frac{1}{(m\beta
N)^{d}}N_{\psi_{0}}^{2}.$ (130)

Eqs. (129) and (130) give

$N_{\psi_{0}}=O\left((\beta N)^{d}\right).$ (131)

We recover the upper bound in the limit $\beta=O(1)$. The case of $M>1$ will
yield the same estimate since it is just the average of $M$ different initial
data of the same type. ∎

The proof of Lemma 20 proceeds in exactly the same way for $n_{g}$ and the
lower bound of $n_{\Phi_{0}}$. For more general cases, there is a small
modification in the proof.

###### Proof.

The estimate of $N^{\prime}_{\psi_{0}}$ is pretty much the same as that in
Lemma 14, except here $p\in\mathbb{R}^{1}$. Thus

$\displaystyle\int_{\Omega_{\beta}}\int|\psi(0,x,p)|^{2}dp\,dx\approx\int\int\delta_{\omega}(p-u^{[k]})^{2}dp\,dx\leq\frac{1}{\omega}\int\int\delta_{\omega}(p-u^{[k]})dp\,dx$
$\displaystyle=$
$\displaystyle\frac{1}{\omega}\int_{\Omega_{\beta}}1\,dx=\frac{\beta^{d}}{\omega},$
(132)

while its numerical quadrature approximation is

$\int\int|\psi(0,x,p)|^{2}dp\,dx\approx\frac{1}{m(\beta
N)^{d}}\sum_{\boldsymbol{j}}^{N}\sum_{\boldsymbol{l}}^{N}|\psi_{0,\boldsymbol{j},\boldsymbol{l}}|^{2}=\frac{1}{m(\beta
N)^{d}}(N^{\prime}_{\psi_{0}})^{2}.$ (133)

Eqs. (132) and (133) give

$N_{\psi_{0}}=O\left(\beta^{d}N^{(d+1)/2}\right),$ (134)

so the upper bound when $\beta=O(1)$ is
$N_{\psi_{0}}=O\left(N^{(d+1)/2}\right)$. ∎

The proof of Lemma 28 for estimating $n_{A}$ and the lower bound of
$n_{\Phi_{0}}$ is identical to the estimation of $n_{G}$ and the lower bound
to $n_{\psi_{0}}$, respectively, in Lemma 14. The proof for the upper bound of
$n_{\Phi_{0}}$ is the following.

###### Proof.

The estimate of $N_{\Phi_{0}}$ is pretty much the same as that in Lemma 14,
except that here there is no $x$ and $q\in\mathbb{R}^{D}$ thus

$\int|\Phi(0,q)|^{2}dq\approx\int\delta_{\omega}(q-X_{0}^{[k]})^{2}dq\leq\frac{1}{\omega^{D}}\int\delta_{\omega}(q-u^{[k]})dq=\frac{1}{\omega^{D}},$
(135)

while the quadrature rule becomes

$\int|\Phi(0,q)|^{2}dq\approx\frac{1}{m^{D}}\sum_{\boldsymbol{j}}^{N}|\Phi_{0,\boldsymbol{j}}|^{2}\\\
=\frac{1}{m^{D}}N_{\Phi_{0}}^{2}.$ (136)

Consequently $N_{\Phi_{0}}=O(N^{D/2})$. ∎

##  Appendix F Algorithm for computing $\Upsilon$ and proof of Lemma 16

Here we prove Lemma 16 which is used for computing
$\Upsilon\equiv\langle\psi_{0}|(\mathcal{M}^{-1})^{\dagger}\mathcal{G}\mathcal{M}^{-1}|\psi_{0}\rangle$.
We attempt to keep this as self-contained as possible, so we include in full
similar arguments and steps made in [36], but with modifications required due
to our assumption of beginning with sparse access
$(s,\|\mathcal{M}\|_{max},O_{M},O_{F})$ to $\mathcal{M}$, instead of block
access to $\mathcal{M}$. We also do not begin with block access to
$\mathcal{G}$. However, block access to $\mathcal{M}$ can be constructed from
sparse access from the following lemma.

######  Lemma 32.

([36, 42]) Given $m$ qubits and sparse access
$(s,\|\mathcal{M}\|_{max},O_{M},O_{F})$ to a $m$-qubit Hermitian matrix
$\mathcal{M}$, block access
$(\alpha_{\mathcal{M}}=s\|\mathcal{M}\|_{max},n_{\mathcal{M}}=2,\delta_{\mathcal{M}},U_{\mathcal{M}})$
to $\mathcal{M}$ can be implemented by making $\mathcal{O}(1)$ queries to
sparse access oracles for $\mathcal{M}$ and
$\mathcal{O}(m+\log^{2.5}(s\|\mathcal{M}\|_{max}/\delta_{\mathcal{M}}))$
additional $2$-qubit gates.

###### Proof.

See Lemma II.5 in [36] and also Lemma 6 in [42] for an explicit construction.
∎

This means one cannot construct a block access to $\mathcal{M}$ with
$\delta_{\mathcal{M}}=0$ using a finite number of $2$-qubit gates if one
begins with sparse access to $\mathcal{M}$. Thus one cannot apply Lemma 8 from
[36] directly at this point and needs to extend to the case where one is
instead given block access
$(\alpha_{\mathcal{M}},n_{\mathcal{M}},\delta_{\mathcal{M}}>0,U_{\mathcal{M}})$
to $\mathcal{M}$ and block access
$(\alpha_{\mathcal{G}}=1,n_{\mathcal{G}},0,U_{\mathcal{G}})$ to $\mathcal{G}$.
In our scenario, $\mathcal{G}$ is a pure density matrix, which means it is
simple to construct a $\delta_{\mathcal{G}}=0$ block access if only given
access to unitaries that construct the state, using the following lemma.

######  Lemma 33.

([54]) Given $\mathcal{G}=|G_{n,\boldsymbol{j}}\rangle\langle
G_{n,\boldsymbol{j}}|$ is a $m$-qubit density matrix and $L(n,\boldsymbol{j})$
is a $m$-qubit unitary operator such that
$L(n,\boldsymbol{j})|0^{m}\rangle=|G\rangle$, then
$U_{\mathcal{G}}=(L(n,\boldsymbol{j})^{\dagger}\otimes\mathbf{1}^{\otimes
m+1})(\mathbf{1}^{\otimes m+1}\otimes
SWAP_{m})(L(n,\boldsymbol{j})\otimes\mathbf{1}^{m+1})$ is a
$(1,n_{\mathcal{G}}=2m+1,\delta_{\mathcal{G}}=0,U_{\mathcal{G}})$ block
encoding of $\mathcal{G}$.

###### Proof.

Apply Lemma 45 in [54] to the case of a pure state $\mathcal{G}$. Note that
the SWAPm gate can be constructed from $O(m)$ 2-qubit gates. ∎

The block access
$(\alpha_{\mathcal{M}},n_{\mathcal{M}},\delta_{\mathcal{M}},U_{\mathcal{M}})$
can then be used to construct a block-access
$(\alpha_{\mathcal{M}^{-1}},n_{\mathcal{M}^{-1}},\delta,U_{\mathcal{M}^{-1}})$
to $\mathcal{M}^{-1}$ as seen from Lemma 38. The block access
$(\alpha_{\mathcal{G}}=1,n_{\mathcal{G}},0,U_{\mathcal{G}})$ can be used to
construct block access $(1,n_{\mathcal{G}}+1,0,U_{\mathcal{G}^{\prime}})$ to
$\mathcal{G}^{\prime}=|0^{n_{\mathcal{M}^{-1}}}\rangle\langle
0^{n_{\mathcal{M}^{-1}}}|\otimes\mathcal{G}$, which will be shown in Lemma 34.

######  Lemma 34.

([36]). Given block-access $(1,n_{\mathcal{G}},0,U_{\mathcal{G}})$ to
$\mathcal{G}$, block-access $(1,n_{\mathcal{G}}+1,0,U_{\mathcal{G}^{\prime}})$
to $\mathcal{G}^{\prime}=|0^{n_{\mathcal{M}^{-1}}}\rangle\langle
0^{n_{\mathcal{M}^{-1}}}|\otimes\mathcal{G}$ can be constructed where
$U_{\mathcal{G}}$ is queried once with $\mathcal{O}(n_{\mathcal{M}^{-1}}^{2})$
additional $2$-qubit gates.

###### Proof.

See Lemma IV.16 in [36] for an explicit construction. ∎

Given these ingredients, we are now ready to estimate $\Upsilon$ to error