Problem-aware sampling for bias and variance reduction in randomised Hamiltonian simulation

The paper in detail

With \(H=\sum_jc_jP_j\), each qDRIFT step samples \(j\) with probability \(q_j\) and applies \(e^{-i\theta_jP_j}\), \(\theta_j=c_jdt/q_j\). With \(M\) trajectories the error is \(B_N^2+V_N/M\), and only the bias does not shrink with \(M\). One step is exact: \(\hat O_1(j)=\langle O\rangle+\sin^2\theta_j\,\gamma_j+\tfrac12\sin2\theta_j\,\beta_j\), with \(\gamma_j=\langle P_jOP_j\rangle-\langle O\rangle\) and \(\beta_j=-i\langle[O,P_j]\rangle\). In the toy \(H=3Z_0Z_1+X_0+X_1\), vanilla qDRIFT spends 60% of samples on \(Z_0Z_1\), which commutes with the magnetization and has \(\gamma=0\). The state-aware rule moves weight to the transverse terms and cuts the root-mean-square error by 2–7 times.

For \(N\) steps, \(B_N=N(b_{\gamma,2}-\tfrac12\varphi_{H,H})dt^2+O(dt^3)\), with \(b_{\gamma,2}=\sum_jc_j^2\gamma_j/q_j\) and \(\varphi_{H,H}=d^2\langle O\rangle/dt^2|_0\). The variance starts at \(Nv_{\beta^2,2}dt^2\). Cross-step terms carry higher powers of \(N\) but still fall as \(1/N\) at fixed \(T\). If the state, observable and terms are real in some basis, as for real-orbital chemistry with Hartree–Fock references, then \(\beta_j=0\) and odd orders vanish. When the \(\gamma_j\) have mixed signs, a distribution can make \(b_{\gamma,2}=\tfrac12\varphi_{H,H}\) and cancel the \(dt^2\) bias.

Densities of the form \(\sum_jc_j^2f_j/q_j\) are minimized by \(q_j\propto|c_j|\sqrt{f_j}\), a Neyman allocation. This gives \(q_j\propto|c_j|\sqrt{|\gamma_j|}\) when \(\beta_j=0\), \(q_j\propto|c_j||\beta_j|\) otherwise, and the family \(q_j\propto|c_j|\sqrt{|\gamma_j|+\alpha\beta_j^2}\) between them. Along this family the bias density grows and the variance density shrinks, so for sign-definite \(\gamma\) some \(\alpha\) beats vanilla in both. Deeper circuits use per-term weights on five densities, optimized for mean-squared error with neither bias nor variance more than 5% worse than vanilla. Earlier state-aware rules are special cases with \(O=|\psi\rangle\langle\psi|\), for which \(\gamma_j\le0\) and the bias cannot be cancelled.

For random product states under five lattice models of 6–18 qubits, most single-step instances improve both bias and variance. For quenches between molecular geometries, from H\(_3\) to CH\(_4\), the steps needed to resolve the change in a weighted occupation fall 1.7–26 times. A validity ratio \(R^{(4)}_{\rm val}>2\) separated wins from losses in every case tested. For deployment, a Hartree–Fock surrogate misjudges correlated LiH and does worse than vanilla, while an MPS with \(\chi=4\) suffices. Matchgate shadows need only the 3-RDM for every distribution-dependent input and the 5-RDM for one fixed term, with \(O(n^5\log n/\varepsilon^2)\) snapshots.

Transferable methods

1. Neyman allocation for importance-sampled channels

   The leading bias and variance are sums \(\sum_jc_j^2f_j/q_j\). Cauchy–Schwarz gives the minimizer in closed form, and the same argument applies to any weighted sum of such densities, which yields a one-parameter family trading bias against variance.

   For \(f_j\ge0\) and \(\sum_jq_j=1\), \[\Bigl(\sum_j|c_j|\sqrt{f_j}\Bigr)^2=\Bigl(\sum_j\tfrac{|c_j|\sqrt{f_j}}{\sqrt{q_j}}\sqrt{q_j}\Bigr)^2\le\sum_j\tfrac{c_j^2f_j}{q_j}\sum_jq_j,\] with equality iff \(q_j\propto|c_j|\sqrt{f_j}\). With \(f_j=|\gamma_j|\) of one sign, the single-step variance \(v_\gamma=\operatorname{Var}_q[(c_j/q_j)^2\gamma_j]\) vanishes at the same \(q\), because the random variable becomes constant.

   The bias \((b_{\gamma,2}-\tfrac12\varphi_{H,H})dt^2\) is minimized at that \(q\) but cancels only if \(\tfrac12\varphi_{H,H}\) lies beyond \(b_\gamma^\star=\pm(\sum_j|c_j|\sqrt{|\gamma_j|})^2\). Mixed signs allow exact cancellation, which is how the staggered magnetization of the XY chain reaches \(|B_1|\approx10^{-4}\) of vanilla.

   Editorial check: for the toy model with \(|\psi_a\rangle\), \(O=Z_0+Z_1\) and \(T=0.2\), \(\gamma=(0,-1.732,-1)\). The 20%-mixed rule \(q=(0.12,0.49,0.39)\) lowers the exact single-step bias from \(-0.498\) to \(-0.339\) and the variance from 0.251 to 0.0072. The root-mean-square error at \(M=1\) falls from 0.707 to 0.349, within the 2–7 times range of Fig. 2.

   How to reuse it

   Any randomized compiler or importance-sampled estimator (randomized measurements, LCU sampling, Pauli-path sums) whose cost is \(\sum_jw_j^2f_j/q_j\) can be tuned with \(q_j\propto|w_j|\sqrt{f_j}\), provided \(f_j\) is cheap to estimate in advance. Mix with the baseline so that no probability falls below \(\tau\) times its baseline value, which bounds every reweighting by \(1/\tau\).

2. Majorana degree counting fixes the RDM order

   To estimate the distribution inputs on a fermionic state, count how many Majorana operators the nested commutators can contain. This fixes which reduced density matrix is needed and therefore the matchgate-shadow cost.

   Majorana monomials \(\zeta_{\mathcal I}\) and \(\zeta_{\mathcal I'}\) with \(|\mathcal I'|\) even commute unless \(|\mathcal I\cap\mathcal I'|\) is odd, in which case \([\zeta_{\mathcal I},\zeta_{\mathcal I'}]=\pm2\zeta_{\mathcal S}\) with \(\mathcal S\) the symmetric difference and \(|\mathcal S|\le|\mathcal I|+|\mathcal I'|-2\). Hamiltonian terms have degree 2 or 4, so each commutator raises the degree by at most 2.

   A one-body \(O\) has degree 2. Conjugation by \(P_j\) flips the sign of pairs with exactly one index in \(\mathcal I_j\), so \(\gamma_j=-2i\sum_{\mu\in\mathcal I_j,\nu\notin\mathcal I_j}A_{\mu\nu}\langle\zeta_\mu\zeta_\nu\rangle\) needs only the 1-RDM. \(\beta_j\) needs the 2-RDM and \(\varphi_{j,b}\) the 3-RDM. The fixed term \(\langle[H,[H,[H,[H,O]]]]\rangle\) has degree up to 10 and needs the 5-RDM.

   Matchgate shadows estimate a degree-\(2k\) correlator with single-snapshot variance at most \(\binom{2n}{2k}/\binom nk\approx(2n)^k/(2k-1)!!\), polynomial in \(n\) at fixed \(k\), whereas local Pauli shadows scale as \(3^w\) in the Pauli weight.

   How to reuse it

   Before measuring nested-commutator quantities on fermionic states (Trotter error bounds, Krylov moments, perturbative corrections), count Majorana degrees to find the highest RDM order needed. One matchgate-shadow dataset then serves every such quantity, provided the Hamiltonian has at most quartic fermionic terms.

Source: Section II (Eqs. 1–22); Section III (Eqs. 25–39, Lemma 3, Corollary 2); Section VI (Eqs. 53–58). Based on the linked paper version. Derivations marked as editorial, and suggested uses in other problems, are editorial guidance rather than results of the paper.

Context & limitations

Why it matters. It is an ancilla-free change to randomized simulation that targets the quantity actually estimated. Its inputs can come from a bond-dimension-4 matrix-product state surrogate (for LiH) or from matchgate shadows with polynomially many snapshots, without full classical simulation.

Limitations. The expansion is fourth order. Its truncation bound scales as \((\lambda T)^3\) and is loose, so the method targets short times: water fails at \(T=1\), as the validity metric predicts. Most benchmarks use exact inputs on classically simulable sizes. A Hartree–Fock surrogate performs worse than vanilla qDRIFT for LiH. The distribution is fixed at \(t=0\), and gate savings depend on the cost model.

Description based on the full text of the linked arXiv version, read 2026-10-25; equations follow the paper’s notation.