Sign outOpen workspaceSign in

MethodLayer 1

Truncated-Taylor LCU simulation

Truncate the Taylor series of e−iHte^{-iHt} over short segments and implement the truncated sum as a linear combination of unitaries — PREPARE loads the coefficients, SELECT applies the terms, PREPARE†† unprepares — with the all-zeros ancilla flag amplified.

Takes

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

Returns

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

Same contract as the slot it fills.

This one, drawn

Truncated-Taylor LCU simulationProduct-formula (Trotter-Suzuki) simulationTruncated-Taylor LCU simulation · 2 parts inside, openTruncated-Taylor LCU simulation · 2 parts inside, openTruncated-Taylor LCU simulation · 2 parts inside, open — click the line to close itBlock-encode a matrix · 4 ways throughBlock-encode a matrix · 4 ways through — click the line to open it hereTruncated-Taylor LCU simulationQubitization walk simulation · 2 parts insideQubitization walk simulation · 2 parts inside — click the line to open it hereProduct-formula (Trotter-Suzuki) simulation — click the name to read about itTrotter-Suzuki simulationTruncated-Taylor LCU simulation · 2 parts inside, open — click the name to close itTruncated-Taylor LCU simulationBlock-encode a matrix · 4 ways through — click the name to read about itBlock-encode a matrixTruncated-Taylor LCU simulation — click the name to read about itsum the truncated series, then amplifyQubitization walk simulation · 2 parts inside — click the name to read about itQubitization walk simulationHamiltonian you can query — you start hereCircuit for e^{-iHt} — you finish hereBlock-encoding — what one part hands to the next

Drag to pan. Pinch, or hold ctrl and scroll, to zoom. Arrow keys pan, plus and minus zoom, zero resets the view.

From Hamiltonian you can query to Circuit for e^{-iHt}

A circle is an object you are holding. This method is drawn heavier, opened into its own steps; the other lines between the same two ends are the alternatives recorded for the same slot. Circles are named on hover, and each one is a link.

2 lines here have ways through that this figure does not open. The map opens them in place. See it on the map

What it fills

  • Simulate Hamiltonian evolution

    Implement e−iHte^{-iHt} to error ε\varepsilon given some access model for HH. It is an application in its own right and also the engine inside phase estimation and several linear-system solvers.

When it applies

Requires HH 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∣∣||H|| — a decomposition with many comparable terms is expensive even for a benign HH.

Requires

These do not move the route along. The method needs each of them alongside its own work, and the cost of getting them is part of what the method costs.

  • Prepare an input state

    Map ∣0…0⟩|0\ldots0⟩ to a state whose amplitudes are proportional to a specified vector bb, to within ε\varepsilon. The cost is set by which description of bb 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⟩=(1/s)∑jβj∣j⟩B|0\rangle = (1/\sqrt{s})\sum_j \sqrt{\beta_j}|j\rangle with s:=∑jβjs := \sum_j \beta_j, over the register indexing j=(k,ℓ1,…,ℓk)j = (k,\ell_1,\ldots,\ell_k), whose amplitudes are the LCU coefficients β(k,ℓ1,…,ℓk)=(t/r)kk! αℓ1⋯αℓk\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 BB is a tensor product over the K+1K+1 registers: the order register goes from 0K0^K to a normalized superposition of the unary strings 1k0K−k1^k0^{K-k} in O(K)O(K) gates, by a cascade of singly-controlled rotations; each of the KK index registers goes from ∣0⟩|0\rangle to the normalized ∑ℓ=1Lαℓ∣ℓ⟩\sum_{\ell=1}^{L}\sqrt{\alpha_\ell}|\ell\rangle at a gate cost of O(L)O(L). The inverse is wanted as often as the preparation: the paper prices one amplification step at two instances of select(V)\mathrm{select}(V), one of select(V)†\mathrm{select}(V)^\dagger, and three each of BB and B†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 aa, 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†⊗I) select(V) (B⊗I)W = (B^\dagger \otimes I)\,\mathrm{select}(V)\,(B \otimes I), whose all-zeros ancilla branch reads PW∣0⟩∣ψ⟩=(1/s)∣0⟩U~∣ψ⟩PW|0\rangle|\psi\rangle = (1/s)|0\rangle\tilde{U}|\psi\rangle with P:=∣0⟩⟨0∣⊗IP := |0\rangle\langle 0| \otimes I and U~\tilde{U} the order-KK truncation of one segment's evolution Ur=e−iHt/rU_r = e^{-iHt/r} — and nothing at all about ∣ψ⟩|\psi\rangle. Hence the realization has to be the oblivious one: the ancilla reflection I−2PI - 2P that AA interleaves with WW and W†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 U~\tilde{U}, A∣0⟩∣ψ⟩=∣0⟩U~∣ψ⟩A|0\rangle|\psi\rangle = |0\rangle\tilde{U}|\psi\rangle, unmeasured, once per segment and rr times over the evolution — once, because s=2s = 2 leaves nothing for a second round. assumption: the robust form actually invoked returns PA∣0⟩∣ψ⟩=∣0⟩(3sU~−4s3U~U~†U~)∣ψ⟩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(δ)|s - 2| = O(\delta) and ∥U~−Ur∥=O(δ)\|\tilde{U} - U_r\| = O(\delta), before ∥PA∣0⟩∣ψ⟩−∣0⟩Ur∣ψ⟩∥=O(δ)\|PA|0\rangle|\psi\rangle - |0\rangle U_r|\psi\rangle\| = O(\delta) follows; a normalization near 22 without U~\tilde{U} near UrU_r buys nothing here.

    assumption

Example

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

Cost, as the source states it

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

Implementations

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

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

    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.

    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.

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

    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!`.

What it needs

  1. Block-encode a matrix 4 methods

    Wrap an operator AA inside a larger unitary UU so that A/αA/α sits in UU's top-left block, giving every routine above it one uniform way to touch the matrix. The subnormalization αα and the ancilla count are outputs of this layer, not free parameters.

Other ways to fill the same slot

Different approaches

  • Product-formula (Trotter-Suzuki) simulation

    Split HH into terms that can each be exponentiated directly and alternate short evolutions of them — the Lie-Trotter formula and its higher-order generalizations. No block-encoding is built and there is no all-zeros flag to amplify.

  • Qubitization walk simulation

    From a block-encoding pair (U,∣G⟩)(U, |G⟩) with H=(⟨G∣⊗I)U(∣G⟩⊗I)H = (⟨G|⊗I)U(|G⟩⊗I), build a walk operator WW that splits the Hilbert space into invariant two-dimensional SU(2)SU(2) subspaces, one per eigenvalue of HH, with eigenvalues e±i⋅arccos(H/α)e^{±i·arccos(H/α)}. Quantum signal processing phases applied to WW then produce e−iHte^{-iHt}.

In the Atlas

Sources