Chebyshev spectral method, global collocation

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

**Problem 1** fixes the access model and one smoothness assumption before any theorem runs: $A(t) \in \mathbb{C}^{d\times d}$ is $s$-sparse for every $t \in [0,T]$, and "we assume that $A_{ij}, f_i \in C^{\infty}(0,T)$ for all $i,j \in [d]$". Everything is reached through oracles — one computing "the locations and values of the nonzero entries of $A(t)$" in a given row or column, and oracles preparing $|\gamma\rangle$ and $|f(t)\rangle$ that "also compute $\|\gamma\|$ and $\|f(t)\|$" — and the deliverable is a quantum state, "a quantum state $\varepsilon$-close to $|x(T)\rangle$", never an explicit vector. **Theorem 1** adds the spectral hypothesis: "$A(t)$ can be diagonalized as $A(t) = V(t)\Lambda(t)V^{-1}(t)$ where $\Lambda(t) = \mathrm{diag}(\lambda_1(t),\ldots,\lambda_d(t))$ with $\mathrm{Re}(\lambda_i(t)) \leq 0$ for each $i \in [d]_0$ and $t \in [0,T]$". The paper gives the reason and the comparison in one breath: the non-positive real parts are required "so that the solution cannot grow exponentially", and this "is essentially the same as in the time-independent case considered in [10] and improves upon the analogous condition in [5] (which requires an additional stability condition)" — [5] being Berry's multistep route on this same slot, whose $A(\alpha)$-stability is the extra condition being improved upon. **The hypothesis the abstract does not carry is smoothness of the solution, and the headline rests on it.** Section 2: "if the solution is in $C^{\infty}$, the spectral method approximates the solution to within $\varepsilon$ using only $n = \mathrm{poly}(\log(1/\varepsilon))$ terms in the Chebyshev series", whereas "for a solution in $C^{r+1}$, the spectral method approximates the solution with $n = \mathrm{poly}(1/\varepsilon)$". The Discussion lists it as an open problem rather than a technicality: "our algorithm must assume that the solution is smooth. If the solution is in $C^r$, the solution error is $O(1/n^{r-2})$ by Lemma 1. Can we improve the complexity to $\mathrm{poly}(\log(1/\varepsilon))$ under such weaker smoothness assumptions?" Diagonalisability is the softer of the two hypotheses and the paper says so: "our algorithm can produce approximate solutions for non-diagonalizable $A(t)$, although the dependence on $\varepsilon$ degrades to $\mathrm{poly}(1/\varepsilon)$".

Requires

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

Example

Read in full: arXiv:1901.00961, the only paper this record cites, reports no run of the method. There is no figure, no table, no computed value, no dataset, no simulator and no hardware anywhere in it — the word "simulation" occurs only in the sense of Hamiltonian simulation and in reference titles. Its ten sections are a construction and its analysis end to end: the spectral method, the linear system, the solution error, the condition number, the success probability, state preparation, the main result and the boundary value problem. Appendix B, "An example of the quantum spectral method", is a structural display and not a run — $d = 1$, $m = 3$, $n = 2$, $p = 1$, with $A(t)$, $f(t)$ and $\gamma$ left symbolic. So there is no run to transcribe here, and no nearby numerics to point at instead.

given  oracles for A(t) (s-sparse), f(t), gamma; horizon T; error budget e

choose m so that  ||A||*T / (2*m) <= 1
    # NOT an accuracy choice. Rescaling [0,T] onto [-1,1] sends A -> -(T/2)A,
    # which 'can dramatically increase the spectral norm'; subdividing avoids it
divide [0,T] into [0,G_1], [G_1,G_2], ... , [G_{m-1},T]
    map each [G_h, G_{h+1}] onto [-1,1] by K_h, with K_h(G_h)=1, K_h(G_{h+1})=-1

choose the series length n from Eq. (8.6)
    n = (e/2) * max( floor(log(Omega)/log(log(Omega))),
                     floor(log(omega)/log(log(omega))) )
    Omega = g'*e*m*(1+eps)/(g*eps),   omega = (g'/||gamma||)*(m+1)
    # poly-logarithmic in 1/eps only because the solution is assumed C^infinity

for each subinterval, hand the layer above's collocation to this one:
    build the rows of  (L_1 + L_2(A_h)) |X> = |B_h>
    L_3 carries the coefficients c_{i,l}(G_{h+1}) into the next block
        as its initial condition
assemble one system  L |X> = |B>          # Eq. (3.12)
    L is (m+p+1)*d*(n+1) square, with O(n*s) nonzeros per row
    pad the final state p = O(m) = O(||A||*T) times via L_4, L_5

hand L |X> = |B> to the layer below          # the high-precision QLSA

measure the first and third registers
    keep outcomes in S = {m, m+1, ... , m+p}
    # success probability (p+1)(n+1) / (pi*m*q^2 + (p+1)(n+1)) by Lemma 5
amplitude-amplify O(q/sqrt(n)) times to reach success probability Omega(1)

return the second register: a state e-close to x(T)/||x(T)||

Cost, as the source states it

**Theorem 1**, the main result: $O(\kappa_V s \|A\| T q\,\mathrm{poly}(\log(\kappa_V s \|A\| g' T/\varepsilon g)))$ queries to the oracles $O_A(h,l)$, $O_x$ and $O_f(h,l)$, with a gate complexity "larger than the query complexity by a factor of $\mathrm{poly}(\log(\kappa_V d s \|A\| g' T/\varepsilon))$". Eq. (8.2) names every parameter rather than leaving them to a reader: $\|A\| := \max_{t\in[0,T]}\|A(t)\|$; $\kappa_V := \max_t \kappa_V(t)$, where $\kappa_V(t)$ is the condition number of the eigenvector matrix $V(t)$; $g := \|\hat{x}(T)\|$; $g' := \max_{t\in[0,T]}\max_{n\in\mathbb{N}}\|\hat{x}^{(n+1)}(t)\|$; and $q := \max_{t\in[0,T]}\|\hat{x}(t)\|/\|x(T)\|$. **Corollary 1** specialises to time-independent equations, replacing $g'$ by $\|\gamma\| + 2\tau\|f\|$: $O(\kappa_V s\|A\|Tq\,\mathrm{poly}(\log(\kappa_V s\gamma\|A\|\|f\|T/\varepsilon g)))$. **The boundary value problem is much worse in the horizon** — Theorem 2 gives $O(\kappa_V s\|A\|^{4}T^{4}q\,\mathrm{poly}(\log(\kappa_V s\|A\|g'T/\varepsilon g)))$, a fourth power where the initial value problem has a first. Two of these parameters carry the paper's own warnings and neither is decoration. Of $g'$: "in general, $g'$ could be unbounded above as $n \to \infty$", and removing it is left open — "the query complexity of Hamiltonian simulation is independent of that parameter. Can we develop quantum algorithms for general differential equations with query complexity independent of $g'$?" Of $q$, which "characterizes the decay of the final state relative to the initial state": "it is unlikely that the dependence on $q$ can be significantly improved, since renormalization of the state effectively implements postselection and an efficient procedure for performing this would have the unlikely consequence $\mathsf{BQP} = \mathsf{PP}$". The horizon dependence, by contrast, is close to the floor: "our algorithm has nearly optimal dependence on $T$, scaling as $O(T\,\mathrm{poly}(\log T))$", against the no-fast-forwarding bound that "the complexity must be at least linear in $T$".

Implementations

• The global Chebyshev linear system $L|X\rangle = |B\rangle$ (sections 2 and 3)

  • Quantum spectral methods for differential equations

  Why a global approximation, in the paper's own order. Section 1 sets out the obstacle first: "most of the aforementioned algorithms use a local approximation: they discretize the differential equations into small time intervals", and even at high order "when solving an equation over the interval $[0,T]$, the number of iterations is $T/h = \Theta(\varepsilon^{-1/k})$ for fixed $k$, giving a total complexity that is $\mathrm{poly}(1/\varepsilon)$ even using high-precision methods for the QLSA or Hamiltonian simulation". Two escapes already existed and both need structure the general case does not have: "when $A(t)$ is anti-Hermitian and $f(t) = 0$, we can directly apply Hamiltonian simulation; if $A$ and $f$ are time-independent, then [10] uses a Taylor series to achieve complexity $\mathrm{poly}(\log(1/\varepsilon))$. However, the case of general time-dependent linear ODEs had remained elusive." Section 2's answer is the one this record is about: represent the components "as linear combinations of basis functions $\phi_j(t)$ expressing the time dependence", so that no step count stands between the horizon and the error.

  The approximation is Eq. (2.1), $x_i(t) = \sum_{k=0}^{n}c_{i,k}T_k(t)$, and the coefficients "are determined by demanding that $x(t)$ satisfies the ODE and initial conditions at a set of interpolation nodes $\{t_l\}_{l=0}^{n}$ (with $1 = t_0 > t_1 > \cdots > t_n = -1$), where $x(t_0)$ and $x(t_n)$ are the initial and final states" — Eq. (2.2) and Eq. (2.3). The nodes are Chebyshev-Gauss-Lobatto, $t_l = \cos(l\pi/n)$, "since these nodes achieve the highest convergence rate among all schemes with the same number of nodes", and because "these nodes also have the convenient property that $T_k(t_l) = \cos(kl\pi/n)$". Differentiation stays inside the basis: from $2T_k(t) = T'_{k+1}(t)/(k+1) - T'_{k-1}(t)/(k-1)$ comes the upper-triangular matrix $D_n$ of Eq. (2.8), $[D_n]_{kj} = 2j/\sigma_k$ for $k+j$ odd and $j > k$, with $\sigma_0 = 2$ and $\sigma_k = 1$ otherwise, giving Eq. (2.10). Section 3 assembles this into $L|X\rangle = |B\rangle$, Eq. (3.12), over $m$ rescaled subintervals with initial conditions chained by $x_h(1) = x_{h-1}(-1)$, and pads the final state $(p+1)(n+1)$ times through $L_4$ and $L_5$ so the readout succeeds with constant probability. **The subdivision is the part most easily mistaken for a step count.** It is not chosen against the error budget at all: rescaling $[0,T]$ onto $[-1,1]$ sends $A \mapsto -(T/2)A$, "which can dramatically increase the spectral norm", so the subintervals are sized by Eq. (3.8), $\tau_h \leq 2/\max_{t\in[\Gamma_h,\Gamma_{h+1}]}\|A(t)\|$. The error is controlled by $n$, and $n$ is logarithmic in $1/\varepsilon$ by Lemma 2.

  No dataset and no instantiated matrix. $A(t)$, $f(t)$ and $\gamma$ stay symbolic throughout and are reached only through oracles, exactly as on `berry-multistep`. The sole concrete display is Appendix B, "An example of the quantum spectral method", which fixes $d = 1$, $m = 3$, $n = 2$ and $p = 1$ and writes out the resulting block structure of $L$ — drawn to make that structure legible, with $A(t)$ and $f(t)$ still symbolic and no value computed. The quantitative choices the construction itself fixes are the node set $t_l = \cos(l\pi/n)$, the subinterval bound $\tau_h \leq 2/\max_t\|A(t)\|$, and $p = O(m) = O(\|A\|T)$.

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.

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

• 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

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 spectral methods for differential equations↗Andrew M. Childs, Jin-Peng Liu · 2019theory yes · simulation no · hardware no