LCHS with the improved kernel

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.

A narrower version of LCHS — linear combination of Hamiltonian simulation

When it applies

Carries the same requirement as the original LCHS: the Hermitian part $L(t) = (A(t)+A(t)^†)/2$ must be positive semi-definite throughout the interval. The authors describe this as the first approach enabling quantum algorithms to solve linear differential equations with both optimal state preparation cost and near-optimal scaling in matrix queries on all parameters, which is why it is the current reference point for this layer.

Requires

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

Example

The numerics in arXiv:2312.03916 measure the kernel's truncation, not a run of the solver. Section 2.2 with Figure 2 approximates the single matrix function e^{-(L+iH)} for random 8 × 8 Hermitian matrices and plots how the truncation error of the k-integral falls as the cutoff grows, against the original Cauchy kernel. No differential equation is solved.

given  A(t) = L(t) + iH(t) with L(t) positive semi-definite on [0,T],  u_0,  b(t),  target error eps
       beta in (0,1) fixed
       kernel  f(z) = 1 / ( C_beta * exp( (1+iz)^beta ) ),   C_beta = 2*pi*exp(-2^beta)
       g(k)   = f(k) / (1 - ik)
       U(T,k) = T exp( -i * int_0^T ( k L(s) + H(s) ) ds )

# 1. recast: truncate and discretize the k-integral of the LCHS identity   [Sec. 3.1]
K  = O( (log(1/eps))^(1/beta) )                  # truncate to [-K,K]              [Lemma 10]
h1 = 1 / ( e * T * max_t ||L(t)|| )              # composite quadrature step size  [Lemma 11]
Q  = O( log(1/eps) )                             # Gauss nodes per subinterval     [Lemma 11]
M  = 2*K*Q/h1 = O( T * max_t ||L(t)|| * (log(1/eps))^(1+1/beta) )   # unitaries    [Lemma 11]
c_{q,m} = w_q * g(k_{q,m})  at the Gauss nodes k_{q,m};  sum_{q,m} |c_{q,m}| = O(1)  [Lemma 12]

# 2. simulate: one Hamiltonian simulation problem per node                  [Sec. 4.2.1]
build HAM-T for k_j L + H from HAM-T_L, HAM-T_H by an inner LCU:
    O_k loads k_j, controlled rotation c-R, controlled HAM-T_L / HAM-T_H, then c-R^T
    a second rotation makes the block-encoding factor uniform at alpha_L*K + alpha_H
SEL = sum_j |j><j| (x) W_j,   W_j block-encodes V_j with ||V_j - U(T,k_j)|| <= eps_1
    implemented by the truncated Dyson series method
    cost: Otilde( alpha_A * K * T * log(1/eps_1) ) queries to HAM-T_L, HAM-T_H  [Thm. 15 proof]

# 3. combine by LCU and post-select                                          [Sec. 4.2.1]
apply (O_{c,l}^dagger (x) I) SEL (O_{c,r} (x) I)  to  |0>|u_0>
    = |0> (1/||c||_1) ( sum_j c_j V_j ) |u_0>  +  |perp>
the inhomogeneous term is prepared the same way over nodes s_{j'} and combined with it  [Sec. 4.2.2]
measure the ancillas and keep the all-zeros branch
O(||c||_1) rounds of amplitude amplification for a constant-level success probability   [Lemma 24]

Cost, as the source states it

$\tilde{O}( ((||u_0|| + ||b||_{L^1})/||u(T)||) \alpha_A T (\log(1/\varepsilon))^{1+1/\beta} )$ matrix queries, with $\alpha_A \geq \max_t ||A(t)||$, $T$ the evolution time and $\beta ∈ (0,1)$; this improves to $(\log(1/\varepsilon))^{1/\beta}$ for time-independent $A$. State preparation costs $O( (||u_0|| + ||b||_{L^1})/||u(T)|| )$ queries, independent of both $T$ and $\varepsilon$.

Implementations

• Truncated Dyson series select oracle combined by LCU

  • Quantum algorithm for linear non-unitary dynamics with near-optimal dependence on all parameters

  The paper's main construction, and the one Theorem 15 prices. The problem is the general inhomogeneous system $du(t)/dt = -A(t)u(t) + b(t)$ with $u(0) = u_0$, and the deliverable is an $\varepsilon$-approximation of the normalized solution state $|u(T)\rangle$, prepared with $\Omega(1)$ probability and a flag indicating success. The solution representation is rewritten by the improved LCHS formula as a weighted integral of unitary evolutions, both integrals are discretized by composite Gaussian quadrature, each surviving evolution $U(T,s,k_j)$ is implemented by the truncated Dyson series method, and the sums are taken by LCU. The choice of simulation subroutine is not forced and the authors say so: any time-dependent Hamiltonian simulation algorithm can implement $U(T,s,k)$, and they name the product formula as an alternative with worse asymptotic scaling that works under weaker oracle assumptions and may give a simpler circuit. The truncated Dyson series is used because it gives the lowest asymptotic complexity.

  The select oracle is assembled in two layers. Starting from $|j\rangle |l\rangle |0\rangle_k |0\rangle_R |0\rangle_a |\psi\rangle$, where $j$ indexes the quadrature node and $l$ the time step inside the Dyson interval, the oracle $O_k$ writes $k_j$ into the $k$ register; a controlled rotation $c\text{-}R$ splits the $R$ register with amplitudes $\sqrt{\alpha_L k}$ and $\sqrt{\alpha_H}$ over $\sqrt{\alpha_L|k| + \alpha_H}$; controlled HAM-T oracles for $L$ and $H$ act on the two branches; and $c\text{-}R^\top$ recombines them into a block-encoding of $(k_j L + H)/(\alpha_L |k_j| + \alpha_H)$. That factor still varies with $k_j$, so one further ancilla qubit and the controlled rotation $|k\rangle|0\rangle \mapsto |k\rangle \left( \frac{\alpha_L|k|+\alpha_H}{\alpha_L K + \alpha_H}|0\rangle + \sqrt{1 - \left(\frac{\alpha_L|k|+\alpha_H}{\alpha_L K + \alpha_H}\right)^2}|1\rangle \right)$ pin it at the worst case $\alpha_L K + \alpha_H$ for every node, after which $O_k^\dagger$ clears the $k$ register; the ancilla count grows by $O(\log K) + 2$. The truncated Dyson series method then gives $\mathrm{SEL} = \sum_j |j\rangle\langle j| \otimes W_j$, where $W_j$ block-encodes some $V_j \approx U(T,k_j)$ with $||V_j|| \le 1$, and $(O_{c,l}^\dagger \otimes I)\,\mathrm{SEL}\,(O_{c,r} \otimes I)$ block-encodes the homogeneous propagator. The inhomogeneous part reuses the same machinery: the varying lower limit is absorbed by writing $U(T,s,k) = \mathcal{T}e^{-i\int_0^T (k \mathbf{1}_{s' \ge s}L(s') + \mathbf{1}_{s' \ge s}H(s'))\,ds'}$, whose HAM-T oracles cost one application of the originals plus a comparison written onto an extra ancilla qubit, and a doubly-indexed select oracle covers both $j$ and $j'$ with $O_b$ loading $|b(s_{j'})\rangle$. The two parts are finally joined by one more LCU on a single extra qubit through $R: |0\rangle \mapsto (\sqrt{||u_0||}\,|0\rangle + \sqrt{||c'||_1}\,|1\rangle)/\sqrt{||u_0|| + ||c'||_1}$, and postselecting all the ancilla registers on 0 yields the desired state.

  What the construction consumes is oracles, not matrices. For the coefficient matrix it assumes HAM-T access — a simultaneous block-encoding of the matrix at the $M_D$ time points inside each short interval $[qh,(q+1)h]$ — namely $\mathrm{HAM\text{-}T}_{A,q}$ with $\alpha_A \ge \max_t ||A(t)||$, from which $\mathrm{HAM\text{-}T}_{L,q}$ and $\mathrm{HAM\text{-}T}_{H,q}$ with $\alpha_L, \alpha_H = O(\alpha_A)$ follow by one LCU with a controlled rotation and one extra ancilla; the paper adds that these two may instead be assumed directly, which can give better block-encoding factors in some applications. The quadrature coefficients arrive as a state-preparation pair $(O_{c,l}, O_{c,r})$ loading $\sqrt{c_j}$ and $\sqrt{\overline{c_j}}$ with gate complexity $O(\alpha_L T (\log(1/\varepsilon))^{1+1/\beta})$, and the inhomogeneous weights as a matching pair $(O_{c',l}, O_{c',r})$; the square root is taken by principal argument, which is what makes it well defined for the improved kernel's complex coefficients. The nodes themselves come from $O_k: |j\rangle|0\rangle \mapsto |j\rangle|k_j\rangle$, built by classical arithmetic. The initial condition enters as $O_u: |0\rangle \mapsto |u_0\rangle$ and the source term as $O_b: |j'\rangle|0\rangle \mapsto |j'\rangle|b(s_{j'})\rangle$.

  Theorem 15 is what came out, and it is a proof rather than a measurement. Taking $M = \tilde{O}(\alpha_A T (\log((||u_0||+||b||_{L^1})/(||u(T)||\varepsilon)))^{1+1/\beta})$ nodes for the homogeneous part and $M' = \tilde{O}(T(\Lambda+\Xi)(\log((1+||b||_{L^1})/(||u(T)||\varepsilon)))^{1+1/\beta})$ for the inhomogeneous one, where $\Lambda$ and $\Xi$ collect the time derivatives of $A$ and $b$, the algorithm prepares the $\varepsilon$-approximation with $\tilde{O}(((||u_0||+||b||_{L^1})/||u(T)||)\,\alpha_A T (\log(1/\varepsilon))^{1+1/\beta})$ queries to the HAM-T oracle and $O((||u_0||+||b||_{L^1})/||u(T)||)$ queries to $O_u$ and $O_b$. Inside one run the simulation costs $\tilde{O}(\alpha_A K T \log(1/\min\{\varepsilon_1,\varepsilon_2\}))$ HAM-T queries; the outer factor is the $O(||c||_1(||u_0||+||c'||_1)/||v||)$ rounds of amplitude amplification, which the proof reduces to $O((||u_0||+||b||_{L^1})/||u(T)||)$ using $||v|| \ge ||u(T)||(1-\varepsilon/2)$. Corollary 16 restates the homogeneous case at $||b||_{L^1} = 0$. Table 1 sets the result beside the spectral method, the linear-system truncated Dyson method, time marching and the original LCHS; the improved-LCHS row reads $\tilde{O}((||u_0||/||u(T)||)\alpha_A T(\log(1/\varepsilon))^{1+1/\beta})$ matrix queries at $O(||u_0||/||u(T)||)$ state preparations, where the original LCHS's matrix dependence on precision is linear in $1/\varepsilon$ and time marching is quadratic in $T$ and $||A(t)||$. No circuit was executed and none was simulated; the paper reports no run of this construction.

• Time-independent specialization simulated by QSP/QSVT

  • Quantum algorithm for linear non-unitary dynamics with near-optimal dependence on all parameters

  A specialization for $A(t) \equiv A$. Corollary 16 already covers this case, but the time-independent setting admits better Hamiltonian simulation algorithms than the truncated Dyson series, and the point of the specialization is one property of them: in QSP and QSVT the Hamiltonian norm and the precision enter additively rather than multiplicatively, so a cutoff $K$ that grows with $\log(1/\varepsilon)$ stops multiplying the precision cost. The authors credit Dominic Berry with pointing out the complexity improvement for time-independent inhomogeneous equations.

  The HAM-T family collapses to a standard block-encoding $\mathrm{BE}_A$ with $(\langle 0|_a \otimes I)\,\mathrm{BE}_A\,(|0\rangle_a \otimes I) = A/\alpha_A$. The select-oracle construction is unchanged up to the block-encoding of $(k_j L + H)/(\alpha_L K + \alpha_H)$; from there QSVT rather than the truncated Dyson series produces $\mathrm{SEL}$, and the LCU combination follows as before. For a time-dependent source term over a time-independent $A$, the differing evolution times $T - s_{j'}$ are handled by one further controlled rotation $|j'\rangle|0\rangle \mapsto |j'\rangle\left((1 - s_{j'}/T)|0\rangle + \sqrt{1-(1-s_{j'}/T)^2}\,|1\rangle\right)$, which turns the oracle into a HAM-T for the rescaled Hamiltonian $(1-s/T)(kL+H)$. That rescaled Hamiltonian only has to be simulated for the single time $T(\alpha_L K + \alpha_H)$, which is independent of both $j$ and $j'$, so one set of QSVT phase factors — those of the polynomial approximation of $e^{-iT(\alpha_L K + \alpha_H)x}$ — serves every node.

  Only the matrix input model changes: a standard block-encoding of the single matrix $A$ with $\alpha_A \ge ||A||$, in place of the HAM-T family at $M_D$ time points. The state-preparation oracles $O_u$ and, in the inhomogeneous variant, $O_b$ are unchanged.

  Corollary 17: $\tilde{O}((||u_0||/||u(T)||)\,\alpha_A T(\log(1/\varepsilon))^{1/\beta})$ queries to the block-encoding of $A$ and $O(||u_0||/||u(T)||)$ queries to $O_u$. The precision dependence falls from $(\log(1/\varepsilon))^{1+1/\beta}$ to $(\log(1/\varepsilon))^{1/\beta}$, which almost matches the $\tilde{\Omega}(\log(1/\varepsilon))$ lower bound as $\beta \to 1$; the mechanism is that each run costs $O((\alpha_L K + \alpha_H)T + \log(1/\varepsilon_1)) = O(\alpha_A K T + \log(1/\varepsilon_1))$ queries, additive in the precision rather than multiplicative. The inhomogeneous variant keeps the same $\tilde{O}((\log(1/\varepsilon))^{1/\beta})$ dependence, the Hamiltonian simulation step still dominating. This is the last row of Table 1. Proved, not measured — no run is reported.

• Gibbs state preparation from the time-independent algorithm

  • Quantum algorithm for linear non-unitary dynamics with near-optimal dependence on all parameters

  An application rather than a new solver: preparing the Gibbs state $\frac{1}{Z_\gamma}e^{-\gamma L}$ of a positive semi-definite Hamiltonian $L$, where $\gamma$ is the inverse temperature and $Z_\gamma = \mathrm{Tr}(e^{-\gamma L})$. The connection is that $e^{-\gamma L/2}$ is exactly the kind of non-unitary operator the time-independent homogeneous algorithm block-encodes, so the ODE machinery is used unchanged and only the reading of its output differs.

  Prepare the purified Gibbs state $|\psi\rangle = \sqrt{N/Z_\gamma}\,(I \otimes e^{-\gamma L/2})\,\frac{1}{\sqrt{N}}\sum_{j=0}^{N-1}|j\rangle|j\rangle$ by constructing a block-encoding of $e^{-\gamma L/2}$ with the improved LCHS method in its time-independent form and applying it to the maximally entangled state, then trace out the first register to leave the Gibbs state itself. The complexity is read straight off Corollary 17 on observing that the norm of the unnormalized solution $(I \otimes e^{-\gamma L/2})\frac{1}{\sqrt{N}}\sum_j |j\rangle|j\rangle$ is $\sqrt{Z_\gamma/N}$, which is what turns the ODE's $||u_0||/||u(T)||$ into $\sqrt{N/Z_\gamma}$.

  One block-encoding of $L$ with factor $\alpha_L \ge ||L||$, and nothing more: no block-encoding of the shifted Hamiltonian $I - L/\alpha_L$, none of $\sqrt{L}$, and no assumption that $L$ is strictly positive definite rather than positive semi-definite. The maximally entangled state is prepared rather than supplied.

  Corollary 18: an $\varepsilon$-approximation of the purified Gibbs state, with $\Omega(1)$ probability and a flag indicating success, using $\tilde{O}\left(\sqrt{N/Z_\gamma}\,\gamma\alpha_L(\log(1/\varepsilon))^{1/\beta}\right)$ queries to the block-encoding of $L$. The authors place this against prior work: several earlier algorithms require $L$ strictly positive definite, or require a block-encoding of $I - L/\alpha_L$ or of $\sqrt{L}$, and the recent algorithm for a general Hamiltonian scales super-linearly in $\alpha_L$ and $\gamma$ and exponentially in $\mathrm{poly}(\log(1/\varepsilon))$ when applied to a positive semi-definite $L$. The LCHS route needs only the ability to simulate the dynamics of $L$, and is almost linear in $\gamma$ and $\alpha_L$ and poly-logarithmic in $1/\varepsilon$; Appendix B carries the longer comparison. Proved, not measured — no state was prepared on hardware or in simulation.

• Adaptive choice of the kernel parameter β

  • Quantum algorithm for linear non-unitary dynamics with near-optimal dependence on all parameters

  Everywhere else in the paper $\beta$ is a fixed constant in $(0,1)$ and the constants hidden inside the big-$O$ may depend on it. This variant lets $\beta$ be a function of the evolution time, the target error and the other parameters of the problem. The two effects pull against each other: as $\beta$ approaches 1 the asymptotic query complexity improves, but the constant factor blows up. The authors credit Robin Kothari with suggesting the adaptive choice.

  The proof keeps the $\beta$ dependence that the rest of the paper hides, through the three places it enters — the truncation $K$, and the coefficient norms, where $||c||_1 = O(1/(1-\beta))$ and $||c'||_1 = O(||b||_{L^1})$ — and then minimizes the total query complexity over $\beta$. The choice is $\beta = 1 - O\!\left(1/\log\log\frac{||u_0||+||b||_{L^1}}{||u(T)||\varepsilon}\right)$, at which the two diverging factors resolve: $1/(1-\beta)^{1+1/\beta} = O\!\left(\left(\log\log\frac{||u_0||+||b||_{L^1}}{||u(T)||\varepsilon}\right)^2\right)$ and $\left(\frac{1}{1-\beta}\log\frac{||u_0||+||b||_{L^1}}{||u(T)||\varepsilon}\right)^{1+1/\beta} = O\!\left(\left(\log\frac{||u_0||+||b||_{L^1}}{||u(T)||\varepsilon}\right)^2\right)$. Nothing in the circuit changes; what changes is which constant the kernel is instantiated at.

  Corollary 19: with that adaptive $\beta$, the algorithm prepares the $\varepsilon$-approximation with $\Omega(1)$ probability and a success flag using $\tilde{O}\left(\frac{||u_0||+||b||_{L^1}}{||u(T)||}\alpha_A T \left(\log\frac{1}{\varepsilon}\right)^2\right)$ queries to the HAM-T oracle and $\tilde{O}\left(\frac{||u_0||+||b||_{L^1}}{||u(T)||}\log\log\frac{1}{\varepsilon}\right)$ queries to the state preparation oracles. The matrix-query exponent drops from $1+1/\beta$, which is always greater than 2, to exactly 2. The paper is plain that this is a trade rather than a free improvement: the state preparation cost picks up the extra $\log\log(1/\varepsilon)$ factor, so this variant does not have the $T$- and $\varepsilon$-independent state preparation cost that the fixed-$\beta$ algorithm attains. Proved, not measured.

• Hybrid importance-sampling implementation for early fault-tolerant devices

  • Quantum algorithm for linear non-unitary dynamics with near-optimal dependence on all parameters

  An alternative for the early fault-tolerant regime. The coherent implementation uses several extra ancilla registers and coherently controlled Hamiltonian simulation, and the paper says both are computational challenges where there are few logical qubits and limited capability to handle complicated control structures. This variant follows the same hybrid quantum-classical approach as the original LCHS paper, with one nuance the improved kernel forces: its coefficients $c_j$ are no longer real and positive, so the real and imaginary parts must be split, estimated as two separate sums, and their signs carried. It is stated for the homogeneous case with $||u_0|| = 1$, and it answers a different question from the quantum algorithm — it estimates the observable $u(T)^\dagger O u(T)$ of the possibly unnormalized solution for a known Hermitian $O$, rather than preparing the solution state.

  The discretized LCHS gives $u(T) \approx \sum_j c_j U(T,k_j)|u_0\rangle$, so the observable becomes the double sum $\sum_{l,j} \overline{c_l}\,c_j \langle u_0 | U(T,k_l)^\dagger O\, U(T,k_j) | u_0 \rangle$, split into a real-weighted and an imaginary-weighted part. Four steps. Classically sample $(j,l)$ with probability proportional to $|\mathrm{Re}(\overline{c_l}c_j)|/\Gamma$ and $(j',l')$ with probability proportional to $|\mathrm{Im}(\overline{c_{l'}}c_{j'})|/\Gamma'$, where $\Gamma = \sum_{j,l}|\mathrm{Re}(\overline{c_l}c_j)|$ and $\Gamma' = \sum_{j',l'}|\mathrm{Im}(\overline{c_{l'}}c_{j'})|$. For the sampled indices, estimate $o_{j,l} = \langle u_0|U(T,k_l)^\dagger O\, U(T,k_j)|u_0\rangle$ on the quantum computer, by the Hadamard test for non-unitary matrices together with amplitude estimation. Classically average the sampled $\mathrm{sgn}(\mathrm{Re}(\overline{c_l}c_j))\,o_{j,l}$ and $\mathrm{sgn}(\mathrm{Im}(\overline{c_l}c_j))\,o_{j',l'}$ into $\sigma$ and $\sigma'$. Output $\Gamma\sigma + i\Gamma'\sigma'$. Each circuit uses one additional ancilla qubit and needs no $k$-controlled implementation of $U(T,k)$. The authors note that the importance-sampling strategy is general to LCU and not specific to LCHS.

  Access to a state preparation oracle $O_u$ for $|u_0\rangle$, a block-encoding $U_O$ of the observable with $\alpha_O \ge ||O||$, and a quantum circuit $\tilde{U}(T,k)$ approximating $U(T,k)$ to error at most $\varepsilon_{HS}$ for any $k$. The coefficients $c_j$ stay on the classical side as numbers to be sampled from, rather than being loaded into an ancilla state by a prepare pair.

  Proposition 20, which rephrases Theorem 9 of the original LCHS paper: $u(t)^{*}Ou(t)$ is estimated to error at most $\varepsilon$ with probability at least $1-\delta$ by choosing $\varepsilon_{HS} = O(\varepsilon/||O||)$, using $O\!\left(\left(\frac{||O||}{\varepsilon}\right)^2\log\frac{1}{\delta}\right)$ samples, each circuit making $O\!\left(\frac{\alpha_O}{\varepsilon}\log\frac{\alpha_O}{\varepsilon}\log\frac{||O||}{\delta\varepsilon}\right)$ queries to $O_u$ and to controlled versions of $U_O$ and $\tilde{U}(T,k)$. The paper states the trade in both directions: fewer ancilla qubits and no coherently controlled simulation, but an overall query complexity that becomes worse. Proved, not measured — no execution of the hybrid on hardware or in simulation is reported.

What it needs

1. Recast a non-Hermitian generator as Hamiltonian evolution 2 methods

   Given a generator $A(t)$ whose evolution is not unitary, produce a Hermitian generator — or a quadrature-indexed family of them — on a space at least as large, whose unitary evolution reproduces the original dynamics, together with the map that recovers the original solution. Both halves are required: a construction that reaches a Hamiltonian and cannot get back is not a route.

2. Simulate Hamiltonian evolution 3 methods

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

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.

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

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.

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

• 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

• Linear combination of unitaries

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

Sources

• 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
• Linear combination of Hamiltonian simulation for nonunitary dynamics with optimal state preparation cost↗Dong An, Jin-Peng Liu, Lin Lin · 2023theory yes · simulation no · hardware no