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

YUTARO IIYAMA^✉, International Center for Elementary Particle Physics, The University of Tokyo, Japan

A scalable ground-state solver to be used in sample-based quantum diagonalization (SQD) is presented. The solver is built on JAX API, and is intended to be run on GPU platforms and possibly many-GPU high-performance computing facilities. It is capable of subspace diagonalization in up to $2^{30}\text{-}$dimensional subspace, with possible upgrades to enable even larger-scale operations. The two main components, a general-purpose eigensolver and an SQD-specific matrix-vector operation function, are described in detail. Efficiency of the eigensolver with respect to the GPU memory is demonstrated using a simple physics model. Furthermore, good strong scaling of the eigensolver in many-GPU environment is observed.

CCS Concepts: • Applied computing $\to$ Physics; • Software and its engineering $\to$ Software libraries and repositories.

Additional Key Words and Phrases: Subspace quantum diagonalization, Quantum-centric supercomputing, Open source

ACM Reference Format:

Yutaro Iiyama. 2026. A scalable sample-based quantum diagonalization solver for general Pauli-sum Hamiltonians. In Proceedings of Make sure to enter the correct conference title from your rights confirmation email (Conference acronym ’XX). ACM, New York, NY, USA, 10 pages. https://doi.org/XXXXXXX.XXXXXXX

1 Introduction

Sample-based quantum diagonalization (SQD) has emerged as one of the most promising application of near-term quantum computers[6, 9–12]. Regarded as a superior alternative to the variational quantum eigensolver[8], SQD has been applied to a wide variety of ground-state problems of quantum many-body systems, demonstrating energy approximation accuracies comparable to the best classical methods.

In SQD, a many-body Hamiltonian of an intractable dimensionality is projected onto a small subspace, in which the ground-state eigenpair of the reduced Hamiltonian is found. The projection subspace dimensions are identified through sampling quantum states realized on a quantum computater. In the original formulation, which is also known as quantum selected configuration interaction (QSCI), these states are the approximations of the ground state prepared with ansatz circuits. In the variant dubbed sample-based Krylov quantum diagonalization (SKQD), final states of time evolutions under the target Hamiltonian for multiple durations are sampled instead. Crucially, in all versions of SQD, the only role quantum computers play is to suggest the subspace dimensions. Everything else, including correcting potential errors in the obtained samples, is performed as classical computation. This separation of responsibility between quantum and classical computation, with a strong emphasis on the latter, is the innovation that allows SQD to be a candidate framework in which quantum advantage is achieved even with pre-fault-tolerant quantum computers.

It is of critical importance that “small” is a relative term in this context; in practice, the number of subspace dimensions can reach up to tens of billions. Even when the number of measurement outcomes from the quantum computer are less than a million, classical postprocessing can drastically enlarge the subspace size by e.g. identifying additional relevant dimensions through symmetry considerations. The implication is that the Hamiltonian projection and diagonalization is not performed sequentially, because it is impossible to store the reduced Hamiltonian matrix in memory unless it is known to be extremely sparse. Instead, a matrix-free method must be employed for diagonalization, in which the matrix-vector operation, which computes the result of multiplying a vector with the reduced Hamiltonian, must perform the projection in situ. Additionally, this computation would likely be executed in distributed parallel environments such as high-performance computing (HPC) facilities.

The central software components for large-scale SQD are therefore a matrix-free eigensolver and a specialized matrix-vector operation. Since both happen to involve many vector-vector and matrix-vector multiplications, they are amenable to GPU acceleration. In fact, as modern HPC architectures provide the majority of their compute through GPUs, it is essentially necessary that these components are implemented in a GPU-native manner.

There already exist at least two such implementations that are open source. The first is SBD-G[4], which is a header-only C++ library utilizing Thrust[3] and MPI to enable multi-GPU parallelism. This library specializes in SQD of chemistry and hard-core boson Hamiltonians, and has achieved state-of-the-art performance of diagonalizing a $1.36\times 10^{18}\text{-}$dimensional Hamiltonian projected onto a $2.34\times 10^{10}\text{-}$dimensional subspace using 128 GPU nodes in Supercomputer Miyabi[5].

The second implementation, which is more similar in scope and design to what is presented in this paper, is in Ref. [2] as a part of the demonstration of the functionalities of the CUDA Quantum library. Unlike SBD-G, the Hamiltonian it works with is expressed as a general Pauli sum of up to 64 qubits. The eigensolver is based on the thick-restarted Lanczos method, which requires retention of many large vectors in memory.

Our SQD solver, which is also open source, is written in Python and built on the JAX API[1]. Like the CUDA Quantum implementation, it targets general Pauli-sum Hamiltonians, but there is no limit to the number of qubits the Hamiltonians act on. The eigensolver employs the locally optimal conjugate gradient method, and is optimized for GPU memory footprint. By taking advantage of the built-in parallelization capabilities in JAX, the full solver can be run in many-GPU HPC as well as a single workstation, with minimal change in the user code.

We designed the library to make the eigensolver a standalone general utility which can be used to find the lowest few eigenpairs of any Hermitian operator, as long as its operation on a vector is implemented as a JAX just-in-time compiled function. This eigensolver can be scaled up to eigenvalue problems of ${10}^{10}$ or higher dimensions. The SQD solver is a downstream user of this eigensolver, providing the SQD-specific matrix-vector operator constructed from the problem Hamiltonian and the subspace dimensions. This matrix-vector operator currently throttles the scaling of the solver, as it has a hard limit of $2^{31}$ for the number of subspace dimensions that can be handled. An upgrade to remove this restriction is planned.

This paper presents our SQD solver with the goal of broadening its user base by clarifying its internal workings and demonstrating its capabilities. We describe the conjugate gradient eigensolver and matrix-vector operation in Section 2. The performance and scaling of the standalone eigensolver and the overall SQD solver are discussed in Section 3. Section 4 concludes with comments on the planned extensions.

2 SQD components

2.1 Eigensolver

The locally optimal conjugate gradient method for an eigenvalue problem

computes a series of solution estimates $\{v^{(0)},v^{(1)},...\}$ iteratively, with progressively better approximations. Here, $A$ is an $D\times D$ symmetric or Hermitian matrix and ${\lambda }_0$ is its lowest eigenvalue. At each iteration, the new estimate is obtained as the minimizer of the Rayleigh quotient of $A$ within the three-dimensional subspace spanned by the current and previous estimates and the current residual. Derivation and detailed discussion of the algorithm, which has been known for several decades, can be found elsewhere[7]. We also provide our implementation as a pseudocode in Appendix A.

This algorithm can be regarded as a single-vector unpreconditioned version of the better-known locally optimal block preconditioned conjugate gradient (LOBPCG) method, which can identify multiple lowest or highest eigenpairs of $A$ simultaneously. We chose to create our own eigensolver, even though there exist multiple public implementations of LOBPCG, in order to minimize the GPU memory footprint. By focusing on the lowest eigenpair only and aggressively reusing memory, our eigensolver can complete its calculation with as low as five $D\text{-}$dimensional vectors simultaneously materialized in the GPU RAM. This means that Equation (18) with $D={10}^9$ and $v$ represented as an array of 128-bit floating point complex numbers may be solved on a single GPU with a $100\mathrm{G}\mathrm{B}$ memory. Note however that the actual maximum number of vectors in RAM can be larger, depending on the matrix-vector operation function that encodes $v\mapsto Av\text{.}$

When tackling a higher-dimensional problem or using a matrix-vector function that requires a large RAM allocation, the eigensolver can be executed in parallel computing environments by splitting the $D\text{-}$dimensional vectors into chunks that are distributed among the parallel nodes. There is thus in principle no limit to the problem dimension the eigensolver can handle, although in reality the maximum number of parallel GPU nodes available at HPC centers are typically at a few hundreds, implying a dimensionality limit of about $10^{10}\text{.}$ For a fixed problem size, the eigensolver is expected to show a good strong scaling, as the only operations it performs on the vectors outside the matrix-vector function are element-wise additions, scalar multiplications, and vector inner products, all of which entail minimal data transfer between the parallel nodes.

It should be mentioned that we do also provide an option to compute the second and higher eigenpairs by passing to the eigensolver a list of vectors that the solution should be orthogonalized against. This crude approach introduces minimal additional memory footprint but is known to suffer from numerical instability, which is exactly why LOBPCG was invented. For SQD, this functionality is not used at any rate.

2.2 Reduced Hamiltonian matrix-vector operation

2.2.1 Algorithm.

We now discuss the SQD-specific matrix-vector operation function. First, we must clarify what exactly correspond to the $A$ and $v$ of Equation (18) in the case of SQD. Let $B=\{{\mathbf{b}}_n{\}}_{n=1}^N$ be the $q\text{-}$bit binary indices (bit strings) of the subspace dimensions identified by quantum states sampling and postprocessing. Using the projector onto this subspace

the subspace diagonalization of the $q\text{-}$qubit Hamiltonian $H$ is formulated as the process of identifying $\alpha :=({\alpha }_n{)}_{n=1}^N$ and ${\tilde{E}}_0$ in

The eigenvalue problem to solve is therefore

which means that the matrix-vector operation must compute

The Hamiltonian $H$ is a Pauli sum

where each Pauli string $P_k\in \{I,X,Y,Z{\}}^{\otimes q}$ is a tensor product of $q$ Pauli operators. We break down Equation (5) using (6) and compute each $⟨{\mathbf{b}}_m|P_k|{\mathbf{b}}_n⟩$ separately. To do so, $P_k$ is expressed in the symplectic representation

where ${\mathbf{x}}_k=(x_{k,j}{)}_{j=0}^{q-1},{\mathbf{z}}_k=(z_{k,j}{)}_{j=0}^{q-1}\in \{0,1{\}}^q$ are the X and Z signatures of $P_k\text{.}$ The phase factor $(-i{)}^{{\mathbf{x}}_k\cdot {\mathbf{z}}_k}$ is static and can be combined with the Pauli coefficient into

Looking at the matrix elements of the operators, because $Z$ is diagonal in the computational basis and $X$ maps one basis state to another, we have

Therefore, the projection of the Hamiltonian boils down to identifying a ${\mathbf{b}}_{n_m}$ for each ${\mathbf{b}}_m$ such that

There is at most one such index in $B\text{.}$

We use Algorithm 1 (X source index search) to find $n_{k,m}$ for each $m\text{.}$ This is an application of a search algorithm designed for fixed tensor size architectures (e.g. GPU). The returned vector of indices $(n_{k,m}{)}_{m=1}^N$ may contain invalid entries, indicating that bit strings satisfying Equation (10) do not exist in $S$ for those slots.

Algorithm 1 X source index search. In the lines, StableSort refers to a sorting that keeps the original relative order of elements that have equal values, and $\mathrm{a}\mathrm{r}\mathrm{g}\text{StableSort}$ returns the indices into the original array which would result in the sorted array.

Require: Sorted list of binary indices $S=({\mathbf{b}}_m{)}_{m=1}^N\text{,}$ binary X signature $\mathbf{x}$

Ensure: List of indices $(n_m{)}_{m=1}^N\text{,}$ taking values in $\{0,...,N\}\text{,}$ such that ${\mathbf{b}}_{n_m}=\mathbf{x}\oplus {\mathbf{b}}_m$ if $n_m/=0\text{.}$

$$ Concatenate $S$ and $\mathbf{x}\oplus S$ and take the sorting indices

• $T\leftarrow (\mathbf{x}\oplus {\mathbf{b}}_m{)}_{m=1}^N$

• $U\leftarrow S\Vert T$

• $P\leftarrow \text{StableSort}(U)$

• $Q\leftarrow \mathrm{a}\mathrm{r}\mathrm{g}\text{StableSort}(U)$

$$ Collect $(n_m{)}_{m=1}^N$ from $Q$ by checking repetitions in $P$

• for all $m\in \{1,...,N\}$ do

• $$ $k_m\leftarrow k\in \{1,...,2N\}\text{s.t.}Q[k]=m$

• $$ if $P[k_m+1]=P[k_m]$ then

• $$ $n_m\leftarrow Q[k_m+1]-N$

• $$ else

• $$ $n_m\leftarrow 0$

• $$ end if

• end for

• return $(n_m{)}_{m=1}^N$

Plugging Equation (6)-(9) back into Equation (5), we get

where

In practice, the “take” function is used on $\alpha$ for indices $n_{k,m}$ to obtain the right-hand side of Equation (12), with a fill-in value of 0 for invalid indices.

Because X source index search and the “take” function are both expensive operations, the latter especially when the $\alpha$ vector is distributed, we seek to minimize their usage by grouping the Hamiltonian terms by the X signature. Equation (6) and (7) combined is rewritten as

where $\{{\mathbf{x}}_{(\ell )}{\}}_{s=1}^L$ are the $L$ distinct X signatures found in the Hamiltonian, and $G_{\ell }$ is the set of indices for the term in the original breakdown (6) with the X signature ${\mathbf{x}}_{(\ell )}\text{.}$ With this grouping, we compute the X source indices $\{n_{(\ell ),m}{\}}_{m=1}^N$ for each ${\mathbf{x}}_{(\ell )}\text{.}$ Equation (11) is then rewritten as

with

2.2.2 Caching static data.

The matrix-vector operation of Section 2.2.1 is called hundreds of times during the eigensolver execution in some cases. Only the $\alpha$ vector changes between the invocations, while the diagonals

and the source indices $\{n_{(\ell ),m}{\}}_{\ell =1,m=1}^{L,N}$ stay constant. Therefore, computing these quantities once and reusing them can speed up the diagonalization process, with a potential cost in GPU memory. We thus provide such caching options in our matrix-vector function.

If each $d_{(\ell ),m}$ is represented with a 128-bit complex number, storing the diagonals in memory would require $16LN$ bytes. For the source indices, the memory requirement would be $4LN$ bytes, if 32-bit integers are used (which is actually necessary, as explained in the following). Depending on the number of distinct X signatures $L\text{,}$ these costs can be prohibitive. Given that X source finding is significantly expensive compared to all other operations in the diagonalization process, the default setting is to cache the X source indices but compute the diagonals on the fly at each invocation of the matrix-vector function. Note, though, that if the diagonals are also cached, the list of dimension indices $S$ will no longer be necessary to be kept in memory, freeing up $\lceil \frac{q}{8}\rceil N$ bytes. Therefore, if $q>128L\text{,}$ it is strictly advantageous to cache also the diagonals.

2.2.3 Limitation on the number of subspace dimensions.

As evident from Algorithm 1, X source search is heavily reliant on stable sorting of subspace dimension indices. Since stable sorting of more than $2^{31}-1$ elements is currently not implemented in JAX, while Algorithm 1 involves sorting a concatenated array of $2N$ elements, the maximum number of subspace dimensions $N$ is currently $2^{30}-1\text{.}$

However, global sorting of the binary indices is not a strict necessity. When the $\alpha$ vector is distributed among parallel compute nodes, we can foresee defining separate X source indices for each node, in which case Algorithm 1 would be run on chunks of $S$ in parallel. Although using such indices can substantially slow down the function execution, we will be able to work around the current limitation this way.

3 Performance of the SQD solver

3.1 Eigensolver

To validate the conjugate gradient eigensolver, we compute the two lowest energy eigenvalues $E_0,E_1$ of the one-dimensional transverse-field Ising model (TFIM) with $q=30$ sites and a periodic boundary condition. Explicitly, the Hamiltonian of the model is

where $Z_q$ is identified with $Z_0$ and $h$ is the transverse field parameter. This model is known to exhibit a phase transition at $h=1$ in the thermodynamic limit, in which the ground state is doubly degenerate in the ordered $\text{(}h<1\text{)}$ phase and is unique in the disordered $\text{(}h>1\text{)}$ phase.

Figure 1a shows the energy gap $E_1-E_0$ as a function of $h\text{.}$ Lifting of the degeneracy at $h=1$ is evident already at this system size. These calculations were performed on a single NVIDIA H100 $\text{(}96\mathrm{G}\mathrm{B}$ RAM) node, thus demonstrating that the eigensolver has the expected memory footprint, and is capable of solving the ground and first excited states of nontrivial Hamiltonians (TFIM with $h\sim 1\text{).}$

Next, we verify that the eigensolver can run with distributed vectors, and investigate its time-to-solution scaling as the number of parallel nodes is increased. We consider the same TFIM Hamiltonian at $h=0.5$ with $q=32$ sites, which cannot be diagonalized on a single H100 GPU any more. The problem is then solved using 8, 16, 32, 64, and 128 nodes in parallel. Exactly identical energy values are obtained in all instances. Figure 1b plots the time to completion for computing the ground eigenpair for different parallelizations. Good strong scaling is observed.

(a) 30-site TFIM energy gap as a function of $h\text{.}$

(b) Time for 32-site TFIM ground state energy determination.

Fig. 1. Validation of the conjugate gradient eigensolver

4 Conclusion

We presented a scalable, memory-efficient SQD solver for general Pauli Hamiltonians based on JAX. The solver consists of a general-purpose conjugate gradient eigensolver and an SQD-specific matrix-vector multiplication function that performs Hamiltonian projection internally. Both the eigensolver and the matrix-vector function support distributed vectors, allowing the overall solver to run in massively parallel GPU environments to leverage modern HPC architectures. The matrix-vector function can operate in several different caching modes, allowing control over the tradeoff between the execution time and memory footprint.

The SQD solver is currently limited to $2^{30}$ dimensions due to the need to sort the states array globally. To support the cutting-edge SQD use case with tens of billions of subspace dimensions, an upgrade to block-wise sorting is foreseen.

Acknowledgments

This work was supported by Japan Society for the Promotion of Science KAKENHI Grant Numbers JP24K07042 and JP24H00689.

Code availability

The SQD solver code is available at https://github.com/UTokyo-ICEPP/rqutils/blob/main/rqutils/sqd.py with a corresponding documentation at https://rqutils.readthedocs.io/en/latest/apidoc/rqutils.sqd.html.

References

[1] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Yash Katariya, Chris Leary, Dougal Maclaurin, George Necula, Adam Paszke, Jake VanderPlas, Skye Wanderman-Milne, and Qiao Zhang. 2018. JAX: composable transformations of Python+NumPy programs. http://github.com/jax-ml/jax

[2] NVIDIA Corporation. 2026. CUDA-Q Applications: Sample-Based Krylov Quantum Diagonalization (SKQD). https://nvidia.github.io/cuda-quantum/latest/applications/python/skqd.html.

[3] NVIDIA Corporation. 2026. Thrust. https://developer.nvidia.com/thrust.

[4] Jun Doi, Tomonori Shirakawa, Yukio Kawashima, Seiji Yunoki, and Hiroshi Horii. 2026. GPU-accelerated selected Basis Diagonalization with Thrust for SQD-based algorithms. arXiv:2601.16637 [cs.DC] doi:10.48550/arXiv.2601.16637

[5] Toshihiro Hanawa, Kengo Nakajima, Yohei Miki, Takashi Shimokawabe, Kazuya Yamazaki, Shinji Sumimoto, Osamu Tatebe, Taisuke Boku, Daisuke Takahashi, Akira Nukada, Norihisa Fujita, Ryohei Kobayashi, Hiroto Tadano, and Akira Naruse. 2024. Preliminary performance evaluation of Grace-hopper GH200. In 2024 IEEE International Conference on Cluster Computing Workshops (CLUSTER Workshops) (2024-09-24). IEEE, 184–185. doi:10.1109/clusterworkshops61563.2024.00050

[6] Keita Kanno, Masaya Kohda, Ryosuke Imai, Sho Koh, Kosuke Mitarai, Wataru Mizukami, and Yuya O Nakagawa. 2026. Quantum-selected configuration interaction: Classical diagonalization of Hamiltonians in subspaces selected by quantum computers. 8 (2026), 023268. Issue 2. doi:10.1103/dmn4-snfx

[7] Andrew V Knyazev. 1990. A preconditioned conjugate gradient method for eigenvalue problems and its implementation in a subspace. In Numerical Treatment of Eigenvalue Problems Vol. 5 / Numerische Behandlung von Eigenwertaufgaben Band 5: Workshop in oberwolfach, February 25 - march 3, 1990 / tagung in oberwolfach, 25. Februar - 3. März 1990 (1990-02-25) (International Series of Numerical Mathematics), J Albrecht, L Collatz, P Hagedorn, and W Velte (Eds.). Springer, 143–154. doi:10.1007/978-3-0348-6332-2

[8] Alberto Peruzzo, Jarrod McClean, Peter Shadbolt, Man-Hong Yung, Xiao-Qi Zhou, Peter J Love, Alán Aspuru-Guzik, and Jeremy L O’Brien. 2014. A variational eigenvalue solver on a photonic quantum processor. 5 (2014), 4213. Issue 1. doi:10.1038/ncomms5213

[9] Javier Robledo-Moreno, Mario Motta, Holger Haas, Ali Javadi-Abhari, Petar Jurcevic, William Kirby, Simon Martiel, Kunal Sharma, Sandeep Sharma, Tomonori Shirakawa, Iskandar Sitdikov, Rong-Yang Sun, Kevin J Sung, Maika Takita, Minh C Tran, Seiji Yunoki, and Antonio Mezzacapo. 2025. Chemistry beyond the scale of exact diagonalization on a quantum-centric supercomputer. 11 (2025), eadu9991. Issue 25. doi:10.1126/sciadv.adu9991

[10] Tomonori Shirakawa, Javier Robledo-Moreno, Toshinari Itoko, Vinay Tripathi, Kento Ueda, Yukio Kawashima, Lukas Broers, William Kirby, Himadri Pathak, Hanhee Paik, Miwako Tsuji, Yuetsu Kodama, Mitsuhisa Sato, Constantinos Evangelinos, Seetharami Seelam, Robert Walkup, Seiji Yunoki, Mario Motta, Petar Jurcevic, Hiroshi Horii, and Antonio Mezzacapo. 2025. Closed-loop calculations of electronic structure on a quantum processor and a classical supercomputer at full scale. arXiv:2511.00224 [quant-ph] doi:10.48550/arXiv.2511.00224

[11] Kenji Sugisaki, Shu Kanno, Toshinari Itoko, Rei Sakuma, and Naoki Yamamoto. 2025. Hamiltonian simulation-based quantum-selected configuration interaction for large-scale electronic structure calculations with a quantum computer. 27 (2025), 20869–20884. Issue 38. doi:10.1039/d5cp02202a

[12] Jeffery Yu, Javier Robledo Moreno, Joseph T Iosue, Luke Bertels, Daniel Claudino, Bryce Fuller, Peter Groszkowski, Travis S Humble, Petar Jurcevic, William Kirby, Thomas A Maier, Mario Motta, Bibek Pokharel, Alireza Seif, Amir Shehata, Kevin J Sung, Minh C Tran, Vinay Tripathi, Antonio Mezzacapo, and Kunal Sharma. 2025. Quantum-Centric Algorithm for Sample-Based Krylov Diagonalization. arXiv:2501.09702 [quant-ph]

A Locally-optimal conjugate gradient

Given a symmetric or Hermitian matrix $A\text{,}$ the locally optimal conjugate gradient method for the standard minimum-eigenvalue problem

is an iterative procedure where the estimate solutions $\{x^{(0)},x^{(1)},...\}$ are formally refined through

with search (conjugate) directions $\{p^{(i)}{\}}_i$ and optimal step sizes $\{{\alpha }^{(i)}{\}}_i\text{.}$ The search directions are linear combinations of the gradient of the Rayleigh quotient

evaluated at $x\in \{x^{(i)}{\}}_i\text{,}$ and are recursively expressed as

where

and

using some $\{{\beta }^{(i)}{\}}_i\text{.}$

In Equation (19), the step size ${\alpha }^{(i)}$ is optimal when $\rho (x^{(i+1)})$ is minimized. Furthermore, from Equations (19) and (21), we find

Therefore, the solution update at iteration $i$ is equivalently obtained as

The right-hand side of Equation (25) is found by the Rayleigh-Ritz method.

Algorithm 2 is the pseudocode representation of our eigensolver implementation. It includes the details on the treatment of the first two iterations, which is ommitted in the discussion above.

Algorithm 2 Locally optimal conjugate gradient eigensolver

Require: Hermitian operator $A$ (via $x\mapsto Ax\text{),}$ initial vector $x_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}\in {\mathbb{C}}^n\text{,}$ maximum iterations $N_{\max }\text{,}$ tolerance $\tau$

Ensure: Lowest eigenvalue $\theta \text{,}$ eigenvector $x\text{,}$ number of iterations $k$

• function $\text{ProjectOut(}(b_1,...,b_m),v\text{)}$

• $$ for $i=0,1,2,3$ do

• $$ $c_j\leftarrow b_j^{†}v$ for $j=1,...,m$

• $$ $v\leftarrow v-{\sum }_j{c}_j{b}_j$

• $$ if $i<2$ and $\Vert v\Vert /=0$ then

• $$ $v\leftarrow v/\Vert v\Vert$

• $$ end if

• $$ end for

• $$ if $\Vert v\Vert <0.99$ then

• $$ $v\leftarrow 0$ $▷$ Discard numerically unreliable direction

• $$ end if

• $$ return $v$

• end function

• function $\text{RayleighRitz(}{w}_1,...,w_m\text{)}$ $▷$ $m\in \{2,3\}$

• $$ $S_{jk}\leftarrow w_j^{†}A{w}_k$ for $j,k=1,...,m$ $▷$ $m\times m$ Hermitian

• $$ $(\theta ,\kappa )\leftarrow$ lowest eigenpair of $S\text{,}$ $\Vert \kappa \Vert =1$ $▷$ Analytic for $2\times 2$ / Cardano for $3\times 3$

• $$ return $(\theta ,\kappa )$

• end function

$$ Iteration 0: steepest-descent residual

• $x\leftarrow x_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}/\Vert x_{\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}}\Vert$

• $\rho \leftarrow x^{†}Ax$

• $r\leftarrow Ax-\rho x$

$$ Iteration 1: two-vector Rayleigh–Ritz

• $p\leftarrow r/\Vert r\Vert$ if $\Vert r\Vert /=0$ else $r$

• $(\theta ,\kappa )\leftarrow \text{RayleighRitz}(x,p)$

• $t\leftarrow {\kappa }_{0}p-{\kappa }_{1}x$

• $u\leftarrow {\kappa }_{0}x+{\kappa }_{1}p$

• $x\leftarrow u/\Vert u\Vert$

• $y\leftarrow t/\Vert t\Vert$

• $r\leftarrow Ax-\theta x$

$$ Main LOCG loop

• $k\leftarrow 0\text{,}$ $\mathit{c}\mathit{o}\mathit{n}\mathit{v}\mathit{e}\mathit{r}\mathit{g}\mathit{e}\mathit{d}\leftarrow \mathbf{\text{false}}\text{,}$ $\theta \leftarrow 0$

• while $k<N_{\max }$ and not $\mathit{c}\mathit{o}\mathit{n}\mathit{v}\mathit{e}\mathit{r}\mathit{g}\mathit{e}\mathit{d}$ do

• $$ $p\leftarrow \text{ProjectOut}((x,y),r)$

• $$ $(\theta ,\kappa )\leftarrow \text{RayleighRitz}(x,y,p)$

• $$ $s\leftarrow {\kappa }_{1}y+{\kappa }_{2}p$

• $$ $t\leftarrow ({\kappa }_0/\Vert s\Vert )s-\Vert s\Vert x$ $▷$ Orthogonal to new $x\text{,}$ spans old $\{x,s\}$

• $$ $u\leftarrow {\kappa }_{0}x+s$

• $$ $x\leftarrow u/\Vert u\Vert$

• $$ $y\leftarrow t/\Vert t\Vert$

• $$ $v\leftarrow Ax$

• $$ $r\leftarrow v-\theta x$

• $$ $\delta \leftarrow 10n(\Vert v\Vert -\theta )$ $▷$ Relative tolerance scale

• $$ $\mathit{c}\mathit{o}\mathit{n}\mathit{v}\mathit{e}\mathit{r}\mathit{g}\mathit{e}\mathit{d}\leftarrow \Vert r\Vert <\tau \delta$

• $$ $k\leftarrow k+1$

• end while

• return $(\theta ,x,k)$

Author’s Contact Information: Yutaro Iiyama, iiyama@icepp.s.u-tokyo.ac.jp, International Center for Elementary Particle Physics, The University of Tokyo, Bunkyo, Tokyo, Japan.

Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than the author(s) must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and /or a fee. Request permissions from permissions@acm.org.

© 2026 Copyright held by the owner/author(s). Publication rights licensed to ACM.
Manuscript submitted to ACM

About this conversion

“A scalable sample-based quantum diagonalization solver for general Pauli-sum Hamiltonians” by Yutaro Iiyama. Original arXiv submission, 2610.13907v1. Shared under Creative Commons Attribution 4.0 International, as identified by the exact-version arXiv record.

Linneon LLC converted the original LaTeX submission to continuous Markdown and HTML using LEC on October 14, 2026. Pagination and page furniture are removed; tables, notes and figures are arranged for continuous reading. HTML mathematics uses MathJax. These reading copies are not the official publication and do not imply endorsement by the authors, their institutions or arXiv. Original notices are retained. Dates generated by the source may reflect recompilation.