Linear multistep method, 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

Theorem 9 assumes the multistep method is of order $p$ and $A(\alpha)$-stable, that $A$ is diagonalisable as $A = VDV^{-1}$, and that its eigenvalues all satisfy $|\arg(-\lambda_i)| \leq \alpha$ — so the eigenvalues must sit in a wedge in the left half-plane, and where $A$ has an eigenvalue with positive real part the paper's remedy is to "simply subtract a multiple of the identity, and rescale the solution". Two further hypotheses: $\max_{t \in [t_0, t_0+\Delta t]} ||x(t)|| = O(||x(t_0+\Delta t)||)$, which rules out a solution that decays away from its own final value, and $\varepsilon = o(||x_{in}||)$. $A(\alpha)$-stability rather than $A$-stability is forced rather than chosen: "the second Dahlquist barrier is that an $A$-stable multistep method must be of order $p \leq 2$", and the point of the paper is $p > 2$. **The rigorous result is for constant coefficients**: the Conclusions state that the approach "can also be used to solve linear differential equations with time-dependent coefficients, though the error analysis will be more difficult", so unlike the Dyson and LCHS routes this one does not meet the slot's $A(t)$ contract with a proof behind it.

Requires

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

Example

Read in full: 1010.2745 carries no numerics at all. It is a construction and a complexity analysis end to end — Definitions 1 and 2, the stability polynomials, the block system of Eq. (8), and Theorem 9 — with no instance, no simulation and no figure. The register row records that as a full-text read rather than an abstract one.

given  A (sparse, constant coefficients), b, x_in, horizon dt, error budget e

choose a linear k-step method of order p that is A(alpha)-stable
    # A-stability would cap p at 2 (the second Dahlquist barrier), and the
    # whole point of going high-order is p > 2
    # alpha must be wide enough that every eigenvalue of A has |arg(-lambda_i)| <= alpha

set N_t = 2*dt/h  time steps
    # twice the horizon on purpose: x is held constant on [dt, 2dt], so that
    # measuring the clock register lands on the answer with probability Omega(1)
    # instead of 1/(N_t + 1)

assemble the sparse linear system  script_A x_vec = b_vec  whose rows are:
    row 0:            x_0 = x_in
    stepping rows:    sum_l alpha_l x_{j+l}  =  h * sum_l beta_l ( A x_{j+l} + b )
    holding rows:     x_{j+1} - x_j = 0        # the constant tail

hand the system to the layer below
    # the solution comes back as one state  sum_j |t_j> |x_j>  over a clock
    # register: the whole history at once, not the final time alone

measure the clock register; keep the outcomes in [dt, 2dt]
    # what is left is a state encoding x(dt), to within trace distance e

Cost, as the source states it

Theorem 9, the main result: $\tilde{O}( \log(N_x) s^{9/2} (||A||\Delta t)^{2+2/p} \kappa_V^5 (||x_{in}|| + ||b||/||A||)/\varepsilon^2 )$ calls to the oracles for $A$, $b$ and $x_{in}$, where $s$ is the sparsity, $N_x$ the dimension and $\kappa_V$ the condition number of the eigenvector matrix. Dropping the term for the error in *starting* the multistep method improves the other quantities without changing the $||A||\Delta t$ scaling: $\tilde{O}( \log(N_x) s^{9/2} (||A||\Delta t)^{2+2/p} \kappa_V^{2+4/p} (||x_{in}|| + ||b||/||A||)^{1/p}/\varepsilon^{1+2/p} )$. The Conclusions quote the headline in $||A||$ and $\Delta t$ alone: $\tilde{O}( (||A||\Delta t)^2 )$. The paper says outright that this is probably not the last word — "it is likely that this is suboptimal, because the lower bound is linear scaling", the no-fast-forwarding bound $O(||A||\Delta t)$ — and every later route on this slot is an attempt on that gap.

Implementations

• The clock-register linear system $\mathcal{A}\vec{x} = \vec{b}$ (sections III and IV)

  • High-order quantum algorithm for solving linear differential equations

  Why a linear system rather than a simulation, in the paper's own order. Section II sets up the alternative first and rejects it: for a general $A$, split into Hermitian and antiHermitian parts and combine the two evolutions by a Lie-Trotter formula, and "it appears to give a complexity that increases exponentially with the time interval $\Delta t$" — because the amplitude "must therefore decay exponentially, because we must allow for the maximum eigenvalue of $A_H$". That approach also "does not enable simulation of inhomogeneous differential equations". Section III's answer is the trick the rest of the construction rests on: a "Feynman's clock", an additional register encoding the time, so that one state $|\psi\rangle := \sum_{j=0}^{N_t}|t_j\rangle|x_j\rangle$ holds the solution at all times rather than the final time alone.

  The system is $\mathcal{A}\vec{x} = \vec{b}$ of Eq. (7), where "each entry of $\mathcal{A}$ is a block of the dimension of $A$, and each entry of $\vec{x}$ and $\vec{b}$ is a block of the dimension of $x$". Eq. (8) writes it out for the Euler case with $A$ and $b$ time-independent: the first block row is $\mathbb{1}$ against $x_{in}$ and "sets the initial value, $x_0 = x_{in}$"; the stepping rows are $-(\mathbb{1}+Ah)$ on the subdiagonal against $\mathbb{1}$ on the diagonal, so that "$x_{j+1} - (x_j + Ax_jh) = bh$"; and "the final rows indicate equations where $x_{j+1} - x_j = 0$", the tail that holds $x$ constant. Section IV replaces the stepping rows with a general linear multistep method, Eq. (9): $\sum_{\ell=0}^{k}\alpha_{\ell}x_{j+\ell} = h\sum_{\ell=0}^{k}\beta_{\ell}[A(t_{j+\ell})x_{j+\ell} + b(t_{j+\ell})]$, whose stability is read off the generating polynomials $\rho(\zeta) = \sum_j \alpha_j\zeta^j$ and $\sigma(\zeta) = \sum_j \beta_j\zeta^j$ through the roots of $\rho(\zeta) - \mu\sigma(\zeta) = 0$: the stability domain is $S := \{\mu \in \mathbb{C};\ \text{all roots } \zeta_j(\mu) \text{ satisfy } |\zeta_j(\mu)| \leq 1,\ \text{multiple roots satisfy } |\zeta_j(\mu)| < 1\}$, and the method has order $p$ if and only if $\rho(e^h) - h\sigma(e^h) = O(h^{p+1})$. Definition 2 is the property the algorithm needs: $A(\alpha)$-stable for $0 < \alpha < \pi/2$ means $S \supset S_{\alpha} = \{\mu;\ |\arg(-\mu)| < \alpha,\ \mu \neq 0\}$. **The doubling is load-bearing and easy to miss.** Reading the answer off the clock register succeeds with probability only $1/(N_t+1)$, so the paper solves to $t_0 + 2\Delta t$ with $x$ held constant on the second half — "then any measurement result for the time in this interval will give the state corresponding to the solution" — which is what makes $p_{ex} = \Omega(1)$ in Eq. (69) and costs nothing in the scaling of $N_t$.

  No dataset and no instantiated matrix. $A$, $b$ and $x_{in}$ stay symbolic throughout, reached only through oracles, and the sole concrete display is the five-block-row Euler example of Eq. (8) — drawn to make the block structure legible, with $A$ and $b$ still symbolic. The one quantitative choice the construction fixes is $N_t = 2\Delta t/h$.

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

• Taylor propagator, all-at-once encoding

  Encode a truncated Taylor series of the propagator into a single sparse linear system approximating the whole evolution, then solve it with a quantum linear system algorithm. This is what brought the precision dependence down to polynomial in $\log(1/\varepsilon)$.

• 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$.

• 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.

• 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.

In the Atlas

• Quantum algorithms for linear differential equations

  Given a first-order linear differential equation d/dt x = A(t)x + b(t) with N-dimensional vectors x and b and an N×N matrix A, and given an initial condition x(0), produce the solution x(t) at a later time t to precision ε, in the sense that the normalized vector x(t)/‖x(t)‖ returned is at distance at most ε from the exact solution.

Sources

• High-order quantum algorithm for solving linear differential equations↗Dominic W. Berry · 2010theory yes · simulation no · hardware no