Efficient quantum algorithm for dissipative nonlinear differential equations

Liu, Jin-Peng; Kolden, Herman Øie; Krovi, Hari K.; Loureiro, Nuno F.; Trivisa, Konstantina; Childs, Andrew M.

doi:10.1073/pnas.2026805118

Efficient quantum algorithm for dissipative nonlinear differential equations

Jin-Peng Liu https://orcid.org/0000-0002-1885-8879, Herman Øie Kolden https://orcid.org/0000-0003-4121-389X, Hari K. Krovi, +2 , Nuno F. Loureiro, Konstantina Trivisa, and Andrew M. Childs https://orcid.org/0000-0002-9903-837X [email protected]Authors Info & Affiliations

Edited by Anthony Leggett, University of Illinois at Urbana–Champaign, Urbana, IL, and approved July 19, 2021 (received for review March 6, 2021)

August 26, 2021

118 (35) e2026805118

https://doi.org/10.1073/pnas.2026805118

PDF/EPUB

Significance

Nonlinear differential equations appear in many domains and are notoriously difficult to solve. Whereas previous quantum algorithms for general nonlinear differential equations have complexity exponential in the evolution time, we give the first quantum algorithm for dissipative nonlinear differential equations that is efficient provided the dissipation is sufficiently strong relative to nonlinear and forcing terms and the solution does not decay too rapidly. We also establish a lower bound showing that differential equations with sufficiently weak dissipation have worst-case complexity exponential in time, giving an almost tight classification of the quantum complexity of simulating nonlinear dynamics. Furthermore, numerical results for the Burgers equation suggest that our algorithm may potentially address complex nonlinear phenomena even in regimes with weaker dissipation.

Abstract

Nonlinear differential equations model diverse phenomena but are notoriously difficult to solve. While there has been extensive previous work on efficient quantum algorithms for linear differential equations, the linearity of quantum mechanics has limited analogous progress for the nonlinear case. Despite this obstacle, we develop a quantum algorithm for dissipative quadratic

n

-dimensional ordinary differential equations. Assuming

R < 1

, where

R

is a parameter characterizing the ratio of the nonlinearity and forcing to the linear dissipation, this algorithm has complexity

T^{2} q poly (\log T, \log n, \log 1 / ϵ) / ϵ

, where

T

is the evolution time,

ϵ

is the allowed error, and

q

measures decay of the solution. This is an exponential improvement over the best previous quantum algorithms, whose complexity is exponential in

T

. While exponential decay precludes efficiency, driven equations can avoid this issue despite the presence of dissipation. Our algorithm uses the method of Carleman linearization, for which we give a convergence theorem. This method maps a system of nonlinear differential equations to an infinite-dimensional system of linear differential equations, which we discretize, truncate, and solve using the forward Euler method and the quantum linear system algorithm. We also provide a lower bound on the worst-case complexity of quantum algorithms for general quadratic differential equations, showing that the problem is intractable for

R \geq \sqrt{2}

. Finally, we discuss potential applications, showing that the

R < 1

condition can be satisfied in realistic epidemiological models and giving numerical evidence that the method may describe a model of fluid dynamics even for larger values of

R

.

Models governed by both ordinary differential equations (ODEs) and partial differential equations (PDEs) arise extensively in the natural and social sciences, medicine, and engineering. Such equations characterize physical and biological systems that exhibit a wide variety of complex phenomena, including turbulence and chaos. We focus here on differential equations with nonlinearities that can be expressed with quadratic polynomials, which include many archetypal models in biology, fluid dynamics, and plasma physics.

Quantum algorithms offer the prospect of rapidly characterizing the solutions of high-dimensional systems of linear ODEs (1–3) and PDEs (4–10). Such algorithms can produce a quantum state proportional to the solution of a sparse (or block-encoded)

n

-dimensional system of linear differential equations in time

poly (\log n)

using the quantum linear system algorithm (QLSA) (11).

Early work on quantum algorithms for differential equations already considered the nonlinear case (12). This work gave a quantum algorithm for ODEs that simulates polynomial nonlinearities by storing multiple copies of the solution. The complexity of this approach is polynomial in the logarithm of the dimension but exponential in the evolution time, scaling as

O (1 / ϵ^{T})

. While some heuristic quantum algorithms for nonlinear ODEs have been studied more recently (13, 14), these works do not make precise statements about concrete implementations or running times of quantum algorithms. The recent preprint (15) also describes a quantum algorithm to solve a nonlinear ODE by linearizing it using a different approach from the one taken here. However, a proof of correctness of their algorithm involving a bound on the condition number and probability of success is not given. The authors also do not describe how barriers such as those of ref. 16 could be avoided in their approach.

Constructing efficient quantum algorithms for general classes of nonlinear dynamics has been considered largely out of reach since the linearity of quantum mechanics makes it challenging to efficiently represent nonlinear dynamics. Furthermore, nonlinear modifications of quantum mechanics generically enable quickly solving hard computational problems such as performing unstructured search among

n

items in time

poly (\log n)

. Indeed, the search lower bound (17) can be used to show that nonlinear dynamics is exponentially difficult to simulate in general (16, 18, 19).

In this article, we design and analyze a quantum algorithm that overcomes this limitation using Carleman linearization (20–22). This approach embeds polynomial nonlinearities into an infinite-dimensional system of linear ODEs and then truncates it to obtain a finite-dimensional linear approximation. We discretize the finite ODE system in time using the forward Euler method and solve the resulting linear equations with the QLSA (11, 23). Subject to the condition

R < 1

, where the quantity

R

(defined in Problem 1 below) characterizes the relative strength of the nonlinear and dissipative linear terms, we show that the total complexity of this quantum Carleman linearization (QCL) algorithm is

s T^{2} q poly (\log T, \log n, \log 1 / ϵ) / ϵ

, where

s

is the sparsity,

T

is the evolution time,

q

is the ratio of the norm of the initial vector to that of the solution vector at time

T

,

n

is the dimension, and

ϵ

is the allowed error (Theorem 1). In the regime

R < 1

, this is an exponential improvement over ref. 12, which has complexity exponential in

T

.

Note that the solution cannot decay exponentially in

T

for the algorithm to be efficient, as captured by the dependence of the complexity on

q

—a known limitation of quantum ODE algorithms (2). For homogeneous ODEs with

R < 1

, the solution necessarily decays exponentially in time (SI Appendix, Eq. S30), so the algorithm is not asymptotically efficient. However, even for solutions with exponential decay, we still find an improvement over the best previous result

O (1 / ϵ^{T})

(12) for sufficiently small

ϵ

. Thus, our algorithm might provide an advantage over classical computation for studying evolution for short times. More significantly, our algorithm can handle inhomogeneous quadratic ODEs, for which it can remain efficient in the long-time limit since the solution can remain asymptotically nonzero (for an explicit example, see the discussion just before the proof of Lemma 1 in SI Appendix, Section 3.A.1) or can decay slowly [i.e.,

q

can be

poly (T)

]. Inhomogeneous equations arise in many applications, including, for example, the discretization of PDEs with nontrivial boundary conditions.

We also show that the requirement

R < 1

cannot be significantly improved in general (Theorem 2). Specifically, we provide a quantum lower bound showing worst-case hardness of simulating strongly nonlinear dynamics for

R \geq \sqrt{2}

, establishing a limitation on the ability of quantum computers to simulate general nonlinear dynamics.

Our quantum algorithm could potentially be applied to study models governed by quadratic ODEs arising in biology and epidemiology as well as in fluid and plasma dynamics. In particular, the celebrated Navier–Stokes equation with linear damping, which describes many physical phenomena, can be treated by our approach if the Reynolds number is sufficiently small. While the formal validity of our arguments assumes

R < 1

, we find in numerical experiments that our proposed approach remains valid for larger

R

in certain cases.

We emphasize that as in related quantum algorithms for linear algebra and differential equations, instantiating our approach requires an implicit description of the problem that allows for efficient preparation of the initial state and implementation of the dynamics. Furthermore, since the output is encoded in a quantum state, readout is restricted to features that can be revealed by efficient quantum measurements. More work remains to understand how these methods might be applied, as we discuss further in Discussion.

Problem Statement

We focus on an initial value problem described by the

n

-dimensional quadratic ODE

\frac{d u}{d t} = F_{2} u^{\otimes 2} + F_{1} u + F_{0} (t), u (0) = u_{i n} .

[1]

Here

u = {[u_{1}, \dots, u_{n}]}^{T} \in R^{n}

,

u^{\otimes 2} = {[u_{1}^{2}, u_{1} u_{2}, \dots, u_{1} u_{n}, u_{2} u_{1}, \dots, u_{n} u_{n - 1}, u_{n}^{2}]}^{T} \in R^{n^{2}}

, each

u_{j} = u_{j} (t)

is a function of

t

on the interval

[0, T]

for

j \in [n] ≔ {1, \dots, n}

,

F_{2} \in R^{n \times n^{2}}

,

F_{1} \in R^{n \times n}

are time-independent matrices, and the inhomogeneity

F_{0} (t) \in R^{n}

is a

C^{1}

continuous function of

t

.

The main computational problem we consider is as follows.

Problem 1. In the quantum quadratic ODE problem, we consider an

n

-dimensional quadratic ODE as in [1]. We assume

F_{2}

,

F_{1}

, and

F_{0}

are

s

-sparse (i.e., have at most

s

nonzero entries in each row and column),

F_{1}

is diagonalizable, and the eigenvalues

λ_{j}

of

F_{1}

satisfy

Re (λ_{n}) \leq \dots \leq Re (λ_{1}) < 0

. We parametrize the problem in terms of the quantity

R ≔ \frac{1}{| Re (λ_{1}) |} (‖ u_{i n} ‖ ‖ F_{2} ‖ + \frac{‖ F_{0} ‖}{‖ u_{i n} ‖}) .

[2]

For some given

T > 0

, we assume the values

Re (λ_{1})

,

‖ F_{2} ‖

,

‖ F_{1} ‖

,

‖ F_{0} (t) ‖

for each

t \in [0, T]

, and

‖ F_{0} ‖ ≔ max_{t \in [0, T]} ‖ F_{0} (t) ‖

,

‖ F_{0}^{'} ‖ ≔ max_{t \in [0, T]} ‖ F_{0}^{'} (t) ‖

are known, and that we are given oracles

O_{F_{2}}

,

O_{F_{1}}

, and

O_{F_{0}}

that provide the locations and values of the nonzero entries of

F_{2}

,

F_{1}

, and

F_{0} (t)

for any specified

t

, respectively, for any desired row or column. We are also given the value

‖ u_{i n} ‖

and an oracle

O_{x}

that maps

|00 \dots 0⟩ \in C^{n}

to a quantum state proportional to

u_{i n}

. Our goal is to produce a quantum state proportional to

u (T)

for some given

T > 0

within some prescribed

ℓ^{2}

error tolerance

ϵ > 0

.

When

F_{0} (t) = 0

(i.e., the ODE is homogeneous), the quantity

R = \frac{‖ u_{i n} ‖ ‖ F_{2} ‖}{| Re (λ_{1}) |}

characterizes the ratio of the nonlinearity to the linear dissipation. It is qualitatively similar to the Reynolds number in fluid dynamics.

Our algorithm uses a rescaling of Problem 1 that we introduce in SI Appendix, Section 1.

Main Results

For Problem 1 with

R < 1

, we present a quantum algorithm that improves the complexity from an exponential dependence on

T

to a nearly quadratic dependence. This algorithm is based on the concept of Carleman linearization (20–22). Carleman linearization is a method for converting a finite-dimensional system of nonlinear differential equations into an infinite-dimensional linear one. This is achieved by introducing powers of the variables into the system, allowing it to be written as an infinite sequence of coupled linear differential equations. We then truncate the system to

N

equations, where the truncation level

N

depends on the allowed approximation error, giving a finite linear ODE system.

Given a system of quadratic ODEs [1], we apply the Carleman procedure to obtain the system of linear ODEs

\frac{d ŷ}{d t} = A (t) ŷ + b (t), ŷ (0) = ŷ_{i n}

[3]

with the tridiagonal block structure

\frac{d}{d t} (\begin{matrix} ŷ_{1} \\ ŷ_{2} \\ ⋮ \\ ŷ_{N} \end{matrix}) = (\begin{matrix} A_{1}^{1} & A_{2}^{1} \\ A_{1}^{2} & A_{2}^{2} & ⋱ \\ ⋱ & ⋱ & A_{N}^{N - 1} \\ A_{N - 1}^{N} & A_{N}^{N} \end{matrix}) (\begin{matrix} ŷ_{1} \\ ŷ_{2} \\ ⋮ \\ ŷ_{N} \end{matrix}) + (\begin{matrix} F_{0} (t) \\ 0 \\ ⋮ \\ 0 \end{matrix}),

[4]

where

ŷ_{j} = u^{\otimes j} \in R^{n^{j}}

,

ŷ_{i n} = [u_{i n}; u_{i n}^{\otimes 2}; \dots; u_{i n}^{\otimes N}]

, and

A_{j + 1}^{j} \in R^{n^{j} \times n^{j + 1}}

,

A_{j}^{j} \in R^{n^{j} \times n^{j}}

,

A_{j - 1}^{j} \in R^{n^{j} \times n^{j - 1}}

for

j \in [N]

satisfying

A_{j + 1}^{j} = F_{2} \otimes I^{\otimes j - 1} + I \otimes F_{2} \otimes I^{\otimes j - 2} + \dots + I^{\otimes j - 1} \otimes F_{2},

[5]

A_{j}^{j} = F_{1} \otimes I^{\otimes j - 1} + I \otimes F_{1} \otimes I^{\otimes j - 2} + \dots + I^{\otimes j - 1} \otimes F_{1},

[6]

A_{j + 1}^{j} = F_{0} (t) \otimes I^{\otimes j - 1} + I \otimes F_{0} (t) \otimes I^{\otimes j - 2} + \dots + I^{\otimes j - 1} \otimes F_{0} (t) .

[7]

To construct a system of linear equations, we next divide

[0, T]

into

m = T / h

time steps and apply the forward Euler method on [3], letting

y^{k + 1} = [I + A (k h) h] y^{k} + b (k h)

[8]

where

y^{k} \in R^{n + n^{2} + \dots + n^{N}}

approximates

ŷ (k h)

for each

k \in {[m + 1]}_{0} ≔ {0,1, \dots, m}

, with

y^{0} = y_{i n} ≔ ŷ (0) = ŷ_{i n}

, and letting all

y^{k}

be equal for

k \in {[m + p + 1]}_{0} \ {[m + 1]}_{0}

, for some sufficiently large integer

p

. (It is unclear whether another discretization could improve performance, as discussed further in Discussion.)

In the above system, the first

n

components of

y^{k}

for

k \in {[m + p + 1]}_{0}

(i.e.,

y_{1}^{k}

) is

ϵ

-close to the exact solution

u (T)

, up to normalization. We apply the high-precision QLSA (23) to [8] and postselect on

k

to produce

y_{1}^{k} / ‖ y_{1}^{k} ‖

for some

k \in {[m + p + 1]}_{0} \ {[m]}_{0}

. SI Appendix, Section 2, gives more details about the system of linear equations and the measurement procedure.

We control the additive approximation error

ϵ

by combining a convergence theorem for Carleman linearization that improves previous analysis (22) with a bound for the global error of the Euler method. To analyze the complexity of this approach, we upper bound the condition number of the linear system and lower bound the success probability of the final measurement. We state our main algorithmic result as follows.

Theorem 1. (Informal version.) Consider an instance of the quantum quadratic ODE problem as defined in Problem 1, with its Carleman linearization as defined in [3]. Assume

R < 1

. Let

q ≔ ‖ u_{i n} ‖ / ‖ u (T) ‖

. There exists a quantum algorithm producing a state that approximates

u (T) / ‖ u (T) ‖

with error at most

ϵ \leq 1

, succeeding with probability

Ω (1)

, with query and gate complexity

\frac{s T^{2} q}{ϵ} \cdot poly (\log T, \log n, \log 1 / ϵ) .

[9]

SI Appendix, Section 2, states the formal version of Theorem 1 with more details.

We also provide a quantum lower bound for the worst-case complexity of simulating strongly nonlinear dynamics, showing that the problem is intractable for

R \geq \sqrt{2}

. Following the approach of refs. 16, 18, we construct a protocol for distinguishing two states of a qubit driven by a certain quadratic ODE. Provided

R \geq \sqrt{2}

, this procedure distinguishes states with overlap

1 - ϵ

in time

poly (\log (1 / ϵ))

. Since nonorthogonal quantum states are hard to distinguish, this implies a lower bound on the complexity of the quantum ODE problem. This gives the following limitation on quantum algorithms for nonlinear ODEs.

Theorem 2. Assume

R \geq \sqrt{2}

. Then there is an instance of the quantum quadratic ODE problem defined in Problem 1 such that any quantum algorithm for producing a quantum state approximating

u (T) / ‖ u (T) ‖

with bounded error must have worst-case time complexity exponential in

T

.

Algorithm Analysis

In this section we describe the main ingredients used to prove Theorem 1.

First, we provide an upper bound for the error from Carleman linearization for arbitrary evolution time

T

. To the best of our knowledge, the first and only explicit bound on the error of Carleman linearization appears in ref. 22. However, the authors only consider homogeneous quadratic ODEs, and to bound the error for arbitrary

T

, the logarithmic norm of

F_{1}

is assumed to be negative (theorems 4.2 and 4.3 of ref. 22), which is too strong for our case. Instead, we give an analysis under milder conditions, providing a convergence guarantee for general inhomogeneous quadratic ODEs.

Lemma 1. Consider an instance of the quadratic ODE [1], with its corresponding Carleman linearization as defined in [3]. As in Problem 1, assume that the eigenvalues

λ_{j}

of

F_{1}

satisfy

Re (λ_{n}) \leq \dots \leq Re (λ_{1}) < 0

. Then for any

j \in [N]

, the error

η_{j} (t) ≔ u^{\otimes j} (t) - ŷ_{j} (t)

satisfies

‖ η_{j} (t) ‖ \leq t N ‖ F_{2} ‖ {‖ u_{i n} ‖}^{N + 1}

[10]

where we require

R < 1

as defined in [2].

We present the proof of Lemma 1 (and an improved version for homogeneous equations) and explain the necessity of our assumptions in SI Appendix, Section 3.A.1.

Next, we provide an upper bound for the error incurred by approximating [3] with the forward Euler method.

Lemma 2. Consider an instance of the quantum quadratic ODE problem as defined in Problem 1, with

R < 1

as defined in [2]. Choose a sufficiently small time step

h

. Suppose the error from Carleman linearization satisfies

‖ η (t) ‖ \leq ‖ u (T) ‖ / 4

. Then the global error from the forward Euler method [8] on the interval

[0, T]

for [3] satisfies

‖ ŷ_{1} (T) - y_{1} (T) ‖ \leq 3 N^{2.5} T h [{(‖ F_{2} ‖ + ‖ F_{1} ‖ + ‖ F_{0} ‖)}^{2} + ‖ F_{0}^{'} ‖] .

[11]

SI Appendix, Section 3.A.2, gives the proof of Lemma 2.

Next we upper bound the condition number of the linear system.

Lemma 3. Consider an instance of the quantum quadratic ODE problem as defined in Problem 1. Apply the forward Euler method [8] to the Carleman linearization [3]. Then the condition number of the linear system satisfies

κ \leq 3 (m + p + 1) .

[12]

SI Appendix, Section 3.B, gives the proof of Lemma 3.

We next describe a procedure for preparing the right-hand side of the linear system.

Lemma 4. Assume we are given the value

‖ u_{i n} ‖

, and let

O_{x}

be an oracle that maps

|00 \dots 0⟩ \in C^{n}

to a normalized quantum state

|u_{i n}⟩

proportional to

u_{i n}

. Assume we are also given the values

‖ F_{0} (t) ‖

for each

t \in [0, T]

, and let

O_{F_{0}}

be an oracle that provides the locations and values of the nonzero entries of

F_{0} (t)

for any specified

t

. Then a quantum state encoding the right-hand side of the linear system can be prepared using

O (N)

queries to

O_{x}

and

O (m)

queries to

O_{F_{0}}

, with gate complexity larger by a

poly (\log N, \log n)

factor.

The proof of Lemma 4 appears in SI Appendix, Section 3.C, along with a discussion of why the encoding of the initial state suffices.

After applying the QLSA to [8], we perform a measurement to extract a final state of the desired form. We now consider the probability that this measurement succeeds.

Lemma 5. Consider an instance of the quantum quadratic ODE problem defined in Problem 1, with the QLSA applied to the linear system using the forward Euler method [8]. Suppose the error from Carleman linearization satisfies

‖ η (t) ‖ \leq ‖ u (T) ‖ / 4

, and the error from Euler method satisfies

‖ ŷ (T) - y^{m} ‖ \leq ‖ u (T) ‖ / 4

. Then the probability of measuring a state

|y_{1}^{k}⟩

for

k = {[m + p + 1]}_{0} \ {[m + 1]}_{0}

satisfies

P_{m e a s u r e} \geq \frac{p + 1}{9 (m + p + 1) N q^{2}} .

[13]

We prove Lemma 5 in SI Appendix, Section 3.D.

The proof of the main algorithmic result (Theorem 1) follows by applying the above lemmas with appropriate parameters and carefully analyzing the complexity. The detailed analysis appears in SI Appendix, Section 3.E.

Lower Bound

In this section, we establish a limitation on the ability of quantum computers to solve the quadratic ODE problem when the nonlinearity is sufficiently strong. We quantify the strength of the nonlinearity in terms of the quantity

R

defined in [2]. Whereas there is an efficient quantum algorithm for

R < 1

(as shown in Theorem 1), we show here that the problem is intractable for

R \geq \sqrt{2}

. We establish this result by showing how the nonlinear dynamics can be used to distinguish nonorthogonal quantum states, a task that requires many copies of the given state.

Previous work on the computational power of nonlinear quantum mechanics shows that the ability to distinguish nonorthogonal states can be applied to solve unstructured search (and other hard computational problems) (16, 18, 19). Here we show a similar limitation using an information-theoretic argument.

Lemma 6. Let

| ψ 〉, | ϕ 〉

be states of a qubit with

| ⟨ ψ | ϕ ⟩ | = 1 - ϵ

. Suppose we are given either a black box that prepares

|ψ⟩

or a black box that prepares

|ϕ⟩

. Then any bounded-error protocol for determining whether the state is

|ψ⟩

or

|ϕ⟩

must use

Ω (1 / ϵ)

queries.

The proof of Lemma 6 is presented in SI Appendix, Section 4.A.

Lemma 6 can be used to establish limitations on the ability of quantum computers to simulate general nonlinear dynamics since nonlinear dynamics can be used to distinguish nonorthogonal states. Whereas previous work considers models of nonlinear quantum dynamics [such as the Weinberg model (18, 19) and the Gross-Pitaevskii equation (16)], here we aim to show the difficulty of efficiently simulating more general nonlinear ODEs—in particular, quadratic ODEs with dissipation—using quantum algorithms.

Lemma 7. There exists an instance of the quantum quadratic ODE problem as defined in Problem 1 with

R \geq \sqrt{2}

, and two states of a qubit with overlap

1 - ϵ

(for

0 < ϵ < 1 - 3 / \sqrt{10}

) as possible initial conditions, such that the two final states after evolution time

T = O (\log (1 / ϵ))

have an overlap no larger than

3 / \sqrt{10}

.

SI Appendix, Section 4.B, presents the proof of Lemma 7.

The proof of our main lower bound result (Theorem 2) follows by combining Lemma 6 and Lemma 7, and is presented in SI Appendix, Section 4.C.

Note that exponential time is achievable since our QCL algorithm can solve the problem by taking

N

to be exponential in

T

, where

N

is the truncation level of Carleman linearization. [The algorithm of Leyton and Osborne also solves quadratic differential equations with complexity exponential in

T

but requires the additional assumptions that the quadratic polynomial is measure-preserving and Lipschitz continuous (12).]

Applications

Our quantum algorithm could be applied to study models governed by inhomogeneous quadratic ODEs with linear dissipation. Such equations arise in biology and epidemiology as well as in fluid and plasma dynamics.

One concrete example is the so-called SEIR model of an epidemic (24). The model

\frac{d P_{S}}{d t} = - Λ \frac{P_{S}}{P} - r_{vac} P_{S} + Λ - r_{tra} P_{S} \frac{P_{I}}{P}

[14]

\frac{d P_{E}}{d t} = - Λ \frac{P_{E}}{P} - \frac{P_{E}}{T_{lat}} + r_{tra} P_{S} \frac{P_{I}}{P}

[15]

\frac{d P_{I}}{d t} = - Λ \frac{P_{I}}{P} + \frac{P_{E}}{T_{lat}} - \frac{P_{I}}{T_{inf}}

[16]

\frac{d P_{R}}{d t} = - Λ \frac{P_{R}}{P} + r_{vac} P_{S} + \frac{P_{I}}{T_{inf}}

[17]

describes the population of

P

individuals composed of four components: susceptible (

P_{S}

), exposed (

P_{E}

), infected (

P_{I}

), and recovered (

P_{R}

). Here

r_{tra}

is the rate of transmission from an infected to a susceptible person,

r_{vac}

is the vaccination rate,

T_{lat}

is the latent time until an exposed person becomes infectious,

T_{inf}

is the infectious time that an infectious person can infect others, and

Λ

is a flux of individuals constantly refreshing the population. This type of model has been used to describe the early spread of the COVID-19 virus (24). Realistic parameters from that study give a value

R < 1

(as discussed in more detail in SI Appendix, Section 5), showing that the assumptions of our algorithm can correspond to some real-world problems that are only moderately nonlinear. This example can be generalized to a high-dimensional system of ODEs that models the early spread over a large number of cities with interaction (25, 26).

We also consider the celebrated Navier–Stokes equation

\partial_{t} u + (u \cdot \nabla) u + β u = ν \nabla^{2} u + f

[18]

with linear damping and a forcing term. Equations of the form [18] are ubiquitous in fluid mechanics (27), and related models such as those studied in refs. 28–30 are used to describe the formation of large-scale structure in the universe. Similar equations also appear in magnetohydrodynamics (31) and in models that describe the motion of free particles that stick to each other upon collision (32). In the inviscid case,

ν = 0

, the resulting Euler-type equations with linear damping are also of interest, both for modeling micromechanical devices (33) and for their intimate connection with viscous models (34).

As a specific example, consider the one-dimensional forced viscous Burgers equation

\partial_{t} u + u \partial_{x} u = ν \partial_{x}^{2} u + f,

[19]

which is the one-dimensional case of Eq. 18 with

β = 0

. This is often used as a simple model of convective flow (35). For concreteness, let the initial condition be

u (x, 0) = U_{0} \sin (2 π x / L_{0})

on the domain

x \in [- L_{0} / 2, L_{0} / 2]

, and use Dirichlet boundary conditions

u (- L_{0} / 2,0) = u (L_{0} / 2,0) = 0

. We force this equation using a localized off-center Gaussian with a sinusoidal time dependence, given by

f (x, t) = U_{0} \exp (- \frac{{(x - L_{0} / 4)}^{2}}{2 {(L_{0} / 32)}^{2}}) \cos (2 π t)

. To solve this equation using the Carleman method, we discretize the spatial domain into

n_{x}

points and use central differences for the derivatives to get a system of quadratic ODEs. This equation is of the form [1] and can thus generate the Carleman system [4]. The resulting linear ODE can then be integrated using the forward Euler method, as shown in Fig. 1. Even though the parameters used in this example result in

R \approx 44

, which violates the requirement

R < 1

of the QCL algorithm, we see from the absolute error plot in Fig. 1 that the maximum absolute error over time decreases exponentially as the truncation level

N

is incremented (in this example, the maximum Carleman truncation level considered is

N = 4

). Surprisingly, this suggests that in this example, the error of the classical Carleman method converges exponentially with

N

, even though

R > 1

.

Fig. 1.

Integration of the forced viscous Burgers equation using Carleman linearization on a classical computer (source code available in ref. 38). The viscosity is set so that the Reynolds number $Re = U_{0} L_{0} / ν = 20$ . The parameters $n_{x} = 16$ and $n_{t} = 4,000$ are the number of spatial and temporal discretization intervals, respectively. The corresponding Carleman convergence parameter is $R = 43.59$ . (*Left*) Initial condition, source, and solutions plotted at a third of the nonlinear time $\frac{1}{3} T_{n l} = \frac{L_{0}}{3 U_{0}}$ . The direct solution is obtained from the ode45 solver in MATLAB. (*Center*) $l_{2}$ norm of the absolute error between the Carleman solutions at various truncation levels $N$ and the direct solution. (*Right*) Convergence of the time-maximum error.

Discussion

In this paper we have presented a QCL algorithm for a class of quadratic nonlinear differential equations. Compared to the previous approach of ref. 12, our algorithm improves the complexity from an exponential dependence on

T

to a nearly quadratic dependence, under the condition

R < 1

as defined in [2]. Qualitatively, this means that the system must be dissipative and that the nonlinear and inhomogeneous effects must be small relative to the linear effects. We have also provided numerical results suggesting the classical Carleman method may work on certain PDEs that do not strictly satisfy the assumption

R < 1

. Furthermore, we established a lower bound showing that for general quadratic differential equations with

R \geq \sqrt{2}

, quantum algorithms must have worst-case complexity exponential in

T

. We also discussed several potential applications arising in biology and fluid and plasma dynamics.

It is natural to ask whether the result of Theorem 1 can be achieved with a classical algorithm, i.e., whether the assumption

R < 1

makes differential equations classically tractable. Clearly, a naive integration of the truncated Carleman system [4] is not efficient on a classical computer since the system size is

Θ (n^{N})

. However, furthermore, it is unlikely that any classical algorithm for this problem can run in time polylogarithmic in

n

. Indeed, if we consider Problem 1 with

‖ λ_{1} ‖ T ≪ 1

but let the nonlinearity and forcing be even smaller such that

R < 1

, then in the asymptotic limit we have a linear differential equation with no dissipation. Hence, any classical algorithm that could solve Problem 1 could also solve nondissipative linear differential equations, which is a

BQP

-hard problem even when the dynamics are unitary (36). In other words, an efficient classical algorithm for this problem would imply efficient classical algorithms for any problem that can be solved efficiently by a quantum computer, which is considered unlikely.

Our upper and lower bounds leave a gap in the range

1 \leq R < \sqrt{2}

, for which we do not know the complexity of the quantum quadratic ODE problem. We hope that future work will close this gap and determine for which

R

the problem can be solved efficiently by quantum computers in the worst case.

Furthermore, the complexity of our algorithm has nearly quadratic dependence on

T

, namely,

T^{2} poly (\log T)

. It is unknown whether the complexity for quadratic ODEs must be at least linear or quadratic in

T

. Note that sublinear complexity is impossible in general because of the no-fast-forwarding theorem (37). However, it should be possible to fast-forward the dynamics in special cases, and it would be interesting to understand the extent to which dissipation enables this.

The complexity of our algorithm depends on the parameter

q

defined in Theorem 1, which characterizes the decay of the final solution relative to the initial condition. This restricts the utility of our result since we must have a suitable initial condition and terminal time such that the final state is not exponentially smaller than the initial state. However, it is unlikely that such a dependence can be significantly improved since renormalization of the state can be used to implement postselection, which would imply the unlikely consequence

B Q P = P P

(see section 8 of ref. 2 for further discussion). As discussed in the introduction, the solution of a homogeneous dissipative equation necessarily decays exponentially in time, so our method is not asymptotically efficient. However, for inhomogeneous equations the final state need not be exponentially smaller than the initial state even in a long-time simulation, suggesting that our algorithm could be especially suitable for models with forcing terms.

The quantum part of the algorithm might also be improved. In this paper we limit ourselves to the first-order Euler method to discretize the linearized ODEs in time. This is crucial for the analysis in Lemma 2, which states that the global error increases at most linearly with

T

(see the formal version of Lemma 2 in SI Appendix, Section 3.A.2). However, it is unclear how to give a similar bound for higher-order numerical schemes. If this obstacle could be overcome, the error dependence of the complexity might be improved.

It is also natural to ask whether our approach can be improved by taking features of particular systems into account. Since the Carleman method has only received limited attention and has generally been used for purposes other than numerical integration, it seems likely that such improvements are possible. In fact, the numerical results discussed in Applications (see in particular Fig. 1) suggest that the condition

R < 1

is not a strict requirement for the viscous Burgers equation, since we observe convergence even though

R \approx 44

. This raises the possibility that our approach could be applicable in a broader context than our analytical results might suggest. We leave a detailed investigation of this for future work.

A related question is whether our algorithm can efficiently simulate systems exhibiting dynamical chaos. The condition

R < 1

might preclude chaos, but we do not have a proof of this. More generally, the presence or absence of chaos might provide a more fine-grained picture of the algorithm’s efficiency.

When contemplating applications, it should be emphasized that our approach produces a state vector that encodes the solution without specifying how information is to be extracted from that state. Simply producing a state vector is not enough for an end-to-end application since the full quantum state cannot be read out efficiently. In some cases it may be possible to extract useful information by sampling a simple observable, whereas in other cases, more sophisticated postprocessing may be required to infer a desired property of the solution. Our method does not address this issue but can be considered as a subroutine whose output will be parsed by subsequent quantum algorithms. We hope that future work will address this issue and develop end-to-end applications of these methods.

Finally, the algorithm presented in this paper might be extended to solve related mathematical problems on quantum computers. Obvious candidates include initial value problems with time-dependent coefficients and boundary value problems. Carleman methods for such problems are explored in ref. 21, but it is not obvious how to implement those methods in a quantum algorithm. It is also possible that suitable formulations of problems in nonlinear optimization or control could be solvable using related techniques.

Data Availability

Source code data have been deposited in GitHub (38).

Acknowledgments

We thank Paola Cappellaro for valuable discussions and comments. J.-P.L. thanks Aaron Ostrander for inspiring discussions. H.Ø.K. thanks Bernhard Paus Graesdal for introducing him to the Carleman method. H.K.K. and N.F.L. thank Seth Lloyd for preliminary discussions on nonlinear equations and the quantum linear systems algorithm. We also thank Dong An and anonymous referees for pointing out the exponential decay of solutions to dissipative homogeneous equations. A.M.C. and J.-P.L. did part of this work while visiting the Simons Institute for the Theory of Computing in Berkeley and gratefully acknowledge its hospitality. A.M.C. and J.-P.L. acknowledge support from the Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Quantum Algorithms Teams, and Accelerated Research in Quantum Computing programs and from the NSF (CCF-1813814). H.Ø.K., H.K.K., and N.F.L. were partially funded by award DE-SC0020264 from the Department of Energy.

Supporting Information

Appendix (PDF)

Download
384.02 KB

References

1

D. W. Berry, High-order quantum algorithm for solving linear differential equations. J. Phys. A Math. Theor. 47, 105301 (2014).

Crossref

Google Scholar

2

D. W. Berry, A. M. Childs, A. Ostrander, G. Wang, Quantum algorithm for linear differential equations with exponentially improved dependence on precision. Commun. Math. Phys. 356, 1057–1081 (2017).

Crossref

Google Scholar

3

A. M. Childs, J. P. Liu, Quantum spectral methods for differential equations. Commun. Math. Phys. 375, 1427–1457 (2020).

Crossref

Google Scholar

4

B. D. Clader, B. C. Jacobs, C. R. Sprouse, Preconditioned quantum linear system algorithm. Phys. Rev. Lett. 110, 250504 (2013).

Crossref

PubMed

Google Scholar

5

Y. Cao, A. Papageorgiou, I. Petras, J. Traub, S. Kais, Quantum algorithm and circuit design solving the Poisson equation. New J. Phys. 15, 013021 (2013).

Crossref

Google Scholar

6

A. Montanaro, S. Pallister, Quantum algorithms and the finite element method. Phys. Rev. A (Coll. Park) 93, 032324 (2016).

Crossref

Google Scholar

7

P. C. S. Costa, S. Jordan, A. Ostrander, Quantum algorithm for simulating the wave equation. Phys. Rev. A (Coll. Park) 99, 012323 (2019).

Crossref

Google Scholar

8

A. M. Childs, J. P. Liu, A. Ostrander, High-precision quantum algorithms for partial differential equations. arXiv [Preprint] (2020) https://arxiv.org/abs/2002.07868 (Accessed 18 February 2020).

Google Scholar

9

A. Engel, G. Smith, S. E. Parker, Quantum algorithm for the Vlasov equation. Phys. Rev. A (Coll. Park) 100, 062315 (2019).

Crossref

Google Scholar

10

N. Linden, A. Montanaro, C. Shao, Quantum vs. classical algorithms for solving the heat equation. arXiv [Preprint] (2020). https://arxiv.org/abs/2004.06516 (Accessed 14 April 2020).

Google Scholar

11

A. W. Harrow, A. Hassidim, S. Lloyd, Quantum algorithm for linear systems of equations. Phys. Rev. Lett. 103, 150502 (2009).

Crossref

PubMed

Google Scholar

12

S. K. Leyton, T. J. Osborne, A quantum algorithm to solve nonlinear differential equations. arXiv [Preprint] (2008). https://arxiv.org/abs/0812.4423 (Accessed 5 March 2020).

Google Scholar

13

I. Joseph, Koopman-von Neumann approach to quantum simulation of nonlinear classical dynamics. arXiv [Preprint] (2020). https://arxiv.org/abs/2003.09980 (Accessed 22 March 2020).

Google Scholar

14

I. Y. Dodin, E. A. Startsev, On applications of quantum computing to plasma simulations. arXiv [Preprint] (2020). https://arxiv.org/abs/2005.14369 (Accessed 29 May 2020).

Google Scholar

15

S. Lloyd et al., Quantum algorithm for nonlinear differential equations. arXiv [Preprint] (2020). https://arxiv.org/abs/2011.06571 (Accessed 12 November 2020).

Google Scholar

16

A. M. Childs, J. Young, Optimal state discrimination and unstructured search in nonlinear quantum mechanics. Phys. Rev. A (Coll. Park) 93, 022314 (2016).

Crossref

Google Scholar

17

C. H. Bennett, E. Bernstein, G. Brassard, U. Vazirani, Strengths and weaknesses of quantum computing. SIAM J. Comput. 26, 1510–1523 (1997).

Crossref

Google Scholar

18

D. S. Abrams, S. Lloyd, Nonlinear quantum mechanics implies polynomial-time solution for NP-complete and #P problems. Phys. Rev. Lett. 81, 3992 (1998).

Crossref

Google Scholar

19

S. Aaronson, NP-complete problems and physical reality. ACM SIGACT News 36, 30–52 (2005).

Crossref

Google Scholar

20

T. Carleman, Application de la théorie des équations intégrales linéaires aux systèmes d’équations différentielles non linéaires. Acta Math. 59, 63–87 (1932).

Crossref

Google Scholar

21

K. Kowalski, W. H. Steeb, Nonlinear Dynamical Systems and Carleman Linearization (World Scientific, 1991).

Crossref

Google Scholar

22

M. Forets, A. Pouly, Explicit error bounds for Carleman linearization. arXiv [Preprint] (2017). https://arxiv.org/abs/1711.02552 (Accessed 5 March 2020).

Google Scholar

23

A. M. Childs, R. Kothari, R. D. Somma, Quantum algorithm for systems of linear equations with exponentially improved dependence on precision. SIAM J. Comput. 46, 1920–1950 (2017).

Crossref

Google Scholar

24

C. Wang et al., Evolving epidemiology and impact of non-pharmaceutical interventions on the outbreak of coronavirus disease 2019 in Wuhan, China. JAMA 323, 1915–1923 (2020).

PubMed

Google Scholar

25

D. Bichara, Y. Kang, C. Castillo-Chávez, R. Horan, C. Perrings, SIS and SIR epidemic models under virtual dispersal. Bull. Math. Biol. 77, 2004–2034 (2015).

Crossref

PubMed

Google Scholar

26

R. Qurratu Aini, D. Aldila, K. Sugeng, Basic reproduction number of a multi-patch SVI model represented as a star graph topology. AIP Conf. Proc. 2023, 020237 (2018).

Crossref

Google Scholar

27

P. G. Lemarié-Rieusset, The Navier-Stokes Problem in the 21st Century (CRC Press, 2018).

Crossref

Google Scholar

28

L. Kofman, D. Pogosian, S. Shandarin, Structure of the universe in the two-dimensional model of adhesion. Mon. Not. R. Astron. Soc. 242, 200–208 (1990).

Crossref

Google Scholar

29

S. F. Shandarin, Y. B. Zeldovich, The large-scale structure of the universe: Turbulence, intermittency, structures in a self-gravitating medium. Rev. Mod. Phys. 61, 185–220 (1989).

Crossref

Google Scholar

30

M. Vergassola, B. Dubrulle, U. Frisch, A. Noullez, Burgers’ equation, devil’s staircases and the mass distribution for large-scale structures. Astron. Astrophys. 289, 325–356 (1994).

Google Scholar

31

P. A. Davidson, An Introduction to Magnetohydrodynamics (Cambridge Texts in Applied Mathematics) (Cambridge University Press, 2001).

Crossref

Google Scholar

32

Y. Brenier, E. Grenier, Sticky particles and scalar conservation laws. SIAM J. Numer. Anal. 35, 2317–2328 (1998).

Crossref

Google Scholar

33

M. H. Bao, Micro Mechanical Transducers: Pressure Sensors, Accelerometers and Gyroscopes (Elsevier, 2000).

Google Scholar

34

C. M. Dafermos, Hyperbolic Conservation Laws in Continuum Physics (Springer, 2005), vol. 3.

Crossref

Google Scholar

35

J.M. Burgers, A mathematical model illustrating the theory of turbulence. Adv. Appl. Mech. 1, 171–199 (1948).

Crossref

Google Scholar

36

R. P. Feynman, Quantum mechanical computers. Opt. News 11, 11–20 (1985).

Crossref

Google Scholar

37

D. W. Berry, G. Ahokas, R. Cleve, B. C. Sanders, Efficient quantum algorithms for simulating sparse Hamiltonians. Commun. Math. Phys. 270, 359–371 (2007).

Crossref

Google Scholar

38

H. Ø. Kolden, Carleman solution of the viscous Burgers equation. GitHub. https://github.com/hermankolden/CarlemanBurgers. Deposited 1 March 2021.

Google Scholar

Information & Authors

Information

Published in

Proceedings of the National Academy of Sciences

Vol. 118 | No. 35
August 31, 2021

PubMed: 34446548

Classifications

Copyright

Data Availability

Source code data have been deposited in GitHub (38).

Submission history

Published online: August 26, 2021

Published in issue: August 31, 2021

Keywords

Acknowledgments

We thank Paola Cappellaro for valuable discussions and comments. J.-P.L. thanks Aaron Ostrander for inspiring discussions. H.Ø.K. thanks Bernhard Paus Graesdal for introducing him to the Carleman method. H.K.K. and N.F.L. thank Seth Lloyd for preliminary discussions on nonlinear equations and the quantum linear systems algorithm. We also thank Dong An and anonymous referees for pointing out the exponential decay of solutions to dissipative homogeneous equations. A.M.C. and J.-P.L. did part of this work while visiting the Simons Institute for the Theory of Computing in Berkeley and gratefully acknowledge its hospitality. A.M.C. and J.-P.L. acknowledge support from the Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Quantum Algorithms Teams, and Accelerated Research in Quantum Computing programs and from the NSF (CCF-1813814). H.Ø.K., H.K.K., and N.F.L. were partially funded by award DE-SC0020264 from the Department of Energy.

Notes

This article is a PNAS Direct Submission.

Authors

Affiliations

Jin-Peng Liu https://orcid.org/0000-0002-1885-8879

Joint Center for Quantum Information and Computer Science, University of Maryland, College Park, MD 20742;

Institute for Advanced Computer Studies, University of Maryland, College Park, MD 20742;

Department of Mathematics, University of Maryland, College Park, MD 20742;

View all articles by this author

Herman Øie Kolden https://orcid.org/0000-0003-4121-389X

Department of Physics, Norwegian University of Science and Technology, NO-7491 Trondheim, Norway;

Plasma Science and Fusion Center, Massachusetts Institute of Technology, Cambridge, MA 02139;

View all articles by this author

Hari K. Krovi

Quantum Engineering and Computing, Raytheon BBN Technologies, Cambridge, MA 02138;

View all articles by this author

Nuno F. Loureiro

Plasma Science and Fusion Center, Massachusetts Institute of Technology, Cambridge, MA 02139;

View all articles by this author

Konstantina Trivisa

Department of Mathematics, University of Maryland, College Park, MD 20742;

Institute for Physical Science and Technology, University of Maryland, College Park, MD 20742;

View all articles by this author

Andrew M. Childs¹ https://orcid.org/0000-0002-9903-837X [email protected]

Joint Center for Quantum Information and Computer Science, University of Maryland, College Park, MD 20742;

Institute for Advanced Computer Studies, University of Maryland, College Park, MD 20742;

Department of Computer Science, University of Maryland, College Park, MD 20742

View all articles by this author

Notes

1

To whom correspondence may be addressed. Email: [email protected].

Author contributions: J.-P.L., H.Ø.K., H.K.K., N.F.L., K.T., and A.M.C. designed research; J.-P.L., H.Ø.K., H.K.K., N.F.L., K.T., and A.M.C. performed research; J.P.L. led the theoretical analysis; H.Ø.K. led the numerical experiments; and J.-P.L., H.Ø.K., H.K.K., N.F.L., K.T., and A.M.C. wrote the paper.

Competing Interests

The authors declare no competing interest.

Metrics & Citations

Metrics

Note: The article usage is presented with a three- to four-day delay and will update daily once available. Due to ths delay, usage data will not appear immediately following publication. Citation information is sourced from Crossref Cited-by service.

Altmetrics

Citations

Cite this article

J. Liu,
H.Ø. Kolden,
H.K. Krovi,
N.F. Loureiro,
K. Trivisa,
& A.M. Childs,

Efficient quantum algorithm for dissipative nonlinear differential equations, Proc. Natl. Acad. Sci. U.S.A. 118 (35) e2026805118, https://doi.org/10.1073/pnas.2026805118 (2021).

Copied!

Copying failed.

Export the article citation data by selecting a format from the list below and clicking Export.

Cited by

View Options

View options

PDF format

Download this article as a PDF file

DOWNLOAD PDF

Login options

Check if you have access through your login credentials or your institution to get full access on this article.

Personal login Institutional Login

Recommend to a librarian

Recommend PNAS to a Librarian

Save for later

Purchase options

Purchase this article to access the full text.

Single Article Purchase

Efficient quantum algorithm for dissipative nonlinear differential equations

Proceedings of the National Academy of Sciences

Vol. 118
No. 35

Restore content access

Restore content access for purchases made as a guest

Featured Topics

Articles By Topic

Featured Topics

Articles By Topic

Featured Topic

Articles By Topic

Significance

Abstract

Sign up for PNAS alerts.

Problem Statement

Main Results

Algorithm Analysis

Lower Bound

Applications

Discussion

Data Availability

Acknowledgments

Supporting Information

References

Information

Published in

Classifications

Copyright

Data Availability

Submission history

Keywords

Acknowledgments

Notes

Authors

Affiliations

Notes

Competing Interests

Metrics

Altmetrics

Citations

Cite this article

Cited by

View options

PDF format

Login options

Recommend to a librarian

Purchase options

Restore content access

Share

Share article link

Share on social media

Further reading in this issue

Developmental HCN channelopathy results in decreased neural progenitor proliferation and microcephaly in mice

Interaction between decision-making and interoceptive representations of bodily arousal in frontal cortex

Replication across space and time must be weak in the social and environmental sciences

X and Y gene dosage effects are primary contributors to human sexual dimorphism: The case of height

Increased flood exposure in the Pacific Northwest following earthquake-driven subsidence and sea-level rise

Multiple sclerosis and gut microbiota: Lachnospiraceae from the ileum of MS twins trigger MS-like disease in germfree transgenic mice—An unbiased functional study

Sign up for thePNAS Highlights newsletter