#  Towards Quantum Advantage in Financial Market Risk using Quantum Gradient
Algorithms

Nikitas Stamatopoulos  Goldman, Sachs & Co., New York, NY  Guglielmo Mazzola
IBM Quantum, IBM Research – Zurich  Stefan Woerner  IBM Quantum, IBM Research
– Zurich  William J. Zeng  Goldman, Sachs & Co., New York, NY

###### Abstract

We introduce a quantum algorithm to compute the market risk of financial
derivatives. Previous work has shown that quantum amplitude estimation can
accelerate derivative pricing quadratically in the target error and we extend
this to a quadratic error scaling advantage in market risk computation. We
show that employing quantum gradient estimation algorithms can deliver a
further quadratic advantage in the number of the associated market
sensitivities, usually called _greeks_. By numerically simulating the quantum
gradient estimation algorithms on financial derivatives of practical interest,
we demonstrate that not only can we successfully estimate the greeks in the
examples studied, but that the resource requirements can be significantly
lower in practice than what is expected by theoretical complexity bounds. This
additional advantage in the computation of financial market risk lowers the
estimated logical clock rate required for financial quantum advantage from
Chakrabarti et al. [Quantum 5, 463 (2021)] by a factor of 50, from 50MHz to
1MHz, even for a modest number of greeks by industry standards (four).
Moreover, we show that if we have access to enough resources, the quantum
algorithm can be parallelized across 30 QPUs for the same overall runtime as
the serial execution if the logical clock rate of each device is $\sim 30$kHz,
same order of magnitude as the best current estimates of feasible target clock
rates of around 10kHz. Throughout this work, we summarize and compare several
different combinations of quantum and classical approaches that could be used
for computing the market risk of financial derivatives.

##  I Introduction

Recently, quantum algorithms have been proposed to accelerate the pricing and
risk analysis of financial derivatives Rebentrost  _et al._ (2018); Woerner
and Egger (2019); Egger  _et al._ (2020); Stamatopoulos  _et al._ (2020);
Chakrabarti  _et al._ (2021). These algorithms use quantum amplitude
estimation to achieve quadratic advantage compared to the classical Monte
Carlo methods that are used in practice for most computationally expensive
pricing. Let $\epsilon_{p}$ be the error in pricing. The quantum advantage
stems from the runtime of a classical Monte Carlo simulation scaling as
$\mathcal{O}\left(1/\epsilon_{p}^{2}\right)$, while the quantum algorithms
have scaling $\mathcal{O}\left(1/\epsilon_{p}\right)$ Montanaro (2015).

A related and important financial application is the computation of the
sensitivity of derivative prices to model and market parameters. This amounts
to computing gradients of the derivative price with respect to input
parameters. A primary business use of calculating these gradients is to enable
hedging of the market risk that arises from exposure to derivative contracts.
Hedging this risk is of critical importance to financial firms Hull (2006). In
some cases, gradients can be computed analytically, e.g. when the derivative
price has an analytical form. In this work, we consider computing gradients
where there is not a closed form for the price and where the classical
comparison is to Monte Carlo simulation. Gradients of financial derivatives
are typically called  _greeks_ , as these quantities are commonly labeled
using Greek alphabet letters. For $k$ greeks, i.e. for a $k$-dimensional
gradient, _classical finite difference_ methods compute the price at multiple
points in each parameter dimension for a scaling of
$\mathcal{O}\left(k/\epsilon_{p}^{2}\right)$.

One approach to quantum acceleration of greek computation is to construct a
finite difference approximation of the financial derivative’s price on a
quantum computer and perform amplitude estimation on that quantity instead of
its price. In Sec. II.4 we show that this approach, called the  _semi-
classical_ method in Gilyén  _et al._ (2019), scales as
$\mathcal{O}\left(k/\epsilon\right)$, where $\epsilon$ is the target error for
the gradients. However, one can further improve the scaling in the number of
greeks $k$, using quantum algorithms for computing the gradients Jordan
(2005); Gilyén  _et al._ (2019). A quantum algorithm to compute the gradient’s
components in superposition was originally introduced by Jordan in Ref. Jordan
(2005) and more recently revisited by Gilyen, Arunachalam, and Wiebe (GAW)
Gilyén  _et al._ (2019).

Gilyen et al. Gilyén  _et al._ (2019) perform a rigorous analysis of Jordan’s
original quantum gradient algorithm Jordan (2005) to show that in general it
scales as $\mathcal{O}\left(\sqrt{k}/\epsilon^{2}\right)$. They then
generalize that algorithm to arbitrary higher order $m$ and show that an order
$m=\log(\sqrt{k}/\epsilon)$ quantum gradient algorithm has runtime of
$\mathcal{O}\left(\sqrt{k}/\epsilon\right)$ for a specific class of smooth
functions. For the sake of simplicity, in the rest of the manuscript we call
this result as _GAW quantum gradient_ algorithm, from the name of the authors
of Ref. Gilyén  _et al._ (2019). Quantum gradient algorithms have also been
previously studied as a subroutine for accelerated convex optimization
Chakrabarti  _et al._ (2020). Since, in this optimization context, subtleties
connected with the error’s scaling and parameters setting were found to be
important, it is timely to put forward a comprehensive study on the advantage
that could be found in financial applications.

In this work, we apply the algorithmic framework from Gilyén  _et al._ (2019)
to compute financial greeks and study several varieties of quantum gradient
algorithms. We use the quantum pricing method of Chakrabarti et al.
Chakrabarti  _et al._ (2021) as a subroutine to the GAW algorithm to
numerically estimate the resource requirements of computing the greeks of two
types of option contracts (a) a European call option, which we use as a
benchmark to establish the validity of the algorithm and estimate the
corresponding resource requirements, and (b) a path-dependent basket option
whose pricing profile is representative of typical financial derivative
contracts of practical interest.

Then, we introduce a second-order accurate quantum gradient algorithm, the
$m=1$ version of Gilyén  _et al._ (2019), for which we can give an explicit
and compact quantum implementation in the financial derivative case, and which
does not rely on block encoding or Hamiltonian simulation. We use this
_Simulation-Free Quantum Gradient_ (SFQG) method to compute the greeks of a
path-dependent derivative and find that it is significantly cheaper to
construct than the Hamiltonian-based method. Additionally, we show that we can
improve the overall performance of quantum gradient estimation algorithms by
employing a maximum likelihood (MLE) method to extract the most likely
estimate of the gradients with concrete confidence intervals. With these
tools, we calculate that quantum advantage for calculating risk may be
achievable with quantum computers whose clock rates are 50 times slower than
that required for pricing itself. We discuss these implications in more detail
in Sec. VII.

Finally, we perform a comparison between quantum, classical and semi-classical
gradient estimation algorithms in the context of financial derivatives,
summarized in Table 4.

We highlight the new contributions in this work:

  * •

We numerically study the quantum gradient estimation algorithms from Jordan
(2005); Gilyén  _et al._ (2019) for functions of practical interest to
financial market risk and compare the observed oracular cost to theoretical
expectations. (Section IV)

  * •

We devise a method to construct a second-order accurate oracle for quantum
gradient estimation for functions computed using quantum amplitude estimation
that is cheaper in required resources compared to existing methods. (Section
V)

  * •

We introduce a way to improve gradient estimation algorithms using classical
maximum likelihood estimation (MLE). (Section VI)

  * •

We propose a technique to employ automatic differentiation (AD) methods on
quantum computers which can enhance the quantum gradient estimation
performance in certain cases. (Appendix A)

  * •

We update the resource estimates for quantum advantage in financial derivative
pricing from prior research. (Section VII)

###  I.1 Quantum Amplitude Estimation

Our gradient calculation algorithms extend the quantum accelerated method for
derivative pricing that is based on quantum amplitude estimation (QAE)
Brassard  _et al._ (2002). We review QAE here. Let the $k$-dimensional vector
$\mathbf{x}$ denote the set of market data parameters of the derivative
represented as a basis state $\ket{\mathbf{x}}$ and $f(\mathbf{x})$ its price,
rescaled to satisfy $f(\mathbf{x})\in[0,1]$. We assume the existence of a
unitary operator $\mathcal{A}(\mathbf{x})$ which produces the state

$\mathcal{A}:\ket{\vec{0}}\ket{\mathbf{x}}\rightarrow\left(\sqrt{1-f(\mathbf{x})}\ket{\psi_{0}(\mathbf{x})}\ket{0}+\sqrt{f(\mathbf{x})}\ket{\psi_{1}(\mathbf{x})}\ket{1}\right)\ket{\mathbf{x}},$
(1)

where $\ket{\psi_{0}(\mathbf{x})}$ and $\ket{\psi_{1}(\mathbf{x})}$ are
arbitrary, normalized quantum states. The state $\ket{\mathbf{x}}$
representing the market data parameters acts only as input to $\mathcal{A}$
and we can ignore it for the remainder of this section. When an operator
$\mathcal{U}$ depends on some input $\ket{\mathbf{x}}$, we will henceforth
write it as $\mathcal{U}(\mathbf{x})$ to highlight the dependence.

Explicit constructions of $\mathcal{A}(\mathbf{x})$ with resource estimates
for specific path dependent derivative instances are given in Chakrabarti et
al. Chakrabarti  _et al._ (2021). QAE estimates $f(\mathbf{x})$ with repeated
applications of the operator
$\mathcal{Q}(\mathbf{x})=\mathcal{A}(\mathbf{x})S_{0}\mathcal{A}^{\dagger}(\mathbf{x})S_{\psi_{0}}$,
where $S_{0}=\mathbb{I}-2\ket{\vec{0}}\bra{\vec{0}}$ and
$S_{\psi_{0}}=\mathbb{I}-2\ket{\psi_{0}(\mathbf{x})}\ket{0}\bra{0}\bra{\psi_{0}(\mathbf{x})}$.
The QAE algorithm relies on the fact that, by the construction of
$\mathcal{Q}(\mathbf{x})$, the state produced after the application of
$\mathcal{A}(\mathbf{x})$ in Eq. (1) can be written in the eigenbasis of
$\mathcal{Q}(\mathbf{x})$ Brassard  _et al._ (2002)

$\frac{-i}{\sqrt{2}}\left(e^{i\theta(\mathbf{x})}\ket{\psi_{+}(\mathbf{x})}-e^{-i\theta(\mathbf{x})}\ket{\psi_{-}(\mathbf{x})}\right)\equiv\ket{\Psi(\mathbf{x})},$
(2)

where $\ket{\psi_{+}(\mathbf{x})}$ and $\ket{\psi_{-}(\mathbf{x})}$ are the
eigenvectors of $\mathcal{Q}(\mathbf{x})$ with eigenvalues $e^{\pm
2i\theta(\mathbf{x})}$, and $f(\mathbf{x})=\sin^{2}(\theta(\mathbf{x}))$. Let
integer $m>0$ denote the bits of precision with which we want to estimate the
amplitude. Applying $H^{\otimes m}$ on an $m$-qubit register initialized at
$\ket{0}_{m}$ and using the register to control different powers of
$\mathcal{Q}(\mathbf{x})$, a phase kickback of the eigenvalues of
$\mathcal{Q}(\mathbf{x})$ is induced onto the control register

$\frac{1}{\sqrt{M}}\sum_{j=0}^{M-1}\ket{j}\mathcal{Q}^{2^{j}}(\mathbf{x})\ket{\Psi(\mathbf{x})}=\frac{1}{\sqrt{M}}\sum_{j=0}^{M-1}e^{2ij\theta(\mathbf{x})}\ket{j}\ket{\Psi_{+}(\mathbf{x})}-\frac{1}{\sqrt{M}}\sum_{j=0}^{M-1}e^{-2ij\theta(\mathbf{x})}\ket{j}\ket{\Psi_{-}(\mathbf{x})},$
(3)

where $\ket{\Psi_{\pm}(\mathbf{x})}=-ie^{\pm
i\theta(\mathbf{x})}\ket{\psi_{\pm}(\mathbf{x})}/\sqrt{2}$ and $M=2^{m}$.
Finally, applying an inverse Quantum Fourier Transform on the first register
and measuring in the computational basis gives an $m$-bit approximation of
either $\tilde{\theta}_{+}(\mathbf{x})=M\theta(\mathbf{x})/\pi$ or
$\tilde{\theta}_{-}(\mathbf{x})=-M\theta(\mathbf{x})/\pi$. This can then be
classically mapped to an estimate for $f(\mathbf{x})$ using

$\tilde{f}(\mathbf{x})=\sin^{2}\left(\frac{\pi\tilde{\theta}_{\pm}(\mathbf{x})}{M}\right).$
(4)

In order to approximate $\tilde{f}(\mathbf{x})$ within additive error
$\epsilon_{p}$ with high probability, QAE requires
$M=\mathcal{O}\left(1/\epsilon_{p}\right)$ invocations to
$\mathcal{Q}(\mathbf{x})$. We measure the query complexity of the QAE
algorithm by counting the number of calls to the oracles
$\mathcal{A}(\mathbf{x})$ and $\mathcal{A}^{\dagger}(\mathbf{x})$ and this
scales as $\mathcal{O}\left(2/\epsilon_{p}\right)$. In the case of financial
derivatives, the function $f$ is the expectation value of the payoff of a
derivative contract. We explore this context in the following section.

##  II Gradient Methods for Financial Derivatives

In this section we describe several methods for computing the gradients of
financial derivatives. The complexities of these approaches are summarized in
Table 1.

Classical Finite Difference | CFD with CRN | Semi-classical | Jordan’s | GAW Quantum Gradient
---|---|---|---|---
$\mathcal{O}\left(k/\epsilon^{3}\right)$ | $\mathcal{O}\left(k/\epsilon^{2}\right)$ | $\mathcal{O}\left(k/\epsilon\right)$ | $\mathcal{O}\left(\sqrt{k}/\epsilon^{2}\right)$ | $\mathcal{O}\left(\sqrt{k}/\epsilon\right)$
Table 1: Summary of the complexity scaling of the different market risk
algorithms. Here $\epsilon$ is the absolute error in gradient estimation, and
$k$ is the dimension of the gradient. For the CFD and semi-classical methods,
which depend additionally on the discretization step $h$, we assume the
optimal choices, which are derived in Section II.2.

###  II.1 Greeks

Let $\mathbf{S}^{t}\in\mathbb{R}_{+}^{d}$ be a vector of values for $d$
underlying assets at time $t$. Let
$(\mathbf{S}^{0},...,\mathbf{S}^{T})=\bar{\omega}\in\bar{\Omega}$ be a path of
a discrete time multivariate stochastic process describing the values of those
assets. We use both notations for a path in the text. The corresponding
probability density function is denoted by $\bar{p}(\bar{\omega})$. Let
$g(\bar{\omega})=g\left(\mathbf{S}^{0},...,\mathbf{S}^{T}\right)\in\mathbb{R}$
be the discounted payoff of some derivative on those assets. To price the
derivative we calculate

$\displaystyle\mathbb{E}(g)=\int_{\bar{\omega}\in\bar{\Omega}}\bar{p}(\bar{\omega})g(\bar{\omega})d\bar{\omega}\approx\sum_{\omega\in\Omega}p(\omega)g(\omega),$
(5)

where we have removed the bar notation to indicate that we have switched from
a continuous to a discrete model of prices. Gradients of this price are known
as greeks.

######  Example 1 (Delta)

Consider a single underlying ($d=1$). We define the gradient of the underlying
with respect to the spot price $S^{0}$ as

$\displaystyle\Delta=\frac{\partial\mathbb{E}(g)}{\partial S^{0}}.$ (6)

Other commonly used greeks are gradients of the price with respect to the time
$T$ (_theta_), the volatility of the underlying model for $S$ (_vega_), the
correlation between assets, or other parameters. In general one calculates
$\mathbb{E}(g,\mathbf{x})$ for some model and/or market parameters
$\mathbf{x}$ and then wishes to compute the set
$\\{\frac{\partial\mathbb{E}(g)}{\partial\theta_{i=1...k}}\\}$. Remarkably,
$k$ can be on the order of hundreds or thousands in practical cases. For this
reason, a $k$ scaling improvement represents an important advantage in this
context.

###  II.2 Classical Finite Difference

Let $f(\mathbf{x})=\mathbb{E}(g,\mathbf{x})$ be the pricing function for a
fixed payoff $g$. Classically, we can compute the gradients of a function
$f:\mathbb{R}^{k}\rightarrow\mathbb{R}$ using finite-difference methods by
sampling the function $f$ over a sufficiently small region $h$ so that
expanding $f$ to first order gives a good approximation to
$f(\mathbf{x})\approx f(\mathbf{a})+(\mathbf{x}-\mathbf{a})\cdot\nabla f$. The
simplest _forward_ finite-difference classical scheme to compute the gradients
at a point $\mathbf{x}_{0}$ requires $k+1$ evaluations of $f$, one at
$\mathbf{x}_{0}$ and $k$ evaluations displaced from $\mathbf{x}_{0}$ by $h$ in
each dimension. The gradient is then approximated using $\partial f/\partial
x_{i}\approx[f(\mathbf{x}_{0}+h\hat{e}_{i})-f(\mathbf{x}_{0})]/h+\mathcal{O}(h)$,
where $\hat{e}_{i}$ is the $i$th normalized basis vector. For a more accurate
approximation, we can instead use a second-order scheme which requires $2k$
function evaluations, and approximate the gradients with

$\partial f/\partial
x_{i}\approx\frac{f(\mathbf{x}_{0}+\frac{h}{2}\hat{e}_{i})-f(\mathbf{x}_{0}-\frac{h}{2}\hat{e}_{i})}{h}+\mathcal{O}(h^{2}).$
(7)

Now suppose that the function $f(x)$ is evaluated with a finite accuracy
$\delta\geq 0$. The associated error for the forward or central finite
differences formulas reads

$\displaystyle\frac{f(x+h)-f(x-(p-1)h)+\mathcal{O}(\delta)}{ph}+\mathcal{O}(h^{p})$
$\displaystyle=$ $\displaystyle\partial_{x}f(x)+\mathcal{O}(\delta/h+h^{p}),$
(8)

where $p=1$ for forward and $p=2$ for central finite differences. Suppose now
that we want to achieve an overall estimation error of $\epsilon>0$, i.e.

$\displaystyle\delta/h+h^{p}\leq\epsilon,$ (9)

where we focus on the asymptotic scaling. While the step size $h$ can be
freely chosen, improving the accuracy $\delta$ is usually related to
increasing computational costs. Thus, we want to maximize $\delta$ by setting
$h$ while still achieving the target accuracy $\epsilon$. This leads to

$\displaystyle\delta=(\epsilon-h^{p})h=\epsilon h-h^{p+1}.$ (10)

Setting the first derivative of the right-hand-side with respect to $h$ to
zero leads to

$\displaystyle h=(\epsilon/(p+1))^{(1/p)},$ (11)

which leads to the optimal $\delta$ for a target $\epsilon$ given by

$\displaystyle\delta$ $\displaystyle=$
$\displaystyle\mathcal{O}(\epsilon^{1+1/p}).$ (12)

When $f(x)$ is approximated by algorithms relying on sampling, the number of
samples for a target approximation error scales as
$\delta=\mathcal{O}(1/M^{q})$, where $M$ denotes here the number of samples
and $q$ depends on the convergence rate of the algorithm i.e., $q=1/2$ for
classical Monte Carlo simulation and $q=1$ for QAE. Then, combining everything
together leads to

$\displaystyle M=\mathcal{O}(\epsilon^{-(1+1/p)/q}).$ (13)

Thus, using $p=2$ and $q=1/2$, we calculate that the complexity of computing
$k$ greeks using a central-difference method with Monte Carlo is
$\mathcal{O}(k/\epsilon^{3})$. We call this approach the  _classical finite
difference_ (CFD) method.

###  II.3 Finite-Difference with Common Random Numbers

Another approach for computing greeks, when derivative pricing is done
classically with Monte Carlo, is to use the second-order central-difference
method of Eq. (7), but to perform correlated sampling by using the same random
numbers in the Monte Carlo evaluation of both
$f(\mathbf{x}_{0}+\frac{h}{2}\hat{e}_{i})$ and
$f(\mathbf{x}_{0}-\frac{h}{2}\hat{e}_{i})$. This way the statistical
fluctuations present in both terms cancel out, effectively removing the
overall error dependence on the discretization step $h$. The complexity of
evaluating $k$ greeks in this case is $\mathcal{O}(k/\epsilon^{2})$, which is
also classically optimal Glasserman and Yao (1992). We call this method the
_classical finite difference with common random numbers_ (CFD-CRN).

###  II.4 Semi-classical Quantum Gradients

A straightforward approach to improve the finite-difference method using
quantum computation is to use a central-difference formula with quantum
amplitude estimation algorithm for each pricing step Rebentrost  _et al._
(2018); Stamatopoulos  _et al._ (2020); Chakrabarti  _et al._ (2021). This
_semi-classical quantum gradient_ (SQG) method Gilyén  _et al._ (2019) then
scales as $\mathcal{O}(k/\epsilon^{1.5})$, which we get by substituting $p=2$
and $q=1$ in Eq. (13). This approach can be improved upon in a similar manner
to using CRN for classical finite difference. Instead of computing
$\mathbb{E}\left(f(\mathbf{x}_{0}+\frac{h}{2}\hat{e}_{i})\right)$ and
$\mathbb{E}\left(f(\mathbf{x}_{0}-\frac{h}{2}\hat{e}_{i})\right)$ separately
using amplitude estimation, we compute
$\mathbb{E}\left((f(\mathbf{x}_{0}+\frac{h}{2}\hat{e}_{i})-f_{\omega}(\mathbf{x}_{0}-\frac{h}{2}\hat{e}_{i}))/h\right)$.
We can do this by computing
$\ket{(f_{\omega}(\mathbf{x}_{0}+\frac{h}{2}\hat{e}_{i})-f_{\omega}(\mathbf{x}_{0}-\frac{h}{2}\hat{e}_{i}))/h}$,
where $f_{\omega}$ denotes the payoff for each path $\omega$ (cf. Sec. II.1)
and then using amplitude estimation on this quantity. In this case the output
of amplitude estimation is

$\sum_{\omega\in\Omega}p(\omega)\left(f_{\omega}(\mathbf{x}_{0}+(h/2)\hat{e}_{i})-f_{\omega}(\mathbf{x}_{0}-(h/2)\hat{e}_{i})\right)/h\approx\mathbb{E}(\partial_{i}f(\mathbf{x})),$
(14)

giving us the expectation value of the finite difference approximation of the
$i$th derivative. The advantage of this method is that there are no separate
statistical fluctuations associated individually with
$f_{\omega}(\mathbf{x}_{0}\pm(h/2)\hat{e}_{i})$ and we recover the standard
$\mathcal{O}(1/\epsilon)$ scaling of amplitude estimation. For $k$ gradients,
we then get an overall complexity of $\mathcal{O}(k/\epsilon)$.

###  II.5 Quantum Gradients

There are other quantum approaches that use a quantum Fourier transform to
compute the gradient with an improved scaling in the dimension $k$ Jordan
(2005); Gilyén  _et al._ (2019). These algorithms require access to a
fractional _phase oracle_ for the target function
$f:\mathbb{R}^{k}\mapsto\mathbb{R}$, in our case the pricing function of Sec.
I.1. This oracle, given a point $\mathbf{x}\in\mathbb{R}^{k}$ and $S>0$,
performs the operation

$O_{Sf}:\ket{\mathbf{x}}\rightarrow e^{2\pi iSf(\mathbf{x})}\ket{\mathbf{x}}.$
(15)

In order to estimate the $k$-dimensional gradient of $f$ at point
$\mathbf{x}_{0}$, the oracle is evaluated over a uniform superposition of $N$
points $\bm{\delta}$ in a sufficiently small $k$-dimensional discretized
hypercube $G_{\mathbf{x}_{0}}^{k}$ of edge length $l$ around $\mathbf{x}_{0}$,
such that

$\frac{1}{\sqrt{N^{k}}}\sum_{\bm{\delta}\in
G_{\mathbf{x}_{0}}^{k}}O_{Sf}\ket{\bm{\delta}}=\frac{1}{\sqrt{N^{k}}}\sum_{\bm{\delta}\in
G_{\mathbf{x}_{0}}^{k}}e^{2\pi
iSf(\mathbf{x_{0}}+\bm{\delta})}\ket{\bm{\delta}}\approx e^{2\pi
iSf(\mathbf{x}_{0})}\frac{1}{\sqrt{N^{k}}}\sum_{\bm{\delta}\in
G_{\mathbf{x}_{0}}^{k}}e^{2\pi iS\cdot\nabla
f_{\mathbf{x}_{0}}\cdot\bm{\delta}}\ket{\bm{\delta}},$ (16)

assuming $f(\mathbf{x}_{0}+\bm{\delta})\approx f(\mathbf{x}_{0})+\nabla
f_{\mathbf{x}_{0}}\cdot\bm{\delta}$ for $||\bm{\delta}||\ll 1$. Applying the
$k$-dimensional inverse Quantum Fourier Transform on the resulting state with
a appropriately chosen value of $S$ and measuring in the computational basis
gives us an approximation of $\nabla f_{\mathbf{x}_{0}}$ with accuracy
$\epsilon=\mathcal{O}(1/N)$ and high probability. The complexity of these
gradient estimation algorithms depends on the resources required to construct
the oracle of Eq. (15) and the size of $G_{\mathbf{x}_{0}}^{k}$ required to
make the approximation in Eq. (16) sufficiently accurate.

##  III Higher-order methods for quantum gradients

Gilyen et al Gilyén  _et al._ (2019) showed that in order to estimate $k$
gradients of a function $f$ with accuracy $\epsilon$, the oracle in Eq. (15)
has to be evaluated $S=N/l\sim D_{2}\sqrt{k}/\epsilon^{2}$ times, where
$D_{2}$ is an upper bound on the magnitude of the second-order derivatives of
$f$. While this gives a quadratically improved scaling in $k$, it does not
improve on the error scaling of the semi-classical method. To improve also the
scaling in accuracy $\epsilon$, Ref. Gilyén  _et al._ (2019) introduces
higher-degree central-difference schemes.

The $2m$-point central-difference approximation for the gradient of $f$ at
$\mathbf{0}$ is given by

$\nabla
f_{(2m)}(\mathbf{0})=\sum_{\ell=-m}^{m}a_{\ell}^{(2m)}f(\ell\mathbf{x}),$ (17)

where the coefficients $a_{\ell}^{(2m)}$ depend on the choice of $m$ and for
uniform grid spacing we have $a_{\ell}^{(2m)}=-a_{-\ell}^{(2m)}$ Fornberg
(1988). A phase oracle for the general $2m$-point scheme of Eq. (17) can be
constructed by composing individual fractional phase oracles for each of the
$2m$ terms, scaled by the appropriate coefficient $a_{\ell}^{(2m)}$. This
leads to a family of algorithms with different phase oracles at different
orders $m$

$O_{Sf}^{m}:\ket{\mathbf{x}}\rightarrow e^{2\pi
iS\sum_{\ell=-m}^{m}a_{\ell}^{(2m)}f(\ell\mathbf{x})}\ket{\mathbf{x}},$ (18)

with $S=N/l$. For example, the phase oracle for the two-point approximation

$O_{f}^{1}\ket{\mathbf{x}}=e^{2\pi
i(f(\mathbf{x})-f(\mathbf{-x}))/2}\ket{\mathbf{x}},$ (19)

can be constructed as the product of oracles $O_{f}^{+}{\mathbf{x}}=e^{\pi
if(\mathbf{x})}\ket{\mathbf{x}}$ and $O_{f}^{-}{\mathbf{x}}=e^{-\pi
if(-\mathbf{x})}\ket{\mathbf{x}}$.

With these higher-order methods, they show that estimating the gradients of a
class of smooth functions111Also known as Gevrey class $G^{\frac{1}{2}}$
functions Gevrey (1918) to accuracy $\epsilon$ using Jordan’s algorithm, can
instead be done with $S=\mathcal{O}(\sqrt{k}/\epsilon)$ phase oracle
applications by picking a large enough value for $m$ in the central-difference
approximation used. Therefore, for this family of smooth functions, gradient
estimation using high-order central-difference methods scales quadratically
better in the desired accuracy $\epsilon$ compared to the original Jordan’s
algorithm. This is described formally in the following Theorem:

######  Theorem 1 (5.4 from Gilyén  _et al._ (2019))

Let $\mathbf{x}\in\mathbb{R}^{k},\epsilon<c\in\mathbb{R}_{+}$ be fixed
constants and suppose $f:\mathbb{R}^{k}\mapsto\mathbb{R}$ is analytic and
satisfies the following: for every $j\in\mathbb{N}$ and $\alpha\in[k]^{j}$

$|\partial_{\alpha}f(\mathbf{x})|\leq c^{j}j^{\frac{j}{2}}.$ (20)

Using the GAW Algorithm and setting $m=\log(c\sqrt{k}/\epsilon)$ and rescaling
$R=\Theta(cm\sqrt{k})$ we can compute an $\epsilon$-approximate gradient
$\tilde{\nabla}f(\mathbf{x})\in\mathbb{R}^{k}$ such that

$\|\nabla f(\mathbf{x})-\tilde{\nabla}f(\mathbf{x})\|_{\infty}\leq\epsilon,$
(21)

with probability at least $1-\delta$, using
$\widetilde{\mathcal{O}}\left(\frac{c\sqrt{k}}{\epsilon}\log\left(\frac{k}{\delta}\right)\right)$
queries to a probability or (fractional) phase oracle of $f$.

In addition to the desired $\mathcal{O}(\sqrt{k}/\epsilon)$ scaling, the
complexity of the GAW algorithm in Theorem 1 includes a factor of
$\log(k/\delta)$ stemming from the fact that we need to extract the medians
after the application of the quantum gradient estimation algorithm, with
probability of $1-\delta$. While in the following sections we focus on the
$\mathcal{O}(\sqrt{k}/\epsilon)$ factor in order to establish the dominant
scaling of the complexity with respect to $k$, in Sec. VI we show how we can
eliminate the $\log(k/\delta)$ factor by using classical maximum likelihood
estimation. This approach not only decreases the overall computational
complexity of the algorithm, but also provides us with concrete confidence
intervals for a given confidence level.

We summarize in Table 1 the scaling in $k$ and $\epsilon$ of the algorithms
discussed in this Section.

###  III.1 Creating phase oracles from probability oracles

Because the function whose gradients we would like to compute is only
accessible through a _probability oracle_ in the form of Eq. (1), we need to
create a corresponding _phase oracle_ of Eq. (15) in order to use the quantum
gradient method described in this section. To do so, we can use the block
encoding technique from Ref. Gilyén  _et al._ (2019):

######  Definition 1 (4.4 from Gilyén  _et al._ (2019))

Suppose that $A$ is an operator on a Hilbert space $\mathcal{H}$, then we say
that the unitary $U$ acting on $\mathcal{H}_{\text{aux}}\otimes\mathcal{H}$ is
a block-encoding of $A$ if

$A=(\bra{\vec{0}}\otimes\mathbb{I})U(\ket{\vec{0}}\otimes\mathbb{I}).$

Intuitively, the block-encoding $H$ is a unitary whose top-left block contains
$A$.

$U=\begin{bmatrix}\mathcal{A}&\text{ . }\\\ .&\text{ . }\end{bmatrix}$

For the probability oracle $\mathcal{A}(\mathbf{x})$ of Eq. (1), the authors
of Ref. Gilyén  _et al._ (2019) observe that

$(\bra{\vec{0}}\otimes\mathbb{I})(\mathcal{A}^{\dagger}(\mathbf{x})(Z\otimes\mathbb{I})\mathcal{A}(\mathbf{x}))(\ket{\vec{0}}\otimes\mathbb{I})=\text{diag}(1-2f(\mathbf{x})),$
(22)

and from Definition 1,
$U\equiv\mathcal{A}^{\dagger}(\mathbf{x})(Z\otimes\mathbb{I})\mathcal{A}(\mathbf{x})$
is a block-encoding of a diagonal matrix $H$ with diagonal entries
$(1-2f(\mathbf{x}))$. With access to this block-encoding, the Hamiltonian
simulation method from Low and Chuang (2019); Gilyén  _et al._ (2019) allows
us to implement an $\epsilon_{\textsubscript{phase}}$-approximation of the
unitary $e^{itH}$, through repeated applications of $U$ and $U^{\dagger}$,
which scales as $\mathcal{O}(|t|+\ln(1/\epsilon_{\textsubscript{phase}}))$.
Therefore, for appropriately chosen values of $t$, we can use block-encoding
and Hamiltonian simulation to produce an
$\epsilon_{\textsubscript{phase}}$-approximation of the phase oracle $O_{Sf}$
from Eq. (15). Note that while the block-encoding of Eq. (22) allows us to
create a phase oracle with phase of $1-2f(\mathbf{x})$ instead of
$f(\mathbf{x})$ in Eq. (15), the first term ends up as a global phase which
can be ignored, and the factor of $-2$ can be absorbed in the factor $t$ of
the Hamiltonian simulation.

While references Low and Chuang (2019); Gilyén  _et al._ (2019) describe how
to perform Hamiltonian simulation and the resources required to realize the
necessary unitary evolution, the methods presented therein require post-
selection, which needs to be factored in for an end-to-end resource
estimation. On the other hand, Ref. Martyn  _et al._ (2021) introduces a
coherent Hamiltonian simulation method which does not require post-selection
and instead succeeds with arbitrarily high probability $1-\delta$, scaling as
$\mathcal{O}(|t|+\ln(1/\epsilon_{\textsubscript{phase}})+\ln(1/\delta))$. For
a target $t$, approximation error $\epsilon_{\textsubscript{phase}}$ and
$\delta=2\epsilon_{\textsubscript{phase}}$, this coherent Hamiltonian
simulation algorithm queries $U$ and its inverse a total number of times

$N_{U}(t,\epsilon_{\textsubscript{phase}},\beta)=2\left\lfloor\frac{1}{2}r\left(\frac{e|t|}{2\beta},\frac{5\epsilon_{\textsubscript{phase}}}{24}\right)\right\rfloor+\gamma\left(\frac{\epsilon_{\textsubscript{phase}}}{3},1-\beta\right)+1,$
(23)

where

  * •

$\beta\in(0,1)$ is a user-chosen parameter,

  * •

$r(\tau,\epsilon)=|\tau|e^{W(\ln(1/\epsilon)/|\tau|)}$, where $W(x)$ is the
Lambert-$W$ function,

  * •

$\gamma(\epsilon,\Delta)=2\cdot\left\lceil\max\left(\frac{e}{\Delta}\sqrt{W\left(\frac{8}{\pi\epsilon^{2}}\right)W\left(\frac{512}{e^{2}\pi\epsilon^{2}}\right)},\sqrt{2}W\left(\frac{8\sqrt{2}}{\sqrt{\pi}\Delta\epsilon}\sqrt{W\left(\frac{8}{\pi\epsilon^{2}}\right)}\right)\right)\right\rceil$
\+ 1.

Using this method, the total number of oracle (Eq. 1) calls required to
construct an $\epsilon_{\textsubscript{phase}}$-approximation of the $m$-order
phase oracle of Eq. (18) is then given by

$N_{o}=\sum_{\ell=-m}^{m}N_{U}\left(2\pi\frac{N}{l}|a_{\ell}^{(2m)}|,\epsilon_{\textsubscript{phase}},\beta\right).$
(24)

The optimal value of $\beta$ depends on the target approximation error
$\epsilon_{\textsubscript{phase}}$. For the cases studied in this manuscript,
we find that $\beta=0.5$ is the optimal choice and we fix that value whenever
we use Eq. (24).

##  IV Resource Estimation of Quantum Gradient Methods

In this section, we perform an asymptotic resource estimation for the gradient
methods described previously. We choose representative parameters from the
financial domain for gradient estimation problems, targeting $k=1000$ greeks
and an approximation error of $\epsilon=10^{-3}$. For the resource estimation
of the GAW method, we assume that the smoothness conditions of Theorem 1, with
smoothness parameter $c=1$, apply for the problem at hand and that the
Hamiltonian simulation phase error is
$\epsilon_{\textsubscript{phase}}=10^{-4}$. Here we choose $c=1$ in order to
estimate the possible usefulness of the algorithm in a best-case scenario of
smoothness from Theorem 1. While we treat the phase error
$\epsilon_{\textsubscript{phase}}$ as a free parameter at this point, our
particular choice of $10^{-4}$ is motivated by numerical simulations which
show this to be a good choice. We discuss the numerical simulations and the
impact of the phase error in more detail in subsequent sections. Then, we set
the finite-difference approximation degree to

$m=\log(c\sqrt{k}/\epsilon)$ (25)

and the spacing parameter to222See proof of Theorem 5.4 from Gilyén  _et al._
(2019)

$l^{-1}=9cm\sqrt{k}\left(81\times 8\times 42\pi
cm\sqrt{k}/\epsilon\right)^{1/(2m)}.$ (26)

Using the proofs from Gilyén  _et al._ (2019) we can estimate the number of
oracle calls that will be required to achieve the target error. Notice that,
in these estimates, we consider only the asymptotic scaling (assuming
remaining constant factors are 1) and assume that
$\epsilon_{\textsubscript{phase}}$ is sufficiently small to have no impact on
performance. This results in $6.3\times 10^{7}$ oracle calls.

We can compare this asymptotically to the performance of both classical finite
difference (with common random numbers and without) and semi-classical quantum
finite difference. These results are summarized in Table 2. We choose the
optimal values of the discretization step $h$ as calculated in Section II.2 to
minimize the gradient estimation error $\epsilon$, and pick $\epsilon=10^{-3}$
for these benchmarks.

Classical Finite Difference | CFD with CRN | Semi-classical | GAW Quantum Gradient
---|---|---|---
$10^{12}$ | $10^{9}$ | $10^{6}$ | $10^{7}$
Table 2: Table of asymptotic oracle calls required to compute quantum
gradients for financially relevant greeks. Parameters for this benchmark are
described in the text. This asymptotic analysis indicates that quantum
gradient method has potential to outperform classical finite difference
methods but is on par with semi-classical methods.

In Table 2, we notice that using the parameters from the proofs in Gilyén  _et
al._ (2019), the GAW algorithm fails to deliver an advantage compared to the
semi-classical method for $k$ as high as $10^{3}$. However, while this
analysis gives us an idea of how the different methods compare in theory, the
estimates are based on asymptotic bounds with parameters which may be loose in
practice and are highly dependent on the smoothness of the functions
considered. In the following section we study the performance of the GAW
method numerically on small examples that are representative of some practical
cases in finance.

###  IV.1 GAW Numerical estimates

The GAW method described in Sec. III gives a quantum algorithm for gradient
estimation which scales as $\mathcal{O}(\sqrt{k}/\epsilon)$ when the target
function satisfies the conditions of Theorem 1. However, in most relevant
financial models of interest, we do not have access to closed-form solutions
that we can examine to check whether they satisfy the smoothness conditions
required of the theorem. As such, in this section we numerically examine the
behavior of the high-order methods of the GAW gradient estimation algorithm
for two financial use cases: (a) A simple (_vanilla_) European call option for
which we have an analytical closed form solution and can benchmark the
performance of the algorithm against the exact gradients, and (b) a path-
dependent basket option with a  _knock-in_ feature which has no known
analytical solution and is in practice classically evaluated using Monte Carlo
methods. In particular, we examine the central-difference approximation order
$m$ and spacing $l$ required for adequately precise gradient estimation and
compare it to the theoretical values of Eq. (25) and Eq. (26) respectively. We
focus on these two parameters because they determine the overall complexity of
the algorithm, which we can then compare to the other methods of Table 2.

Because the resulting quantum circuits are prohibitively wide and deep for
numerical simulation in practice, we adopt the following practical method to
emulate the algorithm’s performance: to estimate $k$ greeks using $n$ bits of
precision we initialize a $k*2^{n}$-dimensional array with the amplitudes of
Eq. (18) computed classically for the chosen derivative order $m$. We then
perform a $k*2^{n}$-dimensional classical inverse Fourier transform of the
array to get the resulting probability distribution which is the output of the
GAW algorithm before measurement. In order to account for the phase error
$\epsilon_{\textsubscript{phase}}>0$ from the Hamiltonian simulation, we add a
random number to each encoded phase in Eq. (18), uniformly picked from the
interval
$[-\epsilon_{\textsubscript{phase}},\epsilon_{\textsubscript{phase}}]$.

####  IV.1.1 Vanilla Options

The simplest example in derivative pricing is a European call option whose
price depends on the performance of a single asset at a pre-determined future
time (the _expiration_ date), relative to a reference price (the _strike_). In
the Black-Scholes-Merton model Black and Scholes (1973), where the asset
undergoes Geometric Brownian Motion (GBM), a call option on a non dividend-
paying asset has a closed form solution given by

$C=SN(d_{1})-Ke^{-rT}N(d_{2}),\quad\text{ with
}d_{1}=\frac{\ln(S/K)+(r+\sigma^{2}/2)T}{\sigma\sqrt{T}},\;d_{2}=d_{1}-\sigma\sqrt{T},$
(27)

where $S$ is the asset price today, $K$ the strike of the option, $r$ is a
risk-free rate of return, $\sigma$ the annualized volatility of the asset, $T$
is the time until the option’s expiration date and $N(x)$ denotes the CDF of
the standard normal distribution. We test the quantum gradient estimation
algorithm for the four greeks of this option/model:
$\textit{delta}=\partial{C}/{\partial S}$,
$\textit{rho}=\partial{C}/\partial{r}$,
$\textit{vega}=\partial{C}/{\partial{\sigma}}$, and
$\textit{theta}=\partial{C}/{\partial T}$. We numerically simulate the GAW
algorithm for increasing $k$ (the number of greeks we compute simultaneously)
and central-difference approximation order $m\in[1,4]$. In each case we search
for the largest value of the spacing $l$ in for which the algorithm produces
an estimate $\epsilon$-close to the exact value with probability $\geq 85\%$
for each greek. We target a gradient error of $\epsilon=2\times 10^{-2}$ which
requires $n=\lceil\log(1/\epsilon)\rceil=6$ qubits in each dimension. 333Our
numerical simulations scale exponentially in this parameter as are emulating
the quantum circuit on these qubits. This limits us to this size.

For each value of $k$, we then compare what we numerically find as the optimal
values for $(m,l)$ — that minimize the total number of oracle calls $N_{o}$
from Eq. (24) while maintaining a success probability of $\geq 85\%$ — to
those used in the proof of Theorem 1, given by Eq. (25) and Eq. (26), assuming
the smoothest possible parameter $c=1$. We set the approximation error from
the Hamiltonian simulation to $\epsilon_{\textsubscript{phase}}=10^{-4}$ and
include it as an error source in our numerical simulations. The gradients of
the vanilla option in Eq. (27) are evaluated at the point
$(S,r,\sigma,T)=(99.5,1\%,20\%,0.1)$, with $K=100$, chosen so that the
parameter values are reasonably realistic from a finance point of view, but at
the same time probing a domain where the function is as non-linear as
possible.444Where the function is (approximately) linear then simple finite
difference methods perform well already. In Table 3 we show the results of the
numerical simulation and the corresponding theoretical estimates from Theorem
1 and in Fig. 1 we show the resulting probability distribution from the
quantum gradient estimation algorithm for each greek, for the case $k=4$ in
Table 3.

| Numerical | Theoretical
---|---|---
$k$ | $m$ | $l$ | $N_{o}$ | $m$ | $l$ | $N_{o}$
2 | 1 | 0.65 | 1976 | 5 | 0.0028 | 570592
3 | 3 | 0.65 | 4664 | 5 | 0.0022 | 712008
4 | 3 | 0.58 | 4904 | 5 | 0.0019 | 833296
Table 3: In this table we show (a) numerical estimates of the query complexity
$N_{o}$ (Eq. (24)) and parameter values $(m,l)$ required by the GAW gradient
estimation algorithm in order to estimate $k$ greeks for the vanilla option of
Eq. (27) within $\epsilon=2\times 10^{-2}$ with probability $\geq 85\%$ and
(b) the corresponding values used in the proof of Theorem 1. We notice that
for a vanilla option, we can reduce the central-difference order and increase
the spacing by more than two orders of magnitude compared to the theoretical
values. As such, the query complexity we estimate numerically is $\sim 200$
times smaller than the theoretical.

Interestingly, we notice from Table 3 that the query complexity required to
estimate the gradients of the vanilla option with high probability in
practice, is orders of magnitude smaller than what is expected from the
parameters used in proof of Theorem 1. Because the types of vanilla options
explored in this section are primarily a motivating example of relevance to
finance that are simpler to analyze, we now turn our attention to more complex
derivatives for which gradient estimation is required for business use in
practice.

Figure 1: The probability distribution resulting from the numerical simulation
of the GAW quantum gradient estimation algorithm with $m=3$ and $l=0.58$ for
the four greeks (delta, rho, vega, theta) of the vanilla option of Eq. (27),
and the corresponding exact values shown as dashed vertical lines. Measurement
of the registers corresponding to each greek will result in a value at most
$\epsilon=1/N=2\times 10^{-2}$ away from the exact value with probability
$\geq 85\%$.

####  IV.1.2 Path-dependent Basket Options

In the previous section, we numerically examined the performance of the GAW
algorithm for vanilla options for which we have analytical solutions and we
can benchmark the algorithm’s performance compared to the exact gradients of
the model. We now perform a similar analysis to a path-dependent option on
multiple underlying assets, which has no known closed form solution, and
compare the query complexity from the numerical simulation to the theoretical
complexity from Theorem 1 as well as to that of the SQG method (Sec. II.4)
which calculates gradients using finite-difference using values estimated
using QAE. Similarly to the previous section, we set the approximation error
from the Hamiltonian simulation to $\epsilon_{\textsubscript{phase}}=10^{-4}$.

The option we consider in this section is defined on three underlying assets
undergoing GBM with volatilities
$\sigma_{1}=20\%,\sigma_{2}=20\%,\sigma_{3}=10\%$ and spot prices
$\vec{S}(t=0)=(S_{1}(t=0),S_{2}(t=0),S_{3}(t=0))=(2.0,2.0,2.0)$. The risk free
rate is set to $r=1\%$ and the option expires in $T=3$ years. The weighted sum
of the asset prices $\vec{w}\cdot\vec{S}(t)$ with weights
$\vec{w}=(w_{1},w_{2},w_{3})=(0.5,0.3,0.2)$ is observed on five days
$t^{B}=[T/5*i]$ for $i\in[1,5]$ across the duration of the contract and the
option’s payoff is given by

$f(\vec{S}(t),K,B)=\begin{cases}\max(\vec{w}\cdot\vec{S}(T)-K,0),&\text{if}\
\vec{w}\cdot\vec{S}(t)>B\text{ for any}\ t\in t^{B}\\\
0,&\text{otherwise,}\end{cases}$ (28)

where we set the strike $K=1.0$ and $B=2.5$. This option is said to have a
_knock-in_ feature, because it only pays off (i.e is _knocked-in_) if the
observed weighted sum exceeds the knock-in barrier $B$ at any of the pre-
defined times before the contract’s expiration. We simulate the GAW algorithm
to estimate four gradients of this option contract’s price $V$: the three
deltas $(\partial V/\partial S_{i})$ and the vega with respect to the first
asset $\partial V/\partial\sigma_{1}$. Because this option does not have a
closed-form solution, we price the option using classical Monte Carlo with
$10^{6}$ paths and use finite-difference to compute the expected gradients
which we use as benchmarks for the quantum algorithm. The quantum algorithm is
simulated using $n=4$ qubits for each gradient register, which sets the target
error of the algorithm to $\epsilon\leq 1/2^{4}=0.0625$. Using $m=1$ and
$l=0.25$ the simulated GAW algorithm gives us an $\epsilon$-close estimate
with probability $\geq 85\%$ for each greek. The SQG method requires
$1/\epsilon$ calls to the $\mathcal{Q}$ operator of QAE, and each
$\mathcal{Q}$ includes two calls to the unitary $\mathcal{A}$ of Eq. (1) and
its inverse, therefore the complexity for $k$ greeks is $2k/\epsilon$. Because
with the SQG method we need to compute the payoff twice in order to construct
the finite difference (Eq. (14)), the $\mathcal{A}$ operator will be
approximately twice as large as the regular pricing oracle. In order to
compare the query complexity more accurately with the other quantum methods
methods we thus include a factor of two in the complexity, for a total of
$4k/\epsilon$.

In Table 4 we show the query complexity and parameters from the numerical
simulation of the GAW method for this path-dependent basket option, along with
(a) the asymptotic estimates from Theorem 1, (b) the query complexity of the
SQG method and (c) the query complexity of the CFD and CFD-CRN methods, all
for the same target approximation error. For the CFD and CFD-CRN methods, the
reported query complexity is the total number of Monte Carlo paths required
for the evaluation of the $k$ greeks within the target approximation error
$\epsilon$ with probability $\geq 85\%$, computed numerically. The table also
includes the parameters and resources required for the same calculation using
the _Simulation-Free Quantum Gradient_ (SFQG) algorithm described in the next
section. The resulting probability distribution for all greeks along with the
MC-estimated “true” values is shown in Fig. 2.

Figure 2: The probability distribution resulting from the numerical simulation
of the GAW quantum gradient estimation algorithm with $m=1$ and $l=0.25$ for
the four greeks of the path-dependent basket option of Sec. IV.1.2, along with
the values estimated using classical Monte Carlo (MC) with $10^{6}$ paths
(vertical dashed lines). Measurement gives us estimates for each greek with
error $\epsilon\leq 0.0625$ with probability $\geq 85\%$.

##  V Simulation-Free Quantum Gradient Method

In this section, we construct a second-order accurate quantum gradient
algorithm, corresponding to $m=1$ in Eq. (17). This is higher order than
Jordan’s algorithm, which has unfavorable scaling, but lower order than
required from the analysis in Gilyen et al. to guarantee the
$\mathcal{O}(\sqrt{k}/\epsilon)$ scaling. A benefit is that, in the case of
derivative pricing, we are able to give an explicit construction of the phase
oracle without needing to appeal to block encoding or Hamiltonian simulation.
We call this the _Simulation-Free Quantum Gradient_ (SFQG) method.

To describe this construction, we first show how to turn the derivative
pricing setting of Sec. I.1 into a first-order phase oracle $O^{1}_{Sf}$ and
then how to build a second-order extension to construct $O^{2}_{Sf}$. We then
simulate this algorithm to estimate the greeks of the basket option from Sec.
IV.1.2, allowing us to compare its performance to the corresponding second-
order accurate GAW method.

###  V.1 First-Order Pricing Phase Oracle

In order to apply the quantum gradient algorithms, we need to construct a
phase oracle of the form of Eq. (15). We evaluate this oracle on a
superposition of points $\bm{\delta}$ in a $k$-dimensional hypercube of edge
length $l$, $G^{k}=[-l/2,l/2]^{k}$ around $\mathbf{x}_{0}$, with $l$ chosen
small enough so that $f(\mathbf{x}_{0}+\bm{\delta})$ is approximately linear
in that region. Let $\mathbf{x}\equiv\mathbf{x}_{0}+\bm{\delta}$ and construct
the following unitary $\mathcal{A}_{\theta}(\mathbf{x})$

$\displaystyle\mathcal{A}_{\theta}(\mathbf{x}):\frac{1}{\sqrt{N^{k}}}\sum_{\mathbf{x}}\ket{\mathbf{x}}\ket{0}_{q+1}$
$\displaystyle\xrightarrow{\mathbb{I}^{\otimes
nk}\otimes\mathcal{A}(\mathbf{x})}$
$\displaystyle\frac{1}{\sqrt{N^{k}}}\sum_{\mathbf{x}}\ket{\mathbf{x}}\left(\sqrt{1-f(\mathbf{x})}\ket{\psi_{0}}_{q}\ket{0}+\sqrt{f(\mathbf{x})}\ket{\psi_{1}}_{q}\ket{1}\right)$
(29) $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N^{k}}}\sum_{\mathbf{x}}\ket{\mathbf{x}}\left(\ket{\Psi_{+}(\mathbf{x})}-\ket{\Psi_{-}(\mathbf{x})}\right)$

where $\ket{\Psi_{\pm}(\mathbf{x})}=-ie^{\pm
i\theta(\mathbf{x})}\ket{\psi_{\pm}(\mathbf{x})}/\sqrt{2}$ and
$\theta(\mathbf{x})=\sin^{2}(f(\mathbf{x}))$, similarly to Eq. (3). Now,
define a modified Grover operator
$\mathcal{Q}_{\theta}(\mathbf{x})=\mathcal{A}_{\theta}(\mathbf{x})\tilde{S}_{0}\mathcal{A}_{\theta}^{\dagger}(\mathbf{x})\tilde{S}_{\psi_{0}}$,
where $\tilde{S}_{0}=\mathbb{I}^{\otimes
nk}\otimes\left(\mathbb{I}-2|0\rangle_{q+1}\langle 0|_{q+1}\right)$ and
$\tilde{S}_{\psi_{0}}=\mathbb{I}^{\otimes
nk}\otimes\left(\mathbb{I}-2|\psi_{0}\rangle_{q}|0\rangle\langle
0|\langle\psi_{0}|_{q}\right)$. Let $D$ be an upper bound on
$|\partial\theta(\mathbf{x})/\partial x_{i}|$ for all $i\in[1,k]$ and apply
$\mathcal{Q}_{\theta}^{\pi S}(\mathbf{x})$ to the state in Eq. (29) with
$S=N/Dl$ to get

$\mathcal{Q}_{\theta}^{\pi
S}(\mathbf{x})\mathcal{A}_{\theta}(\mathbf{x}):\frac{1}{\sqrt{N^{k}}}\sum_{\mathbf{x}}\ket{0}_{q+1}\ket{\mathbf{x}}\rightarrow\frac{1}{\sqrt{N^{k}}}\sum_{\mathbf{x}}e^{2\pi
iS\theta(\mathbf{x})}\ket{\mathbf{x}}\ket{\Psi_{+}(\mathbf{x})}-\frac{1}{\sqrt{N^{k}}}\sum_{\mathbf{x}}e^{-2\pi
iS\theta(\mathbf{x})}\ket{\mathbf{x}}\ket{\Psi_{-}(\mathbf{x})}.$ (30)

Note that the application of
$O_{S\theta}(\mathbf{x})=\mathcal{Q}_{\theta}^{\pi
S}(\mathbf{x})\mathcal{A}_{\theta}(\mathbf{x})$ in Eq. (30) is close to the
phase oracle we are looking for in quantum gradient estimation, but (i) for
$\theta$ rather than $f$, and (ii) with a superposition over the positive and
negative phases that we seek.

An inverse Quantum Fourier Transform
$\left(\mathcal{F}^{-1}_{n}\right)^{\otimes k}$ on the first register then
gives us an estimate of the derivatives of $\theta(\mathbf{x})$ at
$\mathbf{x}_{0}$

$\Ket{\frac{N}{D}\frac{\partial\theta}{\partial
x_{1}}}\Ket{\frac{N}{D}\frac{\partial\theta}{\partial
x_{2}}}\dots\Ket{\frac{N}{D}\frac{\partial\theta}{\partial
x_{k}}}\quad\text{or}\quad\Ket{-\frac{N}{D}\frac{\partial\theta}{\partial
x_{1}}}\Ket{-\frac{N}{D}\frac{\partial\theta}{\partial
x_{2}}}\dots\Ket{-\frac{N}{D}\frac{\partial\theta}{\partial x_{k}}}.$ (31)

This method gives us an oracle for $\theta$ and not for $f$ directly. However,
because we know that $\theta$ and $f$ are related through
$\theta(\mathbf{x})=\sin^{2}(f(\mathbf{x}))$, we can compute the derivatives
of $\theta$ and then use the chain rule to get the gradients of $f$:

$\left.\pm\frac{\partial f}{\partial
x_{i}}\right|_{\mathbf{x}_{0}}=\left.\pm\frac{\partial\theta}{\partial
x_{i}}\right|_{\mathbf{x}_{0}}\times\sin\left(2\theta(\mathbf{x}_{0})\right).$
(32)

This requires knowledge of $\theta(\mathbf{x}_{0})$ which can be calculated
separately through standard QAE via the operator of Eq. (1).

The positive and negative cases can be distinguished by adding a dummy
dimension to $f$ with a gradient that has a known sign. For example one could
change
$\left(f:\mathbb{R}^{k}\mapsto\mathbb{R}\right)\mapsto\left(f+0.5x_{k+1}:\mathbb{R}^{k+1}\mapsto\mathbb{R}\right)$.
Inspecting the sign of the $k+1$-th derivative tells us if we are in the
positive or negative case.

###  V.2 Second-Order Pricing Phase Oracle

The simplest high-order extension we can do is to construct an oracle which
encodes the two-point approximation shown in Eq. (19), which corresponds to
$m=1$ in Eq. (17). Because in the QAE setting we are computing the gradients
of $\theta(\mathbf{x})=\sin^{-1}\sqrt{f(\mathbf{x})}$ instead of
$f(\mathbf{x})$, we need to construct an oracle which performs

$O_{\theta}\ket{\mathbf{x}}=e^{2\pi
iS(\theta(\mathbf{x})-\theta(-\mathbf{x}))/2}\ket{\mathbf{x}},$ (33)

where $\theta(-\mathbf{x})=\sin^{-1}\sqrt{f(-\mathbf{x})}$. In Gilyén  _et
al._ (2019), the authors suggest that high-order oracles of this form can be
constructed as the product of two separate oracles with opposite phases, which
in this case would be the QAE gradient oracles $\mathcal{Q}_{\theta}^{\pi
S}(\mathbf{x})$ and $\mathcal{Q}_{\theta}^{\pi S}(-\mathbf{x})$ of Eq. (30).
However, because these oracles encode both the positive and negative phases in
superposition, their product generates the state

$\displaystyle\frac{1}{\sqrt{N^{k}}}\sum_{\mathbf{x}}e^{2\pi
iS(\theta(\mathbf{x})-\theta(-\mathbf{x}))}\ket{\mathbf{x}}\ket{\psi_{1}(\mathbf{x})}+\frac{1}{\sqrt{N^{k}}}\sum_{\mathbf{x}}e^{-2\pi
iS(\theta(\mathbf{x})-\theta(-\mathbf{x}))}\ket{\mathbf{x}}\ket{\psi_{2}(\mathbf{x})}$
(34) $\displaystyle+\frac{1}{\sqrt{N^{k}}}\sum_{\mathbf{x}}e^{2\pi
iS(\theta(\mathbf{x})+\theta(-\mathbf{x}))}\ket{\mathbf{x}}\ket{\psi_{3}(\mathbf{x})}\frac{1}{\sqrt{N^{k}}}\sum_{\mathbf{x}}e^{-2\pi
iS(\theta(\mathbf{x})+\theta(-\mathbf{x}))}\ket{\mathbf{x}}\ket{\psi_{4}(\mathbf{x})}$
(35)

where $\psi_{i}$ denotes products of eigenstates of $\mathcal{Q}(\mathbf{x})$
and $\mathcal{Q}(-\mathbf{x})$ which we can then be ignored for the rest of
the algorithm. While the first two terms contain the appropriate phase
kickback (up to the sign) for the two-point approximation method of Eq. (33),
the last two terms encode a phase proportional to
$\theta(\mathbf{x})+\theta(-\mathbf{x})=\theta(\mathbf{x}_{0})+\mathcal{O}(\partial^{2}\theta(\mathbf{x})/\partial\mathbf{x}^{2})$,
which create a probability peak around zero instead of the gradient. We now
show a way to construct the operator which encodes the phase as required by
Eq. (33) without garbage terms using the methods from Carrera Vazquez and
Woerner (2021) to construct QAE oracles which compute products and linear
combinations of functions.

Firstly, note that the $\mathcal{A}$ operator of QAE in Eq. (1) can also be
written as

$\mathcal{A}(\mathbf{x})\ket{0}_{q+1}=\cos\left({\theta(\mathbf{x})}\right)\ket{\psi_{0}}_{q}\ket{0}+\sin\left({\theta(\mathbf{x})}\right)\ket{\psi_{1}}_{q}\ket{1},$
(36)

for real $f(\mathbf{x})=\sin^{2}(\theta(\mathbf{x}))$. Constructing
$\mathcal{Q}_{\theta}^{\pi S}(\mathbf{x})$ based on $\mathcal{A}(\mathbf{x})$,
we generate a phase kickback of $e^{\pm 2\pi iS\theta(\mathbf{x})}$. An
operator $\tilde{\mathcal{Q}}_{\theta}^{\pi S}(\mathbf{x})$ based on

$\tilde{\mathcal{A}}_{\theta}(\mathbf{x})\ket{0}_{m+n}=\cos\left(\theta(\mathbf{x})-\theta(-\mathbf{x})\right)\ket{\psi_{0}}_{m}\ket{0}_{n}+\sin\left(\theta(\mathbf{x})-\theta(-\mathbf{x})\right)\ket{\psi_{1}}_{m}\ket{1}_{n},$
(37)

for some appropriate choices of $m,n$ could therefore be used to generate the
phase kickback of $e^{\pm 2\pi
iS\left(\theta(\mathbf{x})-\theta(-\mathbf{x})\right)}$, and can hence serve
as the phase oracle of Eq. (33) we want to construct.

In derivative pricing settings, we can construct a more general variant of Eq.
(1)

$\mathcal{A}(\mathbf{x},c)\ket{0}_{q+1}=\sqrt{1-cf(\mathbf{x})}\ket{\psi_{0}}_{q}\ket{0}+\sqrt{cf(\mathbf{x})}\ket{\psi_{1}}_{q}\ket{1},$
(38)

with $c\in\mathbb{R}$, by scaling the payoff of the derivative with $c$, as
long as $cf(\mathbf{x})\leq 1$. The gradient estimation method in this section
requires first the computation of $f(\mathbf{x}_{0})$ (through
$\theta(\mathbf{x}_{0})$) in order to use the chain rule in Eq. (32),
therefore if $f(\mathbf{x}_{0})>0.25$ we normalize $\mathcal{A}(\mathbf{x})$
so that we have $f(\pm\mathbf{x})\leq 0.25$. We start with the state
$\ket{0}_{2q+2}\frac{1}{\sqrt{2}}\left(\ket{0}+\ket{1}\right)$ and control on
the last qubit to perform
$\mathcal{A}(\mathbf{x},4)\otimes\mathcal{A}(-\mathbf{x},1)$ if the qubit is
$\ket{0}$ and $\mathcal{A}(\mathbf{x},1)\otimes\mathcal{A}(-\mathbf{x},4)$
otherwise. The first $\mathcal{A}$ operator acts on the qubits
$[1,2,\dots,q,q+1]$ and the second on qubits $[q+2,q+3,\dots,2q+2]$. The
resulting state is

$\displaystyle\frac{1}{\sqrt{2}}\left(\left(\sqrt{1-4f(\mathbf{x})}\ket{\psi_{0}}_{q}\ket{0}+\sqrt{4f(\mathbf{x})}\ket{\psi_{1}}_{q}\ket{1}\right)\otimes\left(\sqrt{1-f(-\mathbf{x})}\ket{\psi_{0}}_{q}\ket{0}+\sqrt{f(-\mathbf{x})}\ket{\psi_{1}}_{q}\ket{1}\right)\right)\ket{0}$
(39)
$\displaystyle+\frac{1}{\sqrt{2}}\left(\left(\sqrt{1-f(\mathbf{x})}\ket{\psi_{0}}_{q}\ket{0}+\sqrt{f(\mathbf{x})}\ket{\psi_{1}}_{q}\ket{1}\right)\otimes\left(\sqrt{1-4f(-\mathbf{x})}\ket{\psi_{0}}_{q}\ket{0}+\sqrt{4f(-\mathbf{x})}\ket{\psi_{1}}_{q}\ket{1}\right)\right)\ket{1}.$
(40)

The combined state of the first and last two qubits when the first two of them
are in state $\ket{11}$ is

$\frac{1}{\sqrt{2}}\left(\sqrt{4f(\mathbf{x})}\sqrt{1-f(-\mathbf{x})}\ket{11}\ket{0}+\sqrt{1-f(\mathbf{x})}\sqrt{4f(-\mathbf{x})}\ket{11}\ket{1}\right).$
(41)

Applying a Hadamard on the last qubit, the amplitude of these three qubits
being in state $\ket{111}$ is

$\frac{1}{2}\left(\sqrt{4f(\mathbf{x})}\sqrt{1-f(-\mathbf{x})}-\sqrt{1-f(\mathbf{x})}\sqrt{4f(-\mathbf{x})}\right)=\sqrt{f(\mathbf{x})}\sqrt{1-f(-\mathbf{x})}-\sqrt{1-f(\mathbf{x})}\sqrt{f(-\mathbf{x})}=\sin(\theta(\mathbf{x})-\theta(-\mathbf{x})),$
(42)

where we have used $\sin(a-b)=\sin a\cos b-\cos a\sin b$,
$\sin(\theta(\mathbf{x}))=\sqrt{f(\mathbf{x})}$ and
$\cos(\theta(\mathbf{x}))=\sqrt{1-f(\mathbf{x})}$. We have now constructed the
desired operator of Eq. (37), where the $\ket{111}$ state marks the ”good”
state for QAE. A circuit diagram of this method to construct the
$\tilde{\mathcal{A}}_{\theta}(\mathbf{x})$ operator is shown in Fig. 3.

Now that we have a method to construct the operator of Eq. (33), we use it to
numerically simulate the SFQG algorithm to compute the greeks of the basket
option of Sec. IV.1.2. In this case we use the algorithm to compute the
gradients $\partial\theta/\partial S_{1}$, $\partial\theta/\partial S_{2}$,
$\partial\theta/\partial S_{3}$, $\partial\theta/\partial\sigma_{1}$, where
$\theta=\sin^{-1}\sqrt{V}$, and $V$ is the option price, from which we can
then estimate the option’s greeks using Eq. (32), assuming we have already
priced the contract. Similarly to the simulation of the GAW algorithm in Sec.
IV.1.2, we search for the parameter values which estimate the gradients within
$\epsilon\leq 0.0625$ with probability of success $\geq 85\%$. The probability
distribution for all four gradients after the application of the SFQG
algorithm is shown in Fig. 4, and the parameter values and query complexity of
the algorithm is shown in Table 4 along with all the other quantum, classical
and semi-classical methods studied in this manuscript. We notice that the SFQG
method is the best performing for this option contract, with oracle complexity
$\sim 4$x and $\sim 25$x smaller than semi-classical methods and Hamiltonian-
based GAW respectively.

Figure 3: Circuit diagram of the method used to generate the $\tilde{\mathcal{A}}_{\theta}(\mathbf{x})$ operator of Eq. (37). An open circle in controlled operations denotes that the unitary is performed if the state of the control qubit is $\ket{0}$ and a unitary controlled on a qubit with a filled circle taks place if the control qubit is $\ket{1}$. Note that we can combine the two invocations of $\mathcal{A}(\mathbf{x},c)$ by first encoding the value of $c$ in an ancilla register, controlled on the last qubit, and then compute $\mathcal{A}(\mathbf{x},c)$ once in superposition, controlled on the ancilla register. The same applies to the two invocations of $\mathcal{A}(-\mathbf{x},c)$. |  | $N_{o}$ | m | l
---|---|---|---|---
Quantum | Simulation-Free (SFQG) | 64 | 1 | 0.25
Semi-classical (SQG) | 256 | n/a | n/a
GAW (Numerical) | 1600 | 1 | 0.25
GAW (Theoretical) | 201,528 | 4 | 0.0018
|  Classical
---
Finite-Difference w/ CRN (CFD-CRN) | 32,000 | n/a | n/a
|  | Finite-Difference (CFD)
---
400,000 | n/a | n/a
Table 4: Comparison between (a) the numerical estimates of the query
complexity $N_{o}$ (Eq. (24)) and parameter values $(m,l)$ required by the GAW
gradient estimation algorithm in order to estimate $k$ greeks for the basket
option of Eq. (28) within $\epsilon\leq 0.0625$ with probability $\geq 85\%$,
(b) the corresponding theoretical values used in the proof of Theorem 1, (c)
the parameters and resources required for the same calculation using the
Simulation-Free Quantum Gradient (SFQG) algorithm described in Sec. V, (d) the
query complexity of an SQG method for the same target accuracy and confidence
interval and (e) the total number of classical Monte Carlo paths required for
the same accuracy and confidence interval using CFD and CFD-CRN methods. We
numerically find that for this path-dependent basket option, the $m=1$ GAW
method can estimate the four greeks within $\epsilon$ with probability $\geq
85\%$, with $\sim 125$ times smaller query complexity implied by the proof of
Theorem 1 and 20 times smaller than the best finite-difference-based classical
method (CFD-CRN). The SFQG method far outperforms the query complexity of all
other gradient methods studied here. Figure 4: The probability distribution
resulting from the numerical simulation of the second-order accurate,
Simulation-Free Quantum Gradient (SFQG) algorithm described in Sec. V with
$m=1$ and $l=0.25$ for the four greeks of the path-dependent basket option of
Sec. IV.1.2, along with the exact values estimated using classical Monte Carlo
simulation (vertical dashed lines). This algorithm allows us to estimate the
gradients of the function $\theta=\sin^{-1}\sqrt{V}$, where $V$ is the option
price, from which we can then estimate the option’s greeks using Eq. (32).
Similarly to Fig. 2, measurement gives us estimates for each greek within
error $\epsilon\leq 0.0625$ with probability $\geq 85\%$.

##  VI Quantum Gradient Estimation using Maximum Likelihood Estimation

The success of the quantum gradient estimation algorithms studied so far in
this manuscript, as well as those in Jordan (2005); Gilyén  _et al._ (2019) is
judged by their ability to estimate the gradient within $\epsilon$ with high
probability. In Gilyén  _et al._ (2019) the authors suggest repeating the
gradient estimation algorithm $\mathcal{O}(\log(k/\rho))$ times for $k$
gradients and taking the median to get the estimates within error $\epsilon$
with probability at least $(1-\rho)$. However, this approach has two main
drawbacks: a) In practice we would like to be able to characterize the
correctness of the output with precise confidence intervals, and b) the output
can only be one of the $N$ possible values of the discretized hypercube
$G_{\mathbf{x}_{0}}^{k}$ of Eq. (16).

Recently, approaches using Maximum Likelihood Estimation (MLE) have been
proposed to address these issues for amplitude estimation algorithms Suzuki
_et al._ (2020); Tanaka  _et al._ (2021); Grinko  _et al._ (2021), where the
quantum circuits are sampled more than once, and the results are classically
post-processed with MLE. In this section we show that we can apply the MLE
method from Grinko  _et al._ (2021) to the quantum gradient estimation
algorithms to enhance the final estimate.

Given a probability distribution $p$ with unknown parameter $g$ and data
${x_{i}}$ with $i=1,\dots,M$ sampled from it, MLE is used to obtain an
estimate $\hat{g}$ for $g$. This is done by maximizing the log-likelihood
$\log L$

$\hat{g}=\arg\underset{g^{\prime}}{\max}\log
L(g^{\prime})=\arg\underset{g^{\prime}}{\max}\log\left(\prod_{i=1}^{M}p(x_{i}|g^{\prime})\right)=\arg\underset{g^{\prime}}{\max}\left(\sum_{i=1}^{M}\log
p(x_{i}|g^{\prime})\right),$ (43)

which measures how likely it is to measure the data ${x_{i}}$ if $g^{\prime}$
is the true parameter. The general quantum gradient estimation algorithm from
Sec. II.5 requires a region $l$ where the function whose gradients we are
computing is approximately linear. In this case, the probability distribution
after the inverse Quantum Fourier Transform is applied to Eq. (16) is the same
as that of phase estimation and is given by Brassard  _et al._ (2002)

$p(x)=\frac{\sin^{2}\left(N\Delta\pi\right)}{N^{2}\sin^{2}\left(\Delta\pi\right)},$
(44)

where $\Delta=(x-g^{\prime})$ and $N$ is the number of possible measurements.

Confidence intervals for the MLE estimate $\hat{g}$ can be derived using the
likelihood ratio (LR) Koch (1999). Following the analysis in the supplementary
information of Grinko  _et al._ (2021), the confidence interval at the
$(1-\alpha)$ confidence level is the value $g^{\prime}$ satisfying
$\left\\{g^{\prime}\in[0,1]:\log L(g^{\prime})\geq\log
L(\hat{g})-q_{\chi_{1}^{2}}(1-\alpha)/2\right\\}$, where $q_{\chi_{1}^{2}}$
denotes the $(1-\alpha)$ quantile of the $\chi^{2}$ distribution.

In Fig. 5 we show how MLE can be used to improve the estimate of Vega
($\partial\theta/\partial\sigma$) for the basket option from Fig. 4 calculated
with the simulation-free quantum gradient method. First, in Fig. 5a we show
that the final probability distribution after the application of the quantum
gradient estimation algorithm does indeed fit Eq. (44), and thus the QAE with
MLE results from Grinko  _et al._ (2021) can be used here too. For clarity, we
only plot the measurement outcomes in the interval which contains most of the
probability mass. In Fig. 5b we plot the log-likelihood $\log L(g^{\prime})$
as a function of $g^{\prime}$ across the same interval when we sample the
quantum gradient circuit that produces the probability distribution in Fig. 5a
30 times. The MLE estimate $\hat{g}$ is the global maximum of the function
(green dot), which is very close to the true value $g$ (dashed red line) and
significantly better than the median of the final probability distribution
(yellow arrow).

We note that one significant advantage of the MLE method is that because the
final estimate is not constrained to be one of the $N$ possible discrete
values, we can decrease the value of $N$ at the cost of increasing the number
of samples we take, thus lowering the overall width and depth of the quantum
circuit. For instance, while the distance between the $N$ possible discrete
values in Fig. 5a is $0.0625$, the error $\epsilon$ in the MLE estimate in
Fig. 5b is $\epsilon=4\times 10^{-3}$.

(a)

(b)

Figure 5: a) The discrete probability distribution before measurement of Vega
($\partial\theta/\partial\sigma$) for the basket option from Fig. 4 (blue
bars) is fitted to the theoretical distribution $p(x)$ of Eq.(44) (black
line). The maximum of the distribution occurs at the exact analytical value
(dashed red line). b) The global maximum of the log-likelihood $\log
L(g^{\prime})$ (green dot) gives us a better estimate of the true value
(dashed red line) than the most likely result if we sample from the
probability distribution and take the median (yellow triangle). The log-
likelihood plot was produced by sampling the quantum gradient estimation
circuit 30 times and allows us to estimate this greek to within $\epsilon\leq
2\times 10^{-3}$ with confidence level $1-\alpha=0.68$.

##  VII Updated Estimates for Quantum Advantage

The resource estimation of Chakrabarti et al. Chakrabarti  _et al._ (2021)
established that quantum advantage for derivative pricing with respect to
classical Monte Carlo methods could require a quantum processor that can
execute $\mathcal{O}(10^{7})$ T-gates per second at a code distance that can
support $\mathcal{O}(10^{10})$ logical operations. More specifically, pricing
the autocallable contract studied in Ref. Chakrabarti  _et al._ (2021) using
the reparameterization technique introduced therein to within $\epsilon\leq
2\times 10^{-3}$ with confidence level $1-\alpha=0.68$ requires a T-depth of
$5\times 10^{7}$. As such, in order to match the classical Monte Carlo pricing
time estimated as 1 second, a logical quantum clock rate of $50$Mhz would be
needed. We use the same reparameterization technique from Ref. Chakrabarti
_et al._ (2021) to construct oracle size estimates, and estimate the T-depth
of the SFQG quantum gradient circuit required to compute the four greeks of
the basket option of Sec. IV.1.2 and Table 4. We find that in order to
calculate the greeks to within the same $\epsilon\leq 2\times 10^{-3}$ and
confidence level $1-\alpha=0.68$, we can use the SFQG gradient estimation
method with the parameters of Fig. 4 and the maximum likelihood estimation
(MLE) method from Sec. VI using 30 shots (see Fig. 5b). Using these
parameters, we estimate that the total T-depth of the SFQG circuit multiplied
by the required number of shots is $9\times 10^{6}$, the end-to-end circuit
would require 6k logical qubits, and we would need to execute T-gates at a
code distance that can support $10^{8}$ logical operations. Assuming the
contract can be classically priced using Monte Carlo in 1 second to within the
same error and confidence level Chakrabarti  _et al._ (2021) and that the
greeks can be calculated using a second-order finite-difference method applied
to Monte Carlo pricings, the greeks of the basket option can be classically
estimated in 8 seconds (four greeks, two pricings per greek). Therefore,
quantum advantage in calculating the greeks of this derivative contract would
require executing T-gates at a rate of 1Mhz, $50$ times lower than the
estimate of Chakrabarti et al. for quantum advantage in derivative pricing.
Moreover, because the quantum gradient circuit is sampled 30 times when we use
the MLE method (see Fig. 5), if we have parallel access to 30 QPUs on which
the quantum circuit can be loaded and sampled simultaneously, we can achieve
the same runtime as the serial execution if the logical cock rate of each
device is $\sim 30$kHz, on par with current estimates of feasible logical
clock rates around 10kHz Fowler and Gidney (2019).

While the computational cost required for the numerical simulations of the
quantum gradient estimation algorithm limits our analysis to a maximum of
$k=4$ greeks for practically relevant use cases, we expect that the algorithm
can scale favorably to derivative pricing problems of higher dimensionality,
motivated by the fact that most, if not all, derivative contracts of practical
interest have piecewise-linear payoffs Hull (2006). While Theorem 1 provides a
complexity of $\mathcal{O}(\sqrt{k}/\epsilon)$ for the class of functions
considered, it does not preclude higher-order speedups with respect to $k$ for
smoother functions. For instance, the same gradient estimation algorithm
applied to simple polynomial functions can achieve a complexity of
$\mathcal{O}(\log(k)/\epsilon)$ Gilyén  _et al._ (2019). Derivative contracts
often include market parameters with little or no cross-dependence, i.e.
$\partial^{n}f/(\partial x_{1}^{m}\partial x_{2}^{n-m})\sim 0$ for
$n>1,m\in[1,n-1]$. For example, the basket option from Sec. IV.1.2 without the
knock-in feature satisfies $\partial^{n}V/(\partial S_{i}^{n}\partial
S_{j}^{n-m})=0$. The absence of such higher-order terms allows us to pick
values of $m$ and $l$ in the application of the GAW algorithm which lead to
smaller overall oracular cost than what is required by Theorem 1. An upper
bound on the potential advantage compared to classical Monte Carlo is the case
where the function is at most a second-degree polynomial in all $k$ variables,
yielding an overall speedup of $\mathcal{O}(k)$. Therefore, while we have made
a first step in establishing the relevance and promise of the quantum gradient
estimation methods in the context of financial derivative risk analysis, the
possible extent of quantum speedup will be highly dependent on the nature of
the price function for each derivative.

##  VIII Discussion

We introduce a method to compute gradients of financial derivatives (greeks)
using the gradient algorithms from Jordan (2005); Gilyén  _et al._ (2019).
This method suggests additional quantum advantage is possible in risk
analysis, on top of the quadratic speedup of derivative pricing Rebentrost
_et al._ (2018); Stamatopoulos  _et al._ (2020); Chakrabarti  _et al._ (2021).
Classically computing $k$ greeks with finite-difference methods - when the
underlying derivative is priced using Monte Carlo - has complexity
$\mathcal{O}(k/\epsilon^{2})$ and straightforward extension of finite-
difference methods to derivative pricing using amplitude estimation provides a
quadratic advantage with complexity $\mathcal{O}(k/\epsilon)$. In this work,
we explore an additional quadratic advantage for overall complexity of
$\mathcal{O}(\sqrt{k}/\epsilon)$. The gradient estimation algorithm from
Gilyén  _et al._ (2019) guarantees this quadratic advantage with respect to
the number of greeks when the pricing function satisfies the smoothness
conditions of Theorem 1. Because derivative pricing problems of practical
interest in finance involve numerous diverse multivariate price functions and
generally have no analytical solutions, understanding whether and which
financial derivatives satisfy the aforementioned smoothness conditions is a
challenging task. For this reason, we employ numerical methods to simulate the
gradient estimation algorithm for two example derivatives: a) a European call
option which has a closed-form solution and is used to establish the validity
and benchmarks of the algorithm and b) a path-dependent basket option which
has no known analytical solution and is representative of typical derivative
price functions. We find that the quantum gradient algorithms not only succeed
in estimating the associated greeks for these examples with high probability,
but that the resulting query complexity is in fact significantly smaller than
that suggested by Theorem 1 (Tables 3 and 4), suggesting that the associated
price functions are smoother than those studied in Theorem 1.

Another question we tackle in this work is the rigorous resource estimation of
the quantum oracles involved in the quantum gradient estimation algorithm. Due
to the extra cost associated with the block-encoding and Hamiltonian
simulation required to approximately construct the phase oracle of Eq. (15)
from the probability oracle of Eq. (1) used in derivative pricing, in Sec. V.2
we develop a method to construct a cheaper, second-order ($m=1$) phase oracle
exactly by taking advantage of the structure of amplitude estimation. An
interesting question is whether this method can be extended to construct phase
oracles for higher-order ($m>1)$ gradient methods which would apply more
generally to quantum gradient estimation problems.

In Sec. VI we show that it is possible to enhance quantum gradient algorithms
by employing maximum likelihood estimation (MLE), allowing us to determine the
resources required to estimate gradients with precise confidence intervals and
confidence levels. Using this MLE method, in Sec. VII we estimate the
resources required for quantum advantage in derivative market risk for typical
use cases of practical interest. We find that employing quantum gradient
methods in derivative pricing lowers the logical clock rate estimate for
quantum advantage from Chakrabarti et al. Chakrabarti  _et al._ (2021) by a
factor of 50.

While finite-difference methods are still used in practice to compute greeks,
more recently classical automatic differentiation (AD) methods have been
gathering considerable interest because of their ability to significantly
reduce the associated computational costs, at the cost of increased memory
footprint Homescu (2011); Pages  _et al._ (2016). In particular, the adjoint
mode of automatic differentiation (AAD) in certain cases allows the
computation of all $k$ gradients of a scalar function $f$ at a cost which is
independent of $k$, meaning that the overall classical complexity cost in this
case becomes $\mathcal{O}(\omega/\epsilon^{2})$, for some constant $\omega$
depending on the function $f$ Capriotti (2010). While the complexity of the
GAW quantum gradient estimation algorithm scales as
$\mathcal{O}(\sqrt{k}/\epsilon)$ for the class of smooth functions in Theorem
1, in Sec. IV.1 we saw that for practical use cases in finance, the algorithm
scales 100x-200x times better than theoretical estimate from Theorem 1 for a
given $\epsilon$ (Tables 3 and 4). As such, depending on the practical scaling
of the GAW algorithm to larger values of $k$ for finance use cases, it is
possible that it could also outperform the complexity of AAD methods. It is
also interesting to consider whether a similar construct as that employed by
AD can be applied in a quantum setting. In Appendix A we provide such a
construct and show that in certain settings it can lead to similar performance
profile as classical AD, in that the runtime of the algorithm is independent
of the number of greeks at the expense of increased memory usage. A detailed
comparison of the performance between the quantum gradient algorithms and AD
methods is a worthy study on its own and it is left for future research.

Intuitively, the expensive part of quantum derivative pricing is extracting
the result through quantum amplitude estimation. Therefore it is advantageous
to perform additional calculations involving that value before reading it out.
In this work we consider calculating gradients, but there are other associated
risk metrics that are equally of practical interest, such as the computation
of portfolio value-at-risk (VaR). Based on our results, we are therefore
cautiously optimistic that further quantum advantage in the derivative pricing
subroutine can be leveraged at these higher levels of calculation and
aggregation, and suggests an additional research path going forward.

###### Acknowledgements.

We thank Rajiv Krishnakumar, Shouvanik Chakrabarti and Srinivasan Arunachalam
for useful discussions regarding quantum gradient estimation algorithms, and
Paul Burchard, Graham Griffiths, Alex Hurst, Dunstan Marris and Elmer Tan for
their technical and business insights regarding financial derivatives and
market risk.

## References

  * Rebentrost  _et al._ (2018) Patrick Rebentrost, Brajesh Gupt, and Thomas R. Bromley, “Quantum computational finance: Monte carlo pricing of financial derivatives,” Phys. Rev. A 98, 022321 (2018).
  * Woerner and Egger (2019) Stefan Woerner and Daniel J. Egger, “Quantum risk analysis,” npj Quantum Information 5 (2019), 10.1038/s41534-019-0130-6.
  * Egger  _et al._ (2020) Daniel J Egger, Ricardo Garcia Gutierrez, Jordie Cahue Mestre, and Stefan Woerner, “Credit risk analysis using quantum computers,” IEEE Transactions on Computers (2020), 10.1109/TC.2020.3038063.
  * Stamatopoulos  _et al._ (2020) Nikitas Stamatopoulos, Daniel J. Egger, Yue Sun, Christa Zoufal, Raban Iten, Ning Shen, and Stefan Woerner, “Option pricing using quantum computers,” Quantum 4, 291 (2020).
  * Chakrabarti  _et al._ (2021) Shouvanik Chakrabarti, Rajiv Krishnakumar, Guglielmo Mazzola, Nikitas Stamatopoulos, Stefan Woerner, and William J. Zeng, “A threshold for quantum advantage in derivative pricing,” Quantum 5, 463 (2021).
  * Montanaro (2015) Ashley Montanaro, “Quantum speedup of monte carlo methods,” Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 471 (2015), 10.1098/rspa.2015.0301.
  * Hull (2006) John C. Hull,  _Options, futures, and other derivatives_, 6th ed. (Pearson Prentice Hall, Upper Saddle River, NJ [u.a.], 2006).
  * Gilyén  _et al._ (2019) András Gilyén, Srinivasan Arunachalam, and Nathan Wiebe, “Optimizing quantum optimization algorithms via faster quantum gradient computation,” Proceedings of the Thirtieth Annual ACM-SIAM Symposium on Discrete Algorithms , 1425–1444 (2019).
  * Jordan (2005) Stephen P. Jordan, “Fast quantum algorithm for numerical gradient estimation,” Physical Review Letters 95 (2005), 10.1103/physrevlett.95.050501.
  * Chakrabarti  _et al._ (2020) Shouvanik Chakrabarti, Andrew M. Childs, Tongyang Li, and Xiaodi Wu, “Quantum algorithms and lower bounds for convex optimization,” Quantum 4, 221 (2020).
  * Brassard  _et al._ (2002) Gilles Brassard, Peter Hoyer, Michele Mosca, and Alain Tapp, “Quantum Amplitude Amplification and Estimation,” Contemporary Mathematics 305 (2002), 10.1090/conm/305/05215.
  * Glasserman and Yao (1992) P. Glasserman and D. Yao, “Some guidelines and guarantees for common random numbers,” Management Science 38, 884–908 (1992).
  * Fornberg (1988) B. Fornberg, “Generation of finite difference formulas on arbitrarily spaced grids,” Mathematics of Computation 51, 699–706 (1988).
  * Gevrey (1918) Maurice Gevrey, “Sur la nature analytique des solutions des équations aux dérivées partielles. premier mémoire,” Annales scientifiques de l’École Normale Supérieure 3e série, 35, 129–190 (1918).
  * Low and Chuang (2019) Guang Hao Low and Isaac L Chuang, “Hamiltonian simulation by qubitization,” Quantum 3, 163 (2019).
  * Gilyén  _et al._ (2019) András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe, “Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics,” in  _Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_ (2019) pp. 193–204.
  * Martyn  _et al._ (2021) John M Martyn, Yuan Liu, Zachary E Chin, and Isaac L Chuang, “Efficient fully-coherent hamiltonian simulation,” (2021), arXiv:2110.11327 [quant-ph] .
  * Black and Scholes (1973) Fischer Black and Myron Scholes, “The pricing of options and corporate liabilities,” Journal of Political Economy 81, 637–654 (1973).
  * Carrera Vazquez and Woerner (2021) Almudena Carrera Vazquez and Stefan Woerner, “Efficient state preparation for quantum amplitude estimation,” Physical Review Applied 15 (2021), 10.1103/physrevapplied.15.034027.
  * Suzuki  _et al._ (2020) Yohichi Suzuki, Shumpei Uno, Rudy Raymond, Tomoki Tanaka, Tamiya Onodera, and Naoki Yamamoto, “Amplitude estimation without phase estimation,” Quantum Information Processing 19, 75 (2020).
  * Tanaka  _et al._ (2021) Tomoki Tanaka, Yohichi Suzuki, Shumpei Uno, Rudy Raymond, Tamiya Onodera, and Naoki Yamamoto, “Amplitude estimation via maximum likelihood on noisy quantum computer,” Quantum Information Processing 20, 293 (2021).
  * Grinko  _et al._ (2021) Dmitry Grinko, Julien Gacon, Christa Zoufal, and Stefan Woerner, “Iterative quantum amplitude estimation,” npj Quantum Information 7 (2021), 10.1038/s41534-021-00379-1.
  * Koch (1999) Karl-Rudolf Koch,  _Parameter Estimation and Hypothesis Testing in Linear Models_ (Springer-Verlag Berlin Heidelberg, 1999).
  * Fowler and Gidney (2019) Austin G. Fowler and Craig Gidney, “Low overhead quantum computation using lattice surgery,” (2019), arXiv:1808.06709 [quant-ph] .
  * Homescu (2011) Cristian Homescu, “Adjoints and automatic (algorithmic) differentiation in computational finance,” Risk Management eJournal (2011), 10.2139/ssrn.1828503.
  * Pages  _et al._ (2016) Gilles Pages, Olivier Pironneau, and Guillaume Sall, “Vibrato and automatic differentiation for high order derivatives and sensitivities of financial options,” Journal of Computational Finance 22 (2016), 10.21314/JCF.2018.350.
  * Capriotti (2010) Luca Capriotti, “Fast greeks by algorithmic differentiation,” J. Comput. Financ. 14 (2010), 10.2139/ssrn.1619626.
  * Capriotti and Giles (2010) L. Capriotti and M. Giles, “Fast correlation greeks by adjoint algorithmic differentiation,” ERN: Simulation Methods (Topic) (2010), 10.2139/ssrn.1587822.
  * Bennett (1973) Charles H. Bennett, “Logical reversibility of computation,” IBM Journal of Research and Development 17 (1973), 10.1147/rd.176.0525.

##  Appendix A Automatic Differentiation and Multi-Objective QAE

Another way to compute gradients in the considered setting is automatic
differentiation (AD) Homescu (2011); Pages  _et al._ (2016); Capriotti (2010).
AD repeatedly applies the chain rule to every elementary arithmetic operation
that is used to compute an objective function and keeps track of the
analytical gradient throughout the calculation. Different variants of AD exist
Homescu (2011) and it has been shown that in many practical applications the
gradient can be computed at only a constant overhead, independent of the
dimension Capriotti (2010). In some cases the overhead can even be stated to
be bounded by a factor of four compared to evaluating the function alone at
the expense of larger memory requirements Capriotti and Giles (2010).

Suppose now we want to estimate an expectation value $\mathbb{E}(g(S,x))$ for
a payoff function $g$, a random variable $S$ and some given parameters $x$ as
well as the corresponding gradient $\nabla_{x}\mathbb{E}(g(S,x))$. To
construct the probability oracle $\mathcal{A}(x)$ required by QAE for a fixed
$x$, we usually first create a weighted superposition of all scenarios, then
evaluate the corresponding payoff for each scenario, and last prepare an
objective qubit, i.e., we get

$\displaystyle\sum_{j=0}^{2^{m}-1}\sqrt{p_{j}}\ket{s_{j}}\ket{g(s_{j},x)}\left(\sqrt{1-g(s_{j},x)}\ket{0}+\sqrt{g(s_{j},x)}\ket{1}\right),$
(45)

such that the probability of measuring $\ket{1}$ in the last qubit corresponds
to $\mathbb{E}(g(S,x))$, and where the $s_{j}$ denote the possible
realizations of $S$ represented by $m$ qubits and the $p_{j}$ denote the
corresponding probabilities.

For every scenario $\ket{s_{j}}$, we apply quantum arithmetic to compute the
payoff $\ket{g(s_{j},x)}$. Thus, for each $s_{j}$, we can also use AD in the
same way as classically to compute the gradient $\nabla_{x}g(s_{j},x)$, while
using at most twice the resources required classically due to the need of a
reversible implementation Bennett (1973). Thus, with a constant overhead
compared to the evaluation of the expectation value, this results in the state

$\displaystyle\sum_{j=0}^{2^{m}-1}\sqrt{p_{j}}\ket{s_{j}}\ket{g(s_{j},x)}\bigotimes_{i=1}^{k}\ket{\partial_{i}g(s_{j},x)}.$
(46)

In the following, we show how to use QAE to read out multiple objectives
defined on the same random variables, which then immediately applies to the
gradient as constructed in Eq. (46).

Suppose a random variable $S$ and a set of functions $f_{i}$, $i=1,\ldots,k$,
that map realizations of $S$ to $\mathbb{R}$. Further, suppose we are
interested in estimating the expectation values $\mathbb{E}(f_{i}(S))$ for all
$i$, and that we can construct a state of the form

$\displaystyle\sum_{j=0}^{2^{m}-1}\sqrt{p_{j}}\ket{s_{j}}\bigotimes_{i=1}^{k}\ket{f_{i}(s_{j})}.$
(47)

Then, to estimate the values $\mathbb{E}(f_{i}(S))$, we first introduce $k$
additional $m$-qubit registers $\ket{c_{i}}$, each initialized with some value
$c_{i}$, and second, we use quantum arithmetic to compute the sum

$\displaystyle\sum_{i=1}^{k}c_{i}f_{i}(s_{j})$ (48)

into another register. In other words, we construct an operator that acts as

$\displaystyle\bigotimes_{i=1}^{k}\ket{c_{i}}\ket{0}\bigotimes_{i=1}^{k}\ket{0}\ket{0}\mapsto\bigotimes_{i=1}^{k}\ket{c_{i}}\sum_{j=0}^{2^{m}-1}\sqrt{p_{j}}\ket{s_{j}}\bigotimes_{i=1}^{k}\ket{f_{i}(s_{j})}\ket{\sum_{i=1}^{k}c_{i}f_{i}(s_{j})}.$
(49)

By adding an objective qubit and applying a rotation controlled by the last
register we can also use this to construct a probability oracle
$\mathcal{A}(c)$ that corresponds to the function

$\displaystyle f(c)$ $\displaystyle=$
$\displaystyle\mathbb{E}\left(\sum_{i=1}^{k}c_{i}f_{i}(S)\right).$ (50)

This is a linear function in $c$ and using the quantum gradient algorithm with
respect to $c$ results in

$\displaystyle\nabla_{c}f(c)=\left(\mathbb{E}(f_{1}(S)),\ldots,\mathbb{E}(f_{k}(S))\right)^{T},$
(51)

i.e., in the read out of all $k$ expectation values.

If the values of $f_{i}$ are such that the weighted sum in Eq. (48) satisfies
$\sum_{i=1}^{k}c_{i}f_{i}(s_{j})\leq 1$, since the function by construction is
linear in $c$, the resulting complexity of the quantum gradient algorithm for
a target accuracy $\epsilon>0$ scales as $\mathcal{\tilde{O}}(1/\epsilon)$,
i.e., independent of $k$ (ignoring logarithmic terms), following (Gilyén  _et
al._ , 2019, Thm. 23, arxiv version). The multi-objective QAE requires $k\cdot
m$ additional qubits as well as the weighted sum in Eq. (48), which can be
computed in logarithmic depth by using a divide-and-conquer summation scheme.
If on the other hand, Eq. (48) needs to be normalized by dividing the weighted
sum with a factor $D$, the complexity of the algorithm becomes
$\mathcal{\tilde{O}}(1/\epsilon D)$. When $f_{i}\sim 1,\forall i$, we would
need to choose $D\sim 1/k$, which adds a factor of $k$ back to the complexity
of the algorithm, negating the advantage of this method.

Since Eq. (46) has the required shape, we can immediately apply the multi-
objective QAE to evaluate the gradient that has been evaluated using AD
implemented by quantum arithmetic. If no additional normalization is required,
we can combine AD and (multi-objective) QAE to get the gradient algorithm with
runtime $\mathcal{\tilde{O}}(1/\epsilon)$, i.e., independent of the dimension
and with a quadratic speed-up in the accuracy. Thus, like classically, AD
could represent a promising approach to estimate market risks with a
significant advantage over finite difference schemes. What remains to be
analyzed in more depth, is how this affects the required memory, i.e. qubits,
and how to best automate automatic differentiation for reversible quantum
arithmetic.