⚛️ quantum physics

A scalable sample-based quantum diagonalization solver for general Pauli-sum Hamiltonians

This paper presents a scalable, JAX-based ground-state solver for sample-based quantum diagonalization that leverages GPU and multi-GPU architectures to perform subspace diagonalization in spaces up to $2^{30}$ dimensions with demonstrated memory efficiency and strong scaling.

Original authors: Yutaro Iiyama

Published 2026-10-12

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Yutaro Iiyama

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

In the quest to understand the deepest workings of matter, scientists often turn to quantum mechanics, the set of rules that governs the behavior of atoms and subatomic particles. For systems containing just a few particles, these rules are manageable, but as the number of particles grows, the complexity explodes. Calculating the lowest energy state of a large collection of interacting particles—a task essential for designing new materials or understanding chemical reactions—becomes so difficult that even the most powerful classical supercomputers struggle to solve it. This is where the promise of quantum computing enters the story. The idea is to use a quantum machine, which naturally follows these complex rules, to help solve problems that are otherwise intractable. However, current quantum computers are still in their early stages, prone to errors and limited in size. They cannot yet run the massive, perfect calculations required for real-world applications on their own.

To bridge this gap, researchers have developed a hybrid approach known as sample-based quantum diagonalization. In this method, a small quantum computer does not solve the whole problem. Instead, it acts as a specialized sensor, sampling specific states of the system to suggest which parts of the vast mathematical landscape are most important. The heavy lifting of the actual calculation is then handed over to a classical computer, which uses these suggestions to find the answer. This division of labor allows scientists to leverage the unique strengths of both technologies, potentially achieving useful results even before quantum computers become fully error-free. The challenge, however, lies in the classical side of the equation: as the quantum computer suggests more and more states to consider, the classical computer must be able to handle calculations involving billions of dimensions without running out of memory or time.

A researcher led by Yutaro Iiyama at the University of Tokyo has addressed this bottleneck by creating a new, highly scalable software tool designed to run on modern high-performance computing facilities. Their work focuses on the classical engine that processes the data from the quantum sampler. They have built a solver capable of performing the necessary mathematical operations on subspaces containing up to roughly one billion dimensions, a scale that was previously difficult to manage on a single machine. The software is written to run efficiently on graphics processing units, the powerful chips typically found in video cards but now central to supercomputing, allowing it to harness the parallel power of hundreds of these processors working together.

The core of this new solver is a specialized algorithm that finds the lowest energy state of a system without ever needing to store the entire mathematical map of that system in memory. In traditional approaches, calculating the energy of a complex system often requires creating a massive grid of numbers representing every possible interaction, which quickly fills up a computer's memory. This new method avoids that trap by calculating the necessary numbers only when they are needed, a technique known as a matrix-free approach. It works by repeatedly applying a specific operation to a vector of numbers, refining the estimate of the ground state with each step. Because the researcher optimized this process to use very little memory, they demonstrated that a single graphics card with 100 gigabytes of memory could solve problems involving a billion dimensions, a feat that would have been impossible with older methods.

To test the limits of their system, the researcher ran simulations on a model of a magnetic material known as the transverse-field Ising model. They successfully calculated the energy levels of a system with 30 particles on a single high-end graphics card, confirming that the software could handle non-trivial physics problems with the expected efficiency. They then pushed the system further by distributing the calculation across a cluster of up to 128 graphics cards. In this setup, the problem was split among the machines, and the results showed that the time required to find the answer decreased almost perfectly as more computers were added. This "strong scaling" indicates that the software can effectively utilize the massive parallel power of modern supercomputing centers, making it possible to tackle problems that are far too large for any single machine.

The software also includes a specialized function to handle the unique way the quantum computer suggests which states to include. This function must quickly identify how different quantum states relate to one another under the influence of the system's Hamiltonian, which is a mathematical description of the total energy. The researcher developed a clever search method that groups these states to minimize the number of expensive operations required. They found that by caching certain static data, they could speed up the process, though this required a careful balance between using more memory and saving time. The system is currently limited by the ability of the underlying software library to sort lists of numbers, capping the number of dimensions it can handle at about one billion. However, the researcher has already outlined a plan to upgrade the sorting method to allow for even larger scales in the future.

This work represents a significant step forward in making quantum-centric supercomputing a practical reality. By providing a robust, open-source tool that can run on both single workstations and massive supercomputing clusters, the researcher has lowered the barrier for others to apply these advanced techniques. The solver is designed to be flexible, capable of handling general types of quantum systems described by sums of Pauli operators, which cover a wide range of physical models. While the current version has a hard limit on the size of the subspace it can manage, the architecture is built to be upgraded. The ability to run these calculations efficiently on existing hardware means that scientists can begin exploring complex quantum problems today, using the best of both classical and quantum computing to push the boundaries of what is calculable.

Technical Summary: A Scalable Sample-Based Quantum Diagonalization Solver for General Pauli-Sum Hamiltonians

Problem Statement
Sample-based quantum diagonalization (SQD) is a promising framework for solving ground-state problems of quantum many-body systems using near-term quantum computers. In SQD, a high-dimensional Hamiltonian is projected onto a small subspace defined by quantum states sampled from a quantum device. The challenge lies in the scale of these subspaces; while the quantum computer provides the basis states, classical post-processing can expand the subspace dimension to tens of billions (e.g., $10^{10}$ ).

Standard diagonalization methods fail in this regime because storing the reduced Hamiltonian matrix in memory is impossible unless it is extremely sparse. Consequently, a "matrix-free" approach is required, where the matrix-vector product is computed in situ without explicitly forming the matrix. Furthermore, these computations must be efficient enough to run on modern High-Performance Computing (HPC) facilities, which are increasingly GPU-centric. Existing open-source implementations either specialize in specific Hamiltonian types (e.g., chemistry or hard-core bosons) or rely on eigensolvers (like LOBPCG) that require excessive GPU memory retention, limiting their scalability for general Pauli-sum Hamiltonians.

Methodology
The author presents a scalable SQD solver built on the JAX API, designed to run on single workstations and many-GPU HPC environments. The solver comprises two primary components:

A Memory-Efficient Eigensolver:
- The solver employs the Locally Optimal Conjugate Gradient (LOC) method to find the lowest eigenpair of the reduced Hamiltonian.
- Unlike the more common Locally Optimal Block Preconditioned Conjugate Gradient (LOBPCG) method, which retains multiple large vectors to ensure stability, this implementation focuses solely on the lowest eigenpair.
- By aggressively reusing memory, the solver minimizes the GPU footprint, requiring only five $D$ -dimensional vectors simultaneously in RAM. This allows for the solution of eigenvalue problems with dimensions up to $10^9$ on a single GPU with 100 GB of memory.
- For larger problems, the solver supports distributed vectors across parallel nodes, utilizing JAX's built-in parallelization capabilities.
SQD-Specific Matrix-Vector Operation:
- The core operation computes $\beta = P H P \alpha$ , where $P$ is the projector onto the sampled subspace and $H$ is the general Pauli-sum Hamiltonian.
- The algorithm decomposes the Hamiltonian into Pauli strings and utilizes a symplectic representation to handle the Pauli operators.
- A critical step involves an "X source index search" (Algorithm 1), which identifies the mapping between basis states under the action of Pauli $X$ operators. This is implemented using stable sorting on concatenated arrays of binary indices, optimized for GPU architectures.
- To optimize performance, the algorithm groups Hamiltonian terms by their $X$ signatures, reducing the number of expensive search and "take" operations.
- The implementation includes caching mechanisms for static data (source indices and diagonal coefficients) to trade off memory usage for execution speed.

Key Contributions

General-Purpose Implementation: Unlike previous libraries restricted to specific Hamiltonian types, this solver handles general Pauli-sum Hamiltonians with no inherent limit on the number of qubits.
GPU-Native Scalability: The code is written in Python using JAX, enabling seamless execution on single GPUs and distributed many-GPU HPC clusters with minimal code changes.
Memory Optimization: The custom eigensolver is explicitly designed to minimize GPU memory usage, enabling the handling of dimensions ( $D \approx 10^9$ ) that would exceed the capacity of standard block-preconditioned methods on current hardware.
Open Source Availability: The solver is released as an open-source library (rqutils), providing a general utility for finding eigenpairs of any Hermitian operator where the matrix-vector operation is implemented as a JAX JIT-compiled function.

Results
The paper validates the solver through two primary experiments:

Memory Footprint Validation: The solver successfully computed the ground and first excited states of a 30-site transverse-field Ising model (TFIM) on a single NVIDIA H100 GPU (96 GB RAM). This demonstrated the solver's ability to handle non-trivial Hamiltonians within strict memory constraints.
Strong Scaling Performance: The solver was tested on a 32-site TFIM (which exceeds single-GPU memory capacity) across 8, 16, 32, 64, and 128 GPU nodes. The results showed good strong scaling, with the time-to-solution decreasing effectively as the number of nodes increased, while maintaining identical energy values across all instances.

Significance and Limitations
The paper positions this solver as a critical tool for "quantum-centric supercomputing," enabling the classical post-processing required to extract quantum advantage from pre-fault-tolerant quantum devices. By separating the quantum role (suggesting subspace dimensions) from the classical role (diagonalization and error correction), the framework leverages massive classical compute power to overcome current quantum hardware limitations.

However, the author explicitly notes a current limitation: the X source index search relies on global stable sorting, which is currently limited by JAX to arrays of size $2^{31}-1$ . Consequently, the solver currently supports a maximum subspace dimension of $N = 2^{30}-1$ . The author acknowledges that while this covers many current use cases, future upgrades to block-wise sorting are planned to support the "tens of billions" of dimensions anticipated in cutting-edge SQD applications. The paper does not claim to solve the fault-tolerance problem but rather to maximize the utility of near-term quantum computers through efficient classical scaling.

Technical Summary: A Scalable Sample-Based Quantum Diagonalization Solver for General Pauli-Sum Hamiltonians

More like this