Extending Sample-based Quantum Diagonalization to Systems of Linear Equations

Abstract

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\rangle$ 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\times4$ 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.

Can sample-based quantum diagonalization solve linear systems using quantum sampling followed by a classical subspace solve? This paper applies the sample-based Krylov protocol to Hermitian $Ax=b$: sample computational-basis configurations from $e^{-iAt_j}|b\rangle$, include $\operatorname{supp}(b)$, and solve the Galerkin-projected system. The implementation uses Pauli decomposition and first-order Lie–Trotter evolution, sharing its sampling-and-reconstruction structure with recently proposed quantum linear solvers. A Cayley–Hamilton and analytic-function argument establishes that exact evolution has the full Krylov support except at countably many times; capturing that support gives an exact solution when the projected matrix is invertible. For positive-definite $A$, standard Galerkin optimality guarantees nonincreasing energy-norm error under subspace enlargement, including noise-added configurations if ideal configurations remain present. For a screened-Poisson problem with $N=1024$, sampling exact evolution distributions selects 33 configurations and achieves relative error $1.1\times10^{-5}$. Truncated selections match the largest-solution-component oracle, but inexpensive classical Krylov and residual-greedy selection perform equally well. Actual circuits at $N=64,128$ achieve $1.7\times10^{-5}$ error using approximately $2.1\times10^4$ and $4.5\times10^4$ CNOTs per circuit. A two-qubit hardware demonstration recovers the solution despite distributional deviations eleven times the shot-noise level, because all four configurations are sampled. The stated limitations include $L=N$ Pauli terms for the tested tridiagonal family, Trotter-induced subspace inflation, poor compression for delocalized solutions, and sampling weights that need not track solution importance. Gate-noise robustness and hardware performance remain untested in the compressed regime; no advantage over classical selection is demonstrated.