# Quantum algorithms for uncertainty quantification: application to partial
differential equations

Francois Golse  Ecole Polytechnique, CMLS, 91128 Palaiseau Cedex, France.  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

Most problems in uncertainty quantification, despite its ubiquitousness in
scientific computing, applied mathematics and data science, remain formidable
on a classical computer. For uncertainties that arise in partial differential
equations (PDEs), large numbers $M>>1$ of samples are required to obtain
accurate ensemble averages. This usually involves solving the PDE $M$ times.
In addition, to characterise the stochasticity in a PDE, the dimension $L$ of
the random input variables is high in most cases, and classical algorithms
suffer from curse-of-dimensionality. We propose new quantum algorithms for
PDEs with uncertain coefficients that are more efficient in $M$ and $L$ in
various important regimes, compared to their classical counterparts. We
introduce transformations that transfer the original $d$-dimensional equation
(with uncertain coefficients) into $d+L$ (for dissipative equations) or $d+2L$
(for wave type equations) dimensional equations (with certain coefficients) in
which the uncertainties appear only in the initial data. These transformations
also allow one to superimpose the $M$ different initial data, so the
computational cost for the quantum algorithm to obtain the ensemble average
from $M$ different samples is then independent of $M$, while also showing
potential advantage in $d$, $L$ and precision $\epsilon$ in computing ensemble
averaged solutions or physical observables.

##  I Introduction

Most physical models are not based on first-principles and are thus subject to
uncertainties due to modeling or measurement errors. Examples include
constitutive relations or equations of state in continuum mechanics,
scattering coefficients in wave propagation and transport, initial or boundary
data, forcing or source terms, diffusivity in porous or heterogeneous media.
Quantifying these uncertainties are important to validate, verify and
calibrate the models, to conduct risk management, and to control the
uncertainties [1, 2]. In the last two decades the area of uncertainty
quantification (UQ) was one of the most active areas of research in scientific
computing, applied mathematics and data sciences.

Uncertainties are typically modelled stochastically. In the case of ordinary
and partial differential equations (ODEs and PDEs), one often uses stochastic
ODEs and PDEs. To solve such models the frequently-used method is the Monte-
Carlo method, which does not suffer from the curse-of-dimensionality. Its
drawback, however, is its low–in fact only halfth- order accuracy and
statistical noise. Therefore a large number of samples $M$ is needed for
simulation so their ensemble averages are accurate. Higher order methods have
been introduced more recently, for example the stochastic Galerkin and
stochastic collocation methods [3, 4], that are attractive when the solution
has good regularity in the random space.

In most cases the dimension of the random space is high. For example, the
stochastic process is often approximated by the Karhunen-Loeve expansion [5],
which is a linear combination of a large number of i.i.d. random variables
(referred to as random inputs). Hence the dimension of the random input
variables could be high and one needs to needs to solve the underlying
uncertain differential equations in very high-dimensional parameter space, on
top of the possibly already high dimensionality in the physical space. Sparse
grid [6] or greedy algorithms [7] have been used but they are not effective in
high dimensions. They also need sufficient regularity requirements on the
parameter space and the accuracy could still be just marginally higher than
the Monte-Carlo methods. Therefore, most of UQ problems remain formidable for
a classical computer.

On the other hand, quantum algorithms based on computations on qubits, could
offer possibly up to exponential speedup over their classical counterparts,
although these are also subject to many caveats. While quantum computers that
can solve practical problems are still out of reach in the near future, the
development of quantum algorithms that show quantum advantages over classical
algorithms have seen increasing activities in recent decades, for instance in
linear algebra problems [8, 9, 10]. Quantum algorithms for certain ODEs and
PDEs (linear and nonlinear) have been proposed which, in certain regimes,
could also demonstrate quantum advantages in dimension, precision, and the
number of simulations [11, 12, 13, 14, 15, 16, 17, 18]. Here one usually
discretizes the equations in space and time, and then formulate them into
linear algebra problems which are then solved by quantum algorithms for linear
systems of equations.

In this paper we propose quantum algorithms for some of the most important
linear PDEs with uncertain coefficients. We aim to develop efficient quantum
algorithms that allow a large number of simulations for multiple initial
and/or boundary data, in particular those arising in sampling methods, such as
Monte-Carlo methods or stochastic collocation methods for uncertain PDEs. In
these simulations one has to run repeated simulations for a large number of
samples, or initial and boundary data in order to obtain ensemble averages of
the solutions.

Our main idea is to use transformations that transfer the equation to phase
space, so the random coefficients appear only in the initial data, not in the
equation itself. This allows us to run the simulation with multiple initial
data just once to get the ensemble average of the solution. This means that we
gain speedup with respect to $M$ – the (possibly very large) number of samples
or initial data. While the original transformation was introduced in [19], it
was used for the homogenization of PDEs with oscillatory coefficients. Here we
use it for a completely different purpose and we also generalize the
transformation for our purpose in the case of more general heterogeneous
uncertain coefficients. The resulting phase space equations are then defined
in higher dimensions. This means it can no longer be efficiently solved on a
classical device, but can be efficiently solved by quantum linear PDE solvers
with quantum advantage with respect to several critical parameters including
physical dimension. Most importantly, by transferring the stochastic
parameters into the initial conditions, one only needs to solve a single
instance of the PDE to arrive at the ensemble average of the solutions or
quantities of interest.

Throughout this paper, when we say quantum advantage, we refer to the
comparison to canonical classical algorithms that are used in practice – which
have good numerical accuracy and high resolution – to solve these problems.
Here we use $r^{\text{th}}$-order finite difference approximations.

We will use some of the most important linear PDEs to illustrate our approach,
and then compare with the cost of both classical algorithms and standard
quantum algorithms that simulate directly the original problems to demonstrate
the new advantage gained. Among the equations we will study include the heat
equation, the linear Boltzmann equation, the linear advection equation, and
the linear Schrödinger equation. The treatment of the first two equations–both
are dissipative–will be different than the last two which are wave type
equations. The coefficients that contain uncertainties is assumed to have the
following form:

$a(x,z)=\sum_{i=1}^{L}a_{i}(z)b_{i}(x)$ (1)

where $x\in\mathbb{R}^{d}$ is the space variable, and $z$ is the (possibly
high-dimensional) random or deterministic parameters that model uncertainties.
Without loss of generality, we assume $a_{i}(z)>0$ for all $i$ (since
otherwise one can absorb the negative sign into $b_{i}(x)$), and

$\sum_{i=1}^{L}a^{2}_{i}(z)\leq C$ (2)

for some $C$ independent of $L$. This already covers a fairly general class of
inhomogeneous and uncertain coefficients. This is because one can often
approximate a general inhomogeneous coefficient $a(x,z)$ by a linear
combination of basis functions in a suitable space, for example in $x$ with an
orthogonal basis in $L^{2}$ (like in the Galerkin approximation in space), or
in $z$ through polynomial chaos expansion or Karhunen-Loeve expansion of a
random or stochastic process [5, 2, 3], which gives an approximation in the
form of Eq. (1).

In each of these cases we will estimate the query and gate complexities of the
quantum algorithms to compute ensemble averaged physical observables, where
the cost can be independent of $M$. For large $M$, advantage with respect to
$L$ is possible. In certain regimes, the new quantum algorithms also offer
potential speedup with respect to dimension $d$ and precision $\epsilon$.

In Section II we briefly review quantum algorithms for the system of linear
equations which are important subroutines for solving ODEs and PDEs. In
Section III we introduce new transformations for linear PDEs with uncertainty
– shown explicitly for the heat equation, the linear Boltzmann equation, the
linear advection equation, and the linear Schrödinger equation – so the
equations become deterministic in higher dimensional phase space. For each
equation, we show the end-to-end quantum costs (including initial state
preparation and final measurement costs) for computing physical quantities of
interest in these examples and show when quantum advantage can be expected.
Our main results are summarised in Table 1.

Table 1: Quantum ($\mathcal{Q}$) and classical ($\mathcal{C}$) cost comparison
for $r^{\text{th}}$-order approximations in computing ensemble averaged
solutions, at $\Lambda$ final meshpoints, over $M$ samples, when
$n_{0}^{2}\Lambda=O(N^{b})$, where $n_{0}^{2}$ is a normalisation of the
initial state. Quantum advantages are possible when $\gamma_{i}>0$ for the
corresponding parameters. We also give the sufficient range of $b$ where
quantum advantage is possible. In the table $M_{heat}\equiv
O(L^{2+(d+L+3)/c}(d/\epsilon)^{(L+1)/3})$, $M_{Boltz}\equiv
O(L^{(2d+L+1)/c}\max(L,d)(d/\epsilon)^{L/c}/d)$, $M_{adv}\equiv
O(L^{(d+2L+3)/c+2}(d/\epsilon)^{(2L+2)/c})$ and $M_{Schr}\equiv
O((d+L)^{(d+2L+2)/c+2}/(d^{(d+2)/c+2}\epsilon^{2L/r}))$. See text in Section
III for details.

$\displaystyle\mathcal{O}\left(\frac{\mathcal{C}}{\mathcal{Q}}\right)=\tilde{O}\left(M^{\gamma_{1}}d^{\gamma_{2}}L^{\gamma_{3}}(d+L)^{\gamma_{4}}\left(\frac{1}{\epsilon}\right)^{\gamma_{5}}\right)$

$(d+1)$-dim PDE $\gamma_{1}$ $\gamma_{2}$ $\gamma_{3}$ $\gamma_{4}$
$\gamma_{5}$ $b$ range Parameters ($M$ initial data) with advantage (possible)
Linear heat $M<M_{heat}$ $1$ $\frac{d-7-b}{r}-1$ $-4-\frac{9+b}{r}$ $0$
$\frac{d-7-b}{r}$ $[0,d-7-r]$ $M,d,\epsilon$ $M>M_{heat}$ $0$
$\frac{d+L-6-b}{r}-2$ $\frac{d+L-6-b}{3}$ $0$ $\frac{d+L-6-b}{r}-1$
$[0,d+L-6-2r]$ $L,d,\epsilon$ Linear Boltzmann $M<M_{Boltz}$ $1$
$\frac{d-2-b}{r}+1$ $-1-\frac{3+d+b}{r}$ $-3$ $\frac{d-2-b}{r}-1$ $[0,d-2-2r]$
$M,d,\epsilon$ $M>M_{Boltz},L>d$ $0$ $\frac{d+L-2-b}{r}$ $\frac{d+L-2-b}{r}$
$-3$ $\frac{d+L-2-b}{r}-1$ $[0,d+L-2-2r]$ $L,d,\epsilon$ $M>M_{Boltz},L<d$ $0$
$\frac{d+L-2-b}{r}+1$ $\frac{d+L-2-b}{r}-1$ $-3$ $\frac{d+L-2-b}{r}-1$
$[0,d+L-2-2r]$ $L,d,\epsilon$ Linear advection $M<M_{adv}$ $1$
$\frac{d-8-b}{r}-2$ $-4-\frac{9}{r}$ $0$ $\frac{d-8-b}{r}-1$ $[0,d-8-2r]$
$M,d,\epsilon$ $M>M_{adv}$ $0$ $\frac{d+2L-6-b}{r}-2$ $\frac{d+2L-6-b}{r}-2$
$0$ $\frac{d+2L-6-b}{r}-1$ $[0,d+2L-6-2r]$ $L,d,\epsilon$ Schrödinger
$M<M_{Schr}$ $1$ $\frac{d+2}{r}+2$ $0$ $-4-\frac{6+b}{r}$ $\frac{d-4-b}{r}-1$
$[0,d-4-r]$ $M,d,\epsilon$ $M>M_{Schr}$ $0$ $0$ $0$ $\frac{d+2L-4-b}{r}-2$
$\frac{d+2L-4-b}{r}-1$ $[0,d+2L-4-2r]$ $L,d,\epsilon$

##  II Quantum algorithms for systems of linear equations

Numerical methods to solve linear ODEs and PDEs, with solutions at the mesh
points and/or time steps given by a vector $x$, can be written as a system of
linear algebraic equations $\mathcal{M}x=y$. Without loss of generality, we
can assume the matrix $\mathcal{M}$ to be Hermitian (since a general square
matrix can be made Hermitian through dilation), with the $(i,j)$-th entry
denoted $\mathcal{M}_{ij}$. Assume
$\|\mathcal{M}\|_{max}=\max_{i,j}(|\mathcal{M}_{ij}|)\leq 1$ (otherwise one
can divide $\mathcal{M}$ by a constant so this condition is satisfied). Its
condition number is $\kappa$, and sparsity (the number of non-zero entries in
each row and column) is $s$. Although quantum subroutines do not output all
the solutions $x$, many quantum subroutines exist to solve an alternative
problem. This is the quantum linear systems problem (QLSP) that outputs the
quantum state $|x\rangle$ instead of the full classical solution $x$.

Assume vectors $x$ and $y$ have elements $\\{x_{i}\\}$, $\\{y_{i}\\}$. One can
then define the following $m$-qubit quantum states
$|x\rangle\equiv\frac{1}{\mathcal{N}_{x}}\sum_{i}x_{i}|i\rangle$,
$|y\rangle\equiv\frac{1}{\mathcal{N}_{y}}\sum_{i}y_{i}|i\rangle$ where
$\mathcal{N}_{x}=\sqrt{\sum_{i}|x_{i}|^{2}}=\|x\|$,
$\mathcal{N}_{y}=\sqrt{\sum_{i}|y_{i}|^{2}}=\|y\|$ are normalisation
constants.

######  Problem 1.

(QLSP) Let $\mathcal{M}$ be a $2^{m}\times 2^{m}$ Hermitian matrix with
spectral norm $\|\mathcal{M}\|\leq 1$ with condition number $\kappa$. Assume
vectors $x$ and $y$ solve $\mathcal{M}x=y$. 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$.

In Problem 1, ‘access to $\mathcal{M}$’ is defined with respect to an oracle
that can access the entries of $\mathcal{M}$. There are primarily two types of
oracles considered, called sparse-access and block-access respectively,
defined below. The query complexity of the algorithm is defined with respect
to the number of times such an oracle is used during the protocol.

######  Definition 1.

Sparse access to a Hermitian matrix $\mathcal{M}$ with sparsity $s$ is defined
to be the $4$-tuple $(s,\|\mathcal{M}\|_{max},O_{\mathcal{M}},O_{F})$. 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$ (3)

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

For instance, the common subroutines for QLSP that use sparse-access are the
HHL algorithm [8] and the linear combination of unitaries (LCU) algorithms
[15], including those versions that use the variable-time amplitude
amplification algorithm (VTAA) [20]. The HHL algorithm has query complexity
$\tilde{\mathcal{O}}(s\kappa^{2}\/\eta)$ and the LCU algorithms have the
complexity $\tilde{\mathcal{O}}(s\kappa^{2}\text{poly}\log(1/\eta))$, where
VTAA can be used to reduce each complexity by a factor of $\kappa$.

On the other hand, more recent algorithms like those preparing $|x\rangle$ via
the quantum adiabatic theorem [21, 22], which demonstrate improved scaling
with respect to conditions number $\kappa$ of $\mathcal{M}$ and error $\eta$,
often use block-access. The optimal scaling is achieved in [22] with scaling
$\mathcal{O}(\kappa\log(1/\eta))$. Since sparse-access to $\mathcal{M}$ can be
used to construct block-access to $\mathcal{M}$ [23, 24], query complexity
with respect to block-access can be converted into sparse-access results with
an extra factor of up to $s$ and some small overhead.

######  Definition 2.

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}}.$
(4)

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.

However, for a quantum algorithm to be useful for PDE problems in real
applications, it should output values of meaningful physical quantities of
interest, instead of quantum states $|x\rangle$. Thus it is important to solve
instead expectations values with respect to $|x\rangle$.

######  Problem 2.

Given a Hermitian matrix $\mathcal{G}$, which is of the same size as
$\mathcal{M}$, access to $\mathcal{M}$ and $U_{initial}$, the aim is to
compute the expectation value $\boldsymbol{x}^{T}\mathcal{G}\boldsymbol{x}$ to
precision $\epsilon$.

We note that here the error $\epsilon$ in the expectation value is different
to $\eta$, which is the error in the normalised quantum state itself.

With respect to sparse access, Problem 2 can be solved using the quantum
singular value decomposition [10] algorithm, which has the following query and
gate complexities.

######  Lemma 3.

[18, 24] 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_{K},\boldsymbol{j}/N)$
where $L(t_{K},\boldsymbol{j}/N)|0\rangle=|G_{K,\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=|y\rangle$. The algorithm then
returns with probability at least $2/3$ an $\epsilon$-additive approximation
to $\boldsymbol{x}^{T}\mathcal{G}\boldsymbol{x}$ where
$\mathcal{G}\equiv|G_{K,\boldsymbol{j}}\rangle\langle G_{K,\boldsymbol{j}}|$.
This algorithm makes
$\mathcal{O}(\kappa^{2}\mathcal{N}_{y}^{2}/(\|\mathcal{M}\|\epsilon))$ queries
to $U_{\mathcal{G}}$ and $U_{initial}$,
$\tilde{\mathcal{O}}(s\|\mathcal{M}\|_{max}\kappa^{3}\mathcal{N}_{y}^{2}/(\|\mathcal{M}\|\epsilon))$
queries to sparse oracles for $\mathcal{M}$ and
$\tilde{\mathcal{O}}(\kappa^{2}\mathcal{N}_{y}^{2}(m+s\|\mathcal{M}\|_{max}\kappa)/(\|\mathcal{M}\|\epsilon))$
additional $2$-qubit gates.

This scaling can be further improved with respect to $\kappa$, by employing an
alternative set of algorithms based on quantum adiabatic computation. Given
block-access to $\mathcal{M}$, an optimal scaling in terms of both $\kappa$
and $\epsilon$ can for instance be achieved by combining quantum algorithms
for computing observables via block-encoding [25] (optimal scaling in
$\epsilon$) and the optimal linear systems solver via the discrete adiabatic
theorem [22] (optimal scaling in $\kappa$ and $\epsilon$).

######  Lemma 4.

A quantum algorithm can be constructed that takes the following inputs: (i)
$(1,n_{\mathcal{M}},0,U_{\mathcal{M}})$ block access to a $2^{m}\times 2^{m}$
invertible Hermitian matrix $\mathcal{M}$; (ii)
$(1,n_{\mathcal{G}},0,U_{\mathcal{G}})$ block access to $2^{m}\times 2^{m}$
density matrix $\mathcal{G}$; (iii) an accuracy $\epsilon$; and (iv) $m$-qubit
unitary black box $U_{initial}$ where $U_{initial}|0\rangle=|y\rangle$. The
algorithm then returns with constant probability an $\epsilon$-additive
approximation to $\boldsymbol{x}^{T}\mathcal{G}\boldsymbol{x}$ by making
$\tilde{\mathcal{O}}(\mathcal{N}_{x}^{2}\kappa/(||\mathcal{M}||\epsilon))$
queries to $U_{\mathcal{M}}$, $U_{\mathcal{G}}$ and $U_{initial}$.

###### Proof.

Suppose we are given a $(1,n_{\mathcal{M}},0,U_{\mathcal{M}})$ block-encoding
of the matrix $\mathcal{M}$. Then results of [21, 22] show that with
$\tilde{\mathcal{O}}(\kappa/||\mathcal{M}||\log(1/\eta))$ query calls to
$U_{\mathcal{M}}$ and $U_{initial}$, one can prepare the state
$|\tilde{x}\rangle$ that is $\eta$-close to the state
$|x\rangle\propto\mathcal{M}^{-1}|y\rangle$. One can then adapt the algorithm
from [25] (Lemma 5) to compute the expectation $\langle
x|\mathcal{G}|x\rangle$ with amplitude estimation. Here one can produce, with
probability $1-\delta$, an $\epsilon^{\prime}$ estimate of
$\langle\tilde{x}|\mathcal{G}|\tilde{x}\rangle$ while using
$\mathcal{O}(\log(1/\delta)/\epsilon^{\prime})$ queries to the preparation of
$|\tilde{x}\rangle$ with quantum adiabatic computation and calls to
$U_{\mathcal{G}}$. Combining these results, we obtain a
$\epsilon^{\prime}$-additive approximation to $\langle x|\mathcal{G}|x\rangle$
with $\tilde{\mathcal{O}}(\kappa/\|\mathcal{M}\|\epsilon^{\prime})$ queries to
$U_{\mathcal{M}}$, $U_{\mathcal{G}}$ and $U_{initial}$, where one can choose
any constant $\delta>1/2$. Since
$\boldsymbol{x}^{T}\mathcal{G}\boldsymbol{x}=\mathcal{N}_{x}^{2}\langle
x|\mathcal{G}|x\rangle$, then $\epsilon^{\prime}=\mathcal{N}_{x}^{2}\epsilon$
and we have our result. Since $\mathcal{N}_{x}$ may not be assumed to be known
in advance, unlike $\mathcal{N}_{y}$, another quantum algorithm with cost
$\tilde{\mathcal{O}}(\kappa\mathcal{N}_{x}/\epsilon^{\prime\prime})$ [26, 27]
is required to compute $\mathcal{N}_{x}$ to precision
$\epsilon^{\prime\prime}$. ∎

######  Lemma 5.

[25] Given a density matrix $\mathcal{G}$ and let
$|\rho\rangle=U_{\rho}|0\rangle$ where $U_{\rho}$ consists of $R$ elementary
gates. Then for every $\epsilon$, $\delta>0$, there exists a quantum algorithm
that estimates $\langle\rho|\mathcal{G}|\rho\rangle$ to precision $\epsilon$
with probability at least $1-\delta$. This algorithm has gate complexity
$O((R/\epsilon)\log(1/\delta))$.

Unlike the protocol in Lemma 3, the protocol in Lemma 4 requires the
normalisation of the solution state $\mathcal{N}_{x}$ instead of the initial
state normalisation $\mathcal{N}_{y}$. Since the classical values of the
initial state is known, we can assume $\mathcal{N}_{y}$ to be given. We can
bound $\mathcal{N}_{x}$ using
$\mathcal{N}_{x}=\|\boldsymbol{x}\|=\|\mathcal{M}^{-1}\boldsymbol{y}\|\leq\|\mathcal{M}^{-1}\|\|\boldsymbol{y}\|=\kappa\mathcal{N}_{y}/\|\mathcal{M}\|$
since $\|\mathcal{M}\|\|\mathcal{M}^{-1}\|=\kappa$. This means that by taking
the worst-case scenario Lemma 4 still requires a cost that scales like Lemma 3
with respect to $\kappa$ and $\epsilon$. Thus, to give the most conservative
costs in this paper, and also to avoid implementing an extra quantum algorithm
to compute $\mathcal{N}_{x}$, it is sufficient for us to use the protocol in
Lemma 3.

We note that the protocols above all assume the preparation of the initial
state $|y\rangle$. The total gate complexity in an end-to-end quantum
algorithm must also include the gate complexity in the preparation of the
initial state. The gate complexity for preparing $|y\rangle$ can be greatly
reduced when the state is sparse. So if $|y\rangle$ has sparsity $\sigma$
(i.e. vector $y$ has $\sigma$ non-zero entries), its deterministic preparation
can have the following gate complexity

######  Lemma 6.

[28] A circuit producing an $m$-qubit state $|y\rangle$ from $|0\rangle$ with
given classical entries $\\{y_{i}\\}$ can be implemented using
$\mathcal{O}(m\sigma)$ CNOT gates and $\mathcal{O}(\sigma(\log\sigma+m))$ one-
qubit gates, where the specification of the circuit can be found with a
classical algorithm with run-time $\mathcal{O}(m\sigma^{2}\log\sigma)$.

We note that the above protocol uses a constant number of ancilla qubits.
Alternative preparation strategies are proposed in [29], which has a reduced
gate complexity $\Theta(\log(m\sigma))$, but require
$\mathcal{O}(m\sigma\log\sigma)$ ancilla qubits.

##  III Linear PDEs with uncertainty

In this section, we demonstrate our new algorithm for the heat equation, the
linear Boltzmann equation, the linear advection equation and the linear
Schrödinger equation, all with uncertain coefficients of the form in Eq. (1).

###  III.1 Linear heat equation

Consider the following initial value problem of the linear heat equation

$\begin{cases}\partial_{t}u-a(x,z)\Delta u=0\,\\\
u(0,x,z)=u_{0}(x,z)\end{cases}$ (5)

where $u=u(t,x,z)$, $x\in\mathbb{R}^{d}$ is the position, $t\geq 0$ is the
time, $a(x,z)>0$ is heat conductivity given by Eq. (1), due to the
heterogeneity of the media or background (for example in porous media [30]).
To perform a Monte-Carlo simulation one needs to select a large number (say
$M\gg 1$) of samples in $z$, solve the PDE system $M$ times and then take the
ensemble average. However, solving the equation $M$ times for large $M$ is
computationally expensive. The computation is also costly for large $L$.

We aim instead to find a transformation to another linear PDE with certain
coefficients and only a single initial condition, that still enables us to
compute ensemble averaged observables corresponding to the original PDE.

Let $p=(p_{1},\cdots,p_{L})^{T}$, with $p_{i}\in(-\infty,\infty)$ for all
$i=1,...,L$. We introduce the transformation

$\displaystyle
U(t,x,z,p)=\frac{1}{2}\prod_{i=1}^{L}a_{i}(z)e^{-a_{i}(z)|p_{i}|}u(t,x,z),$
(6)

from which one can recover $u$ from $U$ via

$\displaystyle
u(t,x,z)=\int_{(-\infty,\infty)^{L}}U(t,x,z,p)dp=2\int_{(0,\infty)^{L}}U(t,x,z,p)dp.$
(7)

A simple computation shows that $U$ solves

$\partial_{t}U+\sum_{i=1}^{L}\text{sign}(p_{i})b_{i}(x)\Delta\partial_{p_{i}}U=0,$
(8)

in which the coefficients of the equation are independent of $z$!

To make sure the solution to new Eq. (8) is mathematically well-defined, we
first assume $b_{i}$ is independent of $x$ for simplicity and take a Fourier
Transform on $x$

$\partial_{t}\hat{U}-|\xi|^{2}\sum_{i=1}^{L}{\text{sign}}(p_{i})b_{i}\partial_{p_{i}}{\hat{U}}=0,$
(9)

where the Fourier variable in $x$ is denoted by $\xi\in\mathbb{R}^{d}$. This
is a linear transport equation for $\hat{U}$ with a discontinuous coefficient.
Often one needs to make sense of the solution at $p=0$, depending on the
physical background of the problem. One possibility is to provide a physically
relevant jump condition at $p=0$ [31]. However, from the definition in Eq.
(6), $U$ is continuous at $p_{i}=0$, so there is no need to impose any
interface condition for $U$ at $p_{i}=0$. This can also be justified from the
solution to Eq. (9). For clarity, consider the case of $L=1$ and $b=1$. Then
Eq. (9) is

$\partial_{t}\hat{U}-\xi^{2}{\text{sign}}(p)\partial_{p}{\hat{U}}=0.$ (10)

By the method of characteristics,

$\hat{U}(t,\xi,z,0^{-})=\hat{U}(0,\xi,z,-\xi^{2}t),\qquad\hat{U}(t,\xi,z,0^{+})=\hat{U}(0,\xi,z,\xi^{2}t).$
(11)

From the definition in Eq. (6) one clearly sees that
$\hat{U}(0,\xi,z,-\xi^{2}t)=\hat{U}(0,\xi,z,\xi^{2}t)$, therefore

$\hat{U}(t,\xi,z,0^{-})=\hat{U}(t,\xi,z,0^{+})$ (12)

and $\hat{U}$ is continuous at $p=0$.

By method of characteristics, it is also easy to check that $\hat{U}$ is an
even function for each $p_{i}$, since its initial data is also an even
function.

To work with $M$ samples $\\{z_{m}\\}$, $m=1,...,M$, we now define

$V(t,x,p)=\frac{1}{M}\sum_{m=1}^{M}U(t,x,z_{m},p)$

for $p_{i}\in(-\infty,\infty),\,i=1,\cdots,L$, which solves

$\displaystyle\begin{cases}\partial_{t}V+\sum_{i=1}^{L}{\text{sign}}(p_{i})b_{i}(x)\Delta\partial_{p_{i}}V=0\,\\\
V(0,x,p)=\frac{1}{M}\sum_{m=1}^{M}\prod_{i=1}^{L}a_{i}(z_{m})e^{-a_{i}(z_{m})|p_{i}|}u(0,x,z_{m})\,.\end{cases}$
(13)

This is the linear PDE we will solve, which has certain coefficients and a
single initial condition. Now the average of the solutions of the original
problem, Eq. (5), can be recovered from $V(t,x,p)$ using

$\displaystyle\overline{u}(t,x)$
$\displaystyle=\frac{1}{M}\sum_{m=1}^{M}u(t,x,z_{m})=\int V(t,x,p)\,dp\,.$
(14)

Thus, in solving for $V$, the computational cost is clearly independent of
$M$!

To solve Eq. (13) numerically to, say $t=1$, we use $N_{t}$ steps in time:
$t=n\Delta t$ for $k=0,...,N_{t}-1$ where $N_{t}\Delta t=1$. Throughout the
paper we assume the spatial computational domain to be $[0,1]^{d}$, with $N$
spatial mesh points for each dimension. For example, in one space dimension
$x=j\Delta x$ for $j=0,...,N-1$ where $N\Delta x=1$. Since $V$ decays
exponentially in $|p|$, one can truncate the computational domain in finite
$p$, say at $p_{i}=\pm 1$, namely $p_{i}\in[-1,1]$ for all $1\leq i\leq L$.
Define the discrete $p$ as $p=-1+k\Delta p$ for $k=0,...,N_{p}-1$, where
$N_{p}\Delta p=2$. In $d$ spatial dimensions, the total number of spatial
meshpoints for $x$ is $N^{d}$ and for $p$ is $N_{p}^{L}$ since $p$ is an
$L$-dimensional vector.

We use the center finite difference in space, the upwind scheme in $p$ (so the
overall spatial error is of $O(1/N)$) and the forward Euler scheme in time (so
the time error is of $O(\Delta t)$) to solve the phase space heat equation in
Eq. (13). See Appendix A.1 for details of the discretisation scheme and the
corresponding properties of matrix $\mathcal{M}$.

We use the quadrature rule for the integration in using the discrete values
$V_{N_{t},k}=V(t=N_{t}\Delta t,p=k\Delta p)$. If the numerical solution vector
is denoted
$\boldsymbol{V}=\sum_{n}^{N_{t}}\sum_{j}^{N^{d}}\sum_{k}^{N_{p}^{L}}V_{n,j,k}|n\rangle|j\rangle|k\rangle$,
then the ensemble average of $u$ at point $(t=N_{t}\Delta t,x=J\Delta x)$ is

$\displaystyle\bar{u}(t,x)=\int
V(t,x,p)dp\approx\bar{u}_{N_{t},J}\equiv\frac{1}{N_{p}^{L}}\sum_{k}^{N_{p}^{L}}V_{N_{t},J,k}$
(15)

and
$\boldsymbol{V}^{T}\mathcal{G}\boldsymbol{V}=N_{p}^{L}|\bar{u}_{N_{t},J}|^{2}$,
where $\mathcal{G}=|G_{N_{t},J}\rangle\langle G_{N_{t},J}|$ with
$|G_{N_{t},J}\rangle=(1/\sqrt{N_{p}^{L}})\sum_{k}^{N_{p}^{L}}|N_{t}\rangle|J\rangle|k\rangle$.
The quantum state embedding of the initial condition of the state is defined
as

$\displaystyle|V_{0}\rangle=\frac{1}{\mathcal{N}_{V_{0}}}\sum_{j}^{N^{d}}\sum_{k}^{N_{p}^{L}}V_{0,j,k}|0\rangle|j\rangle|k\rangle,$
(16)

where
$\mathcal{N}^{2}_{V_{0}}=\sum_{j}^{N^{d}}\sum_{k}^{N_{p}^{L}}|V_{0,j,k}|^{2}$
and the state is assumed to have sparsity $\sigma_{V_{0}}$.

######  Lemma 7.

The constant $n^{2}_{V_{0}}\equiv\mathcal{N}^{2}_{V_{0}}/N_{p}^{L}$ has range
$O(1)\leq n^{2}_{V_{0}}\leq O(N^{d})$. Different $n_{\psi_{0}}$ corresponds to
different initial data. If we assume the initial data has support in a box of
size $\beta\in(0,1)$, then $n_{V_{0}}^{2}=O((\beta N)^{d})$.

###### Proof.

We first prove that $\int dp\,|V(0,x,p)|^{2}\leq C=O(1)$. First, it is easy to
check that $\int
dp\,|U(0,x,p,z)|^{2}=\frac{1}{4}\prod_{i=1}^{L}a_{i}(z)^{2}|u_{0}(x,z)|^{2}\leq
C=O(1)$. Here the boundedness of $\prod_{i=1}^{L}a_{i}(z)^{2}$ is a result of
the boundedness of $\sum_{i=1}^{L}a_{i}(z)^{2}$ (from Eq. (2)). Then

$\int dp\,|V(0,x,p)|^{2}=\int
dp\left(\frac{1}{M}\sum_{m=1}^{M}U(0,x,p,z_{m})\right)^{2}\leq\frac{1}{M}\sum_{m=1}^{M}\int
dp\,U(0,x,p,z_{m})^{2}\leq C=O(1).$

Then using the quadrature rule we have $O(1)\geq\int dx\int
dp\,|V(0,x,p)|^{2}\approx(1/(\beta N)^{d})\sum_{j}^{N^{d}}\int
dp\,|V(0,x=j\Delta x,p)|^{2}$. Then $O(1)\geq\int dx\int
dp\,|V(0,x,p)|^{2}\approx\mathcal{N}^{2}_{V_{0}}/((\beta N)^{d}N_{p}^{L})$.
Thus $O(1)\leq n^{2}_{V_{0}}\equiv\mathcal{N}_{V_{0}}^{2}/N_{p}^{L}\leq
O((\beta N)^{d})$. ∎

Then the quantum algorithm needed to recover $\bar{u}$ has the following
complexity.

######  Theorem 8.

A quantum algorithm that takes sparse access to $\mathcal{M}$ (using an
$r^{\text{th}}$-order method with $r\geq 1$) and access to classical values of
the initial conditions $V_{0,J,l}$, is able to estimate the density
$|\bar{u}_{N_{t},J}|^{2}$ to precision $\epsilon$ at $\Lambda$ meshpoints,
with an upper bound on the gate complexity

$\displaystyle\mathcal{Q}=\tilde{\mathcal{O}}(n^{2}_{V_{0}}\Lambda(Ld)^{3}(Ld/\epsilon)^{1+9/r})$
(17)

and a smaller query complexity.

###### Proof.

The sparsity and condition number of the corresponding matrix $\mathcal{M}$
for Eq. (13) is $s=O(Ld)$ and $\kappa=O(LdN^{3})=O(Ld(Ld/\epsilon)^{3/r})$
respectively, where $N=O((Ld/\epsilon)^{1/r})$ from Appendix B.1. From Lemma 3
the cost in recovering
$\boldsymbol{V}^{T}\mathcal{G}\boldsymbol{V}=N_{p}^{L}\bar{u}_{N_{t},J}^{2}$
to precision $\epsilon^{\prime}$ is
$\tilde{\mathcal{O}}(s\kappa^{3}\mathcal{N}_{V_{0}}^{2}/\epsilon^{\prime})=\tilde{\mathcal{O}}(s\kappa^{3}n_{V_{0}}^{2}/\epsilon)=\tilde{\mathcal{O}}(n_{V_{0}}^{2}(Ld)^{3}(Ld/\epsilon)^{1+9/r})$
where $\epsilon^{\prime}=N_{p}^{L}\epsilon$, $\epsilon$ is the error in
$\bar{u}_{N_{t},J}^{2}$. Since $\Lambda$ different states
$|G_{N_{t},J}\rangle$ are required for each mesh point $(N_{t},J)$, so the
total cost must be multiplied by $\Lambda$. ∎

######  Lemma 9.

When $M<M_{heat}\equiv O(L^{2+(d+L+3)/r}(d/\epsilon)^{(L+1)/3})$, the
classical algorithm has minimal cost
$\mathcal{C}=O(Md^{2}(d/\epsilon)^{(d+2)/r})$. If $M>M_{heat}$, the classical
algorithm has minimal cost
$\mathcal{C}=O(L^{2+(d+L+3)/r}d^{2}(d/\epsilon)^{(d+L+3)/r})$.

###### Proof.

See Appendix B.1 for details. ∎

Now we define what we mean by quantum advantage with respect to classical
algorithms with cost $\mathcal{C}$.

######  Definition 10.

We say there is a quantum advantage in estimating the quantity of interest
when $\mathcal{Q}={\scriptstyle\mathcal{O}}(\mathcal{C})$.

######  Corollary 11.

To attain a quantum advantage when $M<M_{heat}$, it is sufficient for the
following condition to hold

$\displaystyle{\scriptstyle\mathcal{O}}\left(\frac{M}{n^{2}_{V_{0}}\Lambda
L^{4+9/r}d}\left(\frac{d}{\epsilon}\right)^{(d-7)/r}\right)=\tilde{O}(1).$
(18)

When $M>M_{heat}$, it is sufficient that the following is satisfied

$\displaystyle{\scriptstyle\mathcal{O}}\left(\frac{1}{n^{2}_{V_{0}}\Lambda
Ld}\left(\frac{Ld}{\epsilon}\right)^{(d+L-6)/r-1}\right)=\tilde{O}(1).$ (19)

###### Proof.

Combining Lemma 9 and Theorem 8. ∎

We see from Eq. (18) that quantum advantage with respect to $M$ is always
possible in the range $M^{\prime}_{heat}<M<M_{heat}$ where
$M^{\prime}_{heat}=\tilde{\mathcal{O}}(n^{2}_{V_{0}}\Lambda
L^{4+9/r}d(\epsilon/d)^{(d-7)/r})$, with up to exponential quantum advantage
in $d$ and $\epsilon$. As a simple example, if we begin with a point source
initial condition and only require the final solution at $O(1)$ points, then
$n^{2}_{V_{0}}=O(1)=\Lambda$. When $d>7$, the range of possible $M$, captured
by $M_{heat}-M^{\prime}_{heat}$, is very large. A necessary condition for such
an $M$ to exist in more general cases is
$O(1)<n^{2}_{V_{0}}\Lambda<\tilde{\mathcal{O}}(L^{(d+L-6)/r-2}(d/\epsilon)^{(L+1)/3-(d-7)/r}/d)$,
which allows for a wide range of possibilities.

In the case when $M>M_{heat}$, we see there is no quantum advantage with
respect to $M$, but there is quantum advantage in $L,d,\epsilon$ when
$n^{2}_{V_{0}}\Lambda<\tilde{\mathcal{O}}((Ld/\epsilon)^{(d+L-6)/r-1}/(Ld))$.

###  III.2 The linear Boltzmann equation

Now consider the linear Boltzmann equation with isotropic scattering [32]

$\begin{cases}\partial_{t}f+v\cdot\nabla_{x}f=a(x,z)\left[\frac{1}{\Omega}\int_{\mathbb{R}^{d}}f\,dv-f\right]\,,\\\
f(0,x,z)=f_{0}(x,z)\,,\end{cases}$ (20)

with suitable boundary conditions. Here $f(t,x,v,z)>0$ is the particle density
distribution, $x\in\mathbb{R}^{d}$ is the position, $v\in\mathbb{R}^{d}$ is
the velocity of the particles, $a(x,z)>0$, defined in (1), is the scattering
rate and $\Omega$ is the volume of the domain for $v$. Here $a(x,z)$ is often
uncertain due to experimental or modeling errors [33]. The left hand side of
the equation models particle transport, while the right hand side models
scattering collision of particles with the background.

We introduce the following transformation

$F^{\prime}(t,x,v,z,p)=\frac{1}{2}\prod_{i=1}^{L}a_{i}(z)e^{-|p_{i}|a_{i}(z)}f(t,x,v,z)$
(21)

where $p_{i}\in(-\infty,\infty),i=1,\cdots,L$. Then $f$ can be recovered from
$F$ via

$f(t,x,v,z)=\int_{(-\infty,\infty)^{L}}F^{\prime}(t,x,v,z,p)=2\int_{[0,\infty)^{L}}F^{\prime}(t,x,v,z,p)\,dp,$

and $F^{\prime}$ solves

$\partial_{t}F^{\prime}+v\cdot\nabla_{x}F^{\prime}=-\sum_{i=1}^{L}{\text{sign}}(p_{i})\left[\frac{1}{\Omega}\int_{\mathbb{R}^{d}}b_{i}(x)\partial_{p_{i}}F^{\prime}\,dv-
b_{i}(x)\partial_{p_{i}}F^{\prime}\right]\,.$ (22)

Then the following stability and conservation results hold.

######  Theorem 12.

Assume vanishing boundary condition for $F^{\prime}$ and
$\sum_{i=1}^{L}b_{i}(x)\geq 0$. Then

$\displaystyle\partial_{t}\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}F^{\prime}\,dv\,dx=0$
(23)
$\displaystyle\partial_{t}\int_{\mathbb{R}^{d}}\int_{\mathbb{R}^{d}}\int_{[0,\infty)^{L}}(F^{\prime})^{2}\,dp\,dv\,dx\leq
0\,.$ (24)

###### Proof.

Eq. (23) is easy to see by integrating over $v$ and $x$ in (22) and using the
vanishing boundary condition. Since $F^{\prime}$ is an even function in
$p_{i}$, we will prove (24) for the case of $p_{i}>0$ for all $i$. Multiplying
(22) by 2F’ gives

$\partial_{t}(F^{\prime})^{2}+v\cdot\nabla_{x}(F^{\prime})^{2}=-\sum_{i=1}^{L}b_{i}(x)\left[\frac{2}{\Omega}F^{\prime}\int_{\mathbb{R}^{d}}\partial_{p_{i}}F^{\prime}\,dv-\partial_{p_{i}}(F^{\prime})^{2}\right]\,.$
(25)

Integrating over $v$, one gets

$\partial_{t}\int_{\mathbb{R}^{d}}(F^{\prime})^{2}dv+\nabla_{x}\cdot\int_{\mathbb{R}^{d}}v(F^{\prime})^{2}dv=-\sum_{i=1}^{L}b_{i}(x)\left[\frac{1}{\Omega}\partial_{p_{i}}\left(\int_{\mathbb{R}^{d}}F^{\prime}dv\right)^{2}-\partial_{p_{i}}\int_{\mathbb{R}^{d}}(F^{\prime})^{2}\,dv\right]\,.$
(26)

Now integrating over $p$. Denoting
$p_{i}^{-}=(p_{1},\cdots,p_{i-1},p_{i+1},\cdots,p_{L})$ and using
$F^{\prime}(t,x,v,p)|_{p_{i}=0}$, gives

$\displaystyle\partial_{t}\int_{[0,\infty)^{L}}\int_{\mathbb{R}^{d}}(F^{\prime})^{2}\,dv\,dp+\nabla_{x}\cdot\int_{[0,\infty)^{L}}\int_{\mathbb{R}^{d}}v(F^{\prime})^{2}\,dv\,dp$
$\displaystyle=$
$\displaystyle\sum_{i=1}^{L}b_{i}(x)\int_{[0,\infty)^{L-1}}\left[\frac{1}{\Omega}\left(\int_{\mathbb{R}^{d}}F^{\prime}(t,x,v,z,p)|_{p_{i}=0}\,dv\right)^{2}-\int_{\mathbb{R}^{d}}F^{\prime}(t,x,v,z,p)^{2}|_{p_{i}=0}\,dv\right]dp_{i}^{-}\leq
0\,,$ (27)

where the inequality is obtained by Jensen’s inequality. Now integrating over
$x$, (24) is proved. ∎

Suppose one is interested in computing the ensemble average of $M$ different
initial data resulting from $M$ samples of $z$. Then similarly to the heat
equation case, we can define

$\displaystyle
F(t,x,v,p)=\frac{1}{M}\sum_{m=1}^{M}F^{\prime}(t,x,v,z_{m},p)=\frac{1}{M}\sum_{m=1}^{M}\prod_{i=1}^{L}a_{i}(z_{m})e^{-|p_{i}|a_{i}(z_{m})}f(t,x,v,z_{m})$
(28)

to obtain the following problem:

$\begin{cases}\partial_{t}F+v\cdot\nabla_{x}F=-\sum_{i=1}^{L}{\text{sign}}(p_{i})\left[\frac{1}{\Omega}\int_{\mathbb{R}^{d}}b_{i}(x)\partial_{p_{i}}F\,dv-
b_{i}(x)\partial_{p_{i}}F\right]\,\\\
F(0,x,v,p)=\frac{1}{M}\sum_{m=1}^{M}\Pi_{i=1}^{L}a_{i}(z_{m})e^{-|p_{i}|a_{i}(z_{m})}f_{0}(x,v,z_{m})\,.\end{cases}$
(29)

Then the average of the solutions of the original problem, Eq. (20), can be
recovered from $F(t,x,p)$ using

$\displaystyle\overline{f}(t,x,v)$
$\displaystyle=\frac{1}{M}\sum_{m=1}^{M}f(t,x,v,z_{m})=\int F(t,x,v,p)\,dp\,.$
(30)

In practical applications one is interested in the moments of $f$, which give
rises to physical quantities of interest including, for example, density
$\rho(t,x)=\int f(t,x,v)\,dv=\int\int F(t,x,v,p)\,dp\,dv$, flux $J(t,x)=\int
vf(t,x,v)\,dv=\int\int vF(t,x,v,p)\,dp\,dv$, and kinetic energy
$E(t,x)=\int\frac{v^{2}}{2}f(t,x,v)\,dv=\int\int\frac{v^{2}}{2}F(t,x,v,p)\,dp\,dv$.
The ensemble averages of these physical quantities can then be computed using
Eq. (31).

Again, notice that $F$ decays exponentially in $p$, which means we can
truncate the computational domain in some finite point, without loss of
generality, for instance at at $p_{i}=\pm 1$. We consider $N_{t}$ steps in
time $t_{n}=n\Delta t$ for $n=0,...,N_{t}-1$ where $N_{t}\Delta t=1$. We also
use $N$ mesh points for each space dimension. For example, in one dimension
$x_{j}=j\Delta x$ for $j=0,...,N-1$ where $N\Delta x=1$, and $p_{k}=-1+k\Delta
p$ for $k=0,...,N_{p}-1$, where $N_{p}\Delta p=2$. To simplify the notation we
let $N=N_{p}$. The discrete-ordinate method (namely, quadrature rules for the
collision term) is used [32]) with quadrature points $v_{l}=l\Delta v$, for
$n=0,...,N_{v}-1$ where $N_{v}\Delta v=1$. In $d$ spatial dimensions, the
total number of spatial mesh points for $x$ is $N^{d}$, for $v$ is $N_{v}^{d}$
and for $p$ is $N_{p}^{L}$ since $p$ is $L$-dimensional. To solve Eq. (29), we
use the upwind scheme for $x$ and $w$, and forward Euler method in time. For
details of the numerical scheme and the corresponding matrix $\mathcal{M}$ for
the problem, see Appendix A.2. The overall spatial error is of $O(1/N)$ and
temporal error is of $O(\Delta t)$.

Let $F_{n,j,l,k}=F(t=n\Delta t,x=j\Delta x,v=l\Delta v,p=k\Delta p)$,
$\rho_{n,j}=\rho(t=n\Delta t,x=j\Delta x)$ and $J_{n,j}=J(t=n\Delta
t,x=j\Delta x)$, $v_{l}=l\Delta v$, where $n=1,...,N_{t}$, $j=1,...,N^{d}$,
$l=1,...,N_{v}^{d}$, $k=1,...,N_{p}^{L}$. Then one can numerically approximate
the density, flux and energy through the following quadrature rule:

$\displaystyle\rho(t,x)=\int\int
F(t,x,v,p)dvdp\approx\rho_{n,j}\equiv\frac{1}{N_{v}^{d}N_{p}^{L}}\sum_{l}^{N_{v}^{d}}\sum_{k}^{N_{p}^{L}}F_{n,j,l,k},$
(31) $\displaystyle J(t,x)=\int\int vF(t,x,v,p)dvdp\approx
J_{n,j}\equiv\frac{1}{N_{v}^{d}N_{p}^{L}}\sum_{l}^{N^{d}_{v}}\sum_{k}^{N_{p}^{L}}v_{l}F_{n,j,l,k},$
(32) $\displaystyle E(t,x)=\int\int(v^{2}/2)F(t,x,v,p)dvdp\approx
E_{n,j}\equiv\frac{1}{2N_{v}^{d}N_{p}^{L}}\sum_{l}^{N_{v}^{d}}\sum_{k}^{N_{p}^{L}}v^{2}_{l}F_{n,j,l,k}.$
(33)

Denoting the solution vector to Eq. (29) by
$\boldsymbol{F}=\sum_{n}^{N_{t}}\sum_{j}^{N^{d}}\sum_{l}^{N_{v}^{d}}\sum_{k}^{N_{p}^{L}}F_{n,j,l,k}|n\rangle|j\rangle|l\rangle|k\rangle$,
one can write

$\displaystyle\boldsymbol{F}^{T}\mathcal{G}^{density}\boldsymbol{F}=N_{v}^{d}N_{p}^{L}|\rho_{n,j}|^{2},$
(34)
$\displaystyle|G^{density}_{n,j}\rangle=\frac{1}{\sqrt{N_{v}^{d}N_{p}^{L}}}\sum_{l}^{N_{v}^{d}}\sum_{k}^{N_{p}^{L}}|n\rangle|j\rangle|l\rangle|k\rangle,$
$\displaystyle\boldsymbol{F}^{T}\mathcal{G}^{flux}\boldsymbol{F}=N^{L}_{p}\left(\frac{N^{d}_{v}}{\mathcal{N}_{V}}\right)^{2}|J_{n,j}|^{2}=O(N_{v}^{d}N_{p}^{L}|J_{n,j}|^{2}),$
(35)
$\displaystyle|G^{flux}_{n,j}\rangle=\frac{1}{\mathcal{N}_{V}\sqrt{N_{p}^{L}}}\sum_{l}^{N_{v}^{d}}\sum_{k}^{N_{p}^{L}}v_{l}|n\rangle|j\rangle|l\rangle|k\rangle,$
$\displaystyle\boldsymbol{F}^{T}\mathcal{G}^{energy}\boldsymbol{F}=4N^{L}_{p}\left(\frac{N^{d}_{v}}{\mathcal{N}_{V^{2}}}\right)^{2}|E_{n,j}|^{2}=O(N_{v}^{d}N_{p}^{L}|E_{n,j}|^{2}),$
(36)
$\displaystyle|G^{energy}_{n,j}\rangle=\frac{1}{\mathcal{N}_{V^{2}}\sqrt{N_{p}^{L}}}\sum_{l}^{N_{v}^{d}}\sum_{k}^{N_{p}^{L}}v^{2}_{l}|n\rangle|j\rangle|l\rangle|k\rangle,$

where
$\mathcal{N}_{V}^{2}=\sum_{l}^{N_{v}^{d}}|v_{l}|^{2}=\sum_{l}^{N_{v}^{d}}(l\Delta
v)^{2}=O(N_{v}^{d})$ and
$\mathcal{N}^{2}_{V^{2}}=\sum_{l}^{N_{v}^{d}}|v^{2}_{l}|^{2}=\sum_{l}^{N_{v}^{d}}(l\Delta
v)^{4}=O(N_{v}^{d})$. The initial condition can be encoded in $|F_{0}\rangle$
defined by

$\displaystyle|F_{0}\rangle=\frac{1}{\mathcal{N}_{F_{0}}}\sum_{j}^{N^{d}}\sum_{l}^{N_{v}^{d}}\sum_{k}^{N_{p}^{L}}F_{0,j,l,k}|0\rangle|j\rangle|l\rangle|k\rangle,$
(37)

where the normalisation
$\mathcal{N}^{2}_{F_{0}}=\sum_{j}^{N^{d}}\sum_{l}^{N_{v}^{d}}\sum_{k}^{N_{p}^{L}}|F_{0,j,l,k}|^{2}$
and the initial state is assumed to have sparsity $\sigma_{F_{0}}$. The
kinetic energy can be computed similarly and we omit its details.

######  Lemma 13.

The constant $n^{2}_{F_{0}}\equiv\mathcal{N}^{2}_{F_{0}}/(N_{p}^{L}N_{v}^{d})$
has range $O(1)\leq n^{2}_{F_{0}}\leq O(N^{d})$. 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^{2}_{F_{0}}=O((\beta N)^{d})$.

###### Proof.

First, similar to the proof in Lemma 7, $\int dv\int dp|F(0,x,v,p)|^{2}\leq
O(1)$. Then using the quadrature rule we have $\int dx\int dv\int
dp|F(0,x,v,p)|^{2}\approx(1/(\beta N)^{d})\sum_{j}^{N^{d}}\int dv\int
dp|F(0,x=j\Delta x,v,p)|^{2}\leq O(1)$. Then $O(1)\geq\int dx\int dv\int
dp|F(0,x,v,p)|^{2}\approx\mathcal{N}^{2}_{F_{0}}/((\beta
N)^{d}N_{v}^{d}N_{p}^{L})$. Thus
$n^{2}_{F_{0}}\equiv\mathcal{N}_{F_{0}}^{2}/(N_{p}^{L}N_{v}^{d})\leq O((\beta
N)^{d})$. However, the quantum inner product $\Upsilon\equiv\langle
F_{0}|(\mathcal{M}^{-1})^{\dagger}\mathcal{G}\mathcal{M}^{-1}|F_{0}\rangle$ is
upper bounded by $1$, where
$1\geq\Upsilon=\boldsymbol{F}^{T}\mathcal{G}\boldsymbol{F}/\mathcal{N}^{2}_{F_{0}}=O(N_{v}^{d}N_{p}^{L}/\mathcal{N}^{2}_{F_{0}})=O(1/n^{2}_{F_{0}})$,
since the observables $|\rho|^{2},|J|^{2}=O(1)$ for
$\mathcal{G}=\mathcal{G}^{density},\mathcal{G}^{flux}$. ∎

Remark: The definition
$n_{F_{0}}\equiv\mathcal{N}_{F_{0}}/(N_{p}^{L}N_{v}^{d})^{1/2}$ can be an
underestimate in situations where the $N_{p}^{L}N_{v}^{d}$ term in the
quadrature sum overestimates the number of non-zero factors in the summation.

With this, we can compute the total cost for an ‘end-to-end’ quantum algorithm
for computing both the density and flux.

######  Theorem 14.

A quantum algorithm that takes sparse access to $\mathcal{M}$ (using an
$r^{\text{th}}$-order method with $r\geq 1$) is able to estimate
$|\rho_{N_{t},J}|^{2}$, $|J_{N_{t},J}|^{2}$ and $|E_{N_{t},J}|^{2}$ at
$\Lambda$ mesh points to precision $\epsilon$ with an upper bound on the gate
complexity

$\displaystyle\mathcal{Q}=\tilde{\mathcal{O}}\left(n^{2}_{F_{0}}\Lambda(L/\epsilon)(d+L)^{3}\left(\frac{Ld}{\epsilon}\right)^{(3+d)/r}\right)$
(38)

with a smaller query complexity.

###### Proof.

The matrix $\mathcal{M}$ corresponding to the discretisation of Eq. (29) has
sparsity $s=O(LN^{d}_{v})=O(L(Ld/\epsilon)^{d/r})$ and condition number
$\kappa=O((d+L)N)=O((d+L)(Ld/\epsilon)^{1/c})$, where
$N=N_{p}=N_{v}=O((Ld/\epsilon)^{1/r})$ from Appendix B.2. To compute
$\boldsymbol{F}^{T}\mathcal{G}^{density}\boldsymbol{F}$,
$\boldsymbol{F}^{T}\mathcal{G}^{flux}\boldsymbol{F}$ and
$\boldsymbol{F}^{T}\mathcal{G}^{energy}\boldsymbol{F}$ to precision
$\epsilon^{\prime}$, from Lemma 3 the maximum query complexity is
$\tilde{\mathcal{O}}(s\kappa^{3}\mathcal{N}^{2}_{F_{0}}/\epsilon^{\prime})=\tilde{\mathcal{O}}(s\kappa^{3}n^{2}_{F_{0}}/\epsilon)=\tilde{\mathcal{O}}(n^{2}_{F_{0}}(d+L)^{3}(L/\epsilon)(Ld/\epsilon)^{(3+d)/r})$
where $\epsilon=\epsilon^{\prime}/(N_{p}^{L}N_{v}^{d})$ is the error in
$|\rho_{N_{t},J}|^{2}$, $|J_{N_{t},J}|^{2}$ and $|E_{N_{t},J}|^{2}$. The same
order of additional 2-qubit gates are required. If one wants the density at
$\Lambda$ meshpoints, then total gate complexity must be multiplied by
$\Lambda$. The same argument applies to computing the flux. ∎

######  Lemma 15.

When $M<M_{Boltz}\equiv O(L^{(2d+L+1)/r}\max(L,d)(d/\epsilon)^{L/r}/d)$, the
classical algorithm for the problem has the cost
$\mathcal{C}=O(Md(d/\epsilon)^{(2d+1)/r})$. When $M>M_{Boltz}$, the classical
algorithm has cost $\mathcal{C}=O((Ld/\epsilon)^{(2d+L+1)/r}\max(L,d))$.

###### Proof.

See Appendix B.2 for details. ∎

######  Corollary 16.

To attain a quantum advantage when $M<M_{Boltz}$, it is sufficient for the
following condition to hold

$\displaystyle{\scriptstyle\mathcal{O}}\left(\frac{Md^{1+(d-2)/r}}{n^{2}_{F_{0}}\Lambda
L^{1+(3+d)/r}(d+L)^{3}}\left(\frac{1}{\epsilon}\right)^{(d-2)/r-1}\right)=\tilde{O}(1).$
(39)

When $M>M_{Boltz}$, it is sufficient that the following is satisfied

$\displaystyle{\scriptstyle\mathcal{O}}\left(\frac{\max(L,d)}{n^{2}_{F_{0}}\Lambda(L/\epsilon)(d+L)^{3}}\left(\frac{Ld}{\epsilon}\right)^{(d+L-2)/r}\right)=\tilde{O}(1).$
(40)

One can see from Eq. (39) that quantum advantage with respect to $M$ is always
possible in the range $M^{\prime}_{Boltz}<M<M_{Boltz}$ where
$M^{\prime}_{Boltz}=\tilde{\mathcal{O}}(n^{2}_{F_{0}}\Lambda
L^{1+(3+d)/r}(d+L)^{3}\epsilon^{(d-2)/r-1}/d^{(d-2)/r+1})$, with up to
exponential quantum advantage in $d$ and $\epsilon$. As a simple example, if
we begin with a point source initial condition and only require the final
solution at $O(1)$ points, then $n^{2}_{F_{0}}=O(1)=\Lambda$. When $d>2+r$,
the range of possible $M$, captured by $M_{Boltz}-M^{\prime}_{Boltz}$, is very
large. A necessary condition for such an $M$ to exist in more general cases is
$O(1)<n^{2}_{F_{0}}\Lambda<\tilde{\mathcal{O}}(L^{(d+L-2)/r-1}d(d/\epsilon)^{(d+L-2)/r-1}\max(L,d)/(d+L)^{3})$,
which allows for a wide range of possibilities.

In the case when $M>M_{Boltz}$, we see from Eq. (40) there is no quantum
advantage with respect to $M$, but there is quantum advantage in
$L,d,\epsilon$ when
$n^{2}_{V_{0}}\Lambda<\tilde{\mathcal{O}}((Ld/\epsilon)^{(d+L-2)/r}\max(L,d)/((L/\epsilon)(d+L)^{3}))$.
Depending on $n^{2}_{V_{0}}\Lambda$, up to exponential quantum advantage in
$L,d,\epsilon$ is possible.

###  III.3 Linear advection equation

Now consider the linear advection equation

$\begin{cases}\partial_{t}u\pm a(x,z)\sum_{j=1}^{d}\partial_{x_{j}}u=0,\\\
u(0,x,z)=u_{0}(x,z).\end{cases}$ (41)

Here $a(x,z)$ in uncertain due to the heterogeneity of the media [34]. For the
case of $L=1$, applying the transformation $U(t,x,z,p)=a(z)e^{-pa(z)}u(t,x,z)$
as before, one can easily show that $U$ solves

$\begin{cases}\partial_{t}U\mp\sum_{j=1}^{d}\partial_{x_{j}p}\,U=0,\\\
U(0,x,z,p)=a(z)e^{-pa(z)}u(0,x,z),\end{cases}$ (42)

where $U$ satisfies a second order equation in the phase space. However, this
problem is ill-posed, which can be easily seen by applying a Fourier transform
on $x$. Denote the Fourier variable by $\xi\in\mathbb{R}^{d}$ in $x$ and
$\iota=\sqrt{i}$. Then in the Fourier space the equation becomes

$\partial_{t}{\hat{U}}\pm\iota\left(\sum_{j=1}^{d}\xi_{j}\right)\partial_{p}{\hat{U}}=0$
(43)

This is a convection equation with imaginary wave speed, which is ill-posed.

Instead we introduce a new transformation

$U(t,x,z,p)=\prod_{i=1}^{L}\sqrt{a_{i}(z)}\,e^{-\sqrt{a_{i}(z)}p_{i}}u(t,x,z)$
(44)

and one arrives at the following problem:

$\begin{cases}\partial_{t}U\pm\sum_{i=1}^{L}b_{i}(x)\sum_{j=1}^{d}\partial_{x_{j}p_{i}p_{i}}\,U=0,\\\
U(0,x,p)=\prod_{i=1}^{L}\sqrt{a_{i}(z)}e^{-p_{i}\sqrt{a_{i}(z)}}u(0,x,z),\\\
\partial_{p_{i}}U(t,x,z,p)+\sqrt{a_{m}(z)}\,U(t,x,z,p)\Big{|}_{p_{i}=0}=0,\qquad{\text{for
all}}\,i\,,\end{cases}$ (45)

where is a third order equation in the phase space.

The well-posedness of the initial value problem can be seen by applying a
Fourier transform on $x$, for the case of constant $b_{i}(x)$. Assume the
Fourier variable is $\xi=(\xi_{1},\cdots,\xi_{d})^{T}\in\mathbb{R}^{d}$. Then
taking a Fourier transform for $U$ on $x$ gives

$\partial_{t}{\hat{U}}\mp\iota\sum_{i=1}^{L}b_{i}\left(\sum_{j=1}^{d}\xi_{j}\right)\,\partial_{p_{i}p_{i}}{\hat{U}}=0,$
(46)

This is the free Schrödinger equation which is a good equation (well-posed)!

However, since the Robin boundary condition in Eq. (45) still depends on
$a(z)$, one will not be able to linearly superimpose solutions for different
sample $z_{m}$ to obtain the ensemble average. To remedy this problem, we
introduce a new transformation

${W}(t,x,p,q)=\frac{1}{M}\sum_{m=1}^{M}\prod_{i=1}^{L}\sqrt{a_{i}(z_{m})}\,e^{-\sqrt{a_{i}(z_{m})}(p_{i}+q_{i})}u(t,x,z_{m})$
(47)

where $q\in\mathbb{R}^{L}$, which then solves

$\begin{cases}\partial_{t}{W}\pm\sum_{i=1}^{L}b_{i}(x)\sum_{j=1}^{d}\partial_{x_{j}p_{i}p_{i}}\,{W}=0,\\\
{W}(0,x,p,q)=\frac{1}{M}\sum_{m=1}^{M}\prod_{i=1}^{L}\sqrt{a_{i}(z_{m})}e^{-(p_{i}+q_{i})\sqrt{a_{i}(z_{m})}}u(0,x,z_{m}),\\\
\partial_{p_{i}}{W}-\partial_{q_{i}}{W}\Big{|}_{p_{i}=0}=0\quad{\text{for
all}}\,\,i\,.\end{cases}$ (48)

Now the coefficients in both the equation and the boundary condition no longer
depend on $a(z)$!

The average of the solutions of the original problem, Eq. (41), can be
recovered from $W(t,x,p,q)$ using

$\displaystyle\overline{u}(t,x)$
$\displaystyle=\frac{1}{M}\sum_{m=1}^{M}u(t,x,z_{m})=\int\int{W}(t,x,p,q)\,dp\,dq\,.$
(49)

Due to the exponential decay in $p$ and $q$ we also use the homogeneous
boundary condition at the right hand side of the domain for $p$ and $q$ in a
suitably truncated computational domain.

Remark: Unlike for the linear heat equation and Boltzmann equations, the
transformation where we extend the domain of $p_{i}$ to $(-\infty,\infty)$
does not work for the advection equation. For instance, if one uses
$U(t,x,z,p)=(1/2)\sqrt{a(z)}\,e^{-\sqrt{a(z)}|p|}u(t,x,z)$, then one
transforms the advection equation $\partial_{t}u+a(z)\partial_{x}u=0$ into
$\partial_{t}U+\partial_{xpp}U=-\sqrt{a(z)}\delta(p)\partial_{x}U$, which
clearly still has $a(z)$ dependence in its coefficient.

Remark: To compute the variance, we can multiply (41) by $u$ to derive

$\begin{cases}\partial_{t}u^{2}\pm\frac{1}{2}a(x,z)\sum_{j=1}^{d}\partial_{x_{j}}u^{2}=0,\\\
u^{2}(0,x,z)=u^{2}_{0}(x,z).\end{cases}$ (50)

Once can then solve this equation similarly as the original advection equation
to obtain the ensemble average of $u^{2}$.

Among the quantities of physical interest include $|u|$ (amplitude), $u^{2}$
(energy) and $\nabla u$ (flux). To compute the ensemble average of the
amplitude and energy, we can use

$\displaystyle\bar{u}(t,x)=\frac{1}{M}\sum_{m=1}^{M}u(t,x,z_{m})=\int\int{W}(t,x,p,q)\,dp\,dq\approx\bar{u}_{n,j}\equiv\frac{1}{N_{p}^{L}N_{q}^{L}}\sum_{k}^{N_{p}^{L}}\sum_{l}^{N_{q}^{L}}W_{n,j,k,l}$
(51)

where $W_{n,j,k,l}=W(t=n\Delta t,x=j\Delta x,p=k\Delta p,q=l\Delta q)$. The
discretization meshes and corresponding notations in $x,t$ are the same as for
the linear heat equation. The discretisation of $p$ and $q$, in one dimension
for example, is $p_{k}=k\Delta p$ and $q_{l}=l\Delta q$ for $k=0,...,N_{p}$,
$l=0,...,N_{q}$, where $N_{p}\Delta p=1=N_{q}\Delta q$. In $d$ spatial
dimensions, the total number of spatial mesh points for $N_{p}$ and $N_{q}$
are respectively $N_{p}^{L}$ and $N_{q}^{L}$ since there are both
$L$-dimensional. We will use the forward Euler discretisation for time, center
finite difference in $x$, and upwind scheme in $w$ (so the overall spatial
error is of $O(1/N)$ and temporal error is of $O(\Delta t)$. See Appendix A.3
for details of the discretisation scheme and the corresponding matrix
$\mathcal{M}$ for the problem.

Denote the solution vector of Eq. (48) by
$\boldsymbol{W}=\sum_{n}^{N_{t}}\sum_{j}^{N^{d}}\sum_{k}^{N_{p}^{L}}\sum_{l}^{N_{q}^{L}}W_{n,j,k,l}|n\rangle|j\rangle|k\rangle|l\rangle$.
Then
$\boldsymbol{W}^{T}\mathcal{G}^{density}\boldsymbol{W}=N_{p}^{L}N_{q}^{L}|\bar{u}_{N_{t},J}|^{2}$,
where $\mathcal{G}^{density}=|G^{density}_{N_{t},J}\rangle\langle
G^{density}_{N_{t},J}|$ and
$|G^{density}_{N_{t},J}\rangle=(1/\sqrt{N_{p}^{L}N_{q}^{L}})\sum_{k}^{N_{p}^{L}}\sum_{l}^{N_{q}^{L}}|N_{t}\rangle|J\rangle|k\rangle|l\rangle$.
The quantum state that embeds the initial condition of $\boldsymbol{W}$ with
sparsity $\sigma_{W_{0}}$ is defined as

$\displaystyle|W_{0}\rangle=\frac{1}{\mathcal{N}_{W_{0}}}\sum_{j}^{N^{d}}\sum_{k}^{N_{p}^{L}}\sum_{l}^{N_{q}^{L}}W_{0,j,k,l}|0\rangle|j\rangle|k\rangle|l\rangle$
(52)

where the normalisation constant is
$\mathcal{N}_{W_{0}}^{2}=\sum_{j}^{N^{d}}\sum_{k}^{N_{p}^{L}}\sum_{l}^{N_{q}^{L}}|W_{0,j,k,l}|^{2}$.

######  Lemma 17.

The upper bound to the normalisation constant
$n^{2}_{W_{0}}\equiv\mathcal{N}^{2}_{W_{0}}/(N_{p}^{L}N_{q}^{L})$ has range
$O(1)\leq n^{2}_{W_{0}}\leq O(N^{d})$

###### Proof.

First, similar to the proof in Lemma 7, $\int dp\int dq|W(0,x,p,q)|^{2}\leq
O(1)$. Then using the quadrature rule one gets $\int dx\int dp\int
dq|W(0,x,p,q)|^{2}\approx(1/N^{d})\sum_{j}^{N^{d}}\int dp\int dq|W(0,x=j\Delta
x,p,q)|^{2}\leq O(1)$. Then $O(1)\geq\int dx\int dp\int
dq|W(0,x,p,q)|^{2}\approx\mathcal{N}^{2}_{W_{0}}/(N^{d}N_{p}^{L}N_{q}^{L})=n^{2}_{W_{0}}/N^{d}$.
Thus $O(1)\leq n^{2}_{W_{0}}\leq O(N^{d})$. However, the quantum inner product
$\Upsilon$ is upper bounded by $1$, where
$1\geq\Upsilon=\boldsymbol{W}^{T}\mathcal{G}^{density}\boldsymbol{W}/\mathcal{N}^{2}_{W_{0}}=O(N_{p}^{L}N_{q}^{L}/\mathcal{N}^{2}_{W_{0}})=O(1/n^{2}_{W_{0}})$
since $\bar{u}^{2}(t,x)=O(1)$. ∎

We can then proceed as before using Lemmas 3 and 6 to find the total cost in
computing the density $|\bar{u}_{N_{t},J}|^{2}$ at $\Lambda$ mesh points.

######  Theorem 18.

A quantum algorithm that takes sparse access to $\mathcal{M}$ (using an
$r^{\text{th}}$-order method with $r\geq 1$) is able to estimate the energy
$|\bar{u}_{N_{t},J}|^{2}$ at $\Lambda$ points to precision $\epsilon$ with an
upper bound on the gate complexity

$\displaystyle\mathcal{Q}=\mathcal{O}\left(n_{W_{0}}^{2}\Lambda(Ld)^{3}\left(\frac{Ld}{\epsilon}\right)^{1+9/r}\right)$
(53)

with a smaller query complexity.

###### Proof.

Here $\mathcal{M}$ has sparsity $s=O(Ld)$ and condition number
$\kappa=O(LdN^{3})=O(Ld(Ld/\epsilon)^{3/r})$ where $N=O((Ld/\epsilon)^{3/r})$
from Appendix B.3. Then Lemma 3 can be applied to show a query complexity
$\tilde{\mathcal{O}}(s\kappa^{3}\mathcal{N}^{2}_{W_{0}}/(N_{p}^{L}N_{q}^{L}\epsilon))=\tilde{\mathcal{O}}(n^{2}_{W_{0}}(Ld)^{3}(Ld/\epsilon)^{1+9/r})$
with the same order of additional $2$-qubit gates, where $\epsilon$ is the
error in $|\bar{u}_{N_{t},J}|^{2}$. ∎

######  Lemma 19.

When $M<M_{adv}\equiv O(L^{(d+2L+3)/r+2}(d/\epsilon)^{(2L+2)/r})$ the
classical algorithm for the problem has the cost
$\mathcal{C}=O(Md^{2}(d/\epsilon)^{(d+1)/r})$. When $M>M_{adv}$, the classical
algorithm has cost $\mathcal{C}=O((Ld)^{2}(Ld/\epsilon)^{(d+2L+3)/r})$.

###### Proof.

See Appendix B.3 for details. ∎

######  Corollary 20.

To attain a quantum advantage $M<M_{adv}$, it is sufficient for the following
condition to hold

$\displaystyle{\scriptstyle\mathcal{O}}\left(\frac{M}{\Lambda
n^{2}_{W_{0}}L^{4+9/r}d}\left(\frac{d}{\epsilon}\right)^{(d-8)/r-1}\right)=\tilde{O}(1).$
(54)

When $M>M_{adv}$ it is sufficient for the following condition to hold

$\displaystyle{\scriptstyle\mathcal{O}}\left(\frac{1}{n^{2}_{W_{0}}\Lambda
Ld}\left(\frac{Ld}{\epsilon}\right)^{(d+2L-6)/r-1}\right)=\tilde{O}(1).$ (55)

One can see from Eq. (54) that quantum advantage with respect to $M$ is always
possible in the range $M^{\prime}_{adv}<M<M_{adv}$ where
$M^{\prime}_{adv}=\tilde{\mathcal{O}}(n^{2}_{W_{0}}\Lambda
L^{4+9/r}d(\epsilon/d)^{(d-8)/r-1})$, with up to exponential quantum advantage
in $d$ and $\epsilon$. As a simple example, if we begin with a point source
initial condition and only require the final solution at $O(1)$ points, then
$n^{2}_{W_{0}}=O(1)=\Lambda$. When $d>8+r$, the range of possible $M$,
captured by $M_{adv}-M^{\prime}_{adv}$, is very large. A necessary condition
for such an $M$ to exist in more general cases is
$O(1)<n^{2}_{W_{0}}\Lambda<\tilde{\mathcal{O}}(L^{(d+L-6)/r-2}(d/\epsilon)^{(d+L-6)/r-1}/d)$,
which allows for a wide range of possibilities.

In the case when $M>M_{adv}$, one sees from Eq. (55) there is no quantum
advantage with respect to $M$, but there is quantum advantage in
$L,d,\epsilon$ when
$n^{2}_{W_{0}}\Lambda<\tilde{\mathcal{O}}((Ld/\epsilon)^{(d+2L-6)/r-1}/(Ld))$.
Depending on $n^{2}_{W_{0}}\Lambda$, up to exponential quantum advantage in
$L,d,\epsilon$ is possible.

###  III.4 The Schrödinger equation

Consider the following Schrödinger equation with uncertainty

$\begin{cases}i\hbar\,\partial_{t}u=-\frac{\hbar^{2}}{2}\Delta u+a(x,z)u\,,\\\
u(0,x,z)=u_{0}(x,z)\,,\end{cases}$ (56)

where $u=u(t,x,z)\in\mathbb{C}$ is the complex-valued wave function,
$a(x,z)>0$ given in (1) is the potential, which could be uncertain due to
modeling or computational errors, or random media [35]. $\hbar$ is Planck’s
constant. For $L=1$, if we use the transformation

$U(t,x,z,p)=a(z)e^{-pa(z)}u(t,x,z)$ (57)

where $p\geq 0$, then $U$ solves

$\begin{cases}i\hbar\,\partial_{t}U=-\frac{\hbar^{2}}{2}\Delta
U-b(x)\partial_{p}U\,,\\\ U(0,x,z,w)=a(z)e^{-pa(z)}u(0,x,z)\,.\end{cases}$
(58)

However, this problem is again ill-posed. For example if $b(x)$ is a constant,
by taking a Fourier transform on both $x$ and $p$ and assume the Fourier
variable is $\xi\in\mathbb{R}^{d}$ in $x$ and $\eta\in\mathbb{R}$ in $p$, then

$i\hbar\partial_{t}\hat{\hat{U}}=\left[\frac{1}{2}\hbar^{2}|\xi|^{2}+ib\eta\right]\hat{\hat{U}},$
(59)

which gives

$\partial_{t}\hat{\hat{U}}=\left[-\frac{i}{2}\hbar|\xi|^{2}+\frac{\eta}{\hbar}\right]\hat{\hat{U}}.$
(60)

Clearly this equation is unstable since $\eta$ may be positive.

Like for the advection equation, for the Schrödinger equation we also cannot
use the transformation applied to the linear heat and Boltzmann equations by
extending $p_{i}$ to $(-\infty,\infty)$. Like for the advection equation, we
define

${W}(t,x,p,q)=\frac{1}{M}\sum_{m=1}^{M}\prod_{i=1}^{L}\sqrt{a_{i}(z_{m})}\,e^{-\sqrt{a_{i}(z_{m})}(p_{i}+q_{i})}u(t,x,z_{m})\,.$
(61)

Then $W$ satisfies

$\displaystyle\begin{cases}i\hbar\partial_{t}{W}=-\frac{\hbar^{2}}{2}\Delta{W}+\sum_{i=1}^{L}b_{i}(x)\partial_{p_{i}p_{i}}{W},\\\
{W}(0,x,p,q)=\frac{1}{M}\sum_{m=1}^{M}\prod_{i=1}^{L}\sqrt{a_{i}(z_{m})}e^{-(p_{i}+q_{i})\sqrt{a_{i}(z_{m})}}u(0,x,z_{m}),\\\
\partial_{p_{i}}{W}-\partial_{q_{i}}{W}\Big{|}_{p_{i}=0}=0\quad{\text{for
all}}\,\,i\,.\end{cases}$ (62)

Note Eq. (62) is the free Schrödinger equation in the $(x,p)$ space, which is
a good equation and the well-posedness of the initial value problem is
classical. The average of the solutions of the original problem, Eq. (56), can
be recovered from $W(t,x,p,q)$ using

$\displaystyle\overline{u}(t,x)$
$\displaystyle=\frac{1}{M}\sum_{m=1}^{M}u(t,x,z_{m})=\int\int{W}(t,x,p,q)\,dp\,dq\,.$
(63)

The quantities of interest in real applications are the important physical
observables of the original Schrödinger equation, which include position
density $|u|^{2}$, flux or moment density $\hbar\,\text{Im}(\overline{u}\nabla
u)$, kinetic energy $\frac{1}{2}\hbar^{2}|\nabla u|^{2}$, and total energy
$\frac{1}{2}\hbar^{2}|\nabla u|^{2}+a(x,z)|u|^{2}$.

The simplest method to extract the ensemble average of the density
$|\bar{u}(t,x)|^{2}=(1/M^{2})|\sum_{m=1}^{M}u(t,x,z_{m})|^{2}$ is from
observing $\bar{u}(t,x)=\int\int W(t,x,p,q)dpdq$. This integral can be
approximated using the quadrature rule

$\displaystyle|\bar{u}(t=N_{t}\Delta t,x=J\Delta
x)|\approx|\bar{u}_{N_{t},J}|=\frac{1}{N_{p}^{L}N_{q}^{L}}\sum_{k}^{N_{p}^{L}}\sum_{l}^{N_{q}^{L}}W_{N_{t},J,k,l}$
(64)

where $W_{N_{t},J,k,l}=W(t=N_{t}\Delta t,x=J\Delta x,p=k\Delta p,q=l\Delta
q)$. The discretization meshes and corresponding notations in $x,t$, $p$ and
$q$ are the same as for the linear advection equation. We use the center
difference scheme in $x$ and $p,q$, and forward Euler in $t$ to solve Eq. (62)
to $t=1$. See Appendix A.4 for details of the discretisation scheme and the
corresponding matrix $\mathcal{M}$ for the problem.

Denote the solution vector of Eq. (62) by
$\boldsymbol{W}=\sum_{n}^{N_{t}}\sum_{j}^{N^{d}}\sum_{k}^{N_{p}^{L}}\sum_{l}^{N_{q}^{L}}W_{n,j,k,l}|n\rangle|j\rangle|k\rangle|l\rangle$.
Then
$\boldsymbol{W}^{T}\mathcal{G}^{density}\boldsymbol{W}=N_{p}^{L}N_{q}^{L}|\bar{u}_{N_{t},J}|^{2}$,
where $\mathcal{G}^{density}=|G^{density}_{N_{t},J}\rangle\langle
G^{density}_{N_{t},J}|$ and
$|G^{density}_{N_{t},J}\rangle=(1/\sqrt{N_{p}^{L}N_{q}^{L}})\sum_{k}^{N_{p}^{L}}\sum_{l}^{N_{q}^{L}}|N_{t}\rangle|J\rangle|k\rangle|l\rangle$.
The quantum state that embeds the initial condition of $\boldsymbol{W}$ with
sparsity $\sigma_{W_{0}}$ is defined as

$\displaystyle|W_{0}\rangle=\frac{1}{\mathcal{N}_{W_{0}}}\sum_{j}^{N^{d}}\sum_{k}^{N_{p}^{L}}\sum_{l}^{N_{q}^{L}}W_{0,j,k,l}|0\rangle|j\rangle|k\rangle|l\rangle$
(65)

where the normalisation constant is
$\mathcal{N}_{W_{0}}^{2}=\sum_{j}^{N^{d}}\sum_{k}^{N_{p}^{L}}\sum_{l}^{N_{q}^{L}}|W_{0,j,k,l}|^{2}$.

######  Lemma 21.

The upper bound to the normalisation constant
$n^{2}_{W_{0}}\equiv\mathcal{N}^{2}_{W_{0}}/(N_{p}^{L}N_{q}^{L})$ has range
$O(1)\leq n^{2}_{W_{0}}\leq O(N^{d})$

###### Proof.

First, similar to the proof in Lemma 7, $\int\int dpdq|W(0,x,p,q)|^{2}\leq
O(1)$. Then using the quadrature rule we have $\int dx\int dp\int
dq|W(0,x,p,q)|^{2}\approx(1/N^{d})\sum_{j}^{N^{d}}\int dp\int dq|W(0,x=j\Delta
x,p,q)|^{2}\leq O(1)$. Then $O(1)\geq\int dx\int dp\int
dq|W(0,x,p,q)|^{2}\approx\mathcal{N}^{2}_{W_{0}}/(N^{d}N_{p}^{L}N_{q}^{L})$.
Thus $O(1)\leq\mathcal{N}^{2}_{W_{0}}\leq O(N_{p}^{L}N_{q}^{L}N^{d})$. ∎

We can then proceed as before using Lemmas 3 and 6 to find the total cost in
computing the density $|\bar{u}_{N_{t},J}|^{2}$ at $\Lambda$ mesh points.

######  Theorem 22.

A quantum algorithm that takes sparse access to $\mathcal{M}$ (using an
$r^{\text{th}}$-order method with $r\geq 1$) is able to estimate the density
$|\bar{u}_{N_{t},J}|^{2}$ at $\Lambda$ points to precision $\epsilon$ with an
upper bound on the gate complexity

$\displaystyle\mathcal{Q}=\tilde{\mathcal{O}}\left(n_{W_{0}}^{2}\Lambda(d+L)^{3}\left(\frac{d+L}{\epsilon}\right)^{1+6/r}\right)$
(66)

with a smaller query complexity.

###### Proof.

Here $\mathcal{M}$ has sparsity $s=O(d+L)$ and condition number
$\kappa=O((d+L)N^{2})=O((d+L)((d+L)/\epsilon)^{2/r})$ where
$N=O(((d+L)/\epsilon)^{1/r})$ from Appendix B.4. Then Lemma 3 can be applied
to show a query complexity
$\tilde{\mathcal{O}}(s\kappa^{3}\mathcal{N}^{2}_{W_{0}}/(N_{p}^{L}N_{q}^{L}\epsilon))=\tilde{\mathcal{O}}(n^{2}_{W_{0}}(d+L)^{3}((d+L)/\epsilon)^{1+6/r})$
with the same order of additional $2$-qubit gates, where $\epsilon$ is the
error in $|\bar{u}_{N_{t},J}|^{2}$. ∎

######  Lemma 23.

When $M<M_{Schr}\equiv
O((d+L)^{(d+2L+2)/r+2}/(d^{(d+2)/r+2}\epsilon^{2L/2}))$, the classical
algorithm for the problem has the cost
$\mathcal{C}=O(Md^{2}(d/\epsilon)^{(d+2)/r})$. When $M>M_{Schr}$, the
classical algorithm has cost
$\mathcal{C}=O((d+L)^{2}((d+L)/\epsilon)^{(d+2L+2)/r})$.

###### Proof.

See Appendix B.4 for details. ∎

######  Corollary 24.

To attain a quantum advantage when $M<M_{Schr}$, it is sufficient for the
following condition to hold

$\displaystyle{\scriptstyle\mathcal{O}}\left(\frac{Md^{(d+2)/r+2}}{n^{2}_{W_{0}}\Lambda(d+L)^{4+6/r}\epsilon^{(d-4)/r-1}}\right)=\tilde{O}(1).$
(67)

When $M>M_{Schr}$ it is sufficient for the following to hold

$\displaystyle{\scriptstyle\mathcal{O}}\left(\frac{(d+L)^{(d+2L-4)/r-2}}{n^{2}_{W_{0}}\Lambda\epsilon^{(d+2L-4)/r-1}}\right)=\tilde{O}(1).$
(68)

One sees from Eq. (67) that quantum advantage with respect to $M$ is always
possible in the range $M^{\prime}_{Schr}<M<M_{Schr}$ where
$M^{\prime}_{Schr}=\tilde{\mathcal{O}}(n^{2}_{W_{0}}\Lambda(d+L)^{4+6/r}\epsilon^{(d-4)/r-1}/d^{(d+2)/r+2})$,
with up to exponential quantum advantage in $d$ and $\epsilon$. As a simple
example, if one begins with a point source initial condition and only requires
the final solution at $O(1)$ points, then $n^{2}_{W_{0}}=O(1)=\Lambda$. When
$d>c+4$, the range of possible $M$, captured by $M_{Schr}-M^{\prime}_{Schr}$,
is very large. A necessary condition for such an $M$ to exist in more general
cases is
$O(1)<n^{2}_{W_{0}}\Lambda<\tilde{\mathcal{O}}((d+L)^{(d+L-4)/r+2}/(\epsilon^{L/2+(d-4)/r}))$,
which allows for a wide range of possibilities.

In the case when $M>M_{Schr}$, one sees from Eq. (68) there is no quantum
advantage with respect to $M$, but there is quantum advantage in
$L,d,\epsilon$ when
$n^{2}_{W_{0}}\Lambda<\tilde{\mathcal{O}}((d+L)^{(d+2L-4)/r-2}/(\epsilon^{(d+2L-4)/r-1}))$.
Depending on $n^{2}_{W_{0}}\Lambda$, up to exponential quantum advantage in
$L,d,\epsilon$ is possible.

##  IV Discussion

We remark that it is also possible to compute ensemble averaged observables
from the original PDEs directly with a corresponding quantum algorithm, where
the PDE needs to be solved $M$ times for $M$ samples of $z$. In this case,
however, different oracle assumptions are required. To solve the original PDE
with $M$ initial conditions, there are $M$ different corresponding
$\mathcal{M}$ matrices, since $a(x,z)$ appears explicitly in the PDE to to
solved. Therefore $M$ different sparse oracles are required. This is in
contrast to solving the phase space PDEs, where only a single oracle is
needed, which is furthermore independent of the details of the stochastic
model $a(x,z)$! This can be relevant for scenarios where either the sparse
oracles for $\mathcal{M}$ are not easily programmable to high precision or one
is not given access to be able to tune the sparse oracles for different values
of $z$. To solve the phase space problem, only a single type of oracle to
access $\mathcal{M}$ is required for each type of PDE, rather than a different
oracle for every different instance of the PDE.

However, in cases that multiple sparse oracles for $\mathcal{M}$, one for each
value of $z$, is assumed to be given without any cost, then solving the
original PDE with a quantum algorithm with cost $\mathcal{Q}_{orig}$ is only
advantageous compared to our new algorithms only for relatively smaller values
of $M$, given in Table 2. For modest sizes of $M$, we see that our new quantum
algorithms are still preferred, i.e. $\mathcal{Q}<\mathcal{Q}_{orig}$. See
Appendix C for details.

We note that the conditions in Table 1 and Table 2 do not assume costs in
preparing the initial quantum state. As discussed in Section III it is
straightforward to include the initial state preparation cost in Table 1 by
multiplying the quantum cost by the sparsity $\sigma_{0}$ of the initial state
and to include another factor of $(d+L)$, following Lemma 6. Here if we define
the input-output factor (IOF) as
$\Gamma\equiv\sigma_{0}n_{0}^{2}\Lambda=N^{b}$, then the minimum $b=L$,
instead of $b=0$ in the current table. We see that this still leaves room for
quantum advantage usually when $d$ is larger than $L$.

It is crucial to emphasise that here the initial state preparation costs (and
hence the IOF) is independent of $M$. Neither the sparsity, the dimension nor
the normalisation of the initial quantum state depends on $M$. This is because
each amplitude of the initial quantum state is already an average over $M$
terms associated with each initial condition, and this amplitude is first
computed classically before being embedded in a quantum state. This means that
the end-to-end cost of our quantum algorithms in Table 1 are all independent
of $M$.

Suppose we now include initial state preparation costs in comparing
$\mathcal{Q}$ with $\mathcal{Q}_{orig}$ in Table 2. We note that their
respective IOF factors $\Gamma\equiv\sigma_{0}n_{0}^{2}\Lambda$ are different
because the sparsities of their initial condition is different up to a factor
of $N^{L}$. For instance, for a point source in the original problem, where
the sparsity of of the initial state is $1$, the sparsity of the initial
quantum state for the phase space equation is $N^{L}$. If we can assume that
in solving the original PDE $M$ times we are given all the corresponding
oracles for $\mathcal{M}$ for free – which is not always possible – and if we
choose the same initial state preparation scheme in Lemma 6 as before, the
quantum algorithm to solve the original PDE $M$ times is preferred unless we
have very large $M$. Alternatively, if we choose a different initial
preparation scheme like in [29], then if we count only gate complexity, we
only need to multiply the quantum cost by a logarithmic factor of
$\sigma_{0}$, although order $\sigma_{0}(d+L)$ ancilla qubits are now
necessary. In this case, the conditions in Table 2 still hold up to an order
of $\tilde{\mathcal{O}}(L(d+L))$ when including initial state preparation.

Table 2: Here $\mathcal{Q}_{orig}$ and $\mathcal{Q}$ are the respective quantum costs for $r^{\text{th}}$-order methods in computing ensemble averaged observables, over $M$ samples, based on solving the original equation ($\mathcal{Q}_{orig}$) versus the phase space representation ($\mathcal{Q}$). Here $s$ and $\kappa$ are the sparsity and condition numbers corresponding to $\mathcal{M}$ of the original $(d+1)$-dimensional PDE and the phase space PDE. However, there are more oracle assumptions needed for $\mathcal{Q}_{orig}$ and when these assumptions are not obeyed, the phase space method is always preferable. See main text for discussion. PDE |  $s$ |  $\kappa$ |  $\mathcal{Q}<\mathcal{Q}_{orig}$
---|---|---|---
Linear heat |  $d$ |  $(d/\epsilon)^{2/r}$ |  $M>\tilde{\mathcal{O}}(L^{4+9/r}d^{3}(d/\epsilon)^{3/r})$
Phase space linear heat |  $Ld$ |  $Ld(Ld/\epsilon)^{3/r}$ |
Linear Boltzmann |  $(d/\epsilon)^{d/r}$ |  $d(d/\epsilon)^{1/r}$ |  $M>\tilde{\mathcal{O}}(L^{1+(3+d)/r}(d/\epsilon)^{2/r}/(d\epsilon))$
Phase space linear Boltzmann |  $L(Ld/\epsilon)^{d/r}$ |  $(d+L)(Ld/\epsilon)^{1/r}$ |
Linear advection |  $d$ |  $d(d/\epsilon)^{1/r}$ |  $M>\tilde{\mathcal{O}}(L^{4+9/r}d^{2}(d/\epsilon)^{8/r})$
Phase space linear advection |  $Ld$ |  $Ld(Ld/\epsilon)^{3/r}$ |
Schrödinger |  $d$ |  $(d/\epsilon)^{2/r}$ |  $M>\tilde{\mathcal{O}}((d+L)^{3}(1+L/d)^{1+2/r}((L+d)/\epsilon)^{4/r})$
Phase space Schrödinger |  $d+L$ |  $(d+L)((d+L)/\epsilon)^{2/r}$ |

We observe that while classical random walk methods can also in principle be
used to solve the phase space versions of these PDEs and these algorithms do
not suffer from the curse of dimensionality, we don’t choose to compare our
quantum algorithms to these classical costs. Firstly, quantum advantage
compared to classical random walk methods are still possible, albeit are at
most polynomial and not exponential, for example [27]. However, only advantage
in $d$, $\epsilon$ are usually discussed, whereas here our main focus is on a
different type of advantage with respect to $M$ and $L$, while still having
potential advantage in $d$ and $\epsilon$. Secondly, it is well-known that
classical random walk methods can be prone to a lot of statistical noise and
also have a lower order of accuracy. We consider it more appropriate therefore
to compare the quantum algorithm with its direct classical counterpart, i.e.
using finite difference methods, which can offer highly accurate and high
resolution solutions if competitive high-order finite difference methods are
used.

##  V Conclusion

We have introduced new quantum algorithms that can compute ensemble averaged
observables of PDEs with uncertainty. The end-to-end cost of these algorithms
(including initial state preparation and final measurement) can be independent
of the number of initial conditions $M$ and in certain regimes can also show
advantage in $L$, $d$ and $\epsilon$.

The key idea of ours is to develop new transformations that map PDEs with
uncertainty to deterministic PDEs. We transform the original PDEs with $M$
initial conditions into phase space PDEs in higher dimensions, so the
transformed equations are independent of the uncertain coefficients with a
single initial condition. Although solving these new PDEs classically – while
also guaranteeing high accuracy and resolution – would be inefficient in $d$
and $\epsilon$, quantum algorithms are efficient in $d$, $\epsilon$ and also
provide advantage in $M$ and $L$ in different regimes.

These transformations can also be extended to certain nonlinear PDEs using
techniques in [18], which introduced another way of embedding the ensemble
average over $M$ in the initial condition–to deal with the nonlinearity rather
than the uncertainty in this case. This will be explored in future work.

These new algorithms offer exciting opportunities to begin exploring how
quantum methods can be applied to uncertainty quantification.

## Acknowledgement

SJ was partially supported by the NSFC grant No. 12031013, the Shanghai
Municipal Science and Technology Major Project (2021SHZDZX0102), and the
Innovation Program of Shanghai Municipal Education Commission (No.
2021-01-07-00-02-E00087). NL acknowledges funding from the Science and
Technology Program of Shanghai, China (21JC1402900), the Shanghai Pujiang
Talent Grant (no. 20PJ1408400) and the NSFC International Young Scientists
Project (no. 12050410230).

## Appendix

##  Appendix A Numerical discretisation of PDEs, $\mathcal{M}$ matrices and
their properties

###  A.1 Linear heat equation

####  A.1.1 Numerical discretizations

As an example we use center finite difference scheme to solve the phase space
heat equation (8). For convenience we use homogeneous boundary condition in
$x$, and the case $b_{i}(x)>0$, so the zero boundary condition in $u$ is given
on the right boundary. We also just consider $p>0$, since the inclusion of the
case of $p<0$ will not affect the computation of condition number and
sparsity. We first consider spatially one-dimensional (for both $x$ and $w$)
problem. Consider $N_{t}+1$ steps in time $0=t_{0}<t_{1}<\cdots<t_{N_{t}}=1$
and $N_{x}+1$ spatial mesh points $0<x_{0}<x_{1}<\cdots<x_{N_{x}}=1$,
$0<p_{0}<p_{1}<\cdots<p_{N_{x}}=1$ by setting $t_{n}=n\tau$ and
$x_{j}=p_{j}=jh$, where $\tau=1/N_{t}$ and $h=1/N_{x}$. Let $u^{n}_{jk}$
denote the numerical approximation of $u$ at $(t_{n},x_{j},p_{k})$. We use the
forward Euler method in time, center difference in $x$, and upwind scheme in
$p$ (so the overall spatial error is of $O(h)$):

$\displaystyle\frac{u^{n+1}_{jk}-u^{n}_{jk}}{\tau}+\frac{1}{h^{3}}\left[(u_{j-1,k+1}^{n}-2u_{j,k+1}^{n}+u_{j+1,k+1}^{n})-(u_{j-1,k}^{n}-2u_{j,k}^{n}+u_{j+1,k}^{n})\right]=0,$
(69) $\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad
j=1,\cdots,N_{x}-1;\quad k=0,\cdots,N_{x}-1.$

Define

$L_{h}=\begin{bmatrix}-2&1&&&\\\ 1&-2&\ddots&&\\\ &\ddots&\ddots&\ddots&\\\
&&\ddots&-2&1\\\ &&&1&-2\\\ \end{bmatrix}_{(N_{x}-1)\times(N_{x}-1)},\qquad
b(k)=\begin{bmatrix}u_{0}(k+1)-u_{0}(k)\\\ 0\\\ \vdots\\\ 0\\\
u_{N_{x}}(k+1)-u_{N_{x}}(k)\\\ \end{bmatrix},$

and set $\lambda=\tau/h^{3}$,
$u_{k}^{n}=(u_{1,k}^{n},\cdots,u_{N_{x},k}^{n})^{T}$, then (69) can be written
as

$u_{k}^{n+1}-\left(I+\lambda L_{h}\right)u_{k}^{n}+\lambda
L_{h}u_{k+1}^{n}=\lambda b(k),$ (70)

which can be marched forward in time by solving a linear system for a
classical computer.

Let $u^{n}=(u_{0}^{n},\cdots,u_{N_{x}-1}^{n})^{T}$. Using the boundary
condition $u^{n}_{N_{x}}=0$, then (70) can be further written as

$-Bu^{n}+u^{n+1}=f^{n+1}=:\lambda(b_{0},\cdots,b_{N_{x}-1})^{T},$ (71)

where

$B=\begin{bmatrix}I+\lambda L_{h}&-\lambda L_{h}&&&\\\ &I+\lambda
L_{h}&-\lambda L_{h}&&\\\ &&\ddots&\ddots&\\\ &&&\ddots&-\lambda L_{h}\\\
&&&&I+\lambda L_{h}\\\ \end{bmatrix}.$

By introducing the notation $U=[u^{1},\cdots,u^{N_{t}}]^{T}$, one obtains the
following linear system

$LU=F,$ (72)

where

$L=\begin{bmatrix}I&&&\\\ -B&I&&\\\ &\ddots&\ddots&\\\ &&-B&I\\\
\end{bmatrix},\qquad F=\begin{bmatrix}f^{1}+Bu^{0}\\\ f^{2}\\\ \vdots\\\
f^{N_{t}}\\\ \end{bmatrix}.$

Note that the $\mathcal{M}$ matrix is denoted

$H=\begin{bmatrix}O&L\\\ L^{T}&O\end{bmatrix}.$

####  A.1.2 Estimation of the eigenvalues

Let $\lambda$ be an eigenvalue of $H$, i.e.,

$\mbox{det}(\lambda I-H)=\mbox{det}(\lambda^{2}I-LL^{T})=:\mbox{det}(\mu
I-LL^{T}),$

where $\mu=\lambda^{2}$ is the eigenvalue of $LL^{T}$, or $\mu^{1/2}$ is the
singular value of $L$. A direct calculation gives

$LL^{T}=\begin{bmatrix}I&-B&&\\\ -B&I+B^{2}&\ddots&\\\ &\ddots&\ddots&-B\\\
&&-B&I+B^{2}\\\ \end{bmatrix}.$

We establish the upper and lower bounds of the eigenvalues of $LL^{T}$ by
using the Gershgorin circle theorem.

######  Lemma 25.

Let $\lambda=\tau/h^{3}$. Then the minimum and maximum eigenvalues of $H$
satisfy

$|\lambda|_{\min}\geq\tau,\qquad|\lambda|_{\max}\leq 2.$

Consequently, for matrix $H$,
$s=7,\kappa\leq\frac{2}{\tau}=O(N_{t})=O(N_{x}^{3})$.

###### Proof.

Let $P^{-1}L_{h}P=\Lambda$, where $\Lambda$ is the diagonal matrix consisting
of the eigenvalues of $L_{h}$. Then

$B\sim\Pi_{B}:=\begin{bmatrix}I+\lambda\Lambda&-\lambda\Lambda&&&\\\
&I+\lambda\Lambda&-\lambda\Lambda&&\\\ &&\ddots&\ddots&\\\
&&&\ddots&-\lambda\Lambda\\\ &&&&I+\lambda\Lambda\\\ \end{bmatrix}$

and

$B^{2}\sim\Pi_{B^{2}}=\begin{bmatrix}(I+\lambda\Lambda)^{2}&-2\lambda\Lambda(I+\lambda\Lambda)&\lambda^{2}\Lambda^{2}&&\\\
&\ddots&\ddots&\ddots&\\\
&&(I+\lambda\Lambda)^{2}&-2\lambda\Lambda(I+\lambda\Lambda)&\lambda^{2}\Lambda^{2}\\\
&&&(I+\lambda\Lambda)^{2}&-2\lambda\Lambda(I+\lambda\Lambda)&\\\
&&&&(I+\lambda\Lambda)^{2}\\\ \end{bmatrix}.\\\ $

Let $\widetilde{P}=\mbox{diag}(P,\cdots,P)$. Then

$\widetilde{P}^{-1}(LL^{T})\widetilde{P}=\begin{bmatrix}I&-\Pi_{B}&&&\\\
-\Pi_{B}&I+\Pi_{B^{2}}&-\Pi_{B}&&\\\ &\ddots&\ddots\ddots&\\\
&&\ddots\ddots&-\Pi_{B}\\\ &&-\Pi_{B}&I+\Pi_{B^{2}}\\\ \end{bmatrix}.$

The eigenvalues of $L_{h}$ are

$\nu_{l}=-4\sin^{2}\frac{l\pi}{2N_{x}}=-4\sin^{2}\frac{l\pi h}{2},\quad
l=1,\cdots,N_{x}-1$ (73)

Noting that the similarity transformation does not change the eigenvalues,
thus one can apply the Gershgorin circle theorem to estimate the maximum
eigenvalues of this matrix. In particular,

$\displaystyle\mu_{\max}\leq\max_{l}\\{1+(1+\lambda\nu_{l})^{2}+2(|1+\lambda\nu_{l}|+\lambda|\nu_{l}|)|+2\lambda|\nu_{l}||1+\lambda\nu_{l}|+\lambda^{2}\nu_{l}^{2}\\}\leq
4,$ (74)

where we need the assumption $\lambda\leq 1/4$.

By definition,

$\sigma_{\min}(L)=\frac{1}{\sigma_{\max}(L^{-1})}.$

By simple algebra, one has

$L^{-1}=\begin{bmatrix}I&&&&\\\ B&I&&&\\\ B^{2}&B&I&&\\\
\vdots&\ddots&\ddots&\ddots&\\\ B^{N_{t}-1}&\cdots&&B&I\\\
\end{bmatrix}=\begin{bmatrix}I&&&&\\\ &I&&&\\\ &&\ddots&&\\\ &&&\ddots&\\\
&&&&I\\\ \end{bmatrix}+\begin{bmatrix}&&&&\\\ B&&&&\\\ &B&&&\\\ &&\ddots&&\\\
&&&B\end{bmatrix}+\cdots,$ (75)

which gives

$\displaystyle\sigma_{\max}(L^{-1})$
$\displaystyle=\|L^{-1}\|_{2}\leq\|I\|_{2}+\|B\|_{2}+\cdots+\|B^{N_{t}-1}\|_{2}$
$\displaystyle\leq\|I\|_{2}+\|B\|_{2}+\|B\|_{2}^{2}+\cdots+\|B\|_{2}^{N_{t}-1},$

Moreover,

$BB^{T}\sim\Pi_{B_{2}}=\begin{bmatrix}(I+\lambda\Lambda)^{2}+\lambda^{2}\Lambda^{2}&-\lambda\Lambda(I+\lambda\Lambda)&&&\\\
-\lambda\Lambda(I+\lambda\Lambda)&(I+\lambda\Lambda)^{2}+\lambda^{2}\Lambda^{2}&-\lambda\Lambda(I+\lambda\Lambda)&&\\\
&\ddots&\ddots&\ddots&\\\
&&-\lambda\Lambda(I+\lambda\Lambda)&(I+\lambda\Lambda)^{2}+\lambda^{2}\Lambda^{2}&-\lambda\Lambda(I+\lambda\Lambda)\\\
&&&-\lambda\Lambda(I+\lambda\Lambda)&(I+\lambda\Lambda)^{2}+\lambda^{2}\Lambda^{2}\\\
&&&&(I+\lambda\Lambda)^{2}+\lambda^{2}\Lambda^{2}\\\ \end{bmatrix}.$

Since $\|B\|_{2}$ is the square root of the maximum eigenvalue of $BB^{T}$,
again by Gershgorin theorem,

$\|B\|_{2}^{2}\leq\max_{l}\\{(1+\lambda\nu_{l})^{2}+\lambda^{2}\nu_{l}^{2}+2\lambda|\nu_{l}||1+\lambda\nu_{l}|\leq
1$

by using $\lambda\leq 1/4$. Hence,

$\displaystyle\sigma_{\max}(L^{-1})\leq N_{t}=1/\tau.$

Thus, $\sigma_{\min}(L)\geq\tau$.

To conclude, the matrix $H$ has order $N=2N_{t}(N_{x}-1)^{2}N_{x}$ and
sparsity number $s\leq 7$. From the above eigenvalue estimations, one knows
that the condition number

$\kappa\leq\frac{2}{\tau}\sim N_{t}\sim N_{x}^{3}$

using the condition $\lambda=\tau/h^{3}$. ∎

####  A.1.3 Higher space dimension

Now we consider the case of $x\in[0,1]^{d}$, with homogeneous Dirichlet
boundary condition in $x$, Other boundary conditions can be similarly treated
with possible different $\kappa$. Let $p\in[0,1]^{L}$, with zero Dirichlet
boundary condition at $p_{i}=1$ for all $i$. We will use a uniform mesh size
$h$ with the same number of mesh points $N_{x}$ in all space dimensions of $x$
and $p$.

Let $D=dL$. The matrix $L_{h}$ of the $D$-dimensional problem will be replaced
by

$L_{h,D}=\underbrace{L_{h}\otimes I\otimes\cdots\otimes
I}_{D~{}\mbox{matrices}}+I\otimes L_{h}\otimes\cdots\otimes I+\cdots+I\otimes
I\otimes\cdots\otimes L_{h},$

and everything else remains the same. We can still carry out the similar proof
in Lemma 25, except that the eigenvalues $\nu_{l}$ of $L_{h}$ will replaced by
the corresponding results of the new matrix. By the properties of tensor
products, the eigenvalues of $L_{h,D}$ can be represented by the sum of the
eigenvalues of $L_{h}$ as

$\nu^{D}_{i_{1},i_{2},\cdots,i_{D}}=\nu_{i_{1}}+\nu_{i_{2}}+\cdots+\nu_{i_{D}},\quad\nu_{l}=-4\sin^{2}(l\pi
h/2)$

for $i_{l}=1,\cdots,N_{x}-1,l=1,\cdots,D$. Like was done in (74), one gets

$\displaystyle\mu^{D}_{\max}\leq\max_{l}\\{1+(1+\lambda\nu^{D}_{i_{1},i_{2},\cdots,i_{D}})^{2}+2(|1+\lambda\nu^{D}_{i_{1},i_{2},\cdots,i_{D}}|+\lambda|\nu^{D}_{i_{1},i_{2},\cdots,i_{D}}|)|+2\lambda|\nu^{D}_{i_{1},i_{2},\cdots,i_{D}}||1+\lambda\nu^{D}_{i_{1},i_{2},\cdots,i_{D}}|$
$\displaystyle\qquad+\lambda^{2}(\nu^{D}_{i_{1},i_{2},\cdots,i_{D}})^{2}\\}\leq
4$ (76)

if

$\lambda=\tau/h^{3}\leq 1/(4D).$ (77)

Like what was done in Lemma 25,
$\nu^{D}_{max}(L^{-1})=O(N_{t})=O(DN_{x}^{3})$. Hence
$\kappa=O(DN_{x}^{3})=O(LdN_{x}^{3})$.

One can easily check that the sparsity $s=O(Ld)$.

###  A.2 Linear Boltzmann equation

We now consider the numerical approximation of the phase space linear
Boltzmann equation with isotropic scattering in Eq. (29). For simplicity we
only consider homogeneous boundary condition in $x$ and the case of
$b_{i}(x)>0$ so only the homogeneous boundary condition on the right boundary
of $p$ is needed. We also just consider $p_{i}\in[0,1]$.

We start with one-dimensional $x\in[0,1],v\in(-1,1)$ and $p\in[0,1]$. Then Eq.
(29) is

$\begin{cases}\partial_{t}F+v\partial_{x}F=-\left[\frac{1}{2}\int_{0}^{1}\partial_{p}F\,dv-\partial_{p}F\right]\,,\\\
F(0,x,v,z,p)=a(z)e^{-pa(z)}f_{0}(x,v,z).\end{cases}\,\\\ $ (78)

We use the upwind scheme for $x$ and $p$, and discrete-ordinate method for
$v$, with $(\omega_{m},v_{m}),m=\pm 1,\cdots,\pm M$ the pair of quadrature
points and weights to approximate the integration in $v$. Here we first define
quadrature points on $(0,1)$ and then use symmetry to define the quadrature
points over $(-1,0)$, therefore $v_{-m}=-v_{m},\omega_{-m}=\omega_{m}$ and $0$
is not a quadrature point. The weights $\omega_{m}$ satisfy

$\frac{1}{2}\sum_{|m=1|}^{M}\omega_{m}=1\,.$ (79)

Let $F^{n}_{j,m,l,k}$ be the approximation of $F(t,x,v,z,p)$ at
$(t^{n},x_{j},v_{m},z_{l},p_{k})$. Then the discrete system of (78) is

$\displaystyle\begin{cases}&\frac{1}{\tau}(F^{n+1}_{j,m,l,k}-F^{n}_{j,m,l,k})+\frac{1}{2h}(v_{m}+|v_{m}|)(F^{n}_{j,m,l,k}-F^{n}_{j-1,m,l,k})+\frac{1}{2h}(v_{m}-|v_{m}|)(F^{n}_{j+1,m,l,k}-F^{n}_{j,m,l,k})\\\
&=-\frac{1}{h}\frac{1}{2}\sum_{|m^{\prime}|=1}^{M}\omega_{m^{\prime}}(F^{n}_{j,m^{\prime},l,k+1}-F^{n}_{j-1,m^{\prime},l,k})+\frac{1}{h}(F^{n}_{j,m,l,k+1}-F^{n}_{j-1,m,l,k})\\\
&\\\
&F^{0}_{j,m,l,k}=a(z_{l})e^{-p_{k}a(z_{l})}f_{0}(x_{j},v_{m},z_{l})\end{cases}\,$
(80)

Since the evolution equation does not depend on $z$, we will drop the index
$l$ from $F_{j,m,l,k}$ in the sequel as long as there is no confusion.

Introduce $f_{j,k}=(f_{j,-M,k},\cdots,f_{j,M,k})^{T}$. Define

$V^{-}=\begin{bmatrix}v_{-M}&&&&&\\\ &\ddots&&&&\\\ &&v_{-1}&&&\\\ &&&0&&\\\
&&&&\ddots&\\\ &&&&&0\end{bmatrix},\qquad V^{+}=\begin{bmatrix}0&&&&&\\\
&\ddots&&&&\\\ &&0&&&\\\ &&&v_{1}&\\\ &&&&\ddots&\\\
&&&&&v_{M}\end{bmatrix},\qquad I_{M}=\begin{bmatrix}1&\cdots&1\\\
\cdots&\cdots&\cdots\\\ 1&\cdots&1\end{bmatrix}_{2M\times 2M}.$

and $W=\frac{1}{2}\text{dial}(\omega_{-M},\cdots,\omega_{M})^{T}$. Then scheme
(80) can be written in vector form as

$f^{n+1}_{j,k}-f^{n}_{j,k}+\lambda{v^{-}}\,(f^{n}_{j+1,k}-f^{n}_{j,k})+\lambda{v^{+}}\,(f^{n}_{j,k}-f^{n}_{j-1,k})-\lambda
I_{M}W(f^{n}_{j,k+1}-f^{n}_{j,k})+\lambda(f^{n}_{j,k+1}-f^{n}_{j,k})=0,$ (81)

where $\lambda=\tau/h$. Denote $f^{n}=(f^{n}_{1},\cdots,f^{n}_{N_{x}-1})^{T}$,
and

$L_{h}=\begin{bmatrix}-I&I&&&\\\ &-I&I&&\\\ &&\ddots&\ddots&\\\
&&&-I&I\end{bmatrix}_{2M(N_{x}-1)\times
2M(N_{x}-1)}\qquad\mathbb{V}^{\pm}={\text{diag}}(V^{\pm},\cdots,V^{\pm}),\qquad\mathbb{W}={\text{diag}}(I_{M}W,\cdots,I_{M}W)$

where $I$ is the $2M\times 2M$ identity matrix, and $\mathbb{I}$ is the
$2M\times(N_{x}-1)$ identity matrix. Then the above scheme can be written as

$f^{n+1}_{k}-f^{n}_{k}+\lambda\mathbb{V}^{-}L_{h}\,f^{n}_{k}+\lambda\mathbb{V}^{+}\,L_{h}^{T}f^{n}_{k}-\lambda\mathbb{W}(f^{n}_{k+1}-f^{n}_{k})+\lambda(f^{n}_{k+1}-f^{n}_{k})=\lambda
b(k)$ (82)

with $b(k)=(v^{+}f_{0,k},\cdots,-v^{-}f_{N_{x},k})$.

Let $u^{n}=(u_{0}^{n},\cdots,u_{N_{x}-1}^{n})^{T}$. Using the boundary
condition $u^{n}_{N_{x}}=0$, then (70) can be further written as

$-Bu^{n}+u^{n+1}=f^{n+1}=:\lambda(b_{0},\cdots,b_{N_{x}-1})^{T},$ (83)

where

$B=\begin{bmatrix}B_{1}&B_{2}&&&\\\ &B_{1}&B_{2}&&\\\ &&\ddots&\ddots&\\\
&&&\ddots&B_{2}\\\ &&&&B_{1}\\\ \end{bmatrix},$

with
$B_{1}=I-\lambda\mathbb{V}^{-}L_{h}-\lambda\mathbb{V}^{+}L_{h}^{T}-\lambda\mathbb{W}+\lambda
I$ and $B_{2}=\lambda\mathbb{W}-\lambda\mathbb{I}$. By introducing the
notation $U=[u^{1},\cdots,u^{N_{t}}]^{T}$, one obtains the following linear
system

$LU=F,$ (84)

where

$L=\begin{bmatrix}I&&&\\\ -B&I&&\\\ &\ddots&\ddots&\\\ &&-B&I\\\
\end{bmatrix},\qquad F=\begin{bmatrix}f^{1}+Bu^{0}\\\ f^{2}\\\ \vdots\\\
f^{N_{t}}\\\ \end{bmatrix}.$

Note that $\mathcal{M}$ matrix is denoted

$H=\begin{bmatrix}O&L\\\ L^{T}&O\end{bmatrix}.$

Let $\lambda$ be an eigenvalue of $H$, i.e.,

$\mbox{det}(\lambda I-H)=\mbox{det}(\lambda^{2}I-LL^{T})=:\mbox{det}(\mu
I-LL^{T}),$

where $\mu=\lambda^{2}$ is the eigenvalue of $LL^{T}$, or $\mu^{1/2}$ is the
singular value of $L$. A direct calculation gives

$LL^{T}=\begin{bmatrix}I&-B&&\\\ -B&I+B^{2}&\ddots&\\\ &\ddots&\ddots&-B\\\
&&-B&I+B^{2}\\\ \end{bmatrix}.$

We again establish the upper and lower bounds of the eigenvalues of $LL^{T}$
by using the Gershgorin circle theorem.

######  Lemma 26.

The minimum and maximum eigenvalues of $H$ satisfy

$|\lambda|_{\min}=O(\tau),\qquad|\lambda|_{\max}=O(1).$

###### Proof.

Note

$B^{2}=\begin{bmatrix}B_{1}^{2}&B_{1}B_{2}+B_{2}B_{1}&B_{2}^{2}&&&\\\
&B_{1}^{2}&B_{1}B_{2}+B_{2}B_{1}&B_{2}^{2}&&\\\ &&\ddots&\ddots&\ddots&\\\
&&&B_{1}^{2}&B_{1}B_{2}+B_{2}B_{1}&B_{2}^{2}\\\
&&&&B_{1}^{2}&B_{1}B_{2}+B_{2}B_{1}\\\ &&&&&B_{1}^{2}\end{bmatrix}.\\\ $

If $\lambda\leq 1$, and also due to the condition (79), one has $W^{2}=W$ and
the sums of the each rows of $W$ is bounded by $1$, it is rather easy to see
that Gershgorin’s theorm applied to $LL^{T}$ implies that

$\mu_{max}\leq C$

where $C$ is a constant independent of $M,N_{x}$ and $N_{t}$.

To estimate the maximum eigenvalues of $LL^{T}$, by definition, the spectral
radius of $L$

$\sigma_{\min}(L)=\frac{1}{\sigma_{\max}(L^{-1})}.$

By simple algebra, one has

$L^{-1}=\begin{bmatrix}I&&&&\\\ B&I&&&\\\ B^{2}&B&I&&\\\
\vdots&\ddots&\ddots&\ddots&\\\ B^{N_{t}-1}&\cdots&&B&I\\\
\end{bmatrix}=\begin{bmatrix}I&&&&\\\ &I&&&\\\ &&\ddots&&\\\ &&&\ddots&\\\
&&&&I\\\ \end{bmatrix}+\begin{bmatrix}&&&&\\\ B&&&&\\\ &B&&&\\\ &&\ddots&&\\\
&&&B\end{bmatrix}+\cdots,$ (85)

which gives

$\displaystyle\sigma_{\max}(L^{-1})$
$\displaystyle=\|L^{-1}\|_{2}\leq\|I\|_{2}+\|B\|_{2}+\cdots+\|B^{N_{t}-1}\|_{2}$
$\displaystyle\leq\|I\|_{2}+\|B\|_{2}+\|B\|_{2}^{2}+\cdots+\|B\|_{2}^{N_{t}-1}$

Note

$BB^{T}=\begin{bmatrix}B_{1}B_{1}^{T}+B_{2}B_{2}^{T}&B_{2}B_{1}^{T}&&&\\\
B_{1}B_{2}^{T}&B_{1}B_{1}^{T}+B_{2}B_{2}^{T}&B_{2}B_{1}^{T}&B_{2}B_{1}^{T}&&\\\
&\ddots&\ddots&\ddots&\\\
&&B_{1}B_{2}^{T}&B_{1}B_{1}^{T}+B_{2}B_{2}^{T}&B_{2}B_{1}^{T}&\\\
&&&B_{1}B_{2}^{T}&B_{1}B_{1}^{T}+B_{2}B_{2}^{T}&\\\ &&&&B_{2}B_{1}^{T}\\\
\end{bmatrix},$

Since $\|B\|_{2}$ is the square root of the maximum eigenvalue of $BB^{T}$,
again by Gershgorin’s theorem, under the condition that $\lambda\leq 1$, one
has

$\|B\|_{2}^{2}\leq C,$

where $C$ is independent of $M,N_{x}$ and $N_{t}$. Hence,

$\displaystyle\sigma_{\max}(L^{-1})\leq CN_{t}=O(1/\tau)$

Thus, $\sigma_{\min}(L)=O(\tau)=O(1/N_{x})$. The second equality is due to the
CFL condition $\lambda=\tau/h=O(1)$. ∎

The above estimates of eigenvalues lead to $k=O(N_{x})$. The sparsity of
number $s=O(M)$ due to the integral collision term.

For $x,v\in\mathbb{R}^{d}$ and $w\in\mathbb{R}^{L}$, similar to the multi-
dimensional case for the linear heat equation, it is easy to get

$\kappa=O((d+L)N_{x}),\quad s=O(LM^{d}).$

We omit the details of the proof.

###  A.3 Linear advection equation

For the advection equation in phase space (48), we use the forward Euler
method in time, center difference in $p$ and upwind discretization in $x$.
This equation contains third order derivatives, so the computational costs
will be similar to those of the phase space heat equation. The difference here
is that we have the boundary condition at $p_{i}=0$.

We first discuss the implementation of the boundary condition in (48). Since
it is a convection equation, a natural discretization is the upwind
discretization:

$\frac{u^{n}_{j1l}-u^{n}_{j0l}}{\tau}-\frac{u^{n}_{j0,l+1}-u^{n}_{j0l}}{h}=0,\quad
n=1,\cdots,N_{t},\quad j=1,\cdots,N_{x}-1;\quad l=0,\cdots,N_{x}-1$ (86)

or

$(1-\tau/h)u^{n}_{j0l}+u^{n}_{j0,l+1}=u^{n}_{j1l}$ (87)

This involves the coupling of two grids in $q$, which will make the the matrix
structure that defines the QLSP to be solved more complicated. Thus computing
the condition number of the matrix to be inverted becomes difficult. To
simplify the analysis, we first choose $\tau=h$ in (87). This corresponds to
the method of characteristics for the boundary condition. Then (87) becomes

$u^{n}_{j0,l+1}=u^{n}_{j1l}$ (88)

To avoid the coupling between different grids in $q$ we further use

$u^{n}_{j0,l+1}=u^{n}_{j0,l}+O(h)$

Now the boundary condition becomes, after ignoring the $O(h)$ term,

$u^{n}_{j0,l+1}=u^{n}_{j0,l}$

which is just the numerical implementation of the Neumann boundary condition.
At $p_{i}=1$ the Dirichet boundary condition is used. Since the largest
discrete eigenvalues of the discrete Laplacian with mixed Neumann-Dirichlet
boundary condition is the same order of that with the Dirichlet boundary
condition [36], a similar analysis as in section A.1 gives

$\kappa=O(LdN_{x}^{3}),\qquad s=O(Ld)\,.$

###  A.4 The Schrödinger equation

For the Schrödinger equation in phase space (62), the equations have second
order derivatives in both $x$ and $s$, but the highest order derivative
remains second order, so the equation is similar to that of the original
Schrödinger equation (56) except the total dimension increases from $d$ to
$d+L$, and one has the boundary condition at $w_{i}=0$. We use forward Euler
method in time and center difference in $x$ and $w$. We use the same
apporximation of the boundary condition as the case of the convection
equation, then a similar analysis as in the case of heat equation [37] gives:

$\kappa=O((d+L)N_{x}^{2}),\,\,s=O(d+L).$

##  Appendix B Numerical discretisations and classical cost

In this section we always assume that the spatial derivatives are approximated
by $r^{\text{th}}$ order method and explicit time discretization is used
(although in the previous section our analysis was done for given $r$’s.

###  B.1 Linear heat equation

For the original linear heat equation Eq. (5), the classical cost is
$\mathcal{C}_{can}=O(M\bar{N}^{d}\bar{N}_{t}d)$. To compute how $\bar{N}$ and
$\bar{N}_{t}$ depend on $d$ and $\epsilon$, we require the following error and
CFL conditions

$\displaystyle\begin{cases}Error:\qquad O(d(\Delta
x)^{r})=O(d/\bar{N}^{r})=O(\epsilon)\\\ CFL:\qquad O(1/\bar{N}_{t})=\Delta
t=O((\Delta x)^{2}/d)=O(1/(\bar{N}^{2}d))\end{cases}$ (89)

for $r^{\text{th}}$-order methods, where $c\geq 1$. Then Eq. (89) implies

$\displaystyle\bar{N}=O((d/\epsilon)^{1/r}),\qquad\bar{N}_{t}=O(d(d/\epsilon)^{2/r}).$
(90)

Therefore

$\displaystyle\mathcal{C}_{can}=O(M\bar{N}^{d}\bar{N}_{t}d)=O(Md^{2}(d/\epsilon)^{(d+2)/r}).$
(91)

For the modified linear heat equation in Eq. (8), the classical cost is
$\mathcal{C}_{mod}=O(LN^{d}N_{p}^{L}N_{t}d)$. We define $N=N_{p}$. To compute
how $N=N_{p}$ and $N_{t}$ depend on $d$ and $\epsilon$, we require the
following error and CFL conditions

$\displaystyle\begin{cases}Error:\qquad O(Ld(\Delta
x)^{r})=O(Ld/N^{r})=O(\epsilon)\\\ CFL:\qquad O(1/N_{t})=\Delta t=O((\Delta
x)^{3}/(L(d+1)))=O(1/(LdN^{3}))\end{cases}$ (92)

for $r^{\text{th}}$-order methods. Eq. (92) then implies

$\displaystyle N=O((Ld/\epsilon)^{1/r}),\qquad
N_{t}=O((Ld)^{1+3/r}/\epsilon^{3/r}).$ (93)

Therefore

$\displaystyle\mathcal{C}_{mod}=O(LN^{d}N_{p}^{L}N_{t}d)=O((Ld)^{2+(d+L+3)/r}/\epsilon^{(d+L+3)/r}).$
(94)

We want to compare quantum costs to the minimal classical cost. Here, the
minimal classical cost $\mathcal{C}=\mathcal{C}_{can}<\mathcal{C}_{mod}$ when
$M<M_{heat}\equiv O(L^{2+(d+L+3)/r}(d/\epsilon)^{(L+1)/3})$. If $M$ is very
large and $d,L$ are relatively small so we can have the condition
$M>M_{heat}$, then $\mathcal{C}=\mathcal{C}_{mod}<\mathcal{C}_{can}$.

###  B.2 Linear Boltzmann equation

For the original linear Boltzmann equation Eq. (20), the classical cost is
$\mathcal{C}_{can}=O(M\bar{N}^{d}\bar{N}_{v}^{d}\bar{N}_{t})$. To compute how
$\bar{N}$ and $\bar{N}_{t}$ depend on $d$ and $\epsilon$, we require the
following error and CFL conditions

$\displaystyle\begin{cases}Error:\qquad O(d\Delta x^{r}+d(\Delta
v)^{r})=O(d/\bar{N}^{r}+d/\bar{N}_{v}^{r})=O(\epsilon)\\\ CFL:\qquad
O(1/\bar{N}_{t})=\Delta t=O(\Delta x/d)=O(1/(\bar{N}d))\end{cases}$ (95)

We define $\bar{N}=\bar{N}_{p}$. Then Eq. (95) implies

$\displaystyle\bar{N}=O((d/\epsilon)^{1/r}),\qquad\bar{N}_{t}=O(d(d/\epsilon)^{1/r}).$
(96)

Therefore

$\displaystyle\mathcal{C}_{can}=O(M\bar{N}^{d}\bar{N}_{v}^{d}\bar{N}_{t})=O(Md(d/\epsilon)^{(2d+1)/r}).$
(97)

For the modified linear Boltzmann equation in Eq. (29), the classical cost is
$\mathcal{C}_{mod}=O(N^{d}N_{p}^{L}N_{v}^{d}N_{t})$. We define
$N=N_{v}=N_{p}$. To compute how $N$ and $N_{t}$ depend on $d$ and $\epsilon$,
we require the following error and CFL conditions

$\displaystyle\begin{cases}Error:\qquad O(d(\Delta x)^{r}+L(d(\Delta
v)^{r}+(\Delta w)^{r}))=O((d+L+Ld)/N^{r})=O(\epsilon)\\\ CFL:\qquad
O(1/N_{t})=\Delta t=\min(\Delta x/d,\Delta w/L)=\min(1/L,1/d)/N\end{cases}$
(98)

Then Eq. (98) implies

$\displaystyle N=O(((d+L+Ld)/\epsilon)^{1/r})=O((Ld/\epsilon)^{1/r}),\qquad
N_{t}=O(\max(L,d)(Ld/\epsilon)^{1/r})$ (99)

Therefore

$\displaystyle\mathcal{C}_{mod}=O(N^{d}N_{p}^{L}N_{v}^{d}N_{t})=O((Ld/\epsilon)^{(2d+L+1)/r}\max(L,d)).$
(100)

We want to compare the quantum cost to the minimal classical cost. Here the
minimal classical cost is $\mathcal{C}=\mathcal{C}_{can}<\mathcal{C}_{mod}$
when $M<M_{Boltz}\equiv O(L^{(2d+L+1)/r}\max(L,d)(d/\epsilon)^{L/r}/d)$. If
$M>M_{Boltz}$, then $\mathcal{C}=\mathcal{C}_{mod}<\mathcal{C}_{can}$.

###  B.3 Linear advection equation

For the original linear heat equation Eq. (41), the classical cost is
$\mathcal{C}_{can}=O(M\bar{N}^{d}\bar{N}_{t}d)$. To compute how $\bar{N}$ and
$\bar{N}_{t}$ depend on $d$ and $\epsilon$, we require the following error and
CFL conditions

$\displaystyle\begin{cases}Error:\qquad O(d(\Delta
x)^{r})=O(d/\bar{N}^{r})=O(\epsilon)\\\ CFL:\qquad O(1/\bar{N}_{t})=\Delta
t=O(\Delta x/d)=O(1/(\bar{N}d))\end{cases}$ (101)

Eq. (101) implies

$\displaystyle\bar{N}=O((d/\epsilon)^{1/r}),\qquad\bar{N}_{t}=O(d(d/\epsilon)^{1/r}).$
(102)

Therefore

$\displaystyle\mathcal{C}_{can}=O(M\bar{N}^{d}\bar{N}_{t}d)=O(Md^{2}(d/\epsilon)^{(d+1)/r}).$
(103)

For the modified linear advection equation in Eq. (48), the classical cost is
$\mathcal{C}_{mod}=O(LN^{d}N_{p}^{L}N_{q}^{L}N_{t}d)$. We define
$N=N_{p}=N_{q}$. To compute how $N=N$ and $N_{t}$ depend on $d$ and
$\epsilon$, we require the following error and CFL conditions

$\displaystyle\begin{cases}Error:\qquad O(Ld(\Delta
x)^{r})=O(Ld/N^{r})=O(\epsilon)\\\ CFL:\qquad O(1/N_{t})=\Delta t=O((\Delta
x)^{3}/(Ld))=O(1/(LN^{3}d))\end{cases}$ (104)

Eq. (104) implies

$\displaystyle N=O((Ld/\epsilon)^{1/r}),\qquad
N_{t}=O(Ld(Ld/\epsilon)^{3/r}).$ (105)

Therefore

$\displaystyle\mathcal{C}_{mod}=O(LN^{d}N_{p}^{L}N_{q}^{L}N_{t}d)=O((Ld)^{2}(Ld/\epsilon)^{(d+2L+3)/r}).$
(106)

We want to compare quantum costs to the minimal classical cost. Here, the
minimal classical cost $\mathcal{C}=\mathcal{C}_{can}<\mathcal{C}_{mod}$ when
$M<M_{adv}\equiv O(L^{(d+2L+3)/r+2}(d/\epsilon)^{(2L+2)/r})$. If $M$ is very
large and $d,L$ are relatively small so we can have the condition $M>M_{adv}$,
then $\mathcal{C}=\mathcal{C}_{mod}<\mathcal{C}_{can}$.

###  B.4 The Schrödinger equation

For the original linear heat equation Eq. (56), the classical cost is
$\mathcal{C}_{can}=O(M\bar{N}^{d}\bar{N}_{t}d)$. To compute how $\bar{N}$ and
$\bar{N}_{t}$ depend on $d$ and $\epsilon$, we require the following error and
CFL conditions

$\displaystyle\begin{cases}Error:\qquad O(d(\Delta
x)^{r})=O(d/\bar{N}^{r})=O(\epsilon)\\\ CFL:\qquad O(1/\bar{N}_{t})=\Delta
t=O((\Delta x)^{2}/d)=O(1/(N^{2}d))\end{cases}$ (107)

Eq. (107) implies

$\displaystyle\bar{N}=O((d/\epsilon)^{1/r}),\qquad\bar{N}_{t}=O(d(d/\epsilon)^{2/r}).$
(108)

Therefore

$\displaystyle\mathcal{C}_{can}=O(M\bar{N}^{d}\bar{N}_{t}d)=O(Md^{2}(d/\epsilon)^{(d+2)/r})$
(109)

For the modified Schrödinger equation in Eq. (62), the classical cost is
$\mathcal{C}_{mod}=O(N^{d}N_{p}^{L}N_{q}^{L}N_{t}(d+L))$. We define
$N=N_{p}=N_{q}$. To compute how $N$ and $N_{t}$ depend on $d$ and $\epsilon$,
we require the following error and CFL conditions

$\displaystyle\begin{cases}Error:\qquad O(d(\Delta x)^{r}+L(\Delta
w)^{r})=O((d+L)/N^{r})=O(\epsilon)\\\ CFL:\qquad O(1/N_{t})=\Delta t=(\Delta
x)^{2}/(d+L)=1/(N^{2}(d+L))\end{cases}$ (110)

Eq. (110) implies

$\displaystyle N=O(((d+L)/\epsilon)^{1/r}),\qquad
N_{t}=O((d+L)((d+L)/\epsilon)^{2/r}).$ (111)

Therefore

$\displaystyle\mathcal{C}_{mod}=O(N^{d}N^{2L}N_{t}(d+L))=O((d+L)^{2}((d+L)/\epsilon)^{(d+2L+2)/r}).$
(112)

We want to compare quantum costs to the minimal classical cost. Here, the
minimal classical cost $\mathcal{C}=\mathcal{C}_{can}<\mathcal{C}_{mod}$ when
$M<M_{Schr}\equiv O((d+L)^{(d+2L+2)/r+2}/(d^{(d+2)/r+2}\epsilon^{2L/r}))$. If
$M$ is very large and $d,L$ are relatively small so we can have the condition
$M>M_{Schr}$, then $\mathcal{C}=\mathcal{C}_{mod}<\mathcal{C}_{can}$.

##  Appendix C Quantum cost of solving original PDEs

To solve for ensemble averaged observables of the original PDEs (not in phase
space) with $M$ initial conditions, each PDE needs to be solved $M$ times
separately. This means for each initial condition, a different $\mathcal{M}$,
hence a different sparse-access is required. There may be cases where this can
be either costly to obtain or is not available, hence solving the phase space
equations is preferred.

Suppose we assume that the $M$ different sparse access to $\mathcal{M}$ are
available and can be accessed without extra cost. The quantum cost to compute
the observables for the linear heat, linear Boltzmann and linear advection
equations can be obtained from Lemma 3. Here the sparsity and condition
nubmers of the corresponding $\mathcal{M}$ are found in Appendix A.

For the linear heat equation, the sparsity $s$ and conditions number $\kappa$
of $\mathcal{M}$ are $O(d)$ and $O(N^{2})$ respectively, where
$N=O(d/\epsilon)^{1/r}$ for $r^{\text{th}}$-order methods. Lemma 3 can be used
directly to arrive at
$\mathcal{Q}_{orig}=\tilde{\mathcal{O}}(Ms\kappa^{3}\mathcal{N}_{y}^{2}/\epsilon)$
to compute the observable crorresponding to $M$ initial conditions to
precision $\epsilon$, where the normalisation $\mathcal{N}_{y}$ is the
normalisation corresponding to the initial condition of the original PDE,
which is a factor of order $O(N_{p}^{L})$ less than the normalisation of the
initial state for the phase space equation. This means that
$\mathcal{N}^{2}_{y}$ term is comparable in size to the normalisation
$n^{2}_{W_{0}}$ corresponding to the phase space equation. Then comparing
$\mathcal{Q}_{orig}$ with $\mathcal{Q}$ from Theorem 8, we have
$\mathcal{Q}<\mathcal{Q}_{orig}$ when
$M>\tilde{\mathcal{O}}(L^{4+9/r}d^{3}(d/\epsilon)^{3/r})$.

The analysis is the same for the linear Boltzmann equation ($s=O(N_{v}^{d})$
and $\kappa=O(dN)$) and the linear advection equation ($s=O(d)$ and
$\kappa=O(dN)$) is the same, with $N=N_{v}=O(d/\epsilon)^{1/r}$. For the
Schrödinger equation, we can still use Lemma 3, but potentially better scaling
is possible when one uses either Lemma 4 or quantum simulation methods. For
consistency in also using a matrix inversion method, we use Lemma 4 where
$s=O(d)$ and $\kappa=O(N^{2})$ for $N=(d/\epsilon)^{1/r}$. Quantum simulation
methods will differ only by small polynomial factors, which will not
significantly affect the final comparison to $\mathcal{Q}$. We omit details of
the proof since it is straightforward.

## References

  * Smith [2013] R. C. Smith,  _Uncertainty quantification: theory, implementation, and applications_ , Vol. 12 (Siam, 2013).
  * Ghanem  _et al._ [2017] R. Ghanem, D. Higdon, H. Owhadi, _et al._ ,  _Handbook of uncertainty quantification_ , Vol. 6 (Springer, 2017).
  * Xiu [2010] D. Xiu, Numerical methods for stochastic computations, in  _Numerical Methods for Stochastic Computations_ (Princeton university press, 2010).
  * Gunzburger  _et al._ [2014] M. D. Gunzburger, C. G. Webster, and G. Zhang, Stochastic finite element methods for partial differential equations with random input data, Acta Numerica 23, 521 (2014).
  * Ghanem and Spanos [2003] R. G. Ghanem and P. D. Spanos,  _Stochastic finite elements: a spectral approach_ (Courier Corporation, 2003).
  * Bungartz and Griebel [2004] H.-J. Bungartz and M. Griebel, Sparse grids, Acta numerica 13, 147 (2004).
  * Cohen and DeVore [2015] A. Cohen and R. DeVore, Approximation of high-dimensional parametric pdes, Acta Numerica 24, 1 (2015).
  * Harrow  _et al._ [2009] A. W. Harrow, A. Hassidim, and S. Lloyd, Quantum algorithm for linear systems of equations, Phys. Rev. Lett. 103, 150502, 4 pp. (2009).
  * Childs  _et al._ [2017] A. M. Childs, R. Kothari, and R. D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision, SIAM J. Comput. 46, 1920 (2017).
  * Gilyén  _et al._ [2019] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics, in  _Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_ (2019) pp. 193–204.
  * Berry [2014] D. W. Berry, High-order quantum algorithm for solving linear differential equations, J. Phys. A: Math. Theor. 47, 105301, 17 pp. (2014).
  * Joseph [2020] I. Joseph, Koopman–von neumann approach to quantum simulation of nonlinear classical dynamics, Physical Review Research 2, 043102 (2020).
  * Dodin and Startsev [2021] I. Y. Dodin and E. A. Startsev, On applications of quantum computing to plasma simulations, Physics of Plasmas 28, 092101 (2021).
  * Lloyd  _et al._ [2020] S. Lloyd, G. De Palma, C. Gokler, B. Kiani, Z.-W. Liu, M. Marvian, F. Tennie, and T. Palmer, Quantum algorithm for nonlinear differential equations, arXiv preprint arXiv:2011.06571 (2020).
  * Childs  _et al._ [2021] A. M. Childs, J. P. Liu, and A. Ostrander, High-precision quantum algorithms for partial differential equations, Quantum 5, 574 (2021).
  * Liu  _et al._ [2021] J. Liu, H. Ø. Kolden, H. K. Krovi, N. F. Loureiro, K. Trivisa, and A. M. Childs, Efficient quantum algorithm for dissipative nonlinear differential equations, Proc. Natl. Acad. Sci. U. S. A. 118 (2021).
  * Jin  _et al._ [2022a] S. Jin, N. Liu, and Y. Yu, Time complexity analysis of quantum difference methods for linear high dimensional and multiscale partial differential equations, arXiv preprint arXiv:2202.04537 (2022a).
  * Jin and Liu [2022] S. Jin and N. Liu, Quantum algorithms for computing observables of nonlinear partial differential equations, arXiv preprint arXiv:2202.07834 (2022).
  * Bernard  _et al._ [2010] E. Bernard, F. Golse, and F. Salvarani, Homogenization of transport problems and semigroups, Mathematical methods in the applied sciences 33, 1228 (2010).
  * Ambainis [2012] A. Ambainis, Variable time amplitude amplification and quantum algorithms for linear algebra problems, LIPIcs. Leibniz Int. Proc. Inform. 14, 636 (2012).
  * Lin and Tong [2020] L. Lin and Y. Tong, Optimal polynomial based quantum eigenstate filtering with application to solving quantum linear systems, Quantum 4, 361 (2020).
  * Costa  _et al._ [2021] P. Costa, D. An, Y. R. Sanders, Y. Su, R. Babbush, and D. W. Berry, Optimal scaling quantum linear systems solver via discrete adiabatic theorem, arXiv preprint arXiv:2111.08152 (2021).
  * Low and Chuang [2019] G. H. Low and I. L. Chuang, Hamiltonian simulation by qubitization, Quantum 3, 163 (2019).
  * Alase  _et al._ [2021] A. Alase, R. R. Nerem, M. Bagherimehrab, P. Høyer, and B. C. Sanders, Tight bound for estimating expectation values from a system of linear equations, arXiv preprint arXiv:2111.10485 (2021).
  * Rall [2020] P. Rall, Quantum algorithms for estimating physical quantities using block encodings, Physical Review A 102, 022408 (2020).
  * Chakraborty  _et al._ [2019] S. Chakraborty, A. Gilyén, and S. Jeffery, The power of block-encoded matrix powers: Improved regression techniques via faster hamiltonian simulation, Leibniz international proceedings in informatics 132 (2019).
  * Linden  _et al._ [2020] N. Linden, A. Montanaro, and C. Shao, Quantum vs. classical algorithms for solving the heat equation, arXiv preprint arXiv:2004.06516 (2020).
  * Gleinig and Hoefler [2021] N. Gleinig and T. Hoefler, An efficient algorithm for sparse quantum state preparation, in  _2021 58th ACM/IEEE Design Automation Conference (DAC)_ (IEEE, 2021) pp. 433–438.
  * Zhang  _et al._ [2022] X.-M. Zhang, T. Li, and X. Yuan, Quantum state preparation with optimal circuit depth: Implementations and applications, arXiv preprint arXiv:2201.11495 (2022).
  * Zhang [2001] D. Zhang,  _Stochastic methods for flow in porous media: coping with uncertainties_ (Elsevier, 2001).
  * Jin [2009] S. Jin, Numerical methods for hyperbolic systems with singular coefficients: well-balanced scheme, Hamiltonian preservation, and beyond, in  _Hyperbolic problems: theory, numerics and applications_, Proc. Sympos. Appl. Math., Vol. 67 (Amer. Math. Soc., Providence, RI, 2009) pp. 93–104.
  * Lewis and Miller [1984] E. E. Lewis and W. F. Miller,  _Computational methods of neutron transport_ (John Wiley and Sons, Inc., New York, NY, 1984).
  * Fichtl [2009] E. D. Fichtl,  _Stochastic methods for uncertainty quantification in radiation transport_ , Ph.D. thesis, The University of New Mexico (2009).
  * Mishra  _et al._ [2016] S. Mishra, C. Schwab, and J. Šukys, Multi-level monte carlo finite volume methods for uncertainty quantification of acoustic wave propagation in random heterogeneous layered medium, Journal of Computational Physics 312, 192 (2016).
  * Ryzhik  _et al._ [1996] L. Ryzhik, G. Papanicolaou, and J. B. Keller, Transport equations for elastic and other waves in random media, Wave motion 24, 327 (1996).
  * Kuo and Shieh [2022] Y.-C. Kuo and S.-F. Shieh, Boundary effects on eigen-problems of discrete laplacian in lattices, Taiwanese Journal of Mathematics 1, 1 (2022).
  * Jin  _et al._ [2022b] S. Jin, N. Liu, and Y. Yu, Time complexity analysis of quantum difference methods for linear high dimensional and multiscale partial differential equations, arXiv preprint arXiv:2202.04537 (2022b).