Time-marching with uniform singular value amplification

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

It surpasses existing QLSA-based solvers in three respects stated by the authors: $A(t)$ need not be diagonalizable; $A(t)$ may be non-smooth and is only required to be of bounded variation; and it can use fewer queries to the initial state. It needs no global linear solve, which is what makes the bypass real rather than rhetorical. The trade is in matrix queries: An, Childs and Lin's comparison table gives time-marching a worse dependence on $\alpha_A T$ than the LCHS family while crediting it with the same initial-state-preparation advantage.

Requires

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

Example

The only numerical experiment in arXiv:2208.06941 — Section 2.3 with Figure 2, and Appendix C — measures the degree of the odd polynomial that uniform singular value amplification hands to QSVT, and nothing else. No differential equation is integrated, no quantum state is prepared, and no value of the problem's own parameters enters anywhere.

given  A(t) on [0,T] with alpha = sup_t ||A(t)||,  Uinit preparing |psi(0)>,  tolerance eps

choose the mesh 0 = t_0 < ... < t_L = T  with  t_l - t_{l-1} <= 1/(2 alpha)
    # so L = Theta(alpha T) segments

for l = 1 ... L:                    # once per time step
    # the time-discretization slot, invoked per segment, not once for the whole horizon
    U_l <- block encoding of  Xi_l = T-ordered exp( int_{t_{l-1}}^{t_l} A(t) dt )
           # truncated Dyson series (Sec 3.1), or the first-order Magnus step (Sec 5)
           # consistent to ||Xi_l - Xibar_l|| <= eps_l ||Xi_l||, with sum_l eps_l <= 1/2

    # the amplification, at every turn: subnormalization alpha_l -> ||Xi_l||/(1-delta)
    Utilde_l <- uniform singular value amplification of U_l, delta = 1/L
                # oblivious to the state, so no earlier step is re-prepared

# no all-at-once system is assembled and no quantum linear solve is called

counter <- |L>  on ceil(log2 L) + 1 qubits
apply Utilde_1, ..., Utilde_L in order, each decrementing counter by controlled ADD-dagger
    # coherent: nothing is measured between turns

post-select counter = |0>              # succeeds with probability Omega(Q^-2)
amplitude-amplify, O(Q) rounds        # Q = prod_l ||Xi_l|| / |||psi(T)>||
return the state on the flagged branch

Cost, as the source states it

The complexity depends linearly on the amplification ratio, which quantifies the deviation from a unitary dynamics; that linear dependence is proved to attain the query complexity lower bound and thus cannot be improved in the worst case.

Implementations

• Time-marching paired with the high-order truncated Dyson series (Theorem 8)

  • Time-marching based quantum solvers for time-dependent linear differential equations

  The paper's headline algorithm, and the one its optimality result is about. The problem is Equation (1): prepare the solution of $\frac{d}{dt}|\psi(t)\rangle = A(t)|\psi(t)\rangle$ at time $T$, for a coefficient matrix $A$ of bounded variation on $[0,T]$, given a unitary circuit $U_{\mathrm{init}}$ with $U_{\mathrm{init}}|0^n\rangle = |\psi(0)\rangle$. Fang, Lin and Tong present it as a design route alternative to solvers built on quantum linear systems algorithms rather than as an improvement inside one.

  The mesh satisfies $t_l - t_{l-1} \le (2\alpha)^{-1}$ for every segment, and each segment carries a time-dependent matrix encoding $\mathrm{MAT}_l$, an $(n_q, m, t_{l-1}, t_l, \alpha, \varepsilon'')$-MAT in the sense of the paper's Definition 1 with $\varepsilon'' < \varepsilon/(2TQ)$, out of which a block encoding of the high-order truncated Dyson series approximation to $\Xi_l$ is built. Every step's block encoding is then rebuilt by uniform singular value amplification at $\delta = 1/L$, the whole chain is made coherent by the compression gadget of Lemma 3 so that the ancilla register need not be duplicated $L$ times, and $O(Q)$ rounds of amplitude amplification lift the $\Omega(Q^{-2})$ success probability to at least $2/3$.

  What the theorem proves, rather than anything measured: $O\!\left(\alpha^2 T^2 Q \log(\alpha TQ)\, \frac{\log(\alpha TQ\varepsilon^{-1})}{\log\log(\alpha TQ\varepsilon^{-1})}\right)$ queries to all $\mathrm{MAT}_l$ and $O(Q)$ applications of controlled-$U_{\mathrm{init}}$ and its inverse, using $O(n + m + \mathrm{polylog}(V_0^T(A)\alpha TQ\varepsilon^{-1}))$ qubits and $\tilde{O}(\alpha^2 T^2 Q (m + \mathrm{polylog}(V_0^T(A)\alpha TQ\varepsilon^{-1})))$ additional elementary gates, preparing a state within $O(\varepsilon)$ of the normalized solution with probability at least $2/3$, success flagged by the measurement result of a single qubit. The authors report no execution of this algorithm on hardware or on a simulator.

• Time-marching paired with the first-order truncated Magnus integrator (Section 5, Theorem 14)

  • Time-marching based quantum solvers for time-dependent linear differential equations

  Section 5's simplified implementation. Implementing the high-order truncated Dyson series requires complicated quantum control logic for handling time-ordering operators, so the authors pair the same marching with the cheapest integrator that avoids that logic entirely: $\bar{\Xi} = \exp\left(\int_{t_j}^{t_{j+1}} A(s)\,ds\right)$, the first-order truncated Magnus series, which is the time-ordered propagator with the time-ordering operator simply dropped. In Hamiltonian simulation the same step is the quantum highly oscillatory protocol of An, Fang and Lin, and it carries commutator scaling in the high-precision limit. When $A(t)$ commutes with itself at any pair of times the step is exact, and the precision dependence of the query complexity becomes polylogarithmic.

  The integral over a segment is approximated by a Riemann sum with $M = 2^{n_q}$ quadrature points, $\int_a^b A(s)\,ds \approx \frac{b-a}{M}\sum_{k=0}^{M-1} A(a + k(b-a)/M)$, and that sum is block encoded by applying a Hadamard on each of the $n_q$ quadrature qubits, then $\mathrm{MAT}_j$, then uncomputing the Hadamards. The resulting time-independent matrix exponential is then applied through the contour integral $e^A = \frac{1}{2\pi i}\oint_{\Gamma} e^z (z - A)^{-1}\,dz$, discretized by numerical quadrature into a linear combination of resolvents $(z_j - A)^{-1}$ each implemented with QSVT. The contour route is what the non-normal case forces: for a general $A$ the singular value transformation no longer agrees with the eigenvalue transformation, so the exponential cannot be applied by transforming singular values directly. Everything around the step, the amplification at each turn, the compression gadget and the $O(Q)$ rounds of amplitude amplification, is unchanged from the Dyson pairing.

  Theorem 14 states the general query count in terms of $\max\{\alpha_{\mathrm{comm}}T^2Q/\varepsilon,\ \alpha T\}$, where $\alpha_{\mathrm{comm}} := \sup_{s,\tau \in [0,T]}\|[A(s), A(\tau)]\|$. In the high-precision limit with a non-zero $\alpha_{\mathrm{comm}}$ it simplifies to $O\!\left(\frac{\alpha_{\mathrm{comm}}^2 T^4 Q^3}{\varepsilon^2}\log^2\frac{\max\{\alpha_{\mathrm{comm}}, \alpha\}TQ}{\varepsilon}\right)$ queries to all $\mathrm{MAT}_l$, together with $O(Q)$ applications of controlled-$U_{\mathrm{init}}$ and its inverse, on $O(n + m + \mathrm{polylog}(V_0^T(A)\alpha TQ\varepsilon^{-1}))$ qubits, with success flagged by the measurement result of a single qubit. When $\alpha_{\mathrm{comm}} = 0$ the authors note this first-order implementation has the same complexity scaling as the high-order truncated Dyson series. These are proved bounds; no run is reported.

• Sparse-matrix input model (Corollary 9)

  • Time-marching based quantum solvers for time-dependent linear differential equations

  The Theorem 8 algorithm re-costed on sparse access rather than on a supplied time-dependent matrix encoding. Fang, Lin and Tong give it because that input model is the one other quantum-linear-systems-based differential equation solvers and time-dependent Hamiltonian simulation algorithms are stated in, which is what makes a comparison against them a like-for-like one.

  Gilyén, Su, Low and Wiebe's Lemma 48 builds the $(n_q, m, a, b, \alpha, \varepsilon)$-MAT that Theorem 8 needs out of those three oracles, at a single query each to $U_{\mathrm{row}}$ and $U_{\mathrm{col}}$ and two queries to $U_{\mathrm{val}}$. The precision $\varepsilon$ can be made arbitrarily small, and holding the error below it costs $O(n + \log^{5/2}(d\varepsilon^{-1}))$ additional elementary gates and $O(n_b + \log^{5/2}(d\varepsilon^{-1}))$ additional ancilla qubits, where $n_b$ is the number of bits encoding $A_{jk}(t)$; those ancillas can be reused. Nothing in the marching itself changes.

  $A(t)$ is $d$-sparse with $\|A(t)\| \le 1$, and the locations of its non-zero elements do not move with time. Access is three oracles: $U_{\mathrm{row}}|j,s\rangle = |j, \mathrm{row}(j,s)\rangle$, $U_{\mathrm{col}}|j,s\rangle = |j, \mathrm{col}(j,s)\rangle$ and $U_{\mathrm{val}}|t,j,k,z\rangle = |t,j,k,z \oplus A_{jk}(t)\rangle$, where $\mathrm{row}(j,s)$ is the row index of the $s$-th non-zero element in the $j$-th column and $\mathrm{col}(j,s)$ the column index of the $s$-th non-zero element in the $j$-th row.

  Corollary 9's proved cost: $O\!\left(d^2 T^2 Q \log(dTQ)\, \frac{\log(dTQ\varepsilon^{-1})}{\log\log(dTQ\varepsilon^{-1})}\right)$ applications of controlled-$U_{\mathrm{row}}$, $U_{\mathrm{col}}$, $U_{\mathrm{val}}$ and their inverses, and $O(Q)$ applications of controlled-$U_{\mathrm{init}}$ and its inverse, on $O(n + \mathrm{polylog}(V_0^T(A)dTQ\varepsilon^{-1}))$ qubits with $\tilde{O}(d^2 T^2 Q(n + \mathrm{polylog}(V_0^T(A)dTQ\varepsilon^{-1})))$ additional elementary gates, at success probability at least $2/3$, flagged by the measurement result of a single qubit. Proved, not measured.

• Convex-optimization construction of the amplifying polynomial (Appendix C)

  • Time-marching based quantum solvers for time-dependent linear differential equations

  A replacement for the polynomial the amplification lemma is quoted from. Gilyén, Su, Low and Wiebe's construction approximates $(1-\delta)\,\mathrm{XRect}(x)$, and because $\mathrm{XRect}(x) = \gamma' x\,\mathrm{Rect}(x)$ is discontinuous at $x = \gamma'^{-1}$ the Gibbs phenomenon forces the $(1-\delta)$ factor and leaves a large preconstant, hence high degrees even at moderate $\alpha$, $\delta$ and $\varepsilon$. Fang, Lin and Tong observe that the amplification never needed that: it only needs a polynomial agreeing with $\gamma' x$ on $I = [-\gamma'^{-1}, \gamma'^{-1}]$, and outside $I$ the value may be anything the norm constraint $|p(x)| \le 1$ on $[-1,1]$ allows. Dropping the requirement to vanish outside $I$ is the whole of the idea.

  The odd polynomial of degree $d$ is written in odd Chebyshev polynomials as $F(x) = \sum_{k=0}^{(d-1)/2} T_{2k+1}(x) c_k$, the interval $[-1,1]$ is discretized at $M$ grid points $x_j = -\cos\frac{j\pi}{M-1}$, and the coefficient matrix is $A_{jk} = T_{2k+1}(x_j)$. The coefficients then solve $\min_{\{c_k\}} \max_{x_j \in [0,\gamma'^{-1}]} |F(x_j) - (1-\delta)\gamma' x_j|$ subject to $|F(x_j)| \le c$ at every grid point, which is a convex problem. The norm constraint is imposed at the sampled points only, so $|F(x)| \le 1$ is relaxed to $|F(x_j)| \le c$ with $c$ chosen close to $1$, for instance $c = \max\{0.9999,\ 1 - 0.1\delta\}$, to absorb the slight overshoot possible between grid points. Because the program approximately solves a min-max problem, it is near-optimal in the $L^\infty$ sense by definition, in the asymptotic and the pre-asymptotic regime alike.

  The convex program is solved with CVX, Grant and Boyd's Matlab software for disciplined convex programming, version 2.1. The resulting Chebyshev coefficients are then fed to QSPPACK, at https://github.com/qsppack/QSPPACK, which finds the phase factors using the optimization-based method of Dong, Meng, Whaley and Lin. Both are named in Appendix C; the paper publishes no repository of its own for the solver.

  At $\gamma' = 5$, $\delta = 0.05$, $\varepsilon = 0.01$ the construction returns degree 21, against 2001 for the polynomial approximating $(1-\delta)\,\mathrm{XRect}(x)$ at the same parameters. Figure 2 draws both polynomials in panels (a) and (b) and their errors on $I$ in panels (c) and (d); panel (e) sweeps the degree at fixed $\gamma'$ and $\delta$ and reports the $L^\infty$ error $\max_{x \in I}|p(x) - (1-\delta)\,\mathrm{XRect}(x)|$, showing both converge exponentially in the degree while the convex-optimization rate is significantly faster, which the authors say reduces the number of queries to each $U_l$ by orders of magnitude. This is a classical numerical construction and no hardware or quantum simulator is named for it; it measures the polynomial, not the solver.

What it needs

1. Choose a time discretization or propagator approximation 6 methodsruns once per time step, with an amplification at every one

   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.

   The loop stays coherent: nothing is measured between turns. The preparation may still be reapplied every turn — what the loop never pays is a readout and a restart from classical data. The price is depth, and a success probability that multiplies down the chain. Fang, Lin and Tong's method is the repetition stated as the design: propagate one step at a time, and defeat the exponentially vanishing success probability by repeatedly invoking uniform singular value amplification. Nothing is measured between turns — the decay is a coherent one and it is bought back coherently — which is why the cost lands on the amplification ratio rather than on a shot count. The authors prove that linear dependence attains the query-complexity lower bound, so this repetition cannot be made cheaper in the worst case; it can only be avoided by not marching.

Slots it makes unnecessary

These do not fill the slot. They replace the span it belongs to, so this layer is not on their path at all.

• Quantum linear solve

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

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

• 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

• Time-marching based quantum solvers for time-dependent linear differential equations↗Di Fang, Lin Lin, Yu Tong · 2022theory yes · simulation yes · hardware no
• Quantum algorithm for linear non-unitary dynamics with near-optimal dependence on all parameters↗Dong An, Andrew M. Childs, Lin Lin · 2023theory yes · simulation yes · hardware no