Variational quantum evolution equation solver.

Variational quantum algorithms offer a promising new paradigm for solving partial differential equations on near-term quantum computers. Here, we propose a variational quantum algorithm for solving a general evolution equation through implicit time-stepping of the Laplacian operator. The use of encoded source states informed by preceding solution vectors results in faster convergence compared to random re-initialization. Through statevector simulations of the heat equation, we demonstrate how the time complexity of our algorithm scales with the Ansatz volume for gradient estimation and how the time-to-solution scales with the diffusion parameter. Our proposed algorithm extends economically to higher-order time-stepping schemes, such as the Crank-Nicolson method. We present a semi-implicit scheme for solving systems of evolution equations with non-linear terms, such as the reaction-diffusion and the incompressible Navier-Stokes equations, and demonstrate its validity by proof-of-concept results.

Introduction

Partial differential equations (PDEs) are fundamental to solving important problems in disciplines ranging from heat and mass transfer, fluid dynamics and electromagnetics to quantitative finance and human behavior. Finding new methods to solve PDEs more efficiently—including making use of new algorithms or new types of hardware—has been an active area of research.

Recently, the advent of quantum computers and the invention of new quantum algorithms have provided a novel paradigm for solving PDEs. A cornerstone of many of these quantum algorithms is the seminal Harrow-Hassidim-Lloyd (HHL) algorithm[1] for solving linear systems, which can be utilized to solve PDEs by discretizing the PDE and mapping it to a system of linear equations. Compared to classical algorithms, the HHL algorithm can be shown to exhibit an exponential speedup. Unfortunately, attractive as it may sound, the HHL algorithm works only in an idealized setting, and a list of caveats must be addressed before it can be used to realize a quantum advantage[2]. Moreover, implementing HHL and many other quantum algorithms would require the use of a fault-tolerant quantum computer, which may not be available in the near future[3]. Instead, the machines we have today are imperfect, noisy intermediate-scale quantum (NISQ) devices[4] with both coherent and incoherent errors limiting practical circuit depths.

Over the last few years, variational quantum algorithms (VQAs) have emerged as a leading strategy to realize a quantum advantage on NISQ devices. Specifically, VQAs employ shallow circuit depths to optimize a cost function, expressed in terms of an Ansatz with tunable parameters, through iterative evaluations of expectation values[5]. Applications of VQAs include the variational quantum eigensolver (VQE) for finding the ground or excited states of a system Hamiltonian[6-8], the quantum approximate optimization algorithm (QAOA) for solving combinatorial optimization problems[9], and solvers for linear[10-12] and non-linear[13] systems of equations.

Here, we are interested in variational quantum algorithms for solving differential equations[14], such as the Black–Scholes equation[15,16], the Poisson equation[17,18], and the Helmholtz equation[19]. Specifically, the Poisson equation can be solved efficiently through explicit decomposition of the coefficient matrix derived from finite difference discretization[17] using minimal cost function evaluations[18] and shallower circuit depth compared to other non-variational quantum algorithms[14,20-22]. A natural question to ask, then, is whether such variational algorithms for Poisson equations can be extended to solving evolution equations, i.e. partial differential equations including a time domain. McArdle et al.[23] proposed a variational quantum algorithm which simulates the real (imaginary) time evolution of parametrized trial states via forward Euler time-stepping of the Wick rotated Schrödinger equation, thereby solving the Black–Scholes equation, and by extension, the heat equation[15,16]. Besides issues of Ansatz selection and quantum complexity, time-stepping based on an explicit Euler method may be unstable, a limiting condition exacerbated by noise. With existing variational quantum algorithms[10], an implicit scheme for evolution equations is expected to preserve any quantum advantage[1] over classical algorithms, with reduced time complexity[18].

This paper is organized as follows. In section “Theory”, we outline general implicit time-stepping schemes for solving evolution partial differential equations and propose the use of a variational quantum solver to resolve the Laplacian operator iteratively. In section “Applications to the heat/diffusion equation”, we apply the variational quantum algorithm to solving a heat or diffusion equation without source terms as a proof of concept. With that, we explore potential applications to more general evolution problems with non-linear source terms, including the reaction–diffusion (section “Applications to the reaction–diffusion equations”) and the Navier–Stokes equations (section “Applications to the Navier–Stokes equations”), where variables can be coupled through semi-implicit schemes.

Theory

Consider the second-order homogeneous evolution equation defined on the set , where denotes a d-dimensional bounded spatial domain and , where denotes a bounded temporal domain, aswhere is a function of spatial vector and time t, is the diffusion coefficient and f is an unspecified source term. For now, Dirichlet and Neumann boundary conditions are applicable on the boundary , respectively,where is the outward normal derivative on boundary .

For a two-dimensional rectangular domain , partitioning the space-time domain yields the space-time grid pointswhere , and are prescribed positive integers, such that , , , , and , where and . The discrete domain grid is denoted by and boundary grid by .

The finite difference (FD) approximation for the second-order spatial derivative (5-point) of the Laplacian operator taken at iswhere and are diffusion parameters.

Using first-order FD for temporal derivative weighted by , the evolution Eq. (1) can be expressed in vector shorthand aswhere is the identity matrix of the same size, and .

Depending on the choice of parameter , actual time-stepping may follow an explicit (forward Euler) method (), an implicit (backward Euler) method (), a semi-implicit (Crank–Nicolson) method () or a variable- method[24]. The explicit method () is efficient for each time-step but is only stable if it satisfies the stability condition . The implicit (backward Euler) method () is unconditionally stable and first-order accurate in time (), which readsThe semi-implicit Crank–Nicolson (CN) method () is popular as it is not only stable, but also second-order accurate in both space and time (), which readswhere . However, the CN method may introduce spurious oscillations to the numerical solution for non-smooth data unless the algorithm parameters satisfy the maximum principle[25].

Variational quantum solver

Here, we explore a variational quantum approach towards the solution of the evolution equation (1). In addition to potential quantum speedup, a variational quantum algorithm could also benefit from data compression, where a matrix of dimension N can be expressed by a quantum system with only qubits, where N is the size of the problem. Consider the Poisson equation, which is a time-independent form of Eq. (1), expressed asThe Laplacian operator in one dimension can be discretized using the finite difference method in the x direction into an N N coefficient matrix aswhere refers to either the Dirichlet (D) or Neumann (N) boundary condition, and and . This extends naturally to higher dimensions, for instance in the y direction.

A variational quantum solution is to prepare a state such that is proportional to a state in a way that satisfies Eq. (10). To do that, a canonical approach[10-12] is to first decompose the matrix A over the Pauli basis (where , , and are the Pauli matrices and is the identity matrix) aswhere are the coefficients of A in the Pauli basis. Using simple operators , , the number of terms in the decomposition can be reduced to [17]. A more efficient approach, however, is to express A as a linear combination of unitary transformations of simple Hamiltonians[18]. Accordingly, the decomposition of A in one dimension can be written as[19]where and , as before except here, and . Here, S is the n-qubit cyclic shift operator defined as

The expectation values of a Hamiltonian H including the shift operator S are evaluated by applying the unitary shift operator to the quantum state[18],where is an arbitrary n-qubit state and . Note that Eq. (13) can be re-written asSince expectation values of the identity operator are equal to 1, i.e. , evaluating the expectation value of the operator requires only the evaluation of expectation values of the simple Hamiltonians ( for Dirichlet boundary condition). The required number of quantum circuits is therefore limited to a constant [18]. Similar decomposition expressions apply to problems of higher dimensions, including in the y direction[19].

Once the matrix A is decomposed, a parameterized quantum state is prepared using an Ansatz represented by a sequence of quantum gates parameterized by applied to a basis state , such that . Here, we use a hardware-efficient Ansatz consisting of multiple layers of gates across n qubits entangled by controlled-X gates (see Fig. 1). For the source term f in (10), a quantum state is prepared by encoding a real vector with the unitary , such that . Depending on the actual input, conventional amplitude encoding methods[26,27] may introduce a global phase that must be corrected by its complex argument for computing in the real space.

Figure 1

Schematic of the quantum circuit hardware-efficient Ansatz used in this study.

With and , the cost function E can be optimized in terms of A as[18]where . The norm of the state vector is represented by , where

The quantum circuit required for the numerators of (17) and (18) consists of an encoding unitary and a parameterized Ansatz (Fig. 1), both oppositely controlled by an ancilla qubit placed in superposition[17]. As for the denominators, the number of quantum circuits required corresponds to the number of decomposed terms of the Hamiltonian (13), each paired with the Ansatz . Finally, the resulting states of these circuits are measured in the computational basis.

Using classical optimization tools, the cost function (17) is minimized with updated iteratively until convergence is reached. The optimization process follows either a gradient-based or gradient-free approach, depending on how the gradient of the cost function is evaluated. A gradient-free optimizer is guided by an estimate of the inverse Hessian matrix, whereas a gradient-based optimizer by the partial derivative of the cost function E with respect to parameters , i.e. , which can be evaluated by a quantum computer (for details, see[18,19]). Regardless of the choice of gradient optimizer used, the optimization routine halts when the cost error falls under a convergence threshold () whence the parameters are at optimum . The converged solution vector satisfies[10]where is the norm of the solution to the Poisson equation (10).

In this study, we propose to solve the evolution equation (1) through successive time-stepping of the quasi-steady Poisson equation using a variational quantum algorithm. Using a parameter set obtained at time-step k, we encode a normalized source state from and seek an implicit solution towhere is the parameter set and is the norm at next time-step . This process is then iterated in time up to number of time-steps as desired (see Algorithm 1).

The cost function may be evaluated on a quantum computer by computing each of the inner products in the expression in Eq. (17) separately. Using the decomposition provided in Eq. (16), these inner products may be expressed in terms of expectation values for preparable states and simple Hermitian operators . Each expectation value is evaluated on a quantum computer by preparing the state using the quantum circuits described above and then measuring the operator in the state [18].

In this study, the variational quantum algorithm is implemented in Pennylane (Xanadu)[28] using a statevector simulator with the Qulacs[29] plugin as a backend for quantum simulations, and the L-BFGS-B optimizer for parametric updates. Amplitude encoding is carried out via the standard Mortonnen state preparation template[30] with custom global phase correction. For hardware emulation via the QASM simulator (Qiskit), we refer the reader to the excellent cost-sampling analysis of Sato et al. [18].

Applications to the heat/diffusion equation

Consider the following one-dimensional heat or diffusion equation without a source term

Dirichlet conditions are applied on the boundaries of a 1D domain , where and , such that the boundary vector is known for all t.

To solve Eq. (22), the variational quantum evolution algorithm (Algorithm 1) can be employed with a suitable time-stepping scheme (7). For the implicit Euler (IE) method (8), the matrix A and source state can be decomposed intowhere .

For the Crank-Nicolson (CN) method (9), it follows that the term, which carries a small but non-trivial evaluation cost, can be eliminated using the source state of the previous time-step , leading towhere . Here, the presence of a term in is not unexpected due to temporal finite differencing at second-order accuracy.

For a space-time domain , let the number of time-steps be and the spatial intervals be , where n is the number of qubits, and is the diffusion parameter. We employ the Dirichlet boundary condition with boundary values and initial values . With initial random parameters , we run a limited-memory Broyden-Fletcher-Goldfarb-Shanno boxed (L-BFGS-B) optimizer[31-34] to optimize with absolute and gradient tolerances set at .

Figure 2a compares solutions obtained from the variational quantum solver (23) and classical methods to a 1D heat or diffusion problem in time-increments of 0.1, where the number of qubits and Ansatz layers expressed as a set n- l, are 3- 3 and 4- 4. Here we define the time-averaged trace error aswhere is the normalized classical solution vector at time k. The trace errors of solutions shown in Fig. 2a are 0.0008 and 0.0025 for n- l sets of 3-3 and 4-4 respectively.

Figure 2

(a) Implicit variational quantum solutions to a 1D heat conduction or diffusion problem in time-increments of 0.1, with boundary values , initial values and diffusion parameter . (Left) Qubit–layer count n--3 and time-averaged trace error ; (Right) n--4, . Both classical and quantum solutions overlap with vanishingly small trace errors. (b) Cost function vs. number of optimization steps for 10 runs. Inset: Input parameters are re-initialized randomly, , before each time-step for 5 runs.

Figure 2b shows how the cost function E depends on the number of optimization steps for n-l of 3-3 and 4-4 (10 sampled runs each). Each distinct step in E represents sequential optimization from solution at time-step k towards the solution at . For small time-step , provides a good initial parameter set for solving optimization step . If the Ansatz parameters were re-initialized randomly before each time-step, significantly more optimization steps would be required on average for convergence for each run (see Fig. 2b, inset).

Time complexity

Here we briefly examine the time complexity of the quantum algorithm excluding the classical computing components. Following the analysis of the variational Poisson solver[18], the time complexity of the proposed variational evolution equation solver per time-step readswhere the terms within the inner parentheses indicate the time complexity of the state preparation scaling as , which consists of the Ansatz depth l, the encoding depth [35], the depth of the circuit needed to implement the n-qubit cyclic shift operator , and that of the number of shots required for estimation of expectation values up to a mean squared error of . The required number of quantum circuits depends on the boundary conditions applied (3 for periodic, 4 for Dirichlet and 5 for Neumann conditions), scaling only as . is the time-averaged number of function evaluations,where is the sum of function evaluations required for a run with time-steps. Using a gradient-based optimizer, the time complexity for gradient estimation via quantum computing would scale as the Ansatz volume representing the number of quantum circuits required for parameter shifting. Otherwise, with a gradient-free optimizer, the time complexity simply contributes towards as additional function evaluations required to evaluate the Hessian for gradient descent.

To see if the time complexity for gradient-free optimization scales as , we plot the time-averaged number of function evaluations against the number of parameters nl (Fig. 3a). Indeed, we found that scales reasonably with nl (see trendline of slope 1), despite apparent tapering at higher l. Figure 3b shows that the time-averaged trace error decreases with circuit depth l, even for over-parameterized quantum circuits where the number of layers exceeds the minimum required for convergence, [36]. For low grid resolution , the trace error is limited to a minimum of Since the time complexity for solving the Poisson equation classically is , where , quantum advantage could be realized with the proposed algorithm with linear time scaling by at sub-exponential time complexity[18].

Figure 3

Logarithmic plots of time-averaged (a) number of function evaluations vs. number of parameters nl, (b) trace error vs. number of layers l, (c) number of iterations and (d) trace error vs. diffusion parameter for up to . Each data point and error bar represents, respectively, the mean and standard deviation out of 25 runs.

For deep and wide quantum circuits, the increase in optimization time is exacerbated by the presence of barren plateaus, or vanishingly small gradients in the energy landscape, where re-initialization can leave one trapped at a position far removed from the minimum[37-39]. Conversely, short time-steps lead to efficient solution trajectories that remain close to the local cost minima, leading to significant reduction in optimization times. To verify this, we conduct numerical simulations varying the diffusion parameter with up to time . Figure 3c shows that the number of iterations, or required optimization steps, per time-step increases linearly with . Close inspection of the time-averaged trace distance shows bimodal distributions at higher , which separates success and failure during convergence towards the global minimum (see Fig. 3d, dotted lines), resembling local minima traps due to poor optimization or expressivity of Ansätze[19,40].

Discretization error

Time evolution can be at a higher order, specifically for the Crank-Nicolson method. The problem statement is identical to the previous one, except with Dirichlet boundary values and the initial condition ), where we use as the spatial length of the domain. This admits an exact analytical solution,Figure 4 compares variational quantum and exact solutions using implicit Euler and Crank-Nicolson (CN) schemes. The discretization error for the higher-order CN scheme is reduced significantly, especially at lower grid resolution (). Although the complexity costs for both methods (23) and (24) are similar, note however that the CN method may introduce spurious oscillations for non-smooth data[24], an issue which may be exacerbated by quantum noise[41].

Figure 4

Variational quantum solutions to a 1D heat conduction or diffusion problem for qubit-layers (a) n--3 and (b) n--4 in time-increments of 0.1 using implicit Euler and Crank–Nicolson schemes, with Dirichlet boundary values , initial values and diffusion parameter . Dashed lines denote exact solutions.

Higher dimensions

The preceding analysis can be extended to higher dimensions. Consider the following two-dimensional heat or diffusion equation in , where :Under the implicit Euler scheme (8), the matrix A and source state can be decomposed into

Dirichlet conditions are applied on the boundaries, where and . Let the number of spatial grid intervals be and , where is the number of qubits allocated to the x grid, to the y grid, and is the total number of qubits. Accordingly, A is decomposed in terms of simple Hamiltonians in x and y aswhere are simple Hamiltonians to be evaluated ( for Dirichlet boundary condition).

Figure 5 shows solution snapshots to a 2D heat conduction or diffusion problem taken at time with Dirichlet boundary values and , initial values , and diffusion parameters . Results obtained from variational quantum solver agree with classical solutions with time-averaged trace errors of up to .

Figure 5

(a) Contour solution plot of 2D heat conduction or diffusion problem at time on a x-y square grid (qubit-layer ---3-6) with Dirichlet boundary values and , initial values , and diffusion parameters . (b) Solution vectors in time-increments of 0.1.

Applications to the reaction–diffusion equations

Here, we extend applications of our variational quantum solver to evolution equations with non-trivial source terms. Consider a two-component homogeneous reaction-diffusion system of equationswhere is a concentration tensor, is a diffusion tensor and is a coupled reaction source term. First proposed by Turing[42], the reaction-diffusion equations are useful for understanding pattern formation and self-organization in biological and chemical systems[43,44], such as morphogenesis[45] and autocatalysis[46].

Here, we propose a semi-implicit time-stepping scheme, whereby the coupled, non-linear source term is solved at the current time-step k. With explicit source term , the implicit Euler scheme (8) readsThe two-component tensor and source state can then be decomposed intowhere is the two-component diffusion parameter vector. With a linear Hermitian source matrix f, a fully implicit time-stepping scheme becomes available (Appendix A).

Implementation

The semi-implicit variational quantum solver solves for the Laplacian for each component using a quantum computer and the solution vectors are explicitly coupled through source terms prior to re-encoding in preparation for the next time-step (see Algorithm 2).

1D Gray–Scott model

The Gray-Scott model[46] was originally conceived to model chemical reactions of the type , , where U, V and P are chemical species with reaction termwhere and are kinetic rate constants.

An interesting class of Gray-Scott solutions involves periodic splitting of chemical wave pulses[24,47]. Here, we conduct a pulse splitting numerical experiment under limited spatial and temporal resolutions, using input parameters , and , with a mid-pulse initial condition and Dirichlet boundary conditions aswhere time t extends up to on .

Figure 6a shows how an initial mid-pulse can spontaneously and periodically split in space and time, a phenomenon captured using variational quantum diffusion reaction solver (see Algorithm 2) even on relatively low spatial resolutions.

Figure 6

Space-time solutions of (a) mid-pulse wave-splitting in 1D two-component Gray-Scott model obtained using semi-implicit variational quantum reaction-diffusion solver on grid points up to , for chemical species (left) and (right). Parameters include , , and . (b) Traveling wave solutions for 1D Brusselator model on grid points up to on Neumann boundary conditions. Parameters include , , and .

1D Brusselator model

So far, we have been looking at only Dirichlet boundary conditions. Here, we demonstrate a test example for Neumann boundary conditions in a diffusion-reaction model, namely, the Brusselator model[48], which was developed by the Brussels school of Prigogine to model the behavior of non-linear oscillators in a chemical reaction system. The model reaction term readsUsing , and , with initial conditionswhere time t extends up to on .

Figure 6b shows how a chemical pulse can be spontaneously created, which continually travels leftwards in time, creating traveling waves that appear as striped patterns in time despite low spatial resolutions.

Applications to the Navier–Stokes equations

The Navier-Stokes equations are a set of non-linear partial differential equations that describes the motion of fluids across continuum length scales. There are several studies aimed at applying quantum algorithms to computational fluid dynamics (see review[49]), ranging from reduction of partial differential equations to ordinary differential equations[50] and quantum solutions of sub-steps of the classical algorithm[51,52] to the quantum Lattice Boltzmann scheme[53].

Here, we look into the potential use of variational quantum algorithms to evolve the fluid momentum equations in time. Consider the incompressible Navier-Stokes equations in non-dimensional formwhere is the velocity vector and p is the fluid pressure. The ratio is the Reynolds number, where is the characteristic flow velocity across a characteristic length-scale and is the fluid kinematic viscosity.

Unlike other temporal evolution equations, the incompressible Navier-Stokes equations cannot be time-marched directly as the resultant velocities do not satisfy the continuity constraint (45), and hence are not divergence-free. To resolve this, the projection method[54], also known as the predictor-corrector or fractional step method, separates the solution time-step into velocity and pressure sub-steps, also known as the predictor and corrector steps.

Projection method

Predictor step

The predictor step first approximates an intermediate velocity by solving the fluid momentum equation (44) in the absence of pressure, i.e. the Burgers’ equations[55], of the formThrough a semi-implicit scheme, the viscous terms are handled implicitly using the variational quantum evolution equation solver and the non-linear inertial terms explicitly as source terms using classical computation. For quantum algorithms for non-linear problems, the reader is referred to separate works on quantum ordinary differential equation solvers[50,55], Carleman linearization[56] and a variational quantum nonlinear processing unit (QNPU)[13].

On a two-dimensional domain with Dirichlet boundary conditions, the tensor and source state can be decomposed aswhere and . is an operator which approximates the non-linear inertial term, where are diagonal matrices with velocity vectors along the diagonals and B is a divergence matrix discretized through center differencing, for instance in the x direction, aswhere refers to either Dirichlet (D) or Neumann (N) boundary condition. Here, and .

Corrector step

The second corrector step solves for the velocity by correcting the intermediate velocities using the pressure gradient as a Lagrange multiplier to enforce continuity. Applying divergence to the correction equations yields the pressure Poisson equation for the pressure field at half-stepwhich can be solved implicitly in two dimensions (x, y) via the following decomposition:Note the addition of a simple Hermitian to the pressure matrix , which would otherwise be singular (corank 1) under Neumann boundary conditions for the pressure field.

With the new pressure , the velocities are updated at the time-step aswhere are the gradient operators.

Overall, the variational quantum solver for Navier-Stokes equations using the projection method (see Algorithm 3) involves two sequential steps, the first requiring a number of Algorithm 1 iterations equal to the number of velocity components, and the second for the pressure Poisson step. For two-dimensional flows, the number of velocity components to be solved can be effectively reduced by one through the vorticity stream-function formulation (Appendix B). In computational fluid dynamics, these implicit systems of linear equations are often the most computationally expensive parts to solve in classical algorithms, providing incentives for potential speedup via quantum computing[49,52].

2D cavity flow

The lid-driven cavity flow is a standard benchmark for testing incompressible Navier–Stokes equations[57]. Consider a two-dimensional square domain with only one wall sliding tangentially at a constant velocity. For simplicity, we employ a fixed collocated grid, instead of a staggered grid which helps avoid spurious pressure oscillations but at the cost of increased mesh and discretization complexity. No-slip boundary conditions apply on all walls, so that zero velocity applies on all wall boundaries except one moving at .

Figure 7a shows a snapshot of a test case conducted on a grid at up to , with the central vortex shown by normalized velocity quivers in white. In terms of time complexity, we note that the pressure correction step requires a greater number of function evaluations for convergence compared to an implicit velocity step (Fig. 7b). This is due to the additional quantum circuits for evaluating the Hamiltonians (32, 33) for Neumann boundary conditions and one for specifying the reference pressure (50), leading to a total of 9 evaluation terms compared to 6, a ratio which corroborates with the apparent increase in function evaluations shown in Fig. 7b.

Figure 7

(a) 2D lid-driven cavity flow () on a grid at up to and upper boundary sliding in the x direction at . Color map indicates velocity magnitude with normalized velocity quivers in white indicating direction of flow. (b) Plots of cumulative number of function evaluations (Nfeval) vs. time for intermediate velocities (u,v) and pressure (p).

While not directly comparable to classical computational fluid dynamics in numerical accuracy, this exercise, nevertheless, roadmaps potential applications of the variational quantum method towards more complicated flow problems[52].

Conclusion

In this study, we proposed a variational quantum solver for evolution equations which include a Laplacian operator to be solved implicitly. For short time-steps , the use of initial parameter sets encoded from prior solution vectors results in faster convergence compared to random re-initialization. The overall time complexity scales with the Ansatz volume for gradient estimation and with the number of time-steps for temporal discretization. Our proposed algorithm extends naturally to higher-order time-stepping and higher dimensions. For evolution equations with non-trivial source terms, the semi-implicit scheme can be applied, where non-linear source terms are handled explicitly. Using statevector simulations, we demonstrated that variational quantum algorithms can be useful in solving popular partial differential equations, including the reaction-diffusion and the incompressible Navier-Stokes equations. Together, our proposed algorithm extends the use of quantum Poisson solvers to solve time-dependent problems with reduced time complexity from variational quantum algorithms over classical computation.

The present work aims at bridging the gap between variational quantum algorithms and practical applications. Our work has assumed that the state preparations, unitary transformations and measurements are implemented perfectly, and does not consider the effects of quantum noise from actual hardware or any potential amplification from iterative time-stepping. In our implementation, we only considered the hardware-efficient Ansatz with rotation gates and controlled-NOT entanglers, and thus leave open the question about the performance of other Ansätze. Future work can include noise mitigation[58-60], quantum random access memory[61-63], tensor networks[64], Ansatz architecture, non-linear algorithms and cost-efficient encoding.