Truncated Dyson series of the propagator

What it fills

• Choose a time discretization or propagator approximation

  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.

When it applies

Theorem 4.1 is stated for $\dot{x}(t) = A(t)x(t) + b(t)$ with $A(t)$ of non-positive logarithmic norm, the equation's parameters provided through unitaries $U_A$, $U_b$, $U_x$ with known normalisations $\lambda_A$, $\lambda_b$, $\lambda_x$. No smoothness condition is required: the oracle counts are independent of derivatives of the parameters.

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) are theorem and proof throughout: no figure, no table, no computed value, and no hardware or simulator named. What stands where a worked instance would be is the symbolic three-step block matrix of Eq. (19), whose entries are symbols — no generator, no inhomogeneity and no dimension is ever chosen.

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

truncate the Dyson series for the propagator at order k:
    sum_{j=0..k}  int ... int  A(t_1) ... A(t_j)  dt_j ... dt_1
                  over  T >= t_1 >= ... >= t_j >= 0
    # the expansion that stands in for the propagator once the generator
    # varies with time

encode the surviving terms as rows of one system of linear equations:
    for j = 1 ... k:  rows for the j-th time-ordered integral, discretized
    plus the rows tying the truncated propagator to the solution at T

hand the system to the layer below
    # no smoothness condition is required: the oracle counts are independent
    # of derivatives of the parameters

Cost, as the source states it

Berry and Costa's Theorem 4.1 gives, for $\dot{x}(t) = A(t)x(t) + b(t)$ with $A(t)$ of non-positive logarithmic norm, an average of $O(R \lambda_A T \log(1/\varepsilon))$ calls to the state-preparation oracles and $O(R \lambda_A T \log(1/\varepsilon) \log(\lambda_{Ax} T/\varepsilon))$ calls to the matrix oracle — $T$ the evolution time, $\varepsilon$ the allowed error, $\lambda_A$ the matrix oracle's normalisation, $\lambda_{Ax} = \max(\lambda_A, b_{max}/x_{max})$, and $R$ a rescaling constant the theorem bounds explicitly, growing when the solution decays. Gate counts depend on the parameters' first derivatives only through a logarithm, and on no higher derivative. That is a statement about the full algorithm, not a standalone cost for the discretization.

Implementations

• Block encoding of the truncated propagator by sorted time registers (section 3.1)

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

  How the blocks $V_m = W_K(m\Delta t,(m-1)\Delta t)$ are actually produced, which the paper flags twice as the hard part and as the point of departure from the Taylor route: the blocks are multiple integrals over time, so the block encoding is not trivial, and this is a major departure from the method for time-independent differential equations of Berry, Childs, Ostrander and Wang, where it was possible to encode the terms of the sum via extra lines in the block matrix.

  The input is the intermediate matrix of Eq. (55), block diagonal with $A(\delta t)$, $A(\Delta t + \delta t)$, $A(2\Delta t + \delta t)$ and so on down the evolving rows and zeros below, where $\delta t$ indexes an offset between $0$ and $\Delta t$; it is obtained from the block encoding of $A(t)$ by using the qubits encoding $\delta t$ together with the qubits encoding the row of the block matrix as the time input, taking the number of intervals per $\Delta t$ as a power of two, and flipping one ancilla to eliminate the lower-right zeros, most easily when $M$ is also a power of two. That matrix is fed to the truncated series of Eq. (56) by the procedure of Kieferova, Scherer and Berry or of Low and Wiebe: $K$ time registers are prepared in equal superposition and sorted by a sorting network, at an optimal $O(K \log K)$ steps on registers of $\log M$ bits, so $O(K \log K \log M)$ gates, and the block encoding of $A(\delta t)$ is used $K$ times, which is $K$ calls to the block encoding of $A(t)$. The register holding $k$ is prepared in unary with controlled rotations; the truncated series needs precision $O(\varepsilon/r)$ so each of the $K$ rotations needs $O(\varepsilon/(Kr))$, giving $O(K\log K \log(r/\varepsilon))$. Two alternatives are named and neither is adopted: the inequality-testing preparation of the time registers costs $O(K \log M)$ from $K-1$ tests on registers of $\log M$ bits, and the inequality-testing preparation of $k$ due to Su, Berry, Wiebe, Rubin and Babbush needs $O(K \log K)$ bits for similar order scaling, while an equal superposition over $k$ alone would be $O(\log K)$ — but in practice the sorting approach is preferred because it gives a better constant factor. With $\lambda_A\Delta t \le 1$ the $\lambda$-value of the result is $\sum_{k=0}^{K}(\lambda_A\Delta t)^k/k! < e^{\lambda_A\Delta t} \le e$, checked twice over, once against the normalisation factor of Kieferova, Scherer and Berry's Eq. (55) and once by the paper's own rule that a block-encoded polynomial takes the same arithmetic on its $\lambda$-values as on its operators. $M$ is then fixed by requiring the discretisation error of Eq. (102), $O((T/r)^2\max_t\lVert A'(t)\rVert/M)$, to stay below $\varepsilon/r$ per segment.

  No dataset. The only input is the block encoding of $A(t)$ accepting the time register as a quantum input, with known normalisation $\lambda_A$; the paper assumes it rather than constructing it, saying explicitly that it does not consider how it would specifically be done for a particular encoding such as sparse matrix oracles. The bound $\max_t \lVert A'(t) \rVert$ enters only through the choice of $M$ and hence only through the constant $D$.

  No numerics. What the construction yields is counted rather than run: $K$ calls to the block encoding of $A(t)$ per block, which is where Theorem 4.1's factor $\log(\lambda_{Ax}T/\varepsilon)$ in the $U_A$ count comes from, and $O(K\log K\log M)$ gates, whose $\log M$ becomes the $\log(TD/(\lambda_A\varepsilon))$ term in the theorem's gate count. Nothing is executed and no simulator or device is named.

• Block encoding of the driven term and preparation of the right-hand side (section 3.1)

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

  The other half of the same discretisation: the blocks $v_m = v_K(m\Delta t,(m-1)\Delta t)$ of Eq. (59), which are the same time-ordered integrals with the innermost $A(t_k)$ replaced by $b(t_k)$. The paper credits one of its improvements over the spectral method of Childs and Liu directly to this shape — the complexity in calls to block encodings of $b(t)$ is only the $\log(1/\varepsilon)$ factor from the linear equation solver, because $b(t)$ is only needed once in the Dyson series solution.

  The block encoding is applied in almost an identical way to the one for the propagator, except that the innermost block encoding of $A(t_k)$ is replaced by the preparation of $\lvert b(t_k) \rangle$: the same preparation and sorting of time registers at $O(K\log K\log M)$ gates, $K-1$ calls to the block encoding of $A(t)$ and a single call to that of $b(t)$. With $\lambda_A\Delta t \le 1$ the $\lambda$-value is $(\lambda_b/\lambda_A)\sum_{k=1}^{K}(\lambda_A\Delta t)^k/k! < (\lambda_b/\lambda_A)(e^{\lambda_A\Delta t}-1) \le (e-1)\lambda_b\Delta t$. The whole right-hand side is then assembled by Eqs. (61)-(62): a state on the time register weighting $\lambda_x\lvert 0 \rangle$ against $(e-1)\lambda_b\Delta t\sum_{m=1}^{r}\lvert m \rangle$, then the preparation of $\lvert x(t_0) \rangle$ controlled on $\lvert 0 \rangle$ and of $\lvert v_m \rangle$ controlled on $\lvert m \rangle$. The time register itself is built by a rotation on one ancilla to set the weight between $m = 0$ and $m = 1,\ldots,r$, then Hadamards on the $\log r$ time qubits when $r$ is a power of two (or the $O(\log r)$ non-power-of-two preparation of Sanders and co-authors otherwise), CNOTs from the ancilla onto the remaining time qubits, and a modular addition of $1$ with the ancilla as most significant bit, so that the ancilla being $\lvert 1 \rangle$ collapses to all zeros and being $\lvert 0 \rangle$ gives $1$ through $r$; this costs $O(\log(1/\varepsilon))$ for the rotation and $O(\log r)$ for the rest, and is dwarfed by the rest of the procedure. The success amplitude of Eq. (63) is bounded below by $\min_m\lVert v_m \rVert/((e-1)\lambda_b\Delta t)$ with $\lVert v_m \rVert \ge \lVert v(m\Delta t,(m-1)\Delta t) \rVert - \varepsilon x_{\max}/r$; because the right-hand side is prepared twice inside the walk step of the linear equation solver, amplitude amplification first boosts the preparation to $\Theta(1)$ at the count of Eq. (65), and the amplitude, usually unknown, is estimated at the beginning at logarithmic overhead — far cheaper than the solve, which carries an extra factor of the condition number. The discretisation error of these integrals, Eq. (104), is $O((T/r)^2\max_t\lVert b'(t)\rVert/M + (T/r)^3 b_{\max}\max_t\lVert A'(t)\rVert/M)$, the first term from $k=1$ and the second from $k=2$, and it is the only place a derivative of $b$ enters the algorithm at all.

  No dataset. The inputs are the block encoding of $b(t)$ taking the time register as a quantum input with a time-independent $\lambda_b$, the block encoding of $A(t)$, and the unitary $U_x$ preparing $\lvert x_0 \rangle/\lambda_x$. The quantity $\min_m \lVert v(m\Delta t,(m-1)\Delta t) \rVert$ is treated as a given bound rather than computed; the paper remarks that if the norm of $b$ does not vary too much and there are no cancellations in the integral this is proportional to $b_{\max}\Delta t$, so that with $\lambda_b$ comparable to $b_{\max}$ the amplitude would be at least a constant, and says it keeps the full expression only for generality.

  No numerics. The construction contributes the $O(R\lambda_A T\log(1/\varepsilon))$ count of calls to $U_b$ and $U_x$ in Theorem 4.1 and the second of the two factors in that theorem's constant $R$, and it is never executed: the paper reports no numerical experiment and names no hardware and no simulator.

• Collapse to a truncated Taylor block when the generator is constant (sections 1.3 and 3.2)

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

  What this discretisation becomes when $A$ does not vary. The paper reaches it not by a separate derivation but by substituting a constant $A$ into the solution for the time-varying $A$, and frames it as applying a similar principle to the Taylor series solution of Berry, Childs, Ostrander and Wang, except that it will encode into the block matrix in a simpler way.

  The truncated propagator becomes $W_K = \sum_{k=0}^{K}(A\Delta t)^k/k!$ and the driven block becomes $v_K = \sum_{k=1}^{K}A^{k-1}\Delta t^k b/k!$, both independent of $m$, with the truncation error of Eq. (18), $O((\lVert A \rVert\Delta t)^{K+1}\lVert x(t_0) \rVert/(K+1)! + \lVert A \rVert^K\Delta t^{K+1}\lVert b \rVert/(K+1)!)$. The block encoding proceeds exactly as in the time-dependent case except that there is no need for the time integrals, so no time registers are prepared and no sorting network is used, and only the register holding $k$ has to be prepared for the two series; with $\lambda_A\Delta t \le 1$ the Taylor series is block encoded with a factor of at least $1/e$. The cancellation that the time-dependent driven block has to guard against cannot occur here: Eq. (69) bounds $\lVert v \rVert \ge (3-e)\lVert b \rVert\Delta t$ under $\lVert A \rVert\Delta t \le \lambda_A\Delta t \le 1$.

  No dataset. The inputs lose their time register entirely: an ordinary block encoding of $A$ with normalisation $\lambda_A$, and unitary preparations of $\lvert b \rangle$ and $\lvert x_0 \rangle$, with $\lambda_b$ now simply $\lVert b \rVert$.

  No numerics. The saving is stated rather than measured: with no arithmetic on time registers, the $\log(TD/(\lambda\varepsilon))$ term disappears from the gate count of Theorem 4.2, and with the state-preparation amplitude now at least a constant, that factor disappears from the theorem's constant $R$ as well.

What it needs

Nothing below this — it bottoms out here.

Other ways to fill the same slot

Different approaches

• Forward (explicit) Euler

  First-order explicit stepping, $u_{k+1} = (I + hA)u_k + h b_k$, assembled into a banded all-at-once linear system. Liu et al. use it inside the Carleman route because its structure is simple enough to bound explicitly.

• Backward (implicit) Euler

  First-order implicit stepping: each step solves $(I - hA)u_{k+1} = u_k + h b_{k+1}$. $A$-stability is the classical reason to prefer it for stiff generators, since it removes the explicit method's step-size restriction.

• Trapezoidal rule (Crank-Nicolson)

  Second-order implicit stepping that averages the generator at the two ends of each step, and is $A$-stable. As a rational approximation of $e^{hA}$ it is the $(1,1)$ diagonal Padé approximant.

• Truncated Taylor series of the propagator

  Rather than approximating the derivative, approximate the propagator $e^{hA}$ itself by $k$ Taylor terms and encode those terms as extra rows of a sparse linear system. Truncation error falls factorially in $k$, so accuracy is bought by adding rows rather than by shrinking $h$ and lengthening the system.

• Chebyshev pseudospectral collocation

  A **global** discretization, and the only one on this slot. Write each component of the solution as a truncated Chebyshev series $x_i(t) = \sum_{k=0}^{n} c_{i,k}T_k(t)$ with the coefficients undetermined, then fix them by demanding that the differential equation hold exactly at the $n+1$ Chebyshev-Gauss-Lobatto nodes $t_l = \cos(l\pi/n)$. What comes out is one linear system in the coefficients. Childs and Liu put the difference from every stepping scheme plainly: "instead of locally approximating the ODE at discretized times, these linear equations use the behavior of the differential equations at the $n+1$ times $\{t_l\}$ to capture their behavior over the entire interval $[-1,1]$".

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