Dyson 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 time-dependent linear differential equations, with a simplified approach in the time-independent case. It reduces to a quantum linear solve. An, Childs and Lin's comparison table places truncated Dyson in the family whose initial-state-preparation cost grows with the evolution, not in the bypass family.

Requires

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

Example

Berry and Costa (arXiv:2212.03544) carry no numerical work at all — no figure, no table, no computed value from the first section through the appendix, and no hardware or simulator named. An, Childs and Lin, cited beside it, do compute, but every number there measures the LCHS kernel, which is the competing route rather than this one.

given  A(t) with non-positive logarithmic norm, b(t), x_0, horizon T, budget e
       the parameters provided through the unitaries U_A, U_b, U_x
       with known normalisations lambda_A, lambda_b, lambda_x
       x_max >= max_t ||x(t)||

split [0,T] into r steps of length dt = T / ceil(lambda_A T)
    # dt taken as large as lambda_A dt <= 1 allows, which is what makes r = lambda_A T

# encode the Dyson series in a system of linear equations
for each step m = 1 ... r:
    V_m = W_K(m dt, (m-1) dt)    # the Dyson series for the propagator, truncated at K
    v_m = v_K(m dt, (m-1) dt)    # the same integrals with b(t_k) in place of A(t_k)
    # forming these two is this route's time discretization

assemble the lower bidiagonal block system  AA X = B  over R = 2r rows:
    row 0       :  x~(0)                            = x_0     # the initial condition
    row m <= r  :  x~(m dt) - V_m x~((m-1) dt)      = v_m     # one row per time step
    row m > r   :  x~(m dt) -     x~((m-1) dt)      = 0       # solution held constant
    # the trailing rows do not evolve anything; they are there to boost the success
    # probability of the final amplitude amplification
    # AA is the paper's block matrix; kappa_AA = ||AA|| ||AA^-1|| = O(R) = O(lambda_A T)

prepare B, and amplitude-amplify its state preparation to amplitude Theta(1)
    # B is prepared twice inside the solver's walk step; without this boost the
    # success amplitude would enter squared

solve  AA X = B  with the optimal quantum linear equation solver
    # O(kappa_AA log(1/e)) calls to the block encodings of AA and B

amplitude-amplify the component of X over rows r+1 ... R
    # O(x_max / ||x(T)||) steps on average

return the state approximating |x(T)>

Cost, as the source states it

Logarithmic dependence of the complexity on the error and derivative, with the usual exponential improvement over classical approaches in the scaling with the dimension, the solution being encoded in the amplitudes of a quantum state.

Implementations

• Bidiagonal block system with a Dyson-series block encoding (sections 2.1, 3.1 and 4.1)

  • Quantum algorithm for time-dependent differential equations using Dyson series

  The construction Theorem 4.1 is proved for, and the paper says plainly why it is not the earlier one: "We use a block matrix similar to [3], but we do not use extra lines of the block matrix to implement terms in the series as in that work. Instead we construct this matrix via a block encoding using a Dyson series, in an analogous way as the block encodings in [12, 15]." The conclusion names the payoff — this approach makes the analysis of the complexity simpler than in Berry, Childs, Ostrander and Wang, because it is not necessary to account for the extra lines in the encoding. The problem it is built for is $\dot{x}(t) = A(t)x(t) + b(t)$ with $x(0) = x_0$ and $A(t)$ an $N \times N$ coefficient matrix of non-positive logarithmic norm, the deliverable being a state $\hat{x}$ with $d_{BW}(\hat{x}, \lvert x(T) \rangle \langle x(T) \rvert) \le \varepsilon x_{\max}$ in the Bures-Wasserstein metric, which is adopted so that an approximate state that is not pure can still be scored and which reduces to the 2-norm distance for pure states.

  The block matrix $\mathcal{A}$ of Eq. (26) is never held as a matrix. Eq. (54) splits it as the identity minus the product of an increment on the time register with a block-diagonal matrix carrying the $V_m$. The increment is block encoded by incrementing the register non-modularly and using the carry qubit, so that the projection onto $\lvert 0 \rangle$ deletes the top row. The $V_m$ diagonal comes from the intermediate matrix of Eq. (55), which holds $A(\delta t), A(\Delta t + \delta t), A(2\Delta t + \delta t), \ldots$ on the evolving rows and zeros on the held rows; it is built from $U_A$ by feeding the qubits encoding the offset $\delta t$ together with the qubits encoding the row of the block matrix as the time input, with the number of subintervals per $\Delta t$ and $M$ both taken as powers of two so that a single ancilla qubit flags the zero block. That matrix is fed to the truncated Dyson series of Eq. (56), block encoded as in the Hamiltonian-simulation constructions of Kieferova, Scherer and Berry and of Low and Wiebe. With $\lambda_A \Delta t \le 1$ that series has $\lambda$-value $\sum_{k=0}^{K}(\lambda_A\Delta t)^k/k! < e^{\lambda_A \Delta t} \le e$, and the whole of $\mathcal{A}$ no larger than $e+1$. The right-hand side $\mathcal{B}$ is $x_0$ from $U_x$ together with the $v_m$ of Eq. (59), whose block encoding is the same machinery with the innermost $A(t_k)$ replaced by the preparation of $\lvert b(t_k) \rangle$, at $\lambda$-value below $(e-1)\lambda_b\Delta t$; the composite state of Eqs. (61)-(62) is prepared and then amplitude-amplified to $\Theta(1)$ before the solve, because $\mathcal{B}$ is prepared twice inside the solver's walk step and would otherwise contribute its amplitude squared. The conditioning is read off the bidiagonal structure: $\lVert \mathcal{A} \rVert \le 1 + \max_m \lVert V_m \rVert = O(1)$, and the explicit inverse of a lower bidiagonal matrix, Eqs. (30)-(31), expanded as a sum over its $R$ diagonals gives $\lVert \mathcal{A}^{-1} \rVert \le R(1+\varepsilon) = O(R)$, so $\kappa_{\mathcal{A}} = O(R)$, and with an equal number of evolving and held steps, $R = 2r$ with $r = \lambda_A T$, this is $O(\lambda_A T)$. The system is solved by the optimal quantum linear equation solver of Costa, An, Sanders, Su, Babbush and Berry at $O(\kappa_{\mathcal{A}} \log(1/\varepsilon_{\mathrm{QLSP}}))$ calls, with $\varepsilon_{\mathrm{QLSP}}$ taken proportional to $\varepsilon$ because $R = 2r$; the measurement that solver uses to confirm its filtering is omitted so that the final amplitude amplification over rows $r+1$ to $R$ can be performed jointly with the filtering, which the paper notes introduces no further factors.

  No dataset. The inputs are the three block encodings of Eqs. (44)-(46): $U_A$ with the time register as a quantum input, giving $A(t)\lvert n_t \rangle/\lambda_A$; $U_b$ likewise giving $\lvert n_t \rangle \lvert b(t) \rangle_s/\lambda_b$, defined as a block encoding rather than a unitary preparation precisely so that $\lambda_b$ can be taken independent of time even though $\lVert b(t) \rVert$ varies; and $U_x$ preparing $\lvert x_0 \rangle_s/\lambda_x$. The normalisations $\lambda_A$, $\lambda_b$, $\lambda_x$ are assumed known, and so are the constants $x_{\max} \ge \max_t \lVert x(t) \rVert$, $b_{\max} \ge \max_t \lVert b(t) \rVert$, $D$ of Eq. (88) and $R$ of Eq. (87). The paper's reason for taking these as given rather than computed is that it may be too difficult to determine the maxima and minima required exactly.

  The construction is carried as far as counted complexity and no further. Theorem 4.1 gives an average of $O(R\lambda_A T \log(1/\varepsilon))$ calls to $U_b$ and $U_x$, $O(R\lambda_A T \log(1/\varepsilon)\log(\lambda_{Ax}T/\varepsilon))$ calls to $U_A$, and $O(R\lambda_A T \log(1/\varepsilon)\log(\lambda_{Ax}T/\varepsilon)[\log(TD/(\lambda_A\varepsilon)) + \log(\lambda_A T/\varepsilon)])$ additional gates, where $\lambda_{Ax} = \max(\lambda_A, b_{\max}/x_{\max})$. The gate count is stated for a fixed gate set such as Toffoli or $T$ plus Clifford rather than arbitrary-precision rotations, which the paper equates up to a constant with the usual non-Clifford count. The two factors of $R$ are the two amplitude amplifications, $O(x_{\max}/\lVert x(T) \rVert)$ for extracting the final time and Eq. (95) for the preparation of $\mathcal{B}$, and the paper expects $R = O(1)$ in practice when the solution does not decay significantly and $b(t)$ does not vary so as to cancel. Nothing here is executed: the paper reports no numerical experiment, and names no hardware and no simulator.

• Simplified encoding for a time-independent generator (sections 1.3, 2.2, 3.2 and 4.2)

  • Quantum algorithm for time-dependent differential equations using Dyson series

  The same encoding when the coefficients do not vary — the abstract's "simplified approach in the case of time-independent differential equations", which the conclusion calls a significantly simplified result for the complexity in the time-independent case. The problem is $\dot{x}(t) = Ax(t) + b$ with $x(0) = x_0$ and $A$ of non-positive logarithmic norm. It is deliberately not the encoding Berry, Childs, Ostrander and Wang gave for that same problem: there the orders of the Taylor sum occupied extra lines of the block matrix, and here, as the conclusion puts it, the block matrix has each successive line encoding a new time step, with the whole sum inside one block encoding. The paper reaches it by substituting a constant $A$ into the solution for the time-varying case.

  Every $V_m$ collapses to the single $V = \sum_{k=0}^{K}(A\Delta t)^k/k!$ of Eq. (40) and every $v_m$ to $v = \sum_{k=1}^{K}A^{k-1}\Delta t^k b/k!$ of Eq. (41), so the block matrix of Eq. (39) carries one repeated subdiagonal block. The oracles lose their time register: $\langle 0 \rvert_A U_A \lvert 0 \rangle_A = A/\lambda_A$, $U_b \lvert 0 \rangle = \lvert b \rangle/\lambda_b$ with $\lambda_b = \lVert b \rVert$, and $U_x \lvert 0 \rangle = \lvert x_0 \rangle/\lambda_x$, all of them ordinary unitaries on one register. The block encoding proceeds through the same Eq. (54) split except that there is no need for the time integrals, so only the register holding $k$ has to be prepared, and with $\lambda_A\Delta t \le 1$ the Taylor series is block encoded with a factor of at least $1/e$. The conditioning argument is the time-dependent one word for word, $\lVert \mathcal{A} \rVert \le 1 + \lVert V \rVert = O(1)$ and $\lVert \mathcal{A}^{-1} \rVert = O(R)$, so $\kappa_{\mathcal{A}} = O(R)$ again. One thing genuinely improves rather than merely simplifying: the cancellations in $v_m$ that force the time-dependent case to carry an amplitude-amplification factor cannot occur, because Eq. (69) gives $\lVert v \rVert > \lVert b \rVert \Delta t\,(1 - [\exp(\lVert A \rVert \Delta t) - (1 + \lVert A \rVert \Delta t)]/(\lVert A \rVert \Delta t)) \ge (3-e)\lVert b \rVert\Delta t$ under $\lVert A \rVert \Delta t \le \lambda_A \Delta t \le 1$, putting the state-preparation success amplitude at no less than the constant $(3-e)/(e-1)$, which therefore drops out of the scaling.

  No dataset. The inputs are the three time-independent unitaries of Eqs. (110)-(112) with known $\lambda_A$, $\lambda_b = \lVert b \rVert$ and $\lambda_x$, together with $x_{\max} \ge \max_t \lVert x(t) \rVert$ and the constant $R \ge x_{\max}/\lVert x(T) \rVert$. The bound $D$ that the time-dependent theorem needs does not appear, because no derivative of the parameters enters.

  Theorem 4.2 gives an average of $O(R\lambda T \log(1/\varepsilon))$ calls to $U_b$ and $U_x$, $O(R\lambda T \log(1/\varepsilon)\log(\lambda_{Ax}T/\varepsilon))$ calls to $U_A$, and $O(R\lambda T \log(1/\varepsilon)\log(\lambda_{Ax}T/\varepsilon)\log(\lambda_A T/\varepsilon))$ additional gates, with $\lambda_{Ax} = \max(\lambda_A, \lVert b \rVert/x_{\max})$ and $\lambda$ standing where the time-dependent theorem writes $\lambda_A$. Two things shrink against Theorem 4.1 and the proof says which: the state preparation now succeeds with at least constant probability so that factor leaves $R$, and no arithmetic on time registers is needed so the $\log(TD/(\lambda\varepsilon))$ term leaves the gate count, though the rotations preparing the $k$ register keep the $\log(\lambda T/\varepsilon)$ factor. As with the time-dependent case, this is a specified circuit and not an executed one: no numerical experiment, no hardware and no simulator appears in the paper.

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

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

• 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 algorithm for time-dependent differential equations using Dyson series↗Dominic W. Berry, Pedro C. S. Costa · 2022theory yes · simulation no · 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