Quantum algorithms for linear and non-linear fractional reaction-diffusion equations

Reaction–diffusion equations arise in many areas in science and engineering [33, 9, 35, 17, 22]. In population dynamics models in biology, the reaction term typically accounts for growth, whereas the diffusion term accounts for migration [24, 31]. The classical diffusion term has its origin in a model in physics. Recent research investigations indicate that the classical diffusion equation is inadequate to model many real situations, where a particle plume spreads faster than that predicted by the classical model, and may exhibit significant asymmetry [32]. In a fractional diffusion equation, the classical Laplace operator in the spatial variable is replaced by a fractional Laplacian of order less than two $(-\Delta)^{\alpha/2}$ where $0<\alpha\leq 2$. The fundamental solutions of these equations still exhibit useful scaling properties that make them attractive for applications.

The present article deals with

$$\begin{split}\partial_{t}u(t,x)&=-(-\Delta)^{\alpha/2}u(t,x)-c(t,x)u(t,x)+au(t,x)(1-u(t,x)),\quad t\in[0,T],x\in[0,1]^{d},\\ u(0,x)&=u_{0}(x).\end{split}$$

(1)

Here $(-\Delta)^{\alpha/2}$ is the fractional Laplacian where $0<\alpha\leq 2$, and $c(t,x)$ is the potential function. For the nonlinear term, in this work we only consider the quadratic non-linearity which yields Fisher’s equation, but our results can be generalized to equations with high-order polynomial non-linearity. For the fractional Laplacian operator, there are several different definitions of $(-\Delta)^{\alpha/2}$ on bounded domain, including spectral definition and Riesz definition [28]. This works focuses on spectral fractional Laplacian with periodic boundary condition. Classical numerical algorithms for solving Equation 1 typically require exponential computational resources when the spatial variable is in high dimension: suppose that we use $N$ grid points or basis functions for spatial discretization in each dimension, then the dimension of semi-discretized differential equations becomes as large as $N^{d}$. Therefore we would like to explore the power of quantum algorithms for fractional reaction-diffusion equations and whether quantum algorithms can be efficient in high-dimensional case.

Quantum algorithms for differential equations aim at preparing a quantum state encoding the solutions at discrete grid points in its amplitudes. The first quantum differential equation algorithm was proposed in [6], which transforms differential equations into a linear system of equations using multi-step discretization and then applies quantum linear system algorithms such as HHL algorithm [21] or advanced ones [10, 36, 3, 27, 16]. Since then, there have been remarkable progresses on designing better quantum algorithms, for linear differential equations based on refined discretization [7, 12, 13, 26, 8], time-marching strategy [18], Schrödingerization [25], linear combination of unitaries technique [4], and for nonlinear differential equations using linearization techniques [29, 1]. With the caveat that the output is a quantum state encoding solutions in its amplitudes, these quantum algorithms can achieve exponential speedup in the system size compared to classical algorithms.

However, when directly applied to fractional reaction-diffusion equations, existing generic quantum algorithms are not as efficient as expected due to two major difficulties. First, similar as the standard Laplacian operator, fractional Laplacian operator is an unbounded operator, so its spatially discretized version has a huge spectral norm as the spatial dimension and the number of grid points increase. Most existing quantum differential equation algorithms scale at least linearly on the spectral norm of the coefficient matrix, and thus can be computationally expensive in high-dimensional case for accurate simulation. Second, when the equation is genuinely fractional (i.e., $\alpha<2$), the coefficient matrix for the linear part after spatial discretization is unavoidably dense, because fractional differential operators are global operators that depends on the function evaluated in the entire space. The dense coefficient matrix poses computational difficulties in solving nonlinear equations. This is because all the existing quantum Carleman linearization algorithms require the coefficient matrix to be sparse in order to bypass the difficulty caused by the enlarged Carleman matrix, which is a direct sum of matrices in different dimensions [29].

We remark that the (fractional- or integer-order) Laplacian operator, and more general spatial differential operators, widely appear in various types of partial differential equations. In the contexts other than fractional reaction-diffusion equations, there have been several work managing to overcome the computational difficulty brought by its large spectral norm after spatial discretization, but unfortunately those techniques do not apply to the fractional reaction-diffusion equations. For example, in real-space Schrödinger equation, a poly-logarithmic dependence on the spectral norm of the Laplacian operator can be achieved by simulating the Hamiltonian in the interaction picture [30, 11], which simulates the transformed wavefunction under the rotation associated with the Laplacian operator. The resulting interaction picture Hamiltonian becomes bounded and thus can be efficiently simulated. To avoid the large spectral norm dependence in the rotations, the algorithm takes advantage of an important feature of the Laplacian operator: it can be diagonalized by the quantum Fourier transform (QFT) circuit and thus can be fastforwarded (i.e., Hamiltonian simulation governed by the Laplacian operator can be implemented with cost independent of evolution time and its spectral norm). However, such a technique does not directly work for fractional reaction-diffusion equations. The interaction picture transformation requires both forward and backward time evolution of the Laplacian operator. This is efficient for Schrödinger equations because the dynamics is reversible and both forward and backward time evolution operators are unitary, but fractional reaction-diffusion equation is a dissipative system and implementing its backward time evolution can be prohibitively expensive. The recent work [5] proposes efficient quantum algorithms for various partial differential equations beyond Schrödinger equations. The key technique there is a generalization of the fastforwarding simulation of the Laplacian operator. However, the algorithms in [5] require the entire coefficient matrix of the linear part to be fast-forwardable (i.e., its spectral decomposition has quantumly implementable eigenstates and classically computable eigenvalues), which cannot be satisfied even by linear fractional reaction-diffusion equations with the presence of the potential $c(t,x)$.

In this work, we investigate efficient quantum algorithms for fractional reaction-diffusion equations. The majority of our work is devoted to linear fractional reaction-diffusion equations. We numerically treat the equations by the method of lines, i.e., first discretizing the spatial variable to obtain a system of ODEs

$$\frac{d}{dt}\vec{u}=-B\vec{u}-C(t)\vec{u},$$

(2)

and then solving the resulting ODE system with different quantum ODE algorithms. Here $\vec{u}$ represents the solution evaluated at different spatial grid points, $B$ is the discretized fractional Laplacian operator, and $C(t)$ is the potential matrix. For time evolution, we consider four different algorithms: second-order Trotter formula [14], time-marching method [18], truncated Dyson series method [8], and linear combination of Hamiltonian simulation in the interaction picture (LCHS-IP). We analyze the complexity of these four methods and the results are shown in Table 1. Our main results and contributions are summarized as follows:

1. 1.

   Second-order Trotter formula: unlike Hamiltonian simulation, Equation 2 is not a unitary dynamics, and in the Trotter formula we need to implement non-unitary operators $e^{-Bs}$ and $e^{-C(t)s}$. We discuss efficient construction of these operators via controlled rotations, and implementing their multiplication through a generalization of the compression gadget technique [30, 18] with only poly-logarithmic many ancilla qubits. For complexity analysis, we derive an improved Trotter error bound that avoids the exponential factor in [14] and generalizes to the time-dependent case. Compared to time-marching and truncated Dyson series methods, second-order Trotter has better dependence on the dimension $d$ thanks to its commutator scalings, but has worse dependence on the precision.

2. 2.

   Time-marching method: we directly apply the standard time-marching method in [18] and analyze its complexity for fractional reaction-diffusion equations. It has low state preparation cost and poly-logarithmic dependence on precision, but has worse scalings in the dimension $d$ and the evolution time $T$.

3. 3.

   Truncated Dyson series method: we directly apply the standard truncated Dyson series method in [8] and analyze its complexity for fractional reaction-diffusion equations. It still depends polynomially on the dimension due to its spectral norm dependence and has high state preparation cost due to the usage of quantum linear system algorithms, but can achieve poly-logarithmic scaling in the precision and linear scaling in time simultaneously.

4. 4.

   LCHS-IP: The LCHS-IP method is a novel method that combines the linear combination of Hamiltonian simulation (LCHS) [4] technique with the interaction picture Hamiltonian simulation [30]. The LCHS method first represents the evolution operator of Equation 2 as a linear combination of several Hamiltonian simulation problems associated with the matrices $B$ and $C(t)$. To avoid the computational overhead brought by the discretized fractional Laplacian $B$, we implement each Hamiltonian simulation in the interaction picture by rotating the Hamiltonian with respect to $B$. Therefore, the resulting algorithm only has poly-logarithm dependence on the dimension $d$ and thus is the most preferable algorithm in the high-dimensional case. It also has low state preparation cost, but only linear scaling in precision.

For nonlinear fractional reaction-diffusion equations, we discuss a block-encoding version of the Carleman linearization technique to deal with the dense coefficient matrix. Let $\vec{u}(t)$ denote the solution vector of the fractional reaction-diffusion equation after spatial discretization. Standard quantum Carleman linearization algorithm [29] considers the dynamics of the enlarged vector $[\vec{u}(t);\vec{u}(t)^{\otimes 2};\cdots;\vec{u}(t)^{\otimes M}]$, which approximately satisfies a linear system of ODEs governed by a so-called Carleman matrix. Since $\vec{u}(t)^{\otimes m}$’s are in different sizes for different powers, the Carleman matrix is a direct sum of matrices in different dimensions. In our work, we propose a simple generalization of the Carleman linearization, by extending the Carleman matrix to even higher dimension such that it becomes the direct sum of matrices in the same dimension and the corresponding solution vector $[\vec{w}_{1};\vec{w}_{2};\cdots;\vec{w}_{M}]$ can be exactly mapped to the solution $[\vec{u}(t);\vec{u}(t)^{\otimes 2};\cdots;\vec{u}(t)^{\otimes M}]$ via discarding all the zero entries. Then the extended Carleman matrix can be easily block-encoded through the block-encoding of the original coefficient matrix, and thus the corresponding extended linearized system can be solved via quantum linear differential algorithms with block-encoding input model (e.g., [18, 8, 4]). This makes the quantum Carleman linearization algorithm applicable to the equations with dense coefficient matrices, including the fractional reaction-diffusion equations.

The rest of this paper is organized as follows. Section 2 first discusses the mathematical setup of the fractional reaction-diffusion equations, notations being used throughout the paper and some preliminary results. The main results of this paper start by Section 3 with the simplest case where the equations only involve a single fractional Laplacian operator without potential or nonlinear term, then Section 4 considers the general linear equations with potential. Quantum algorithm for nonlinear equations is presented in Section 5, followed by conclusion and open questions in Section 6.

We start with a more rigorous setup of the fractional reaction-diffusion equation we are interested in and a summary of theoretical tools and technical lemmas being used in our analysis.

We start with the simplest case where $c(t,x)=0$ and $a=0$. In this case, the right hand side of Equation 1 only involves a fractional Laplacian operator. By the spectral definition, this operator has closed-form eigenvalues and known eigenfunctions, so the corresponding time-evolution operator can be implemented fast-forwardly [5] (i.e., query complexity is independent of the spectral norm and the evolution time). Here we present our algorithm in general high dimension and establish rigorous complexity estimate taking into consideration the spatial discretization errors.

We consider general $d$-dimensional case. The eigenvalues and eigenfunctions of $(-\Delta)$ are $4\pi^{2}(k_{0}^{2}+\cdots+k_{d-1}^{2})$ and $e^{2\pi i(k_{0}x_{0}+\cdots+k_{d-1}x_{d-1})}$, where $x=(x_{0},\cdots,x_{d-1})$ denotes the spatial coordinate and $k=(k_{0},\cdots,k_{d-1})$ denotes a set of integers. Let the Fourier series of $u_{0}(x)$ be

$$u_{0}(x)=\sum_{k\in\mathbb{Z}^{d}}\hat{u}_{k}e^{2\pi i(k_{0}x_{0}+\cdots+k_{d-1}x_{d-1})},$$

(7)

where $\hat{u}_{k}$ denotes the Fourier coefficients. Then the solution of Equation 1 has the form

$$u(T,x)=\sum_{k\in\mathbb{Z}^{d}}\hat{u}_{k}e^{-(4\pi^{2}(k_{0}^{2}+\cdots+k_{d-1}^{2}))^{\alpha/2}T}e^{2\pi i(k_{0}x_{0}+\cdots+k_{d-1}x_{d-1})}.$$

(8)

Numerical solutions can be obtained by truncating the Fourier series at a finite order.

Equation 8 can be quantumly implemented using quantum Fourier transform (QFT) and controlled rotations. Recall that our goal is to prepare an approximation of the quantum state encoding the normalized solution at discrete spatial grid points $(j/N)_{j\in[N]^{d}}$. Let $N$ be the number of the grid points in each spatial dimension. Suppose that we are given the oracle $O_{u_{0}}$ that prepares the normalized initial condition

$$\ket{u_{0}}=\frac{1}{\|\vec{u}_{0}\|}\sum_{n\in[N]^{d}}u_{0}(n_{0}/N,\cdots,n_{d-1}/N)\ket{n_{0}}\cdots\ket{n_{d-1}},$$

(9)

where

$$\vec{u}_{0}=\sum_{n\in[N]^{d}}u_{0}(n_{0}/N,\cdots,n_{d-1}/N)\ket{n_{0}}\cdots\ket{n_{d-1}}.$$

(10)

We first compute a quantum state encoding the Fourier coefficients $\hat{u}_{k}$ by QFT. Specifically, let $\omega_{N}=e^{2\pi i/N}$ and $\mathcal{F}$ denote the one-dimensional QFT, i.e., for any computational basis state $\ket{j}$,

$$\mathcal{F}\ket{j}=\frac{1}{\sqrt{N}}\sum_{l=0}^{N-1}\omega_{N}^{jl}\ket{l}.$$

(11)

Then

	$\displaystyle(\mathcal{F}^{-1})^{\otimes d}\vec{u}_{0}$	$\displaystyle=\frac{1}{N^{d/2}}\sum_{n\in[N]^{d}}u_{0}(n_{0}/N,\cdots,n_{d-1}/N)\left(\sum_{m_{0}=0}^{N-1}\omega_{N}^{-n_{0}m_{0}}\ket{m_{0}}\right)\cdots\left(\sum_{m_{d-1}=0}^{N-1}\omega_{N}^{-n_{d-1}m_{d-1}}\ket{m_{d-1}}\right)$		(12)
		$\displaystyle=\frac{1}{N^{d/2}}\sum_{m\in[N]^{d}}\sum_{n\in[N]^{d}}u_{0}(n_{0}/N,\cdots,n_{d-1}/N)\omega_{N}^{-n_{0}m_{0}-\cdots-n_{d-1}m_{d-1}}\ket{m_{0}}\cdots\ket{m_{d-1}}$		(13)
		$\displaystyle=\frac{1}{N^{d/2}}\sum_{m\in[N]^{d}}\sum_{n\in[N]^{d}}\sum_{k\in\mathbb{Z}^{d}}\hat{u}_{k}e^{2\pi i(k_{0}n_{0}/N+\cdots+k_{d-1}n_{d-1}/N)}\omega_{N}^{-n_{0}m_{0}-\cdots-n_{d-1}m_{d-1}}\ket{m_{0}}\cdots\ket{m_{d-1}}$		(14)
		$\displaystyle=\frac{1}{N^{d/2}}\sum_{m\in[N]^{d}}\sum_{k\in\mathbb{Z}^{d}}\hat{u}_{k}\left(\sum_{n_{0}=0}^{N-1}\omega_{N}^{n_{0}(k_{0}-m_{0})}\right)\cdots\left(\sum_{n_{d-1}=0}^{N-1}\omega_{N}^{n_{d-1}(k_{d-1}-m_{d-1})}\right)\ket{m_{0}}\cdots\ket{m_{d-1}}$		(15)
		$\displaystyle=N^{d/2}\sum_{m\in[N]^{d}}\sum_{j\in\mathbb{Z}^{d}}\hat{u}_{m+jN}\ket{m_{0}}\cdots\ket{m_{d-1}}.$		(16)

When the function $u_{0}(x)$ satisfies certain regularity assumption (to be specified later), the summation $\sum_{j\in\mathbb{Z}^{d}}\hat{u}_{m+jN}$ is dominated by the index with smallest absolute value, because the Fourier coefficients decay rapidly with respect to its frequency. Therefore

$$(\mathcal{F}^{-1})^{\otimes d}\ket{\vec{u}_{0}}\approx\frac{N^{d/2}}{\|\vec{u}_{0}\|}\sum_{m\in[N]^{d}}\hat{u}_{i(m)}\ket{m_{0}}\cdots\ket{m_{d-1}},$$

(17)

where $i(m)=(i_{0}(m),\cdots,i_{d-1}(m))$ represents an $d$-dimensional vector with the $j$-th index

$$i_{j}(m)=\begin{cases}m_{j},&\text{ if }0\leq m_{j}\leq N/2\\ m_{j}-N,&\text{ if }N/2+1\leq m_{j}\leq N-1.\end{cases}$$

(18)

Now we append three ancilla registers to Equation 17 and get the (approximate) quantum state

$$\frac{N^{d/2}}{\|\vec{u}_{0}\|}\sum_{m\in[N]^{d}}\hat{u}_{i(m)}\ket{m_{0}}\cdots\ket{m_{d-1}}\ket{0}\ket{0}\ket{0}.$$

(19)

The first two ancilla registers are used to binarily encode $e^{-(4\pi^{2}(i_{0}(m)^{2}+\cdots+i_{d-1}(m)^{2}))^{\alpha/2}T}$ and the third register is used for rotation. Specifically, suppose that we are given the oracle that encodes the eigenvalues as

$$O_{1}:\ket{m_{0}}\cdots\ket{m_{d-1}}\ket{0}\rightarrow\ket{m_{0}}\cdots\ket{m_{d-1}}\ket{(2\pi\|i(m)\|)^{\alpha}},$$

(20)

and the oracle for computing an exponential function as

$$O_{\exp,1}:\ket{x}\ket{0}\rightarrow\ket{x}\ket{e^{-xT}}.$$

(21)

Notice that both functions have closed-form expression so the oracles can be efficiently constructed using classical arithmetic [34]. Applying $O_{1}$ on the index and first ancilla registers and then $O_{\exp,1}$ on the first and second ancilla registers yields

$$\frac{N^{d/2}}{\|\vec{u}_{0}\|}\sum_{m\in[N]^{d}}\hat{u}_{i(m)}\ket{m_{0}}\cdots\ket{m_{d-1}}\ket{\|i(m)\|}\ket{e^{-(2\pi\|i(m)\|)^{\alpha}T}}\ket{0}.$$

(22)

Let $\text{c-}R$ denote the controlled rotation

$$\text{c-}R:\ket{\theta}\ket{0}\rightarrow\ket{\theta}\left(\theta\ket{0}+\sqrt{1-|\theta|^{2}}\ket{1}\right).$$

(23)

Then, by applying $\text{c-}R$ on the last two registers in Equation 22, we obtain

$$\frac{N^{d/2}}{\|\vec{u}_{0}\|}\sum_{m\in[N]^{d}}\hat{u}_{i(m)}e^{-(2\pi\|i(m)\|)^{\alpha}T}\ket{m_{0}}\cdots\ket{m_{d-1}}\ket{\|i(m)\|}\ket{e^{-(2\pi\|i(m)\|)^{\alpha}T}}\ket{0}+\ket{\perp},$$

(24)

where $\ket{\perp}$ represents the orthogonal part with $1$ in the last ancilla register. Uncomputing the first two ancilla registers yields

$$\frac{N^{d/2}}{\|\vec{u}_{0}\|}\sum_{m\in[N]^{d}}\hat{u}_{i(m)}e^{-(2\pi\|i(m)\|)^{\alpha}T}\ket{m_{0}}\cdots\ket{m_{d-1}}\ket{0}\ket{0}\ket{0}+\ket{\perp},$$

(25)

Notice that, by replacing $u_{0}(x)$ with $u(T,x)$ in Equation 17 and using Equation 8, we have

$$(\mathcal{F}^{-1})^{\otimes d}\vec{u}(T,x)\approx N^{d/2}\sum_{m\in[N]^{d}}\hat{u}_{i(m)}e^{-(2\pi\|i(m)\|)^{\alpha}T}\ket{m_{0}}\cdots\ket{m_{d-1}}.$$

(26)

Therefore, applying the QFT $\mathcal{F}^{\otimes d}$ to Equation 25 approximately gives

$$\frac{1}{\|\vec{u}_{0}\|}\sum_{n\in[N]^{d}}u(T,n/N)\ket{n_{0}}\cdots\ket{n_{d-1}}\ket{0}\ket{0}\ket{0}+\ket{\perp}.$$

(27)

Measuring the ancilla registers to get all $0$’s yields an approximation of $\ket{u(T)}$, and the averaged number of repeats for success after amplitude amplification scales $\mathcal{O}(\|\vec{u}_{0}\|/\|\vec{u}(T)\|)$.

We summarize the quantum algorithm in Figure 2 and present its error and complexity estimates in the following result. The main result is that the overall error will decrease as the number of the grid points increases, since more grid points imply larger truncation order in the Fourier series. Detailed proof can be found in Appendix B.

Consider the linear spatial fractional reaction-diffusion equations

$$\begin{split}\partial_{t}u(t,x)&=-(-\Delta)^{\alpha/2}u(t,x)-c(t,x)u(t,x),\quad t\in[0,T],x\in[0,1]^{d},\\ u(0,x)&=u_{0}(x).\end{split}$$

(29)

Throughout this section, we assume the potential function $c(t,x)$ to be non-negative. We will first discuss the validity of this assumption by so-called shifting equivalence. The equation Equation 29 can then be treated by the method of lines, i.e., first discretizing the spatial variable to obtain a system of ODEs, and then solving the ODE with different quantum algorithms. We consider four different quantum algorithms: second-order Trotter formula, time-marching method, truncated Dyson series method, and linear combination of Hamiltonian simulation in the interaction picture.