Variational quantum simulations of stochastic diﬀerential equations

,


I. INTRODUCTION
Stochastic differential equations (SDEs), which describe the time evolution of random variables, are among the most important mathematical tools for modeling uncertain systems in diverse fields, such as finance [1], physics [2], and biology [3].From the expectation values of the simulated random variables, we can often extract information about the system of interest.Since the expectation values rarely admit analytical solutions, they are usually obtained by numerical methods such as the Monte Carlo method [4].However, those numerical methods incur high computational costs, especially in high-dimensional problems such as the SDEs of financial applications [5][6][7].Therefore, a method that can speed-up SDE simulations is urgently demanded.
Such a speed up can be achieved on quantum computers.Throughout the past decade, technological developments have realized a primitive form of quantum computers called noisy intermediate-scale quantum (NISQ) devices [8], which can handle problems that are intractably large for classical computers [9].NISQ devices can operate only a few tens to hundreds of qubits without error correction, so they cannot run quantum algorithms requiring deep and complicated quantum circuits.Although quantum algorithms are expected to outperform classical ones on specific computing tasks [10][11][12][13], they usually exceed the capability of NISQ devices.Accordingly, NISQ devices have been leveraged with heuristic algorithms that solve real-world problems.For example, in quantum chemistry and condensed matter physics, the variational quantum eigensolver (VQE) algorithm [14,15] can calculate the ground-state energies of given Hamiltonians [16,17].Another example is quantum machine learning with variational quantum circuits [18][19][20][21].Both algorithms variationally optimize the tuneable classical parameters in quantum circuits, so the speedups of the computation over classical computers and the precision of the obtained results are not guaranteed in general.
Several quantum-computing-based methods obtain the expectation value of a function that takes an SDE solution as its argument.However, all of these methods require prerequisite knowledge of the SDE solution.In [22], the partial differential equation describing the time evolution of the expectation value was simulated by a variational quantum computation, which requires prederivation of the partial differential equation of the expectation value.In [23] and [24], the probability distribution of the SDE solution was embedded in the quantum state, and the expectation value was calculated by a quantum amplitude estimation algorithm (QAE).In this case, the probability distribution of the SDE solutions must be known in advance.As the solution to the SDE is not found, the partial differential equation of the expectation value must also be derived, or the SDE solved beforehand.
In this study, to solve an SDE with quantum algo-rithms, we apply a tree model approximation [25], and hence obtain a linear differential equation describing the probability distribution of SDE solutions.This differential equation is then solved by a variational quantum simulation (VQS) [26][27][28][29][30].Note that linear differential equations can be solved by a quantum linear solver algorithm (QLSA) [12,31,32], which is expected to be quantumaccelerated.However, the QLSA requires a large number of ancilla qubits and deep circuits and is possibly executable only on quantum computers with error correction.Our proposed method possesses several desirable features.First, the probability distribution is simulated by the tree-model approximation, so the model requires only the SDE.No prior knowledge of the probability distribution or expectation value is required.Therefore, our method is applicable to more general SDEs than previous methods.Second, once the VQS is performed, the variational parameters are obtained as classical information, and the probability distribution of the simulation results can be used to compute various expectation values.We can also compute path-dependent expectation values because the time series of the probability distribution is obtained.Third, the algorithm is less resource-intensive than the QLSA.Since VQS is a variational algorithm, it is difficult to estimate the exact computational cost, but VQS requires only a few ancilla qubits and calculates the expectation value for relatively shallow unitary gates at each time step.The number of qubits and the depth of the circuit are expected to be much smaller than QLSA.As our method uses a new scheme for embedding probability distributions in quantum states, the method for computing expectation values is also new.We additionally found that the expectation values are more simply determined by our method than by the QAE.The proposed method facilitates the application of SDEs in quantum computing simulations and is expected to impact various scientific fields.
The remainder of this paper is organized as follows.Section II reviews the trinomial tree-model approximation and the VQS, before introducing our method.Our main theoretical results are contained in Secs.III, IV.Section III proposes a VQS-based method that simulates the dynamics of the probability distribution of the stochastic process in the trinomial tree model.The quantum circuits and operators that perform the VQS are also constructed in this section.Section IV calculates the expectation value of the random variable using the state obtained by simulating SDE with the VQS.Section V discusses the advantages of our method and compares them with previous studies.Section VI numerically evaluates our algorithm on two SDE prototypes: the geometric Brownian motion and the Ornstein-Uhlenbeck process.Conclusions are presented in Section VII.Appendix A analyses the complexity of calculating the expectation value, and Appendix B generalizes our result to a multiple-variable process.Appendix C evaluates the error of expectation values from piecewise polynomial approximation.
FIG. 1. Lattice of the trinomial tree model.Nodes (circles) at (t, x) represent the events in which X(t) takes the value x.Edges represent the transition probabilities between the nodes.The stochastic process starts at node (t0, x0) and "hops" to the other nodes depending on the transition probabilities.

II. PRELIMINARIES
This section reviews the main ingredients of this paper: the trinomial tree-model approximation of the SDE [25] and the VQS algorithm [26][27][28][29].In Sec.III, we combine both ingredients into a method that simulates the SDE by the VQS.

A. Trinomial tree-model approximation of the stochastic differential equation
Let us consider a random variable X(t) taking values on an interval I ⊂ R. We refer to I as an event space.The SDE of a single process {X(t)} t∈[0,T ] , which is a time-series of random variables from t = 0 to t = T , is defined as [1] where µ(X(t), t), σ(X(t), t) are real valued functions of time t and the variable X(t), and W denotes the Brownian motion.In the main text, our proposal is applied to a single process (extensions to multi-variables cases are described in Appendix B).
The tree model numerically simulates the time evolution of an SDE.Let us consider an SDE simulation of the process with event space [0, x max ] from t = 0 to t = T .We discretize the time as t i ≡ i∆t (i = 0, 1, . . ., N t ) and the event space as x i ≡ i∆x (i = 0, 1, . . ., N x ), where N t ∆t = T and N x ∆x = x max .In this discretization scheme, we define a (N x + 1) × (N t + 1) lattice on which each node (i, j) is associated with a probability Prob[X(t j ) = x i ] and each edge represents a transition between two nodes, as shown in Fig. 1.Here, we adopt the trinomial tree model, which has three transition probabilities as follows: These probabilities were chosen to reproduce the first and second moment (mean and variance, respectively) of the random variable X(t) in Eq. (1).Following the Euler-Maruyama method [7], the SDE is discretized as where z ∼ N (0, 1) and O(∆t 2 ) terms are ignored.The conditional expectation value and variance are respectively expressed as The corresponding moments on the trinomial tree model are Equating these moments and considering the normalization condition p u (x, t) + p m (x, t) + p d (x, t) = 1, we obtain In summary, the trinomial tree-model approximates the original SDE by discretizing it on the lattice and setting the transition probabilities between the nodes to reproduce the first and the second moments of the process.The trinomial tree model simulates the SDE as follows.First, the closest value to x ini in {x i } i∈[0,Nx] is set to x i0 , and the probabilities are set as Prob[X(t 0 ) = x i0 ] = 1, Prob[X(t 0 ) = x i =i0 ] = 0. Next, the probability distribution of X(t 1 = ∆t) is calculated using the transition probabilities given by Eqs.(3)(4) (5).Repeating this step for X(t j )(j = 2, 3, ..., N t − 1) yields all probabilities Prob[X(t j ) = x i ] at node (i, j), from which any properties related to the process X(t), such as the expectation values of X(T ) under some function f , E[f (X(T ))], can be determined.In option-pricing financial problems, the nodes of the tree model denote the prices of the option, and the problems are sometimes to be solved in the backward direction from time t.In such cases, the boundary condition is set at t = T .

B. Variational quantum simulation (VQS)
This subsection introduces the VQS algorithm [26][27][28][29], a quantum-classical hybrid algorithm that simulates both unitary and non-unitary time evolution with possibly shallow quantum circuits.Therefore, the VQS algorithms is especially suitable for NISQ devices.
We are interested in the non-unitary time evolution of an unnormalized quantum state | ψ(t) on an n-qubit system, defined as where L(t) is a time-dependent (possibly non-Hermitian) linear operator.To simulate the dynamics of | ψ(t) , let us introduce the following ansatz quantum state |ṽ(θ(t)) : where α(t) is a real number, θ(t) ≡ (α(t), θ 1 (t)) ≡ (α(t), θ 1 (t), . . ., θ M (t)) are variational parameters of the ansatz, |0 is some reference state, and R(θ The gates depend on their parameters and on other non-parametric gates.In particular, G k is assumed as a multi-qubit Pauli gate {I, X, Y, Z} ⊗n . The VQS algorithm maps the dynamics of the quantum state, Eq. ( 6), to those of the variational parameters θ(t) of the ansatz.The mapping is performed by McLachlan's variational principle [33] min where |ϕ ≡ ϕ|ϕ .This equation reduces to an Euler-Lagrange equation, We define θ 0 (t) ≡ α(t) for notational simplicity.When simulating the dynamic Eq. ( 6), one starts from the initial parameters θ ini corresponding to the initial state | ψ(t = 0) = |ṽ(θ ini ) .The time derivative θ(t = 0) is calculated by Eq. ( 9) with |ṽ(θ ini ) in Eqs.(10) and (11).After a small time step δt, the parameters are obtained as θ(δt) = θ ini + δt • θ(t = 0).Repeating this procedure obtains the dynamics of θ(t) and the state |ṽ(θ(t)) .The terms M k,j and V k can be evaluated by the quantum circuits depicted in Fig. 2 [29].The normalized |0 is actually prepared on quantum computers and is multiplied by the normalization constant α(t) in post-processes of the result of the circuit measurements.Decomposing the operator where U k is an easilyimplementable unitary operator (e.g., a multi-qubit Pauli operator) and λ k is a complex coefficient, we must evaluate O(M 2 ) + O(M k term (t)) distinct quantum circuits.The circuits need an ancilla qubit other than the qubits of the system of interest, along with control operations of G k and U k .Therefore, to ensure a feasible VQS algorithm, both M, k term (t), and the depth of the unitaries U k must be O(poly(n)).

III. SOLVING STOCHASTIC DIFFERENTIAL EQUATIONS BY VARIATIONAL QUANTUM SIMULATION
This section presents one of our main results.The SDE simulated by the above-described trinomial tree model is reformulated as the non-unitary dynamics of a quantum state | ψ(t) embedding the probability distribution of the random variable X(t).We explicitly state for the L(t) operator of the VQS and decompose it by the polynomial number of the sum of easily-implementable unitaries.

A. Embedding the probability distribution into a quantum state
To simulate the trinomial tree model of the target SDE by VQS, we define an unnormalized quantum state containing the discretized probability distribution of the random variable X(t j ): where {|i } Nx i=0 is the computational basis.We call this state a directly embedded state.For simplicity, we assume that N x = 2 n − 1, where n is the number of qubits.
Note that this embedding of the probability distribution into the quantum state differs from most of the lit-erature, in which (aiming for a quantum advantage) the expectation values of a probability distribution are calculated using QAE [34].In the literature, the probability distribution is expressed as a normalized quantum state The expectation value of the distribution, for some function f , is computed by the QAE.In this embedding method, VQS cannot be used because the differential equation describing the time evolution of the probability distribution is nonlinear.There are ways to solve the nonlinear differential equation with a quantum algorithm [35][36][37][38], but they require more complicated quantum circuits.
Because our embedding (12) differs from this embedding scheme, we also developed a method for evaluating its expectation values (see Sec. IV).Note that the normalization constant α in Eq. ( 7) may be exponentially small.In fact, for a uniform distribution Prob[X(t j ) =

B. Reformulating the trinomial tree model and applying the variational quantum simulation
In the trinomial tree model, the probability Prob Substituting the transition probabilities (3), ( 4) and ( 5) into this expression and denoting P (x, t) ≡ Prob[X(t) = x], we get In the limit ∆t → 0, one obtains where P (t) ≡ (P (x 0 , t), P (x 1 , t), . . ., P (x As shown in Eq. ( 16), the time evolution of the state | ψ(t) , or where corresponds to the time evolution of the probability distribution {Prob[X(t) = x i ]} 2 n −1 i=0 .Equation ( 18) is the essence of our proposal to simulate VQS-based SDE simulation: specifically, the VQS algorithm applied to Eq. ( 18) obtains the time-evolved probability distribution as the quantum state | ψ(t) .Hereafter, when the distinction is clear in context, we denote the operator L(t) by L(t) as in Eq. ( 16).

C. Construction of L(t)
As explained in the previous section, in the VQS, we evaluate Eqs.(10) and (11), and decomposes L(t) into a sum of easily-implementable unitaries (composed of single-qubit, two-qubit, and few-qubit gates).These evaluations are important for a feasible VQS.This subsection discusses the explicit decomposition of L(t) given by Eq. (19).
To express the operator L(t) in Eq. ( 19), we define operators These operators can be constructed from the n-qubit cyclic increment/decrement operator CycInc(n) ≡ where |−1 , |2 n are identified with |2 n − 1 , |0 , respectively.These gates are implemented as a product of O(n) Toffoli, CNOT, and X gates with O(n) ancilla qubits [39].
which can be implemented [13] as a product of O(n 2 ) Toffoli, CNOT, and single qubit gates.Using meaning that V ± (n) can be decomposed into a sum of two unitaries composed of O(n 2 ) few-qubit gates.Finally, we define the operator D(n) by where Z i is a Z gate acting on the ith qubit.Therefore, D(n) is a sum of O(n) unitaries composed of a singlequbit gate.It follows that Let us recall Expanding σ 2 (x i , t) and µ(x i , t) as we can decompose L(t) as follows: Therefore, the L(t) decomposition realizes a feasible VQE (Eq.( 18)).

IV. CALCULATION OF EXPECTATION VALUES
In the previous section, we propose a method to simulate the SDE by calculating the dynamics of the probability distribution of a random variable X(t) using VQS.However, in many cases, the goal of the SDE simulation is not the probability distribution of X(t), but the expectation value E[f (X(t))] of X(t) for some function f .In this section, we introduce a means of calculating this expectation value.

A. Problem Setting
Given a function f (x) : R → R, we try to calculate the expectation value E[f (X(T ))] of the SDE (1) at time t = T .The expectation value can be explicitly written as Here, we assume that The additional error from this piecewise polynomial approximation is evaluated in Appendix C. As x is finite, the range of f is also finite.Thus, by shifting the function f by a constant, we can ensure that the range of f is positive and that the expectation value is also positive, i.e.E[f (X(T ))] ≥ 0. In most situations (such as pricing of European call options as we see in Sec.IV C) the number of intervals d does not scale with the number of qubits n.
The imaginary part of the expectation value ℑ ψ(t)| U | ψ(t) is evaluated by the circuit with an S † gate inserted to the left of the second H gate.

B. General formula for calculating expectation values
We now on compute the expectation value (28) using the quantum state | ψ(t) (Eq.( 12)).First, we consider a non-unitary operator satisfying and decompose S f into a sum of easily-implementable unitaries as As |0 0| = I − C n Z • X ⊗n is also a sum of easilyimplementable unitaries as explained in the previous subsection, the Hermitian observable which is again a sum of unitaries.With this decomposition, the left-hand side of Eq. ( 30) is computed by evaluating ψ(t Because we set E[f (X(T ))] ≥ 0, the left hand side of Eq. (30) will determine the expectation value.
There are two options to evaluate the quantities ψ(t . The first one is to use the Hadamard test depicted in Fig. 3.The second one is to use quantum phase estimation [40,41].The former one requires shallower quantum circuits but is inefficient in terms of the number of measurements to determine the quantities with fixed precision.The detailed computational complexity of these methods is given in Sec.V and Appendix A. Next, we explain the construction of the operator S f in Eq. ( 29) and its decomposition.We first define an operator where χ [0,a] (x) is the indicator function valued as 1 for x ∈ [0, a] and 0 for otherwise.Using the binary expansion of a/∆x, we can obtain the decomposition of S χ[0,a] hence the decomposition of S f .As a ∈ [0, x max ], there exists k a ∈ N such that ∆x2 ka−1 ≤ a < ∆x2 ka , 0 < k a ≤ n.The binary expansion of a/∆x is given by a/∆x = ka−1 j=0 s j 2 j , s j ∈ {0, 1}.We define the list of l as l 1 , l 2 , . . ., l B (= k a − 1) satisfying s l = 1 in ascending order, and also define an interval for l ∈ {l 1 , l 2 , . . ., l B }. Using χ a l , we devide [0, a/∆x] into disjoint intervals as follows: The indicator operator S χ [0,a] is obtained by summing the indicator operators on each interval.In binary expansion, the k a th and the lth bit of i ∈ χ a l are 1, and the bit below l is either 0 or 1. Accordingly, X should act on the bit taking 1, and H should act on the bit taking either of {0, 1}.The indicator operator S χ a l on χ a l is defined as follows: where In addition, we define We can construct S χ [0,a] by summing Eqs. ( 35) and (37) on each interval.S χα k on interval α k ≡ [a k , a k+1 ] is which is a sum of at most O(n) unitaries composed of O(n) gates.Using S χα k , we obtain and S f is constructed as In summary, evaluation of the expectation value is calculated by the following steps.3. Decompose S f |0 0| S † f into a sum of unitary terms and calculate each term using the circuits in Fig. 3.

Divide the domain of the target func
When the target function f on each interval is written by a low-degree polynomial (i.e., L is small), especially by a linear function (as in the pricing of European call options shown below), our algorithm can efficiently calculate the expectation value because the number of unitaries O(d 2 n 2L+2 ) gets not so large.When the function f is approximated by the polynomial, we can estimate the error of the expectation value stemming from that approximation.If we want to suppress the error below ǫ, the number of unitaries becomes O(x 2 max ǫ − 2 L+1 n 2L+2 ) (the derivation is presented in Appendix C).Note that as L is increased, ǫ − 2 L+1 becomes smaller while n 2L+2 becomes larger.The number of unitaries, therefore, is not monotonic with respect to L, and there may be an optimal L for the desired accuracy.We note that evaluation of expectation values of those unitaries can be performed completely in parallel by independent quantum devices.

C. Pricing of The European Call Option
As a concrete example, we present the pricing of a European call option with the Black-Scholes (BS) model, which is one of the simplest financial derivatives.The holder of a European call option is entitled to buy the asset at a predetermined strike price at maturity.The price of a European call option with strike price K ≥ 0, interest rate r ≥ 0, and maturity T ≥ 0 is defined by the conditional probability Here, E Q denotes the expectation value under the riskneutral probability measure.Stochastic processes are assumed to follow geometric Brownian motion in the BS model, but are described by more complex mechanisms in other models.Even in these models, the expression Eq. ( 41) of the price of the European call option is the same with the present case.
Setting the probability distribution of X T conditioned by For simplicity, we assume ∆x = 1 and In this case, there are only two intervals [0, K − 1] and [K, 2 n−1 ], and the polynomial in each interval is of firstorder degree at most.Therefore, we can calculate the price of the European call option by Eq. ( 30).

V. POSSIBLE ADVANTAGES OF OUR METHOD
In this section, we discuss the advantages of our method compared to previous studies, as well as the possible quantum advantages.
In general, the SDEs addressed in this paper can be transformed into a partial differential equation (PDE) of the function e f (x, t), where e f (x, t) gives the expectation value E[f (X(T − t))|X(0) = x], by Feynman-Kac formula [1].In fact, the authors of [22] performed a variational quantum computation of a PDE of this function.We point out two advantages of our method compared with this strategy using Feynman-Kac formula.First, the resulting PDE by Feynman-Kac formula must be solved backwardly in time from t = T to t = 0, with the initial condition at t = T being related to the functional form of f (X).It is not trivial to prepare the initial state |ψ(T ) corresponding to the initial condition; the authors of [22] executed an additional VQE to prepare the initial state.Second, when using the Feynman-Kac formula, the initial condition of the PDE is different for each function f for which we want to calculate the expectation value E[f (X(T ))].If we want to calculate a different expectation value E[f ′ (X(T ))], we need to run the whole algorithm simulating the PDE with the different initial state corresponding to f ′ .On the other hand, in our method, once we perform VQS, we obtain the probability distribution of X(T ) as a quantum state and the corresponding variational parameters to reproduce it.We only need to redo the part of the expectation value calculation (Sec.IV) for different f ′ .
The authors of [23] embedded the probability distribution by quantum arithmetic.Their embedding, proposed in [42], requires O(2 n ) gates to embed the probability distribution into an n-qubit quantum state.To moderate the gate complexity, the authors of [24] embedded the probability distribution using a quantum generative adversarial network, which requires only O(Poly(n)) gates.The probability distribution function can also be approximated by a lth-order piecewise polynomial, which can be embedded with O(ln 2 ) gates even in quantum arithmetic [43].However, both methods require prior knowledge of the probability distribution to be embedded.In contrast, our method does not require prior knowledge of the embedding probability distribution since our method simulates the time evolution of a given SDE.
We now compare the computational cost to calculate expectation values with previous studies.In [23] and [24], by employing QAE, the expectation value (Eq.( 28)) was calculated by using an oracle that is complex quantum gate reflecting the functional form of f for O(1/ǫ) times, where ǫ is the precision for the expectation values.The classical Monte Carlo method requires O(1/ǫ 2 ) sampling for precision ǫ, so their methods provide a second-order acceleration.On the other hand, our method measures the expectation value of each term of Eq. ( 31) using the Hadamard test (Fig. 3) or the quantum phase estimation (QPE) [40,41].As shown in Appendix A, the total number of measurements to obtain the expectation value with precision ǫ is O(1/γǫ 2 ) for the Hadamard test and O(log(1/γǫ)) for QPE, where γ is some factor.We note that the depth of the circuit is O(1/γǫ) in QPE, which is in the same order as the QAE whereas our method requires not an complicated oracle but a relatively-small unitary.Hence, when the factor γ is not too small, our method combined with QPE can also exhibit quantum advantage for the evaluation of the expectation values.The factor γ depends on the parameters of the polynomial approximation (a (m) k , d, L), the domain of the approximated function x max , and the probability distribution {Prob[X = x i ]} 2 n −1 i=0 .The detailed evaluation of γ is described in Appendix A.

VI. NUMERICAL RESULTS
In this section, our algorithm is applied to two stochastic processes, namely, geometric Brownian motion and an Ornstein-Uhlenbeck process, which are commonly assumed in financial engineering problems.Geometric Brownian motion simply models the fluctuations of asset prices, and the Ornstein-Uhlenbeck process is a popular model of interest rates.
The ansatz circuit is identical for both models and shown in Fig. 4. As the amplitudes of the quantum state must be real, the ansatz contains only CNOT and RY gates.This depth-k circuit repeats the entangle blocks composed of CNOTs and RY gates k times.The parameters of geometric Brownian motion were r = 0.1, σ = 0.2, ∆x = 1, and t ∈ [0, 4] and those of the Ornstein-Uhlenbeck process were r = 7, σ = 0.5, η = 0.01, ∆x = 1, and t ∈ [0, 4].We simulate the quantum circuits without noise using numpy [44] and jax [45].We set the number of qubits n = 4 and the number of repetitions of entangle blocks k = 2, 3.

B. Results
Panels (a) and (b) of Fig. 5 present the numerical simulations of geometric Brownian motion and the Ornstein-Uhlenbeck process, respectively.In comparison, we also provide a probability density function (PDF) for the solution of the SDE equation obtained by solving the Fokker-Planck equation [46] analytically.We can see that our method describes the time evolution of the probability distribution well.
We calculated the means (Fig. 5(c),(d)) and variances (Fig. 5(e),(f)) of the resulting distributions.We also present the mean and variance obtained from the analytical solution and the solution of ( 16) using the Runge-Kutta method.Because of the approximation with the tree model, even the results of the Runge-Kutta method slightly differ from the analytical solution.In the case of VQS with k = 2, we see that the error from the analytical solution is larger than that of the k = 3 case.This is because the number of VQS parameters is less than the number of lattice points in the event space when k = 2, i.e., the degrees of freedom of ansatz are less than the degrees of freedom of the system, and thus the errors due to the ansatz appear.In the case of k = 3, the number of parameters in ansatz is sufficient, and thus the results are closer to the results of the Runge-Kutta method.

VII. CONCLUSION
This paper proposed a quantum-classical hybrid algorithm that simulates SDEs based on VQS.A continuous stochastic process was discretized in a trinomial tree model and was reformulated as a linear differential equation.The obtained differential equation was solved with VQS, obtaining quantum states representing the probability distribution of the stochastic processes.As our method can embed the probability distribution of the solution of a given SDE into the quantum state, it is applicable to general SDEs.We note that our methods can apply to the Fokker-Plank equation, which also gives the time-evolution of the probability distributions of SDE solutions.
Because the embedding methods of the probability distribution differ in the proposed method and the conventional quantum algorithm, we proposed another method for computing the expectation value.We approximated the functions to calculate expectation values by piecewise polynomials and constructed operators corresponding to the polynomial in each interval.The operators were constructed as sums of unitary operators, which are composed of easily-implementable gates.The expectation value was then computed using the sum of unitary operators.Our algorithm was validated in classical simulations of geometric Brownian motion and the Ornstein-Uhlenbeck process.Both processes were well simulated by the algorithm.Our algorithm is expected to efficiently simulate other stochastic processes, provided that L(t) can be written as a polynomial linear combination of unitary matrices.
Let us summarize the computational cost of our method presented in this work.Our method consists of two parts; one is to perform VQS to simulate the SDE, and the other is to calculate the expectation value of the SDE solution.In the part of running VQS, we decompose matrix L(t) in Eq. ( 16) into a sum of O(n mmax ) different unitaries composed of O(n 2 ) few-qubit gates, where m max is the largest order of the polynomial expansion of µ, σ in Eq. (27).At each time step of VQS, the vector V k in Eq. ( 11) is evaluated as a sum of O(n mmax ) measurement results of the circuits depicted in Fig. 2. As m max is typically finite and small (∼ 1, 2) in most practical applications, the computational cost (i.e., the number of gates in quantum circuits, the number of different circuits to run) of the simulation of SDE is O(Poly(n)).In contrast, QLSA [12,31,32] requires much deeper and more complex quantum circuits and a large number of ancilla qubits because it uses the Hamiltonian simulation and the quantum Fourier transform.This is an advantage of our method leveraging the variational quantum algorithm.
In the part of the expectation value evaluation of the SDE solution, we evaluate it by running different O(d 2 n 2L+2 ) quantum circuits, where d and L are the number of intervals and the order of the piecewise polynomial approximation of the function f in Eq. ( 28), respectively.Each circuit is constructed to compute an expectation value ψ|U |ψ of a unitary U that contains O(n 4 ) quantum gates.When we adopt the Hadamard test (Fig. 3) as such a quantum circuit, the number of measurements to suppress statistical error of the expectation value below ǫ is O(1/γǫ 2 ), where γ is a factor defined in Appendix A. This O(1/ǫ 2 ) scaling is the same as the classical Monte Carlo method to compute the expectation values from the probability distribution of the SDE solution.When we choose the QPE-type circuit to evaluate ψ|U |ψ , the number of measurements becomes O(log(1/γǫ)) while the depth of the circuit in terms of U is O(1/γǫ).This situation can provide a quantum advantage for computing the expectation value of the SDE solution.The error from the piecewise polynomial approximation of f can be made small by increasing d or L, which is detailed in Appendix C.This study focused on computational finance because financial engineering is among the most popular applications of stochastic processes.Pricing of derivatives, and many other problems in financial engineering, satisfy the conditions of the proposed method.However, as the stochastic processes themselves are quite general, the proposed method is expected to contribute to solving problems in various fields.
for d = 1, . . ., D and the covariance of the variables satisfies Cov[X k (t j+1 ) − X k (t j ), X l (t j+1 ) − X l (t j ) |X k (t j ) = x, X l (t j ) = y] = p (k,l) uu − p (B8) As is the same for the case of a single variable we set the transition amplitudes by equating Eqs.(B3),(B4),(B5) with (B6),(B7),(B8).If the solutions of p (k) u,d , p (k,l) uu,ud,du,dd are proportional to ∆t, the linear differential equation can be derived by taking the limit of ∆t → 0 (as in the one-dimensional case Eq. ( 16)).
When D > 1, one should note the numbers of variables and conditional expressions.As the numbers of p m , p (k) u,d , p (k,l) uu,ud,du,dd are 1, 2D, 2D(D − 1), respectively, the number of independent variables is 2D 2 under the normalized probability conditions.On the other hand, the number of equations of the mean, variance, and covariance are D, D, D(D − 1)/2, respectively, so the total number of equations is D(D + 3)/2.When D > 1, the number of variables exceeds the number of conditions, so an infinite number of transition probabilities satisfy the condition.
Here, we show there is indeed a solution of the transition amplitudes which admit taking limit ∆t → 0 and obtain the linear differential equitation of the probability distributions of the SDE.Fixing p (B12) Here, we omit the arguments of µ d and σ d to simplify the notation.

Mapping to VQS and construction of L(t)
In the multivariate case, we can construct L(t) as described in Sec.III.For notational simplicity, we denote Here, we expand σ k (x (k) , t), µ k (x (k) , t) as We also define the operators These operators satisfy the following equations: m and (D(n)) m are composed of the sum of O(n m ) unitaries, each composed of O(n 2 ) few-qubit gates.In typical SDEs, the orders m σ , m µ can be set to small values.For example, geometric Brownian motion case, m = 1 (see Sec. VI).

FIG. 3 .
FIG. 3. Quantum circuit for evaluating the real part of an expectation valueℜ ψ(t)| U | ψ(t) of a unitary operator U = QiQ † i ′ , QiC n Z • X ⊗n Q † i ′ .The imaginary part of the expectation value ℑ ψ(t)| U | ψ(t) is evaluated by the circuit with an S † gate inserted to the left of the second H gate.
B2) where {z d } D d=1 is sampled from the multi-variable Gaussian distribution, E[z d ] = 0, Var[z d ] = 1, Corr[z k , z k ] = ρ kl .The first and second first and second conditional moments satisfy ∆x (k) ∆x (l) .

2 d
the number of variables becomes D(D + 3)/2, which is slightly asymmetric (because we consider only p k uu to be nonzero), but agrees with the number of conditional expressions.In this case, the transition probabilities arep (k,l) uu = σ k σ l ρ kl ∆x (k) ∆x (l) ∆t, k σ d ρ kd ∆x (k) ∆x (d) , k σ d ρ kd ∆x (k) ∆x (d)   − k =l σ k σ l ρ kl ∆x (k) ∆x (l)∆x (d) 2 ∆t.