# Quantum algorithm for stochastic optimal stopping problems with applications
in finance

João F. Doriguello<EMAIL_ADDRESS>Centre for Quantum Technologies, National
University of Singapore, Singapore  Alessandro Luongo<EMAIL_ADDRESS>Centre
for Quantum Technologies, National University of Singapore, Singapore  Jinge
Bao<EMAIL_ADDRESS>Centre for Quantum Technologies, National University of
Singapore, Singapore
Patrick Rebentrost<EMAIL_ADDRESS>Centre for Quantum Technologies, National
University of Singapore, Singapore  Miklos Santha<EMAIL_ADDRESS>Centre for
Quantum Technologies, National University of Singapore, Singapore

###### Abstract

The famous least squares Monte Carlo (LSM) algorithm combines linear least
square regression with Monte Carlo simulation to approximately solve problems
in stochastic optimal stopping theory. In this work, we propose a quantum LSM
based on quantum access to a stochastic process, on quantum circuits for
computing the optimal stopping times, and on quantum techniques for Monte
Carlo. For this algorithm, we elucidate the intricate interplay of function
approximation and quantum algorithms for Monte Carlo. Our algorithm achieves a
nearly quadratic speedup in the runtime compared to the LSM algorithm under
some mild assumptions. Specifically, our quantum algorithm can be applied to
American option pricing and we analyze a case study for the common situation
of Brownian motion and geometric Brownian motion processes.

##  1 Introduction

Within stochastic optimization, optimal stopping theory is a broad area of
applied mathematics that started in the 1940s and 1950s mainly with A. Wald
[101] and is concerned with the problem of deciding the best time to “stop” or
take an action in order to maximize an expected reward [88]. The time at which
the observations are terminated, called _stopping time_ , is a random variable
depending on the data observed so far in the process. A simple example of an
optimal stopping problem is the following: consider a game in which $100$
numbers are written on $100$ pieces of paper without restrictions on the
numbers, except that no number appears more than once. The pieces of paper are
shuffled faced down and you are asked to look at the numbers, without having
seen them before and _once_ at a time, and to stop when you think that you
have found the biggest number. It turns out that there is a stopping rule that
allows you to stop at the biggest number for $1/e$ fraction of the inputs.

Since its conception, optimal stopping theory collected problems from many
disparate areas under a unique umbrella [78], e.g. quickest detection [89],
sequential parameter estimation [71] and sequential hypothesis testing [27].
Probably the most famous optimal stopping problem is the one of option pricing
in finance, especially American options [63]. A central problem in the world
of finance is to assign a monetary value to a hitherto unvalued asset. In the
capital markets, there exists a large variety of financial assets which are
derivative to underlying assets such as stocks, bonds, or commodities. One of
the most well-known examples is the European call option which allows the
buyer to “lock in” a price for buying a stock at some future time (or
“exercise” time). A fair valuation of such an option was first discussed in
the seminal works of Black and Scholes [11] and Merton [68]. In the times
since, the methods proposed in these works have become standard practice in
the financial sector and have been extended and generalized for many financial
derivatives and market models.

American options allow the buyer to exercise the option at any point in time
between the time of purchase and a fixed final time. In contrast to an
European option, there are no known closed formulas for the price of an
American option with finite maturity date even in simple models like the
Black-Scholes-Merton one. Theoretically, an American option can be viewed as a
stochastic optimal stopping problem for the buyer and a super-martingale
hedging problem for the seller [34]. Practical algorithms have been developed
for the pricing of American options [24, 84, 52, 50], an important class being
least squares Monte Carlo (LSM) algorithms originally proposed independently
by Tsitsiklis and Van Roy [96] and by Longstaff and Schwartz [63].

Among the aforementioned classical algorithms for option pricing — and other
topics in finance — is the sub-field of quantum computing of designing quantum
algorithms in the context of financial problems [74, 12, 31], e.g. risk
management [103], financial greeks [92, 2], portfolio optimization [5, 26, 49,
82, 46, 1] and option pricing [66, 100]. A common tool in obtaining a quantum
advantage is amplitude estimation [14] and its generalizations for Monte Carlo
sampling [70, 45, 23, 22]. A few different works devised quantum algorithms
for derivative pricing based on quantum subroutines for Monte Carlo [81, 91,
18], e.g. European [80, 35, 79] and American/Bermudan [69] option pricing, and
option pricing in the local volatility model [55, 2] (of which the Black-
Scholes model is a subcase). Given its versatility and previous cases of
success, it is only natural to explore the applicability of quantum methods
for Monte Carlo to problems in optimal stopping theory. In this work we focus
on tailoring these methods to LSM algorithms.

###  1.1 Problem statement and past results

Optimal stopping theory is concerned with the problem of finding the best
moment to stop a process in order to maximize an expected reward. More
generally, assume a stochastic process (which corresponds to the market model
in financial applications) modelled with a Markov chain $(X_{t})_{t=0}^{T}$
with sample space $\Omega$ and state space $E\subseteq\mathbb{R}^{d}$, i.e.,
the random variables defined as $X_{t}:\Omega\to E$. Each element $X_{t}$ for
$t\in\\{0,\dots,T\\}$ gives rise to an image probability measure $\rho_{t}$ in
$E\subseteq\mathbb{R}^{d}$. Let $L^{2}(E,\rho_{t})=L^{2}(\rho_{t})$ be the set
of squared integrable functions with norm
$\|f\|_{L^{2}(\rho_{t})}=\sqrt{\mathbb{E}_{\rho_{t}}[|f(X_{t})|^{2}]}$. Assume
further a _payoff process_ : another stochastic process $(Z_{t})_{t=0}^{T}$
obtained from $(X_{t})_{t=0}^{T}$ as $Z_{t}=z_{t}(X_{t})$ for square-
integrable real functions $z_{t}\in L^{2}(E)$, for $t\in\\{0,\dots,T\\}$, over
the probability distribution of the Markov chain. The main problem is when to
stop the process and take the payoff so that the expected payoff is maximized.
A stopping time is a random variable $\tau:\Omega\to\\{0,1,\dots,T\\}$ such
that $\\{\omega\in\Omega|\tau(\omega)=t\\}$ does not depend on future
information (for $t^{\prime}>t$) about the stochastic process. The optimal
stopping problem then consists in finding a stopping time $\tau$ that
maximizes the expected value payoff.

######  Problem 1.1 (Optimal stopping problem).

For $\epsilon_{\rm final}>0$, approximate $\mathcal{U}_{0}$ to additive
accuracy $\epsilon_{\rm final}$ with high probability, where

$\mathcal{U}_{0}:=\operatorname*{sup}_{{\tau}~{}{\rm
stopping~{}time}}\mathbb{E}[Z_{\tau}].$

A more precise formulation is given in Problem 2.3. The quantity
$\mathcal{U}_{0}$ is related to a stochastic process called the _Snell
envelope_ [72, 34] and can be expressed as the solution of the dynamic
programming [21, 34]

$\displaystyle\begin{cases}\tau_{T}=T,\\\
\tau_{t}=t\mathbf{1}\\{Z_{t}\geq\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\\}+\tau_{t+1}\mathbf{1}\\{Z_{t}<\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\\},&0\leq
t\leq T-1.\end{cases}$

The stopping time $\tau_{0}$ can be shown to maximize $\mathbb{E}[Z_{\tau}]$
in Problem 1.1. The quantities $\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]$ are called
_continuation values_. In the past, many different approaches were developed
to tackle the dynamic programming above [56, 19, 24, 84, 87, 67, 77, 52, 8,
50, 98, 13, 36, 58, 17]. A famous approach is the least squares Monte Carlo
(LSM) by Longstaff and Schwartz [63]. The LSM algorithm combines linear least-
square regression with Monte Carlo simulation to approximate the hard-to-
compute continuation values arising in the dynamic programming.

####  1.1.1 Least squares Monte Carlo algorithms

The Longstaff-Schwartz algorithm combines linear least-square regression with
Monte Carlo simulation to approximating certain quantities known as
continuation values arising in a discrete-time optimal stopping problem. The
slightly simpler Tsitsiklis-Van Roy algorithm also tackles the problem by
parametric approximations based on stochastic approximation techniques [9].
More specifically, the LSM algorithm from Longstaff and Schwartz [63] solves
the dynamic programming problem by means of two approximations. The first one
is based on approximating the continuation values with a finite set of basis
functions. For each time step $t\in[T-1]$ we choose a set of $m$ expansion
functions $e_{t,k}(x)\in L^{2}(E)$, $k\in[m]$, that is rich enough to
approximate sufficiently well the continuation values
$\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]$ (in the $L^{2}$ norm) but is simple enough
to make the problem tractable. Often a few low-degree polynomials are
sufficient to obtain a good approximation. Then we can linearly decompose the
continuation value function as
$\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\approx\sum_{k=1}^{m}\alpha_{k}e_{t,k}(X_{t})$,
where $\alpha_{k}$ are projection coefficients given by
$\alpha_{t}=A_{t}^{-1}b_{t}$, where $A_{t}\in\mathbb{R}^{m\times m}$ is a
matrix with entries $(A_{t})_{k,l}=\mathbb{E}[e_{t,k}(X_{t})e_{t,l}(X_{t})]$
and $b_{t}\in\mathbb{R}^{m}$ is a vector with entries
$(b_{t})_{k}=\mathbb{E}[Z_{\tau_{t+1}}e_{t,k}(X_{t})]$. The key advantage of
the first approximation is that instead of Monte Carlo estimating the
continuation values $f(x)$ we can Monte Carlo estimate quantities involving
the expansion functions. In particular, we can use Monte Carlo estimation to
approximate the coefficients $\alpha_{t}$.

The second approximation is to solve the dynamic programming for the stopping
times by using Monte Carlo sampling. In other words, we sample $N$ independent
paths of the Markov Chain, compute the associated payoff values and solve the
dynamic programming for each of these paths (using the approximation to the
continuation values, which are easy to evaluate on the sample paths). The end
result is a set of optimal stopping times $\\{\tau^{(n)}\\}_{n=1}^{N}$, from
which the quantity $\sup_{\tau}\mathbb{E}[Z_{\tau}]$ is approximated as
$\frac{1}{N}\sum_{n=1}^{N}Z_{\tau^{(n)}}$. The LSM algorithm [63] takes as
input sampling access to the Markov chain and access to the payoff and the
expansion functions.

While Longstaff and Schwartz [63] offered some theoretical arguments for the
convergence of their algorithm, corroborated by many numerical experiments
with improved performances compared to other methods at the time, they did not
present any explicit convergence result. Tsitsiklis and Van Roy [96], on the
other hand, proved that the estimator in their algorithm converges to an
estimator of an ideal algorithm where sampling error from Monte Carlo
simulation is ignored. However, no convergence rate was provided.

The first theoretical analysis of the LSM algorithm was done by Clément et al.
[21] and focused on finite-dimensional linear approximation spaces. They
proved that the approximation error of the continuation values goes to $0$ as
the number of basis function goes to infinity. Moreover, they also provided an
almost sure convergence of the Monte Carlo procedure via a Central Limit
Theorem, but neither rate of convergence nor sample complexity was analyzed.
The generality of their results is also limited in the sense that nonlinear
approximation spaces were not studied and some assumptions were placed on the
underlying Markov process. The work of Stentoft [93] is placed in the
framework as [21] and claims to have obtained convergence rates in the two-
period case as both the number of basis functions and the number of simulated
paths go to infinity together under some assumptions. However, as noted by
[62], the proof of [93] contains some errors.

Glasserman and Yu [38] were the first to regenerate the sample paths at each
time step in order to simplify the error analysis. Moreover, they also
considered quasi-regression instead of regression. Explicit results were
obtained when taking polynomial functions as basis functions and considering
Brownian motion and geometric Brownian motion as underlying process. Their
analysis was later extended by Gerhold [37] for different Lévy models.

Egloff [32] unified the Longstaff-Schwartz and the Tsitsiklis-Van Roy
algorithms by proposing a new class of algorithms named dynamic look-ahead,
and managed to obtain major advances in the error analysis. By leveraging the
empirical risk minimization framework of Vapnik and Chervonenkis [99], Egloff
provided error estimates and sample complexity results by upper-bounding the
expected sample error by $O(\sqrt{\log(N)/N})$, where $N$ is the number of
sampled paths, for a variety of approximation spaces, including nonlinear
ones. His analysis requires that the payoff functions be bounded, and that the
approximation space be uniformly bounded, closed, and convex in some $L^{p}$
space, $2\leq p\leq\infty$, and possess finite Vapnik–Chervonenkis (VC)
dimension.

Following the results of Egloff [32], a series of works from Zanger slowly
dropped previous assumptions and generalized the error analysis. Starting in
[105], Zanger removed the convexity and closeness conditions on the
approximation space from Egloff [32], but obtained a worse expected $L^{2}$
sample error of $o(1/N^{1/4})$. The uniform boundedness of the payoff and
approximation functions were still presumed, together with finite VC
dimension. In [107], Zanger managed to drop the boundedness assumption on the
approximation functions, extending previous results to arbitrary sets of
$L^{2}$ functions with finite VC dimension, together with an improved
$O(\sqrt{\log(N)/N})$ expected $L^{2}$ sample error. However, like [38, 37],
regenerating new sets of independent sample paths at each time step was
needed. Such assumption was later dropped in [104], but at the expense of
again requiring uniformly bounded approximation functions. Finally, in [106]
Zanger proved an $O(\sqrt{\log(N)/N})$ expected $L^{2}$ sample error by using
a single set of independent sample paths. The approximation space was assumed
to be an arbitrary set of $L^{2}$ functions with finite VC dimension. Further,
the payoff functions were assumed to be bounded in $L^{p}$, $2<p<\infty$, but
not necessarily under the uniform norm. Moreover, this work also considered
the case in which the least-square regression step is solved only
approximately. The results from [106] are the most general ones that we are
aware of which apply to the original version of the LSM algorithm introduced
in [63].

Apart from the aforementioned works, Klimek and Pitera [57] investigated in
the same framework of Clément et al. [21] but considered the case when the
underlying process is non-Markovian and the payoff is path-dependent. Liu et
al. [62] revised the proof of Clément et al. [21] employing quasi-regression
and proved that the mean sample error asymptotically scales as $1/\sqrt{N}$.

####  1.1.2 Applications of LSM algorithm

Among the whole domain of optimal stopping problems, there are many that can
be approached directly with LSM, e.g. the secretary problem [25], modelling
the optimal time to call an election based on data [95], estimating the
solvency of governments with respect to their debt [86], and multi-armed
bandit problems [44]. Another important application of LSM is in the insurance
sector. In fact, LSM can be used to estimate the VaR (Value at Risk) [59] and
life insurance contracts [4] (see also [76] for a comparison of LSM with other
methods). The computational challenges of this domain were further highlighted
by recent European regulatory requirements [30, 29]. LSM is also often used
for solving Backward Stochastic Differential Equations (BSDE). Some numerical
algorithms for BSDE are two-steps stochastic procedures involving a
discretisation step where the solutions obtained at time $t$ of the BSDE are
projected onto a space obtained from the filtrations at time $t-1$. This step
involves a conditional expectation that cannot be calculated analytically, but
must be estimated using some approximation procedure. The idea of applying LSM
to BSDE was first introduced in [39] and further developed in [40, 60].
Recently, this method has been generalized to solve two-dimensional forward-
backward stochastic differential equations [61, 7].

####  1.1.3 Other algorithms for option pricing

LSM is not the only type of algorithm that can be used to price American
options [32]. Besides a few attempts to give an analytical formula under
certain conditions [56], the vast majority of them has been directed towards
giving numerical results, which we briefly discuss in this section. A simple
and well-known way of pricing American options is through the use of binomial
trees. While the origins of this technique are somewhat unclear [19], the
first articles that proposed the idea of binomial trees for pricing options
are considered to be [24, 84], with the first seminal ideas proposed in the
first edition of [87]. McKean [67] realized that the price of an American
option can be cast as a free boundary problem [77], which is a particular
partial differential equation that can be solved numerically. There is a
flurry of other methods to price American options based on partial
differential equations. We name a few approaches such as variational
inequalities [52, 8], linear complementary [50], and those related to free
boundaries [98]. However, as noted in [32], these methods often suffer from
the curse of dimensionality, as they require the computing time and the
storage to grow exponentially with the dimension of the underlying state
space. LSM is also not the only Monte Carlo approach for pricing American
options. One of the first works using Monte Carlo for option pricing is [13].
Reviews of Monte Carlo and other methods for the problem of American option
pricing can be found in [36, 58, 17]. In contrast to giving lower bounds for
the true optimal stopping value — as the LSM algorithm — Rogers [85] proposed
a method which leverages a dual problem, resulting in an upper bound for the
optimal stopping value. Last but not least, semi-analytical approaches for
American option pricing and optimal stopping time are also used [6].

###  1.2 Our results

We propose a quantum version of the LSM algorithm and perform the error
analysis considering a general set of expansion functions
$\\{e_{t,k}\\}_{t,k}$. We achieve a quadratic speedup in the approximation
error compared to the classical counterpart under some assumptions of
smoothness on the continuation values (which are explained below).

Our quantum algorithm assumes quantum sampling access to the Markov chain and
quantum query access to the functions $\\{z_{t}:E\to\mathbb{R}\\}_{t=0}^{T}$
and $\\{e_{t,k}:E\to\mathbb{R}\\}_{t\in[T-1],k\in[m]}$. Given these, it
follows the classical LSM in approximating the continuation values with a
suitable set of expansion functions. However, the dynamic programming is not
solved separately along different sampled paths, but in superposition along
all possible stochastic processes. More specifically, at any given time $t$,
the dynamic programming is solved in a backward fashion from time $T$ to $t$
by constructing a unitary that prepares the stopping times $\tau_{t}$ in
superposition via the mapping $|x\rangle\ket{\bf
0}\mapsto|x\rangle|\tau_{t}(x)\rangle$ for all $x\in E$. This resulting state
is used in the quantum subroutines for Monte Carlo to extract expectation
values $\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]$ that make up the vector $b_{t}$ and
which are required to continue solving the dynamic programming at the next
time step $t-1$. Such procedure is repeated until $t=0$, when the optimal
stopping time $\tau_{0}$ can be computed in superposition and thus the
quantity $\sup_{\tau}\mathbb{E}[Z_{\tau}]$ be finally approximated. We note
that, unlike the classical LSM, our quantum algorithm requires redoing all
previous dynamic programming steps before a given time $t$ in order to
progress into the next time step $t-1$. The final procedure involves
$O(T^{2})$ time steps instead of $O(T)$.

We prove that the classical LSM algorithm and our proposed quantum LSM
algorithm approximate the sought-after quantity $\mathcal{U}_{0}$ up to
additive accuracy with high probability. In the following, define
$\vec{e}_{t}(\cdot):=(e_{t,1}(\cdot),\dots,e_{t,m}(\cdot))^{\top}$.

######  Theorem 1.2 (Informal version of Corollary 3.15).

Consider a set of linearly independent functions
$\\{e_{t,k}:E\to\mathbb{R}\\}_{k=1}^{m}$ for each $t\in[T-1]$ and payoff
functions $\\{z_{t}:E\to\mathbb{R}\\}_{t=0}^{T}$. Then, for $\delta\in(0,1)$
and $\epsilon>0$, the classical and quantum LSM algorithms output
$\widetilde{\mathcal{U}}_{0}$ such that

$\displaystyle\operatorname{Pr}\left[|\widetilde{\mathcal{U}}_{0}-\mathcal{U}_{0}|\geq
5^{T}\left(\epsilon+\operatorname*{\max}_{0<t<T}\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}\right)\right]\leq\delta$

using time, respectively,
$\widetilde{O}\left(\frac{Tm^{6}}{\epsilon^{2}}\right)$ and
$\widetilde{O}\left(\frac{T^{2}m^{4}}{\epsilon}\right)$, up to ${\rm polylog}$
terms.

The error $\epsilon$ arises from the Monte Carlo subroutines and can be made
smaller by increasing the calls to the quantum inputs (or to the number of
sampled paths in the classical counterpart). Compared to the classical
algorithm, the number of oracle calls is quadratically less in the quantum
algorithm. The quantity
$\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}$
appearing in the theorem above is known as _approximation error_. This term
arises from approximating the continuation values by the $m$ expansion
functions and is a deterministic quantity implicitly dependent on $m$ and on
smoothness properties of the continuation values.

In order to obtain a final additive accuracy $\epsilon_{\rm final}$ for
$\widetilde{\mathcal{U}}_{0}$, we must resolve the implicit dependence of the
approximation error on $m$. This is done by considering specific sets of
expansion functions and assuming sufficiently good smoothness properties for
the continuation values. More specifically, for each $t\in[T-1]$ we consider
functions $\\{e_{t,k}:E\to\mathbb{R}\\}_{k=1}^{m}$ that generate the space
$\mathcal{R}_{q}$ of all polynomials of degree at most $q$, so that
$m=\binom{q+d}{d}$. We also assume that $\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\in
C^{n}$, i.e., the continuation values are $n$-differentiable functions. Then
it is possible to bound the approximation error
$\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}$
by using a Jackson-like inequality [51] and obtain the following result.

######  Theorem 1.3 (Informal version of Theorem 4.3).

For each $t\in[T-1]$, consider a set of linearly independent functions
$\\{e_{t,k}:E\to\mathbb{R}\\}_{k=1}^{m}$ that spans the space
$\mathcal{R}_{q}$ with $m=\binom{q+d}{d}$ and consider payoff functions
$\\{z_{t}:E\to\mathbb{R}\\}_{t=0}^{T}$. Assume that
$\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\in C^{n}$ for all $t\in\\{0,\dots,T-1\\}$,
where $n\leq q$. Then, for $\delta\in(0,1)$ and $\epsilon>0$, if
$q=\lceil(5^{T}/\epsilon)^{1/n}\rceil$, the classical and quantum LSM
algorithms output $\widetilde{\mathcal{U}}_{0}$ such that

$\displaystyle\operatorname{Pr}\left[|\widetilde{\mathcal{U}}_{0}-\mathcal{U}_{0}|\geq\epsilon\right]\leq\delta$

using time, respectively,
$\widetilde{O}\left((5^{T}/\epsilon)^{2+6d/n}\right)$ and
$\widetilde{O}\left((5^{T}/\epsilon)^{1+4d/n}\right)$, up to ${\rm polylog}$
terms.

We obtain similar results if the continuation values are Lipschitz continuous,
in which case the time complexities are the ones from the previous theorem
with $n=1$.

If the continuation values are $n$-times differentiable, for
$n=\Theta(\log(5^{T}/\epsilon)/\log\log(5^{T}/\epsilon))$, then we get the
sought-after quadratic improvement from $\widetilde{O}((5^{T}/\epsilon)^{2})$
classical runtime to $\widetilde{O}(5^{T}/\epsilon)$ quantum runtime, up to
polylog terms. We briefly note that such smoothness conditions on the
continuation values are not unusual in areas like finance. Indeed, the
continuation values can even be in $C^{\infty}$ in some models, e.g. Black-
Scholes [37, 94].

Very recently, Miyamoto [69] proposed a quantum LSM algorithm based on
Chebyshev interpolation through Chebyshev nodes and obtained a complexity of
$O(\epsilon^{-1}\log^{d}(1/\epsilon)\text{poly}\log\log(1/\epsilon))$. Our
approach, in contrast, is to project $\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]$ onto a
set of polynomials and is, for this reason, much more general. Moreover, our
final result is a _time_ complexity, while the result from [69] is a _query_
complexity on the number of unitaries called by all quantum routines for Monte
Carlo. Finally, Miyamoto [69] assumes that the continuation values are
analytical functions, i.e., are in $C^{\infty}$, while we only need to assume
$\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\in C^{n}$ for
$n=\Theta(\log(5^{T}/\epsilon)/\log\log(5^{T}/\epsilon))$ in order to recover
$\widetilde{O}(\epsilon^{-1})$ up to polylog factors. One downside of our
approach, though, is the presence of quantities that implicitly depend on the
underlying Markov chain.

As just mentioned, the full results behind the informal theorems above involve
parameters that depend on the underlying Markov chain such as the minimum
singular value of the matrices $A_{t}$. In order to explicitly work these
parameters out, we also study the case when the underlying Markov process
follows Brownian motion or geometric Brownian motion and obtain a simplified
version of our algorithm. In the case of Brownian motion, we choose Hermite
polynomials as the functions $\\{e_{t,k}:E\to\mathbb{R}\\}_{k=1}^{m}$ for each
$t\in[T-1]$, since they are orthogonal under the probability measure
underlying a Brownian motion. This means that the matrices $A_{t}$ are just
the identity. The final result is a mild reduction on the classical and
quantum time complexities to $\widetilde{O}((5^{T}/\epsilon)^{2+4d/n})$ and
$\widetilde{O}((5^{T}/\epsilon)^{1+7d/2n})$, respectively. For the geometric
Brownian motion, we pick suitable monomials that reduce the matrices $A_{t}$
to Vandermonde matrices, whose minimum singular value can be bounded. We
obtain the final classical and quantum complexities
$e^{O((5^{T}/\epsilon)^{2/n})}(5^{T}/\epsilon)^{2+12d/n}$ and
$e^{O((5^{T}/\epsilon)^{2/n})}(5^{T}/\epsilon)^{1+15d/2n}$, respectively. If
the continuation values are again $n$-times differentiable for
$n=\Theta(\log(5^{T}/\epsilon)/\log\log(5^{T}/\epsilon))$, then the classical
and quantum complexities for the Brownian motion setting reduce to the usual
$\widetilde{O}(\epsilon^{-2})$ and $\widetilde{O}(\epsilon^{-1})$,
respectively, while, for the geometric Brownian motion, they reduce to
$e^{O(\log^{c}(5^{T}/\epsilon))}(5^{T}/\epsilon)^{2}$ and
$e^{O(\log^{c}(5^{T}/\epsilon))}(5^{T}/\epsilon)$ for any constant $0<c<1$.
These results for the geometric Brownian motion are slightly weaker than the
usual $\widetilde{O}((5^{T}/\epsilon)^{2})$ and
$\widetilde{O}(5^{T}/\epsilon)$, since the bound on the minimum singular value
of the matrix $A_{t}$ is very sensitive to the degree $q$ of the chosen
monomials.

##  2 Preliminaries

Define $[n]:=\\{1,\dots,n\\}$. Given $b\in\mathbb{R}^{m}$ and
$A\in\mathbb{R}^{m\times m}$ for some $m\in\mathbb{N}$, let
$\|b\|_{2}:=\sqrt{\sum_{i=1}^{m}b_{i}^{2}}$ and $\|A\|_{2}:=\sigma_{\max}(A)$
be the vector and matrix norms, respectively, where $\sigma_{\max}(A)$ is the
maximum singular value of $A$. We shall denote by $\sigma_{\min}(A)$ the
minimum singular value of $A$. Given an event $\mathscr{E}$, we denote by
$\mathbf{1}\\{\mathscr{E}\\}$ the indicator function
$\mathbf{1}\\{\mathscr{E}\\}=1$ if $\mathscr{E}$ is true and $0$ if not. Let
$\omega_{\ast}$ denote the matrix multiplication exponent.

Let $(\Omega,\mathcal{F},\mathbb{P})$ be a probability space and
$\\{0,1,\dots,T\\}$ an index set with $T\in\mathbb{N}$. For every
$t\in\\{0,\dots,T\\}$, let $\mathcal{F}_{t}$ be a sub $\sigma$-algebra on
$\Omega$. If $\mathcal{F}_{k}\subseteq\mathcal{F}_{l}$ for all $k\leq l$, then
$(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t=0}^{T},\mathbb{P})$ is a filtered
probability space with filtration $(\mathcal{F}_{t})_{t=0}^{T}$. Consider a
discrete-time stochastic process $\mathbf{X}=(X_{t})_{t=0}^{T}$, assumed to be
Markovian, defined on a filtered probability space
$(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t=0}^{T},\mathbb{P})$ and with state
space $(E,\mathcal{E})$, where $E\subseteq\mathbb{R}^{d}$. The set $E$ is
equipped with its natural Borel $\sigma$-algebra inherited from
$\mathbb{R}^{d}$, which we denote as $\mathcal{E}$. For $i\in[d]$, let
$X_{t,i}:\Omega\to\mathbb{R}$ be the $i$-th component of $X_{t}$. We shall
assume that $\mathbf{X}$ is _adapted_ with respect to
$(\mathcal{F}_{t})_{t=0}^{T}$, meaning that each $X_{t}$ is
$\mathcal{F}_{t}$-measurable. We shall also assume that $X_{0}=x_{0}$ is
deterministic and, therefore, sometimes write the Markov chain
$(X_{t})_{t=1}^{T}$ as starting from $t=1$. Each element $X_{t}$ for
$t\in\\{0,\dots,T\\}$, called the underlying process at time $t$, gives rise
to an image probability measure (also called pushforward measure) $\rho_{t}$
in $E\subseteq\mathbb{R}^{d}$, i.e.,
$\rho_{t}(Y)=\mathbb{P}[\omega\in\Omega:X_{t}(\omega)\in Y]$ for any
$Y\in\mathcal{E}$ (note that $\rho_{0}$ is the probability measure that
assigns measure $1$ to the singleton set containing $x_{0}$). We denote by
$\mathbf{X}_{t}=(X_{j})_{j=t}^{T}$ the last $T-t+1$ random variables in the
stochastic process. Let $L^{2}(E,\rho_{t})=L^{2}(\rho_{t})$ be the set of
squared integrable functions with norm
$\|f\|_{L^{2}(\rho_{t})}=\sqrt{\mathbb{E}_{\rho_{t}}[|f(X_{t})|^{2}]}$.
Moreover, define the uniform norm $\|f\|_{u}=\sup\\{|f(s)|:s\in E\\}$ for
$f:E\to\mathbb{R}$.

###  2.1 Stochastic optimal stopping problems

The stochastic process $(X_{t})_{t=0}^{T}$ with state space
$E\subseteq\mathbb{R}^{d}$ represents the randomness of the problem under
consideration. In economics, this process can model $d$ factors driving the
economy and in finance this process can represent the values of some $d$
underlying set of securities. For our optimal stopping problem, we are
interested in evaluating some function (which we may call “payoff function”)
of this stochastic process. More formally:

######  Definition 2.1 (Payoff process [34]).

A payoff process is a non-negative adapted process $(Z_{t})_{t=0}^{T}$ on the
filtered probability space
$(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t=0}^{T},\mathbb{P})$. In this work,
we consider an adapted payoff process $(Z_{t})_{t=0}^{T}$ obtained from the
Markov chain $(X_{t})_{t=0}^{T}$ as $Z_{t}=z_{t}(X_{t})$ for Borel functions
$z_{t}(\cdot)\in L^{2}(E,\rho_{t})$, $t\in\\{0,\dots,T\\}$.

An example of a non-negative adapted process is the American contingent claim.
Common examples of American options are the _put_ and _call_ options on a
single asset given by the $Z^{\rm put}=z^{\rm put}(X)\coloneqq\max\\{0,K-X\\}$
and $Z^{\rm call}=z^{\rm call}(X)\coloneqq\max\\{0,X-K\\}$, respectively,
where the constant $K>0$ is called strike price of the option.

Next, we define the notion of a stopping time (also called Markov time). A
stopping time is any random variable that selects one of the possible times
$\\{0,1,\dots,T\\}\cup\\{+\infty\\}$ and satisfies a measureability condition.

######  Definition 2.2 (Stopping time).

A _stopping time_ is a function
$\tau:\Omega\to\\{0,1,\dots,T\\}\cup\\{+\infty\\}$ such that
$\\{\omega\in\Omega|\tau(\omega)=t\\}\in\mathcal{F}_{t}$ for
$t\in\\{0,\dots,T\\}$. The payoff obtained by using $\tau$ is
$Z_{\tau}(\omega):=Z_{\tau(\omega)}(\omega)$. Let
$\mathbb{T}_{t}:=\\{\tau|\tau~{}\text{is a stopping time with}~{}t\leq\tau\leq
T\\}$ be the set of all stopping times taking values in $[t,T]$. A stopping
time $\tau^{\ast}\in\mathbb{T}_{t}$ is called _optimal_ in the interval
$[t,T]$ if

$\displaystyle\mathbb{E}[Z_{\tau^{\ast}}|X_{t}]=\operatorname*{ess~{}sup}_{\tau\in\mathbb{T}_{t}}\mathbb{E}[Z_{\tau}|X_{t}].$

The maximization is expressed via an essential supremum such that the null
sets of the probability measure do not affect the result. For more details on
the essential supremum see [34, Appendix A.5]. The optimal stopping problem
then consists in finding a stopping time $\tau$ that maximizes the expected
value payoff.

######  Problem 2.3 (Precise formulation of Problem 1.1).

Let $(Z_{t})_{t=0}^{T}$ be a payoff process. For $\epsilon_{\rm final}>0$,
approximate the exact value
$\operatorname*{sup}_{\tau\in\mathbb{T}_{0}}\mathbb{E}[Z_{\tau}]$ to additive
accuracy $\epsilon_{\rm final}$ with high probability.

A well-studied solution strategy for the above problem statement is based on
dynamic programming for a set of stopping times. A crucial concept is the
_Snell envelope_ [72, 34].

######  Definition 2.4 (Snell envelope).

Let $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t=0}^{T},\mathbb{P})$ be a
filtered probability space and $(X_{t})_{t=0}^{T}$ a Markov chain with state
space $E\subseteq\mathbb{R}^{d}$. Let $(Z_{t})_{t=0}^{T}$ be a payoff process
where $Z_{t}=z_{t}(X_{t})$ for some $z_{t}\in L^{2}(E,\rho_{t})$,
$t\in\\{0,\dots,T\\}$. The Snell envelope
$\mathcal{U}_{t}:\Omega\to\mathbb{R}$ of $Z_{t}$ is defined as

$\begin{cases}\mathcal{U}_{T}=Z_{T},\\\
\mathcal{U}_{t}=\max\left\\{Z_{t},\mathbb{E}[\mathcal{U}_{t+1}|X_{t}]\right\\},&0\leq
t\leq T-1.\end{cases}$

Define the stopping times $\tau_{t}:=\min\\{u\geq
t~{}|~{}\mathcal{U}_{u}=Z_{u}\\}$. The Snell envelope is related to the
maximal expected payoff according to the next theorem: $\tau_{t}$ maximizes
the expectation of $Z_{\tau}$ among all $\tau\in\mathbb{T}_{t}$, i.e., that
$\tau_{t}$ are optimal stopping times (in their respective intervals).

######  Theorem 2.5 ([34, Theorem 6.18]).

The Snell envelope $\mathcal{U}_{t}$ of $Z_{t}$ satisfies

$\displaystyle\mathcal{U}_{t}=\mathbb{E}[Z_{\tau_{t}}|X_{t}]=\operatorname*{ess~{}sup}_{\tau\in\mathbb{T}_{t}}\mathbb{E}[Z_{\tau}|X_{t}],$

where $\tau_{t}:=\min\\{u\geq t~{}|~{}\mathcal{U}_{u}=Z_{u}\\}$. In
particular,
$\mathcal{U}_{0}=\mathbb{E}[Z_{\tau_{0}}]=\operatorname*{sup}_{\tau\in\mathbb{T}_{0}}\mathbb{E}[Z_{\tau}]=\max\\{Z_{0},\mathbb{E}[Z_{\tau_{1}}]\\}$.

Hence, finding an approximate $\mathcal{U}_{0}$ solves our Problem 2.3. In
order to solve the dynamic programming behind the Snell envelope, it is more
convenient to recast the dynamic programming in terms of the optimal stopping
times $\tau_{t}$ (rather than in terms of value functions) as follows.

######  Theorem 2.6.

The dynamic programming principle in Definition 2.4 can be recast in terms of
the stopping times $\tau_{t}=\min\\{u\geq t~{}|~{}\mathcal{U}_{u}=Z_{u}\\}$ as

$\displaystyle\begin{cases}\tau_{T}=T,\\\
\tau_{t}=t\mathbf{1}\\{Z_{t}\geq\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\\}+\tau_{t+1}\mathbf{1}\\{Z_{t}<\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\\},&0\leq
t\leq T-1.\end{cases}$

###### Proof.

The case $t=T$ is trivial. Assume then $t<T$. Note that
$\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]=\mathbb{E}[\mathbb{E}[Z_{\tau_{t+1}}|X_{t+1}]|X_{t}]$
because of the tower property of the expectation value with the filtration
generated by $X_{t}$. In addition,
$\mathbb{E}[Z_{\tau_{t+1}}|X_{t+1}]=\mathcal{U}_{t+1}$ from Theorem 2.5.
Hence, if $Z_{t}\geq\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]$, then
$Z_{t}\geq\mathbb{E}[\mathcal{U}_{t+1}|X_{t}]$. This latter statement, by the
definition of the Snell envelope, implies $\mathcal{U}_{t}=Z_{t}$ and then
$\tau_{t}=t$. On the other hand, if $Z_{t}<\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]$,
then, $Z_{t}<\mathbb{E}[\mathcal{U}_{t+1}|X_{t}]\implies
Z_{t}\neq\mathcal{U}_{t}$, and so $\tau_{t}=\min\\{u\geq
t~{}|~{}Z_{u}=\mathcal{U}_{u}\\}=\min\\{u\geq
t+1~{}|~{}Z_{u}=\mathcal{U}_{u}\\}=\tau_{t+1}$. ∎

The stopping time $\tau_{0}$ maximizes $\mathbb{E}[Z_{\tau}]$ in Problem 2.3.
The conditional expectations $\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]$ are called
_continuation values_.

###  2.2 The least squares Monte Carlo algorithm

The LSM algorithm consists in solving the dynamic programming in Theorem 2.6
by means of two approximations. The first one is to approximate the
continuation values $\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]$ using a set of
measurable real-valued functions in $L^{2}(E,\rho_{t})$, e.g. by projection
onto a finite-dimensional set of linearly independent polynomials. Let
$\mathscr{H}_{0}\subseteq\mathbb{R}$ and let, for $t\in[T-1]$,
$\mathscr{H}_{t}\subseteq L^{2}(E,\rho_{t})$ be a subset of real-valued
functions on $E$, called _approximation architecture_ or hypothesis class,
that will be used to approximate the continuation values. By approximating
$\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]$ by $f_{t}\in\mathscr{H}_{t}$ for each
$t\in\\{0,\dots,T-1\\}$, we can write the approximate dynamic programming as

$\displaystyle\begin{cases}\widetilde{\tau}_{T}=T,\\\
\widetilde{\tau}_{t}=t\mathbf{1}\\{Z_{t}\geq
f_{t}\\}+\widetilde{\tau}_{t+1}\mathbf{1}\\{Z_{t}<f_{t}\\},&0\leq t\leq
T-1.\end{cases}$ (1)

Note that $\widetilde{\tau}_{t}=\widetilde{\tau}_{t}(f_{t},\dots,f_{T-1})$
depends on the approximation architecture.

The second approximation of the algorithm is to numerically evaluate the
approximations $f_{t}$ in $L^{2}(\rho_{t})$ by a Monte Carlo procedure. We
sample $N$ independent paths
$(X_{t}^{(1)})_{t=0}^{T},\dots,(X_{t}^{(N)})_{t=0}^{T}$ of the Markov chain
$\mathbf{X}=(X_{t})_{t=0}^{T}$ and denote by $Z^{(n)}_{t}=z_{t}(X^{(n)}_{t})$
the associated payoffs conditioned on $X_{t}^{(n)}$, where $z_{t}\in
L^{2}(E,\rho_{t})$, $t\in\\{0,\dots,T\\}$. Write the random variables of the
last $T-t+1$ elements of all the sampled Markov chains by
$\mathbf{X}_{t}^{(N)}=(X_{t}^{(1)},\dots,X_{t}^{(N)},X_{t+1}^{(1)},\dots,X_{t+1}^{(N)},\dots,X_{T}^{(1)},\dots,X_{T}^{(N)})$.
For each path, the dynamic programming in Eq. (1) is solved recursively by
approximating the continuation values in $\mathscr{H}_{t}$ via a least-square
estimator. The result is sampled stopping times $\widetilde{\tau}_{t}^{(n)}$
that $\widetilde{\tau}_{t}$ takes on each random path. We stress that, due to
the recursive nature of Eq. (1), the stopping times
$\widetilde{\tau}_{t}^{(n)}$ will depend on $\mathbf{X}_{t}^{(N)}$, and
consequently also the payoffs
$Z^{(n)}_{\widetilde{\tau}^{(n)}_{t}}=z_{\widetilde{\tau}^{(n)}_{t}(\mathbf{X}_{t}^{(N)})}\big{(}X^{(n)}_{\widetilde{\tau}^{(n)}_{t}(\mathbf{X}_{t}^{(N)})}\big{)}$.
The dependence on $\mathbf{X}_{t}^{(N)}$ should be clear from the context and
therefore we shall simply write $Z^{(n)}_{\widetilde{\tau}^{(n)}_{t}}$. In
summary, combining both the approximation architecture and the Monte Carlo
sampling, at each $t\in[T-1]$ we take $f_{t}\in\mathscr{H}_{t}$, depending on
$\mathbf{X}_{t}^{(N)}$, satisfying

$\displaystyle\frac{1}{N}\sum_{n=1}^{N}\left(Z^{(n)}_{\widetilde{\tau}_{t+1}^{(n)}}-f_{t}(X_{t}^{(n)})\right)^{2}\leq\epsilon+\operatorname*{\inf}_{g\in\mathscr{H}_{t}}\frac{1}{N}\sum_{n=1}^{N}\left(Z^{(n)}_{\widetilde{\tau}_{t+1}^{(n)}}-g(X_{t}^{(n)})\right)^{2}$
(2)

for some given $\epsilon\geq 0$. It might be the case that an exact minimizer
of the above optimization problem does not exist (the infimum does not belong
to $\mathscr{H}_{t}$) or is hard to compute, meaning that an
$\epsilon$-approximation could be used. Given the choice of
$f_{t}\in\mathscr{H}_{t}$, it is then used in Eq. (1) to obtain
$\widetilde{\tau}_{t}^{(n)}$, and so on recursively. At the end of the
recursion we can take
$f_{0}=\frac{1}{N}\sum_{n=1}^{N}Z^{(n)}_{\widetilde{\tau}_{1}^{(n)}}$ as an
exact minimizer, since $X_{0}$ is constant, and obtain the approximation
$\widetilde{\mathcal{U}}_{0}$ to $\mathcal{U}_{0}$ as

$\displaystyle\widetilde{\mathcal{U}}_{0}=\max\left\\{Z_{0},\frac{1}{N}\sum_{n=1}^{N}Z^{(n)}_{\widetilde{\tau}_{1}^{(n)}}\right\\}.$
(3)

In this paper we shall be particularly interested in finite-dimensional linear
approximation architectures, for which an exact minimizer exists in Eq. (2).
Consider then a set $\\{e_{t,k}:E\to\mathbb{R}\\}_{k=1}^{m}$ of $m$ linearly
independent measurable real functions and take the vector space generated by
them as our approximation architecture $\mathscr{H}_{t}$, $t\in[T-1]$.
Therefore, the infimum in Eq. (2) is attained by projecting the continuation
values onto $\mathscr{H}_{t}$ as $\alpha_{t}\cdot\vec{e}_{t}(X_{t})$, where
$\vec{e}_{t}(\cdot):=(e_{t,1}(\cdot),\dots,e_{t,m}(\cdot))^{\top}$ and the
$m$-dimensional vector $\alpha_{t}$, the projection coefficients, is the
least-square estimator given by [21]

$\displaystyle\alpha_{t}=\arg\operatorname*{\min}_{a\in\mathbb{R}^{m}}\mathbb{E}\big{[}(Z_{\widetilde{\tau}_{t+1}}-a\cdot\vec{e}_{t}(X_{t}))^{2}\big{]}.$

Given the assumption that $\\{e_{t,k}\\}_{k=1}^{m}$ are linearly independent
for each $t\in[T-1]$, the vector $\alpha_{t}\in\mathbb{R}^{m}$ has the
explicit expression

$\displaystyle\alpha_{t}=A_{t}^{-1}b_{t}\text{ where
}b_{t}=\mathbb{E}[Z_{\widetilde{\tau}_{t+1}}\vec{e}_{t}(X_{t})]$ (4)

and the $m\times m$ matrix $A_{t}$ has coefficients

$\displaystyle(A_{t})_{k,l}=\mathbb{E}[e_{t,k}(X_{t})e_{t,l}(X_{t})].$ (5)

Often it is hard to compute $\alpha_{t}$ and $A_{t}$ exactly. As previously
mentioned, the LSM algorithm approximates these by Monte Carlo sampling,

$\displaystyle\widetilde{\alpha}_{t}=\widetilde{A}_{t}^{-1}\frac{1}{N}\sum_{n=1}^{N}Z^{(n)}_{\widetilde{\tau}^{(n)}_{t+1}}\vec{e}_{t}(X_{t}^{(n)})$
(6)

and

$\displaystyle(\widetilde{A}_{t})_{k,l}=\frac{1}{N}\sum_{n=1}^{N}e_{t,k}(X^{(n)}_{t})e_{t,l}(X^{(n)}_{t}).$
(7)

More generally, though, any good approximation $\widetilde{\alpha}_{t}$ and
$\widetilde{A}_{t}$ to $\alpha_{t}$ and $A_{t}$, respectively, is valid, and
we shall not restrict the notation $\widetilde{\alpha}_{t}$ and
$\widetilde{A}_{t}$ to only mean the above sampled quantities.

We have introduced the quantities that are important for the LSM algorithm. To
present the algorithm, we first specify the input model, starting with
sampling access to the Markov chain.

######  Definition 2.7 (Sampling access to Markov chain).

Given a Markov chain $(X_{t})_{t=1}^{T}$ on a probability space
$(\Omega,\mathbb{P})$ and with state space $E\subseteq\mathbb{R}^{d}$, we
define sampling access as the ability to draw a sample $\omega\in\Omega$
according to $\mathbb{P}$ and observe the value $X_{t}(\omega)$ for some
$t\in[T]$. One sample costs time $\mathcal{T}_{\rm samp}$.

Furthermore, query access to a function is defined as being able to evaluate a
function for a given input.

######  Definition 2.8 (Query access to function).

Let $E\subseteq\mathbb{R}^{d}$ and $h:E\to\mathbb{R}$ be a function. We define
query access as the ability to observe the value $h(x)$ for any given $x\in
E$. One query costs time $\mathcal{T}_{h}$.

Here, we assume that the functions of the approximation architecture and
functions for the payoff take time $\mathcal{T}_{e}$ and $\mathcal{T}_{z}$,
respectively, to evaluate. Both sampling and function access have natural
quantum extensions, as will be defined in Section 3.

We are now in the position to present the classical LSM algorithm in Algorithm
1. Since we focus on the case where the approximation architectures
$\mathscr{H}_{t}$ are finite-dimensional and linear, we write Algorithm 1 for
this particular case.

Algorithm 1 Classical LSM algorithm for optimal stopping problem

1:Integer $N\in\mathbb{N}$. Sampling access to Markov chain
$(X_{t})_{t=0}^{T}$ defined on a sample space $\Omega$ and with state space
$E\subseteq\mathbb{R}^{d}$. Query access to
$\\{z_{t}:E\to\mathbb{R}\\}_{t=0}^{T}$ and
$\\{e_{t,k}:E\to\mathbb{R}\\}_{t\in[T-1],k\in[m]}$, where
$\\{e_{t,k}\\}_{k=1}^{m}$ are linearly independent for each $t\in[T-1]$. Let
$\vec{e}_{t}(\cdot):=(e_{t,1}(\cdot),\dots,e_{t,m}(\cdot))^{\top}$.

2:Sample $N$ independent paths $(X^{(1)}_{t},\dots,X^{(N)}_{t})_{t=0}^{T}$.

3:Query associated payoffs $(Z^{(1)}_{t},\dots,Z^{(N)}_{t})_{t=0}^{T}$ and
values
$(e_{t,k}(X_{t}^{(1)}),\dots,e_{t,k}(X_{t}^{(N)}))_{t\in[T-1],k\in[m]}$.

4:Compute the matrices $\\{\widetilde{A}_{t}\\}_{t=1}^{T-1}$ with entries as
in Eq. (7).

5:Compute the inverses $\\{\widetilde{A}_{t}^{-1}\\}_{t=1}^{T-1}$.

6:Set $\widetilde{\tau}^{(n)}_{T}=T$ for $n\in[N]$.

7:for $t=T-1$ to $1$ do

8: Calculate the vector
$\widetilde{\alpha}_{t}=\widetilde{A}_{t}^{-1}\frac{1}{N}\sum_{n=1}^{N}Z^{(n)}_{\widetilde{\tau}^{(n)}_{t+1}}\vec{e}_{t}(X_{t}^{(n)})$.

9: Calculate
$\widetilde{\tau}^{(n)}_{t}=t\mathbf{1}\\{Z_{t}^{(n)}\geq\widetilde{\alpha}_{t}\cdot\vec{e}_{t}(X_{t}^{(n)})\\}+\widetilde{\tau}_{t+1}^{(n)}\mathbf{1}\\{Z_{t}^{(n)}<\widetilde{\alpha}_{t}\cdot\vec{e}_{t}(X_{t}^{(n)})\\}$,
$n\in[N]$.

10:end for

11:Output
$\widetilde{\mathcal{U}}_{0}:=\max\Big{\\{}Z_{0},\frac{1}{N}\sum_{n=1}^{N}Z^{(n)}_{\widetilde{\tau}_{1}^{(n)}}\Big{\\}}$.

###  2.3 Error analysis of least squares Monte Carlo algorithm

As previously mentioned in the introduction, there has been a considerable
amount of work in understanding the convergence and expected sample error of
the LSM algorithm, starting with Clément et al. [21] two decades ago until the
more recent results of Zanger [106]. We first review a general result of
Zanger [106] for any uniformly bounded approximation architecture with finite
VC dimension.

######  Theorem 2.9 ([106, Theorem 3]).

Consider the LSM algorithm (Algorithm 1) with $N$ random sampled paths, $T$
time steps, and uniformly bounded approximation architectures
$\mathscr{H}_{t}$, $t\in\\{0,\dots,T-1\\}$, with finite VC dimension. Assume
that an exact minimizer is always chosen in Eq. (2). Then the expected error
from the LSM algorithm is

$\displaystyle\mathbb{E}[|\widetilde{\mathcal{U}}_{0}-\mathcal{U}_{0}|]=O\left(6^{T}\left(\sqrt{\frac{\log{N}}{N}}+\max_{0\leq
t<T}\inf_{f\in\mathscr{H}_{t}}\|f-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}\right)\right).$

The term
$\inf_{f\in\mathscr{H}_{t}}\|f-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}$
is known in learning theory as the _approximation error_ of the approximation
architecture, and is a deterministic quantity that conveys the minimum error
that is possible from using $\mathscr{H}_{t}$.

Often in the study of convergence of the LSM algorithm, previous works [21,
32, 106] do not take into consideration the time complexity of finding an
exact (or $\epsilon$-approximation) minimizer in Eq. (2), but instead consider
only the number of sampled paths $N$ as a measure of complexity. It is thus
not immediate to compare, e.g. the above theorem, with our quantum LSM
algorithm in the next section, since sampled paths can be reused in different
parts of an algorithm, but calls to a quantum oracle generally cannot.
Nonetheless, the time complexity of Algorithm 1 can be precisely analyzed, due
to the closed form for the exact minimizer, and it is given in the next
theorem. This result will serve as our proxy when comparing our quantum
algorithm to its classical counterpart. We shall postpone the proof until the
end of next section, since it is very similar to Theorems 3.11 and 3.13 on our
quantum algorithm.

######  Theorem 2.10.

Within the setting of Algorithm 1, consider $N$ independent sample paths with
sample cost $\mathcal{T}_{\rm samp}$ and let the linearly independent
functions $\\{e_{t,k}:E\to\mathbb{R}\\}_{k=1}^{m}$ for each $t\in[T-1]$ be
such that $L:=\max_{t\in[T-1],k\in[m]}\|e_{t,k}\|_{L^{2}(\rho_{t})}$ and have
query cost $\mathcal{T}_{e}$. Also consider
$\\{z_{t}:E\to\mathbb{R}\\}_{t=0}^{T}$ with
$R:=\max_{t\in[T]}\|z_{t}\|_{u}<\infty$ and query cost $\mathcal{T}_{z}$.
Moreover, let $\sigma_{\min}:=\min_{t\in[T-1]}\sigma_{\min}(A_{t})>0$ and
assume that $\sqrt{m}RL/\sigma_{\min}\geq 1$. Then, for any
$\epsilon\in(0,\sigma_{\min}/2]$, Algorithm 1 runs for
$O(Tm^{2}N+Tm^{\omega_{\ast}}+TN(\mathcal{T}_{\rm
samp}+\mathcal{T}_{z}+m\mathcal{T}_{e}))$ time and returns an estimate
$\widetilde{\mathcal{U}}_{0}$ such that

$\displaystyle\operatorname{Pr}\left[\big{|}\widetilde{\mathcal{U}}_{0}-\mathbb{E}[Z_{\tau_{0}}]\big{|}\geq
5^{T}\left(\frac{4\epsilon
mRL^{2}}{\sigma_{\min}^{2}}+\operatorname*{\max}_{0<t<T}\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\tau_{t+1}}|X_{k}]\|_{L^{2}(\rho_{t})}\right)\right]\leq
6m^{2}e^{-2N\epsilon^{2}/m^{2}}.$

If sampling and query costs are constant, then the complexity is
$O(Tm^{2}N+Tm^{\omega_{\ast}})$. The factor $Tm^{2}$ comes from computing
$\\{A_{t}\\}_{t=1}^{T-1}$ and accounts for runtime instead of only number of
samples.

##  3 Quantum algorithm

In this section we shall present our quantum algorithm, which is based on the
classical LSM algorithm (Algorithm 1). Before we discuss it, we review our
computational model, input assumptions, and the quantum algorithm for Monte
Carlo used in this work. In what follows, for simplicity, we suppose that
$\ket{\bf 0}$ describes a register with sufficiently many qubits initialized
in the all-$0$ state.

###  3.1 Computational model

In this subsection, we address our quantum computational model. We work in the
standard circuit model of quantum computation [73]. Aside from these standard
assumptions, we take the following additional assumptions on the computational
model.

##### Arithmetic model

In our work, we perform the arithmetic computations on the quantum computer by
using a fixed point representation for real numbers. We assume that we can
have enough qubits for storing these numbers, represented as bit strings using
the following definition. We also assume to work with enough precision so that
numerical errors in the computation are negligible, and will not impact the
final output of our algorithm.

######  Definition 3.1 (Fixed-point encoding of real numbers [83]).

Let $c_{1},c_{2}$ be positive integers, and $a\in\\{0,1\\}^{c_{1}}$,
$b\in\\{0,1\\}^{c_{2}}$, and $s\in\\{0,1\\}$ be bit strings. Define the
rational number as:

$\mathcal{Q}(a,b,s):=(-1)^{s}\left(2^{c_{1}-1}a_{c_{1}}+\dots+2a_{2}+a_{1}+\frac{1}{2}b_{1}+\dots+\frac{1}{2^{c_{2}}}b_{c_{2}}\right)\in[-R,R],$
(8)

where $R:=2^{c_{1}}-2^{-c_{2}}$. If $c_{1},c_{2}$ are clear from the context,
we can use the shorthand notation for a number $z:=(a,b,s)$ and write
$\mathcal{Q}(z)$ instead of $\mathcal{Q}(a,b,s)$. Given an $n$-dimensional
vector $v\in(\\{0,1\\}^{c_{1}}\times\\{0,1\\}^{c_{2}}\times\\{0,1\\})^{n}$ the
notation $\mathcal{Q}(v)$ means an $n$-dimensional vector whose $j$-th
component is $\mathcal{Q}(v_{j})$, for $j\in[n]$.

The choice of values for $c_{1}$ and $c_{2}$ depends on the choice of input
functions used when running the algorithm. For the purposes of optimizing the
quantum circuit, these constants can be changed dynamically in various steps
of the computation. While analyzing how error propagates and accumulates
throughout the operations in the quantum circuit is essential to ensure a
correct estimation of the final result, this analysis can only be done for a
given choice of input functions. We avoid the analysis of such details by
using the quantum arithmetic model as in Definition 3.2. A standard result is
that any Boolean function can be reversibly computed. Any reversible
computation can be realized with a circuit involving negation and three-bit
Toffoli gates. Such a circuit can be turned into a quantum circuit with
single-qubit NOT gates and three-qubit Toffoli gates. Since most circuits for
arithmetic operations operate with a number of gates of
$O(\text{poly}(c_{1},c_{2}))$ this implies a number of quantum gates of
$O(\text{poly}(c_{1},c_{2}))$ for the corresponding quantum circuit.

######  Definition 3.2 (Quantum arithmetic model).

Given $c_{1},c_{2}\in\mathbb{N}$ specifying fixed-point precision numbers as
in Definition 3.1, we say we use a quantum arithmetic model of computation if
the four arithmetic operations can be performed in constant time in a quantum
computer.

In our computational model we do not include the cost for performing
operations described in our arithmetic model. For instance, a central
computational step of the quantum algorithm is the circuit computing the
stopping times $\widetilde{\tau}_{t}(x)$, but as the circuit depth depends
polynomially on $c_{1}$ and $c_{2}$, we do not take into account this cost
when stating our runtime. For an example of a resource estimation for a
financial problem that takes into account the cost of arithmetic operations in
fixed-point precision, we refer to [18].

##### Quantum input access

We assume that we have quantum oracles for certain input functions. The
classical algorithm assumes access to two different kinds of oracles. The
first is an oracle that allows us to obtain samples from the Markov chain
$(X_{t})_{t=0}^{T}$. The second kind of oracle is evaluating the functions
$\\{z_{t}\\}_{t=0}^{T}$ and $\\{e_{t,k}\\}_{t\in[T-1],k\in[m]}$. We assume
access to the quantum versions of these oracles (formalized below). The first
kind of quantum oracle prepares a quantum state that is in a superposition
over the different outcomes of the Markov chain, weighted by amplitudes which
are square roots of their classical probabilities. A measurement in the
computational basis of such a state obtains a single sample with the
corresponding probability and hence directly recovers a single use of the
classical sampling access. The second kind of quantum oracle evaluates a given
function in superposition over its inputs. While the functions
$\\{e_{t,k}\\}_{t\in[T-1],k\in[m]}$ are usually chosen to be low-degree
polynomials (and thus admit efficient classical and quantum circuits with gate
complexity proportional to the degree of the polynomial), the functions
$\\{z_{t}\\}_{t=0}^{T}$ might be arbitrarily complex. Usually the complexity
of these functions is not discussed in classical literature, and we use
placeholders for their evaluation cost.

######  Definition 3.3 (Quantum sampling access to a Markov chain).

Let $(X_{t})_{t=1}^{T}$ be a Markov chain defined on a filtered probability
space $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t=1}^{T},\mathbb{P})$, for a
finite $\Omega$, assuming values in a finite state space
$E\subseteq\mathbb{R}^{d}$. Given $x\in E^{T}$, let
$p(x):=\mathbb{P}[X_{1}=x_{1}]\prod_{t=1}^{T-1}\mathbb{P}[X_{t+1}=x_{t+1}|X_{t}=x_{t}]$
be the probability that $X_{1}=x_{1},\dots,X_{T}=x_{T}$. Let $\mathcal{H}$ be
a finite-dimensional Hilbert space with basis $\\{|x\rangle\\}_{x\in E^{T}}$.
We say that we have quantum sampling access to $(X_{t})_{t=1}^{T}$ if we are
given an oracle $U$ on $\mathcal{H}$ such that $U|\mathbf{0}\rangle=\sum_{x\in
E^{T}}\sqrt{p(x)}|x\rangle$. One application of $U$ costs
$T\cdot\mathcal{T}_{\rm samp}$ time. If $T=1$, we say that we have quantum
sampling access to a random variable $X$ if we are given an oracle $U$ on
$\mathcal{H}$ such that $U|\mathbf{0}\rangle=\sum_{x\in
E}\sqrt{p(x)}|x\rangle$, where $p(x):=\mathbb{P}[X=x]$ is the probability that
$X=x$.

We note the alternative definition of the unitary $U$ such that
$U|\mathbf{0}\rangle=\sum_{x\in E^{T}}\sqrt{p(x)}|x\rangle|\psi_{x}\rangle$
for unknown garbage unit states $|\psi_{x}\rangle$. Such garbage states do not
change our analysis, so we shall ignore them and work with the unitary $U$
from Definition 3.3.

Even though we assume the existence of the oracle $U$, constructing such
unitary is an important question on its own. A few methods have been proposed
in order to tackle such problem, one of the most famous is due to Grover and
Rudolph [43] (see [65] for recent improvements on the Grover-Rudolph method),
which loads into a quantum computer a discretization of a distribution with
density function $p(x)$. More specifically, it creates the quantum state
$\sum_{i=0}^{2^{n}-1}\sqrt{p_{i}^{(n)}}|i\rangle$ with
$p_{i}^{(n)}=\int_{x_{i}^{(n)}}^{x_{i+1}^{(n)}}p(x)\text{d}x$ by recursively
(on $n$) computing quantities like
$f_{n}(i)=\int_{x_{i}^{(n)}}^{(x_{i}^{(n)}+x_{i+1}^{(n)})/2}p(x)\text{d}x\big{/}\int_{x_{i}^{(n)}}^{x_{i+1}^{(n)}}p(x)\text{d}x$
for $i=0,\dots,2^{n}-1$. It is also possible to perform simple Taylor
approximations on $f_{n}(i)$ when $n$ is sufficiently large (see [55, Equation
(35)]). We briefly note that the issues about the Grover-Rudolph method
recently pointed out by [48] only arise when one needs to _sample_ from the
distribution $p(x)$ in order to compute $f_{n}(i)$, which is not the case in
many situations, e.g. in finance.

######  Definition 3.4 (Quantum access to a function).

Let $E\subseteq\mathbb{R}^{d}$ be a finite set and let $\mathcal{H}$ be a
Hilbert space with basis $\\{|x\rangle\\}_{x\in E}$. Given
$h:E\to\\{0,1\\}^{n}$, we say that we have $(V_{h},\mathcal{T}_{h})$-quantum
access to $h$ if we have access to a quantum circuit $V_{h}$ on
$\mathcal{H}\otimes\mathbb{C}^{2^{n}}$ such that
$V_{h}|x\rangle|b\rangle=|x\rangle|b\oplus h(x)\rangle$ for any bit string
$b\in\\{0,1\\}^{n}$. One application of $V_{h}$ costs time $\mathcal{T}_{h}$.

Computational Model | Classical / Quantum
---|---
Arithmetic operation | $O(1)$
Query functions $z$ | $\mathcal{T}_{z}$
Query functions $e$ | $\mathcal{T}_{e}$
Sampling a Markov chain | $T\cdot\mathcal{T}_{\rm samp}$

The sampling and query costs used in this paper are summarized in the table
above. We briefly note that a given cost might represent different values in
different contexts, e.g. we shall use the notation $\mathcal{T}_{\rm samp}$ in
both classical and quantum contexts, and it is implied that classical and
quantum sampling costs may not be necessarily the same.

##### Access to quantum controlled rotations

Controlled rotations are a central step in the quantum algorithm for Monte
Carlo (Theorem 3.7). The cost of a controlled rotation depends directly on the
number of bits used to specify the angle of rotation [102]. In our
computational model we assume that controlled rotations come with a unit cost.

######  Definition 3.5 (Access to quantum controlled rotations).

We say that we have access to quantum controlled rotations if we have a
quantum circuit $R$ whose application takes constant time and, for all
rational numbers $x\in[0,1]$ defined by a $(1+c_{2})$-bit string in our fixed-
point arithmetic model 3.2, operates as:

$R\ket{x}\ket{0}=\ket{x}\left(\sqrt{1-x}\ket{0}+\sqrt{x}\ket{1}\right).$ (9)

We note that again this definition allows us to neglect terms $O(c_{1}+c_{2})$
in the runtime and to neglect complications arising from the arithmetic
computation of the $\arcsin$. The access to these rotation unitaries leads to
the following fact.

######  Fact 3.6 (Controlled rotations of a function with an interval).

Consider a rational number representation from Definition 3.1 for some
$c_{1},c_{2}\in\mathbb{N}$. Assume access to controlled rotations according to
Definition 3.5. Assume $(V_{h},\mathcal{T}_{h})$-quantum access to a function
$h$ according to Definition 3.4. For any two bit strings
$a,b\in\\{0,1\\}^{n}$, with $0\leq\mathcal{Q}(a)<\mathcal{Q}(b)$, we can
construct a unitary operator $R^{h}_{a,b}$ on
$\mathcal{H}\otimes\mathbb{C}^{2}$, such that, for all $x\in E$,

$\displaystyle
R^{h}_{a,b}|x,0\rangle=\begin{cases}|x\rangle\left(\sqrt{1-\mathcal{Q}(h(x))/\mathcal{Q}(b)}|0\rangle+\sqrt{\mathcal{Q}(h(x))/\mathcal{Q}(b)}|1\rangle\right)\quad&\text{if}~{}\mathcal{Q}(a)\leq\mathcal{Q}(h(x))\leq\mathcal{Q}(b),\\\
|x,0\rangle\quad&\text{otherwise},\end{cases}$

where an application of $R^{h}_{a,b}$ costs $O(\mathcal{T}_{h})$ time.

###### Proof.

The quantum circuits for the division and for checking the interval run in
constant time in the quantum arithmetic model. The quantum circuits allow us
to prepare $\ket{x}\ket{\mathcal{Q}(h(x))/\mathcal{Q}(b)}$ on the interval
using $V_{h}$ two times, where the second register is of size polynomial in
$c_{1}$ and $c_{2}$. Performing a controlled rotation of an ancilla costs
constant time by Definition 3.5. ∎

##### Quantum algorithm for Monte Carlo

Our quantum algorithm requires the computation of several expectation values.
In this work we use the quantum algorithm for Monte Carlo from Montanaro [70],
already adapted to our computational model.

######  Theorem 3.7 (Quantum algorithm for Monte Carlo $\mathtt{QMonteCarlo}$,
[70, Theorem 2.5]).

Let $X$ be a random variable given via quantum sampling access as in
Definition 3.3. Consider a rational number representation from Definition 3.1
for some $c_{1},c_{2}\in\mathbb{N}$. Let $E\subseteq\mathbb{R}^{d}$ be a
finite set and let $\mathcal{H}$ be a Hilbert space with basis
$\\{|x\rangle\\}_{x\in E}$. Consider a function $h:E\to\\{0,1\\}^{n}$ via
quantum access to the controlled rotations as in Fact 3.6, where
$n\in\mathbb{N}$ such that $c_{1}+c_{2}+1=n$, and each access costs time
$\mathcal{T}_{h}$. Assume that the random variable $\mathcal{Q}(h(X))$ has
finite mean $\mu$, and variance upper-bounded by $\sigma^{2}$ for some known
$\sigma>0$. Given $\delta,\epsilon\in(0,1)$, there is a quantum algorithm,
called $\mathtt{QMonteCarlo}(h(X),\epsilon,\delta,\sigma)$, that runs in
$O((\sigma/\epsilon)\log(1/\delta)\log^{3/2}(\sigma/\epsilon)\log\log(\sigma/\epsilon))\times(T\mathcal{T}_{\rm
samp}+\mathcal{T}_{h})$ time and outputs an estimate $\widetilde{\mu}$ such
that
$\operatorname{Pr}\left[|\widetilde{\mu}-\mu]\geq\epsilon\right]\leq\delta$.

The above result will be used to approximate expected values, e.g.
$\mathbb{E}[e_{t,k}(X_{t})e_{t,l}(X_{t})]$ and
$\mathbb{E}[Z_{\widetilde{\tau}_{t+1}}\vec{e}_{t}(X_{t})]$, and was chosen for
its simplicity. It is possible, though, to use other, more complicated quantum
subroutines for Monte Carlo, e.g. [45, 23, 22]. Refs. [23, 22] propose quantum
algorithms for multivariate Monte Carlo estimation, which could be
particularly suitable in our case, since most of our quantities of interest
are vectors and matrices. However, since these more complex and alternative
quantum subroutines for Monte Carlo lead to the same time complexities up to
polylogarithm factors as Theorem 3.7, we decided to use the above result.

###  3.2 Quantum circuits for the stopping times

Recall that Theorem 2.6 allows us to formulate the stochastic optimal stopping
problem with dynamic programming for the optimal stopping times. Having
introduced the quantum computational model, we are now in the position to
construct quantum circuits for various computations related to the dynamic
programming. In particular, we construct a unitary that propagates backwards
the optimal stopping time by one time step according to Eq. (1). In what
follows, given a path $x\in E^{T}$, by $z_{\widetilde{\tau}_{t}(x)}$ we mean
$z_{\widetilde{\tau}_{t}(x)}(x_{\widetilde{\tau}_{t}(x)})$, i.e., the
associated payoff of the $\widetilde{\tau}_{t}(x)$-th time step of $x$.

######  Lemma 3.8 (Quantum circuits for computing the stopping times).

Let $\widetilde{\tau}_{t^{\prime}}$, $t^{\prime}\in\\{0,\dots,T\\}$, be
stopping times defined in Eq. (1). For all $t\in[T]$, given quantum query
access to functions $\\{z_{t^{\prime}}:E\to\mathbb{R}\\}_{t^{\prime}=1}^{T}$
and $\\{e_{t^{\prime},k}:E\to\mathbb{R}\\}_{t^{\prime}\in[T-1],k\in[m]}$ in
time $\mathcal{T}_{z}$ and $\mathcal{T}_{e}$, respectively, and to the
components of real vectors
$\\{\widetilde{\alpha}_{t^{\prime}}\in\mathbb{R}^{m}\\}_{t^{\prime}=t}^{T}$ in
time $O(1)$, the following statements are true (let $e_{0,k}(x_{0}):=1$):

  1. 1.

There is a unitary $W_{t}$ such that, in time
$O(\mathcal{T}_{z}+m\mathcal{T}_{e})$,

$\displaystyle\begin{cases}W_{t}\ket{x}\ket{\widetilde{\tau}_{t+1}(x)}|\mathbf{0}\rangle^{\otimes
3}=\ket{x}\ket{\widetilde{\tau}_{t+1}(x)}\ket{z_{t}(x_{t})}|\widetilde{\alpha}_{t}\cdot\vec{e}_{t}(x_{t})\rangle|\widetilde{\tau}_{t}(x)\rangle&\quad~{}\text{if}~{}t\neq
T,\\\
W_{t}\ket{x}|\mathbf{0}\rangle=|x\rangle|T\rangle&\quad~{}\text{if}~{}t=T.\end{cases}$

  2. 2.

There is a unitary $V^{(k)}_{t}$ such that
$V^{(k)}_{t}|x\rangle|\widetilde{\tau}_{t}(x)\rangle|\mathbf{0}\rangle=|x\rangle|\widetilde{\tau}_{t}(x)\rangle|z_{\widetilde{\tau}_{t}(x)}e_{t-1,k}(x_{t-1})\rangle$,
for $k\in[m]$, in time $O(T\log(T)\mathcal{T}_{z}+\mathcal{T}_{e})$.

  3. 3.

The unitary $C_{t}^{(k)}:=W_{T}^{\dagger}\dots
W_{t+1}^{\dagger}W_{t}^{\dagger}V_{t}^{(k)}W_{t}W_{t+1}\dots W_{T}$ is such
that
$C_{t}^{(k)}|x\rangle|\mathbf{0}\rangle^{\otimes(3(T-t)+2)}=|x\rangle|z_{\widetilde{\tau}_{t}(x)}e_{t-1,k}(x_{t-1})\rangle|\mathbf{0}\rangle^{\otimes(3(T-t)+1)}$,
for $k\in[m]$, in time $O(T(\log(T)\mathcal{T}_{z}+m\mathcal{T}_{e}))$.

###### Proof.

We start with the first statement. The existence of $W_{T}$ is trivial. Assume
$t\in[T-1]$, then, with one query access to function oracle $z_{t}$ and $m$
query accesses to function oracle $e_{t,k}$ for $k\in[m]$, we can perform

$\ket{x}\ket{\bf 0}\ket{\bf 0}\mapsto\ket{x}\ket{z_{t}(x_{t})}\ket{\bf
0}\mapsto\ket{x}\ket{z_{t}(x_{t})}\ket{\vec{e}_{t}(x_{t})}.$

By using access to the $m$ elements of $\widetilde{\alpha}_{t}$, and $O(m)$
multiplications and additions, we can compute the inner product of
$\widetilde{\alpha}_{t}\cdot\vec{e}_{t}(x_{t})$ in superposition over $x_{t}$,
as

$\ket{x}\ket{\vec{e}_{t}(x_{t})}\ket{\bf
0}\mapsto\ket{x}\ket{\vec{e}_{t}(x_{t})}\ket{\widetilde{\alpha}_{t}\cdot\vec{e}_{t}(x_{t})}.$

Comparing between $z_{t}(x_{t})$ and
$\widetilde{\alpha}_{t}\cdot\vec{e}_{t}(x_{t})$ in constant time, we can
compute $\widetilde{\tau}_{t}(x)$ according to Eq. (1), and hence obtain

$\ket{x}\ket{\widetilde{\tau}_{t+1}(x)}\ket{z_{t}(x_{t})}\ket{\widetilde{\alpha}_{t}\cdot\vec{e}_{t}(x_{t})}\ket{\bf
0}\mapsto\ket{x}\ket{\widetilde{\tau}_{t+1}(x)}\ket{z_{t}(x_{t})}\ket{\widetilde{\alpha}_{t}\cdot\vec{e}_{t}(x_{t})}\ket{\widetilde{\tau}_{t}(x)}.$

Uncomputing the intermediate steps leads to the desired operation. The total
runtime is
$O(\mathcal{T}_{z}+m\mathcal{T}_{e}+m+m+1)=O(\mathcal{T}_{z}+m\mathcal{T}_{e})$.

Regarding the second statement, we require a few circuits which can be
constructed once as a pre-processing step. First, we prepare the input for the
payoff functions in an ancillary register, where the input depends on the
content of the register $\ket{\widetilde{\tau}_{t}(x)}$. For this step, we
prepare a conditional copy quantum circuit $V_{\rm copy}$ which operates as
$\ket{x}\ket{\widetilde{\tau}_{t}(x)}\ket{\bf
0}\to\ket{x}\ket{\widetilde{\tau}_{t}(x)}|x_{\widetilde{\tau}_{t}(x)}\rangle$,
where the register $|{x_{\widetilde{\tau}_{t}(x)}}\rangle$ stores the
$\widetilde{\tau}_{t}(x)$-th step of the path $x$. This circuit operates in
time given by the size of the registers of at most $O(T\log(T))$. Second, from
the access to the different payoff functions we construct access to the
functions in superposition of the time parameter. By assumption, we are given
quantum circuits $V_{z_{t^{\prime}}}$ for $t^{\prime}\in[T]$. From these
quantum circuits we construct the controlled circuit $V_{\rm
select}:=\sum_{t^{\prime}=1}^{T}\ket{t^{\prime}}\bra{t^{\prime}}\otimes
V_{z_{t^{\prime}}}$, which consists of the controlled versions of the circuits
$V_{z_{t^{\prime}}}$ and has a runtime of $O(T\log(T)\mathcal{T}_{z})$ [10].
Now, given $\ket{x}\ket{\widetilde{\tau}_{t}(x)}$, with one application of
$V_{\rm copy}$ and one application of $V_{\rm select}$, we obtain
$z_{\widetilde{\tau}_{t}(x)}$, i.e, the payoff evaluated at
$x_{\widetilde{\tau}_{t}(x)}$, as

$\displaystyle\ket{x}\ket{\widetilde{\tau}_{t}(x)}\ket{\bf 0}\ket{\bf 0}$
$\displaystyle\mapsto$
$\displaystyle\ket{x}\ket{\widetilde{\tau}_{t}(x)}|{x_{\widetilde{\tau}_{t}(x)}}\rangle\ket{\bf
0}\mapsto\ket{x}\ket{\widetilde{\tau}_{t}(x)}|{x_{\widetilde{\tau}_{t}(x)}}\rangle|{z_{\widetilde{\tau}_{t}(x)}}\rangle.$

Using $V_{\rm copy}$ again we uncompute the third register. One query to the
function oracle $e_{t-1,k}$ obtains $e_{t-1,k}(x_{t-1})$ and by multiplication
we obtain,

$\displaystyle\ket{x}\ket{\widetilde{\tau}_{t}(x)}|{z_{\widetilde{\tau}_{t}(x)}}\rangle\ket{e_{t-1,k}(x_{t-1})}\ket{\bf
0}$
$\displaystyle\mapsto\ket{x}\ket{\widetilde{\tau}_{t}(x)}|{z_{\widetilde{\tau}_{t}(x)}}\rangle\ket{e_{t-1,k}(x_{t-1})}|{z_{\widetilde{\tau}_{t}(x)}e_{t-1,k}(x_{t-1})}\rangle.$

We uncompute the third and fourth registers using the given circuits to obtain
the desired operation. The total runtime is
$O(T\log(T)\mathcal{T}_{z}+\mathcal{T}_{e})$.

Finally, for the third statement, it is not hard to see that (let
$\widetilde{\tau}_{T}:=T$)

$\displaystyle C_{t}^{(k)}|x\rangle|\mathbf{0}\rangle^{\otimes(3(T-t)+2)}$
$\displaystyle=W_{T}^{\dagger}\dots
W_{t}^{\dagger}V_{t}^{(k)}|x\rangle|T\rangle|\mathbf{0}\rangle\bigotimes_{j=t}^{T-1}\ket{z_{j}(x_{j})}\ket{\widetilde{\alpha}_{j}\cdot\vec{e}_{j}(x_{j})}\ket{\widetilde{\tau}_{j}(x)}$
$\displaystyle=W_{T}^{\dagger}\dots
W_{t}^{\dagger}|x\rangle|T\rangle|{z_{\widetilde{\tau}_{t}(x)}e_{t-1,k}(x_{t-1})}\rangle\bigotimes_{j=t}^{T-1}\ket{z_{j}(x_{j})}\ket{\widetilde{\alpha}_{j}\cdot\vec{e}_{j}(x_{j})}\ket{\widetilde{\tau}_{j}(x)}$
$\displaystyle=|x\rangle|{z_{\widetilde{\tau}_{t}(x)}e_{t-1,k}(x_{t-1})}\rangle|\mathbf{0}\rangle^{\otimes(3(T-t)+1)}.$

From the two previous statements, the runtime of $C_{t}^{(k)}$ is
$2(T-t+1)O(\mathcal{T}_{z}+m\mathcal{T}_{e})+O(T\log(T)\mathcal{T}_{z}+\mathcal{T}_{e})=O(T(\log(T)\mathcal{T}_{z}+m\mathcal{T}_{e}))$.
∎

###  3.3 Quantum least squares Monte Carlo algorithm

We present our quantum LSM algorithm in Algorithm 2. As previously mentioned,
it follows the classical version in Algorithm 1. It computes the expectations
$\mathbb{E}[e_{t,k}(X_{t})e_{t,l}(X_{t})]$ (for the matrices
$\\{A_{t}\\}_{t=1}^{T-1}$),
$\mathbb{E}[Z_{\widetilde{\tau}_{t+1}}\vec{e}_{t}(X_{t})]$ and
$\mathbb{E}[Z_{\widetilde{\tau}_{1}}]$ using Theorem 3.7 instead of drawing
random samples. Recall that by definition
$Z_{\widetilde{\tau}_{t+1}}=z_{\widetilde{\tau}_{t+1}(\mathbf{X}_{t+1})}(X_{\widetilde{\tau}_{t+1}(\mathbf{X}_{t+1})})$,
i.e., both the optimal stopping time and the payoff depend on the path of the
Markov chain.

Algorithm 2 Quantum LSM algorithm for optimal stopping problem

1:Parameters $\delta\in(0,1)$, $\epsilon>0$. Quantum sampling access to Markov
chain $(X_{t})_{t=1}^{T}$ defined on a finite sample space $\Omega$ and with
finite state space $E\subseteq\mathbb{R}^{d}$. Quantum query access to
$\\{z_{t}:E\to\mathbb{R}\\}_{t=1}^{T}$ and linearly independent functions
$\\{e_{t,k}:E\to\mathbb{R}\\}_{k=1}^{m}$ for $t\in[T-1]$. Let
$L:=\max_{t\in[T-1],k\in[m]}\|e_{t,k}\|_{L^{2}(\rho_{t})}$ and
$R:=\max_{t\in[T]}\|z_{t}\|_{u}$.

2:$\delta_{A}:=\delta/(4Tm^{2})$, $\delta_{b}:=\delta/(4Tm)$,
$\epsilon_{A}:=\epsilon/m$ and $\epsilon_{b}:=\epsilon/\sqrt{m}$.

3:Construct quantum access and controlled rotation access to $e_{t,k}e_{t,l}$,
$\forall k,l\in[m],t\in[T-1]$, with quantum query access to $e_{t,j}$, quantum
circuits for multiplication and Fact 3.6.

4:Compute $\\{\widetilde{A}_{t}\\}_{t=1}^{T-1}$ by calling
$\mathtt{QMonteCarlo}(e_{t,k}(X_{t})e_{t,l}(X_{t})$, $\epsilon_{A}$,
$\delta_{A}$, $L^{2}$) for $k,l\in[m]$.

5:Compute the inverses $\\{\widetilde{A}_{t}^{-1}\\}_{t=1}^{T-1}$.

6:Prepare unitary $W_{T}$ s.t.
$W_{T}|x\rangle|\mathbf{0}\rangle=|x\rangle|\widetilde{\tau}_{T}(x)\rangle$,
where $\widetilde{\tau}_{T}(x)=T$ for all $x\in E^{T}$.

7:for $t=T$ to $2$ do

8: if $t\neq T$ then

9: Prepare unitary $W_{t}$ s.t.
$W_{t}\ket{x}\ket{\widetilde{\tau}_{t+1}(x)}|\mathbf{0}\rangle^{\otimes
3}=\ket{x}\ket{\widetilde{\tau}_{t+1}(x)}\ket{z_{t}(x_{t})}|\widetilde{\alpha}_{t}\cdot\vec{e}_{t}(x_{t})\rangle|\widetilde{\tau}_{t}(x)\rangle$
for any $\widetilde{\tau}_{t+1}(x)\in[T]$ (Lemma 3.8).

10: end if

11: Prepare unitaries $\\{V^{(k)}_{t}\\}_{k=1}^{m}$ s.t.
$V^{(k)}_{t}|x\rangle|\widetilde{\tau}_{t}(x)\rangle|\mathbf{0}\rangle=|x\rangle|\widetilde{\tau}_{t}(x)\rangle|z_{\widetilde{\tau}_{t}(x)}e_{t-1,k}(x_{t-1})\rangle$
(Lemma 3.8).

12: Prepare unitary $W_{T}^{\dagger}\dots
W_{t+1}^{\dagger}W_{t}^{\dagger}V_{t}^{(k)}W_{t}W_{t+1}\dots W_{T}$ for
$k\in[m]$ (Lemma 3.8).

13: Construct quantum access to the controlled rotations of the functions
$z_{\widetilde{\tau}_{t}(x)}e_{t-1,k}(x_{t-1})$ (Fact 3.6).

14: Execute $\mathtt{QMonteCarlo}(Z_{\widetilde{\tau}_{t}}e_{t-1,k}(X_{t-1})$,
$\epsilon_{b}$, $\delta_{b}$, $RL$), for all $k\in[m]$, to compute
$\widetilde{b}_{t-1}$.

15: Compute the vector
$\widetilde{\alpha}_{t-1}=\widetilde{A}_{t-1}^{-1}\widetilde{b}_{t-1}$
classically.

16:end for

17:Prepare unitary $W_{1}$ s.t.
$W_{1}\ket{x}\ket{\widetilde{\tau}_{2}(x)}|\mathbf{0}\rangle^{\otimes
3}=\ket{x}\ket{\widetilde{\tau}_{2}(x)}\ket{z_{1}(x)}|\widetilde{\alpha}_{1}\cdot\vec{e}_{1}(x_{1})\rangle|\widetilde{\tau}_{1}(x)\rangle$
for any $\widetilde{\tau}_{2}(x)\in[T]$ (Lemma 3.8).

18:Prepare unitary $V_{1}$ s.t.
$V_{1}|x\rangle|\widetilde{\tau}_{1}(x)\rangle|\mathbf{0}\rangle=|x\rangle|\widetilde{\tau}_{1}(x)\rangle|z_{\widetilde{\tau}_{1}(x)}\rangle$
(Lemma 3.8).

19:Prepare unitary $W_{T}^{\dagger}\dots
W_{2}^{\dagger}W_{1}^{\dagger}V_{1}W_{1}W_{2}\dots W_{T}$ (Lemma 3.8).

20:Construct quantum access to the controlled rotations of the function
$z_{\widetilde{\tau}_{1}}$ (Fact 3.6).

21:Execute $\mathtt{QMonteCarlo}(Z_{\widetilde{\tau}_{1}}$, $\epsilon$,
$\frac{\delta}{2}$, $R$) to compute $\widetilde{Z}_{\widetilde{\tau}_{1}}$.

22:Output
$\widetilde{\mathcal{U}}_{0}:=\max\big{\\{}Z_{0},\widetilde{Z}_{\widetilde{\tau}_{1}}\big{\\}}$.

We explain the three essential steps of the algorithm. First, we compute
$\mathbb{E}[e_{t,k}(X_{t})e_{t,l}(X_{t})]$ to obtain the matrices
$\\{A_{t}\\}_{t=1}^{T-1}$. For this we provide quantum access to the functions
$e_{t,k}(x_{t})e_{t,l}(x_{t})$ from the quantum access to all $e_{t,k}(x_{t})$
and the arithmetic quantum circuit for multiplication. We take the number of
bits such that the numerical error is much smaller than the error obtained
from estimating the random variable. By our arithmetic model, the computation
is performed in constant time. By availability of the quantum controlled
rotations with Fact 3.6, we then satisfy the hypothesis of Theorem 3.7. Next,
for any given time step $t\in\\{2,\dots,T\\}$, we solve the dynamic recursion
$\widetilde{\tau}_{t}(x)=t\mathbf{1}\\{z_{t}(x_{t})\geq\widetilde{\alpha}_{t}\cdot\vec{e}_{t}(x_{t})\\}+\widetilde{\tau}_{t+1}(x)\mathbf{1}\\{z_{t}(x_{t})<\widetilde{\alpha}_{t}\cdot\vec{e}_{t}(x_{t})\\}$
from time $T$ to time $t$ in superposition via unitaries $C_{t}^{(k)}$ built
by a separate subroutine explained in Lemma 3.8. These unitaries provide
$(C_{t}^{(k)},\mathcal{T}_{ze})$-quantum access to the functions
$z_{\widetilde{\tau}_{t}(x)}e_{t-1,k}(x_{t-1})$, where the cost
$\mathcal{T}_{ze}$ is stated in Lemma 3.8. With the corresponding quantum
controlled rotations, this function access can be used with Theorem 3.7 to
estimate the corresponding vector of expected values
$b_{t-1}:=\mathbb{E}[Z_{\widetilde{\tau}_{t}}\vec{e}_{t-1}(X_{t-1})]$ that is
then used to classically compute the approximation $\widetilde{\alpha}_{t-1}$
for the next time step $t-1$. Finally, for $t=1$, with similar steps, we
provide $(C_{1},\mathcal{T}_{z})$-quantum access to the function
$z_{\widetilde{\tau}_{1}(x)}$, where the cost $\mathcal{T}_{z}$ follows again
from Lemma 3.8. This last estimation then provides the approximation
$\widetilde{\mathcal{U}}_{0}$ to the true initial value of the Snell envelope.
We note that the procedure of approximating a matrix $A_{t}$ and a vector
$b_{t}$ entrywise via quantum algorithms for Monte Carlo followed by the
classical computation of
$\widetilde{\alpha}_{t}=\widetilde{A}_{t}^{-1}\widetilde{b}_{t}$ was already
used in [54].

In order to analyze the approximation error and the complexity of Algorithm 2,
we will need the following result from [105, 106] (already modified to our
notation) that bounds the error between the exact continuation values
$\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]$ and their approximation
$\widetilde{\alpha}_{t}\cdot\vec{e}_{t}(X_{t})$ in terms of the error between
the continuation values $\mathbb{E}[Z_{\widetilde{\tau}_{k+1}}|X_{k}]$
evaluated on the approximated stopping times $\widetilde{\tau}_{k+1}$ and
$\widetilde{\alpha}_{k}\cdot\vec{e}_{t}(X_{k})$ for $k\in\\{t,\dots,T-1\\}$.
Recall the image probability measures $\rho_{t}$ in $E\subseteq\mathbb{R}^{d}$
induced by each element $X_{t}$, $t\in\\{0,\dots,T\\}$.

######  Lemma 3.9 ([105, Lemma 2.2]).

For each $t\in\\{0,\dots,T-1\\}$, we have

$\displaystyle\|\widetilde{\alpha}_{t}\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}$
$\displaystyle\leq
2\sum_{k=t}^{T-1}\|\widetilde{\alpha}_{k}\cdot\vec{e}_{k}(X_{k})-\mathbb{E}[Z_{\widetilde{\tau}_{k+1}}|X_{k}]\|_{L^{2}(\rho_{k})},$
$\displaystyle\|\mathbb{E}[Z_{\widetilde{\tau}_{t+1}}|X_{t}]-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}$
$\displaystyle\leq
2\sum_{k=t+1}^{T-1}\|\widetilde{\alpha}_{k}\cdot\vec{e}_{k}(X_{k})-\mathbb{E}[Z_{\widetilde{\tau}_{k+1}}|X_{k}]\|_{L^{2}(\rho_{k})},$

where
$\widetilde{\alpha}_{0}\cdot\vec{e}_{0}(X_{0}):=\widetilde{Z}_{\widetilde{\tau}_{1}}$
approximates $\mathbb{E}[Z_{\widetilde{\tau}_{1}}]$.

We will also need the following technical result on the sensitivity of square
systems.

######  Theorem 3.10 ([41, Theorem 2.6.2]).

Let $Ax=b$ and $\widetilde{A}\widetilde{x}=\widetilde{b}$, where
$A,\widetilde{A}\in\mathbb{R}^{d\times d}$ and
$b,\widetilde{b}\in\mathbb{R}^{d}$, with $b\neq 0$. Suppose that
$\|A-\widetilde{A}\|_{2}\leq\epsilon_{A}$ and
$\|b-\widetilde{b}\|_{2}\leq\epsilon_{b}$. If
$\epsilon_{A}\leq\sigma_{\min}(A)/2$, where $\sigma_{\min}(A)$ is the minimum
singular value of $A$, then

$\displaystyle\|x-\widetilde{x}\|_{2}\leq\frac{2}{\sigma_{\min}(A)}\left(\frac{\epsilon_{A}\|b\|_{2}}{\sigma_{\min}(A)}+\epsilon_{b}\right).$

We are now able to state a central theorem for our quantum LSM algorithm.
Recall the definition $\|f\|_{u}=\sup\\{|f(x)|:x\in E\\}$ of the uniform norm.

######  Theorem 3.11.

Within the setting of Algorithm 2 with input parameters $\delta$ and
$\epsilon$, let $T\mathcal{T}_{\rm samp}$ be the sampling cost of the Markov
chain and consider a set of linearly independent functions
$\\{e_{t,k}:E\to\mathbb{R}\\}_{k=1}^{m}$ for each $t\in[T-1]$ with
$L:=\max_{t\in[T-1],k\in[m]}\|e_{t,k}\|_{L^{2}(\rho_{t})}$ and query cost
$\mathcal{T}_{e}$. Also consider $\\{z_{t}:E\to\mathbb{R}\\}_{t=0}^{T}$ with
$R:=\max_{t\in[T]}\|z_{t}\|_{u}<\infty$ and query cost $\mathcal{T}_{z}$.
Moreover, let $\sigma_{\min}:=\min_{t\in[T-1]}\sigma_{\min}(A_{t})>0$. Assume
that $\sqrt{m}RL/\sigma_{\min}\geq 1$ and define $\mathcal{T}_{\rm
total}:=\mathcal{T}_{\rm samp}+\mathcal{T}_{z}+\mathcal{T}_{e}$. Then, for any
$\delta\in(0,1)$ and $\epsilon\in(0,\sigma_{\min}/2]$, Algorithm 2 outputs
$\widetilde{\mathcal{U}}_{0}$ such that

$\displaystyle\operatorname{Pr}\left[|\widetilde{\mathcal{U}}_{0}-\mathbb{E}[Z_{\tau_{0}}]|\geq\frac{8T\epsilon
mRL^{2}}{\sigma_{\min}^{2}}+2\sum_{t=1}^{T-1}\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\widetilde{\tau}_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}\right]\leq\delta$
(10)

in time

$\displaystyle O\left(\frac{T^{2}m^{3}}{\epsilon}\mathcal{T}_{\rm
total}L(L+R)\log(T)\log(Tm^{2}/\delta)\log^{3/2}(mL(L+R)/\epsilon)\log\log(mL(L+R)/\epsilon)\right).$

###### Proof.

Set $b_{t}:=\mathbb{E}[Z_{\widetilde{\tau}_{t+1}}\vec{e}_{t}(X_{t})]$. Recall
that $\epsilon_{A}:=\epsilon/m$, $\epsilon_{b}:=\epsilon/\sqrt{m}$,
$\delta_{A}:=\delta/(4Tm^{2})$ and $\delta_{b}:=\delta/(4Tm)$. We start by
computing the complexity of the algorithm. We first note the bounds
$\|b_{t}\|_{2}\leq R\|\mathbb{E}[\vec{e}_{t}(X_{t})]\|_{2}\leq\sqrt{m}RL$ and
$\|A_{t}\|_{2}\leq m\max_{k,l}|\mathbb{E}[e_{t,k}(X_{t})e_{t,l}(X_{t})]|\leq
mL^{2}$ for all $t\in[T-1]$. The computation of all the entries of the
matrices $\\{A_{t}\\}_{t=1}^{T-1}$ requires time

$\displaystyle O\left(\frac{Tm^{2}}{\epsilon_{A}}L^{2}(T\mathcal{T}_{\rm
samp}+\mathcal{T}_{e})\log(1/\delta_{A})\log^{3/2}\left(\frac{L^{2}}{\epsilon_{A}}\right)\log\log\left(\frac{L^{2}}{\epsilon_{A}}\right)\right),$

by calling
$\mathtt{QMonteCarlo}(e_{t,k}(X_{t})e_{t,l}(X_{t}),\epsilon_{A},\delta_{A},L^{2})$
from Theorem 3.7 for $t\in[T-1]$ and $k,l\in[m]$. The computation of all
$b_{t}$ uses time

$\displaystyle O\left(\frac{Tm}{\epsilon_{b}}RL(T\mathcal{T}_{\rm
samp}+T(\log(T)\mathcal{T}_{z}+m\mathcal{T}_{e}))\log(1/\delta_{b})\log^{3/2}\left(\frac{RL}{\epsilon_{b}}\right)\log\log\left(\frac{RL}{\epsilon_{b}}\right)\right),$

by calling
$\mathtt{QMonteCarlo}(Z_{\widetilde{\tau}_{t}}e_{t-1,k}(X_{t-1}),\epsilon_{b},\delta_{b},RL)$
from Theorem 3.7 for $t\in\\{2,\dots,T\\}$ and $k\in[m]$. Note that the term
$T(\log(T)\mathcal{T}_{z}+m\mathcal{T}_{e})$ comes from using the unitaries
$C_{t}^{(k)}$ in $\mathtt{QMonteCarlo}$, each with cost
$O(T(\log(T)\mathcal{T}_{z}+m\mathcal{T}_{e}))$ according to Lemma 3.8.
Computing $\mathbb{E}[Z_{\widetilde{\tau}_{1}}]$ requires time

$\displaystyle O\left(\frac{R}{\epsilon}(T\mathcal{T}_{\rm
samp}+T(\log(T)\mathcal{T}_{z}+m\mathcal{T}_{e}))\log(1/\delta)\log^{3/2}\left(\frac{R}{\epsilon}\right)\log\log\left(\frac{R}{\epsilon}\right)\right),$

by calling
$\mathtt{QMonteCarlo}(Z_{\widetilde{\tau}_{1}},\epsilon,\frac{\delta}{2},R)$
from Theorem 3.7 and where the term
$T(\log(T)\mathcal{T}_{z}+m\mathcal{T}_{e})$ again comes from the unitaries
$C_{t}^{(k)}$ in $\mathtt{QMonteCarlo}$. The classical computation of
$\\{\widetilde{A}^{-1}_{t}\\}_{t=1}^{T-1}$ and
$\\{\widetilde{\alpha}_{t}=\widetilde{A}^{-1}_{t}\widetilde{b}_{t}\\}_{t=1}^{T-1}$
requires time $O(Tm^{\omega_{\ast}})$, where $2\leq\omega_{\ast}<3$. Hence, by
keeping the largest terms of each complexity, the final complexity is upper-
bounded by

$\displaystyle O\left(\frac{T^{2}m^{3}}{\epsilon}\mathcal{T}_{\rm
total}L(L+R)\log(T)\log(Tm^{2}/\delta)\log^{3/2}(mL(L+R)/\epsilon)\log\log(mL(L+R)/\epsilon)\right).$

We now move to the error analysis. Fix $t\in[T-1]$. We start by bounding the
error $\|\widetilde{\alpha}_{t}-\alpha_{t}\|_{2}$ between
$\alpha_{t}=A_{t}^{-1}b_{t}$ and
$\widetilde{\alpha}_{t}=\widetilde{A}_{t}^{-1}\widetilde{b}_{t}$. By using
$\mathtt{QMonteCarlo}$ from Theorem 3.7 we approximate each entry of $A_{t}$
and $b_{t}$ as $|(A_{t})_{jl}-(\widetilde{A}_{t})_{jl}|\leq\epsilon/m$ and
$|(b_{t})_{j}-(\widetilde{b}_{t})_{j}|\leq\epsilon/\sqrt{m}$ for all
$j,l\in[m]$. All approximations hold with probability at least $1-\delta/2T$
by the union bound. This means that, with probability at least $1-\delta/2T$,

$\displaystyle\|A_{t}-\widetilde{A}_{t}\|_{2}\leq\sqrt{\sum_{j,l=1}^{m}|(A_{t})_{jl}-(\widetilde{A}_{t})_{jl}|^{2}}\leq\epsilon,\qquad\|b_{t}-\widetilde{b}_{t}\|_{2}=\sqrt{\sum_{j=1}^{m}|(b_{t})_{j}-(\widetilde{b}_{t})_{j}|^{2}}\leq\epsilon.$

According to Theorem 3.10, we obtain

$\displaystyle\|\widetilde{\alpha}_{t}-\alpha_{t}\|_{2}\leq\frac{2\epsilon}{\sigma_{\min}(A_{t})}\left(1+\frac{\|b_{t}\|_{2}}{\sigma_{\min}(A_{t})}\right)\leq\frac{4\epsilon\sqrt{m}RL}{\sigma_{\min}^{2}}$

with probability at least $1-\delta/2T$, using that
$\sqrt{m}RL/\sigma_{\min}\geq 1$. This, in turn, implies that

$\displaystyle\|\widetilde{\alpha}_{t}\cdot\vec{e}_{t}(X_{t})-\alpha_{t}\cdot\vec{e}_{t}(X_{t})\|_{L^{2}(\rho_{t})}\leq\|\widetilde{\alpha}_{t}-\alpha_{t}\|_{2}\big{\|}\sqrt{\vec{e}_{t}(X_{t})\cdot\vec{e}_{t}(X_{t})}\big{\|}_{L^{2}(\rho_{t})}\leq\frac{4\epsilon
mRL^{2}}{\sigma_{\min}^{2}},$

using that
$\big{\|}\sqrt{\vec{e}_{t}(X_{t})\cdot\vec{e}_{t}(X_{t})}\big{\|}_{L^{2}(\rho_{t})}=\sqrt{\sum_{x\in
E}\rho_{t}(x)\sum_{k=1}^{m}e_{t,k}^{2}(x)}\leq\sqrt{m}L$. Next, we bound

$\displaystyle\|\widetilde{\alpha}_{t}\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\widetilde{\tau}_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}$
$\displaystyle\leq\|\widetilde{\alpha}_{t}\cdot\vec{e}_{t}(X_{t})-\alpha_{t}\cdot\vec{e}_{t}(X_{t})\|_{L^{2}(\rho_{t})}+\|\alpha_{t}\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\widetilde{\tau}_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}$
$\displaystyle\leq\frac{4\epsilon
mRL^{2}}{\sigma_{\min}^{2}}+\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\widetilde{\tau}_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})},$
(11)

using that
$\alpha_{t}=\arg\operatorname*{\min}_{a\in\mathbb{R}^{m}}\mathbb{E}\big{[}(Z_{\widetilde{\tau}_{t+1}}-a\cdot\vec{e}_{t}(X_{t}))^{2}\big{]}$
minimizes the least square estimator.

Finally, by the union bound, with probability at least $1-\delta$, Eq. (11)
holds for all $t\in[T-1]$, together with
$|\widetilde{Z}_{\widetilde{\tau}_{1}}-\mathbb{E}[Z_{\widetilde{\tau}_{1}}]|\leq\epsilon$.
Lemma 3.9 then leads to

$\displaystyle|\widetilde{Z}_{\widetilde{\tau}_{1}}-\mathbb{E}[Z_{\tau_{1}}]|$
$\displaystyle\leq\frac{8T\epsilon
mRL^{2}}{\sigma_{\min}^{2}}+2\sum_{t=1}^{T-1}\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\widetilde{\tau}_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})},$

which implies Eq. (10) by using
$|{\max}\\{a_{0},a_{1}\\}-\max\\{a_{0},a_{2}\\}|\leq|a_{1}-a_{2}|$ with
$a_{0},a_{1},a_{2}\in\mathbb{R}$ on the definition of
$\widetilde{\mathcal{U}}_{0}$ in Eq. (3). ∎

######  Remark 3.12.

As previously mentioned, it is possible to use the multivariable quantum
algorithm for Monte Carlo from [23, 22] within Algorithm 2. This leads to the
same complexity of the previous theorem,
$\widetilde{O}\big{(}\frac{T^{2}m^{3}}{\epsilon}\mathcal{T}_{\rm
total}L(L+R)\big{)}$, up to $\operatorname{polylog}$ factors.

As previously mentioned, the quantities
$\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\widetilde{\tau}_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}$
appearing in the theorem above are sometimes known as _approximation errors_
of the approximation architecture $\mathscr{H}_{t}$, here chosen to be finite-
dimensional, linear and generated by a set of linearly independent functions.
We could rewrite these terms alternatively as
$\inf_{f\in\mathscr{H}_{t}}\|f(X_{t})-\mathbb{E}[Z_{\widetilde{\tau}_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}$.

Note that the approximation errors
$\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\widetilde{\tau}_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}$
appearing in Theorem 3.11 depend on the approximated stopping times
$\widetilde{\tau}_{t+1}$, which in turn depend on $\mathscr{H}_{t^{\prime}}$
for $t^{\prime}>t$. It is possible to restate Theorem 3.11 in terms of
$\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}$,
which we do in the next theorem by following a similar approach to [105,
Theorem 6.1]. The downside is that the time dependence now becomes
exponential.

######  Theorem 3.13.

Within the setting of Algorithm 2 with input parameters $\delta$ and
$\epsilon$, let $T\mathcal{T}_{\rm samp}$ be the sampling cost of the Markov
chain and consider a set of linearly independent functions
$\\{e_{t,k}:E\to\mathbb{R}\\}_{k=1}^{m}$ for each $t\in[T-1]$ with
$L:=\max_{t\in[T-1],k\in[m]}\|e_{t,k}\|_{L^{2}(\rho_{t})}$ and query cost
$\mathcal{T}_{e}$. Also consider $\\{z_{t}:E\to\mathbb{R}\\}_{t=0}^{T}$ with
$R:=\max_{t\in[T]}\|z_{t}\|_{u}<\infty$ and query cost $\mathcal{T}_{z}$.
Moreover, let $\sigma_{\min}:=\min_{t\in[T-1]}\sigma_{\min}(A_{t})>0$. Assume
that $\sqrt{m}RL/\sigma_{\min}\geq 1$ and define $\mathcal{T}_{\rm
total}:=\mathcal{T}_{\rm samp}+\mathcal{T}_{z}+\mathcal{T}_{e}$. Then, for any
$\delta\in(0,1)$ and $\epsilon\in(0,\sigma_{\min}/2]$, Algorithm 2 outputs
$\widetilde{\mathcal{U}}_{0}$ such that

$\displaystyle\operatorname{Pr}\left[\big{|}\widetilde{\mathcal{U}}_{0}-\mathbb{E}[Z_{\tau_{0}}]\big{|}\geq
5^{T}\left(\frac{4\epsilon
mRL^{2}}{\sigma_{\min}^{2}}+\operatorname*{\max}_{0<t<T}\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}\right)\right]\leq\delta$

in time

$\displaystyle O\left(\frac{T^{2}m^{3}}{\epsilon}\mathcal{T}_{\rm
total}L(L+R)\log(T)\log(Tm^{2}/\delta)\log^{3/2}(mL(L+R)/\epsilon)\log\log(mL(L+R)/\epsilon)\right).$

###### Proof.

The proof follows the same steps of the proof of Theorem 3.11, with the
further observation in Eq. (11) that if

$\displaystyle\|\widetilde{\alpha}_{\ell}\cdot\vec{e}_{\ell}(X_{\ell})-\mathbb{E}[Z_{\widetilde{\tau}_{\ell+1}}|X_{\ell}]\|_{L^{2}(\rho_{\ell})}\leq\epsilon_{0}+\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{\ell}(X_{\ell})-\mathbb{E}[Z_{\widetilde{\tau}_{\ell+1}}|X_{\ell}]\|_{L^{2}(\rho_{\ell})},$
(12)

for all $\ell\in\\{t,\dots,T-1\\}$, where we defined
$\epsilon_{0}:=\frac{4\epsilon mRL^{2}}{\sigma_{\min}^{2}}$, then

$\displaystyle
2(T-\ell)\epsilon_{0}+2\sum_{k=\ell}^{T-1}\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{k}(X_{k})-\mathbb{E}[Z_{\widetilde{\tau}_{k+1}}|X_{k}]\|_{L^{2}(\rho_{k})}\leq
5^{T-\ell}(\epsilon_{0}+M^{\ast}_{\ell}),$ (13)

for all $\ell\in\\{t,\dots,T-1\\}$, where
$M^{\ast}_{\ell}:=\max_{k=\ell,\dots,T-1}\big{(}{\operatorname*{\min}}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{k}(X_{k})-\mathbb{E}[Z_{\tau_{k+1}}|X_{k}]\|_{L^{2}(\rho_{k})}\big{)}$.
We prove this bound using backward induction as follows. Eq. (13) clearly
holds for $\ell=T-1$. Assume it holds for $\ell=t+1$, we shall prove it is
also true for $\ell=t$. First notice that, by the triangle inequality followed
by Lemma 3.9 and then Eq. (12),

$\displaystyle\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\widetilde{\tau}_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}$
$\displaystyle\leq\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}+\|\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]-\mathbb{E}[Z_{\widetilde{\tau}_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}$
$\displaystyle\leq\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}+2\sum_{k=t+1}^{T-1}\|\widetilde{\alpha}_{k}\cdot\vec{e}_{k}(X_{k})-\mathbb{E}[Z_{\widetilde{\tau}_{k+1}}|X_{k}]\|_{L^{2}(\rho_{k})}$
$\displaystyle\leq\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}+2\sum_{k=t+1}^{T-1}\big{(}\epsilon_{0}+\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{k}(X_{k})-\mathbb{E}[Z_{\widetilde{\tau}_{k+1}}|X_{k}]\|_{L^{2}(\rho_{k})}\big{)}$
$\displaystyle\leq
M^{\ast}_{t}+2(T-t-1)\epsilon_{0}+2\sum_{k=t+1}^{T-1}\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{k}(X_{k})-\mathbb{E}[Z_{\widetilde{\tau}_{k+1}}|X_{k}]\|_{L^{2}(\rho_{k})}.$

Using the above inequality followed by the induction hypothesis,

$\displaystyle 2(T-t)$
$\displaystyle\epsilon_{0}+2\sum_{k=t}^{T-1}\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{k}(X_{k})-\mathbb{E}[Z_{\widetilde{\tau}_{k+1}}|X_{k}]\|_{L^{2}(\rho_{k})}$
$\displaystyle\leq
2\epsilon_{0}+2M^{\ast}_{t}+6(T-t-1)\epsilon_{0}+6\sum_{k=t+1}^{T-1}\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{k}(X_{k})-\mathbb{E}[Z_{\widetilde{\tau}_{k+1}}|X_{k}]\|_{L^{2}(\rho_{k})}$
$\displaystyle\leq 2\epsilon_{0}+2M^{\ast}_{t}+3\cdot
5^{T-t-1}(\epsilon_{0}+M^{\ast}_{t+1})$ $\displaystyle\leq
5^{T-t}(\epsilon_{0}+M^{\ast}_{t}),$

proving the induction statement. The theorem follows by taking $\ell=0$ in Eq.
(13) and using
$|{\max}\\{a_{0},a_{1}\\}-\max\\{a_{0},a_{2}\\}|\leq|a_{1}-a_{2}|$,
$a_{0},a_{1},a_{2}\in\mathbb{R}$, on the definition of
$\widetilde{\mathcal{U}}_{0}$ in Eq. (3). ∎

Given the above theorems, we prove a classical analogue.

######  Theorem 3.14.

Within the setting of Algorithm 1, consider $N$ independent sample paths with
sample cost $\mathcal{T}_{\rm samp}$ and let the linearly independent
functions $\\{e_{t,k}:E\to\mathbb{R}\\}_{k=1}^{m}$ for each $t\in[T-1]$ be
such that $L:=\max_{t\in[T-1],k\in[m]}\|e_{t,k}\|_{L^{2}(\rho_{t})}$ and have
query cost $\mathcal{T}_{e}$. Also consider
$\\{z_{t}:E\to\mathbb{R}\\}_{t=0}^{T}$ with
$R:=\max_{t\in[T]}\|z_{t}\|_{u}<\infty$ and query cost $\mathcal{T}_{z}$.
Moreover, let $\sigma_{\min}:=\min_{t\in[T-1]}\sigma_{\min}(A_{t})>0$ and
assume that $\sqrt{m}RL/\sigma_{\min}\geq 1$. Then, for any
$\epsilon\in(0,\sigma_{\min}/2]$, Algorithm 1 runs for
$O(Tm^{2}N+Tm^{\omega_{\ast}}+TN(\mathcal{T}_{\rm
samp}+\mathcal{T}_{z}+m\mathcal{T}_{e}))$ time and returns an estimate
$\widetilde{\mathcal{U}}_{0}$ such that

$\displaystyle\operatorname{Pr}\left[\big{|}\widetilde{\mathcal{U}}_{0}-\mathbb{E}[Z_{\tau_{0}}]\big{|}\geq
5^{T}\left(\frac{4\epsilon
mRL^{2}}{\sigma_{\min}^{2}}+\operatorname*{\max}_{0<t<T}\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\tau_{t+1}}|X_{k}]\|_{L^{2}(\rho_{t})}\right)\right]\leq
6m^{2}e^{-2N\epsilon^{2}/m^{2}}.$

###### Proof.

The error analysis is very similar to that of Theorems 3.11 and 3.13,
therefore we shall just point out the required changes. Each entry of $A_{t}$
and $b_{t}$ is approximated using a Chernoff bound, i.e.,
$\operatorname{Pr}[|(A_{t})_{jl}-(\widetilde{A}_{t})_{jl}|\geq\epsilon/m]\leq
2e^{-2N\epsilon^{2}/m^{2}}$ and
$\operatorname{Pr}[|(b_{t})_{j}-(\widetilde{b}_{t})_{j}|\geq\epsilon/\sqrt{m}]\leq
2e^{-2N\epsilon^{2}/m}$ for all $j,l\in[m]$. Moreover,
$\operatorname{Pr}[|\widetilde{Z}_{\widetilde{\tau}_{1}}-\mathbb{E}[Z_{\widetilde{\tau}_{1}}]|\geq\epsilon]\leq
2e^{-2N\epsilon^{2}}$. Therefore, by the union bound, all approximations hold
with probability at least
$1-2m^{2}e^{-2N\epsilon^{2}/m^{2}}-2me^{-2N\epsilon^{2}/m}-2e^{-2N\epsilon^{2}}\geq
1-6m^{2}e^{-2N\epsilon^{2}/m^{2}}$.

Regarding the time complexity, the most expensive computational steps are
calculating the matrices $\widetilde{A}_{t}$, which requires $O(Tm^{2}N)$
time, and inverting them and computing the vectors $\widetilde{\alpha}_{t}$,
which requires $O(Tm^{\omega_{\ast}})$ time. Sampling
$(X_{t}^{(1)},\dots,X_{t}^{(N)})_{t=0}^{T}$ and querying
$(Z_{t}^{(1)},\dots,Z_{t}^{(N)})_{t=0}^{T}$ and
$(e_{t,k}(X_{t}^{(1)}),\dots,e_{t,k}(X_{t}^{(N)}))_{t\in[T-1],k\in[m]}$
require $O(NT(\mathcal{T}_{\rm samp}+\mathcal{T}_{z}+m\mathcal{T}_{e}))$ time.
All the other steps, computing
$\frac{1}{N}\sum_{n=1}^{N}Z^{(n)}_{\widetilde{\tau}^{(n)}_{t+1}}\vec{e}(X^{(n)}_{t})$,
$\widetilde{\mathcal{U}}_{0}$ and
$\widetilde{\alpha}_{t}\cdot\vec{e}_{t}(X_{t}^{(n)})$ require $O(TmN)$ or
$O(Tm)$ time. ∎

If $T_{\rm samp},T_{z},T_{e}=O(1)$, then the complexity is
$O(Tm^{2}N+Tm^{\omega_{\ast}})$. The factor $Tm^{2}$ comes from computing
$\\{A_{t}\\}_{t=1}^{T-1}$ and accounts for runtime instead of only number of
samples.

We summarize and compare the results from Theorems 3.14 and 3.13 into a single
corollary.

######  Corollary 3.15.

Within the setting of Algorithm 1 with input parameters $\delta,N$ and
Algorithm 2 with input parameters $\delta,\epsilon_{0}$, let
$T\mathcal{T}_{\rm samp}$ be the sampling cost of the Markov chain and
consider a set of linearly independent functions
$\\{e_{t,k}:E\to\mathbb{R}\\}_{k=1}^{m}$ for $t\in[T-1]$ with query cost
$\mathcal{T}_{e}$ and
$L:=\max_{t\in[T-1],k\in[m]}\|e_{t,k}\|_{L^{2}(\rho_{t})}$. Also consider
$\\{z_{t}:E\to\mathbb{R}\\}_{t=0}^{T}$ with query cost $\mathcal{T}_{z}$ and
$R:=\max_{t\in[T]}\|z_{t}\|_{u}<\infty$. Let
$\sigma_{\min}\leq\min_{t\in[T-1]}\sigma_{\min}(A_{t})$ be known. Assume that
$\sqrt{m}RL/\sigma_{\min}\geq 1$ and define $\mathcal{T}_{\rm
total}:=\mathcal{T}_{\rm samp}+\mathcal{T}_{z}+\mathcal{T}_{e}$. For any
$\delta\in(0,1)$ and
$\epsilon\in\big{(}0,\frac{2mRL^{2}}{\sigma_{\min}}\big{]}$, if
$\epsilon_{0}=\frac{\epsilon\sigma_{\min}^{2}}{4mRL^{2}}$ and
$N=\big{\lceil}\frac{m^{2}}{2\epsilon_{0}^{2}}\log(6m^{2}/\delta)\big{\rceil}$,
then Algorithms 1 and 2 output $\widetilde{\mathcal{U}}_{0}$ such that

$\displaystyle\operatorname{Pr}\left[|\widetilde{\mathcal{U}}_{0}-\mathbb{E}[Z_{\tau_{0}}]|\geq
5^{T}\left(\epsilon+\operatorname*{\max}_{0<t<T}\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}\right)\right]\leq\delta$

using time, respectively,

$\displaystyle
O\left(\frac{Tm^{6}}{\epsilon^{2}}\frac{R^{2}L^{4}}{\sigma_{\min}^{4}}\mathcal{T}_{\rm
total}\log(m^{2}/\delta)\right)$

and (up to $\log\log$ factors)

$\displaystyle
O\left(\frac{T^{2}m^{4}}{\epsilon}\frac{RL^{3}(L+R)}{\sigma_{\min}^{2}}\mathcal{T}_{\rm
total}\log(T)\log\left(\frac{Tm^{2}}{\delta}\right)\log^{3/2}\left(\frac{m^{2}RL^{3}(L+R)}{\epsilon\sigma_{\min}^{2}}\right)\right).$

###### Proof.

The results concerning the quantum Algorithm 2 were already proven in Theorem
3.13, we just use that
$\epsilon_{0}=\frac{\epsilon\sigma_{\min}^{2}}{4mRL^{2}}$. It is left to prove
the results about the classical Algorithm 1. Indeed, by setting
$N=\big{\lceil}\frac{m^{2}}{2\epsilon_{0}^{2}}\log(6m^{2}/\delta)\big{\rceil}$
into Theorem 3.14 with
$\epsilon_{0}=\frac{\epsilon\sigma_{\min}^{2}}{4mRL^{2}}$, we obtain

$\displaystyle\operatorname{Pr}\left[|\widetilde{\mathcal{U}}_{0}-\mathbb{E}[Z_{\tau_{0}}]|\geq
5^{T}\left(\epsilon+\operatorname*{\max}_{0<t<T}\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}\right)\right]\leq\delta$

in time $O(Tm^{2}N+Tm^{\omega_{\ast}}+TN(\mathcal{T}_{\rm
samp}+\mathcal{T}_{z}+m\mathcal{T}_{e}))=O\left(\frac{Tm^{4}}{\epsilon^{2}}(m^{2}+\mathcal{T}_{\rm
samp}+\mathcal{T}_{z}+m\mathcal{T}_{e})\frac{R^{2}L^{4}}{\sigma_{\min}^{4}}\log(\frac{m^{2}}{\delta})\right)$
$=O\left(\frac{Tm^{6}}{\epsilon^{2}}\frac{R^{2}L^{4}}{\sigma_{\min}^{4}}\mathcal{T}_{\rm
total}\log(\frac{m^{2}}{\delta})\right)$. ∎

######  Remark 3.16.

We note that in the previous theorem it is necessary to know a lower bound on
$\min_{t\in[T-1]}\sigma_{\min}(A_{t})$ since it must be inputted into quantum
subroutines for Monte Carlo to bound the variance. It is possible, though, to
get around this detail by using the more complicated algorithm from [45] which
does not require a bound on the variance.

The accuracy in Corollary 3.15 is made up of two parts: the “sampling” or
stochastic error $\epsilon$, and the approximation architecture error
$\max_{0<t<T}\min_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}$
(ignoring the time dependence). If the sample and query complexities
$\mathcal{T}_{\rm samp},\mathcal{T}_{z}$ and $\mathcal{T}_{e}$ are constant,
then we see that the classical and quantum algorithms obtain the same additive
accuracy in approximating $\mathcal{U}_{0}$ in time
$\widetilde{O}(Tm^{6}/\epsilon^{2})$ and $\widetilde{O}(T^{2}m^{4}/\epsilon)$,
respectively, thus constituting a quadratic improvement in $\epsilon$. The
sampling error can be made arbitrary small by increasing the number of path
samples or calls to quantum inputs. The approximation architecture error, on
the other hand, does not depend on $\epsilon$, and is intrinsic to projecting
the continuation values onto a lower-dimensional space generated by
predetermined measurable functions.

##  4 Linear finite-dimensional spaces of polynomials

In this section we shall focus specifically on linear finite-dimensional
spaces of polynomials, in order to bound the approximation errors
$\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}$
from Corollary 3.15 and thus obtain overall error estimates. The bounds will
come from assuming smoothness properties on the continuation values, which is
not an unusual assumption in finance: the continuation values
$\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]$ can even be $C^{\infty}$-smooth in some
models, e.g. Black-Scholes [37].

Consider an open, connected and bounded set $D\subset\mathbb{R}^{d}$, together
with its closure $\overline{D}$. Let $\mathcal{R}_{q}$ be the space of all
polynomials of degree at most $q$ with real coefficients, i.e., the space of
all linear combinations of terms $x_{1}^{a_{1}}x_{2}^{a_{2}}\dots
x_{d}^{a_{d}}$ with integers $a_{i}\geq 0$ such that $\sum_{i=1}^{d}a_{i}\leq
q$. Let $\mathcal{R}_{q}(D)$ be the space of all polynomials in
$\mathcal{R}_{q}$ that are restrictions to $D$, i.e., with domain in $D$.
Moreover, for any non-negative integer $n$, let $C^{n}(\overline{D})$ be the
space of all continuous, real-valued functions $f$ on $\overline{D}$ whose
continuous derivatives $\partial^{a}f$ exist on $D$ for all multi-indices $a$
such that $\|a\|_{1}=\sum_{i=1}^{d}a_{i}\leq n$ and that have continuous
extensions onto $\overline{D}$. Define the norm on this space as

$\displaystyle\|f\|_{C^{n}(\overline{D})}:=\sum_{\|a\|_{1}\leq
n}\operatorname*{\sup}_{x\in D}|\partial^{a}f(x)|.$

If $f\in C^{n}(\overline{D})$, then the next Jackson-type inequality allows
one to bound the error committed by approximating $f$ by a polynomial in
$\mathcal{R}_{q}(D)$.

######  Lemma 4.1 ([28, Chapter 7, Theorem 6.2]).

Let $D\subset\mathbb{R}^{d}$ be the interior of some closed cube. If $q>n\geq
1$, then, for any $f\in C^{n}(\overline{D})$,

$\displaystyle\operatorname*{\inf}_{p\in\mathcal{R}_{q}}\operatorname*{\sup}_{x\in
D}|f(x)-p(x)|\leq C_{D}q^{-n}\|f\|_{C^{n}(\overline{D})},$

where the constant $C_{D}$ depends on $D$ but not on $q$ or $f$.

Another smoothness condition we can consider is that of Lipschitz continuous
functions. Let $Q_{d}(\lambda)$ be the closed cube of side length
$0<2\lambda<\infty$ in a $d$-dimensional Euclidean space centered at the
origin. Let $\operatorname{Lip}(C_{L},Q_{d}(\lambda))$ be the space of
Lipschitz continuous functions $f$ on $Q_{d}(\lambda)$ such that

$\displaystyle|f(x)-f(x^{\prime})|\leq
C_{L}\|x-x^{\prime}\|_{\infty},\quad\forall x,x^{\prime}\in Q_{d}(\lambda),$

with Lipschitz constant $C_{L}$. The Lipschitz condition implies the following
error bound for approximating functions by polynomials in
$\mathcal{R}_{q}(Q_{d}(\lambda))$.

######  Lemma 4.2 ([33, Chapter 9, Theorem 10]).

For all $C_{L}>0$, for all $f\in\operatorname{Lip}(C_{L},Q_{d}(\lambda))$,

$\displaystyle\operatorname*{\inf}_{p\in\mathcal{R}_{q}}\operatorname*{\sup}_{x\in
Q_{d}(\lambda)}|f(x)-p(x)|\leq\frac{88\lambda C_{L}d}{d+q}.$

Given the above theorems, we obtain a corollary of Theorem 3.13 where the
continuation values are assumed to be functions either in
$C^{n}(Q_{d}(\lambda))$ or in $\operatorname{Lip}(C_{L},Q_{d}(\lambda))$.

######  Theorem 4.3.

Let $Q_{d}(\lambda)$ be the closed cube with side length $2\lambda$ in
$\mathbb{R}^{d}$ centered at the origin. Within the setting of Algorithm 2
with parameters $\delta$ and $\epsilon_{0}$, let $T\mathcal{T}_{\rm samp}$ be
the sampling cost of the Markov chain and consider the approximation
architecture $\mathcal{R}_{q}(Q_{d}(\lambda))$ for all $t\in\\{0,\dots,T-1\\}$
with
$L:=\max_{t\in[T]}\sup_{p\in\mathcal{R}_{q}}\|p\|_{L^{2}(Q_{d}(\lambda),\rho_{t})}$
and query cost $\mathcal{T}_{e}$. Also consider
$\\{z_{t}:Q_{d}(\lambda)\to\mathbb{R}\\}_{t=0}^{T}$ with query cost
$\mathcal{T}_{z}$ and $R:=\max_{t\in[T]}\|z_{t}\|_{u}<\infty$. Moreover, let
$\sigma_{\min}\leq\min_{t\in[T-1]}\sigma_{\min}(A_{t})$ be known. Assume that
$\sqrt{m}RL/\sigma_{\min}\geq 1$ and define $\mathcal{T}_{\rm
total}:=\mathcal{T}_{\rm samp}+\mathcal{T}_{z}+\mathcal{T}_{e}$.

_(a)._ Suppose there is a positive integer $n$, $n<q$, such that
$\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\in C^{n}(Q_{d}(\lambda))$ for all
$t\in\\{0,\dots,T-1\\}$. Let $C\geq C_{D}\max_{0\leq
t<T}\|\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{C^{n}(Q_{d}(\lambda))}$ be known.
Then, for $\delta\in(0,1)$,
$\epsilon\in\big{(}0,5^{T}\big{(}\frac{RL^{2}(e^{2n}C)^{d/n}}{\sigma_{\min}}\big{)}^{n/(n+d)}\big{]}$
and $q=\left\lceil(2\cdot 5^{T}C/\epsilon)^{1/n}\right\rceil$, if
$\epsilon_{0}=\frac{\sigma_{\min}^{2}}{4(e^{2n}C)^{d/n}RL^{2}}\big{(}\frac{\epsilon}{2\cdot
5^{T}}\big{)}^{1+d/n}$, Algorithm 2 outputs $\widetilde{\mathcal{U}}_{0}$ such
that

$\displaystyle\operatorname{Pr}\left[|\widetilde{\mathcal{U}}_{0}-\mathbb{E}[Z_{\tau_{0}}]|\leq\epsilon\right]\geq
1-\delta$

in time

$\displaystyle\widetilde{O}\left(T^{2}(e^{2n}C)^{4d/n}\mathcal{T}_{\rm
total}\frac{RL^{3}(R+L)}{\sigma_{\min}^{2}}\left(\frac{5^{T}}{\epsilon}\right)^{1+4d/n}\right),$

where $\widetilde{O}(\cdot)$ hides $\operatorname{polylog}$ terms in
$5^{T},R,L,\sigma_{\min},C,e^{d},\epsilon,\delta$.

_(b)._ Suppose that
$\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\in\operatorname{Lip}(C_{L},Q_{d}(\lambda))$
for all $t\in\\{0,\dots,T-1\\}$, and suppose that $C_{L}$ is known. Then, for
$\delta\in(0,1)$,
$\epsilon\in\big{(}0,5^{T}\big{(}\frac{RL^{2}(176e^{2}\lambda
C_{L}d)^{d}}{\sigma_{\min}}\big{)}^{1/(1+d)}\big{]}$ and $q=\lceil 176\cdot
5^{T}\lambda C_{L}d/\epsilon\rceil$, if
$\epsilon_{0}=\frac{\sigma_{\min}^{2}}{4(176e^{2}\lambda
C_{L}d)^{d}RL^{2}}\big{(}\frac{\epsilon}{2\cdot 5^{T}}\big{)}^{1+d}$,
Algorithm 2 outputs $\widetilde{\mathcal{U}}_{0}$ such that

$\displaystyle\operatorname{Pr}\left[|\widetilde{\mathcal{U}}_{0}-\mathbb{E}[Z_{\tau_{0}}]|\leq\epsilon\right]\geq
1-\delta$

in time

$\displaystyle\widetilde{O}\left(T^{2}\big{(}\lambda
C_{L}d\big{)}^{4d}\mathcal{T}_{\rm
total}\frac{RL^{3}(R+L)}{\sigma_{\min}^{2}}\left(\frac{5^{T}}{\epsilon}\right)^{1+4d}\right),$

where $\widetilde{O}(\cdot)$ hides $\operatorname{polylog}$ terms in
$5^{T},R,L,\sigma_{\min},\lambda C_{L},e^{d},\epsilon$.

###### Proof.

We start with the case when $\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\in
C^{n}(Q_{d}(\lambda))$ for all $t\in\\{0,\dots,T-1\\}$. Since the uniform norm
upper-bounds the $L^{2}$-norm over a probability space, and using Lemma 4.1,
then

$\displaystyle\operatorname*{\inf}_{p\in\mathcal{R}_{q}(Q_{d}(\lambda))}\|p-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}\leq\operatorname*{\inf}_{p\in\mathcal{R}_{q}(Q_{d}(\lambda))}\|p-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{u}\leq
Cq^{-n}.$

Thus the result follows from Theorem 3.13 by noting that $\mathcal{R}_{q}$ has
dimension $m:=\binom{q+d}{d}\leq e^{d}q^{d}$, and so we must take $m$ linearly
independent polynomials in Algorithm 2. By setting $q=\lceil(2\cdot
5^{T}C/\epsilon)^{1/n}\rceil$ and using the accuracy
$\epsilon_{0}=\frac{\sigma_{\min}^{2}}{4(e^{2n}C)^{d/n}RL^{2}}\big{(}\frac{\epsilon}{2\cdot
5^{T}}\big{)}^{1+d/n}$ in Theorem 3.13 (note that
$\epsilon_{0}\in(0,\sigma_{\min}/2]$), we obtain the final error

$\displaystyle
5^{T}\left(\frac{4\epsilon_{0}mRL^{2}}{\sigma_{\min}^{2}}+\frac{C}{q^{n}}\right)=5^{T}\left(\frac{4\epsilon_{0}e^{d}q^{d}RL^{2}}{\sigma_{\min}^{2}}+\frac{C}{q^{n}}\right)\leq(1/2+1/2)\epsilon=\epsilon.$

This in turn implies the complexity

$\displaystyle\widetilde{O}\left(\frac{T^{2}m^{3}}{\epsilon_{0}}L(L+R)\mathcal{T}_{\rm
total}\right)=\widetilde{O}\left(T^{2}(e^{2n}C)^{4d/n}\mathcal{T}_{\rm
total}\frac{RL^{3}(R+L)}{\sigma_{\min}^{2}}\left(\frac{5^{T}}{\epsilon}\right)^{1+4d/n}\right).$

We now prove the case when
$\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\in\operatorname{Lip}(C_{L},Q_{d}(\lambda))$
for all $t\in\\{0,\dots,T-1\\}$. By the same token as before, but now using
Lemma 4.2, we obtain

$\displaystyle\operatorname*{\inf}_{p\in\mathcal{R}_{q}(Q_{d}(\lambda))}\|p-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}\leq\frac{88\lambda
C_{L}d}{q}.$

Once again, we must pick $m\leq e^{d}q^{d}$ linearly independent polynomials
in Algorithm 2. By setting $q=\lceil 176\cdot 5^{T}\lambda
C_{L}d/\epsilon\rceil$ and using the accuracy
$\epsilon_{0}=\frac{\sigma_{\min}^{2}}{4(176e^{2}\lambda
C_{L}d)^{d}RL^{2}}\big{(}\frac{\epsilon}{2\cdot 5^{T}}\big{)}^{1+d}$ in
Theorem 3.13 (note that $\epsilon_{0}\in(0,\sigma_{\min}/2]$), we obtain the
final error

$\displaystyle
5^{T}\left(\frac{4\epsilon_{0}mRL^{2}}{\sigma_{\min}^{2}}+\frac{88\lambda
C_{L}d}{q}\right)=5^{T}\left(\frac{4\epsilon_{0}e^{d}q^{d}RL^{2}}{\sigma_{\min}^{2}}+\frac{88\lambda
C_{L}d}{q}\right)\leq(1/2+1/2)\epsilon=\epsilon.$

This in turn implies the complexity

$\widetilde{O}\left(\frac{T^{2}m^{3}}{\epsilon_{0}}L(L+R)\mathcal{T}_{\rm
total}\right)=\widetilde{O}\left(T^{2}\big{(}\lambda
C_{L}d\big{)}^{4d}\mathcal{T}_{\rm
total}\frac{RL^{3}(R+L)}{\sigma_{\min}^{2}}\left(\frac{5^{T}}{\epsilon}\right)^{1+4d}\right).\qed$

Similar results can be obtained for the classical Algorithm 1 by using
Corollary 3.15. In the case of $\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\in
C^{n}(Q_{d}(\lambda))$, pick $q=\lceil(2\cdot 5^{T}C/\epsilon)^{1/n}\rceil$
and final accuracy $\epsilon_{f}=\frac{\epsilon}{2\cdot 5^{T}}$ in Corollary
3.15, which implies time

$\displaystyle\widetilde{O}\left(\frac{Tm^{6}}{\epsilon_{f}^{2}}\frac{R^{2}L^{4}}{\sigma_{\min}^{4}}\mathcal{T}_{\rm
total}\right)=\widetilde{O}\left(T(e^{n}C)^{6d/n}\frac{R^{2}L^{4}}{\sigma_{\min}^{4}}\left(\frac{5^{T}}{\epsilon}\right)^{2+6d/n}\mathcal{T}_{\rm
total}\right)$

up to the logarithm factor $\log(m^{2}/\delta)$.

In the case of
$\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\in\operatorname{Lip}(C_{L},Q_{d}(\lambda))$,
pick $q=\lceil 176\cdot 5^{T}\lambda C_{L}d/\epsilon\rceil$ and final accuracy
$\epsilon_{f}=\frac{\epsilon}{2\cdot 5^{T}}$ in Corollary 3.15, which implies
time

$\displaystyle\widetilde{O}\left(\frac{Tm^{6}}{\epsilon_{f}^{2}}\frac{R^{2}L^{4}}{\sigma_{\min}^{4}}\mathcal{T}_{\rm
total}\right)=\widetilde{O}\left(T(\lambda
C_{L}d)^{6d}\frac{R^{2}L^{4}}{\sigma_{\min}^{4}}\left(\frac{5^{T}}{\epsilon}\right)^{2+6d}\mathcal{T}_{\rm
total}\right)$

up to the logarithm factor $\log(m^{2}/\delta)$.

By assuming sufficiently smooth conditions on the continuation values, we
obtain the final classical and quantum complexities
$\widetilde{O}((5^{T}/\epsilon)^{2+6d/n})$ and
$\widetilde{O}((5^{T}/\epsilon)^{1+4d/n})$, respectively. If the continuation
values are $n$-differentiable functions for
$n=\Theta(\log(5^{T}/\epsilon)/\log\log(5^{T}/\epsilon))$ — note that $n\leq
q=O((5^{T}/\epsilon)^{1/n})$ —, then we obtain the simpler runtimes
$\widetilde{O}((5^{T}/\epsilon)^{2})$ and $\widetilde{O}(5^{T}/\epsilon)$,
thus constituting a quadratic improvement from classical to quantum
complexity. The exponential time dependence comes from the error precision
required to obtain a final additive accuracy of $\epsilon$, which is exactly
the quantity that is sped up by using quantum algorithm for Monte Carlo, hence
the quadratic improvement on time. Finally, we note that, given the condition
$n\leq q$ and the value $q=O((5^{T}/\epsilon)^{1/n})$, it is not possible to
take $n\to\infty$. The value
$n=\Theta(\log(5^{T}/\epsilon)/\log\log(5^{T}/\epsilon))$ implies
$n\log{n}=\Theta\left(\frac{\log(5^{T}/\epsilon)}{\log\log(5^{T}/\epsilon)}\log\left(\frac{\log(5^{T}/\epsilon)}{\log\log(5^{T}/\epsilon)}\right)\right)=O(\log(5^{T}/\epsilon))$,
and thus does not violate $n\leq q$.

Theorem 4.3 applies to any set of polynomials, e.g. Chebyshev polynomials. We
can thus compare the above result with the one from Miyamoto [69], which
approximates the continuation values by a Chebyshev interpolation through
Chebyshev nodes and obtains a complexity of
$O(\epsilon^{-1}\log^{d}(1/\epsilon)\text{poly}\log\log(1/\epsilon))$. Our
approach, in contrast, is to project $\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]$ onto a
set of polynomials and is, for this reason, much more general. Moreover, our
final result is a _time_ complexity, while the result from [69] is a _query_
complexity on the number of unitaries called by all quantum routines for Monte
Carlo (thus being a weaker result). Finally, Miyamoto [69] assumes that the
continuation values are analytical functions, i.e., are in $C^{\infty}$, while
we only need to assume $\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\in
C^{n}(Q_{d}(\lambda))$ for
$n=\Theta(\log(5^{T}/\epsilon)/\log\log(5^{T}/\epsilon))$ in order to recover
$\widetilde{O}(\epsilon^{-1})$ up to polylog factors. One downside of our
approach, though, is the presence of the (non-zero) minimum singular value
$\sigma_{\min}$ which depends on the probability measures $\rho_{t}$.
Chebyshev interpolation bypasses such direct probability dependence by
shifting it entirely to the approximation of the continuation values on
Chebyshev nodes via (quantum) Monte Carlo.

We can take advantage of our freedom in choosing the approximation
architecture in order to simplify the LSM algorithm, specially the matrices
$A_{t}$, depending on the underlying probability model. In this section we
analyze the common cases of Brownian motion and geometric Brownian motion.

###  4.1 Brownian motion

Consider a Markov chain $(X_{t})_{t=0}^{T}$ where each component of
$X_{t}:\Omega\to\mathbb{R}^{d}$ follows an independent Brownian motion. A
Brownian motion can be described by a stochastic process $W=(W_{t})_{t=0}^{T}$
(called Wiener process) such that:

  1. 1.

$W_{0}=0$;

  2. 2.

$W_{t}$ is almost surely continuous in $t$;

  3. 3.

$W_{t+s}-W_{s}$ for $t\geq 0$ is independent of $W_{s}$, for $s\leq t$;

  4. 4.

$W_{t+u}-W_{t}\sim\mathcal{N}(0,u)$ for $u\geq 0$, where
$\mathcal{N}(\mu,\sigma^{2})$ denotes the normal distribution.

It follows from the above definition that the image probability measure
$\rho_{t}$ of $X_{t}$ reads $\rho_{t}(x)=\prod_{i=1}^{d}\frac{1}{\sqrt{2\pi
t}}e^{-x_{i}^{2}/2t}$ for $t\in[T]$. Given $t\in[T]$, consider the
approximation architecture generated by
$\\{e_{t,\vec{k}}:\mathbb{R}^{d}\to\mathbb{R}\\}_{\vec{k}\in\mathcal{B}}$,
where $\mathcal{B}:=\\{\vec{k}\in\\{0,\dots,q\\}^{d}:\sum_{i=1}^{d}k_{i}\leq
q\\}$ for some integer $q\geq 0$ and

$\displaystyle
e_{t,\vec{k}}(x)=\prod_{l=1}^{d}\frac{1}{\sqrt{k_{l}!2^{k_{l}}}}H_{k_{l}}(x_{l}/\sqrt{2t}),$
(14)

with $x\in\mathbb{R}^{d}$ and the $H_{l}(y)$ are Hermite polynomials defined
by

$\displaystyle
H_{l}(y):=(-1)^{l}e^{y^{2}}\frac{\text{d}^{l}}{\text{d}y^{l}}e^{-y^{2}}$

where $l\in\mathbb{N}$ and specifically $H_{0}(x)=1$. The choice for Hermite
polynomials comes from the fact that they form an orthogonal basis under the
measure $\rho_{t}(x)$, i.e.,

$\displaystyle\int_{-\infty}^{\infty}H_{k}(x/\sqrt{2t})H_{l}(x/\sqrt{2t})\frac{e^{-x^{2}/2t}}{\sqrt{2\pi
t}}\text{d}x=\int_{-\infty}^{\infty}H_{k}(x)H_{l}(x)\frac{e^{-x^{2}}}{\sqrt{\pi}}\text{d}x=k!2^{k}\delta_{kl}.$

This means that, for each $t\in[T-1]$,
$\\{e_{t,\vec{k}}\\}_{\vec{k}\in\mathcal{B}}$ are orthonormal under
$\rho_{t}(x)$ and thus the matrices $A_{t}$ become simply the identity matrix.
However, in order to invoke approximation results due to smoothness
assumptions, e.g. the Jackson-like inequality in Lemma 4.1, we must restrict
the domain to a close set like $Q_{d}(\lambda)$. We shall then consider the
approximation architecture generated by $\\{e^{\rm
trunc}_{t,\vec{k}}:\mathbb{R}^{d}\to\mathbb{R}\\}_{\vec{k}\in\mathcal{B}}$
with $e^{\rm trunc}_{t,\vec{k}}$ defined as in Eq. (14) if $x\in
Q_{d}(\lambda)$ and $e^{\rm trunc}_{t,\vec{k}}(x)=0$ otherwise. The implied
matrix $A^{\rm trunc}_{t}$ is no longer equal to the identity matrix, but is
close to it. Such closeness will be approximated with the aid of the following
result whose proof is postponed to Appendix A.

######  Lemma 4.4.

Given $\lambda>0$ and $k,l\in\mathbb{N}$ with $l\leq k$, then

$\displaystyle\left|\int_{\lambda}^{\infty}H_{k}(x)H_{l}(x)e^{-x^{2}}\text{d}x\right|\leq
2^{2+(k+l)/2}\sqrt{(k+1)!(l+1)!}e^{(\sqrt{2(k+1)}+\sqrt{2(l+1)})\lambda}e^{-\lambda^{2}}.$

This means that, if we approximate the matrices $A^{\rm trunc}_{t}$ by the
identity matrix, the error incurred will be small for sufficiently large
$\lambda$ and, therefore, the complexity of Algorithms 1 and 2 can be reduced
since the matrices $\widetilde{A}_{t}$ no longer need to be computed and
inverted. The vector $b_{t}$, $t\in[T-1]$, will still be approximated using
quantum algorithm for Monte Carlo, which, in turn, must use a discretization
of the Brownian motion. In order to formalize these points in a theorem, we
first define the truncation operator $\mathscr{T}_{\beta}$ such that, for any
function $f:\mathbb{R}^{d}\to\mathbb{R}$,

$\displaystyle\mathscr{T}_{\beta}f(x)=\begin{cases}f(x)\quad&\text{if}~{}|f(x)|\leq\beta,\\\
\beta\operatorname{sgn}(f(x))\quad&\text{otherwise}.\end{cases}$

Define also the quantity

$\displaystyle
r_{p}:=T\sqrt{\frac{2M_{p}}{p-2}}+\operatorname*{\max}_{t\in[T-1]}\sum_{i=1}^{d}\mathbb{E}[|X_{t,i}|]=T\sqrt{\frac{2M_{p}}{p-2}}+d\sqrt{\frac{2(T-1)}{\pi}}$

for $2<p<\infty$, where
$M_{p}:=\max\big{\\{}\|z_{1}\|_{L^{p}(\rho_{1})}^{p},\dots,\|z_{T}\|_{L^{p}(\rho_{T})}^{p}\big{\\}}$
and we used $\mathbb{E}[|X_{t,i}|]=\sqrt{\frac{2t}{\pi}}$. For simplicity, the
next theorem just considers the case when the continuation values are in
$C^{n}(Q_{d}(\lambda))$. Moreover, we stress that
$L^{2}(\rho_{t})=L^{2}(\mathbb{R}^{d},\rho_{t})$ considers the whole space
$\mathbb{R}^{d}$.

######  Theorem 4.5.

Consider a Markov chain $(X_{t})_{t=0}^{T}$ where each component of
$X_{t}:\Omega\to\mathbb{R}^{d}$ follows an independent Brownian motion. Let
$Q_{d}(\lambda)$ be the closed cube with side length $2\lambda$ in
$\mathbb{R}^{d}$ centered at the origin. Suppose that
$M_{p}:=\max\big{\\{}\|z_{1}\|_{L^{p}(\rho_{1})}^{p},\dots,\|z_{T}\|_{L^{p}(\rho_{T})}^{p}\big{\\}}<\infty$
for some $2<p<\infty$. Within the setting of Algorithm 2 with parameters
$\delta$ and $\epsilon_{0}$, let $T\mathcal{T}_{\rm samp}$ be the sampling
cost of the Markov chain and set all matrices $\widetilde{A}_{t}$,
$t\in[T-1]$, to be the identity matrix and use the payoffs
$z_{0},\mathscr{T}_{\beta}z_{1},\dots,\mathscr{T}_{\beta}z_{T}$, where
$\beta:=\lambda^{2/p}$, with query cost $\mathcal{T}_{z}$ and
$R:=\max_{t\in[T]}\|\mathscr{T}_{\beta}z_{t}\|_{u}$. Also, for $t\in[T-1]$,
consider the approximation architectures $\mathscr{H}_{t}=\\{e^{\rm
trunc}_{t,\vec{k}}:\mathbb{R}^{d}\to\mathbb{R}\\}_{\vec{k}\in\mathcal{B}}$
where $\mathcal{B}:=\\{\vec{k}\in\\{0,\dots,q\\}^{d}:\sum_{i=1}^{d}k_{i}\leq
q\\}$ for some $q\geq 0$ and $e^{\rm trunc}_{t,\vec{k}}(x)$ is defined as in
Eq. (14) if $x\in Q_{d}(\lambda)$ and $e^{\rm trunc}_{t,\vec{k}}(x)=0$
otherwise. Let $\mathcal{T}_{e}$ be the query cost of $\mathscr{H}_{t}$,
$t\in[T]$. Suppose there is $n\in\mathbb{N}$, $n<q$, such that the restriction
$\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]|_{Q_{d}(\lambda)}$ to $Q_{d}(\lambda)$ of
each $\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]$, $t\in\\{0,\dots,T-1\\}$, is in
$C^{n}(Q_{d}(\lambda))$. Let $C\geq C_{d,\lambda}\max_{0\leq
t<T}\|\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]|_{Q_{d}(\lambda)}\|_{C^{n}(Q_{d}(\lambda))}$
be known, where $C_{d,\lambda}$ is a constant depending only on
$Q_{d}(\lambda)$. Then, for $\delta\in(0,1)$,
$\epsilon\in\big{(}0,5^{T}(e^{2d}C^{d/n}R)^{n/(n+d)}\big{]}$,
$q=\left\lceil(3\cdot 5^{T}C/\epsilon)^{1/n}\right\rceil$ and

$\displaystyle\lambda=O\left(\max\left\\{\sqrt{T}\left(\frac{5^{T}}{\epsilon}\right)^{1/2n},T\log\left((e^{2n}C)^{2d/n}Rd\left(\frac{5^{T}}{\epsilon}\right)^{1+2d/n}\right),\left(\frac{r_{p}}{\epsilon}\right)^{p/(p-2)}\right\\}\right),$
(15)

if $\epsilon_{0}=\frac{1}{3(e^{2n}C)^{d/n}R}\big{(}\frac{\epsilon}{3\cdot
5^{T}}\big{)}^{1+d/n}$, Algorithm 2 outputs
$\widetilde{\mathcal{U}}_{0}:=\max\big{\\{}Z_{0},\widetilde{(\mathscr{T}_{\beta}Z)}_{\widetilde{\tau}_{1}}\big{\\}}$
such that

$\displaystyle\operatorname{Pr}\left[|\widetilde{\mathcal{U}}_{0}-\mathbb{E}[Z_{\tau_{0}}]|\leq\epsilon\right]\geq
1-\delta$

in time

$\displaystyle\widetilde{O}\left(T^{2}\mathcal{T}_{\rm
total}(e^{2n}C)^{7d/2n}R^{2}\left(\frac{5^{T}}{\epsilon}\right)^{1+7d/2n}\right),$

where $\mathcal{T}_{\rm total}:=\mathcal{T}_{\rm
samp}+\mathcal{T}_{z}+\mathcal{T}_{e}$ and $\widetilde{O}(\cdot)$ hides
$\operatorname{polylog}$ terms in $5^{T},R,C,e^{d},\epsilon,\delta$.

###### Proof.

First we prove that Algorithm 2, when using $m:=|\mathcal{B}|=\binom{q+d}{d}$
different Hermite polynomials from Eq. (14) as approximation architecture and
parameters $\epsilon_{b}=\epsilon_{0}/\sqrt{m}$ and $\delta_{b}=\delta/(2Tm)$
(there are no $\epsilon_{A}$ and $\delta_{A}$ parameters), outputs
$\widetilde{(\mathscr{T}_{\beta}Z)}_{\widetilde{\tau}_{1}}$ such that, for any
$\epsilon_{0}\in(0,1/2]$,

$\displaystyle\operatorname{Pr}\left[\big{|}\widetilde{(\mathscr{T}_{\beta}Z)}_{\widetilde{\tau}_{1}}-\mathbb{E}[(\mathscr{T}_{\beta}Z)_{\tau_{1}}]\big{|}\geq
5^{T}\left(3m\epsilon_{0}R+\operatorname*{\max}_{0<t<T}\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbb{E}[(\mathscr{T}_{\beta}Z)_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}\right)\right]\leq\delta,$

where $\vec{e}_{t}(\cdot):=(e^{\rm
trunc}_{t,\vec{k}}(\cdot))_{\vec{k}\in\mathcal{B}}$, in time

$\displaystyle O\left(\frac{T^{2}m^{5/2}}{\epsilon_{0}}\mathcal{T}_{\rm
total}R\log(T)\log(Tm/\delta)\log^{3/2}(\sqrt{m}R/\epsilon_{0})\log\log(\sqrt{m}R/\epsilon_{0})\right).$
(16)

Indeed, the proof is very similar to the ones from Theorems 3.11 and 3.13.
According to Lemma 4.4, approximating the matrices $A^{\rm trunc}_{t}$ by the
identity matrix leads to an error in any diagonal entry $\mathbb{E}[e^{\rm
trunc}_{t,\vec{k}}(X_{t})e^{\rm trunc}_{t,\vec{k}}(X_{t})]$ of at most

$\displaystyle
1-\prod_{i=1}^{d}\left(1-\int_{(-\infty,-\lambda]\cup[\lambda,\infty)}\frac{H_{k_{i}}(x/\sqrt{2t})^{2}}{k_{i}!2^{k_{i}}}\frac{e^{-x^{2}/2t}}{\sqrt{2\pi
t}}\text{d}x\right)$ $\displaystyle\leq
2\sum_{i=1}^{d}\left|\int_{\lambda}^{\infty}\frac{H_{k_{i}}(x/\sqrt{2t})^{2}}{k_{i}!2^{k_{i}}}\frac{e^{-x^{2}/2t}}{\sqrt{2\pi
t}}\text{d}x\right|$ (17)
$\displaystyle\leq\frac{8e^{-\lambda^{2}/2t}}{\sqrt{\pi}}\sum_{i=1}^{d}(k_{i}+1)e^{2\lambda\sqrt{(k_{i}+1)/t}}$
$\displaystyle\leq\frac{8e^{-\lambda^{2}/2t}}{\sqrt{\pi}}\sum_{i=1}^{d}e^{4\lambda\sqrt{(k_{i}+1)/t}}$
(18)
$\displaystyle\leq\frac{8d}{\sqrt{\pi}}e^{4\lambda\sqrt{(q+1)/t}}e^{-\lambda^{2}/2t}.$

where Eq. (17) follows from convexity and Eq. (18) uses that $(k_{i}+1)\leq
e^{2\lambda\sqrt{(k_{i}+1)/t}}$ if $\lambda\geq 16\sqrt{(q+1)T}$. If
$\lambda\geq\max\big{\\{}16\sqrt{(q+1)T},4T\log(5md/\epsilon_{0})\big{\\}}$,
the above error is at most
$\frac{8d}{\sqrt{\pi}}e^{-\lambda^{2}/4t}\leq\frac{\epsilon_{0}}{m}$. A
similar bound can be obtained for off-diagonal terms $\mathbb{E}[e^{\rm
trunc}_{t,\vec{k}}(X_{t})e^{\rm trunc}_{t,\vec{l}}(X_{t})]$:

$\displaystyle\prod_{i=1}^{d}\left|\int_{(-\infty,-\lambda]\cup[\lambda,\infty)}\frac{H_{k_{i}}(x/\sqrt{2t})H_{l_{i}}(x/\sqrt{2t})}{\sqrt{k_{i}!l_{i}!2^{k_{i}+l_{i}}}}\frac{e^{-x^{2}/2t}}{\sqrt{2\pi
t}}\text{d}x\right|\leq\left(\frac{8}{\sqrt{\pi}}e^{4\lambda\sqrt{(q+1)/t}}e^{-\lambda^{2}/2t}\right)^{d}\leq\frac{\epsilon_{0}}{m}$

using the bounds on $\lambda$. The above bounds imply the $L^{2}$-error
between $A^{\rm trunc}_{t}$ and the identity is $\|A^{\rm
trunc}_{t}-\mathbb{I}\|_{2}\leq\epsilon_{0}$. The $L^{2}$-error between
$b_{t}$ and $\widetilde{b}_{t}$ is still
$\|b_{t}-\widetilde{b}_{t}\|_{2}\leq\epsilon_{0}$ with probability at least
$1-\delta/2T$. According to Theorem 3.10, we obtain

$\displaystyle\|\widetilde{\alpha}_{t}-\alpha_{t}\|_{2}=\|\widetilde{b}_{t}-(A^{\rm
trunc}_{t})^{-1}b_{t}\|_{2}\leq\|\widetilde{b}_{t}-b_{t}\|_{2}+\|b_{t}-(A^{\rm
trunc}_{t})^{-1}b_{t}\|_{2}\leq\epsilon_{0}+2\epsilon_{0}\|b_{t}\|_{2}\leq
3\epsilon_{0}\sqrt{m}R,$

where we used that $\sigma_{\min}(\mathbb{I})=1$ and that $\|b_{t}\|_{2}\leq
R\|\mathbb{E}[\vec{e}_{t}(X_{t})]\|_{2}\leq\sqrt{m}R$ since, by the definition
of $e^{\rm trunc}_{t,\vec{k}}$, we have
$\max_{\vec{k}\in\mathcal{B},t\in[T]}\|e^{\rm
trunc}_{t,\vec{k}}\|_{L^{2}(\rho_{t})}\leq 1$. This leads, with probability at
least $1-\delta/2T$, to

$\displaystyle\|\widetilde{\alpha}_{t}\cdot\vec{e}_{t}(X_{t})-\alpha_{t}\cdot\vec{e}_{t}(X_{t})\|_{L^{2}(\rho_{t})}\leq\|\widetilde{\alpha}_{t}-\alpha_{t}\|_{2}\big{\|}\sqrt{\vec{e}_{t}(X_{t})\cdot\vec{e}_{t}(X_{t})}\big{\|}_{L^{2}(\rho_{t})}\leq
3m\epsilon_{0}R,$

using
$\big{\|}\sqrt{\vec{e}_{t}(X_{t})\cdot\vec{e}_{t}(X_{t})}\big{\|}_{L^{2}(\rho_{t})}=\sqrt{\sum_{x\in
Q_{d}(\lambda)}\rho_{t}(x)\sum_{\vec{k}\in\mathcal{B}}\big{(}e^{\rm
trunc}_{t,\vec{k}}(x)\big{)}^{2}}\leq\sqrt{m}$. The rest of the error analysis
is the same and we omit it here, the only difference being that the payoffs
are now $z_{0},\mathscr{T}_{\beta}z_{1},\dots,\mathscr{T}_{\beta}z_{T}$.
Regarding the complexity, the most expensive step is the computation of all
$b_{t}$, which requires time

$\displaystyle O\left(\frac{Tm}{\epsilon_{b}}R(T\mathcal{T}_{\rm
samp}+T(\log(T)\mathcal{T}_{z}+m\mathcal{T}_{e}))\log\left(\frac{1}{\delta_{b}}\right)\log^{3/2}\left(\frac{R}{\epsilon_{b}}\right)\log\log\left(\frac{R}{\epsilon_{b}}\right)\right),$

by calling
$\mathtt{QMonteCarlo}(Z_{\widetilde{\tau}_{t}}e_{t-1,\vec{k}}(X_{t-1}),\epsilon_{b},\delta_{b},R)$
from Theorem 3.7 for $t\in\\{2,\dots,T\\}$ and $\vec{k}\in\mathcal{B}$. The
term $T(\log(T)\mathcal{T}_{z}+m\mathcal{T}_{e})$ comes from the unitaries
$C_{t}^{(\vec{k})}$ in $\mathtt{QMonteCarlo}$. Hence, by keeping the largest
terms of each complexity and using that $\epsilon_{b}=\epsilon_{0}/\sqrt{m}$
and $\delta_{b}=\delta/(4Tm)$, the final complexity is the one from Eq. (16).

With this partial result in hands, we again use the Jackson-type inequality
from Lemma 4.1 as in Theorem 4.3. In order to do so, we first need the
following technical result that limits the continuation values to being a
function in the domain $Q_{d}(\lambda)$. More specifically, by [107,
Proposition 3.14],

$\displaystyle\|\mathbf{1}_{Q_{d}(\lambda)}\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]-\mathbb{E}[(\mathscr{T}_{\beta}Z)_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}\leq
r_{p}\lambda^{2/p-1},$ (19)

where $\mathbf{1}_{Q_{d}(\lambda)}$ is the indicator random variable
corresponding to the set $Q_{d}(\lambda)$. Then, since $\\{e^{\rm
trunc}_{t,\vec{k}}\\}_{\vec{k}\in\mathcal{B}}$ spans
$\mathcal{R}_{q}(Q_{d}(\lambda))$, the Jackson-type inequality from Lemma 4.1
reads

$\displaystyle\operatorname*{\min}_{a\in\mathbb{R}^{m}}\|a\cdot\vec{e}_{t}(X_{t})-\mathbf{1}_{Q_{d}(\lambda)}\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]\|_{L^{2}(\rho_{t})}=\operatorname*{\inf}_{p\in\mathcal{R}_{q}}\|(p-\mathbb{E}[Z_{\tau_{t+1}}|X_{t}])\mathbf{1}_{Q_{d}(\lambda)}\|_{L^{2}(\rho_{t})}\leq
Cq^{-n},$

for all $t\in\\{0,\dots,T-1\\}$. Combining both results, we finally get that

$\displaystyle\big{|}\widetilde{(\mathscr{T}_{\beta}Z)}_{\widetilde{\tau}_{1}}-\mathbb{E}[(\mathscr{T}_{\beta}Z)_{\tau_{1}}]\big{|}\leq
5^{T}\left(3m\epsilon_{0}R+Cq^{-n}+r_{p}\lambda^{2/p-1}\right),$

with probability at least $1-\delta$. Before we plug values for $m$, $q$,
$\epsilon_{0}$ and $\epsilon$, we get rid of the truncation operator
$\mathscr{T}_{\beta}$ in the inequality above by using that [32, Proposition
5.2]

$\displaystyle\left\|\mathbb{E}[Z_{\tau_{t+1}}|X_{t}]-\mathbb{E}[(\mathscr{T}_{\beta}Z)_{\tau_{t+1}}|X_{t}]\right\|_{L^{2}(\rho_{t})}\leq
r_{p}\lambda^{1-p/2}$

(which is quite similar to Eq. (19)). By using that
$|{\max}\\{a_{0},a_{1}\\}-\max\\{a_{0},a_{2}\\}|\leq|a_{1}-a_{2}|$ with
$a_{0},a_{1},a_{2}\in\mathbb{R}$, then

$\displaystyle|\widetilde{\mathcal{U}}_{0}-\mathbb{E}[Z_{\tau_{0}}]|\leq\big{|}\widetilde{(\mathscr{T}_{\beta}Z)}_{\widetilde{\tau}_{1}}-\mathbb{E}[Z_{\tau_{1}}]\big{|}\leq
5^{T}\left(3m\epsilon_{0}R+Cq^{-n}+2r_{p}\lambda^{2/p-1}\right),$

since $\lambda^{1-p/2}\leq\lambda^{2/p-1}$.

We use that $m\leq e^{d}q^{d}$, $q=\lceil(3\cdot
5^{T}C/\epsilon)^{1/n}\rceil$,
$\lambda\geq(\frac{6r_{p}}{\epsilon})^{p/(p-2)}$ and take
$\epsilon_{0}=\frac{1}{3C^{d/n}e^{2d}R}\big{(}\frac{\epsilon}{3\cdot
5^{T}}\big{)}^{1+d/n}$ (note that $\epsilon_{0}\in(0,1/2]$) in order to obtain
the final error

$\displaystyle
5^{T}\left(3m\epsilon_{0}R+\frac{C}{q^{n}}+2r_{p}\lambda^{2/p-1}\right)\leq
5^{T}\left(3\epsilon_{0}e^{d}q^{d}R+\frac{C}{q^{n}}+2r_{p}\lambda^{2/p-1}\right)\leq(1/3+1/3+1/3)\epsilon=\epsilon.$

The above in turn implies the complexity, up to $\operatorname{polylog}$
factors,

$\widetilde{O}\left(\frac{T^{2}m^{5/2}}{\epsilon_{0}}\mathcal{T}_{\rm