Approximating the ground state of many-body systems is a key computational bottleneck underlying important applications in physics and chemistry. The most widely known quantum algorithm for ground state approximation, quantum phase estimation, is out of reach of current quantum processors due to its high circuit-depths. Subspace-based quantum diagonalization methods offer a viable alternative for pre- and early-fault-tolerant quantum computers. Here, we introduce a quantum diagonalization algorithm which combines two key ideas on quantum subspaces: a classical diagonalization based on quantum samples, and subspaces constructed with quantum Krylov states. We prove that our algorithm converges in polynomial time under the working assumptions of Krylov quantum diagonalization and sparseness of the ground state. We then demonstrate the scalability of our approach by performing the largest ground-state quantum simulation of impurity models using a Heron quantum processors and the Frontier supercomputer. We consider both the single-impurity Anderson model with 41 bath sites, and a system with 4 impurities and 7 bath sites per impurity. Our results are in excellent agreement with Density Matrix Renormalization Group calculations.
Quantum algorithms based on classical processing of individual samples have recently emerged as a promising route to perform large-scale quantum-chemical calculations on pre-fault-tolerant and early-fault-tolerant quantum devices. In these algorithms, the quantum computer acts as a sampling engine that generates the subspace in which the Hamiltonian is classically diagonalized. The recently proposed Sample-based Krylov Quantum Diagonalization (SKQD) uses quantum Krylov states as circuits from which samples are collected. Convergence guarantees can be derived for SKQD under assumptions similar to those of quantum phase estimation -- namely, a sufficiently large overlap between the reference and the target ground state -- with the additional requirement that the ground-state wave function be well approximated by a polynomial subset of the full Hilbert space. However, implementations of SKQD for complex many-body Hamiltonians, such as quantum chemistry ones, are limited by the depths of time-evolution circuits needed to generate Krylov vectors. In this work, we introduce a method that combines SKQD with a qDRIFT randomized compilation of the Hamiltonian propagator. The resulting algorithm, termed SqDRIFT, enables executing quantum chemistry experiments on present-day quantum processors while retaining provable convergence guarantees, which we extend to account for the error of the randomized compilation and of finite sampling. We demonstrate its viability by applying SqDRIFT to calculate the electronic ground-state energy of several polycyclic aromatic hydrocarbons, up to system sizes beyond the reach of exact diagonalization.
Sampling-based quantum diagonalization method exploits Quantum-centric supercomputing platforms to sample bitstrings for Hamiltonian projection on a quantum computer, and then classically diagonalize the Hamiltonian to estimate the eigenvalues and eigenvectors. In current quantum devices an algorithm is useful when shallow quantum circuits with error mitigation support can discover better results while having either a proof of convergence or some method to explain trust in experiment. In this paper, we introduce DF-SQD, a hybrid algorithm that derives deterministic auxiliary-field circuits from selected double-factorization leaves of the two-electron tensor. The circuits propose occupation-number configurations, while selected configuration interaction evaluates the original active-space Hamiltonian and can recentre subsequent proposal rounds. On N2 (32 qubits; 6-31G basis) and a 40-qubit [Fe2S2(SCH3)4]2- active-space Hamiltonian, we show that DF-SQD improves the energy obtained from sampled determinant spaces while using shallow number-preserving circuits in both simulator and hardware runs. For N2, DF-SQD is 45x more accurate with a 11.23\% smaller subspace, and due to its ability to sample better bitstrings at lesser shots it is 2.93x faster than SQD in quantum devices. For the iron-sulfur cluster, DF-SQD generated a subspace dimension of 221M with 400K shots, while SQD needed 1.5M shots to generate a 238M subspace, thus we have better subspace recovery evident from the hardware at 3.75x reduced shots. At a matched 50M subspace dimension, DF-SQD is 1.32x more accurate (achieves a 24.5\% relative error reduction over standard SQD). So overall, our method is able to discover better results with shallower circuits, is sample efficient, uses configuration recovery (so has targeted error mitigation) and we have empirical convergence observation.
We apply the recently proposed Sample-based Krylov Quantum Diagonalization (SKQD) method to lattice gauge theories, using the Schwinger model with a
θ-term as a benchmark. SKQD approximates the ground state of a Hamiltonian, employing a hybrid quantum-classical approach: (i) constructing a Krylov space from bitstrings sampled from time-evolved quantum states, and (ii) classically diagonalizing the Hamiltonian within this subspace. We study the dependence of the ground-state energy and particle number on the value of the
θ-term, accurately capturing the model's phase structure. The algorithm is implemented on trapped-ion and superconducting quantum processors, demonstrating consistent performance across platforms. We show that SKQD substantially reduces the effective Hilbert space, and although the Krylov space dimension still scales exponentially, the slower growth underscores its promise for simulating lattice gauge theories in larger volumes.
Solving systems of linear equations,
Ax=b, is a ubiquitous task in computational science. The Harrow-Hassidim-Lloyd algorithm requires phase estimation, controlled rotations, and post-selection, impractical on noisy intermediate-scale quantum hardware. We extend sample-based quantum diagonalization (SQD), developed for quantum chemistry, to linear systems. Following the sample-based Krylov variant, the Hermitian matrix
A is decomposed in the Pauli basis, the Trotterized time evolution of
∣b⟩ is sampled at several times to select computational-basis configurations, and the system is solved classically by a Galerkin projection onto the selected subspace. We prove that the subspace solution is exact when the sampled configurations contain the Krylov support of
b and optimal in the energy norm for symmetric positive-definite matrices, so that noise-induced configurations cannot degrade it as long as they do not displace the ideal ones. On IBM Quantum hardware, the sampled distribution of a
4×4 example deviates from the exact one by eleven times the shot-noise level, yet the solution is recovered exactly. For a localized problem of dimension
N=1024, SQD selects 33 configurations, reaches a relative error of
10−5, and its truncated subspaces coincide with the optimal choice; readout errors enlarge the subspace without degrading the accuracy. The tests also show the limitations: the selected subspace coincides with that of inexpensive classical selection, banded matrices require exponentially many Pauli terms in the binary encoding, and Trotter errors populate spurious configurations unless the circuits are deep. An MIT-licensed Python package and notebooks that reproduce the numerical results are provided.
To solve linear systems
Ax=b, we propose a sample-based quantum linear solver combining quantum Krylov-state sampling with classical variational diagonalization. Rather than simulating the cost Hamiltonian
Hb(A)=A†(I−∣b⟩⟨b∣)A, we measure quantum Krylov states
e−iAT∣b⟩ generated by the original Hermitian operator
A in the computational basis. We then classically diagonalize the projected cost Hamiltonian to reconstruct an approximate solution. Using a representation of
A−1∣b⟩ as a linear combination of quantum Krylov states, we derive lower bounds on computational-basis sampling probabilities from solution amplitudes. We establish support-discovery guarantees for finite nonnegative-time grid sampling and signed-time sampling, which permits both positive and negative evolution times. The former bounds the maximum evolution time, whereas the latter improves ideal shot complexity at the cost of occasional long evolutions. For a condition-number bound
κ on
A, a lower bound
βL on the solution probability weight of each of
L target basis states, and failure probability
δ, signed-time sampling achieves
O(κ2βL−1log(L/δ)) shots, matching a worst-case lower bound up to constants within an evolution-only observation model in a specified parameter regime. Finally, we analyze the effects of finite-length qDRIFT compilation and circuit reuse on sampling and derive an infidelity bound for the reconstructed solution.