Sample-based Krylov Quantum Diagonalization (SKQD)

• SKQD is a hybrid method that constructs quantum-prepared Krylov subspaces to estimate low-energy eigenstates and eigenvalues of many-body Hamiltonians.
• The algorithm employs shifting and coefficient splitting techniques to drastically reduce sampling overhead, achieving cost reductions up to 500× in benchmark cases.
• Demonstrated across diverse systems like electronic structures, spin models, and lattice gauge theories, SKQD offers rigorous convergence guarantees and robustness against sampling errors.

Sample-based Krylov Quantum Diagonalization (SKQD) is a hybrid quantum–classical algorithm for estimating low-lying eigenvalues and eigenstates of large quantum many-body Hamiltonians. By leveraging quantum-prepared Krylov subspaces and sampling-based estimation of matrix elements, SKQD combines the rapid convergence of classical Krylov methods with hardware feasibility on pre-fault-tolerant and early fault-tolerant quantum computers. Sampling error and ill-conditioning in generalized eigenvalue problems (GEVP) are key challenges, addressed via statistical analysis and targeted error-reduction strategies. Recent implementations provide rigorous convergence guarantees and demonstrate utility across diverse systems, including electronic structure, spin models, and lattice gauge theories.

1. Formal Algorithmic Structure

SKQD algorithmically projects the Hamiltonian $\hat H$ onto a Krylov subspace generated by repeated application of a short-time propagator or Hamiltonian powers to an initial state $|\phi_0\rangle$. For discrete-time evolution, one constructs the Krylov subspace

$$K_m = \mathrm{span}\big\{\,|\phi_0\rangle,\,|\phi_1\rangle,\,\dots,\,|\phi_{m-1}\rangle\,\big\}, \quad |\phi_k\rangle = \hat B^k|\phi_0\rangle,\,\,\hat B = e^{-i \Delta t\,\hat H}$$

where $m$ is the subspace dimension. Matrix elements of the projected (Toeplitz) overlap and Hamiltonian matrices,

$$S_{ij} = \langle\phi_i|\phi_j\rangle, \qquad H_{ij} = \langle\phi_i|\hat H|\phi_j\rangle$$

are empirically estimated—typically via repeated projective quantum measurements (“shots”) on the quantum device.

The low-energy spectrum is then extracted by solving the generalized eigenvalue problem (GEVP),

$H \mathbf{w} = \lambda S \mathbf{w}$

where the smallest eigenvalue $\lambda$ approximates the ground-state energy $E_0$, converging exponentially in $m$ when $|\phi_0\rangle$ has nonzero overlap with the true ground state (Lee et al., 2024, Firt et al., 19 Dec 2025).

2. Sampling Error and Its Impact

Finite quantum measurement leads to statistical error in each matrix element $|\phi_0\rangle$0. Under standard Hamiltonian decompositions (linear combination of unitaries, diagonalizable fragments), the variance of $|\phi_0\rangle$1 scales as

$|\phi_0\rangle$2

where $|\phi_0\rangle$3 is the total shots, $|\phi_0\rangle$4 are decomposition coefficients, $|\phi_0\rangle$5 is the shot allocation, and $|\phi_0\rangle$6 is the variance of measured observables.

Sampling cost to achieve error $|\phi_0\rangle$7 scales as $|\phi_0\rangle$8, where $|\phi_0\rangle$9 is the 1-norm of Hamiltonian coefficients in the chosen decomposition. Sampling errors become particularly problematic as the GEVP becomes ill-conditioned for large $$K_m = \mathrm{span}\big\{\,|\phi_0\rangle,\,|\phi_1\rangle,\,\dots,\,|\phi_{m-1}\rangle\,\big\}, \quad |\phi_k\rangle = \hat B^k|\phi_0\rangle,\,\,\hat B = e^{-i \Delta t\,\hat H}$$0, so small statistical perturbations can lead to large eigenvalue errors (Lee et al., 2024, Lee et al., 2023).

3. Error-Reduction Techniques: Shifting and Coefficient Splitting

To mitigate sampling costs, two principal strategies have been devised (Lee et al., 2024):

(A) Shifting Technique

Redundant Hamiltonian fragments $$K_m = \mathrm{span}\big\{\,|\phi_0\rangle,\,|\phi_1\rangle,\,\dots,\,|\phi_{m-1}\rangle\,\big\}, \quad |\phi_k\rangle = \hat B^k|\phi_0\rangle,\,\,\hat B = e^{-i \Delta t\,\hat H}$$1 that annihilate all relevant Krylov vectors can be eliminated. A Hermitian shift operator $$K_m = \mathrm{span}\big\{\,|\phi_0\rangle,\,|\phi_1\rangle,\,\dots,\,|\phi_{m-1}\rangle\,\big\}, \quad |\phi_k\rangle = \hat B^k|\phi_0\rangle,\,\,\hat B = e^{-i \Delta t\,\hat H}$$2 with $$K_m = \mathrm{span}\big\{\,|\phi_0\rangle,\,|\phi_1\rangle,\,\dots,\,|\phi_{m-1}\rangle\,\big\}, \quad |\phi_k\rangle = \hat B^k|\phi_0\rangle,\,\,\hat B = e^{-i \Delta t\,\hat H}$$3 is introduced, yielding a shifted Hamiltonian

$$K_m = \mathrm{span}\big\{\,|\phi_0\rangle,\,|\phi_1\rangle,\,\dots,\,|\phi_{m-1}\rangle\,\big\}, \quad |\phi_k\rangle = \hat B^k|\phi_0\rangle,\,\,\hat B = e^{-i \Delta t\,\hat H}$$4

which preserves the spectrum up to a constant. The optimal shift minimizes the effective 1-norm,

$$K_m = \mathrm{span}\big\{\,|\phi_0\rangle,\,|\phi_1\rangle,\,\dots,\,|\phi_{m-1}\rangle\,\big\}, \quad |\phi_k\rangle = \hat B^k|\phi_0\rangle,\,\,\hat B = e^{-i \Delta t\,\hat H}$$5

Numerical studies show this reduces $$K_m = \mathrm{span}\big\{\,|\phi_0\rangle,\,|\phi_1\rangle,\,\dots,\,|\phi_{m-1}\rangle\,\big\}, \quad |\phi_k\rangle = \hat B^k|\phi_0\rangle,\,\,\hat B = e^{-i \Delta t\,\hat H}$$6 by $$K_m = \mathrm{span}\big\{\,|\phi_0\rangle,\,|\phi_1\rangle,\,\dots,\,|\phi_{m-1}\rangle\,\big\}, \quad |\phi_k\rangle = \hat B^k|\phi_0\rangle,\,\,\hat B = e^{-i \Delta t\,\hat H}$$7–$$K_m = \mathrm{span}\big\{\,|\phi_0\rangle,\,|\phi_1\rangle,\,\dots,\,|\phi_{m-1}\rangle\,\big\}, \quad |\phi_k\rangle = \hat B^k|\phi_0\rangle,\,\,\hat B = e^{-i \Delta t\,\hat H}$$8, effecting up to $$K_m = \mathrm{span}\big\{\,|\phi_0\rangle,\,|\phi_1\rangle,\,\dots,\,|\phi_{m-1}\rangle\,\big\}, \quad |\phi_k\rangle = \hat B^k|\phi_0\rangle,\,\,\hat B = e^{-i \Delta t\,\hat H}$$9 cost reduction in small-molecule benchmarks.

(B) Coefficient Splitting (ICS)

Pauli terms appearing in multiple measurement circuits are optimally split. Each coefficient $m$0 is decomposed as $m$1 across fragments. Overall variance is then a quadratic function of the splits, and alternately optimizing the shot allocations and coefficient splits yields optimal resource usage. ICS alone provides $m$2–$m$3 shot reduction for fragments, and $m$4–$m$5 for LCU, while combining with shifting achieves sampling overhead reductions of $m$6–$m$7 in representative molecules.

Molecule	LCU (SI)	LCU (ICS)	LCU (Shift)	LCU (Shift+ICS)
H$m$8	12.3	8.2	0.13	0.13
H$m$9	80.9	72.9	6.9	7.0
LiH	210	224	2.7	3.1

Shot cost $$S_{ij} = \langle\phi_i|\phi_j\rangle, \qquad H_{ij} = \langle\phi_i|\hat H|\phi_j\rangle$$0 (in Hartree$$S_{ij} = \langle\phi_i|\phi_j\rangle, \qquad H_{ij} = \langle\phi_i|\hat H|\phi_j\rangle$$1), from [(Lee et al., 2024), Table II]. The Shift+ICS strategy yields up to 500$$S_{ij} = \langle\phi_i|\phi_j\rangle, \qquad H_{ij} = \langle\phi_i|\hat H|\phi_j\rangle$$2 sampling cost reduction for H$$S_{ij} = \langle\phi_i|\phi_j\rangle, \qquad H_{ij} = \langle\phi_i|\hat H|\phi_j\rangle$$3O.

4. Theoretical Foundations and Convergence Guarantees

SKQD offers rigorous convergence under two principal hypotheses (Yu et al., 16 Jan 2025, Piccinelli et al., 4 Aug 2025):

• The reference state $$S_{ij} = \langle\phi_i|\phi_j\rangle, \qquad H_{ij} = \langle\phi_i|\hat H|\phi_j\rangle$$4 has non-negligible overlap with the true ground state ($$S_{ij} = \langle\phi_i|\phi_j\rangle, \qquad H_{ij} = \langle\phi_i|\hat H|\phi_j\rangle$$5).
• The ground state $$S_{ij} = \langle\phi_i|\phi_j\rangle, \qquad H_{ij} = \langle\phi_i|\hat H|\phi_j\rangle$$6 is $$S_{ij} = \langle\phi_i|\phi_j\rangle, \qquad H_{ij} = \langle\phi_i|\hat H|\phi_j\rangle$$7-sparse in the measurement basis, i.e. a polynomial number $$S_{ij} = \langle\phi_i|\phi_j\rangle, \qquad H_{ij} = \langle\phi_i|\hat H|\phi_j\rangle$$8 of basis states provide most of its norm.

The sampling complexity to guarantee recovery of all important bitstrings in the Krylov subspace with high probability is

$$S_{ij} = \langle\phi_i|\phi_j\rangle, \qquad H_{ij} = \langle\phi_i|\hat H|\phi_j\rangle$$9

for $H \mathbf{w} = \lambda S \mathbf{w}$0 Krylov steps and $H \mathbf{w} = \lambda S \mathbf{w}$1 measurements per step. The ground-state energy estimate then satisfies

$H \mathbf{w} = \lambda S \mathbf{w}$2

Convergence is polynomial in $H \mathbf{w} = \lambda S \mathbf{w}$3 when these parameters scale polynomially (Yu et al., 16 Jan 2025, Piccinelli et al., 4 Aug 2025).

5. Workflow and Implementation Methodologies

The practical SKQD workflow proceeds as follows (Firt et al., 19 Dec 2025, Yu et al., 16 Jan 2025, Piccinelli et al., 4 Aug 2025):

1. Initial state preparation: Prepare the reference state $H \mathbf{w} = \lambda S \mathbf{w}$4, potentially leveraging problem structure (e.g., singlet product states, W-states, or sector-specific ansätze).
2. Krylov basis generation: Apply $H \mathbf{w} = \lambda S \mathbf{w}$5 steps of Trotterized time evolution to produce $H \mathbf{w} = \lambda S \mathbf{w}$6. Circuit depths scale $H \mathbf{w} = \lambda S \mathbf{w}$7 per step for local models.
3. Sampling: For each Krylov state, perform projective measurement in the computational basis for $H \mathbf{w} = \lambda S \mathbf{w}$8 shots. Collect the union of all observed bitstrings.
4. Subspace construction: Form the subspace spanned by unique measured bitstrings, and classically compute the subspace-projected Hamiltonian matrix.
5. Diagonalization: Solve the associated eigenproblem to recover the lowest eigenvalues and eigenstates.

In molecular systems, block-encoding and Quantum Signal Processing (QSP) can be used to prepare Krylov eigenstates in depth linear in the subspace dimension, with expectation values for observables (e.g., reduced density matrices) measured in constant time per term, independent of $H \mathbf{w} = \lambda S \mathbf{w}$9 (Oumarou et al., 9 Jan 2025).

6. Sampling Error Analysis, Conditioning, and Thresholding

Being sample-based, SKQD is susceptible to noise amplification in ill-conditioned GEVPs. The spectral properties of the overlap matrix $\lambda$0 determine the sensitivity of eigenvalues to statistical perturbations. Nonasymptotic matrix concentration results yield operator-norm error bounds: $\lambda$1 where $\lambda$2 and $\lambda$3 for Toeplitz cases. The resultant eigenvalue error, after optimal basis truncation and thresholding, is sharply bounded as: $\lambda$4 for reduced subspace dimension $\lambda$5 and GEVP condition number $\lambda$6 (Lee et al., 2023). Overly aggressive subspace truncation increases bias, while insufficient thresholding leads to instabilities.

7. Applications and Demonstrated Performance

SKQD has been deployed across domains:

• Heisenberg models and materials science: SKQD recovers ground-state energies and magnetization curves for 1D and 2D Heisenberg systems, maintaining errors $\lambda$7 (easy-axis) and showing scaling consistent with exact diagonalization and DMRG. Implementations on IBM Heron hardware (18–30 qubits) validate algorithmic robustness in the presence of hardware noise (Firt et al., 19 Dec 2025).
• Electronic structure: Chemical systems, including polycyclic aromatic hydrocarbons up to 48 qubits, have been simulated with the SKQD variant incorporating qDRIFT randomized Trotterization (SqDRIFT). State-of-the-art results show convergence within chemical accuracy at modest circuit depths (Piccinelli et al., 4 Aug 2025).
• Lattice gauge theory: The Schwinger model with a $\lambda$8-term has been simulated using both trapped-ion and superconducting architectures, showing that SKQD achieves accurate ground-state and order-parameter estimation while reducing effective Hilbert space by up to $\lambda$9 for $E_0$0 (Rosanowski et al., 30 Oct 2025).

In molecular property computation, the QSP-enhanced SKQD framework offers constant-scaling measurements for 1- and 2-particle reduced density matrices and analytic nuclear gradients, a substantial advantage over traditional $E_0$1 scaling (Oumarou et al., 9 Jan 2025).

8. Outlook and Future Directions

SKQD and its variants (notably, SqDRIFT) establish rigorous, noise-robust quantum subspace diagonalization with sampling complexity that is polynomial under mild ground-state concentration hypotheses. The combination of shifting and coefficient splitting dramatically reduces sampling overhead, with resource scaling suitable for near-term and early-fault-tolerant quantum devices (Lee et al., 2024, Piccinelli et al., 4 Aug 2025). Current limitations involve the required circuit depths for highly entangled systems and the challenge of maintaining ground-state sparsity. Ongoing research is directed toward resource-efficient state preparation, improved basis rotations, and further optimizations for hardware connectivity and noise (Yu et al., 16 Jan 2025, Piccinelli et al., 4 Aug 2025).