Requires $H$ as a linear combination of efficiently implementable unitaries with a cheap PREPARE over the coefficients. The projection onto the wanted block succeeds only on the all-zeros flag, so amplification is part of the method rather than an afterthought, and the cost scales with the 1-norm of the coefficient vector rather than with $||H||$ — a decomposition with many comparable terms is expensive even for a benign $H$.

An access model for $H$ — a sum of efficiently exponentiable terms, sparse-access oracles, or a block-encoding — plus an evolution time $t$ and a target error $\varepsilon$.

1. Hamiltonian you can query → Block-encodingBlock-encode a matrix

   With $H = \sum_{\ell=1}^{L} \alpha_\ell H_\ell$ into unitary $H_\ell$, what this method needs from the slot is access not to $H$ but to the truncated Taylor terms. Indexing $j = (k,\ell_1,\ldots,\ell_k)$, the target unitaries are $V_j$, each of the form $(-i)^k H_{\ell_1}\cdots H_{\ell_k}$, and the access operator is abstracted as $\mathrm{select}(V)|j\rangle|\psi\rangle = |j\rangle V_j|\psi\rangle$. Concretely $\mathrm{select}(V)$ maps $|k\rangle|\ell_1\rangle\cdots|\ell_k\rangle|\psi\rangle$ to $|k\rangle|\ell_1\rangle\cdots|\ell_k\rangle(-i)^k H_{\ell_1}\cdots H_{\ell_k}|\psi\rangle$, and is built from $K$ controlled-$\mathrm{select}(H)$ operations, each mapping $|b\rangle|\ell\rangle|\psi\rangle$ to $|b\rangle|\ell\rangle(-iH_\ell)^b|\psi\rangle$; the $\kappa$th acts on the $\kappa$th qubit of $|k\rangle$, which is encoded in unary, and on the register holding $\ell_\kappa$, applying $-iH_{\ell_\kappa}$ to $|\psi\rangle$ only when $\kappa \le k$. approximation: each $H_\ell$ is applied exactly, by the mechanism the paper assumes for it, so the only approximation entering at this layer is the cut at order $K$ — the index register carries $k \le K$ and nothing beyond it.

   approximation

   Open this step

2. Block-encoding → Circuit for e^{-iHt}sum the truncated series, then amplify

   The evolution is cut into $r$ segments and each segment's $e^{-iHt/r}$ is replaced by its Taylor series truncated at order $K$, which expands into a sum of unitaries $\tilde{U} = \sum_j \beta_j V_j$. Berry, Childs, Cleve, Kothari and Somma implement that sum by preparing $B|0\rangle = (1/\sqrt{s}) \sum_j \sqrt{\beta_j}|j\rangle$ with $s = \sum_j \beta_j$ and sandwiching the selection: $W = (B^\dagger \otimes I)\,\mathrm{select}(V)\,(B \otimes I)$, whose all-zeros ancilla branch carries $\tilde{U}/s$. Segment length is then chosen to make the amplitude exactly right rather than merely large: with $T = (\alpha_1 + \cdots + \alpha_L)t$ and $T/r = \ln 2$ one has $s = 2$, and oblivious amplitude amplification turns the branch into certainty as $A = -WRW^\dagger RW$ with $R = I - 2P$. approximation: Truncating at $K$ makes $\tilde{U}$ close to unitary but not unitary, which is exactly what ordinary oblivious amplitude amplification needs, so the paper proves a robust version instead. Taking $K = O(\log(T/\varepsilon)/\log\log(T/\varepsilon))$ keeps $|s - 2| \le \varepsilon/r$ and the per-segment error $O(\varepsilon/r)$, so the whole evolution lands within $\varepsilon$. This is why the cost is paid in the $1$-norm of the decomposition rather than in $\|H\|$: $s$ is a sum of coefficient magnitudes, and $T$ is built from $\alpha_1 + \cdots + \alpha_L$.

   approximation

A circuit approximating $e^{-iHt}$ to within $\varepsilon$, with a stated query or gate count, an ancilla count, and the norm parameter — sparsity times $\lVert H\rVert_{\max}$, or the LCU 1-norm — that the cost is measured against.

• Prepare an input state

  Map $|0\ldots0⟩$ to a state whose amplitudes are proportional to a specified vector $b$, to within $\varepsilon$. The cost is set by which description of $b$ you hold, not by the algorithm that consumes it.

  What this method asks of the slot is one fixed state, not an arbitrary vector: $B|0\rangle = (1/\sqrt{s})\sum_j \sqrt{\beta_j}|j\rangle$ with $s := \sum_j \beta_j$, over the register indexing $j = (k,\ell_1,\ldots,\ell_k)$, whose amplitudes are the LCU coefficients $\beta_{(k,\ell_1,\ldots,\ell_k)} = \frac{(t/r)^k}{k!}\,\alpha_{\ell_1}\cdots\alpha_{\ell_k}$. Those factorize register by register, so $B$ is a tensor product over the $K+1$ registers: the order register goes from $0^K$ to a normalized superposition of the unary strings $1^k0^{K-k}$ in $O(K)$ gates, by a cascade of singly-controlled rotations; each of the $K$ index registers goes from $|0\rangle$ to the normalized $\sum_{\ell=1}^{L}\sqrt{\alpha_\ell}|\ell\rangle$ at a gate cost of $O(L)$. The inverse is wanted as often as the preparation: the paper prices one amplification step at two instances of $\mathrm{select}(V)$, one of $\mathrm{select}(V)^\dagger$, and three each of $B$ and $B^\dagger$ — so a single segment of evolution runs this slot three times forward and three times back.

• Amplify a success branch

  Take a routine whose output lands in the wanted subspace only with probability $a$, and raise that probability to near 1 using quadratically fewer repetitions than restarting the routine would need.

  What this method hands the slot is $W = (B^\dagger \otimes I)\,\mathrm{select}(V)\,(B \otimes I)$, whose all-zeros ancilla branch reads $PW|0\rangle|\psi\rangle = (1/s)|0\rangle\tilde{U}|\psi\rangle$ with $P := |0\rangle\langle 0| \otimes I$ and $\tilde{U}$ the order-$K$ truncation of one segment's evolution $U_r = e^{-iHt/r}$ — and nothing at all about $|\psi\rangle$. Hence the realization has to be the oblivious one: the ancilla reflection $I - 2P$ that $A$ interleaves with $W$ and $W^\dagger$ is the identity on the system register, so one fixed circuit serves every input, where ordinary amplitude amplification would reflect about $|\psi\rangle$ and so need it known and re-preparable. What returns is an operator identity rather than a raised probability: for exactly unitary $\tilde{U}$, $A|0\rangle|\psi\rangle = |0\rangle\tilde{U}|\psi\rangle$, unmeasured, once per segment and $r$ times over the evolution — once, because $s = 2$ leaves nothing for a second round. assumption: the robust form actually invoked returns $PA|0\rangle|\psi\rangle = |0\rangle\big(\tfrac{3}{s}\tilde{U} - \tfrac{4}{s^3}\tilde{U}\tilde{U}^\dagger\tilde{U}\big)|\psi\rangle$, and needs two closeness conditions at once, $|s - 2| = O(\delta)$ and $\|\tilde{U} - U_r\| = O(\delta)$, before $\|PA|0\rangle|\psi\rangle - |0\rangle U_r|\psi\rangle\| = O(\delta)$ follows; a normalization near $2$ without $\tilde{U}$ near $U_r$ buys nothing here.

  assumption

given  H = sum_j c_j V_j  with each V_j efficiently implementable and a cheap
       PREPARE over the coefficients; evolution time t split into short
       segments

for each segment:
    truncate the Taylor series of exp(-iH t_seg) and collect its terms as a
        linear combination of unitaries

    PREPARE            # load the coefficients into the ancilla register
    SELECT             # apply the terms, controlled on that register
    PREPARE^-1         # unprepare

    amplify the all-zeros ancilla flag
        # the projection onto the wanted block succeeds only on that flag,
        # so this is part of the method rather than an afterthought

# the cost scales with the 1-norm of the coefficient vector rather than with
# ||H||: a decomposition with many comparable terms is expensive even for a
# benign H

Berry, Childs, Cleve, Kothari and Somma: the cost of the algorithm depends logarithmically on the inverse of the desired precision, which the authors state is optimal. Berry, Childs and Kothari, using a linear combination of quantum walk steps with coefficients given by Bessel functions, report a complexity in queries and 2-qubit gates that is "logarithmic in the inverse error, and nearly linear in the product $\tau$ of the evolution time, the sparsity, and the magnitude of the largest entry of the Hamiltonian".

• QuDiffEq.jl's `taylorcircuit`/`taylorsolve` construction

  • Simulating Hamiltonian dynamics with a truncated Taylor series

    Dominic W. Berry, Andrew M. Childs, Richard Cleve, Robin Kothari, Rolando D. Somma · 2014

  About

  QuDiffEq.jl is a Julia package built on the Yao.jl circuit framework; its README describes it as "Quantum algorithms for solving differential equations" and states "This project is part of Julia's Season of Contribution 2019". Its documentation page for this feature names this record's own source paper directly: "Taylor truncation based Hamiltonian simulation (https://arxiv.org/abs/1412.4687) has many clear advantages. It has better complexity dependence on the precision and allows a greater range of Hamiltonians to be simulated." That page states the package "provides circuits for five kinds of problems": "Unitary Taylor simulation", "Non-unitary Taylor simulation", "Unitary QuLDE Problem", "Non-unitary QuLDEProblem", and "Solution by linearising a non-linear differential equation" — the first two are the plain e^{Mt} construction this record's hops describe; the other three add a differential-equation source term the record does not cover and this entry does not draw on. The docs describe both branches in the record's own PREPARE/SELECT/PREPARE-dagger terms without using those names: for the case where the generator itself is unitary, "the VS1 block acts on the ancilla register to generate an appropriate superpostion ... Multiplication of M^j block is controlled by |j> in the ancilla register. VS1', the adjoint of VS1, un-computes the ancilla registers. The desired result is obtained when the resulting state is projected onto |0> ancilla state."; for the non-unitary case, "M is expressed as a linear combination of four (at most) unitary i.e. M = sum_i (1/2) F_i", selected by a second per-order register, and again "un-computed and measured in the zero ancilla state".

  Methods

  The package implements the PREPARE-SELECT-PREPARE-dagger sandwich for a single time window, not the record's full construction: there is no loop splitting the evolution into r segments, and no oblivious-amplitude-amplification stage — `taylorsolve` postselects the ancilla on the all-zeros outcome once and returns the sub-normalized branch together with a classical rescaling factor, rather than boosting the branch to certainty. `TaylorParam`'s constructor branches on `u = isunitary(H/opn)` (opn = opnorm(H)); for `u = false` it sets `c = 2; rs = k; l = 2`, and this is the branch that matches the record's block-encode-matrix hop closely: `calc_vs1` writes the PREPARE amplitudes into rows `VS1[(2^k - 2^(k-j) + 1),1]` for `j` in `0:k` (a k-qubit unary encoding of the order, since `2^k - 2^(k-j)` is the bit pattern `1^j 0^(k-j)`), each set to `sqrt(C(j, x, opn, t, 2))/C_tilda` with `C(m,x,opn,t,c) = norm(x)*(opn*t*c)^m/factorial(m)`, then extends this partial first column to a full unitary via `VS1 = -1*qr(VS1).Q` — applying VS1 to the all-zeros ancilla yields its first column, i.e. exactly the amplitudes just written. `unitary_decompose(H)` supplies the terms SELECT chooses among: it splits `Mu = H/opnorm(H)` into Hermitian and anti-Hermitian parts `B1, B2` and returns four unitaries `F[1] = B1 + i*sqrt(I - B1*B1)`, `F[2] = B1 - i*sqrt(I - B1*B1)`, `F[3] = i*B2 - sqrt(I - B2*B2)`, `F[4] = i*B2 + sqrt(I - B2*B2)` — a fixed decomposition into at most 4 unitaries that the package derives automatically from H, not an arbitrary L-term decomposition supplied by the caller as the record's hop assumes. `circuit_intermediate` (the `u = false` method) then loops `for i in 1:k, for j in 0:2^l-1` and appends a Yao `control` gate applying `F[j+1]` to the target register, controlled on order-register qubit `i` together with a 2-qubit sub-register holding `j` — the per-slot select this record's block-encode-matrix hop describes, with L fixed at (up to) 4 rather than general. PREPARE for this branch is not `VS1` alone: `taylorsolve`'s `rs == k` path also builds `VT = calc_vt(CPType)`, a fixed 4-by-4 unitary whose first column is uniform (`VT[:,1] = 0.5*ones(4)`, then extended to a full unitary by the same `-1*qr(...).Q` construction as `VS1`) — an unconditional preparation of the equal superposition over the (at most) 4 `F_i` terms, applied once per order-index qubit. It then calls the four-argument `taylorcircuit(n, blk, VS1, VT)`, which chains `circuit_ends(n, blk, VS1, VT)` (PREPARE: `VS1` on the order register composed with, for each `i` in `1:k`, an unconditional `lc(...)` placing `VT` on that index's own 2-qubit term-selection sub-register), `circuit_intermediate` (SELECT), and `circuit_ends(n, blk, VS1', VT')` (PREPARE-dagger) in that order — matching the package's own description, "the second register is superposed by `VT`, where each new state corresponds to a `F_i`". The three-argument `taylorcircuit(n, blk, VS1)`, which chains the VT-free `circuit_ends(n, blk, VS1)` alone with `circuit_intermediate` and `circuit_ends(n, blk, VS1')`, is dispatched only for `TaylorParam{CPType, true}` — the `u = true` (unitary-H) branch — and is not what runs for the `u = false` branch this paragraph describes. `taylorsolve`'s docstring reads "Simulates a Hamiltonian using the Taylor truncation method. Returns the state register and inverse probability of finding it."; its body runs `r = apply!(inreg,cir) |> focus!(1:n - nbit...,) |> select!(0)` and returns `r, blk.N` — postselecting the ancilla on 0 once, with no amplification, and handing back `N = blk.N` (the classically-summed Taylor-magnitude normalization) for the caller to rescale by afterward. The `u = true` branch (H itself already proportional to a unitary) uses a different, binary order-register encoding, `rs = log2i(k+1)`, and repeatedly multiplies the same H rather than selecting among several terms; it is a special case the record's general L-term construction does not describe and this entry does not claim matches it structurally.

  Data

  No molecule or physical model: `test/TaylorTrunc_tests.jl` builds one random 1-qubit unitary `Au = rand_unitary(2)` and one random 1-qubit non-unitary matrix `An = rand(ComplexF64, 2, 2)`, plus normalized random vectors `b, x`, with `Random.seed!(2)`, Taylor order `k = 3`, and evolution time `tspan = (0.0, 0.1)`. Each is wrapped in a `QuLDEProblem(A, b, x, tspan)` and then passed to the plain Hamiltonian-simulation entry point `taylorsolve(qprob.A, qprob.u0, k, tspan[2])`, not to a QuLDE-specific solver.

  Code

  `src/TaylorTrunc.jl` in https://github.com/QuantumBFS/QuDiffEq.jl — Julia, MIT License (`LICENSE`: "Copyright (c) 2019 Divyanshu Gupta"), 334 lines on `master` as read on 2026-08-26 at commit 83dd4b4a, which is also `master`'s current HEAD and the last commit to touch this file (2024-08-24). `Project.toml` names the package `QuDiffEq`, version `0.2.0`, depending on `Yao = "0.9"`. The relevant exports from this file are `TaylorParam`, `taylorcircuit`, `taylorsolve`, `circuit_ends`, `circuit_intermediate`; `calc_vs1`, `calc_vs2`, `calc_vt` and `unitary_decompose` are unexported helpers in the same file. Documentation source is `docs/src/man/taylor.md`. Tests are `test/TaylorTrunc_tests.jl`.

  Results

  The only check is a numerical comparison against exact matrix exponentiation, run for both branches. `test/TaylorTrunc_tests.jl` computes `out = N*vec(state(r))` from `r, N = taylorsolve(qprob.A, qprob.u0, k, tspan[2])` and compares it to `r_out = exp(qprob.A*tspan[2])*qprob.u0`, asserting `isapprox.(r_out, out, atol = 1e-3) |> all` — once for the unitary instance `Au` and once for the non-unitary instance `An`, both at `k = 3`, `t = 0.1`, on a single qubit. No gate count, T-count, circuit depth or hardware/simulator backend is reported anywhere in the package for this construction; the only named computational engine is Yao.jl's own state-vector simulation via `apply!`.

• Linear combination of unitaries · Qiskit

  From the repository — run, not written up from a paper · native

  About

  A block-encoding primitive that turns a weighted sum of unitary operations into one larger unitary circuit.

  Methods

  None found yet.

  Data

  None found yet.

  Code

  Qiskit

  Results

  Two-term PREPARE–SELECT–unprepare circuit checked on a one-qubit target · success branch and coefficient normalization recorded

• Linear combination of unitaries

  A block-encoding primitive that turns a weighted sum of unitary operations into one larger unitary circuit.