Solving Fractional Differential Equations on a Quantum Computer:
A Variational Approach

We set up a space-time rectangular domain with spatial extent $L$ and temporal duration $T$, dividing it into a regularly spaced grid with integer size $N\times M$. The spatial grid points are defined as $x_{i}=ih$ where $h=L/N$ and $i=1,2,\dots,N$, and the temporal grid points as $t_{k}=k\tau$ where $\tau=T/M$ and $k=1,2,\dots,M$.

The first-order finite difference of the Caputo derivative of order $(0<\alpha\leq 1)$ is approximated by [5, 22]

where $u_{n}^{k}\equiv u(t_{k},x_{n})$, $g_{\alpha,\tau}:=\tau^{-\alpha}/\Gamma(2-\alpha)$ and $w_{j}^{(\alpha)}:=j^{1-\alpha}-(j-1)^{1-\alpha}$. Rearranging Eq. 3 yields a discrete sum of terms $u_{n}^{j}$ as

where $w_{1}=1$ by definition.

The finite difference of the spatial derivative is approximated by

Together, the difference scheme for the fractional diffusion equation can be iterated in time as

where $a=g_{\alpha,\tau}^{-1}h^{-2}$ and $\delta u_{n}:=u_{n+1}-2u_{n}+u_{n-1}$. In a more compact matrix form, this reads

where $\mathbf{u}^{k}:=u(t_{k},x)$, $\Delta w_{j}:=w_{j+1}-w_{j}$, and $A$ is an $N\times N$ symmetric tridiagonal matrix with $b=1+2a$ on the main diagonal, $-a$ on the adjacent off-diagonals, and terms in the corners of the matrix that are specified by the boundary conditions. More precisely, the matrix $A$ is given by

where $(c,d)=(b,-a)$ for periodic boundary condition, $(c,d)=(b,0)$ for Dirichlet boundary condition and $(c,d)=(1+a,0)$ for Neumann boundary condition. The right-hand side of Eq. 7 carries the history of solutions up to $k-1$.

Using Dirac notation, we rewrite the iterative fractional diffusion equation as

where we approximate the unnormalized quantum state $|\tilde{u}^{k}\rangle$ at time $k$ by an ansatz or trial state formed by a set of parameterized unitaries $R$ and unparameterized entangling unitaries $W$ as,

where $r^{k}$ is the norm of the unnormalized state $|\tilde{u}^{k}\rangle$, $\boldsymbol{\theta}$ is a set of gate parameters and $l$ is the number of layers in the ansatz. The unnormalized source state $|\tilde{f}^{k-1}\rangle$ is a linear combination of history states in $[0,k-1]$.

Following [11], we define a cost function based on potential energy at time $k$,

where

With that, $\langle u^{k}|\tilde{f}^{k-1}\rangle$ expands into a list of $k$ measurable overlap terms $\langle u^{k}|u^{[0,k-1]}\rangle$ [23].

Introducing $|f,u\rangle:=2^{-1/2}(|0\rangle|f\rangle+|1\rangle|u\rangle)$, we have

Here, $\langle f,u\rangle:=\langle f,u|X\otimes\mathbb{I}^{\otimes n}|f,u\rangle$ is used as a notation shorthand, where $X=|1\rangle\langle 0|+|0\rangle\langle 1|$ is the Pauli-X matrix and $\mathbb{I}:=\mathbb{I}_{0}+\mathbb{I}_{1}=|0\rangle\langle 0|+|1\rangle\langle 1|$ is the $2\times 2$ identity matrix. Also, $n=\log_{2}N$ denotes the number of qubits representing $N$ spatial grid points.

Taking the derivative of $\mathcal{C}^{k}$ with respect to $r^{k}$, the optimal norm

can be eliminated from the cost function

which can be measured using quantum circuits (Fig. 12a,b) and optimized classically as

via either gradient-based (evaluate $\partial_{\boldsymbol{\theta}^{k}}\mathcal{C}^{k}$) or gradient-free methods. A schematic of the proposed variational quantum algorithm is shown in Fig. 1.

The Hamiltonian matrix $A$ (Eq. 8) can be decomposed into

Since expectation values of the identity operator $\mathbb{I}$ are equal to 1, i.e. $\langle\phi|\mathbb{I}^{\otimes n}|\phi\rangle=\langle\phi^{\prime}|\mathbb{I}^{\otimes n}|\phi^{\prime}\rangle=1$ [11], evaluating the expectation value of $A$ requires only the evaluation of expectation values of the simple terms ($H_{1-2}$ for periodic boundary condition, $H_{1-3}$ Dirichlet boundary condition and $H_{1-4}$ for Neumann boundary condition). The operator $\mathcal{S}$ denotes the n-qubit cyclic shift operator,

which can be implemented using relative-phase Toffoli gates, Toffoli gates, a cnot gate and an X gate [24].

For problems that admit only real solutions, it is preferable to use a real-amplitude ansatz, represented by Eq. 10 formed by $l$ repeating blocks, each consisting of a parameterized $R$ layer with one $R_{Y}(\theta)$ gate on each qubit, followed by unparameterized $W$ layer with a cnot gate between consecutive qubits [25]:

where we have employed the product notation with an initial index greater than the stopping index to signify that the terms in the product are arranged in decreasing index order [26]. In addition, we consider a real-amplitude ansatz with circular entanglement, such that

which has been shown to be efficient for solving partial differential equations [16]. See Fig. 11 for circuit diagrams of the ansaetze.

To solve Eq. 9, we prepare the initial ansatz state $|\tilde{u}^{0}\rangle=r^{0}|u(\boldsymbol{\theta}^{0})\rangle$ using classical optimization on parameters $\boldsymbol{\theta}=\left(\theta_{i,j}\right)_{(i,j)\in\{1,\ldots,l\}\times\{1,\ldots,n\}}$ such that

where $|\hat{u}^{0}\rangle\equiv\mathbf{u}^{0}/\left\|\mathbf{u}^{0}\right\|$ is the normalized initial classical solution vector and $r^{0}=\left\|\mathbf{u}^{0}\right\|$ is the initial norm.

The gradient of the cost function Eq. 15 can be estimated on quantum computers using the parameter shift rule. Following [11], here we write down the partial derivative of the cost function with respect to parameter $\theta^{k}_{i,j}$ indexed by $(i,j)\in[1,nl]$ at $k$, as

where $\langle f,u\rangle_{A}:=\langle f,u|X\otimes A|f,u\rangle$ is used as a notational shorthand. The parametric shift is implemented by a $\pi$ rotation on the $i$-th $R_{Y}$ gate as

where $i\%n$ is the remainder $r\in\{0,1,\ldots,n-1\}$ obtained when $i$ is divided by $n$ and $\lfloor i/n\rfloor$ is the integer part of $i/n$. Thus, for each index $i$ at time $k$, gradient estimation requires at least $k$ overlap measurements.

Time complexity. — Here we briefly estimate the time complexity of the algorithm based on the number of time steps, quantum circuits, gates and shots required, neglecting classical overheads such as parameter update, iteration, functional evaluation and initial encoding.

For a number of time steps $M$, the overall time complexity is

where $n$ is the number of qubits, $l$ is the number of ansatz layers, $N_{\mathrm{it}}$ is the mean number of iterations required per time step and $N_{g}$ is the mean number of evaluations required for gradient estimation per iteration. Specifically, the terms in square brackets are the gate complexities for the fractional derivative $Ml$ based on $M/2$ overlap circuits containing $2l$ parametric gates in depth and for the Hamiltonian $(l+n^{2})$ based on $\mathcal{O}(1)$ circuits ($2-4$ Eq. 17) containing $l$ parametric gates in depth for the ansatz and $\mathcal{O}(n^{2})$ non-parametric gates for the shift operator. Each circuit is sampled repeatedly for $\mathcal{O}(\epsilon^{-2})$ shots per measurement, where $\epsilon$ is the desired precision. We remark that with quantum amplitude estimation [27], the number of shots required may be quadratically reduced from $\mathcal{O}(\epsilon^{-2})$ to $\mathcal{O}(\epsilon^{-1})$. This decrease, however, entails increased circuit depth, as exemplified by the $\Omega(\epsilon^{-1})$ circuit depths utilized in [27, 28, 29]. Given these challenges in reducing circuit depth in quantum amplitude estimation and its extensions [30, 31, 32, 33], we do not consider this approach in the present study, deferring it for future exploration.

For $n\sim l$, the gradient-based optimizer $N_{g}\sim\mathcal{O}(nl)$ yields an overall time complexity

scaling as $\mathcal{O}(M^{2})$ for overlap circuits, or $\mathcal{O}(M)$ for shift operators, the latter scaling applicable to the non-fractional diffusion/heat equation ($\alpha=1$, [11]). Using simultaneous perturbation stochastic approximation (SPSA) [34], only two circuit executions are required per gradient evaluation $N_{g}\sim\mathcal{O}(1)$, independent of the number of parameters, similar to gradient-free optimization.

The present algorithm is efficient in time complexity with respect to number of spatial grid points $N(=2^{n})$ [11, 35], neglecting initial encoding costs. Quantum advantage with respect to time $T$ may be achievable in higher dimensions (Remark 2 [35], Remark 3 [36]).

Memory scaling. — To solve the Caputo derivative Eq. 4, a classical algorithm assesses $\mathcal{O}(kN)$ parameters from memory at time $k$, compared to $\mathcal{O}(knl)$ parameters for variational quantum algorithm Eq. 12. Assuming $n\sim l$, the memory cost of the Caputo sums scales sub-exponentially as $\mathcal{O}(k\log N)$ (Fig. 1), demonstrating efficient algorithmic space complexity.

We solve for time-fractional $(0<\alpha\leq 1)$ sub-diffusion equation Eq. 1 as an initial value problem with initial condition $u(0,x)=x(1-x)$ and Dirichlet boundary condition $u(t,0)=u(t,1)=0$ [6].

Figure 2 shows that the variational quantum solutions agree with classical sub-diffusion solutions in $32\times 32$ space-time domain with qubit-layer count $(n,l)=(5,4)$, for (a) $\alpha=1$ and (b) $\alpha=0.5$. Figure 3 shows that low fractional-order solutions tend to diffuse faster for short times, but slow down for long times, both characteristics of sub-diffusive behavior.

The trace error at time $k$

can be averaged over the number of time steps $M$ for $0<\alpha\leq 1$, where $|\hat{u}^{k}\rangle\equiv\mathbf{u}^{k}/\left\|\mathbf{u}^{k}\right\|$ is the normalized classical solution vector at time $k$.

The fractional Burgers’ equation [40, 41] has found applications in shallow water waves and waves in bubbly liquids [42]. We consider the following 1D time-fractional Burgers’ equation [43] as an extension to the time-fractional diffusion equation Eq. 1,

where $D_{t}^{\alpha}$ is the Caputo fractional derivative of the $\alpha$-th order, as previously defined, and $\nu$ is the fluid kinematic viscosity. Readers interested in quantum computing algorithms for non-fractional Burgers’ equation ($\alpha=1$) are referred to [44, 45, 13].

A natural extension to the nonlinear differential equations involves exploring coupled systems of differential equations, as exemplified by the susceptible-exposed-infectious-recovered (SEIR) epidemic model [44]. Fractional-order epidemic models have garnered significant attention in the wake of the COVID-19 pandemic [47, 48, 49]. These models include spatio-temporal models that account for geographical spread of the disease [50, 51, 52]. A minimal representation of a fractional diffusive SEIR model is given by [53, Eq. (24)–(29)]:

which partitions a population into four distinct cohorts: susceptible individuals $S$, exposed individuals $E$, infectious individuals $I$, and those who have recovered $R$, without accounting for vaccination and quarantine. The symbols $\alpha_{i}$ and $\nu_{i}$, where $i\in\{1,2,3,4\}$, are the fractional order and the diffusion coefficient corresponding to each respective $\{S,E,I,R\}$ cohort in that order, $\Pi$ is the population influx rate, $\beta$ is the infective rate, $\mu$ is the death rate, $\sigma$ is the progression rate and $\rho$ is the recovery rate.