Brute-forcing spin-glass problems with CUDA

TL;DR

This paper presents an exact, GPU-accelerated brute-force algorithm for computing low energy spectra of small spin-glass instances, serving as a benchmark for quantum and classical optimization methods.

Contribution

It introduces a fast, open-source CUDA-based implementation for exact energy spectrum calculation of spin-glass models on arbitrary graphs, enabling benchmarking and comparison with other algorithms.

Findings

Achieves unprecedented speed on commodity hardware.

Can handle spin-glass instances up to N=50.

Demonstrates the effectiveness of the approach as a benchmark tool.

Abstract

We demonstrate how to compute the low energy spectrum for small ($N\le 50$), but otherwise arbitrary, spin-glass instances using modern Graphics Processing Units or similar heterogeneous architecture. Our algorithm performs an exhaustive (i.e., brute-force) search of all possible configurations to select $S\ll 2^N$ lowest ones together with their corresponding energies. We mainly focus on the Ising model defined on an arbitrary graph. An open-source implementation based on CUDA Fortran and a suitable Python wrapper are provided. As opposed to heuristic approaches, ours is exact and thus can serve as a references point to benchmark other algorithms and hardware, including quantum and digital annealers. Our implementation offers unprecedented speed and efficiency already visible on commodity hardware. At the same time, it can be easily launched on professional, high-end graphics cards virtually at no extra effort. As a practical application, we employ it to demonstrate that the recent Matrix Product State based algorithm-despite its one-dimensional nature-can still accurately approximate the low energy spectrum of fully connected graphs of size $N$ approaching $50$.

Figures

Equations

$$H(\mathbf{s})=-\sum _{⟨i,j⟩}{J}_{ij}{s}_i{s}_j-\sum _{i=1}^N{h}_i{s}_i,$$

$$F(\mathbf{q})=-\sum _{⟨i,j⟩}{a}_{ij}{q}_i{q}_j-\sum _{i=1}^N{b}_i{q}_i,$$

$$a_{ij}=4J_{ij},b_i=2h_i-2\sum _{⟨i,j⟩}{J}_{ij}.$$

$$\text{GPU thread index}=(\mathbf{q}{)}_2,$$

$$\text{thread index}\#13=(00001101{)}_2,$$

$$\text{GPU thread index}+\text{offset}=(\mathbf{q}{)}_2,$$

$$\mathcal{O}[2^N(N^2+1)+S\log (S)]=\mathcal{O}[2^N(N^2+1)].$$

$$e^{-\beta H(\mathbf{s})/2}\approx A^{s_1}{A}^{s_2}\dots A^{s_L}=|\Psi (\beta )⟩,$$

$$\begin{array}{r}\text{includegraphics}\end{array}.$$

$$U_i(d\beta )=e^{-d\beta s_i({\sum }_{j>i}{J}_{ij}{s}_j+h_i)/2},$$

$$\begin{array}{r}\text{includegraphics}\end{array}.$$

Full text

1 Introduction

With increasing complexity and interconnectivity in the modern world, the ability to solve optimization problems becomes indispensable. Notwithstanding, these problems are fundamentally hard to resolve as they often require seeking over enormous spaces of possible solutions [1]. A notable example is the famous spin–glass problem encoded via the Ising model [2], where the low energy spectrum (the ground state in particular) is sought after. The importance of this system is reflected in the fact that many NP-complete [3] optimization problems (i.e. Karp’s $21$ problems [4]) can be mapped onto its Hamiltonian [5]. Furthermore, there is growing hardware support for many spin–glass based models [6, 7, 8]. These cutting edge technologies, when combined with classical neural networks [9], lead to quantum artificial intelligence [10]. A type of artificial intelligence believed to be powerful enough to simulate many–body quantum systems efficiently, which is a holy grail of modern physics [11].

The most promising ideas to overcome mathematical difficulties concerning classical optimization could rely on quantum computers [12]. In particular, on quantum annealers such as the D-Wave $2000$Q chip [13]. In principle, such machines could solve variate of (hard) optimization problems (almost) “naturally” by finding low energy eigenstates [14]. However, current quantum annealers are extremely noisy and thus not powerful enough to tackle large scale optimization challenges [15, 16]. In contrast, heuristic approaches, often offering superior performance, can not typically certify that the solution that has been found is, in fact, optimal [17, 18]. Most heuristic solvers rely on strategies ranging from famous simulated annealing [19], branch and bound approaches [20] their chordal extension [18], various Monte Carlo methods [21] throughout dynamical systems simulations [22] to tensor network analysis [23].

In this work, we focus on yet another class of solvers, namely those that perform exact brute–force search [24]. The idea is to search the entire Hilbert space exhaustively to find configurations with the lowest energies. Such a search can be performed either in the probability or energy space [23]. For all classical Hamiltonians, where all their terms commute, this is essentially equivalent to the exact diagonalization. However, in contrast to the quantum case, the eigenvalue problem for classical models can be executed truly in parallel. An efficient implementation nonetheless is not trivial.

Although practical applications of such solvers are limited to small problem sizes ($N\leq 50$), they can solve the Ising model that is defined on an arbitrary graph. Moreover, with the exhaustive search, one can easily certify the output. All of these features are crucial for testing, benchmarking, and validating new methods [25], strategies, and paradigms (e.g., memcomputing [22]) for solving classical optimization problems [26]. It is worth mentioning that brute–force approaches however limited can still serve as a reference point for today’s quantum supremacy experiments [27, 28].

Our implementation offers excellent flexibility and portability, as well as the significant efficiency and speed. Our solver can be executed on either CPU (Central Processing Unit) or GPU (Graphics Processing Unit) using Nvidia’s CUDA (Compute Unified Device Architecture). The latter architecture is of particular importance due to its massive parallel capabilities [29, 30]. Moreover, we provide a simple Python wrapper that allows users to access both architectures effortlessly [31].

Finally, we employ our solver to test the applicability of a particular tensor network ansatz—based on the Matrix Product State (MPS)—to optimization purposes of a fully connected graph of growing size. For a detailed description of this algorithm, we refer the reader to look at Supplementary Information in Ref. [23]. We have verified that indeed, such an ansatz, despite its inherited one–dimensional structure, can still successfully capture the low energy spectrum for tested graphs up to $N=50$. This indicates that the MPS ansatz should still perform well also for much larger systems having a dominant quasi–one–dimensional nature. At the same time, sparse connections at long–range do not necessary exclude the applicability of the MPS approach.

2 Spin-glass problems

In this work, we mainly focus on the Ising Hamiltonian. However, our approach can easily be extended to include other classical spin–glass models [32, 33]. To begin with, consider a simple undirected graph with $N$ nodes (i.e., vertices) as the one drawn in Fig. 1. We assign a unique spin variable, $s_{i}\pm 1$ (blue and red arrows), to each node. Adjacent nodes labeled as $i$, $j$ are coupled via interaction strength $J_{ij}$, which may be viewed as a weight of the edge connecting those two nodes. Additionally, for every spin, we associate a local magnetic field (bias) $h_{i}$ interacting with it. Then the energy of such a system of spins is defined as

[TABLE]

where $\boldsymbol{s}:=(s_{1},\ldots,s_{L})$. The first sum runs over all adjacent sites, which we denote here as $\langle i,j\rangle$.

In many practical applications, one is typically interested in finding a particular spin configuration, say $\boldsymbol{s}_{0}$, for which $H(\boldsymbol{s}_{0})$ in Eq. (1) admits its minimum value. Such configuration is called the ground state. Naturally, states with energies above the ground state energy are called excited states. Finding the low energy spectrum (consisting of the ground state energy and a number $S\ll 2^{N}$ of excited states) of the Ising model (1) can also be formulated as a Quadratic Unconstrained Optimization Problem (QUBO). Namely,

[TABLE]

where $\boldsymbol{q}=(\boldsymbol{s}+1)/2$ are binary variables whereas

[TABLE]

The energy offset reads $H(\boldsymbol{s})-F(\boldsymbol{q})=\sum_{i=1}^{N}h_{i}-\sum_{\langle i,j\rangle}J_{ij}$. Note, if a given $q_{i}$ vanishes so does any product $q_{i}q_{j}$. Therefore, QUBO formulation (2) effectively reduces the number of multiplications almost by half in comparison to Eq. (1).

Despite its straightforward formulation, the problem of solving spin–glass instances can not be easily tackled using a brute force approach even for a modest number of spin variables. This is since the number of possible spin assignments grows exponentially with the number of nodes in the graph. For instance, when $N=40$, the number of possible states is greater than the number of bits in a $32$GB memory chip. Already when $N=64$, the size of the search space is greater than the estimated age of the Universe in seconds [34]. In fact, the problem of finding the ground state of the Ising model defined on an arbitrary graph is long known to be NP–hard [35]. This means, in particular, that even verifying if a given configuration minimizes the cost function (1) is difficult.

3 Description of the algorithm

A general idea underlying this work is to perform an exhaustive search over the whole state space, taking advantage of massive parallel capabilities of modern GPUs. This requires an efficient strategy to encoding all states, $\boldsymbol{q}=(q_{1},q_{2},\dots,q_{L})$, on a GPU. A naive approach would required storing an array of $N$ integers, $q_{i}=0,1$, for each state $\boldsymbol{q}$. However, this would also lead to excessive use of memory and render this approach inefficient. As an optimal strategy, one should try to reuse information already stored in the GPU memory. Therefore, in our algorithm, we take advantage of the following correspondence

[TABLE]

where $k=(\boldsymbol{q})_{2}$ denotes the binary representation of an integer $k$. For instance, when there is $N=8$ spins, one may associate

[TABLE]

Theoretically, this strategy allows one to store $M=2^{64}\sim 10^{19}$ states with no extra cost, limiting the system size to $N=64$ spins. Nonetheless, this is more than the current architecture, based on the von Neumann paradigm of computation, which can process in a reasonable time [36]. Indeed, we estimated that optimal search among $2^{64}$ states to extract the low energy spectrum consisting of $S=10^{2}$ of them would take $821$ years on an efficient Titan V GPU [37]. In comparison, systems of sizes $N=32$, $49$, $50$ can be solved within $5$ seconds, $12$ and $24$ days, respectively. A detailed benchmark is presented in Sec. 5.

One should stress that the fastest (as of $2018$) supercomputer in the world—Summit—is equipped with $27648>2^{14}$ Nvidia Tesla V$100$ GPUs [38]. Therefore, “only” $2^{14}$ of them (processing chunks of size $2^{50}$ each, simultaneously) should be able to reduce the number of $821$ years (for a single GPU tackling $N=64$) substantially. Perhaps, maybe even down to a couple of months. Nevertheless, a priori, it is hard to estimate the exact numbers due to various communication bottlenecks. This interesting open problem is, however, beyond the scope of the current work.

In theory, one could first compute all $M=2^{N}$ energies in parallel and only then select $S\ll 2^{N}$ lowest ones (and the corresponding states if needed). However, even with an efficient storage strategy, this approach quickly becomes impractical for large systems. It requires an exponentially increasing storage space to encode possible solutions. To overcome this problem, one could iterate over the solution space in manageable chunks, each time extracting the desired number of states [e.g., with the bucket select algorithm [39], for which the mean execution time scales linearly with the size of the input vector]. Sorting the energies is executed only in the final step. Since GPU threads and blocks are labeled in the same way for every chunk, an offset is required to correctly enumerate all states, i.e.,

[TABLE]

Note, the energy calculations are independent and thus can be performed in parallel. The overall parallel speedup is limited by the serial part (Amdahl’s law [40]) consisting of the lowest energy states extraction and merging all local information into the global record. Algorithm 1 in the below listing summarizes the underlying structure of our solver.

4 Implementation details

4.1 Languages and technologies employed

The core components of our implementation has been written in modern Fortran $95/2003$ [41], which we have chosen for its flexibility [42], extensive support for linear algebra [43], performance [44] and native support for CUDA technology [45]. To make our code easier to use, we have wrapped it in a Python package using the f2py [46] utility and numpy’s fork of distutils package [47]. Whereas Fortran is widely used mostly for numerical simulations [48], Python is one of the most popular general–purpose programming language [49].

We have also incorporated the fast $k$–selection algorithm for GPU [39], and the Thrust library [50] into our solver for its parallel implementation of many standard methods such as finding the minimum and maximum of an array or partitioning thereof. The Thrust library is utilized both for the GPU implementation and for the pure CPU implementation with the OMP backend. In order to take advantage of various Thrust’s functions, we have written several small C++ modules to bind them into the Fortran code. The source code of the entire package, together with the comprehensive documentation, can be found on GitHub [31].

The Python wrapper allows one to execute the algorithm both on the CPU and GPU. It was designed with simplicity in mind, and as such, its primary usage does not require any specialized knowledge. A basic understanding of the underlying optimization problem is enough, cf. the listing below. In particular, only the system definition (i.e., the graph or adjacency matrix) and the desired number of states needs to be provided by the user. Nevertheless, other parameters, including the chunk sizes, can also be passed to the wrapper.

On virtually all Linux platforms, it is possible to install the very basic version (i.e., with no GPU support) of our solver directly from the Python Package Index, by issuing

[TABLE]

where ising is the name of the package. However, to assure full compatibility with modern GPUs, CUDA requires a custom build from source which can be initiated via

[TABLE]

from the package source directory. For more details regarding custom installation, including CUDA and various Fortran compilers, we refer the reader to documentation [51]. Note, only PGI and IBM XLF support CUDA Fortran. Our package has only been tested using the former.

4.2 GPU execution scheduling

Programming GPUs often pose a nontrivial endeavor. Among many challenges, one has to design the grid on which kernels are launched [52]. We have tested various grid/blocks launching configurations for the energy computing kernel and have obtained the best results with grids consisting of $G=2^{g-5}$ ($g$ being the current chunk size) blocks of size $B=2^{5}$ each. This particular value maximizes the Local Store Throughput (LST)—reported by the NVIDIA Profiler, nvprof—and thus also the execution time of the energy computing kernel, cf. Fig. 4.

Note also that the Titan V is roughly $6$ times faster than GTX $1080$, which is exactly the ratio between LSTs for these two devices. Such apparent correlation validates, to some extent, the LST as a proper metric (in comparison to a typical Throughput used to benchmark GPUs) to asses the performance of our solver.

We would like to stress, however, that the optimal launching configuration we have used in the present studies may need further (experimental) adjustment depending on the hardware, and possibly the problem size.

4.3 Complexity analysis

Our algorithm performs an exhaustive search over the entire, exponentially large, state space in predefined chunks to find $S$ lowest states (cf. Algorithm 1). Thus, unavoidably its time complexity has to be at least exponential in the system size $N$. Computing the energy (2) for a single state, $\boldsymbol{q}$, requires $\mathcal{O}(N^{2})$ operations. The selection procedure executed on a data chunk of size $2^{k}$, however, requires $\mathcal{O}(2^{k})$ comparisons resulting in $\mathcal{O}[2^{k}(N^{2}+1)]$ operations. Finally, taking into account the total number of chunks, $2^{N-k}$, and adding complexity of the final sorting procedure, $\mathcal{O}[S\log(S)]$, results in total complexity being

[TABLE]

Therefore, essentially the solver’s complexity behaves as $\mathcal{O}(2^{N})$ which we also demonstrate experimentally in Sec. 5 (cf. Fig. 3).

As one can see, the GPU implementation takes $30$ seconds (GeForce $1080$) and $5$ seconds (Titan V) on average to solve the Ising problems with $N=32$ spins. The same problem requires about $1500$ seconds on average on a single CPU core. For GPU, the differences in solution times between single and double precision are close to $10\%$ and are not reported on Fig. 3.

5 Benchmarks

5.1 GPU vs CPU comparison

We have tested our algorithm on the following hardware:

CPU: $10$ Cores ${\rm Intel}^{\rm R}$ ${\rm Core}^{\rm TM}$ i7-$6950$X;

• 2.

GPU(1): Nvidia GeForce GTX $1080$, $8$GB GDDR$5$ global memory, $2560$ CUDA Cores;

• 3.

GPU(2): Nvidia Titan V, $12$GB HBM$2$ global memory, $5120$ CUDA Cores.

For benchmarking purposes, we have executed our algorithm on a fully connected, $K=100$ randomly generated (cf. Ref. [53, 54]), problem instances for systems up to $N=50$ (on Titan V). For each instance, we have calculated the low energy spectrum consisting of $S=10^{2}$ states in a single run. Typical results obtained with a high-end CPU (i7-$6950$X) and both a mid-class (GeForce $1080$) and professional (Titan V) GPU are depicted in Fig. 3. We have also estimated time to solution experimentally, for larger systems (up to $N=50$ spins) for which the low energy spectrum can be obtained in a reasonable time (i.e., one month) on Titan V. The estimate is based on the average time required to process a single chunk of data [of size $M=2^{29}$ (CPU), and $M=2^{27}$ (GTX $1080$)]. Our measurements are consistent with the complexity analysis discussed in Sec. 4.3.

5.2 Validation of MPS algorithm

To demonstrate the capabilities of our solver, we employ it to benchmark a more sophisticated, heuristic, approach based on a Matrix Product States (MPS) technique (see Supplementary Material of Ref. [ 23] for details). Here, we are not interested in time to solution, but instead, we would like to investigate the accuracy of the latter. Heuristic algorithms can often solve large systems ($N\gg 50$). However, they cannot certify solutions.

With the MPS based algorithm one aims at approximating the Boltzmann distribution,

[TABLE]

for a sufficiently large inverse temperature, $\beta$, where each $A^{s_{i}}$ ($i=1,2\ldots,L$) is matrix of limited dimensions $\leq D$ (refereed to as the bond dimensions). The above approximation is usual depicted—using a network of tensors—as

[TABLE]

At each bond, one splits the system into two–halves. The exact decomposition would require the bond dimension $D$ to grow exponentially with the number of spins in one half, interacting with spins in the second half (and arbitrary numerical precision). The limited bond dimension $D$ reflects on the amount of entanglement/correlations (related to a given bipartition), which can be stored in a “quantum system” decomposed as MPS [55]—here we understand a “quantum system” as a superposition over all possible classical spin configurations. Having the approximation in Eq. (10) in the form of MPS, we can efficiently calculate any marginal and conditional probability (at the inverse temperature $\beta$) described by $|\Psi(\beta)\rangle$, and then systematically search for the most probable classical configurations (i.e., the ones with the smallest energies) using branch and bound strategy—building the most probable spin configurations one spin at the time.

Finally, to perform the search one needs to find $|\Psi(\beta)\rangle$, which is obtained by starting from $\beta=0$—for which the MPS decomposition $|\Psi(\beta=0)\rangle$ is trivial—and then subsequently simulating the imaginary time evolution (i.e., the annealing). To that end, we apply the sequence of operators,

[TABLE]

which amount to $\prod_{i=1}^{N}U_{i}(d\beta)=e^{-d\beta H(\boldsymbol{s})/2}$. Applying each gate (11) results in doubling of the affected bond dimensions. Moreover, applying all such operators would result in uncontrollable, exponential growth of the MPS matrices. However, the one–dimensional (and loop–free) structure of the MPS ansatz allows one to systematically, at each step, find its approximation, which effectively compresses the information and maintains the bond dimensions limited to $D$. The whole procedure can be graphically depicted as

[TABLE]

While all the applied operators $U_{i}(d\beta)$ formally commute (independent of $d\beta$), due to the finite numerical precision and finite $D$, it is relevant to reach the final inverse temperature, $\beta$, gradually in a couple of consecutive steps, each with smaller $d\beta$. Otherwise, for larger $d\beta$, $U_{i}(d\beta)$ effectively act as projectors trapping the system in a local minima.

The question then becomes how well the MPS ansatz, which by construction is one–dimensional, is able to encode the structure of low energy spectrum for fully–connected graphs. In general, the bigger the system, the higher $D$ necessary to faithfully capture the structure of low energy spectrum. We observe that already moderate $D$ of $128$ is enough to find all ground states for $100$ considered instances, see Fig. 5a. The inverse temperature $\beta=1$ is large enough to sufficiently zoom–in on the low energy states. At the same time, the importance of small enough time–step (here $d\beta=0.25$) is clearly visible. It is also enough to recover most of the $1000$ configurations with lowest energies for those instances, see Fig. 5b for $N=40$. Note that the exact MPS decomposition would require the bond dimension of $2^{N/2}=2^{20}$. This demonstrate the magnitude of the compression of the relevant information encoded in MPS.

A typical lowest energy spectrum for $N=50$ spins, and consisting of $S=10^{3}$ states, is shown in Fig. 6. Therein, we have also incorporated an approximated low–energy spectrum obtained with the recent Monte Carlo based algorithm, introduced in [56], which can determine the density of states. The numerical data was provided to us by the authors of that paper.

6 Summary

We have demonstrated how to perform an exhaustive (brute–force) search in the solution space of the Ising spin–glass model [2] utilizing modern Graphics Processing Units [57]. Our algorithm can also be adapted for different heterogeneous architectures (e.g., Xeon Phi [58]). The Hamiltonian of this particular model encodes a variety of important optimization problems [5]. Moreover, this model has also been realized experimentally as a commercially available D-Wave quantum annealer [8].

Our implementation with CUDA Fortran [45] offers unprecedented speed and efficiency already visible on commodity hardware (e.g., GeForce $1080$). Furthermore, it can be easily tuned for professional GPUs such as Titan V [37] virtually at no extra effort. To give an example, our algorithm, when tailored for the latter GPU, can extract the low energy spectrum (consisting of $N=10^{3}$ states) in roughly $5$ seconds for the spin system admitting $M=2^{32}\sim 10^{9}$ different configurations. In comparison, a single CPU core takes (on average) $25$ minutes to finish the same task (cf. Sec. 5 for detailed benchmark).

Admittedly, practical applications of brute–force algorithms are constrained to small problem sizes ($N\leq 50$). However, they can not only solve the spin–glass problems for arbitrary topologies and instances but also certify solutions [24, 17, 18]. These two features are crucial for developing and validating new methods and strategies for solving classical optimization problems [26]. We have explicitly exemplified this point by comparing our algorithm to a sophisticated recent Ising solver based on tensor network techniques [23]. In particular, we have demonstrated that despite its one–dimensional nature, the Matrix Product State ansatz is still able to approximate well the relevant part of the Boltzmann distribution for a fully connected graph of $N\leq 50$. Therefore, this suggests that the MPS algorithm should be superior for all problems having a dominant quasi–one–dimensional nature that allows for sparse connections to span the full problem.

Finally, to benefit the community, we have made our code publicly available as an open–source project [31]. Moreover, for those users who lack technical knowledge of Fortran or CUDA, we have provided an easy to install and use Python wrapper [31].

Acknowledgments

We appreciate fruitful discussions with Andrzej Ptok, Jerzy Dajka, and Piotr Gawron and Masoud Mohseni. We thank Paweł Wasiak for his valuable remarks regarding solver’s documentation. We gratefully acknowledge the support of NVIDIA Corporation with the donation of the Titan V GPU used for this research. This work was supported by National Science Center (NCN, Poland) under projects 2015/19/B/ST2/02856 (KJ) 2016/20/S/ST2/00152 (BG) and NCN together with European Union through QuantERA ERA NET program 2017/25/Z/ST2/03028 (MMR). MMR acknowledges receiving Google Faculty Research Award 2017.