Taylor propagator, all-at-once encoding

What it fills

• Solve a linear ODE du/dt = A(t)u + b(t)

  Given block-encoding access to $A(t)$ and $b(t)$ and a preparation unitary for $u_0$, output a normalized state $\varepsilon$-close to $u(T)/||u(T)||$. Matrix-query and state-preparation-query counts are stated separately, because methods here differ in them independently.

When it applies

Stated for systems of possibly inhomogeneous linear ODEs with constant coefficients. Berry, Childs, Ostrander and Wang describe the encoding as a sparse, well-conditioned linear system, and state that unlike with finite difference methods their approach does not require additional hypotheses to ensure numerical stability. It reduces to a quantum linear solve; it does not remove that layer.

Requires

Every step this method names moves its route along, so there is nothing it needs alongside them.

Example

Berry, Childs, Ostrander and Wang report no run of any kind. A full-text read (2026-08-12, which moved the register row for arxiv:1701.03684 from an abstract basis to a full-text one) found no figure, no table, no worked instance with numbers, no dataset, no simulator and no hardware anywhere in the paper — the only occurrences of "simulation" are in the sense of Hamiltonian simulation. There is therefore nothing to write up, and there is no neighbouring numerical work to describe either.

given  A (constant coefficients, diagonalizable A = V D V^-1, Re(lambda_i) <= 0),
       b, x_in with known norms ||x_in||, ||b||, horizon T, error budget e <= 1/2,
       oracles O_A for the entries of A, and controlled O_x, O_b

set    h     = T / ceil(T ||A||)
       m = p = T/h = ceil(T ||A||)
       g     = max_{t in [0,T]} ||x(t)|| / ||x(T)||
       delta = e / (25 sqrt(m) g)                   # the solver's error budget
       Omega = 70 g kappa_V m^{3/2} (||x_in|| + T ||b||) / (e ||x(T)||)
       k     = ceil( 2 log(Omega) / log log(Omega) )   # ensures (k+1)! >= Omega

# --- time-discretization, via truncated-taylor-propagator -------------------
assemble  C_{m,k,p}(Ah) x = |0>|x_in> + h sum_{i=0..m-1} |i(k+1)+1>|b>
    Taylor-term rows :  x_{i,j} = (Ah/j) x_{i,j-1},   2 <= j <= k
                        x_{i,1} = Ah x_{i,0} + h b
    step-closing row :  x_{i,0} = sum_{j=0..k} x_{i-1,j}
    padding rows     :  x_{m,j} = x_{m,j-1},          1 <= j <= p
    # sparse and well-conditioned: O(ks) nonzeros per row or column, and
    # condition number O(kappa_V k m) by Theorem 5.  Unlike with finite
    # difference methods, no additional hypothesis is needed to ensure
    # numerical stability.

repeat O(g) times, coherently -- the amplitude amplification rounds:

    # --- own work: prepare the right-hand side ---
    prepare |0>|x_in> + h sum_{i=0..m-1} |i(k+1)+1>|b>
        # a constant number of calls to O_x and O_b, poly(log(mk)) gates

    # --- quantum-linear-solve ---
    solve C_{m,k,p}(Ah) x = (that state) with a QLSA, to within delta of the
        normalised solution
        # this route reduces to a quantum linear solve; it does not remove
        # that layer

    # --- own work: extract the answer ---
    measure the first register in the standard basis
    accept iff the outcome lies in S = {m(k+1), ..., m(k+1)+p}
        # the last p+1 blocks all hold x_{m,0}, the state at time T
        # one round accepts with probability at least 1/78 g^2  (p = m)

return the second register -- a state e-close to x(T)/||x(T)||, with a flag
       indicating success

Cost, as the source states it

The complexity is polynomial in the logarithm of the inverse error, an exponential improvement over previous quantum algorithms for this problem. No matrix-query or state-preparation query count appears in the abstract; a $\kappa_V$-dependent expression often attached to this method in secondary summaries belongs to the spectral-method row of a later comparison table, so it is not reproduced here.

Implementations

• The sparse-oracle construction of Theorem 9 (sections 2, 6 and 7)

  • Quantum algorithm for linear differential equations with exponentially improved dependence on precision

  The algorithm the paper states and proves as its main result. Section 1 gives the reason it was built: the earlier approach handed a linear multistep method to a quantum linear systems algorithm, and its complexity stayed $\mathrm{poly}(1/\varepsilon)$ even when the improved solver was substituted, because the multistep method is itself a significant source of error. Direct use of a linear combination of unitaries was considered and rejected — for general non-anti-Hermitian $A$ the best implementation the authors are aware of has an exponentially small success probability — so a truncation of the Taylor series of the propagator $\exp(At)$ is encoded into a linear system instead, which the paper describes as conceptually similar to a linear combination of unitaries but with significantly better performance because the combination is built stepwise through a system of linear equations. The problem statement is Theorem 9's: $A = VDV^{-1}$ is an $N \times N$ diagonalizable matrix whose eigenvalues all satisfy $\operatorname{Re}(\lambda_j) \le 0$, with at most $s$ nonzero entries in any row and column and an oracle $O_A$ computing them; $\vec{x}_{\mathrm{in}}$ and $\vec{b}$ have known norms and controlled oracles $O_x$ and $O_b$ preparing states proportional to them; the output is a state $\varepsilon$-close in $\ell_2$ norm to $\vec{x}(T)/\lVert\vec{x}(T)\rVert$ with a flag indicating success, obtained with probability $\Omega(1)$.

  The parameters are all fixed in one place, the Statement of the Algorithm in section 7: $h = T/\lceil T\lVert A\rVert\rceil$, $m = p = T/h = \lceil T\lVert A\rVert\rceil$, $\delta = \varepsilon/(25\sqrt{m}\,g)$ with $\varepsilon \le 1/2$, and $k = \lceil 2\log\Omega/\log\log\Omega\rceil$ where $\Omega = 2.8\kappa_V m(\lVert x_{\mathrm{in}}\rVert + mh\lVert b\rVert)/(\delta\lVert x(T)\rVert) = 70g\kappa_V m^{3/2}(\lVert x_{\mathrm{in}}\rVert + T\lVert b\rVert)/(\varepsilon\lVert x(T)\rVert)$. That choice of $k$ ensures $(k+1)! \ge \Omega$, and because $\Omega \ge 70$ it also ensures $k \ge 5$, which is the hypothesis the conditioning and success-probability theorems need. The system built is $C_{m,k,p}(Ah)\lvert x\rangle = \lvert 0\rangle\lvert x_{\mathrm{in}}\rangle + h\sum_{i=0}^{m-1}\lvert i(k+1)+1\rangle\lvert b\rangle$. Its right-hand side is prepared by Lemma 8 of section 6, in three moves: a rotation on the block register that weights $\lVert\vec{x}_{\mathrm{in}}\rVert$ against $\sqrt{m}h\lVert\vec{b}\rVert$, then $O_x$ and $O_b$ applied in arbitrary order, then a unitary that fixes $\lvert 0\rangle$ and sends $\lvert 1\rangle$ to $\frac{1}{\sqrt{m}}\sum_{j=0}^{m-1}\lvert j(k+1)+1\rangle$, spreading the inhomogeneous part over the time-step blocks — a constant number of calls to $O_x$ and $O_b$ plus $\mathrm{poly}(\log(mk))$ elementary gates. The system is solved with the linear systems algorithm of Childs, Kothari and Somma, which section 1 singles out because it avoids phase estimation and inverts the system directly by a linear combination of unitaries, giving a state within $\delta$ of the normalised solution. The first register is then measured in the standard basis and the outcome accepted when it lies in $S = \{m(k+1), m(k+1)+1, \ldots, m(k+1)+p\}$, the block labels that all carry $\lvert x_{m,0}\rangle$. Theorem 7 with $p = m$ puts one round's success probability at $1/(78g^2)$ for the exact solution vector, and Eq. (127) at $1/(121g^2)$ once the solver's own error $\delta$ is carried through. $O(g)$ rounds of amplitude amplification raise it to $\Omega(1)$; section 2 records classical repetition as the alternative to amplification.

  No dataset. $A$ enters only through the oracle $O_A$ that computes its at most $s$ nonzero entries per row and column, and the two vectors only through controlled oracles $O_x$ and $O_b$ preparing $\vec{x}_{\mathrm{in}}/\lVert\vec{x}_{\mathrm{in}}\rVert$ and $\vec{b}/\lVert\vec{b}\rVert$, with both norms assumed known. Section 6 marks this input model as a deliberate change from the multistep algorithm it is compared against, which assumes $\vec{x}_{\mathrm{in}}$ and $\vec{b}$ are sparse vectors whose entries are given by oracles: instead of assuming sparsity, this paper simply assumes controlled oracles that produce states proportional to the two vectors.

  Counted, never run. The solve alone costs $O(\kappa_V k^2 m s \cdot \mathrm{poly}(\log(\kappa_V kms/\delta)))$ queries, using Theorem 5's $O(\kappa_V km)$ bound on the condition number of $C_{m,k,p}(Ah)$ at these parameters and the fact that the matrix is $(d+1)N \times (d+1)N$ with $O(ks)$ nonzero entries in any row or column. With the $O(g)$ amplification rounds, Theorem 9 gives the whole algorithm $O(\kappa_V sgT\lVert A\rVert \cdot \mathrm{poly}(\log(\kappa_V sg\beta T\lVert A\rVert/\varepsilon)))$ queries to $O_A$, $O_x$ and $O_b$, where $\beta = (\lVert x_{\mathrm{in}}\rVert + T\lVert b\rVert)/\lVert x(T)\rVert$, and a gate complexity larger by a factor $\mathrm{poly}(\log(\kappa_V sg\beta T\lVert A\rVert N/\varepsilon))$. Section 8 argues that neither of the two parameters this is charged in can be improved much: evolution in which the state decays can be used to implement postselection, and the ability to postselect on an exponentially small amplitude would imply BQP = PP, so the dependence on $g$ is close to necessary, while Hamiltonian simulation is the special case $g = 1$ and the no-fast-forwarding theorem then makes the dependence on $T$ optimal up to logarithmic factors. Nothing here was measured. The paper reports no numerical experiment, contains no figure and no table, and names no hardware and no simulator, so this implementation exists as a specified and analysed algorithm rather than an executed one.

• Non-diagonalizable $A$ through a nearby diagonalizable matrix (section 8)

  • Quantum algorithm for linear differential equations with exponentially improved dependence on precision

  Theorem 9 assumes $A$ is diagonalizable, and section 8 states that the algorithm can nonetheless produce approximate solutions when it is not. This is the paper's own variant for that case, offered as a remark rather than as a second theorem.

  Diagonalizable matrices are dense within the set of all complex matrices: for any non-diagonalizable $A$ and any $\delta > 0$ there is a diagonalizable $B$ with $\lVert A - B\rVert < \delta$. The construction is then run unchanged with $B$ in place of $A$, so that Eq. (1) is simulated approximately.

  The paper states the price in the same sentence it offers the variant: this approach can yield a matrix $B$ whose diagonalizing transformation has a condition number polynomial in $1/\delta$, so the complexity would no longer be $\mathrm{poly}(\log(1/\varepsilon))$. The variant keeps the algorithm applicable and gives up the property the paper was written for. No bound is proved for it and no numerics are reported.

• Hamiltonian simulation as the $g = 1$ specialisation (sections 1 and 8)

  • Quantum algorithm for linear differential equations with exponentially improved dependence on precision

  The specialisation the paper names for itself. Section 1 opens on differential equations of the form $\dot{\vec{x}} = A\vec{x}$ with $A$ anti-Hermitian, and section 8 says outright that Hamiltonian simulation is a special case of this algorithm with $g = 1$.

  The inhomogeneity is absent and the generator is anti-Hermitian, so the solution does not decay and the parameter $g = \max_{t \in [0,T]}\lVert\vec{x}(t)\rVert/\lVert\vec{x}(T)\rVert$ of Eq. (111) — the quantity the success probability and the $O(g)$ amplification rounds are charged against — takes the value one.

  The paper uses the specialisation as an argument rather than as an algorithm to run: because Hamiltonian simulation is this case, the no-fast-forwarding theorem implies that the algorithm's dependence on the evolution time $T$ is optimal up to logarithmic factors. Nothing is computed, and no separate construction is written down for it.

What it needs

1. Choose a time discretization or propagator approximation 6 methods

   Reduce continuous evolution over $[0,T]$ to a finite algebraic object — a banded linear system, a product of step propagators, or a spectral coefficient system — with a stated truncation error. When a linear system is formed, a conditioning bound is stated with it.

2. Quantum linear solve 5 methods

   Given access to a matrix $A$ and a unitary that prepares $|b>$, produce a flagged quantum state that is $\varepsilon$-close in $l2$ to the normalised $A^{-1}b$. The deliverable is a state, not a classical vector.

Other ways to fill the same slot

Different approaches

• Linear multistep method, all-at-once encoding

  The first quantum algorithm for general linear ODEs. Discretize with a high-order linear multistep method, lay every time step out at once against a clock register so that one state $\sum_j |t_j\rangle|x_j\rangle$ holds the whole history, and solve the resulting sparse linear system with a quantum linear system algorithm. Going high-order is what buys the $\Delta t$ scaling: plain Euler costs at least $O(\Delta t^4)$ here.

• Chebyshev spectral method, global collocation

  The route that brought $\mathrm{poly}(\log(1/\varepsilon))$ precision to linear ODEs with **time-dependent** coefficients, which is what its abstract says was missing: "no such algorithm was previously known for differential equations with time-dependent coefficients". It does it by not stepping. The solution is approximated globally by a truncated Chebyshev series, the coefficients are fixed by collocating the differential equation at Chebyshev nodes, and the resulting sparse system goes to a high-precision quantum linear system algorithm. The exponential precision is bought by smoothness rather than by the solver: it is the $C^{\infty}$ convergence of the Chebyshev series that makes the series length logarithmic in $1/\varepsilon$.

• Dyson propagator, all-at-once encoding

  Encode the Dyson series in a system of linear equations and solve it via the optimal quantum linear equation solver, extending the all-at-once approach to genuinely time-dependent generators.

• Time-marching with uniform singular value amplification

  Propagate the solution one step at a time and defeat the exponentially vanishing success probability by repeatedly invoking uniform singular value amplification, improved further by a compression gadget lemma. Fang, Lin and Tong present it explicitly as a design path alternative to solvers based on quantum linear systems algorithms.

• LCHS — linear combination of Hamiltonian simulation

  Express a general non-unitary evolution operator as a linear combination of unitary evolution operators, each of which solves a Hamiltonian simulation problem, rather than converting the problem into a dilated linear system. An, Liu and Lin state that the method can achieve optimal cost in terms of state preparation.

• LCHS with the improved kernel a narrower version of LCHS — linear combination of Hamiltonian simulation

  A family of identities expressing non-unitary evolution as a linear combination of unitary evolutions, built on the kernel $f(z) = 1/(C_β e^{(1+iz)^β})$ with $β ∈ (0,1)$ and $C_β = 2π e^{-2^β}$. The kernel decays at a near-exponential rate $e^{-c|k|^β}$, replacing the original Cauchy kernel's quadratic decay and exponentially enhancing accuracy.

• Schrödingerisation (linear PDEs as Schrödinger equations)

  A simple change of variable — the warped phase transformation, which introduces one extra variable — recasts any linear PDE or ODE system into a system of Schrödinger equations in real time, which ordinary Hamiltonian simulation then runs. The original solution is recovered from the auxiliary dimension.

Narrower versions of this one

• Krovi's reanalysis of the all-at-once encoding a narrower version of Taylor propagator, all-at-once encoding

  Reanalyses the all-at-once propagator encoding and shows that the norm of the matrix exponential, rather than the eigenvector condition number, characterizes the run time. It still forms a global linear system and still calls a quantum linear solver.

In the Atlas

No record in the Atlas covers this yet. The catalogue is circuits and primitives; this part of the literature is not in it.

Sources

• Quantum algorithm for linear differential equations with exponentially improved dependence on precision↗Dominic W. Berry, Andrew M. Childs, Aaron Ostrander, Guoming Wang · 2017theory yes · simulation no · hardware no