Equation Planting: A Tool for Benchmarking Ising Machines

TL;DR

This paper presents 'equation planting', a novel method for creating benchmark problems for Ising machines by transforming linear systems into Ising models, enabling better evaluation and understanding of these devices.

Contribution

The paper introduces 'equation planting', a new technique to generate challenging benchmark problems for Ising machines from linear systems, aiding in performance assessment and mechanistic studies.

Findings

Equation planting creates NP-hard-like problems from linear systems.

The method helps evaluate Ising machines as optimizers and samplers.

It provides insights into the operation of Ising devices.

Abstract

We introduce a methodology for generating benchmark problem sets for Ising machines---devices designed to solve discrete optimization problems cast as Ising models. In our approach, linear systems of equations are cast as Ising cost functions. While linear systems are easily solvable, the corresponding optimization problems are known to exhibit some of the salient features of NP-hardness, such as strong exponential scaling of heuristic solvers' runtimes and extensive distances between ground and low-lying excited states. We show how the proposed technique, which we refer to as `equation planting,' can serve as a useful tool for evaluating the utility of Ising solvers functioning either as optimizers or as ground-state samplers. We further argue that equation-planted problems can be used to probe the mechanisms underlying the operation of Ising machines.

Figures

Equations

$$\sum _{j=1}^n{a}_{ij}{x}_j=b_i\text{for}i=1\dots m.$$

$$F=\sum _{i=1}^m{\Vert \sum _{j=1}^n{a}_{ij}{x}_j-b_i\Vert }^2.$$

$$\sum _{j=1}^n{a}_{ij}{x}_j=b_i\mathrm{mod}2\text{for}i=1\dots m,$$

$$\prod _{j:a_{ij}=1}{s}_j=(-1{)}^{b_i}.$$

$${\tilde{F}}_2=\sum _{i=1}^m{\left(\prod _{j:a_{ij}=1}{s}_j-(-1{)}^{b_i}\right)}^2,$$

$$F_2=-\sum _{i=1}^m(-1{)}^{b_i}\prod _{j:a_{ij}=1}{s}_j.$$

$$G_{\text{3X}}$$

Full text

Introduction.— Recent years have witnessed a flourishing in experimental ‘Ising machines’—special-purpose programmable devices engineered to solve discrete optimization problems Korte and Vygen (2012); Papadimitriou and Steiglitz (2013) cast as Ising models—with the tacit promise that their performance is superior to those of algorithms running on standard computers. Analog quantum devices that perform quantum annealing Kadowaki and Nishimori (1998); Farhi et al. (2001) designed to find bit assignments that minimize the cost of Ising Hamiltonians have already been realized on various platforms such as arrays of superconducting flux qubits Tolpygo et al. (2015a, b); Jin et al. (2015); Johnson et al. (2010); Berkley et al. (2010); Harris et al. (2010); Bunyk et al. (Aug. 2014). Other notable technologies that recently emerged are coherent Ising machines based on lasers and degenerate optical parametric oscillators Yamamoto et al. (2017); McMahon et al. (2016), FPGA-based quantum-inspired digital annealers Aramon et al. (2018), and memcomputing devices that operate on terminal-agnostic self-organizing logic gates Di Ventra and Traversa (2018); Traversa and Di Ventra (2018); Traversa et al. (2015).

Touting improved performance over traditional algorithms Haribara et al. (2017); Inagaki et al. (2016), Ising machines have gained a large amount of interest, both in the academe as well as with the general public—and rightfully so. Many problems of theoretical and practical relevance, in areas as diverse as machine learning, materials design, software verification, portfolio management, logistics, and many more, can be cast as searching for the global minima of Ising cost functions Korte and Vygen (2012); Papadimitriou and Steiglitz (2013). It is, however, unclear to date whether indeed any one of the aforementioned devices truly offers genuine advantages over its competitors.

One of the main bottlenecks preventing meaningful benchmarking of Ising machines is the absence of appropriate benchmark problems. Generating problem instances that are on the one hand challenging enough and whose solutions are on the other hand known in advance or are easily verifiable, so as to allow proper benchmarking, are not straightforward to construct since these two requirements are, in many respects, contradictory. To overcome this obstacle, various methods have been devised in recent years to generate problem classes with known minimizing configurations (see, e.g., Refs Hen et al. (2015); King et al. (2019); Marshall et al. (2016)). However these have generally been found to lack the hardness that characterizes NP-hard problems. At the other extreme, problem classes that are challenging to solve but whose ground-state energies are not verifiable, have also been developed Mandrà and Katzgraber (2018); Martín-Mayor and Hen (2015); Venturelli et al. (2015).

To bridge the above gap, in this study we advocate a methodology for generating benchmark problem sets for Ising machines that we argue may serve as a suitable tool for evaluating their performance and in turn also to indirectly probe the mechanisms underlying their operation. The proposed technique, which we refer to as ‘equation planting,’ has several desirable properties. Explicitly, the generated problem sets, while generated from easily solvable problems, possess some of the salient features of NP-hard problems—most notably extensive distances between ground and low-lying excited states.

Our approach is based on the embedding of linear systems of equations in Ising models, motivated by the observation that when linear systems of equations are cast as optimization problems, the latter form often stymies heuristic solvers Franz et al. (2001); Barthel et al. (2002); Jörg et al. (2010); Ricci-Tersenghi (2010); Guidetti and Young (2011); Farhi et al. (2012). In addition to serving as suitable benchmark problems to Ising machines functioning as optimizers, we also show that the proposed technique may be used to test the functionality of these devices as Boltzmann samplers, or more specifically as ground-state samplers Adachi and Henderson (2015); Benedetti et al. (2016); Marshall et al. (2017), thanks to another property of these problem sets—namely, a verifiable number of ground state configurations.

We begin by describing the equation planting technique in general, moving on to focus on a specific class of problem instances using which we demonstrate the effectiveness of the method.

Equation planting.— We begin by considering a linear system of $m$ equations in $n$ variables

[TABLE]

Here $\{x_{1},\ldots,x_{n}\}$ stand for variables over a given field, and $\{a_{ij}\}$ and $\{b_{i}\}$ are the equation coefficients. Every such linear system may be cast as an optimization problem with the corresponding cost function

[TABLE]

Since $F$ is a sum over positive terms, any configuration $\{x^{*}_{1},\ldots,x^{*}_{n}\}$ that yields $F=0$ is a minimizing configuration that also solves the linear system.

As already mentioned, even though the computational problem of solving a linear system of equations is easy—any given instance can always be solved in polynomial time using Gaussian elimination—the corresponding optimization problem is not necessarily easy for heuristic solvers Franz et al. (2001); Barthel et al. (2002); Ricci-Tersenghi (2010); Guidetti and Young (2011); Farhi et al. (2012).

We leverage the above setup towards constructing Ising problems that stymie heuristic solvers. To that aim, we focus henceforth on linear systems of equations modulo 2, that is,

[TABLE]

where now both the variable set and coefficients are Booleans taking on values $\{0,1\}$. The $i$th equation can therefore be written as $\bigoplus_{j:a_{ij}=1}{x_{j}}=b_{i}$ (here $\oplus$ denotes the bitwise XOR operation), or in terms of Ising spins $s_{i}\in\{-1,1\}$,

[TABLE]

Cast in optimization form, the system of equations becomes

[TABLE]

or, after dropping immaterial constants,

[TABLE]

The cost function $F_{2}$ is a multi-spin cost function consisting of a sum of products of spin variables. The localities of the terms comprising the cost function are precisely the number of variables in the corresponding equations. The minima of $F_{2}$ correspond to the solutions of the system given in Eq. (4), provided that solutions exist. In more detail, linear systems may behave in one of three possible ways. (i) The system may have no solutions at all. This scenario corresponds to the cost function $\tilde{F}_{2}$ having a strictly positive minimal value (meaning $F_{2}>-m$) and may happen if the number of equations is greater than the number of variables. (ii) The system may have a unique solution. Here, a single configuration minimizes $F_{2}$ whose minimal value in this case would be $-m$. (iii) Last, the system may have multiple solutions if the dimension of its null space, or nullity Lay (2005), which we denote here by $d_{N}$, is nonzero. In this case, the number of minimizing configurations that yield $F_{2}=-m$ is $N_{\text{GS}}=2^{d_{N}}$.

Ising machines are designed to tackle two-body cost functions. The locality of the multi-spin cost $F_{2}$ may in this case be readily reduced to a two-body Ising model of the general form $\sum_{ij}J_{ij}s_{i}s_{j}+\sum_{i}h_{i}s_{i}$ using standard reduction techniques (see, e.g., Ref. Babbush et al. (2013)), where $\{J_{ij}\}$ and $\{h_{i}\}$ are integer-valued coefficients.

As we demonstrate in the next section, while the solutions of the linear systems of equations used to generate these Ising cost functions are straightforwardly obtained (if there are any), heuristic solvers will generically find this type of problems extremely difficult to solve. In what follows, we illustrate this by studying in detail one specific type of equation system—namely, random 3-regular 3-XORSAT.

3-regular 3-XORSAT as 2-body Ising.— To demonstrate the challenges that equation-planted instances present to heuristic Ising solvers, we consider as a test case linear systems $\mod 2$ (also known as XORSAT equations) wherein each equation contains exactly three randomly chosen Boolean variables and each variable, or bit, $x_{j}$ appears in exactly three equations. This type of problem is also known as ‘3-regular 3-XORSAT’ (3R3X) and was studied previously in various contexts Jörg et al. (2010); Guidetti and Young (2011); Farhi et al. (2012); Albash et al. (2017, 2019). An $n$-bit instance would thus consist of $n$ equations of the form $x_{i_{1}}\oplus x_{i_{2}}\oplus x_{i_{3}}=b_{i}$(see the Appendix for additional details). We note that while the results presented below pertain to 3R3X, they should extend in general to other classes of linear XORSAT systems.

The cost function $F_{2}$ of an $n$-bit 3R3X system is therefore a sum of $n$ three-body terms of the form $-(-1)^{b_{i}}s_{i_{1}}s_{i_{2}}s_{i_{3}}$ defined on $n$ Ising spins. To reduce the locality of a term to two- and one-body we use a gadget that shares its four minimizing configurations (see the Appendix for additional details). For the negatively signed clause ($b_{i}=0$), these are the four configurations whose product is positive, namely, $(1,1,1),(1,-1,-1),(-1,1,-1)$ and $(-1,-1,1)$ and an appropriate gadget is

[TABLE]

where $(h,\tilde{h},J,\tilde{J})$ can be either $(-1,-2,1,2)$ or $(-1,2,1,-2)$, yielding in both cases a minimal cost of $-4$ (we note that other gadgets exist that can equivalently be used). The above gadget is a fully connected 4-spin cost function with $s_{a}$ serving as an auxiliary spin (see Fig. 1). The positively signed clauses are treated similarly.

The above gadget allows us to cast $n$-bit 3R3X linear systems as two-body Ising models with $2n$ spins, as each clause adds one auxiliary spin to the problem and an $n$-bit instance has $n$ clauses.

At this point, we note several of the 3R3X Ising instances properties: (i) First, these instances have a bounded degree, namely 9, as each spin connects in general to three other spins for every clause in which the spin appears. (ii) The specification of Ising parameters $\{J_{ij}\}$ and $\{h_{i}\}$ requires only two bits of precision as both spin-spin couplings and longitudinal fields are integers in the range $[-3,3]$. (iii) As already mentioned, these instances have known ground-state energies. (iv) Finally, these instances have a controllable ground-state degeneracy, determined by the nullity of the generating linear system. As will become evident in the next section, these properties may be leveraged to determine the utility of heuristic Ising solvers functioning either as optimizers or as ground-state samplers.

Results.— In what follows, we examine the performance of heuristic solvers on randomly generated two-body Ising 3R3X instances that have unique solutions. As a representative for state-of-the-art heuristic solvers, we use parallel tempering (PT) Hukushima and Nemoto (1996); Marinari (1998) — a refinement of the celebrated yet somewhat outdated simulated annealing algorithm Kirkpatrick et al. (1983), that finds the ground-state configurations of general discrete-variable cost functions Fernandez et al. (2013); Fernandez and Martin-Mayor (2009); Fernandez et al. (2009); Ballesteros et al. (2000); Katzgraber et al. (2006); Earl and Deem (2005); Martín-Mayor and Hen (2015); Zhang et al. (2017); Hen et al. (2015); Albash et al. (2017). In PT, multiple copies of the problem are equilibrated in parallel at different temperatures and spin configurations at adjacent temperatures are regularly swapped (see the Appendix for additional details). In addition, we consider a variant of PT (which we denote PT-H) that employs global cluster moves due to Houdayer Houdayer (2002) and is known to accelerate PT convergence rates for many problem classes (see, e.g., Refs. Katzgraber et al. (2004); Houdayer and Hartmann (2004); Albash and Lidar (2018)).

We test the performance of PT and PT-H on instances of varying sizes by measuring their typical times to find minimizing configurations. Here, typical runtimes are defined as the median time to reach a minimizing configuration over 100 randomly generated instances of a given size 111Both algorithms were run on one core of a $3.5$ GHz 6-Core Intel Xeon E5 processor..

The results are presented in Fig. 2(top) showing the exponential runtime scaling of both PT and PT-H in problem size, proportional to ${\text{e}}^{\alpha n}$ with $\alpha\approx 0.13$ and $0.14$, respectively, indicating a strong exponential scaling. Moreover, the Houdayer cluster updates do not seem to improve the scaling. Also shown in the figure (as shaded areas) are the interquartile ranges (the difference between 75th and 25th percentiles) of the PT and PT-H runtimes, quantitatively indicating the homogeneity in the hardness of the instances.

To better understand why these problems, while being trivial to solve using Gaussian elimination, present severe challenges to heuristic solvers, we next examine their energy landscapes. We do so by measuring the Hamming distances between local minimum configurations, as found by the solver in the course of the simulation, and the global minimum. These are plotted in Fig. 2(bottom). The low-lying excited states, characterized by small yet positive residual energies, are typically $\sim 0.6n$ spin flips away from the global minimum, implying that reaching a local minimum is no indication of the whereabouts of the ground state. The above property is a hallmark of NP-hard problems Papadimitriou and Steiglitz (2013).

We next illustrate the usefulness of equation planting in determining the bias of Ising machines tasked with sampling the ground-state manifolds of Ising cost functions Bastea and Duxbury (1998); Martín-Mayor and Hen (2015); Matsuda et al. (2009), which in addition to optimization is considered to be one of their main suggested uses Adachi and Henderson (2015); Benedetti et al. (2016); Marshall et al. (2017); Zhang et al. (2017).

We first construct random Ising-encoded 3R3X instances with nullities $d_{N}=1,2,3$ having ground-state degeneracies of $2,4$ and $8$, respectively, at different problem sizes. We then study the functionality of PT-H as a ground-state sampler by simply measuring the fraction of occurrences of each ground state for any given instance [see Fig. 3(left)]. Each instance is solved 100 times, and the ground-state configurations are recorded.

We quantify the bias of PT-H on the various instances by measuring the typical $p$-value obtained from one-sided $\chi^{2}$ tests performed on the tallied ground states of each instance under the hypothesis that the true distribution is uniform. The results are summarized in Fig. 3(middle). The smallness of the $p$-values indicates PT’s strong bias in sampling ground-state manifolds of these instances, implying differently sized ‘basins of attraction’ for each ground state. We also find that, similar to the ground state manifolds of NP-hard problems, the distances between the ground states of any given instance are typically extensive. These are depicted in Fig. 3(right) 222The bias exhibited by PT may seem at first to contradict the Boltzmann distribution which (in accordance with the assumption of equal a priori probability) prescribes the same probability to same-cost configurations for thermalized systems. However, we note that when PT is used as an optimizer (as it is here), the number of sweeps is not large enough for thermalization to take place..

Summary and conclusions.— We introduced a method for generating random problem sets that we argued are optimally suitable for verifiably testing the utility of Ising machines. The proposed technique allows for the construction of random problem sets that have a number of properties that are desirable for this task: (i) verified ground-state energies; (ii) controllable ground-state degeneracy; (iii) ‘NP-hard-like’ energy landscapes; (iv) low-precision problem parameters and (v) underlying connectivity with a bounded degree.

We further illustrated the computational difficulties that heuristic state-of-the-art solvers face when tasked with finding or uniformly sampling the minimizing configurations of problem sets possessing intricate energy landscapes. We showed that equation-planted instances become difficult to solve already at small sizes, revealing rather early on the asymptotic exponential scaling expected from heuristic solvers. In addition, we noted that, as is expected from difficult problem sets, Houdayer cluster moves Houdayer (2002) do not provide any additional advantages.

With the above stated, it should be noted that since the problem sets are generated from trivially solvable systems, special-purpose algorithms could conceivably be tailored to identify the underlying equations within the Ising cost function and utilize the equation structure to solve these Ising costs polynomially fast. On the other hand, the equation structure can be easily obfuscated if one (i) employs a wider range of gadgets to embed the various XORSAT clauses, (ii) constructs instances using a randomly weighted sum of terms, (iii) generates non-regular linear systems of equations with potentially different numbers of variables per equation or (iii) adds to existing XORSAT cost functions non-XORSAT gadgets whose ground states coincide with those of the instance, as these are known in advance. An interesting question that arises in this context is whether one can obfuscate such easily solvable problems in a manner that is provably intractable to ‘reverse engineer’ back to easily solvable form.

In the generation of the random XORSAT instances considered in this study, no specific connectivity constraints were assumed, which in turn led to the generation of Ising instances with random connections among the various spins. These instances could not be directly programmed into Ising machines that have rigid and sparse architectures Marshall et al. (2017). This limitation may be addressed in one of two ways. The first is to convert the randomly connected instances to match the Ising machine hardware graph via embedding schemes at the price of adding auxiliary spins to the problem Cai et al. (2014). Alternatively, one can directly generate equation-planted instances that are native to the device connectivity, which may require devising novel 3-XORSAT gadgets.

Another aspect of equation planting that has not been explored here is the relative hardness of $r$-regular $k$-XORSAT instances for different $k,r$ greater than three — problem classes that are expected to yield yet harder instances. Varying $k,r$ would presumably allow for ‘hardness tunability’, which is often a desired property in the benchmarking of heuristic solvers. Another way for achieving hardness tunability that was not pursued here is the use of non-regular random graphs, such as Poisson graphs, rather than random regular graphs Barthel et al. (2002).

The NP-hard-like energy landscapes characterizing equation-planted instances coupled with the fact that their ground state manifolds can be known in advance, may also be used as a tool for gaining insight into the mechanisms underlying the operation of Ising machines. This is because the machines’ performance on these problems serves as an indicator of the locality of the heuristic being employed. As we demonstrated above, local search approaches are not very successful whereas global approaches such as Gaussian elimination are able to decipher their inherent structure very easily. Another possible use of this type of problems is perhaps also the study of spin glasses, specifically the onset of temperature chaos Fernandez et al. (2013); Aspelmeier et al. (2002); Berthier and Young (2005).

One intriguing future research direction would be the development of heuristic solvers that identify XORSAT-type relations between the spins of Ising cost functions and utilize the ease with which such sub-problems can be minimized in order to efficiently find global minimal costs.

Appendix A Generating 3-regular 3-XORSAT instances

The simplest way to generate a uniform ensemble of $n$-spin 3-regular 3-XORSAT instances is by generating triplets of lists ${\bf v}^{(1)},{\bf v}^{(2)}$ and ${\bf v}^{(3)}$ each of which is a random shuffling of the list $(1,\ldots,n)$. As a next step, the 3-XORSAT clauses, or equations, are constructed from the $n$ triplets $C_{i}=({\bf v}^{(1)}_{i},{\bf v}^{(2)}_{i},{\bf v}^{(3)}_{i})$ with $i=1\ldots n$. If any of the clauses are found to contain repeated indices, the process is reinitialized. The triplets of indices $C_{i}$ can then be used to form XORSAT equations $x_{{\bf v}^{(1)}_{i}}+x_{{\bf v}^{(2)}_{i}}+x_{{\bf v}^{(3)}_{i}}={\bf b}_{i}\mod 2$ where ${\bf b}$ is a randomly generated vector of $n$ bit values. These equations can also be written as $x_{{\bf v}^{(1)}_{i}}\oplus x_{{\bf v}^{(2)}_{i}}\oplus x_{{\bf v}^{(3)}_{i}}={\bf b}_{i}$ where $\oplus$ is the XOR operation (hence the name XORSAT). This routine ensures that every index appears in exactly three clauses, i.e., it is 3-regular. The degeneracy of the Ising instance generated from the above XORSAT system of equations can be found by row-reducing the system to obtain the dimension of its null space. Similarly, $r$-regular $k$-XORSAT instances can be generated for different choices of $k$ and $r$.

Appendix B Gadget construction

The gadgets for this study were constructed using a brute-force search of the parameter space of Ising cost functions, namely sweeps overs ranges of $J_{ij}$ couplings and longitudinal fields $h_{i}$. While surely more efficient (and more intelligent) methods exist for constructing such gadgets, since the gadgets need only be constructed once, such methods are not particularly necessary.

To construct an Ising cost function for a $k$-XORSAT equation using $a$ auxiliary spins we consider fully-connected Ising models with $k^{\prime}=k+a$ spins. Beginning with $a=0$, we systematically scan through sets of integer-valued couplers and external fields (starting with $J_{i},h_{i}=0,\pm 1,\pm 2,\ldots$). For each choice of cost function parameters, we enumerate the ground-state configurations of the resultant model. Specifically, we test whether the ground states satisfy the XORSAT conditions, namely, that there are $2^{k-1}$ of them and that all have the same parity over the non-auxiliary spins, i.e., $\prod_{j=1}^{k}s_{k}=(-1)^{b}$ where $b$ is either 0 or 1. If no choice of parameters is found, we conclude that the number of auxiliary spins $a$ is too small. We then increase $a$ by one and restart the routine.

Appendix C The parallel tempering simulations

We outline here the technical details of the parallel tempering (PT) simulations Hukushima and Nemoto (1996); Marinari (1998) used in this study. In PT, one considers $N_{T}$ independent copies of the system executing a Metropolis algorithm in parallel at different temperatures, $T_{1}<T_{2}<\ldots<T_{N_{T}}$. Copies with neighboring temperatures regularly attempt to swap their temperatures with probabilities that satisfy detailed balance Sokal (1997). In this way, each copy performs a temperature random-walk. At high temperatures, free-energy barriers are easily overcome, allowing for a global exploration off configuration space. On the other hand. at lower temperatures the local minima are explored in more detail.

The temperature grid of our PT simulations consisted of $N_{T}=37$ inverse-temperatures $\beta$ in the range $(0.0166667,3.33333)$ [see Fig. 4(left)]. The maximal inverse temperature was chosen such that $\beta_{\max}|J_{\max}|=10$ to ensure an effectively zero final temperature. Ergodicity was maintained by the temperature-swap part of the PT.

To verify an optimized choice of the temperature mesh, the acceptance probabilities of swaps between neighboring temperatures were closely monitored to make sure that there are no bottlenecks in the flow of configurations in temperature space, i.e., that the acceptance ratios are not too small (or too high). Figure 4(right) shows the swap acceptance rates for several randomly chosen instances, indicating the effectiveness of the method.