Terastate-per-second QUBO Brute-Force on a Single GPU: A Matrix Prefix-Suffix Decomposition Aleksandr Maltsev, Mikhail Remnev, Alexey Kapranov and Ekaterina Krivtsova
arXiv:2607.04857v1 [cs.DS] 6 Jul 2026
Correspondence to: A. Maltsev, Email: [email protected].
Abstract—This paper presents a parallel QUBO exhaustive search algorithm for dense matrices, based on a prefix-suffix decomposition and Gray code ordering. The algorithm achieves O(1) per-state complexity: for the QUBO objective function computation only one arithmetic operation per state is performed. An adjustable energy components cache size enables placement in the fastest available memory tier. This reduces memory bandwidth requirements to a negligible level and transforms the problem from memory-bound to compute-bound. Our CUDA-based implementation achieves a state-of-the-art evaluation rate of 7.5 × 1012 states per second on a single GPU, setting a new performance benchmark for the fullspace-search subclass of exact solvers. Index Terms—Quantum computing, Ising model, Ising spin–glass, QUBO, Combinatorial optimization, GPU, CUDA, Exhaustive search, Bruteforce, Quantum annealers, HPC
I. Introduction The QUBO (Quadratic Unconstrained Binary Optimization) problem involves minimizing the function X E(x) = Qij xi xj , i,j
where x ∈ {0, 1}n , and is one of the central problems in combinatorial optimization [1–3]. This problem is equivalent to finding the ground state of the Ising model [1, 4], described by the Hamiltonian X X H(s) = Jij si sj + hi si . i,j
i
A wide range of combinatorial optimization problems can be reduced to the QUBO formulation, which underlies its practical importance. These include graph-theoretic problems (MaxCut, graph coloring, maximum independent set,
community detection in networks), logistic problems such as the traveling salesman problem (TSP) with N 2 variables for solving an N -city problem, and financial applications (portfolio optimization) [3, 5–7]. Machine learning tasks (inference in graphical models, restricted Boltzmann machines), industrial applications (optimization of metal heat treatment schedules), as well as physical simulations (molecular modeling and ground state search in spin glasses) are also worth mentioning [1, 2, 4, 8–10]. Furthermore, a systematic transformation from general 0/1 linear programming models to QUBO formulations is available [5]. The practical significance of these formulations is further strengthened by the fact that the Ising Hamiltonian is naturally expressible in the language of quantum computation, where binary spin variables are directly associated with qubit degrees of freedom and pairwise couplings can be encoded as quantum interactions. As a result, QUBO and Ising representations provide a broadly applicable route for translating realworld optimization tasks into forms amenable to quantum algorithms and quantum hardware, including adiabatic and gate-based approaches [7]. The main reason for the difficulty of QUBO problems lies in their NP-hardness [1–4]. Even with just one additional bit in problem size, the state space doubles, reaching 2n possible configurations [2]. Without special graph structure, finding an exact solution in polynomial time is impossible [1]. This complexity is worsened by the presence of numerous local minima in the energy landscape [11, 12]. Fully-connected problems have non-zero Qij elements for all variable pairs, creating maximum
computational complexity [13]. For example, a fully-connected problem on 100,000 spins contains about 5 billion edges [11]. Sparse problems contain connections only between nearest neighbors or on specific topologies, such as the D-Wave Chimera graph, where each node connects to 5–6 others [1, 4]. A. Exact Solution Methods Branch and Bound is one of the most powerful classical tools for exactly solving NP-hard problems [4, 14, 15]. The algorithm recursively partitions the solution space into subproblems and uses relaxations to obtain lower bounds on the cost of the optimal solution [4, 15]. For QUBO problems, odd-cycle inequalities, semidefinite programming (SDP), and dynamic programming methods for graphs with bounded treewidth are used [4, 15, 16]. The SDP-based Branch-and-Cut methodology demonstrates successful results in solving Chimera graph problems with guaranteed optimal solutions [4, 15]. Despite its power, Branch-andBound is practically applicable only to problems up to ≈100 variables in the fully-connected case [17, 18]. For sparse matrices, problems with 10,000 nodes can be solved [17]. While QUBO formulations offer universality for solving combinatorial optimization problems on quantum and quantum-inspired hardware, direct problem-specific approaches often prove more efficient for large-scale instances. A notable example is the exact solution of the Traveling Salesman Problem with 85,900 cities [19, 20], achieved through specialized branch-and-cut methods implemented in the Concorde solver, demonstrating that tailored algorithms can handle problem sizes far beyond the reach of current QUBO-based approaches while providing rigorous certification of optimality. Tensor network methods enable exact ground state energy computation by compressing the tensor network using tropical algebra, successfully applied to problems up to 1024 spins on the Chimera graph [9]. Industrial solvers CPLEX and Gurobi are standard tools using integer programming [10, 21]. Specialized exact solvers like BEIT
for Chimera topology can handle problems up to 1024 variables with optimality guarantees [1]. B. Heuristic Methods Simulated Annealing remains a classical method but requires sequential variable updates for its core search trajectory [21, 22], which limits its degree of parallelism. Newer physicsinspired algorithms show advantages: Simulated Bifurcation enables simultaneous updates of all variables [11, 12], while Momentum Annealing achieves 250x speedup compared to Simulated Annealing (SA) on GPU for fully-connected problems on 100,000 spins [13]. GPU-based solvers demonstrate superior performance for QUBO problems through massive parallelization, as exemplified by the DABS (Diverse Adaptive Bulk Search) framework that achieves 100% success rate on benchmark instances including K2000 and scales efficiently across 8 NVIDIA A100 GPUs, solving problems with up to 10,000 variables in seconds—100200 times faster than quantum annealers while maintaining guaranteed optimality [6]. A class of specialized hardware solvers is also developing, using physical processes to find the energy minimum: memristor crossbars with quantum-inspired parallel annealing (QPA), coherent Ising machines (CIM), bifurcation algorithms on FPGA and GPU, oscillator systems, and magnetic devices [23–25]. Classical CMOS implementations of fully connected problems include STATICA, a 512-spin annealing processor fabricated in 65-nm technology that enables parallel spin updates through stochastic cellular automata dynamics and achieves state-of-the-art performance on combinatorial optimization [26]. Coherent Ising machines deserve special mention, as they enable solutions to fully connected problems with up to 100,000 variables [27], though not optimally. Machine learning approaches, particularly graph neural networks and deep reinforcement learning, have demonstrated remarkable potential for solving large-scale QUBO and Ising model problems. Physics-inspired graph neural networks achieve near-optimal solutions for instances with up to one million variables, demonstrating supe-
rior scalability compared to quantum annealers while maintaining approximation ratios above 90% for Max-Cut problems [28], while deep reinforcement learning methods have successfully found ground states of 3D spin glass systems with up to 8,000 spins—requiring 540 times fewer initial configurations than parallel tempering to guarantee exact solutions for smaller systems [29]. Importantly, both the Max-Cut and 3D spin glass problems are equivalent to QUBO formulations, as Max-Cut can be directly expressed as a QUBO.
Gurobi proved fastest with linear scaling and 100% optimality, while D-Wave Leap NL HQA (Nonlinear Hybrid Quantum Annealer) achieved optimality with comparable performance, outperforming Fujitsu DA v4, which had the worst performance. For extra-large proteins, Gurobi maintained optimal solutions, while NL HQA showed suboptimal performance. A critically important result was the identification of rank-1 dominance: guanine-cytosine content quadratic terms form a rank-1 QUBO matrix, enabling classical solvers to efficiently linearize the problem. For reaction netC. Quantum Annealers work pathway analysis, classical MIP/CP solvers D-Wave quantum annealers use qubit archi- (Gurobi, CP-SAT) vastly outperformed quantumtectures organized into Chimera, Pegasus, and inspired approaches, achieving optimality with Zephyr graphs [1, 4]. Quantum computers offer minimal computational effort. potential advantages for QUBO solving due to quantum phenomena like superposition and quan- D. Brute-force Approaches tum tunneling [4, 30]. Brute-force, by definition, explores all 2N posHowever, current quantum annealers have significant limitations: limited qubit count (D-Wave sible states, making its complexity inherently Advantage has ∼5000 qubits, while classical exponential. Despite this, brute-force remains imGPUs can handle problems with millions of portant. It guarantees optimality [1, 2] and serves variables) [1, 4, 21], limited connectivity topology as a benchmark for other methods [2, 33, 34]. requiring complex minor embedding [4], presence Also, it is often faster than heuristics for 30 of noise and other error sources [4], and uncertain variable problems [35] and acceptable for 50 time complexity for practically relevant problem variables [36], requires no parameter tuning [2], and provides the best scalability on GPU clussizes [1, 4]. Direct comparisons show that classical exact ters [34, 36], because traversing the state space methods often outperform modern quantum an- is an "embarrassingly parallel" task. The entire nealers [1, 4, 21]. In a study comparing the D- state space can be partitioned into independent Wave 2000Q quantum annealer with optimal clas- subsets processed by GPU cores simultaneously, sical Branch-and-Cut methods, classical methods without complex inter-node interaction [35, 37]. found optimal solutions faster for most prob- Additionally, brute-force enables not only finding the ground state but also constructing the lowlems [4]. A comprehensive comparison of D-Wave’s mod- energy spectrum [36], which is highly useful for ern hybrid quantum solver with classical meth- quantum algorithms [38]. Moreover, as noted in ods showed quantum advantages for large, dense [35], the accuracy of solutions obtained by heurisQUBO matrices, achieving relative accuracy of tic algorithms and quantum annealers depends on 0.013 and 6561x solution time reduction com- the type of problem, unlike that of brute-force. Brute-force in its naive form has O(N 2 · 2N ) pared to the best classical solver for 10,000complexity, as computing energy for each of 2N variable problems [21]. A recent study by Upadhyay and Jones [31] states takes O(N 2 ) time [35]. systematically compared quantum annealing (DHowever, optimized brute-force algorithms sigWave), digital annealing (Fujitsu DA v4 [32]), and nificantly improve this by using Gray code to classical MIP (Mixed Integer Programming)/CP traverse the state space, enabling incremental en(Constraint Programming) solvers on two indus- ergy updates in O(N ) time based on the previous trial problems. For mRNA codon optimization, value, reducing total complexity to O(N · 2N ) [2].
Nevertheless, algorithms with O(1) per state complexity are not novel. For such complexity, the evaluation rate in states/second is independent of problem size, making it a more practical metric than raw execution time. In the paper [35], constant-time energy updates are achieved by combining a tree-based depth-first traversal of the search space with Gray code ordering. This enables evaluation speeds reaching up to 0.07×1012 states per second on a single NVIDIA GeForce RTX 2080Ti GPU. Utilizing expensive 2×4 NVIDIA H100 GPUs, a brute-force solver can handle problems with up to 60 variables in approximately three days [1], which is equivalent to ∼ 0.56×1012 states per second on a single NVIDIA H100 GPU. However, as the corresponding work [39] is still in preparation, its algorithmic complexity and implementation details are not publicly known. Nevertheless, the near-linear trend of log Runtime [s] vs N in Fig. 1, implies an O(1) per-state complexity for the underlying algorithm.
Runtime [s]
106 105
1 day
104
1 hour
10
II. Mathematical Framework and Notation A. Prefix-Suffix Decomposition Approach Consider an arbitrary set of M binary state vectors arranged as rows in the matrix X ∈ {0, 1}M ×N . Given an upper-triangular QUBO matrix Q ∈ RN ×N , we compute the vector E ∈ RM where each element represents the energy of the corresponding state: E = diag(XQXT ).
3
102
and energy components precomputation reduce overall complexity to O(1) per state. Adjustable energy components decomposition enables aggressive caching of these components directly in GPU register memory. This minimizes memory bandwidth requirements and transforms the problem from memory-bound to compute-bound. The proposed CUDA implementation demonstrates an evaluation speed of 7.5 × 1012 states per second, setting a state-of-the-art performance standard for exhaustive solvers. Furthermore, this implementation achieves tera-scale performance on affordable, low-memory GPUs, making it a cost-effective solution.
1 minute
101
Brute Force AutoTune VeloxQ Custom VeloxQ
100 10−1
To establish the foundation for parallel computation of energy contributions, we introduce the following decomposition. Each binary state vector x ∈ {0, 1}N can be represented as a concatenation of two subvectors: x = [xp , xs ],
10−2
(1)
(2)
where xp ∈ {0, 1}A is the prefix part containing the first A bits, and xs ∈ {0, 1}B is the suffix part containing the last B bits, with A + B = N . Following the partitioning (2), we split the Figure 1. Brute-force runtime plot, reproduced from [1]. state matrix X ∈ {0, 1}M ×N and the QUBO maDeveloping faster and more accurate brute- trix Q ∈ RN ×N into conformal blocks according force solvers along with hybrid quantum-classical to prefix size A and suffix size B: algorithms represents a promising direction for Qpp Qps future research. X = [Xp , Xs ], Q = , 0 Qss E. Our Contribution where Xp ∈ {0, 1}M ×A , Xs ∈ {0, 1}M ×B , and In this work, we present a novel approach to A×A , Qps ∈ RA×B , and exact QUBO solving based on matrix formalism the submatrices Qpp ∈ R B×B . and prefix-suffix decomposition. It naturally maps Qss ∈ R the exhaustive search algorithm to the massively parallel GPU architecture. Gray code ordering 40
44
48 52 QUBO variables
56
60
The block multiplication yields: Qpp Qps XTp T = XQX = Xp Xs 0 Qss XTs = Xp Qpp XTp + Xp Qps XTs + Xs Qss XTs .
(3)
Taking the diagonal of both sides of equation (3), we obtain the energy decomposition given by equation (1): E = Ep + Eps + Es ,
(4)
where Ep = diag(Xp Qpp XTp ) is the prefix selfenergy, Eps = diag(Xp Qps XTs ) is the prefix-suffix interaction energy, and Es = diag(Xs Qss XTs ) is the suffix self-energy. For further use, we define the suffix-dependent energy vector as:
2A
A
B
z }| {
z }| { 2B
.. .
.. . 2B
.. .
.. . 2B
Figure 3. Decomposition of the state space UN into 2A chunks of 2B rows, and then into a prefix block and a suffix block.
xp xp .. = 1⊤ 2B · xp . xp
UB
Figure 4. Chunk structure
Esd = Es + Eps . The decomposition (4) enables independent Further, each chunk is partitioned into a prefix and parallel computation of energy contributions. block of the first A columns and a suffix block of the last B columns, as depicted in Figure 3. B. State Space Partitioning This method of partitioning UN exhibits sevN Let UN ∈ {0, 1}2 ×N denote the lexicograph- eral key properties: ically ordered set of all possible binary vectors 1) The suffix block is identical across all of size N . This set is represented as a matrix, chunks and is precisely the matrix UB . where each row corresponds to one binary vector, 2) The prefix block within any given chunk as illustrated in Figure 2. consists of identical rows. In other words, each chunk contains all possible suffixes N for a single, fixed prefix, as depicted in z }| { Figure 4. 3) For the i-th chunk, the common prefix is 0 0 ··· 0 0 0 0 0 · · · 0 0 1 the i-th row of the matrix UA . 0 0 · · · 0 1 0 4) A chunk can be uniquely identified by its N .. .. . . . . . 2 . . fixed prefix. We denote the chunk corre. .. .. .. 1 1 · · · 1 1 0 sponding to prefix xp as Cxp . 1 1 ··· 1 1 1 5) The prefix block of the chunk Cxp can be constructed as the product 1⊤ 2B · xp , where N 2B Figure 2. The complete state space matrix UN with 2 rows 12B ∈ R is a vector of ones of size 2B . and N columns, ordered lexicographically. Thus, the chunk Cxp can be expressed by 2A ×A We define the matrices UA ∈ {0, 1} and horizontal concatenation as: 2B ×B UB ∈ {0, 1} similarly for sizes A and B, reCxp = [1⊤ (5) 2B · xp , UB ]. spectively. Thus, UA contains all possible prefixes C. Energy Vector Partitioning and UB contains all possible suffixes. A The complete energy vector E(UN ) containing We partition the matrix UN into 2 nonoverlapping chunks, each containing 2B rows. energies of all 2N possible states can be con-
Table I Mathematical Notations and Definitions Notation
Entity
Description
N A B x xp xs X
Problem size Prefix size Suffix size State vector Prefix vector Suffix vector State sequence matrix
Xp
Prefix sequence matrix
Xs
Suffix sequence matrix
UN
Complete state space
UA
Prefix state space
UB
Suffix state space
Q Qpp Qps Qss E(x) E(X) Ep (xp ) Ep (Xp ) Es (xs ) Es (Xs ) Eps (xp , xs ) Eps (Xp , Xs )
QUBO matrix Prefix block Cross block Suffix block State energy States energies vector Prefix term Prefix energies vector Suffix term Suffix energies vector Interaction term Interaction energies vector
Esd (xp , xs ) Esd (Xp , Xs ) Cxp
Suffix dependent term Suffix dependent energies vector State space chunk
Total number of binary variables Number of bits in prefix, A + B = N Number of bits in suffix, A + B = N Complete binary state x = [xp , xs ] ∈ {0, 1}N First A bits of state, xp ∈ {0, 1}A Last B bits of state, xs ∈ {0, 1}B Ordered set of arbitrary M binary states {x1 , x2 , . . . , xM } represented as matrix X ∈ {0, 1}M ×N Ordered set of arbitrary M prefix vectors represented as matrix Xp ∈ {0, 1}M ×A Ordered set of arbitrary M suffix vectors represented as matrix Xs ∈ {0, 1}M ×B Lexicographically ordered set of all possible 2N binary states N represented as matrix UN ∈ {0, 1}2 ×N Lexicographically ordered set of all 2A prefix configurations A represented as matrix UA ∈ {0, 1}2 ×A Lexicographically ordered set of all 2B suffix configurations B represented as matrix UB ∈ {0, 1}2 ×B Upper triangular matrix, Q ∈ RN ×N QUBO submatrix, Qpp ∈ RA×A QUBO submatrix, Qps ∈ RA×B QUBO submatrix, Qss ∈ RB×B Scalar energy value xQxT Energies of states in X, E ∈ RM Scalar prefix energy contribution xp Qpp xTp Vector of prefix energy contributions for states in Xp , Ep ∈ RM Scalar energy contribution xs Qss xTs Vector of suffix energy contributions for states in Xs , Es ∈ RM Scalar energy contribution xp Qps xTs Vector of interaction energy contributions for state pairs, Eps ∈ RM Scalar Es (xs ) + Eps (xp , xs ) Vector Es (Xs ) + Eps (Xp , Xs ) Lexicographically ordered set of states with fixed prefix xp , B represented as matrix Cxp ∈ {0, 1}2 ×N
12B Mint
2B -size vector of ones Prefix-suffix interaction matrix
Matrix Mint = (Qps UTB ), size A × 2B
structed by concatenating the energy vectors of individual chunks: E(UN ) = [E(Cx(0) ), E(Cx(1) ), . . . , E(Cx(2A −1) )]. p p
We isolate the prefix-independent part of the second term and denote it as Mint . Mint = Qps UTB .
p
By substituting the chunk structure (5) into the energy decomposition formula (4), we obtain the energy vector for chunk Cxp : E(Cxp ) = Ep (xp )12B T + diag(1⊤ 2B xp Qps UB ) + Es (UB ),
(6)
where 12B is a 2B -dimensional vector of ones.
We can simplify the diagonal operation: diag(1⊤ 2B xp Mint ) = xp Mint .
(7)
Therefore, by substituting (7) into (6), the final expression for the chunk energy vector is: E(Cxp ) = Ep (xp )12B + xp Mint + Es (UB ).
(8)
Note that both Mint and Es (UB ) are independent The first one computes reusable energy comof the prefix configuration xp and can therefore ponents with the following steps: be precomputed. 1) Compute the all-prefix energy vector E(UA ) = diag(UA Qpp UTA ) D. Permutation Invariance 2) Compute the all-suffix energy vector The lexicographical ordering of state spaces UA E(UB ) = diag(UB Qss UTB ) and UB is not fundamental to the proposed de3) Construct the prefix-suffix interaction macomposition framework. The framework remains trix valid under arbitrary reordering of prefix and Mint = Qps UTB suffix configurations. The overall preprocessing complexity is O(A2 · Let πA and πB be arbitrary permutation op- 2A ) + O(B 2 · 2B ) + O(A · B · 2B ). erators acting on the row indices of UA and The main search loop systematically explores UB respectively, such that U′A = πA (UA ) and all 2A prefix configurations to identify the global U′B = πB (UB ). The permuted state spaces main- optimum through chunk-based processing. For tain completeness: each prefix configuration xp ∈ UA , the algorithm: ′ ′ ′ ′ ′ ′ ′ 1) Computes the suffix-dependent energy vecUN = {[xp , xs ] : xp ∈ UA , xs ∈ UB } = UN . tor Esd (xp ) = xp Mint + Es (UB ) Consequently, the complete energy vector 2) Finds the local minimum energy and correE(U′N ) is a permutation of E(UN ), preserving sponding state index within the chunk using the global minimum: Esd (xp ) ′ 3) Adds the prefix energy Ep (xp ), taken from min E(UN ) = min E(UN ). the precomputed E(UA ), to the local minThe interaction matrix transforms according to imum to obtain the final energy value the applied permutation of suffix states: 4) Updates the global minimum if the current local minimum is smaller M′int = Qps (U′B )T = Qps · πB (UB )T . The overall search loop complexity is O(A · Similarly, the suffix energy vector transforms 2B ) + O(2B ) per prefix configuration, resulting as: in total complexity of 2A [O(A · 2B ) + O(2B )] = Es (U′B ) = Es (πB (UB )). O(A·2N ), which dominates over the preprocessing complexity when A (or B) are not significantly Crucially, these transformed quantities retain smaller than N . Thus, the total algorithm comtheir prefix-independence and remain suitable plexity is O(A · 2N ). for precomputation. The chunk energy compuB. Gray-code optimized algorithm with O(2N ) tation (8) generalizes to: Complexity E′ (Cx′p ) = Ep (x′p )12B + x′p M′int + Es (U′B ). The key insight is to iterate through all 2A This property provides flexibility in state space prefixes in Gray code order [40], where each traversal ordering, which can affect computa- consecutive element differs by only a single bit. tional complexity and parallelization strategies. Decimal
0
1
2
3
4
5
6
7
III. Algorithm description Gray code 000 001 011 010 110 111 101 100 A. Basic brute-force algorithm with O(A · 2N ) Figure 5. 3-bit Gray code sequence Complexity (i) (i+1) Let xp and xp be two such prefixes, differing The complete brute-force QUBO solver algorithm operates by partitioning the search space at position k (note that k follows the same left-tointo chunks based on prefix configurations. The right order as the vector itself). The prefix change algorithm consists of two parts: preprocessing and is: ∆xp = x(i+1) − x(i) main search loop. p p = ±ek ,
Algorithm 1: QUBO Solver Algorithm Input : QUBO matrix Q ∈ RN ×N , partition sizes A, B where A+B =N Output: Optimal solution xopt and minimum energy Emin EA ← ComputePrefixEnergies(Qpp ); EB ← ComputeSuffixEnergies(Qss ); Mint ← Qps UTB ; Emin ← +∞; foreach xp ∈ UA do Ecross ← xp Mint ; Esd ← Ecross + EB ; idxlocal ← arg min(Esd ); Elocal ← min(Esd ) + EA [xp ]; if Elocal < Emin then Emin ← Elocal ; xopt ← ConstructState(xp , idxlocal ); end end return (xopt , Emin )
where ek is the k-th standard basis vector and the sign depends on flip direction. Given that Esd (xp ) = xp Mint + Es (UB ) from (8), the change in suffix-dependent energy vector is: (i+1)
∆Esd = Esd
(i)
− Esd = ∆xp Mint .
Therefore: (k)
∆Esd = ±ek Mint = ±Mint ,
(k)
where Mint denotes the k-th row vector of the matrix Mint . Hence, we can update the energy vector incrementally: (i+1) (i) (k) Esd = Esd ± Mint . (9) This reduces the update operation cost from A·2B to 2B additions/subtractions. Such improvement is particularly significant as it simplifies the dominant computational term in the algorithm’s overall complexity. (0) The algorithm starts with zero prefix xp = (0) 0 and computes initial energy vector Esd = Es (UB ). For each subsequent prefix in Gray code order:
1) Get the flipped bit position k and direction δ ∈ {−1, +1} 2) Update the energy vector using (9): (i+1) (i) (k) Esd = Esd + δMint 3) Find the local minimum and update the global optimum if needed For computing prefix energies E(UA ) and suffix energies E(UB ), the same Gray-code algorithm can be employed. However, instead of searching for minima within each chunk, the energy values are stored to construct the complete energy vectors. This approach reduces the asymptotic complexity of preprocessing from O(A2 · 2A ) + O(B 2 · 2B ) to O(2A ) + O(2B ). Since the main search loop complexity is O(2N ), the preprocessing overhead remains minimal in comparison. Thus, the overall complexity has been reduced from O(A · 2N ) to O(2N ). Algorithm 2: Gray-code QUBO Solver Algorithm Input : QUBO matrix Q ∈ RN ×N , partition sizes A, B where A+B =N Output: Optimal solution xopt and minimum energy Emin EA ← ComputePrefixEnergies(Qpp ); EB ← ComputeSuffixEnergies(Qss ); Mint ← Qps UTB ; Emin ← +∞; xp ← 0; Esd ← EB ; for i = 0 to 2A − 1 do idxlocal ← arg min(Esd ); Elocal ← min(Esd ) + EA [xp ]; if Elocal < Emin then Emin ← Elocal ; xopt ← ConstructState(xp , idxlocal ); end if i < 2A − 1 then (k, δ) ← ComputeGrayFlip(i); xp [k] ← xp [k] ⊕ 1; Esd ← Esd + δ · Mint [k, :]; end end return (xopt , Emin )
C. Parallelization Strategies We describe two primary parallelization approaches. Both strategies partition the entire search space into segments. This enables a mapreduce scheme for global minimum search. Each processing unit computes local minima for its assigned segment independently. The global minimum is found by reducing the local results. Furthermore, these two approaches can be combined for hybrid execution. 1) Straightforward Parallelization: Decompose the full search space of size 2N into 2k independent subproblems of size N − k. Split x = [xf , xr ] with xf ∈ {0, 1}k , xr ∈ {0, 1}N −k . Then partition Q into blocks Qf f , Qf r , Qrr as in II-A. For each fixed xf , form the subproblem by reducing original QUBO matrix as follows: Q′ = Qrr + diag(xf Qf r ). For this xf , the energy becomes E(x) = c(xf ) + E ′ (xr ), where c(xf ) = xf Qf f x⊤ f, E ′ (xr ) = xr Q′ x⊤ r .
A. Baseline Implementations The core logic is encapsulated within a QBF() solver class, which can be configured for either CPU or GPU execution through NumPy [42] and CuPy [43] backends respectively. Both configurations share the same codebase which dynamically selects the appropriate computational backend and memory-management calls at runtime. The implementation supports arbitrary data types while delegating overflow handling and parallelism management to the respective libraries. 1) NumPy-based CPU Baseline: The NumPy implementation serves as a straightforward, hardware-agnostic reference. By operating on whole vector rows, the algorithm achieves efficiency through cache-friendly memory access patterns and processor vector instructions. While currently operating on a single CPU core, the implementation can be readily adapted to multithreading through one of the parallelization techniques III-C. 2) Cupy-based GPU Baseline: The CuPybased implementation provides a GPU-native execution path while maintaining codebase compatibility with the NumPy variant. It leverages vectorized operations and implicit parallelization across CUDA cores. However, the inability to fuse addition and minimum kernels in CuPy leads to multiple memory transfers per iteration, creating a memory-bandwidth bottleneck. This limitation motivated the development of a custom CUDA kernel.
Each of the 2k subproblems is solved independently. The global minimum is obtained over the 2k local minima, each augmented by its corresponding constant c(xf ). 2) Column-wise Parallelization: This approach partitions the suffix space UB into multiple segments UB = [UB1 |UB2 | . . .]. Such partitioning induces a corresponding decomposition of the B. CUDA-based Optimized Implementation interaction matrix Mint = [M1 |M2 | . . .] and enTo overcome the limitations of the baseline ergy vector E(UB ) = [E(UB1 )|E(UB2 )| . . .]. Units implementation, we develop a custom CUDA [44] share the complete E(UA ) vector while storing kernel-based version of the algorithm. This impleonly their respective Mi and E(UBi ) portions. mentation targets peak brute-force performance. IV. Implementation Details Our optimization efforts focus on two distinct We develop two Python-based implementa- GPUs: NVIDIA V100 and NVIDIA H100. For tions: a baseline version for demonstration and clarity, this section focuses on the V100-specific validation, and a high-performance optimized ver- optimizations, while the maximum performance sion designed to achieve maximum performance. was achieved on the more modern H100. Many The baseline implementation is available as design choices are driven by the characteristics of open-source software through the qubo_lib li- the target GPU hardware. Below we describe the brary on GitHub [41]. key optimization approaches:
1) Register-Based Energy Components Storage: The primary optimization focuses on maximizing the utilization of GPU registers, the fastest available memory tier. Since we are able to freely choose A and B partition sizes, we can adjust them to make Mint and Esd small enough. Consequently, both components can reside completely in GPU register memory. Our first target GPU is the NVIDIA Tesla V100 with 80 SMs (streaming multiprocessors) and a 20 MB total register file (256 KB on each SM) [45]. This leads to the optimal suffix size B = 18. For B ≥ 19, insufficient register memory limits the prefix size, and leading to trivial problem sizes. Values of B ≤ 17 underutilized available resources. 2) Column-Wise Parallelization: We decompose Mint and Esd into segments, as shown in section III-C2, using CUDA threads as the computational units. CUDA launch-configuration dimensions are set empirically to 256 blocks and 64 threads per block, creating 16384 total threads. This configuration assigns 3.2 threads per FP32 core across the GPU’s 5120 (64 per SM) cores, which presumably improves warpscheduling efficiency. Column-wise parallelization implies that all threads share the prefix-energy vector. This vector is substantial in size (4 GB for N = 49), requiring storage in GPU global memory. However, at a brute-force rate of 1.1 × 1012 states per second with B = 18, only 4.2 × 106 prefix energy fetches per second are required. This results in a insignificant memory bandwidth of 8 MB/s, confirming the compute-bound nature of our implementation. 3) Coalesced Memory Access Patterns: Natural lexicographical prefix energy computation conflicts with Gray code traversal. This causes scattered global memory accesses, which cannot be fully mitigated by caching. Reordering the energy vector to match Gray code sequence achieves fully coalesced access patterns. This reordering yields substantial performance gains, making the preprocessing overhead worthwhile. 4) Data Types: Storing energy components in GPU register memory directly influences data type selection. For energy computations, we
choose the 16-bit signed integer (int16) data type. The integer type for energy eliminates floating-point accumulation errors. The 16-bit sizeof halves memory usage and doubles performance compared to 32-bit types. For prefix representation we choose signed int32 type, imposing a practical limit of A ≤ 31. This allows solving problems up to size 49 (A = 31, B = 18), which is still large enough to be computationally challenging. Unsigned uint32 type is rejected due to slower computation, 64-bit types are unnecessary. 5) Overflow Check Elimination: Overflow control is critical for integer types but cannot be performed within the CUDA kernel without performance loss. Instead, the kernel assumes external guarantees against overflow. Prior to launch, the QUBO problem is solved using a simulated annealing solver. The obtained minimum energy provides a lower bound for intermediate energy values during brute-force. The upper bound is estimated by solving the negated QUBO matrix. If the resulting energy spectrum fits within int16 range, overflow is deemed unlikely. This method requires safety margins due to its approximate nature. 6) Branch Elimination and Min Reduction: Critical loop optimization focuses on eliminating branching within the innermost computation. Arithmetic operations were restructured to avoid conditional addition/subtraction selection. Minimum finding was implemented using tree-based reductions with intrinsic min() operations rather than iterative searches. This provides measurable performance gains through improved instruction throughput. CUDA architecture complexity made analytical parameter selection impractical. Final configuration emerged through extensive empirical experimentation, balancing register usage, memory bandwidth, and computational throughput. The non-linear optimization process reflected the complex interactions between hardware constraints and algorithm requirements. C. Scaling Challenges As discussed in IV-B4, the CUDA kernel is constrained to problem sizes up to N = 49.
To overcome this limit, we adopt the embarrassingly parallel decomposition introduced in Section III-C1. Rather than modifying the highly optimized CUDA kernel, the problem structure is adapted to fit the kernel constraint. We set |xr | to 49, and the remaining |xf | = N − 49 bits induce a subproblem matrix Q′ as defined in III-C1. For each subproblem, the CUDA kernel finds the local minimum. The global minimum is then reduced over all local minima. Since each node independently constructs E(UA ), Mint , and E(UB ), the approach carries a small computational overhead. We estimate this overhead at under 0.3% of total runtime (Table IV). V. Experimental Evaluation A. Setting All experiments were conducted on a server running Debian 12 (bookworm) with the Linux kernel 5.15.0-130-generic and Python 3.10.14. The experiments were executed within a Docker container running on a virtualized server to ensure a reproducible software environment. The computational hardware consists of a dualsocket Intel Xeon Gold 6348 CPU [46], providing 28 physical cores at 2.60 GHz. The CPU features a multi-level cache hierarchy with 48 KB L1 data cache and 32 KB L1 instruction cache per core, a 1.25 MB L2 cache per core, and a 39 MB shared L3 cache per socket. The processor supports a wide range of SIMD instruction sets, including AVX2, FMA3, and AVX-512 [47], which are leveraged by the numerical libraries. The system is equipped with 499 GB of RAM. The primary computational device is an NVIDIA Tesla V100S-PCIE-32GB GPU, which is based on the Volta architecture [45]. It features 80 Streaming Multiprocessors (SMs) for a total of 5120 CUDA cores. Each SM has access to a 128 KB L1 cache and a shared 6 MB L2 cache. The GPU operates at a boost clock of 1597 MHz, delivers a theoretical 16.4 TFLOPS of singleprecision (FP32) performance, and provides a memory bandwidth of 1130 GB/s. Additionally, an NVIDIA H100 80GB GPU [48] based on the Hopper architecture was employed to assess peak performance. This accelerator
integrates 132 Streaming Multiprocessors (SMs), amounting to 16896 CUDA cores, with a 256 KB L1 cache per SM and a 50 MB shared L2 cache. Running at a boost clock of 1845 MHz, it offers a theoretical 66.9 TFLOPS of FP32 compute and a memory bandwidth of 3350 GB/s. The software stack includes NVIDIA driver version 565.57.01, CUDA Toolkit 12.7, and CUDA Runtime version 12.9. Our implementations leverage NumPy 2.2.6 for the CPU baseline and CuPy 13.6.0 for the GPU baseline and as the interface for our custom CUDA kernel, which was compiled using GCC 12.2.0. The NumPy-based CPU baseline utilizes the scipy-openblas library (version 0.3.29) for its linear algebra operations, which was compiled with support for SIMD extensions up to AVX-512. B. GPU Profiling We profiled the CuPy baseline and the custom CUDA kernel to identify performance bottlenecks with Nsight Compute [49]. The profiling was performed on the V100 GPU only, as we were unable to profile the H100. The CuPy implementation is fundamentally memory-bound. Its inefficiency stems from a poor memory access pattern, specifically repeated global memory accesses. For each prefix, the kernel reads the energy vector Esd and a row of the matrix Mint , writes the updated vector back to global memory, and then immediately reads it again for the minimum reduction. This costly memory traffic could be mitigated by an effective cache, but the Read-Modify-Write memory access pattern [50] induces severe cache pollution, rendering the cache useless. The profiler confirms this diagnosis, showing a cache hit rate below 1%. This explains the severe memory bandwidth saturation we observe, which makes the implementation profoundly memorybound. For the int16 data type, the main search loop transfers 8 bytes of data per state, totaling 8 TB of memory traffic for a problem with N = 40. This establishes a theoretical lower bound of 8 seconds on the runtime, assuming an optimistic memory bandwidth of 1 TB/s. The profiler reports an occupancy of 60-80% for the main kernels. This is combined with CPU
synchronization and the overhead of launching six separate kernels per main loop iteration. These factors explain the actual execution time of 15.7 seconds for N = 40. In contrast, the custom CUDA kernel is compute-bound. Its primary advantage is the elimination of global memory traffic during the main search loop. Profiling confirms this with a negligible memory bandwidth consumption of 87.25 MB/s and a near-perfect cache hit rate. The kernel achieves 75% of the theoretical compute throughput. The primary reason for not achieving the theoretical peak is the compiler’s data placement strategy, which moved data from registers to shared memory, making the number of registers per warp the limiting factor. As a result, the custom CUDA kernel solves the problem for N = 40 in just 0.58 seconds, versus the 15.7 seconds taken by the CuPy baseline.
a direct comparison, we linearly normalized these results to a single GPU. This optimistic scaling assumes perfect parallel performance with zero communication overhead. Second, we group results from the NVIDIA RTX 2080Ti (4352 CUDA cores, 13.5 TFLOPS) and Tesla V100 GPUs (5120 CUDA cores, 16.4 TFLOPS), justified by their relatively small performance gap (15-20%). Third, our CPU benchmark for QBF benefits from the algorithm’s inherent vectorizability, enabling the NumPy backend to explicitly use of AVX-512 instructions for processing 32 integers (16-bit) in parallel. In contrast, the CPU performance of QUBO-ES is reported without details on its use of SIMD instructions, creating an uneven comparison that likely favors our implementation. Despite approximations introducing a margin of error, these minor inaccuracies do not alter the C. Performance Comparison qualitative conclusions of the comparison. We evaluated our QBF solver against leading The high performance of the QBF solver stems brute-force algorithms [1, 35, 39] across three directly from its underlying algorithmic design. hardware tiers: single CPU core, mid-range GPU The method is inherently vectorizable for CPUs (RTX 2080Ti/V100), and high-end GPU (H100). and maps naturally to the massive parallelism We assess both our baseline implementation and of GPUs, resulting in a compute-bound, cachethe optimized CUDA version. For the benchmark friendly implementation. This synergy between we use a uniformly random matrix, since the the algorithm and modern hardware architecture exhaustive search algorithm’s execution trace is pushes the practical boundaries of the exhaustive independent of the specific problem instance. search method. Since all compared solvers achieve O(1) comVI. Discussion plexity per state, their runtime scales linearly with the size of the search space, 2N . This prop- A. Why Not Tensor Cores erty makes the state evaluation rate (states/sec) A natural question is whether the algorithm a primary performance metric instead of raw can be accelerated using Tensor Cores, which on execution times. the V100 GPU offer a theoretical peak of 112.2 The performance results are summarized in TFLOPS compared to 16.4 TFLOPS for standard Table II. For a more intuitive understanding of CUDA cores. The answer is likely no. The high throughput of Tensor Cores is the scale, Table III presents the same data as log2 (states/sec). This metric directly corresponds achieved through the HMMA instruction, a 4x4 to the QUBO problem size solvable in one second. matrix multiply-add operation (D = A × B + C) Finally, Table IV provides the raw execution counted as 128 FLOPs. Our core computational kernel, however, consists of a vector addition foltimes for reference. (k) Our comparison involves several methodolog- lowed by a minimum reduction: Esd = Esd +Mint ical approximations. First, as the work [39] is and min(Esd ). unpublished, we rely on approximate values exWhile the addition could be mapped to the tracted from their performance plot (Fig. 1). The addition part of HMMA by reshaping vectors into results were obtained on a 2x4 H100 rig; to enable 4x4 matrices, this would only utilize 16 of the 128
Table II State evaluation rate, 1012 states/sec QUBO-ESa 0.000053 0.066 —
CPU Mid-range GPU: 2080Ti, V100 High-end GPU: H100
QBFb 0.0051 0.072 —
QBF CUDAc — 2.33 7.47
Bruteforced — — ∼ 0.56
QBFb 32.3 36.0 —
QBF CUDAc — 41.1 42.8
Bruteforced — — ∼ 39.0
a
QUBO-ES bruteforce (tree DFS + Gray code) from [35] this work, baseline c this work, optimized d Bruteforce implementation from [1, 39], approximated and scaled b
Table III One-second QUBO size: log2 (states/sec) QUBO-ESa 25.7 35.9 —
CPU Mid-range GPU: 2080Ti, V100 High-end GPU: H100 a
QUBO-ES bruteforce (tree DFS + Gray code) from [35] this work, baseline c this work, optimized d Bruteforce implementation from [1, 39], approximated and scaled b
Table IV Raw execution time (sec) N
CPU Mid-range GPU: 2080Ti, V100 QUBO-ESa QBF(NumPy)b QUBO-ESa QBF(CuPy)b QBF CUDAc 34 322.7 3.135* (0.053) 0.269 0.571 (0.310) 0.120 (0.111) 35 641.9 6.950 (0.055) 0.518 0.775 (0.298) 0.150 (0.131) 36 1289 13.28 (0.057) 1.021 1.259 (0.297) 0.167 (0.130) 37 2565 26.93 (0.060) 2.102 2.255 (0.315) 0.205 (0.139) 38 5153 54.08 (0.063) 4.082 4.176 (0.313) 0.258 (0.141) 39 10265 107.5 (0.079) 8.513 7.933 (0.298) 0.431 (0.195) 40 20629 211.9 (0.086) 16.59 15.66 (0.298) 0.575 (0.112) 41 — 434.3 (0.074) 32.93 31.08 (0.311) 1.049 (0.112) 42 — 872.3 (0.082) 65.85 61.71 (0.299) 1.970 (0.115) 43 — 1686 (0.080) 134.3 123.5 (0.299) 3.841 (0.117) 44 — 3460 (0.081) 272.6 246.2 (0.299) 7.59 (0.120) 45 — 6752 (0.084) 525.4 491.7 (0.302) 15.22 (0.127) 46 — 13690 (0.119) 1053 984.2 (0.327) 30.49 (0.163) 47 — — 2097 1967 (0.300) 61.49 (0.524) 48 — — 4210 3930 (0.309) 121.3 (0.246) 49 — — 8515 7879 (0.335) 238.3 (0.487) 50 — — 17038 15708 (0.296) 485.7 (2.232) 51 — — — — 968.7 (5.256) 52 — — — — 1945 (10.58) 53 — — — — 3877 (21.00) 54 — — — — 7701 (15.18) 55 — — — — 15464 (29.89) 56 — — — — — * Time in parentheses denotes precomputation and synchronization overhead. a QUBO-ES bruteforce (tree DFS + Gray code) from [35] b this work, baseline c this work, optimized d Bruteforce implementation from [1, 39], approximated and scaled
High-end GPU: H100 QBF CUDAc Bruteforced 0.041 (0.037) — 0.049 (0.043) — 0.051 (0.040) — 0.060 (0.039) — 0.079 (0.042) — 0.110 (0.036) — 0.182 (0.038) ∼ 15.7 0.335 (0.042) — 0.620 (0.041) ∼ 22.7 1.200 (0.037) — 2.372 (0.038) ∼ 25.0 4.756 (0.041) — 9.501 (0.036) ∼ 170 19.01 (0.052) — 37.85 (0.057) ∼ 540 74.36 (0.077) — 151.2 (0.160) ∼ 2040 301.3 (0.328) — 604.7 (0.592) ∼ 8030 1207 (2.558) — 2403 (2.623) ∼ 32000 4826 (5.095) — — ∼ 128000
FLOPs. This reduces the effective throughput to approximately (16/128) × 112.2 ≈ 14.0 TFLOPS. Furthermore, the subsequent minimum reduction cannot be performed on Tensor Cores and would require a separate pass on standard CUDA cores, introducing additional overhead. Therefore, the current CUDA core-based implementation remains the most efficient and suitable architecture for this problem.
with Gray code traversal. This enables registerbased caching, reduces memory bandwidth requirements to negligible levels, and shifts the bottleneck from memory-bound to compute-bound. Our custom CUDA implementation reaches a state evaluation rate of 7.5 × 1012 states per second on a single H100 GPU for the int16 data type. An even more significant result is 2.3×1012 states per second on a low-memory V100 GPU, which is more affordable and accessible. B. Future Work This makes the proposed solver a practical tool The proposed approach opens up several for researchers. The results outperform the best promising directions for future enhancement. known alternatives by more than an order of magFirst, the current implementation is optimized nitude, setting a new state-of-the-art for brutefor dense matrices. Adapting the framework for force exact QUBO solvers. The baseline source sparse matrices represents another promising av- code is publicly available [41]; the optimized enue. This would involve reformulating the matrix solver is also provided for experimental access [? operations and data structures to exploit matrix ]. sparsity could yield performance gains for some References problems. The second future direction is supporting [1] J. Pawlowski, J. Tuziemski, P. Tarasiuk, floating-point QUBO matrices. This requires adA. Przybysz, R. Adamski, K. Hendzel, dressing the challenge of numerical error accumuL. Pawela, and B. Gardas. VeloxQ: A fast lation inherent in the long sequential incremental and efficient QUBO solver. arXiv preprint, Gray-code updates, potentially through mixed2025. arXiv:2501.19221. precision or periodic correction schemes. [2] Sascha Mücke. Faster QUBO brute-force Third direction involves adapting the solver for solving using gray code. arXiv preprint, newer accelerator architectures, such as Ampere 2023. arXiv:2310.19373. or Hopper. However, achieving high performance [3] G. Kochenberger, J.-K. Hao, F. Glover, will require a new round of empirical tuning to M. Lewis, and Z. Lü. The unconstrained re-balance key parameters for the new hardware. binary quadratic programming problem: A This effort is necessary to unlock the full potential survey. Journal of Combinatorial Optimizaof any future accelerator. tion, 28:58–81, 2014. Furthermore, the framework can be extended [4] Michael Jünger, Elisabeth Lobe, Petra Mutzel, Gerhard Reinelt, Franz Rendl, beyond finding the single ground state to generate Giovanni Rinaldi, and Tobias Stollenwerk. the low-energy spectrum. This can be achieved Quantum annealing versus digital computby modifying the kernel to maintain a list of ing: An experimental comparison. ACM J. the k lowest energies. In turn, this provides a Exp. Algorithmics, 26, July 2021. more comprehensive benchmark for quantum and heuristic solvers by revealing the structure of the [5] Fred Glover, Gary Kochenberger, Ryan Hennig, and Yuan Du. Quantum bridge analytics solution landscape. I: a tutorial on formulating and using QUBO VII. Conclusion models. Annals of Operations Research, This paper presents a novel parallel exhaus314:141–183, 2022. tive search algorithm for QUBO problems with [6] Koji Nakano, Shunsuke Tsukiyama, Yasuaki dense matrices, achieving O(1) complexity per Ito, Takashi Yazane, Junko Yano, Takumi state. The approach combines complete state Kato, Shiro Ozaki, Rie Mori, and Ryota space energy vector prefix–suffix decomposition Katsuki. Dual-matrix domain wall: A novel
technique for generating permutations by qubo and ising models with quadratic sizes. Technologies, 11(5), 2023. [7] A. Lucas. Ising formulations of many NP problems. Frontiers in Physics, 2:5, 2014. [8] Ikuto Nakatsukasa, Koji Nakano, Victor Parque, and Yasuaki Ito. Optimizing heat treatment schedules via qubo formulation. Applied Sciences, 15(16), 2025. [9] Jin-Guo Liu, Lei Wang, and Pan Zhang. Tropical tensor network for ground states of spin glasses. Phys. Rev. Lett., 126:090506, Mar 2021. [10] Jinglei Cheng, Ruilin Zhou, Yuhang Gan, Chen Qian, and Junyu Liu. Scalable community detection using quantum hamiltonian descent and qubo formulation. In 2025 62nd ACM/IEEE Design Automation Conference (DAC), pages 1–7, 2025. [11] H. Goto, K. Tatsumura, and A. R. Dixon. Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems. Science Advances, 5(4):eaav2372, 2019. [12] H. Goto, K. Endo, M. Suzuki, Y. Sakai, T. Kanao, Y. Hamakawa, R. Hidaka, M. Yamasaki, and K. Tatsumura. Highperformance combinatorial optimization based on classical mechanics. Science Advances, 7(6):eabe7953, 2021. [13] T. Okuyama, T. Sonobe, K. Kawarabayashi, and M. Yamaoka. Binary optimization via momentum annealing. Physical Review E, 100:012111, 2019. [14] F. Barahona, M. Jünger, and G. Reinelt. Experiments in quadratic 0–1 programming. Mathematical Programming, 44(1):127–137, 1989. [15] F. Liers, M. Jünger, G. Reinelt, and G. Rinaldi. Computing exact ground states of hard ising spin glass problems by branchand-cut. In New Optimization Algorithms in Physics. Wiley-VCH, 2004. [16] Franz Rendl, Giovanni Rinaldi, and Angelika Wiegele. Solving Max-Cut to optimality by intersecting semidefinite and polyhedral relaxations. Mathematical Programming, 121(2):307–335, 2010.
[17] Daniel Rehfeldt, Thorsten Koch, and Yuji Shinano. Faster exact solution of sparse MaxCut and QUBO problems. Mathematical Programming Computation, 15(3):445– 470, 2023. [18] Timotej Hrga and Janez Povh. Madam: a parallel exact solver for max-cut based on semidefinite programming and admm. Computational Optimization and Applications, 80(2):347–375, November 2021. [19] David L. Applegate, Robert E. Bixby, Vašek Chvátal, and William J. Cook. 2007. [20] David L. Applegate, Robert E. Bixby, Vašek Chvátal, William Cook, Daniel G. Espinoza, Marcos Goycoolea, and Keld Helsgaun. Certification of an optimal tsp tour through 85,900 cities. Operations Research Letters, 37(1):11–15, 2009. [21] Seongmin Kim, Sang-Woo Ahn, In-Saeng Suh, Alexander W. Dowling, Eungkyu Lee, and Tengfei Luo. Quantum annealing for combinatorial optimization: a benchmarking study. npj Quantum Information, 11(177):1– 12, 2025. [22] Tingting Zhang, Qichao Tao, Bailiang Liu, and Jie Han. A review of simulation algorithms of classical ising machines for combinatorial optimization. In 2022 IEEE International Symposium on Circuits and Systems (ISCAS), pages 1877–1881, 2022. [23] N. Mohseni, P. L. McMahon, and T. Byrnes. Ising machines as hardware solvers of combinatorial optimization problems. Nature Reviews Physics, 4:363–379, 2022. [24] M. Jiang, K. Shan, C. He, and C. Li. Efficient combinatorial optimization by quantuminspired parallel annealing in analogue memristor crossbar. Nature Communications, 14:5927, 2023. [25] Daisuke OKU, Kotaro TERADA, Masato HAYASHI, Masanao YAMAOKA, Shu TANAKA, and Nozomu TOGAWA. A fullyconnected ising model embedding method and its evaluation for cmos annealing machines. IEICE Transactions on Information and Systems, E102.D(9):1696–1706, 2019. [26] Kasho Yamamoto, Kazushi Kawamura, Kota Ando, Normann Mertig, Takashi Take-
moto, Masanao Yamaoka, Hiroshi Teramoto, Akira Sakai, Shinya Takamaeda-Yamazaki, and Masato Motomura. Statica: A 512spin 0.25m-weight annealing processor with an all-spin-updates-at-once architecture for combinatorial optimization with complete spin–spin interactions. IEEE Journal of Solid-State Circuits, 56(1):165–178, 2021. [27] Toshimori Honjo, Tomohiro Sonobe, Kensuke Inaba, Takahiro Inagaki, Takuya Ikuta, Yasuhiro Yamada, Takushi Kazama, Koji Enbutsu, Takeshi Umeki, Ryoichi Kasahara, Ken ichi Kawarabayashi, and Hiroki Takesue. 100,000-spin coherent ising machine. Science Advances, 7(40):eabh0952, 2021. [28] Martin J. A. Schuetz, J. Kyle Brubaker, and Helmut G. Katzgraber. Combinatorial optimization with physics-inspired graph neural networks. Nature Machine Intelligence, 4(4):367–377, 04 2022. [29] Changjun Fan, Mutian Shen, Zohar Nussinov, Zhong Liu, Yizhou Sun, and Yang-Yu Liu. Searching for spin glass ground states through deep reinforcement learning. Nature Communications, 14(1):725, 02 2023. [30] Michael A. Nielsen and Isaac L. Chuang. Quantum Computation and Quantum Information: 10th Anniversary Edition. Cambridge University Press, 2010. [31] M. Upadhyay and M. N. Jones. Comparative studies of quantum annealing, digital annealing, and classical solvers for reaction network pathway analysis and mRNA codon selection. arXiv preprint, 2025. arXiv:2509.09862. [32] Fujitsu Limited. Fujitsu digital annealer user’s guide. API Documentation, 2025. https://portal.aispf.global. fujitsu.com/apidoc/da/jp/da-guide-en.html. [33] Konrad Jałowiecki and Łukasz Pawela. Omnisolver: An extensible interface to ising spin–glass and qubo solvers. SoftwareX, 24:101559, 2023. [34] V. Mehta, F. Jin, K. Michielsen, and H. De Raedt. On the hardness of quadratic unconstrained binary optimization problems. Frontiers in Physics, 10:956882, 08 2022. [35] Masaki Tao, Koji Nakano, Yasuaki Ito, Ryota Yasudo, Masaru Tatekawa, Ryota Kat-
suki, Takashi Yazane, and Yoko Inaba. A work-time optimal parallel exhaustive search algorithm for the qubo and the ising model, with gpu implementation. In 2020 IEEE International Parallel and Distributed Processing Symposium Workshops (IPDPSW), pages 557–566, 2020. [36] Konrad Jałowiecki, Marek M. Rams, and Bartłomiej Gardas. Brute-forcing spin-glass problems with cuda. Computer Physics Communications, 260:107728, 2021. [37] Konrad Jałowiecki, Marek M. Rams, and Bartłomiej Gardas. Brute-forcing spin-glass problems with cuda. Computer Physics Communications, 260:107728, 2021. [38] A.A. Zhukov, A.V. Lebedev, and W.V. Pogosov. Grover’s search meets ising models: A quantum algorithm for finding low-energy states. Computer Physics Communications, 313:109627, 2025. [39] K. Jałowiecki, L. Pawela, B. Gardas, A. Przybysz, and J. Tuziemski. Gpu based brute-force solver for QUBO and ising instances. In preparation, cited in arXiv:2501.19221, 2025. [40] Torsten Mütze. Combinatorial Gray codes—an updated survey. arXiv preprint, 2024. arXiv:2202.01280 [math.CO]. [41] Cloud.ru. open source github repository. github.com, 2026. https://github.com/ cloud-ru-tech/evo-quantum-qubolib. [42] NumPy Developers. Numpy: the fundamental package for scientific computing with python. numpy.org, 2025. Version 2.2.3. [43] CuPy Developers. Cupy: Numpy & scipy for gpu. cupy.dev, 2025. Version 13.4.1. [44] NVIDIA Corporation. Cuda toolkit documentation. Documentation, 2025. https: //docs.nvidia.com/cuda/. [45] NVIDIA Corporation. Nvidia tesla v100 gpu accelerator. Datasheet, 2017. https://images.nvidia.com/ content/technologies/volta/pdf/ tesla-volta-v100-datasheet-letter-fnl-web. pdf. [46] Intel Corporation. Intel xeon gold 6348 processor specifications. Product Specifications, 2021. https://www.intel.com/
content/www/us/en/products/sku/212456/ intel-xeon-gold-6348-processor-42m-cache-2-60-ghz/ specifications.html. [47] Intel Corporation. Intel avx-512 instruction set. Architecture Instruction Set Reference, 2025. https://www.intel.com/content/ www/us/en/developer/articles/technical/ intel-avx-512.html. [48] NVIDIA Corporation. Nvidia h100 tensor core gpu. Datasheet, 2023. https://www.nvidia.com/content/dam/ en-zz/Solutions/gtcs22/data-center/h100/ PB-11133-001_v01.pdf. [49] NVIDIA Corporation. Nvidia nsight systems profiler. Documentation, 2025. https://docs. nvidia.com/nsight-systems/. [50] D. B. Kirk and W. W. Hwu. Programming massively parallel processors: A hands-on approach. 2010. Chapter 6: Memory Issues in CUDA.