Simulating the non-Hermitian dynamics of financial option pricing with quantum computers

Abstract

The Schrödinger equation describes how quantum states evolve according to the Hamiltonian of the system. For physical systems, we have it that the Hamiltonian must be a Hermitian operator to ensure unitary dynamics. For anti-Hermitian Hamiltonians, the Schrödinger equation instead models the evolution of quantum states in imaginary time. This process of imaginary time evolution has been used successfully to calculate the ground state of a quantum system. Although imaginary time evolution is non-unitary, the normalised dynamics of this evolution can be simulated on a quantum computer using the quantum imaginary time evolution (QITE) algorithm. In this paper, we broaden the scope of QITE by removing its restriction to anti-Hermitian Hamiltonians, which allows us to solve any partial differential equation (PDE) that is equivalent to the Schrödinger equation with an arbitrary, non-Hermitian Hamiltonian. An example of such a PDE is the famous Black-Scholes equation that models the price of financial derivatives. We will demonstrate how our generalised QITE methodology offers a feasible approach for real-world applications by using it to price various European option contracts modelled according to the Black-Scholes equation.

A financial derivative is an options contract whose value derives from an underlying financial asset [1]. An options contract defines an agreement between two parties that entails a right to trade an asset at some specified future date for a fixed price. This right to trade agreement thus creates inherent value, which may in turn be traded in the same manner as the underlying financial asset. Consequently, a financial derivative may be viewed as an instrument, which can be used to either exploit arbitrage opportunities or mitigate risk exposure in the market. For this reason, a fundamental task in quantitative finance is how exactly do we determine the fair price of a financial derivative. Determining the fair price of an option is a highly non-trivial task, which is due in part to the stochastic nature of the parameters that define a derivative.

The famous Black–Scholes model [2, 1] is an effective method for determining the fair price of a derivative, and has become the standard for pricing European style financial options. Given the payoff price for an option at the maturity time, we can determine the present price of the option by solving the linear differential equation

$$\frac{\partial u}{\partial t}=-\frac{1}{2}(\sigma x)^{2}\partialderivative[2]{u}{x}-rx\frac{\partial u}{\partial x}+ru,$$ (1)

for $(x,t)\in[x_{0},x_{N}]\times[0,T]$, where the condition $u(x,T)=p(x)$ denotes the payoff of the option. The price of the option is denoted by $u(x,t)$, while $x$ denotes the value of the underlying asset, $t$ represents time, and $T$ is the maturity time. For simplicity, it is assumed that the volatility of the asset, $\sigma$, and the risk-free interest rate, $r$, are constant with respect to time. For convenience, adopting $\tau=T-t$ transforms the Black-Scholes equation Eq. (1) to the initial value problem

$$\frac{\partial u}{\partial\tau}=\frac{1}{2}(\sigma x)^{2}\partialderivative[2]{u}{x}+rx\frac{\partial u}{\partial x}-ru,$$ (2)

for $(x,\tau)\in[x_{0},x_{N}]\times[0,T]$ with the initial condition $u(x,\tau=0)=p(x)$. To numerically solve the Black-Scholes equation, we must discretise the domain $[x_{0},x_{N}]$ to a finite domain and assign appropriate boundary conditions.

The Schrödinger equation models the evolution of the wave function of a quantum mechanical system, and takes the form

$$i\frac{\partial\psi(\vec{x},t)}{\partial t}=\hat{H}\psi(\vec{x},t),$$ (3)

where the Hamiltonian, $\hat{H}$, is a linear differential operator in $\vec{x}$ acting on the wave function $\psi$. Solutions to the Schrödinger equation are expressed in terms of the time evolution operator,

$$\psi(\vec{x},t)=e^{-i\hat{H}t}\psi(\vec{x},0).$$ (4)

The Black-Scholes equation, Eq. (2), can also be expressed in the form of the Schrödinger equation, where its Hamiltonian is given by

$$\hat{H}_{BS}=i\left[\frac{1}{2}(\sigma x)^{2}\partialderivative[2]{x}+rx\frac{\partial}{\partial x}-r\right].$$ (5)

Note that while the Hamiltonian of the Schrödinger equation is a Hermitian operator, which gives rise to unitary time evolution, the Black-Scholes Hamiltonian, Eq. (5), is non-Hermitian, and induces non-unitary time evolution. However, since quantum computers evolve under unitary time evolution, it is the case that simulating non-Hermitian dynamics is not directly feasible on a quantum computer. It is for this reason that quantum computing approaches for solving the Black-Scholes equation have thus far been based on variational algorithms [3, 4, 5] or require post-selection techniques [6].

Another example of non-Hermitian dynamics can be seen in the imaginary time evolution of a quantum system. Following a Wick rotation, which replaces time with an imaginary number $\beta=it$, the Schrödinger equation drives wave functions to become parallel to the ground state of the system. The Wick-rotated form of the Schrödinger equation also takes the form of Eq. (3), but with an anti-Hermitian Hamiltonian. Although the imaginary time Schrödinger equation induces non-unitary dynamics, the normalised evolution can be simulated with quantum algorithms, including quantum imaginary time evolution (QITE) [7] and variational QITE [8].

Variational QITE (varQITE) is a hybrid quantum-classical algorithm that is well suited for noisy intermediate-scale quantum (NISQ) devices. As a variational quantum algorithm, varQITE considers a system of differential equations linking to the gradients of ansatz parameters in imaginary time, and coefficients that depend on measurements of the ansatz. Variational QITE employs a fixed ansatz, where the time complexity is linear in the number of Hamiltonian terms. However, the choice of ansatz is crucial, as it is possible that the states produced by the true imaginary time evolution may not be generated by the particular parameterised ansatz circuit.

On the other hand, the simulated QITE approach is an alternative technique for simulating imaginary time evolution [7]. The technique works by approximating the normalised time evolution operator with Trotter products via unitaries. Simulated QITE with sufficiently large unitary domains is not plagued by barren plateaus, as is the case with its variational counterpart. Simulated QITE on a $k$-local Hamiltonian requires a number of measurements that is exponential in $k$, with the depth of the associated quantum circuits scaling accordingly. Interestingly, however, recent work has focused on optimising the circuit depth and the number of measurements required in simulated QITE. For instance, Fast QITE provides for an exponential reduction in the circuit depth of each unitary and also reduces the number of measurements required per time step, leading to a quadratic speedup over QITE [9]. A time dependent drifted QITE introduces the concept of randomised compiling, which reduces the unitary circuit depth to be a constant and also reduces the number of measurements needed [10]. We also have an implementation of QITE using nonlocal approximation, which reduces circuit depth and is NISQ-friendly [11].

Figure 1: Black-Scholes option pricing simulations using QNUTE. The figure compares the Black-Scholes option prices calculated using QNUTE with varying number of qubits to the corresponding analytical solutions for the following European option types: (a) Call (b) Put (c) Bull Spread (d) Bear Spread (e) Straddle (f) Strangle. The vertical dashed lines at $x=50,75,$ and $100$ correspond to the strike prices of the option contracts. We simulated the solutions for the asset prices $x\in[0,150]$, with the maturity time $T=3$ years, simulated over $N_{T}=500$ time steps. Our simulations used a risk-free interest rate of $r=0.04$, and the volatility $\sigma=0.2$. The unitaries used to approximate the evolution act on all of the qubits used in the simulation.

Although imaginary time evolution was originally envisioned as a technique for determining the ground state of a Hamiltonian [8, 7], the methodology has been recently used as an approach for solving partial differential equations (PDEs), primarily based on varQITE [12, 13, 5, 3, 14]. However, a simulated QITE approach for solving linear PDEs was recently considered [15], although it is restricted to anti-Hermitian Hamiltonians involving only even-ordered derivatives. This application tracks how the non-unitary time evolution scales the quantum state over time. The approach was used to generate solutions to the isotropic heat equation by combining the scale information with the normalised states obtained from QITE.

In this paper, we further widen the scope of simulated QITE by broadening the methodology to simulations involving arbitrary non-Hermitian dynamics. By removing simulated QITE’s underlying restriction to anti-Hermitian Hamiltonians, we enhance the capabilities of the methodology with an ability to simulate arbitrary linear PDEs involving non-unitary time evolution. We have called this generalisation of simulated QITE to arbitrary Hamiltonians quantum non-unitary time evolution (QNUTE).

Results

Quantum Non-Unitary Time Evolution

QNUTE is a quantum algorithm that simulates the dynamics of the Schrödinger equation with an arbitrary non-Hermitian Hamiltonian ${\hat{H}=\sum_{m=1}^{M}i\hat{h}_{m}}$. The non-unitary time evolution operator generated by $\hat{H}$ is approximated by its first order Trotter product, and takes the form

$$e^{-i\hat{H}T}\approx\left(\prod_{m=1}^{M}e^{\hat{h}_{m}\Delta t}\right)^{N_{T}},$$ (6)

where $N_{T}=T/\Delta t$ [16, 17]. The normalised actions of each Trotter step $e^{\hat{h}_{m}\Delta t}$ acting on a state $\ket{\psi}$ are approximated with unitaries of the form $e^{-i\hat{A}\Delta t}$, and implemented with Trotter products of the form

$$e^{-i\hat{A}\Delta t}\approx\prod_{I=1}^{\mathcal{I}}e^{-ia_{I}\hat{\sigma}_{I}\Delta t}.$$ (7)

In Eq. (7), ${\hat{A}=\sum_{I=1}^{\mathcal{I}}a_{I}\hat{\sigma}_{I}}$ is a Hermitian operator with $\hat{\sigma}_{I}$ denoting Hermitian operators chosen such that each unitary $e^{-i\theta\hat{\sigma}_{I}}$ is efficiently implemented with a quantum circuit parameterised by $\theta$. The real-valued coefficients $a_{I}$ are determined by minimising the expression

$$\left\|\frac{e^{\hat{h}_{m}\Delta t}\ket{\psi}}{\sqrt{\expectationvalue{e^{\hat{h}_{m}^{\dagger}\Delta t}e^{\hat{h}_{m}\Delta t}}{\psi}}}-e^{-i\hat{A}\Delta t}\ket{\psi}\right\|,$$ (8)

up to $O(\Delta t)$, which involves solving a system of linear equations, ${(S+S^{\top})\,\vec{a}=\vec{b}}$, constructed using various measurements on $\ket{\psi}$. In particular, we have

$$\begin{aligned} S_{I,J}&=\expectationvalue{\hat{\sigma}_{I}^{\dagger}\hat{\sigma}_{J}}{\psi},\\ c&=\sqrt{1+2\Delta t\,\text{Re}\expectationvalue{\hat{h}_{m}}{\psi}},\\ b_{I}&=\frac{-2}{c}\,\text{Im}\expectationvalue{\hat{\sigma}_{I}^{\dagger}\,\hat{h}_{m}}{\psi},\end{aligned}$$ (9)

see Supplementary Information for further details on the construction. Simulating each Trotter step involves taking $O(\mathcal{I}^{2})$ measurements to construct the $\mathcal{I}\times\mathcal{I}$ matrix equation and generates a quantum circuit of depth $O(\mathcal{I})$. The full simulation therefore requires $O(N_{T}M\mathcal{I}^{2})$ measurements.

The states generated by QNUTE are determined by the choice of $\hat{\sigma}_{I}$. For example, choosing $\hat{\sigma}_{I}$ to encompass all Pauli strings allows us to capture arbitrary state vector rotations in the state space, whereas restricting $\hat{\sigma}_{I}$ to Pauli strings involving an odd number of $\hat{Y}$ gates significantly reduces the operator decomposition count and allows us to capture those rotations that do not introduce complex phases to the quantum state. Given a choice of $\hat{\sigma}_{I}$, the accuracy of the QNUTE implementation is dependent on the support of $\hat{A}$. Ideally, the support of $\hat{A}$ should cover ${D=O(C)}$ adjacent qubits surrounding the support of $\hat{h}_{m}$, where the correlation length $C$ denotes the maximum distance between interacting qubits in the Hamiltonian. However, our choice to express $\hat{A}$ has been in terms of Pauli strings, which gives rise to an exponential dependence on $D$, ${\mathcal{I}=O(2^{D})}$. For this reason, we have considered an inexact implementation of QNUTE that uses a constant domain size $D<C$.

We will demonstrate that QNUTE can be used to approximate solutions to arbitrary linear PDEs with solutions stored in the qubit state vector. Information relevant to the solution is extracted by taking measurements on the final quantum state. It is expected that the number of distinct measurements required to extract the relevant information should scale polynomially with the number of qubits. Further, if it is known that the solution to a PDE will be real-valued and non-negative, then the normalised solution calculated by QNUTE can be extracted obtained by taking the square root of the probability distribution of computational basis states. We will use QNUTE to simulate the Black-Scholes equation, as it has a non-Hermitian Hamiltonian and has non-negative real-valued solutions.

Figure 2: Average fidelities of inexact QNUTE implementations. The figure shows the fidelities of different implementations of inexact QNUTE used to simulate Black-Scholes dynamics averaged over each time step, with the error bars depicting the standard deviation. These simulations share the same parameters values for $r,\sigma,T$ and $N_{T}$ as with the simulations shown in Fig. 1. $n$ denotes the number of qubits used to store the function samples, and $D$ denotes the maximum number of adjacent qubits targeted by the unitaries. The overall low fidelities shown the by inexact QNUTE, where $D<n$, indicate that the Black-Scholes Hamiltonian with linear boundary conditions has a high correlation length, making it difficult to accurately reproduce its evolution with small unitaries.

Simulating Black-Scholes with QNUTE

To model the dynamics of the Black-Scholes equation, we discretise the domain $[x_{0},x_{N}]$ into $2^{n}$ points equally spaced by a distance of $h=\frac{x_{N}-x_{0}}{2^{n}-1}$. The normalised samples of the option price are stored in an $n$-qubit quantum state given by

$$\ket{\bar{u}(\tau)}=\frac{\sum_{k=0}^{2^{n}-1}u(x_{k},\tau)\ket{k}}{\sqrt{\sum_{k=0}^{2^{n}-1}u^{2}(x_{k},\tau)}},$$ (10)

where $x_{k}=x_{0}+kh$. Following Eq. (5), the discretised Black-Scholes Hamiltonian can be represented in terms of a central finite difference matrix of the form

$$-i\hat{H}_{BS}=\begin{bmatrix}\gamma_{0}&\beta_{0}&&\\ \alpha_{1}&\gamma_{1}&\beta_{1}&\\ &\ddots&\ddots&\ddots&&\\ &&\alpha_{2^{n}-2}&\gamma_{2^{n}-2}&\beta_{2^{n}-2}\\ &&&\alpha_{2^{n}-1}&\gamma_{2^{n}-1}\end{bmatrix},$$ (11)

where

$$\begin{aligned} \alpha_{k}&=\frac{\sigma^{2}x_{k}^{2}}{2h^{2}}-\frac{rx_{k}}{2h},\\ \beta_{k}&=\frac{\sigma^{2}x_{k}^{2}}{2h^{2}}+\frac{rx_{k}}{2h},\\ \gamma_{k}&=-r-\alpha_{k}-\beta_{k}.\end{aligned}$$ (12)

Refer to Supplementary Information for the representation of the discretised Hamiltonian of Eq. (11) in the Pauli operator basis.

Norm correction

The scale factor $c$ given in Eq. (9) approximates how the Trotter step scales $\ket{\psi}$ up to $O(\Delta t)$. These approximations can be stored and multiplied to provide an approximation of how the state vector scales over the course of the evolution. Excluding the scenario of the ideal implementation of QNUTE that records a perfect fidelity, errors associated to each scale factor will compound over multiple Trotter steps, which must be corrected periodically. For an anti-Hermitian Hamiltonian $\hat{H}=i\hat{L}$, it was shown that the correction factor can be calculated using knowledge of the non-degenerate ground state $\ket{\psi_{0}}$ of $\hat{L}$ and its corresponding eigenvalue $\lambda_{0}$ [15]. This correction strategy necessarily exploits the mutual orthogonality of the eigenstates of the associated Hamiltonian.

However, since the discretised Black-Scholes Hamiltonian as given in Eq. (11) is not a normal operator, its eigenvectors are not guaranteed to be mutually orthogonal. This, therefore, rules out the norm correction strategy pursued in Ref. [15]. Interestingly, variational QITE has been employed as a technique for solving the Black-Scholes equation. Under this setting, the normalisation factor was considered either as a variational parameter [5] or was determined with prior knowledge of how, specifically, call option prices evolve at the boundary $x_{N}$ [3]. Since the former is not compatible with QNUTE, we generalise the latter approach to cater to various European option types.

Table 1: Average fidelity data for the QNUTE simulations for pricing various European option types. $\mu_{F}$ denotes the mean fidelity over each time step and $\sigma_{F}$ the standard deviation. Simulations involving $D<n$ represent inexact QNUTE implementations.

$n$	$D$	Call		Put
		$\mu_{F}$	$\sigma_{F}$	$\mu_{F}$	$\sigma_{F}$
2	2	1.000	2.67$\times 10^{-5}$	1.000	4.54$\times 10^{-5}$
3	2	0.994	4.77$\times 10^{-3}$	0.982	1.20$\times 10^{-2}$
	4	1.000	3.26$\times 10^{-6}$	1.000	3.06$\times 10^{-6}$
4	2	0.710	2.78$\times 10^{-1}$	0.904	8.02$\times 10^{-2}$
	4	1.000	3.16$\times 10^{-7}$	1.000	2.79$\times 10^{-7}$
5	2	0.158	2.11$\times 10^{-1}$	0.749	6.79$\times 10^{-2}$
	4	0.111	1.79$\times 10^{-1}$	0.655	1.70$\times 10^{-1}$
	6	1.000	4.15$\times 10^{-8}$	1.000	2.10$\times 10^{-7}$
6	2	0.130	7.27$\times 10^{-2}$	0.724	4.18$\times 10^{-2}$
	4	0.122	7.32$\times 10^{-2}$	0.766	5.20$\times 10^{-2}$
	6	1.000	3.81$\times 10^{-8}$	1.000	2.18$\times 10^{-7}$

$n$	$D$	Bull		Bear
		$\mu_{F}$	$\sigma_{F}$	$\mu_{F}$	$\sigma_{F}$
2	2	1.000	2.67$\times 10^{-5}$	1.000	4.54$\times 10^{-5}$
3	2	0.994	4.77$\times 10^{-3}$	0.982	1.20$\times 10^{-2}$
	4	1.000	3.26$\times 10^{-6}$	1.000	3.06$\times 10^{-6}$
4	2	0.710	2.78$\times 10^{-1}$	0.904	8.02$\times 10^{-2}$
	4	1.000	3.16$\times 10^{-7}$	1.000	2.79$\times 10^{-7}$
5	2	0.158	2.11$\times 10^{-1}$	0.749	6.79$\times 10^{-2}$
	4	0.111	1.79$\times 10^{-1}$	0.655	1.70$\times 10^{-1}$
	6	1.000	4.15$\times 10^{-8}$	1.000	2.10$\times 10^{-7}$
6	2	0.130	7.27$\times 10^{-2}$	0.724	4.18$\times 10^{-2}$
	4	0.122	7.32$\times 10^{-2}$	0.766	5.20$\times 10^{-2}$
	6	1.000	3.81$\times 10^{-8}$	1.000	2.18$\times 10^{-7}$

$n$	$D$	Straddle		Strangle
		$\mu_{F}$	$\sigma_{F}$	$\mu_{F}$	$\sigma_{F}$
2	2	1.000	2.67$\times 10^{-5}$	1.000	4.54$\times 10^{-5}$
3	2	0.994	4.77$\times 10^{-3}$	0.982	1.20$\times 10^{-2}$
	4	1.000	3.26$\times 10^{-6}$	1.000	3.06$\times 10^{-6}$
4	2	0.710	2.78$\times 10^{-1}$	0.904	8.02$\times 10^{-2}$
	4	1.000	3.16$\times 10^{-7}$	1.000	2.79$\times 10^{-7}$
5	2	0.158	2.11$\times 10^{-1}$	0.749	6.79$\times 10^{-2}$
	4	0.111	1.79$\times 10^{-1}$	0.655	1.70$\times 10^{-1}$
	6	1.000	4.15$\times 10^{-8}$	1.000	2.10$\times 10^{-7}$
6	2	0.130	7.27$\times 10^{-2}$	0.724	4.18$\times 10^{-2}$
	4	0.122	7.32$\times 10^{-2}$	0.766	5.20$\times 10^{-2}$
	6	1.000	3.81$\times 10^{-8}$	1.000	2.18$\times 10^{-7}$

Consider the Black-Scholes equation, as given in Eq. (5), with option price $u(x,\tau)$ assumed to be linear in $x$ in the neighbourhood of the boundaries $x_{0}$ and $x_{N}$. We will consider linear boundary conditions, since they are widely used in classical option pricing simulations and are known to be numerically stable [18]. Thus, under linear boundary conditions, the option price takes the form ${u(x,\tau)=a(\tau)x+b(\tau)}$ near the boundaries. Substituting this form into Eq. (5) reduces the Black-Scholes equation to an ordinary differential equation (ODE) at the boundaries

$$x\derivative{a}{\tau}+\derivative{b}{\tau}=-rb(\tau).$$ (13)

Solving Eq. (13) yields $a(\tau)=a(0)$ and $b(\tau)=b(0)e^{-r\tau}$, where $a(0)$, and $b(0)$ can be derived from the initial conditions $p(x)$ at each boundary. If $a(0)$ or $b(0)$ are non-zero on at least one of the boundaries, we can rescale a normalised solution to ensure that the value at that boundary is equal to $a(\tau)x+b(\tau)$.

To guarantee that the linear boundary conditions apply during the QNUTE simulation, they must be encoded into the Black-Scholes Hamiltonian. The first and last rows of the matrix in Eq. (11) are updated with the corresponding forward and backward first-order finite difference coefficients, respectively, with the second-derivative terms vanishing as the function is linear. The Black-Scholes Hamiltonian inclusive of linear boundary conditions takes the form

$$-i\hat{H}_{LBS}=\begin{bmatrix}\gamma_{0}^{\prime}&\beta_{0}^{\prime}&&\\ \alpha_{1}&\gamma_{1}&\beta_{1}&\\ &\ddots&\ddots&\ddots&&\\ &&\alpha_{2^{n}-2}&\gamma_{2^{n}-2}&\beta_{2^{n}-2}\\ &&&\alpha_{2^{n}-1}^{\prime}&\gamma_{2^{n}-1}^{\prime}\end{bmatrix},$$ (14)

where

$$\begin{aligned} \gamma_{0}^{\prime}&=\ -r-\frac{rx_{0}}{h},\\ \beta_{0}^{\prime}&=\ \frac{rx_{0}}{h},\\ \alpha_{2^{n}-1}^{\prime}&=\ -\frac{rx_{N}}{h},\\ \gamma_{2^{n}-1}^{\prime}&=\ -r+\frac{rx_{N}}{h}.\end{aligned}$$ (15)

See Supplementary Information for the Pauli decomposition of this Hamiltonian.

Discussion

In this work, we have generalised the quantum imaginary time evolution algorithm to enable the simulation of arbitrary non-Hermitian dynamics on quantum computers. We demonstrated our QNUTE algorithm’s application via a numerical implementation that simulating the pricing dynamics of European options, as dictated by the Black-Scholes equation.

In undertaking these simulations, we assumed that the underlying financial asset had constant volatility, $\sigma$, and risk-free interest rate, $r$. The time dependence of these variables can be encoded in the Hamiltonian with no extra cost to its construction. Further, the inclusion of these variables does not affect the unitary approximations produced by QNUTE, however, modelling volatility and interest rates as stochastic processes may require smaller time steps for more accurate simulations. Indeed, the time dependence of $r$ gives rise to a different boundary ODE, which necessitates modifications to our rescaling protocol.

As depicted in Fig. 1, our implementations of QNUTE were able to match the analytical solutions of the Black-Scholes equation. For convergence, it is important to choose an asset price domain with boundaries such that linear boundary conditions hold for the option’s payoff $u(x,\tau=T)$. A good level of convergence also depends on having access to enough sample points around the strike prices of the option, a lack of which can be seen in the 2-qubit curves in Fig. 1 (c), (d) and (f).

In our implementations of QNUTE, each term, $\hat{h}_{m}$, in the decomposition of the Black-Scholes Hamiltonian was a linear combination of Pauli strings. Since the number of distinct Pauli strings in our decomposition scales exponentially with the number of qubits, an alternative decomposition is required for the scalability of QNUTE for Black-Scholes. Approaches taken to solve PDEs using varQITE have also required the expectation values of finite difference operators. In particular, Liu et al. proposed a scheme to measure such expectation values with a linear overhead [14]. We conjecture that this scheme may be adopted within our approach, leading to exponentially fewer terms in the Hamiltonian decomposition.

The discretised Black-Scholes Hamiltonian has a high degree of correlation between all the qubits used in the simulation. This was demonstrated by implementing inexact QNUTE, wherein the unitary approximations only act on at most $D$ adjacent qubits. Figure 2 depicts the fidelities of inexact QNUTE simulations, averaged over each time step for the various option payoffs, see Table 1 for the exact values. For an increasing number of qubits and fixed domain size $D$, we have it that inexact QNUTE does not capture the correlations between the qubits, rendering it unable to emulate the true time evolution. For future work, we intend to incorporate recent improvements to the simulated QITE methodology, including Fast QITE  [9], time-dependent drifted QITE [10] and QITE with nonlocal approximation [11], within our QNUTE framework to understand their effect on the accuracy and efficiency for simulated PDE dynamics.

Data availability

All data generated and analysed during this study are included in this published article and its supplementary information files.

Code availability

The code that supports the findings of this study is available from the corresponding authors upon reasonable request.

Author Contributions

C.M.W. conceived the project. S.K. performed the numerical simulations. C.M.W. and S.K. analyzed the data and interpreted the results, and proposed improvements. C.M.W. and S.K. both contributed to the writing and editing of the manuscript.

Competing Interests

The authors declare no competing interests.

Supplementary Information

Calculating the unitaries

We approximate the normalised action of a Trotter step $e^{\hat{h}_{m}\Delta t}$ on the state $\ket{\psi}$ by using the unitary $e^{-i\hat{A}\Delta t}$, where $\hat{A}=\sum_{I}a_{I}\hat{\sigma_{I}}$ is a Hermitian operator expressed as a real linear combination of Pauli operator strings indexed by $I$. To calculate the coefficients $a_{I}$, we minimise the distance between the normalised Trotter step and the unitary acting on $\ket{\psi}$ expanded to $O(\Delta t)$. This is achieved by firstly approximating the norm of the Trotter step acting on $\ket{\psi}$. We have it that

$$\begin{aligned} c^{2}&={\expectationvalue{e^{\hat{h}_{m}^{\dagger}\Delta t}e^{\hat{h}_{m}\Delta t}}{\psi}}\\ &=\expectationvalue{(\hat{\mathbbm{1}}+\hat{h}_{m}^{\dagger}\Delta t)(\hat{\mathbbm{1}}+\hat{h}_{m}\Delta t)}{\psi}+O(\Delta t^{2})\\ &=\expectationvalue{\hat{\mathbbm{1}}+(\hat{h}_{m}^{\dagger}+\hat{h}_{m})\Delta t}{\psi}+O(\Delta t^{2})\\ &=1+2\mathrm{Re}\expectationvalue{\hat{h}_{m}}{\psi}\Delta t+O(\Delta t^{2}),\end{aligned}$$ (16)

from which we have

$$c \approx\sqrt{1+2\Delta t\,\mathrm{Re}\expectationvalue{\hat{h}_{m}}{\psi}}.$$ (17)

Next, we set the unitary $e^{-i\hat{A}\Delta t}$ to map $\ket{\psi}$ to the normalised output of the Trotter step

$$e^{-i\hat{A}\Delta t}\ket{\psi} =\frac{e^{\hat{h}_{m}\Delta t}}{c}\ket{\psi},$$ (18)

from which we have

$$(\hat{\mathbbm{1}}-i\hat{A}\Delta t)\ket{\psi} \approx\left(\frac{\hat{\mathbbm{1}}+\hat{h}_{m}\Delta t}{c}\right)\ket{\psi},$$ (19)

and rewriting

$$-i\hat{A}\ket{\psi} \approx\frac{1-c}{c\,\Delta t}\ket{\psi}+\frac{\hat{h}_{m}}{c}\ket{\psi}.$$ (20)

Now, defining

$$\ket{\Delta_{0}}:=\frac{1-c}{c\,\Delta t}\ket{\psi}+\frac{\hat{h}_{m}}{c}\ket{\psi},$$ (21)

and

$$\ket{\Delta}:=-i\hat{A}\ket{\psi},$$ (22)

we solve for $\hat{A}$ such that the squared distance between $\ket{\Delta}$ and $\ket{\Delta_{0}}$ is minimised. It follows that

$$\begin{aligned} \|\ket{\Delta_{0}}-\ket{\Delta}\|^{2}&=(\bra{\Delta_{0}}-\bra{\Delta})(\ket{\Delta_{0}}-\ket{\Delta})\\ &=\braket{\Delta_{0}}-\braket{\Delta_{0}}{\Delta}-\braket{\Delta}{\Delta_{0}}+\braket{\Delta}\\ &=\braket{\Delta_{0}}+\braket{\Delta}-2\mathrm{Re}\braket{\Delta_{0}}{\Delta}.\end{aligned}$$ (23)

Since $\ket{\Delta_{0}}$ does not depend on $\hat{A}$, $\braket{\Delta_{0}}$, in the minimisation, can be ignored. Noting that $\hat{A}$ is Hermitian, we see that

$$\braket{\Delta}=\expectationvalue{\hat{A}^{2}}{\psi},$$ (24)

while the quantity $\braket{\Delta_{0}}{\Delta}$ can be evaluated as

$$\braket{\Delta_{0}}{\Delta} = \left(\frac{1-c}{c\,\Delta t}\bra{\psi}+\bra{\psi}\frac{\hat{h}_{m}^{\dagger}}{c}\right)(-i\hat{A}\ket{\psi})$$ (25)

$$= -i\left(\frac{1-c}{c\,\Delta t}\right)\expectationvalue{\hat{A}}{\psi}-\frac{i}{c}\expectationvalue{\hat{h}_{m}^{\dagger}\hat{A}}{\psi}.$$ (26)

Using that $\hat{A}$ is Hermitian, the expectation value of $\expectationvalue{\hat{A}}{\psi}$ will always be real, which implies that the real part of $\braket{\Delta_{0}}{\Delta}$ only depends on the second term in Eq. (26)

$$\begin{aligned} \mathrm{Re}(\braket{\Delta_{0}}{\Delta})&=\mathrm{Re}\left(\frac{-i}{c}\expectationvalue{\hat{h}_{m}^{\dagger}\hat{A}}{\psi}\right)\\ &=\mathrm{Im}\left(\frac{\expectationvalue{\hat{h}_{m}^{\dagger}\hat{A}}{\psi}}{c}\right)\\ &=-\mathrm{Im}\left(\frac{\expectationvalue{\hat{A}\hat{h}_{m}}{\psi}}{c}\right).\end{aligned}$$ (27)

The function to minimise, therefore, reduces to

$$\begin{aligned} f(\vec{a})&=-2\mathrm{Re}(\braket{\Delta_{0}}{\Delta}+\braket{\Delta})\\ &=\frac{2}{c}\,\mathrm{Im}(\expectationvalue{\hat{A}\hat{h}_{m}}{\psi})+\expectationvalue{\hat{A}^{2}}{\psi}\\ &=\sum_{I}\frac{2a_{I}}{c}\,\mathrm{Im}(\expectationvalue{\hat{\sigma}_{I}\hat{h}_{m}}{\psi})+\sum_{I,J}a_{I}a_{J}\expectationvalue{\hat{\sigma}_{I}\hat{\sigma}_{J}}{\psi}.\end{aligned}$$ (28)

Defining $\vec{b}$ with components

$$b_{I}:=\frac{-2}{c}\,\mathrm{Im}(\expectationvalue{\hat{\sigma}_{I}\hat{h}_{m}}{\psi}),$$ (29)

and the matrix $S$ with entries

$$S_{I,J}:=\expectationvalue{\hat{\sigma}_{I}\hat{\sigma}_{J}}{\psi},$$ (30)

we can rewrite the objective function as

$$\begin{aligned} f(\vec{a})&=\sum_{I,J}a_{I}S_{I,J}a_{J}-\sum_{I}b_{I}a_{I}\\ &=\vec{a}\cdot S\vec{a}-\vec{b}\cdot\vec{a}.\end{aligned}$$ (31)

To minimise $f$ with respect to $\vec{a}$, we equate its gradient to zero, and observe that

$$\begin{aligned} \partialderivative{f}{a_{I}}&=\sum_{J\neq I}(S_{I,J}+S_{J,I})a_{J}+2S_{I,I}a_{I}-b_{I}\\ &=\sum_{I,J}(S_{I,J}+S_{J,I})a_{J}-b_{I}.\end{aligned}$$ (32)

It follows that $\nabla f=(S+S^{\top})\vec{a}-\vec{b}$, and, hence, to minimise $f$, we solve for $\vec{a}$ in

$$(S+S^{\top})\vec{a}=\vec{b}.$$ (33)

The Black-Scholes Hamiltonian in the Pauli basis

In order to discretise the Black-Scholes equation with the form, as given in Eq. (2),

$$\partialderivative{u}{\tau}=\frac{\sigma^{2}}{2}x^{2}\partialderivative[2]{u}{x}+rx\partialderivative{u}{x}-ru,$$ (34)

we rewrite the derivative operators in terms of finite difference representations and substitute the $x$ and $x^{2}$ terms with the diagonal matrices

$${\hat{X}=\begin{bmatrix}x_{0}&&\\ &\ddots&\\ &&x_{N}\end{bmatrix}}\ \textrm{and}\ {\hat{X}^{2}=\begin{bmatrix}x_{0}^{2}&&\\ &\ddots&\\ &&x_{N}^{2}\end{bmatrix}},$$ (35)

respectively, where $x_{k}=x_{0}+kh$. We are required to represent these matrices as linear combinations of the Pauli operators in order to be able to simulate their dynamics with QNUTE. Firstly, let us define the following one-qubit matrices.

$$\begin{aligned} \hat{A}_{\nwarrow}^{(1)}:=\begin{bmatrix}1&0\\ 0&0\end{bmatrix}&=\frac{\hat{I}+\hat{Z}}{2}\\ \hat{A}_{\searrow}^{(1)}:=\begin{bmatrix}0&0\\ 0&1\end{bmatrix}&=\frac{\hat{I}-\hat{Z}}{2}\\ \hat{A}_{\nearrow}^{(1)}:=\begin{bmatrix}0&1\\ 0&0\end{bmatrix}&=\frac{\hat{X}+i\hat{Y}}{2}\\ \hat{A}_{\swarrow}^{(1)}:=\begin{bmatrix}0&0\\ 1&0\end{bmatrix}&=\frac{\hat{X}-i\hat{Y}}{2}\end{aligned}$$ (36)

Further, consider the $n$-qubit generalisations of Eq. (36) given by

$$\begin{aligned} \hat{A}_{\nwarrow}^{(n)}&:=\hat{A}_{\nwarrow}^{(1)}\otimes\hat{A}_{\nwarrow}^{(n-1)},\\ \hat{A}_{\searrow}^{(n)}&:=\hat{A}_{\searrow}^{(1)}\otimes\hat{A}_{\searrow}^{(n-1)},\\ \hat{A}_{\nearrow}^{(n)}&:=\hat{A}_{\nearrow}^{(1)}\otimes\hat{A}_{\nearrow}^{(n-1)},\\ \hat{A}_{\swarrow}^{(n)}&:=\hat{A}_{\swarrow}^{(1)}\otimes\hat{A}_{\swarrow}^{(n-1)}.\end{aligned}$$ (37)

Tensor product strings comprising $n$ copies of the four 1-qubit operators in Eq. (36) allows us to generate matrices that have a 1 in exactly one entry of a $2^{n}\times 2^{n}$ matrix with 0s elsewhere. We will use these $n$-qubit tensor strings of the operators in Eq. (36) to define the linear boundary conditions for the Black-Scholes Hamiltonian.

Let us construct the $\hat{X}$ matrix for $n$-qubits. We have

$$\hat{X}^{(n)}=\begin{bmatrix}x_{0}&&&\\ &x_{1}&&\\ &&\ddots&\\ &&&x_{N}\end{bmatrix}=x_{0}\hat{\mathbbm{1}}+h\begin{bmatrix}0&&&\\ &1&&\\ &&\ddots&\\ &&&2^{n}-1\end{bmatrix}=x_{0}\hat{\mathbbm{1}}+h\hat{\chi}^{(n)}.$$ (38)

Consider the matrix $\hat{\chi}^{(n)}$. For $n=1$, $\hat{\chi}^{(1)}=\hat{A}_{\searrow}^{(1)}$, while for $n=2$, we have

$$\hat{\chi}^{(2)}=\begin{bmatrix}\hat{\chi}^{(1)}&0\\ 0&2\hat{I}+\hat{\chi}^{(1)}\end{bmatrix}=\hat{A}_{\nwarrow}^{(1)}\otimes\hat{\chi}^{(1)}+\hat{A}_{\searrow}^{(1)}\otimes(2\hat{I}+\hat{\chi}^{(1)}).$$ (39)

Generalising to $n$-qubits, we have it that

$$\hat{\chi}^{(n)}=\begin{bmatrix}\hat{\chi}^{(n-1)}&0\\ 0&2^{n-1}\hat{I}^{\otimes n-1}+\hat{\chi}^{(n-1)}\end{bmatrix}=\hat{I}\otimes\hat{\chi}^{(n-1)}+2^{n-1}\left(\hat{A}_{\searrow}^{(1)}\otimes\hat{I}^{\otimes n-1}\right).$$ (40)

Squaring both sides yields

$$(\hat{\chi}^{(n)})^{2} =\left(\hat{I}\otimes\hat{\chi}^{(n-1)}+2^{n-1}\left(\hat{A}_{\searrow}^{(1)}\otimes\hat{I}^{\otimes n-1}\right)\right)\left(\hat{I}\otimes\hat{\chi}^{(n-1)}+2^{n-1}\left(\hat{A}_{\searrow}^{(1)}\otimes\hat{I}^{\otimes n-1}\right)\right)$$

$$=\hat{I}\otimes(\hat{\chi}^{(n-1)})^{2}+2^{n}\left(\hat{A}^{(1)}_{\searrow}\otimes\hat{\chi}^{(n-1)}\right)+2^{2(n-1)}\left((\hat{A}^{(1)}_{\searrow})^{2}\otimes\hat{I}^{\otimes n-1}\right)$$

$$=\hat{I}\otimes(\hat{\chi}^{(n-1)})^{2}+\hat{A}^{(1)}_{\searrow}\otimes\left(2^{n}\hat{\chi}^{(n-1)}+2^{2(n-1)}\hat{I}^{\otimes n-1}\right).$$ (41)

Next, let us consider the first derivative $\partialderivative{x}$. This operator is approximated by the central difference matrix $\frac{1}{2h}\hat{D}_{1}^{(n)}$, where

$$\hat{D}_{1}^{(n)}=\begin{bmatrix}0&1&&&\\ -1&0&1&&\\ &\ddots&\ddots&\ddots&\\ &&-1&0&1\\ &&&-1&0\end{bmatrix}.$$ (42)

For one qubit case, $n=1$, we have

$$\hat{D}_{1}^{(1)}=\begin{bmatrix}0&1\\ -1&0\end{bmatrix}=i\hat{Y}.$$ (43)

Generalising to $n$-qubits, we have

$$\hat{D}_{1}^{(n)}=\begin{bmatrix}\hat{D}_{1}^{(n-1)}&\hat{A}_{\swarrow}^{(n-1)}\\ -\hat{A}_{\nearrow}^{(n-1)}&\hat{D}_{1}^{(n-1)}\end{bmatrix}=\hat{I}\otimes\hat{D}_{1}^{(n-1)}+\hat{A}_{\nearrow}^{(1)}\otimes\hat{A}_{\swarrow}^{(n-1)}-\hat{A}_{\swarrow}^{(1)}\otimes\hat{A}_{\nearrow}^{(n-1)}.$$ (44)

Next, let us consider the second derivative $\partialderivative[2]{x}$, which is approximated by the central difference matrix $\frac{1}{h^{2}}\hat{D}_{2}^{(n)}$, where

$$\hat{D}_{2}^{(n)}=\begin{bmatrix}-2&1&\\ 1&-2&1&\\ &\ddots&\ddots&\ddots&\\ &&1&-2&1&\\ &&&1&-2\end{bmatrix}.$$ (45)

For $n=1$, we have

$$\hat{D}_{2}^{(1)}=\begin{bmatrix}-2&1\\ 1&-2\end{bmatrix}=-2\hat{I}+\hat{X}.$$ (46)

Generalising to $n$-qubits, it can be shown that

$$\hat{D}_{2}^{(n)}=\begin{bmatrix}\hat{D}_{2}^{(n-1)}&\hat{A}_{\swarrow}^{(n-1)}\\ \hat{A}_{\nearrow}^{(n-1)}&\hat{D}_{2}^{(n-1)}\end{bmatrix}=\hat{I}\otimes\hat{D}_{2}^{(n-1)}+\hat{A}_{\nearrow}^{(1)}\otimes\hat{A}_{\swarrow}^{(n-1)}+\hat{A}_{\swarrow}^{(1)}\otimes\hat{A}_{\nearrow}^{(n-1)}.$$ (47)

Using the $\hat{X}^{(n)}$, $\hat{D}_{2}^{(n)}$ and $\hat{D}_{2}^{(n)}$ representations above, the $n$-qubit Black-Scholes Hamiltonian takes the form

$$-i\hat{H}_{BS}^{(n)}=\frac{\sigma^{2}}{2h^{2}}(\hat{X}^{(n)})^{2}\hat{D}_{2}^{(n)}+\frac{r}{2h}\hat{X}^{(n)}\hat{D}_{1}^{(n)}-r\hat{I}^{\otimes n}.$$ (48)

Expressing Eq. (48) in matrix form over $n$-qubits, we have it that

$$-i\hat{H}_{BS}^{(n)}=\begin{bmatrix}\gamma_{0}&\beta_{0}&\\ \alpha_{1}&\gamma_{1}&\beta_{1}&\\ &\ddots&\ddots&\ddots&\\ &&\alpha_{2^{n}-2}&\gamma_{2^{n}-2}&\beta_{2^{n}-2}\\ &&&\alpha_{2^{n}-1}&\gamma_{2^{n}-1}\end{bmatrix},$$ (49)

where

$$\alpha_{k}=\frac{\sigma^{2}x_{k}^{2}}{2h^{2}}-\frac{rx_{k}}{2h},\quad\beta_{k}=\frac{\sigma^{2}x_{k}^{2}}{2h^{2}}+\frac{rx_{k}}{2h},\quad\text{and}\quad\gamma_{k}=-r-\alpha_{k}-\beta_{k}.$$ (50)

To include linear boundary conditions, we neglect the second derivative terms and substitute the first and last row of Eq. (48) with forward and backward finite difference coefficients, respectively. The revised first row coefficients that take into account the linear boundary conditions are

$$\gamma_{0}^{\prime}=-r-\frac{rx_{0}}{h}\quad\text{and}\quad\beta_{0}^{\prime}=\frac{rx_{0}}{h},$$ (51)

while the revised last row coefficients that into account the linear boundary conditions are

$$\alpha_{2^{n}-1}^{\prime}=-\frac{rx_{N}}{h}\quad\text{and}\quad\gamma_{2^{n}-1}^{\prime}=-r+\frac{rx_{N}}{h}.$$ (52)

Consequently, for $n\geq 2$ qubits, the Black-Scholes Hamiltonian with linear boundary conditions is written as

$$-i\hat{H}_{LBS}^{(n)} = -i\hat{H}_{BS}^{(n)}+(\gamma_{0}^{\prime}-\gamma_{0})\hat{A}_{\nwarrow}^{(n)}+(\beta_{0}^{\prime}-\beta_{0})(\hat{A}_{\nwarrow}^{(n-1)}\otimes\hat{A}_{\nearrow}^{(1)})$$ (53)

$$+\ (\alpha_{2^{n}-1}^{\prime}-\alpha_{2^{n}-1})(\hat{A}_{\searrow}^{(n-1)}\otimes\hat{A}_{\swarrow}^{(1)})+(\gamma_{2^{n}-1}^{\prime}-\gamma_{2^{n}-1})\hat{A}_{\searrow}^{(n)}.$$

The rescaling protocol

Under the assumption of linear boundary conditions, $u(x,\tau)$ is linear with respect to $x$ in the neighbourhood of $x_{0}$ and $x_{N}$, hence, we have

$$u(x,\tau)=a(\tau)x+b(\tau)$$ (54)

for $x\approx x_{0}$ and $x\approx x_{N}$. Substituting Eq. (54) into the Black-Scholes equation, Eq. (34), we obtain the boundary equation

$$\partialderivative{\tau}\left(a(\tau)x+b(\tau)\right) =rx\partialderivative{x}\left(a(\tau)x+b(\tau)\right)-r\left(a(\tau)x+b(\tau)\right),$$ (55)

which implies

$$x\derivative{a}{\tau}+\derivative{b}{\tau}=-rb(\tau).$$ (56)

Equation (56) is an ODE with solutions

$$a(\tau)=a(0)\quad\text{and}\quad b(\tau)=b(0)e^{-r\tau}.$$ (57)

Determining $a(0)$ and $b(0)$ at each boundary, Eq. (54) becomes a closed form expression that explains how $u(x,\tau)$ evolves at the boundaries. Assuming $h$ is sufficiently small, we can approximate $a(0)$ and $b(0)$ at each boundary by assuming that the closest sample point also follows the linear condition. Thus, on the left boundary, using $u(x_{0},0)=a_{0}(0)x_{0}+b_{0}(0)$ and $u(x_{1},0)=a_{0}(0)x_{1}+b_{0}(0)$, we have

$$a_{0}=\frac{u(x_{1},0)-u(x_{0},0)}{h}\quad\text{and}\quad b_{0}(0)=u(x_{0},0)-a_{0}x_{0}.$$ (58)

Similarly, on the right boundary, we have

$$a_{N}=\frac{u(x_{N},0)-u(x_{N-1},0)}{h}\quad\text{and}\quad b_{N}(0)=u(x_{N},0)-a_{N}x_{N}.$$ (59)

Using Eq. (58) and Eq. (59), we define a protocol that rescales the normalised state $\ket{\bar{\psi}_{u}(\tau)}$ obtained from QNUTE by a factor $C_{*}(\tau)$ such that

$$\bra{0}C_{*}(\tau)\ket{\bar{\psi}_{u}(\tau)}=a_{0}x_{0}+b_{0}(0)e^{-r\tau},$$ (60)

or

$$\bra{2^{n}-1}C_{*}(\tau)\ket{\bar{\psi}_{u}(\tau)}=a_{N}x_{N}+b_{N}(0)e^{-r\tau}.$$ (61)

The choice between Eq. (60) and Eq. (61) depends on the preference assigned to the boundaries during simulation. This work used the left boundary when modelling put options, and the right boundary when modelling call options. The protocol fails if $a(0)=b(0)=0$ on both boundaries.