Improved quantum algorithms for linear and nonlinear differential equations

We present substantially generalized and improved quantum algorithms over prior work for inhomogeneous linear and nonlinear ordinary differential equations (ODE). Specifically, we show how the norm of the matrix exponential characterizes the run time of quantum algorithms for linear ODEs opening the door to an application to a wider class of linear and nonlinear ODEs. In Berry et al., (2017), a quantum algorithm for a certain class of linear ODEs is given, where the matrix involved needs to be diagonalizable. The quantum algorithm for linear ODEs presented here extends to many classes of non-diagonalizable matrices. The algorithm here is also exponentially faster than the bounds derived in Berry et al., (2017) for certain classes of diagonalizable matrices. Our linear ODE algorithm is then applied to nonlinear differential equations using Carleman linearization (an approach taken recently by us in Liu et al., (2021)). The improvement over that result is two-fold. First, we obtain an exponentially better dependence on error. This kind of logarithmic dependence on error has also been achieved by Xue et al., (2021), but only for homogeneous nonlinear equations. Second, the present algorithm can handle any sparse, invertible matrix (that models dissipation) if it has a negative log-norm (including non-diagonalizable matrices), whereas Liu et al., (2021) and Xue et al., (2021) additionally require normality.


Introduction
Differential equations lie at the heart of many important problems in several fields including fluid mechanics, plasma physics and quantum physics. Indeed, simulating the Schrödinger equation, a central equation in quantum physics, was one of the first motivations to build a quantum computer and it is a linear homogeneous differential equation. Simulating the Schrödinger equation or, in other words, Hamiltonian simulation can be done efficiently for several classes of Hamiltonians and important techniques were developed in a long line of work [4][5][6][7][8]. The quantum algorithmic techniques developed in those works extend beyond Hamiltonian simulation. Some of those techniques such as block encoded matrix operations will be used in this paper. The technique of block encoded matrix operations is further developed in [9,10] by extending the block encoding to functions of a sparse matrix.
Another important development in quantum algorithms is the algorithm to solve linear systems of equations [11]. This algorithm has been applied to numerous problems such as in machine learning and differential equations [1,2,[12][13][14][15]. The run-time of the original algorithm [11] has an exponentially better dependence on dimension and it scales quadratically as function of the condition number κ of the linear system. This dependence on the condition number has been improved in subsequent work to κ log κ [16]. In [17], the dependence on the solution error has been exponentially improved and made logarithmic. In a recent series of papers, the scaling has been improved to κ log 1/ [18][19][20][21]. The Quantum Linear Systems Algorithm (QLSA) with these improvements will be used here as well.
Quantum algorithms for differential equations have been designed in several papers starting with the early work on nonlinear differential equations in [22]. In [12], a quantum algorithm to solve linear inhomogeneous equations is presented that uses high-order methods. In [1], using truncated Taylor series, an exponential improvement in the dependence of the solution error was obtained. In [14], a quantum algorithm to solve time dependent linear differential equations was presented using spectral methods. These spectral methods are also applied to develop a quantum algorithm to solve linear partial differential equations in [13].
Recently, in [2], a quantum algorithm to solve dissipative nonlinear differential equations is given using the so called Carleman linearization. This algorithm is efficient when the ratio of nonlinearity to dissipation is less than 1. In [15], a quantum algorithm for nonlinear differential equations is given when the time scale of simulation is short. However, this was heuristically derived without a rigorous analysis. Recently, using homotopic perturbation methods, another quantum algorithm for nonlinear differential equations is given in [3], where the dependence on the error is exponentially better than in [2]. However, the ratio of nonlinearity to dissipation is required to be smaller than in [2]. The algorithm of [3] is also only for homogeneous differential equations with the assumption of normality. In the regime of nonlinearity considered in [3], homogeneous differential equations have an exponential decay of solutions leading to a run-time that is exponential in T (the simulation time).
Other works that address aspects of differential equations but with an application to plasma physics are [23], where a quantum algorithm to solve a linearized Vlasov equation is given. In [24,25], a formulation of several interesting problems in plasma physics that could be amenable to quantum algorithms through linearization is given and in [26], a quantum algorithm to solve classical dynamics using the Koopman operator approach is given. In [27], a quantum algorithm to solve cold plasma waves is given. In [28], a connection between Wigner-Weisskopf theory of spontaneous decay and non-Hermitian evolution is made and an implementation using "non-Hermitian" quantum circuits is explored.
In this paper, we first present a quantum algorithm to solve time-independent linear inhomogeneous equations that improves on the results of [1] in two ways. In order to do this, we extend their algorithm to include many classes of non-diagonalizable and even singular matrices that are not covered by their algorithm. As was pointed out in [1], every matrix is close to a diagonalizable matrix (since the latter are dense in the set of all matrices). But using the closest diagonalizable matrix leads to an exponentially worse error. Here, we are able to extend to many classes of non-diagonalizable matrices keeping the logarithmic dependence on error. We give a characterization of matrices for which the algorithm is efficient by bounding the gate and query complexity. Diagonalization was needed in [1] in deriving a bound on the condition number and a bound on the solution error. To circumvent the need for diagonalization, we provide a different and improved analyses for both these parts.
The second main difference between the two algorithms is that even when we restrict to diagonalizable matrices there are cases when the bound derived for the algorithm of [1] is exponentially worse than ours. Specifically, for diagonalizable matrices that have non-positive log-norm but an exponentially high condition number κ V , where V is the matrix of eigenvectors, the gate complexity bound for the algorithm of [1] is exponential since it has at least a linear dependence on κ V . In our algorithm, this dependence is removed leading to an exponential improvement in this case. More generally, for all diagonalizable matrices with bounded norm of the matrix exponential (i.e., exp(At) is bounded) but with a high κ V , our algorithm is exponentially better than the bound in [1]. Numerical simulations for a specific example of the actual condition number of [1] show that the bounds derived there might be very loose for such matrices.
Despite these differences, the algorithm presented here has many common elements with [1]. Like that work, we use truncated Taylor series and construct a quantum linear system whose solution gives a quantum state proportional to the solution of the linear differential equation. We also make use of the "ramp" to boost the success probability using the technique introduced in [12] and used in subsequent papers. However, the linear system is different in our case and we implement it using the techniques of block encoding. The reason we pick a different linear system is that it is easier to analyze using the techniques developed here.
We then apply our linear ODE algorithm to solve nonlinear differential equations by using Carleman linearization. We improve on our prior work [2] by exponentially improving the dependence on the error. The main reason that [2] has a polynomial dependence on the error is the use of Euler method. The linear ODE obtained by Carleman linearization need not be diagonalizable and hence [1] is not immediately applicable. But the Carleman ODE has negative log-norm if the dissipation matrix of the nonlinear system also has that property. This makes it possible to use our present algorithm to achieve an exponential improvement in the error even though the matrix may not be diagonalizable. This also means that our linear ODE algorithm applied to Carleman linearized ODEs can handle non-normal and non-diagonalizable matrices. In principle, we should be able to generalize this further since all we need for our linear ODE algorithm is a general property that the norm of the matrix exponential is bounded. However, it seems hard to characterize dissipation matrices of nonlinear equations that lead to Carleman ODEs with this general property.
Since the first version of this paper appeared, there have been some more results some of which use the characterization of run time in terms of matrix exponential introduced here. In [29], the authors apply ideas from [2] to the case of reaction-diffusion equations which form a special class of nonlinear differential equations. In [30], the authors design a quantum algorithm for homogeneous time-dependent ODEs where the run-time bound depends on the matrix exponential. Recently, in [31], the authors extend to time-dependent inhomogeneous linear ODEs, the characterization that we have derived here in terms of the matrix exponential and in [32], some limitations of quantum algorithms for ODEs is explored. This paper is organized as follows. In Section 2, we precisely state the problem that our algorithm solves as well as some background on notation and results on block encoded matrices that we use later on in the paper. In Section 3, we discuss the behaviour of the norm of the matrix exponential, its dependence on quantities such as spectral abscissa and log-norm. We then give a characterization of matrices (even non-diagonalizable ones) that are well-behaved for our algorithm. Then in Section 4, we present the steps of the algorithm. In Section 5, we present the analysis of the algorithm by deriving bounds on the solution error, condition number of the linear system, the probability of success and the implementation of the algorithm. In Section 6, we prove the correctness of the algorithm and bound its gate and query complexity. In Section 7, we apply the linear ODE algorithm to nonlinear differential equations and show an exponential improvement in error. We also show that more general classes of nonlinear differential equations can be solved efficiently than the ones considered earlier. Finally, in Section 8, we present some conclusions and open questions.

Problem statement and preliminaries
In this section, we set notation and give some background and results from the literature on block encoding. We also precisely define the input model and the problem our algorithm solves. First, we set notation and define the quantities we use in this paper. We also give results from prior work that we use later on in the paper. In this paper, · denotes l 2 norm of a vector or a matrix. For an arbitrary matrix A, an eigenvalue is a root of the characteristic polynomial and singular value is the square root of an eigenvalue of A † A. We use the standard notation g(n) = O(f (n)) for a function f (n) to mean g(n) ≤ cf (n) for some c independent of n. To distinguish this O from the O used for oracles, we always use a subscript for the oracles (defined below).
The specific problem we consider is to obtain the solution of a linear ordinary differential equation as a quantum state. We state it below more precisely. Problem 1. Assume we are given a stable, sparse d × d matrix A, consider the linear ordinary differential equation The input model we consider is the same as the one considered in prior work on differential equations [1,2,[12][13][14]. The matrix A is assumed to be sparse with at most s r nonzero entries in any row and at most s c nonzero entries in any column. Definition 1. We assume that A ∈ C 2 n ×2 n can be accessed through the following oracles.

4)
where r ij (respectively c ji ) is the j th nonzero entry in the i th row (resp. column). If there are fewer than j entries, then r ij (resp. c ji ) is set to j +2 n . Here a ij is a b-bit binary representation of the (i, j) matrix entry of A. Similarly, we have oracles O b and O x to prepare a normalized version of b and x 0 . We assume that their norms are known. Specifically, let O x be any unitary that maps |1 |φ to |1 |φ for any state |φ and |0 |0 to |0 |x 0 , wherex 0 = x 0 / x 0 . Let O b be a unitary that similarly maps |1 |φ to |1 |φ for any state |φ and |0 |0 to |0 |b , whereb = b/ b .
Next, we give here some results on block encoded matrices that will be useful for us. Our main tool to construct a quantum algorithm is to use matrix arithmetic to construct the linear system corresponding to a given ODE. The linear system can be broken up into pieces involving subtraction, multiplication and inversion of matrices and we give a list of results from the literature that bound the complexity of these operations.
First, we define block encoding. The idea behind block encoding is to embed an arbitrary sparse matrix into a unitary matrix so that one can use embedding to perform arithmetic operations (see [10] for more details). Definition 2. Suppose A is an a-qubit operator, α, ∈ R + and b ∈ N, then an a + b qubit unitary U is said to be an (α, b, ) block encoding of A if (2.5)

Definition 3. We denote byĀ the Hermitian complement of A defined below
Next, we state the relevant results from the literature that we use in our algorithm. First, we need the following lemma that describes a way to block-encode a sparse matrix. Next, we need the following result from [8] to invert a block encoded matrix. Lemma 2 ([8], Lemma 9). Let A be a matrix with condition number κ ≥ 2. Let H be the Hermitian complement of A and suppose that I/κ ≤ H ≤ I. Let

AB has an
Next, we need the following theorem to implement the QLSA. Theorem 1 ([21], Theorem 19). Let A be such that A = 1 and A −1 = κ. Given a oracle block encoding of A and an oracle for implementing |b , there exists a quantum algorithm which produces the normalized state A −1 |b to within an error using calls to the oracles. Remark 2. In using Lemma 2 and Theorem 1, we will apply them to matrices with norms that are O(1) rather than 1. This does not change the asymptotic scaling of the gate complexities.

Norm of the matrix exponential
In the analysis of our algorithm, we will need to bound the norm of the matrix exponential e At for some matrix A and t ≥ 0. We explain the known results that bound this quantity and also the relation to stability theory of differential equations. For an arbitrary matrix A, in addition to the norm of the matrix, the following quantities are useful to bound the norm of functions of A. Definition 4.
We drop the dependence on A in the above quantities when it is clear from context. The quantity µ(A) is not a norm (even though it is called the log-norm) and can, in fact, be negative. The norm of e At depends on µ and α in the short and long timescales respectively. If we want the norm of e At to be bounded for all t, we need that α(A) < 0. But we will see below that this alone is not enough and we have to make additional assumptions depending on the sign of µ. The quantities α, µ and ρ are all real numbers. For arbitrary matrices, these quantities are related as follows.
In the theory of differential equations [33], a stable matrix is defined as a matrix whose eigenvalues (i.e., roots of the characteristic polynomial) have negative real parts. It is called semi-stable if the eigenvalues are allowed to have zero real parts. The Lyapunov condition, can be used to determine if a matrix is stable. This can be stated as follows. If there exist a positive definite matrix P and a negative definite matrix N such that AP + P A † = N , then A is stable. The Routh-Hurwitz criterion [33] can also be used to check for stability. In this section, we are interested in bounds on the norm of the matrix exponential for stable matrices A. Next, we summarize known results on these bounds. For an arbitrary matrix A, the exponential of At for t ≥ 0 converges to zero in the limit of long time if and only if A is a stable matrix [34], i.e., lim t→∞ exp(At) = 0 ⇐⇒ A is stable i.e., α < 0 . (3.6) The following bounds are well-known [35].
While this is a bound for all time, it turns out that the quantity exp(At) grows (or decays) like the function exp(µt) initially following it closely for very short times. In fact, another way to define the log norm is as the rate of change in exp(At) near t = 0. This behavior can be seen in the following simple example [36]. Consider the following two 2 × 2 matrices.
The plot in Figure 1 shows the difference in short term behavior of the norm of the matrix exponential depending on the sign of µ. We would like to point out that restricting attention to upper triangular matrices is without loss of generality since by the Schur decomposition [37], any matrix is unitarily equivalent to an upper triangular matrix and the l 2 norm is unitarily invariant. Based on the above example, we can consider the two cases depending on the log-norm: when µ ≤ 0 and µ > 0.

Non-positive log-norm
In this case, the matrix exponential is easily bounded using the following inequality.
This bound in useful when the log-norm is known to be negative. For instance, consider the following tridiagonal matrix.
Matrices of this type are known as twisted Toeplitz or Berezin-Toeplitz matrices. Matrices like these appear in certain plasma physics models [38]. Matrices with a negative log-norm appear in high energy physics e.g., in modeling particle decay [28,39]. The above matrix is diagonalizable but is not a normal matrix. Therefore, the diagonalizing matrix V is not unitary. In this case, V has a condition number κ V that is exponential in log d (see Figure 2a). However, it can be seen that (A + A † )/2 is a diagonal matrix with negative entries and therefore µ < 0. In fact, it can be seen that µ = −1/d. This also means that α < 0 since α < µ. If κ V grows exponentially as in this case (see Figure 2a), the bound   is exponentially worse than the bound (3.9).
Using this example, we check here numerically if the condition number of the linear system from [1] is indeed exponential or if the bound is just loose. The matrix C in the linear system from [1] is defined as We compare the actual condition number of the linear system C in [1] to the bound derived in [1] as well as the condition number of the linear system L defined here. In Figure 2a, we plot the relevant quantities for the matrix in (3.10) for dimensions ranging from 15 to 100. We find that the condition number κ C of the linear system C is far smaller than the bound derived in [1]. While the bound, which depends on the condition number κ V of the diagonalizing matrix V , turns out to be exponential, κ C seems to scale polynomially for this example. The condition number κ L of the linear system L that we define here is lower than κ C and the gap between them grows with dimension as shown in Figure 2b. It is possible that for the matrix in (3.10) as well as other matrices A, the condition number κ C scales polynomially if the norm of the exponential of A is bounded. If this is true, then it seems likely that the techniques developed in the present work can be used to prove this.

Positive log-norm
For the case of positive log-norm, it turns out that exp(At) can increase initially since the log norm is positive. If we assume that α < 0, the long time behavior is exponential decay. The initial increase in norm is bounded and the Kreiss matrix theorem gives a bound on its supremum.
The Kreiss matrix theorem states the following.
Theorem 2 (Kreiss matrix theorem [36]). For an arbitrary d×d matrix A, the following bounds are optimal.
where e is the Euler constant, d is the dimension of A and the Kreiss constant K(A) is defined as While this bound is useful for arbitrary matrices, the fact that, as a uniform bound, it is optimal and can be as large as the dimension of A shows that we need to restrict the class of matrices to focus on.
To restrict the class of matrices, we give a couple of bounds on exp(At) . These bounds and others based on pre-conditioning are discussed in [40]. The two main ways to derive bounds on the matrix exponential are to use the Jordan and Schur decompositions [37]. The following bound using Jordan blocks was shown in [34].
Note that when all the Jordan blocks are one dimensional i.e., the diagonalizable case, this reduces to the bound from [1]. The next bound was also shown in [34] and it uses the Schur decomposition.

16)
where d is the dimension of A and  [41]. The reason we present several bounds here is that no one bound completely subsumes the others, especially when the log-norm is positive. When µ(A) > 0, it is possible that the bound of (3.9) is better than the ones in Lemma 4 and Lemma 5 for some t. This is because for the latter bounds, the exponential exp(αt) needs to overcome the constants κ V and polynomials p d−1 ( N t) or t β before falling below a reasonable value, whereas for small positive µ and t, the former could be small.
Finally, we define here a quantity that we use in our algorithm that may be satisfied by matrices in the above classes e.g., in the remark above and more. The reason we are interested in this definition is that it encompasses even non-stable matrices. For instance, if T A is a constant (independent of n = log d), then C(A) ≤ exp(T A ) is also acceptable as a bound where A can be an arbitrary matrix. To verify that C(A) is bounded is not always an easy task. This is not uncommon in quantum or classical algorithms. For instance, verifying that the condition number is bounded for a linear system is not generally easy or checking if a matrix is diagonalizable as in the some prior algorithms. However, in some cases as in the example given above in (3.10), it is easy to verify that the matrix has a non-positive log norm. In general, when A is a block matrix such as in the application for non-linear differential equations, one can bound C(A) without additional assumptions on the component matrices such as e.g., normality or diagonalizability.

Description of the quantum algorithm
We assume that we have a time-independent linear system of the form.
where A and b are time independent. The solution x(T ) at time T of this matrix differential equation is which is valid even when A is a singular matrix [42]. We can re-write it in a form as follows that will be useful later. Proof. Expanding the exponential, we have for the second term This can be seen to be equal to (4.2).
To design a quantum algorithm to find the solution at time T , we use the method from [1] of truncated Taylor series. Using their notation, we get the following discrete equation at time step i.

5)
where h is the time interval of a single step and T k (z) and S k (z) are the following functions.
Here we have m = T /h time steps in the time interval [0, T ]. We will encode this equation into a linear system and use the quantum linear system algorithm (QLSA) to solve it. Our implementation also approximates via a truncated Taylor series similar to [1], but we use a different linear operator that we find easier to analyze to provide a better analysis and generalization.
To solve the linear inhomogeneous equation, we use the linear operator L, where Here the first register, indexed by i, in the sum contains the time step. The second register, indexed by j, is the Taylor step and the third register contains states on which A acts. The quantity p is the number of additional steps needed to boost the success probability (this has become a standard trick that was first used in [12] and subsequently in several other papers [2,3,14]). We give an example of M 1 and M 2 for k = 3 below. (4.12) The above operator L is used to solve the following linear system. where is the unnormalized version of the initial state. We now show that implementing L −1 on the initial state gives us an entangled state from which we can extract the solution. Applying L −1 to the initial state gives since N is nilpotent and N m+p+1 = 0. To examine the above state, let us consider the action of (I − M 1 ) −1 on an arbitrary basis state of the form | , v . We have Note that T 0,k (Ah) = T k (Ah). Therefore, the action of this on the initial state is where y i are given by (4.5). This becomes the initial state for the next power of N giving us The overall state is then Now if one measures the time-step register and post-selects on the outcomes m, . . . , m + p − 1, the state collapses to a state |ψ which is y m / y m . We show in subsequent sections that this state is close to the normalized solution, that the measurement probability is significant, that the condition number is well-behaved and give bounds on the query and gate complexities of implementing the algorithm.

Analysis of the quantum algorithm
In this section, we first derive bounds on the solution error and then bound the condition number and probability of success of the algorithm. For the solution error, we are interested in deriving a bound on the number of terms to keep for the Taylor series for the error between the actual and approximate solutions of the ODE to be small. In [1], the same thing was done assuming that the matrix is diagonalizable. The results here are valid for an arbitrary sparse matrix which can be non-diagonalizable and even singular. Standard books such as [43] do not seem to have the specific results we need.

Solution error
Before we prove our bounds, let us define the following quantities. From Lemma 6, the solution after evolving for time T can be written as

8)
After m time steps, the solution is We need the following bounds later. Lemma 7. When Ah ≤ 1, we have Proof. In terms of the remainder, we have since Ah ≤ 1.
We have a similar bound on l 1 and l 1 . Lemma 8. When Ah ≤ 1, we have Proof. We have Lemma 9. When Ah ≤ 1 and me 2 /(k + 1)! ≤ 1, Proof. Since l 0 and l 0 commute, we have Using Lemma 7 and the triangle inequality, we have Denote r = 1 + e 2 (k+1)! . Using this, we get where the last line follows by using the inequality since the r.h.s is less than unity. This proves the lemma.

Remark 5.
The proof of the above lemma goes through for any j ≤ m i.e., we have when Ah ≤ 1 and me 2 /(k + 1)! ≤ 1.
We also need the following lemma. Finally, this can be written as Now we have Using the above equations as well as Lemma 9 and Lemma 8 and plugging into (5.36), we get where we used the fact that m = T A + 1.

Lemma 11.
We have the following relations.
Proof. Consider Al 1 . We have We can similarly show that We can now prove a bound on the solution error. Theorem 3. Suppose x T is the solution at time T and y m is the solution given by the truncated Taylor series at k terms, where and Ah ≤ 1, then we have the bound Proof. From Lemma 6, the solution after evolving for time T is Now the final state after m time steps of the truncated Taylor series is Now, let x h = L 0 x 0 be the homogeneous part of the solution, then the error is given by where in the second line, we used Lemma 11, the second-to-last line follows from the triangle inequality and the last line from Lemma 9 and Lemma 10. Now, using the value of k defined in the statement of this theorem, we have

Condition number
In this subsection, we prove a theorem that gives a bound on the condition number.

Lemma 12. Recall that the error in the truncated Taylor series T k (At) of exp(At) is
Then for Ah ≤ 1, the norm of this error satisfies Proof. We can write the error norm as Next we prove a lemma on the norm of powers of the truncated Taylor series. Lemma 13. For Ah ≤ 1 and k as in Theorem 3, we have that the truncated Taylor series Proof. The condition number of L is κ L = L L −1 . We bound each of these below. The norm of L is To bound the norm of M 2 (I − M 1 ) −1 , we can use its action on a state of the form | , v for some normalized basis state |v . From (4.18), this is We can see that for all 0 ≤ ≤ k (since Ah ≤ 1), we have the following bound.
Proof. First, let |y g be defined as follows.
The probability of a successful measurement is where y is the state from (4.22). The squared norm of y is Let i 0 be the index that maximizes y i 2 . Then, using Theorem 3, we have x(t) .
(5.86) Therefore, we have Since the solution error is x T − y m ≤ δ x T from Theorem 3, we have Choosing m = p and δ ≤ 1/2, we get P meas ≥ 1 18g 2 . (5.92)

Initial state preparation and circuit implementation
The preparation scheme for the initial state is very similar to the one in [1]. Proof. Recall that the initial state is To prepare this state from |0, 0, 0 , first apply the following rotation on the second register.
This step takes O(1) gates since it is a unitary on a single qubit. Next, we apply the oracles O x and O b to the third register conditioned on the second register being in |0 and |1 respectively to get Finally, conditioned on the second register being |1 , we implement the following rotation on the first register that takes This step takes polylog(m) gates. All this gives us the normalized initial state we want with the complexity stated in the theorem.
We now compute the query and gate complexity of implementing the QLSA on the linear system (4.8). 2. Using Lemma 2, the inverse of (I − M 1 ) −1 can be implemented as an (α 1 , a 1 , 1 ) block unitary with complexity where s = max{s r , s c }, α 1 = 2k and a 1 = log(dk) + O(log(2k/ 1 )). The fact that α 1 = 2k comes from the condition number of I − M 1 .

Using
where T ψ is the complexity of preparing the initial state.
For the total error of the implementation to be bounded by , we first choose 1 = 2 . We then need For this to be satisfied, we need Plugging this into expressions for α 3 , T 4 and a, we get  and gate complexity which is greater by a factor of at most where from (5.91), g is defined as Proof. We pick the parameters of the algorithm as follows. where This choice of k makes (k + 1)! > Ω. Let the output of the algorithm described in Section 4 be the state |ψ . In Section 4, we showed that the quantum state |ψ produced by the algorithm is y m / y m . In Theorem 3, it was shown that This means, we have From Theorem 5, the success probability is proportional to O(1/g 2 ). To make the success probability a constant, we need to perform amplitude amplification [44]. Taking this into account the overall complexity increases by a factor of g. Using the bound on the condition number from Theorem 4 and bounds on query and gate complexities from Theorem 6, after plugging in the values from (6.4), we get the query and gate complexities of the algorithm as in the theorem statement.
7 Application to nonlinear differential equations Consider the nonlinear differential equation from [2], which is of the form [2], we assume that F 0 (and F 1 and F 2 ) are time-independent. However, we do not assume that F 1 is diagonalizable. We only assume that F 1 has a negative log-norm µ(F 1 ) < 0.
The procedure of Carleman linearization takes a quadratic (or higher order) ODE and converts it into a linear equation of higher dimension. For the quadratic ODE above, we get with the tri-diagonal block structure

4)
A j j = F 1 ⊗ I ⊗j−1 + I ⊗ F 1 ⊗ I ⊗j−2 + · · · + I ⊗j−1 ⊗ F 1 , (7.5) Note that A is a (3N s)-sparse matrix, where s is the sparsity of A. The dimension of (7.2) is We now define a quantity similar to the one that was used in [2] to quantify the amount of nonlinearity in the quadratic ODE (7.1). The difference is that we use log norm here instead of max λ Re(λ). Definition 6.
The following lemmas are needed to give a bound on the error in truncating the linearization procedure at level N . Lemma 15. For the quadratic ODE in (7.1), when R < 1 we have Proof. The proof of this lemma is given in the appendix.
Define the error due to truncation of Carleman procedure as follows.
where η 1 (t) is the error between the solution of (7.1) and x 1 (t) (which will be close to the output of the quantum algorithm). Let η(t) be the vector whose entries are η j for j = 1 to N . Lemma 16. For the quadratic ODE from (7.1), normalize the vector u(t) by the norm of the initial condition u 0 i.e.,ū = u u in , (7.11) which leads to the redefinition of F 0 and F 2 asF 0 = F 0 / u in andF 2 = F 2 u in (F 1 The exponential can be bounded as follows. giving us C(A) ≤ 1 , (7.22) when R ≤ 1 and µ(F 1 ) < 0.

Lemma 17. When R < 1 and
N ≥ 2 log(T F 2 /δ u(T ) ) log(1/ u in ) , (7.23) we have that the norm of the error η 1 (t) = x(t) − u(t) coming from truncating Carleman linearization to N steps can be bounded as Proof. First, note that the exact solution u(t) of the original quadratic ODE (7.1) satisfies (7.25) and the approximated solutionŷ j (t) satisfies (7.3). Comparing these equations, we have the tri-diagonal block structure which we write compactly as dη dt = Aη +b(t), η(0) = 0 (7.27) This can be integrated to give The norm of the error satisfies When R < 1 and µ(F 1 ) < 0, using Lemma 15, we have u ≤ u in and using Lemma 16, we have We can re-write the error as follows.
η ≤ N F 2 T u in N +1 . (7.32) This means that when N ≥ 2 log(T F 2 /δ u(T ) ) log(1/ u in ) , (7.33) we get η ≤ δ u(T ) . (7.34) We are now ready to prove the main result in this section. Theorem 8. There is an efficient quantum algorithm to solve the quadratic ODE from (7.1) i.e., to produce a quantum state close to a state proportional to the solution with query complexity (  1 ), log(T A ) , (7.35) and gate complexity larger by a factor of at most Proof. We apply the truncated Carleman linearization with N as above in (7.37).
We get a linear ODE dx dt = Ax + b , (7.40) with the following parameters.

Dimension of
(7.45) 6. Recall that the quantity g is defined in (5.91) as follows (7.46) Using (7.45) above and that max t∈[0,T ] u(t) = u in when R < 1, we get The initial state preparation is the same as in Lemma 5 of [2], where it is shown how to prepare the normalized version of |j ⊗ |u ⊗j in ⊗ |0 N −j , (7.50) where j labels the truncation level. This state is x 0 in the algorithm from Section 4 and is used to prepare the initial state for that algorithm. From Theorem 7, we can use the algorithm from Section 4 to get a quantum state y such that y − x(T ) ≤ δ x(T ) ≤ δ (1 + δ) √ N u(T ) . Now, measuring the quantum register containing the truncation level of Carleman linearization gives outcome 1 with probability (7.53) and the measured state is y 1 (T )/ y 1 (T ) . This is state is close to u(T )/ u(T ) . To see this, first note that Now using Theorem 7, we get the query and gate complexities as in the theorem statement.

Conclusions
We have presented a quantum algorithm to solve linear inhomogeneous ODEs. The gate and query complexity of our algorithm are bounded. Classes of matrices (namely those with nonpositive log-norm but with a large condition number for the diagonalizing matrix) where the bounds derived here are exponentially better than in previous work [1,12] are discussed. Our algorithm is also able to solve the linear ODE for many classes of non-diagonalizable matrices and even singular matrices. We construct a different linear system here because it is easier to apply our analysis to it, but it is interesting to see if the linear system in [1] has a small condition number (even when κ V is large). Our condition involving the matrix exponential will be particularly useful when considering block matrices such as in nonlinear ODEs. Bounding the matrix exponential of a block matrix in terms of the matrix exponentials of individual matrices is easier than checking for diagonalizability or other spectral conditions. An interesting open question is whether the dependence on A can be improved for stable matrices. For example, could A be replaced by a quantity proportional to the spectral radius ρ(A) of A, which can be considerably smaller than A for some matrices. It is known that ρ(A) ≤ A for all matrices. But for stable matrices, Gelfand's theorem says that for large enough k, there exists and such that A k ≤ (ρ + ) k . To be able to use this to pick the step size h, we need to find a quantitative bound relating and k. The theory of pseudospectra [36] gives bounds on norms of powers of matrices. It would be useful to explore this further to see a relationship between the -pseudosprectral radius ρ (A) and the norm of A.
In the realm of nonlinear differential equations, we have exponentially improved the dependence on error of previous work [2]. Our algorithm is also efficient for nonlinear differential equations when F 1 is non-normal (and even non-diagonalizable) if it has a negative log-norm. The linear ODE algorithm here can handle even more general cases (i.e., whenever C(A) is bounded and A is sparse). It would be interesting to extend our algorithm for nonlinear differential equations to these more general cases. For this, one needs to extend the definition of R and prove convergence of Carleman linearization to these more general matrices.