# From Near to Eternity: Spin-glass planting, tiling puzzles, and constraint satisfaction problems

###### Abstract

We present a methodology for generating Ising Hamiltonians of tunable complexity and with a priori known ground states based on a decomposition of the model graph into edge-disjoint subgraphs. The idea is illustrated with a spin-glass model defined on a cubic lattice, where subproblems, whose couplers are restricted to the two values , are specified on unit cubes and are parametrized by their local degeneracy. The construction is shown to be equivalent to a type of three-dimensional constraint satisfaction problem known as the tiling puzzle. By varying the proportions of subproblem types, the Hamiltonian can span a dramatic range of typical computational complexity, from fairly easy to many orders of magnitude more difficult than prototypical bimodal and Gaussian spin glasses in three space dimensions. We corroborate this behavior via experiments with different algorithms and discuss generalizations and extensions to different types of graphs.

###### pacs:

75.50.Lk, 75.40.Mg, 05.50.+q, 03.67.Lx## I Introduction

Hard optimization problems are ubiquitous throughout the sciences and engineering, and have consequently been the subject of intense theoretical and practical study. While it is generally considered highly unlikely that problems classified as NP-hard can be solved efficiently for all members of their class, numerous algorithms have been devised that either strive for an approximate solution (e.g., simulated annealing or stochastic local searches) or solve the problems exactly in exponential time, but through mathematical insight (e.g., branch-and-bound or branch-and-cut), the algorithms increase the practically feasible sizes over what can be achieved using brute force. Evaluation of such algorithms often requires access to benchmarking problems of various types; ideally their difficulty should also be a controllable (or tunable) property. Ideas of generating test instances with planted solutions, that is, whose optimizing values are known to the problem constructor, have been explored in various fields for decades Bach (1983); Pilcher and Rardin (1987, 1992). A method of planting solutions to hard random Boolean satisfiability (SAT) problems based on statistical mechanics was first proposed in Ref. Barthel et al. (2002). In contrast, doing so for topologically structured problems is considerably less charted territory, as the correspondence between random SAT and the diluted spin glass disappears; thus replica symmetry breaking analysis Monasson and Zecchina (1996) no longer applies. As such, generating hard problems with a known solution for nearest-neighbor Ising models on a hypercube is a relatively new field of study.

While the task of creating such problems is of theoretical interest due to their potential assistance in answering open questions about the nature of spin glasses, the pace of research has been hastened in recent years by practical motivations, in particular, the availability of quantum annealing Johnson et al. (2011); Kadowaki and Nishimori (1998) and related analog devices (e.g., optical Wang et al. (2013)) physically implementing Ising Hamiltonian minimization heuristics. Such Hamiltonians are usually constrained by device manufacturing considerations to having short-ranged interaction terms Choi (2008); while more general objectives can typically be encoded onto their native graphs, this may require many system variables for each objective function variable, limiting the problem sizes that can be studied. Hence, a way of encoding topology-native problems is desirable.

In this work, we present a methodology for generating short-range Ising problems with a known ground state. Our primary focus is on the three-dimensional lattice with periodic boundary conditions, primarily because the three-dimensional spin glass is a prototypical complex system and of tremendous interest to condensed matter physicists, but also because the regular structure over the lattice allows a description of our key contributions in a natural and transparent manner. We stress, however, that the method is adaptable to different graph topologies.

By carefully decomposing the underlying graphical structure and selecting Ising subproblems over the resultant components, whose spin-spin interactions are restricted to being bimodal, from classes designed to have specific features influencing hardness, we obtain a factor graph representation having the topology of a body-centered cubic lattice (see Fig. 1). Following a Voronoi tessellation of this lattice, we note that the resultant problem is equivalent to a specific type of three-dimensional tile matching puzzle, a generalization of the two-dimensional edge-matching puzzle Demaine and Demaine (2007). In their two-dimensional form, these puzzles are constraint-satisfaction problems requiring placement of tiles from a given set in permissible spatial locations so that patterns on the tile edges match those of their neighbors. A well-known example is Eternity II ete (), a to date unsolved problem for whose solution a $2 million prize was once offered; the prize deadline expired in 2010. Our construction yields polyhedral rather than square puzzle pieces; in particular the units are truncated octahedra whose facets are of various color patterns (see Fig. 6 below). By varying the pattern properties via manipulation of the underlying Ising subproblems, we quickly and efficiently induce a tremendous variation in problem difficulty observed over all examined heuristic algorithms; hundreds of problems at any hardness level within our attainable range can be obtained in seconds with modest computing resources. Hard problems within our class are observed to be many orders of magnitude more difficult than randomly generated spin-glass problems with bimodal disorder, i.e., where the spin-spin interactions are randomly drawn from . Despite our restriction to using problem parameters in , we can directly generate problems with a known solution that are vastly more difficult than the Gaussian spin glass, which is itself known to be considerably harder than the random bimodal class. We thus consider our construction technique as giving rise to as yet unencountered types of disordered system and discuss our results within the context of hardness phase transitions for combinatorial problems.

The paper is structured as follows. In Sec. II we present an overview of recently-introduced planting techniques, followed by an outline of the planting approach presented in this work in Sec. III. Sections IV and V discuss the relationship of the presented problems to constraint satisfaction problems. Section VI presents numerical experiments using different algorithms highlighting the tunability of typical computational complexity for the problems with planted solutions, followed by discussions, generalizations to other graphs, and concluding remarks.

## Ii Related Work

Our methodology addresses several shortcomings of recently proposed techniques for solution planting in short-range Ising models. In Ref. Hen et al. (2015), a construction method was presented in which frustrated loops with known solution were employed as subproblem units. While the idea is appealing, the resultant problems typically have coupler strengths spanning a large range, posing serious issues for physical devices with precision limitations. An ad hoc limited-precision variant was introduced in Ref. King et al. (2015). Both approaches, however, bear the more serious and systematic disadvantage of tending to generate sets of problems whose members are not of sufficient hardness; we contend that this is a consequence of adding subproblem couplers. It is not hard to see that when partial Hamiltonians having consistent ground states are added, couplers in common which in the individual minima are either both satisfied or both violated will increase in magnitude. In a ground state, however, satisfied couplers are more common than violated ones, hence for any given shared coupler, it is relatively likely that the resultant (added) coupler tends to more strongly prefer alignment with the shared ground state. When this occurs over many edges, the problem’s degrees of freedom are “softly” constrained. The effect is the introduction of “hints” facilitating an optimization algorithm’s task.

In Ref. Marshall et al. (2016), a method was presented for generating difficult instances with bimodal couplers by iteratively changing the couplers’ signs to maximize the time required by a given solver to reach the (hypothesized) ground state. This is a promising approach as it both avoids the coupler addition issue and makes modest precision demands; however, it generally does not allow the ground state to be known with certainty and hence is difficult to scale to larger systems. Furthermore, the approach requires a Monte Carlo-like bond moving algorithm that will likely not easily scale up to much larger problem sizes. The authors do discuss a variant enabling solution planting, but as it again relies on adding couplers, the precision issue returns and, they report, the gains in hardness are modest compared to the mode not controlling the ground state. While the authors provide a useful analysis of post hoc empirical correlates of hardness, for example the parallel tempering mixing time, first proposed in Ref. Yucesoy et al. (2013), the generation procedure is ultimately a local search technique with respect to solver computational time, and hence does not yield or use physical or algorithmic insight at the problem level about what tunes difficulty. In other words, one could not, merely by reference to the resultant Hamiltonian, predict whether the problem is easy or hard. Finally, to design problems based on the time to solution of a given solver, assuming this will carry over to other optimization techniques requires theoretical backing that is missing to date.

A recent paper Wang et al. (2017) introduced the method of patch planting, in which subgraphs with known solution are coupled to each other, satisfying all interactions and thus planting an overall ground state. An advantage is that low-precision and tunable problems can be readily created for a given application domain. The authors report, however, that the resultant problems are observed to be less difficult than random problems, simply because frustrated loops typically do not extend beyond the subgraph building blocks, making the ground state less frustrated than in the fully random case. This problem can be slightly alleviated by post processing mining of the data Wang et al. (2017).

## Iii Problem Construction

Consider a cubic lattice of linear size with periodic boundary conditions. Let the underlying graph be with and its vertex and edge sets, respectively. Define to be the vertices of whose coordinates consist of integers of the same parity, for example, or . For each , define to be the vertices of the cubic unit cell implied when is opposite to the vertex with coordinates , where the addition is modulo to account for the periodic boundary. Let be the collection of all such cell vertex sets. This construction is partially illustrated in Fig. 1, where each colored unit cell represents a subgraph induced by a vertex set as defined. It should be clear that first, each vertex in appears in two and only two unit cells and second, that the unit cell subgraphs do not share any edges. In graph terminology, the family partitions into edge-disjoint induced subgraphs.

We are concerned here with constructing Ising problems on the lattice with zero field, i.e., whose Hamiltonians are of the form

(1) |

Due to the disjointness of the cell edge sets, we can regroup the Hamiltonian into terms each dependent only on the couplings within a unit cell in ,

(2) |

where

(3) |

The terms are called the unit cell subproblem Hamiltonians. Specifying each subproblem’s Hamiltonian is sufficient to imply a full Hamiltonian over the lattice. A straightforward consequence of the relation between summation and minimization, exploited in Ref. Hen et al. (2015), is that if the subproblems share the same minimizing configuration , then is also minimized at . This in turn suggests a natural construction procedure for a problem with a known ground state. We note a seemingly small but in fact deep and far-reaching difference between prior methodologies and what we propose is that subproblem couplers are never added, avoiding the issue discussed in Sec. II.

It may appear surprising at first that any interesting behavior can result from linking such apparently simple unit-cell subproblems, but this impression turns out to be quite false. Indeed, by restricting our attention to simple classes of subproblems whose couplers belong to the set , in other words, representable with a single bit of precision, we can effect dramatic changes in problem complexity, from trivial to many orders of magnitude harder than, e.g., Gaussian spin glasses.

Without loss of generality, we focus on planting the ferromagnetic ground state, i.e., and its image, because once such a problem has been generated, the structure of the construction procedure can be concealed by gauge randomization from a would-be adversary seeking the solution. Specifically, to translate to an arbitrarily chosen ground state , one would transform the initially determined couplings via

(4) |

We thus turn our attention to the set of unit-cell Hamiltonians having couplers of magnitude and ground state . Here we disregard the trivial subproblem comprised uniquely of ferromagnetic bonds; those remaining can be naturally partitioned into three types according to their number of frustrated facets. We call these types , , and , as they contain unit cells with two, four, and six frustrated facets, respectively. These types can be further partitioned into members equivalent under action of the cube graph automorphism group. More precisely, problems within a given equivalence class can be arrived at from one another via transformation from , the octahedral symmetries comprised of rotation and reflection leaving a (generic) cube invariant. It turns out that and contain two such classes (the orbits of acting on and , called and ) while all members of are equivalent under . Figure 2 illustrates an arbitrarily-chosen set of equivalence class representatives for the three problems types, while in Figure 3, a few members of type are shown.

While the subproblems can all readily be verified to have as a ground state, they clearly differ in the number of other minimizers they possess. Up to symmetry, members of have a unique ground state while those of have four. All problems in and have two and eight global minima, respectively. Figure 4 shows the four ground states of a specific member of . As we will see, this variation in subproblem ground-state degeneracy turns out to exert a crucial effect on the computational difficulty of the final planted problem.

We are now ready to straightforwardly specify our problem construction procedure. First, a probability distribution is chosen over the problem classes . A naive example, which in fact turns out to be a poor choice for generating hard problems, is to sample them uniformly. Next, a distribution over members of the classes is specified. In our examples this was always uniform. The randomness in generating the problems is thus decomposed into what we call subproblem disorder (over the classes) and disorder (over the members), the latter so called because it arises equivalently by sampling symmetries of the octahedral group. Subproblems over the unit-cell set are then sampled from the defined distributions.

Note that disorder is sufficient to induce considerable richness in the types of final problems, even when the cells in are all assigned from the same subproblem class. For example, it may appear at first that restricting all subproblems to the class (a regime which turns out to be highly interesting) will result in a so-called fully frustrated model Villain (1977), i.e., where all chordless cycles (or plaquettes) are frustrated in the final Hamiltonian, because problems in are themselves fully frustrated. That, however, would be incorrect; it is easy to check that cycles forming the facets of unit cubes not in may or may not be frustrated depending on the random orientation of the bounding -cell problems.

It should be pointed out that while the method does ensure that the planted state is in fact a ground state, as the couplers are selected from the class, it is in principle possible for other ground states to exist com (a). Our experimental validation in Sec. VI is thus concerned with the difficulty of finding any ground state. Methods of removing the degeneracy are certainly conceivable and are discussed briefly in Sec. VII.

We proceed next to a further exploration of the problem structure, discussing its connection to tiling puzzles.

## Iv The puzzle of spin glasses

The proposed problem generation formalism enables translation to a specific type of geometrically appealing constraint-satisfaction problem (CSP) called a tiling puzzle. While these problems are usually defined over two-dimensional domains, we will see that our three-dimensional technique yields a polyhedral generalization of an edge-matching puzzle, an example of which is Eternity II ete (). The idea is to abstract the Ising problem into a generic CSP representation known as a factor graph Kschischang et al. (2001). While the Ising Hamiltonian itself can always be interpreted as a factor graph, under our planting technique additional structure is present, namely, the restriction to certain subsets of states. Following construction of the CSP exploiting this structure, the natural connection to a tiling puzzle is observed.

Figure 5 illustrates the necessary steps for the mapping. For clarity, a direct two-dimensional analog of the three-dimensional case is shown, in which the unit cells induced by the decomposition into are the cycles forming a checkerboard pattern rather than unit cubes as displayed in Fig. 1.

We first consider the construction of a graph with respect to , which we call the subproblem lattice, denoted by . Suppose a vertex is placed at the centroid of each unit cell. The subproblem lattice is defined by considering two such vertices to be adjacent if and only if their corresponding unit cells share a corner vertex. Clearly, this lattice may be -colored, that is, partitioned into two subsets, none of whose members are neighbors. This in turn implies a partitioning of the cells in into two disjoint sets and , whose elements are labeled with colors and , respectively. Figure 5(a) demonstrates the two-dimensional construction of and its vertex coloring, where yellow (light shading) and blue (dark shading) may represent labels and or vice versa. In three space dimensions, each yellow (light shading) cell is surrounded by eight blue (dark shading) cells instead of four.

We next construct the factor graph CSP representation equivalent to the problem of minimizing . The idea is to derive an objective consisting of independent variables subject to equality constraints along the edges of . First, each vertex in , mapping to cell in the original lattice, is identified with a variable or , over eight Ising variables (four in two space dimensions) depending on its label in the -coloring. Because each vertex of the initial lattice is a shared corner of two unit cells in of opposite color, in the concatenation of the states variables and will occur exactly once for each vertex . The domain of each subproblem lattice variable is simply the ensemble of its subproblem ground states, the ground-state set , so that the full configuration space of the new problem is the Cartesian product

(5) |

To be a feasible solution to the original problem however, equality of neighboring variables must be enforced, that is, . If we define the problem

(6) |

where if and zero otherwise, we readily see that the solution is obtained for the overall planted ground state . To make the graphical connection, we note that each function in Eq. (6) can be interpreted as a factor over neighboring subproblem lattice variables (identified with) , i.e., with . The final maximization objective is then the product of all factors over subproblem lattice variables

(7) |

where the product index refers to edges of by the cells in mapping to their endpoint vertices. The factor graph corresponding to the two-dimensional lattice in Fig. 5 (a) appears in Fig. 5 (b), with variables represented by circular vertices (whose colors identify their labels as or ) and factors by the squares along the edges. In three dimensions, each variable is involved in eight factors rather than four.

Consider now the subproblem lattice with its vertices lying at their “natural” points in Euclidean space, i.e., at the centroids of the unit cells, rather than as a generic graph. A Voronoi tessellation Torquato (2013) is a partitioning of the space surrounding the vertices into convex polytopes, usually called tiles, but for reasons that will soon be clear, we refer to as tiling locations such that all points within a polytope are closest (typically in the -norm sense) to the enclosed vertex. It is easy to see in Fig. 5(a) that in two dimensions, the subproblem lattice forms a periodic pattern, sometimes called the quincunx, of squares of length , each with a central vertex. The appropriate tessellation of is a partitioning into tilted square regions centered at each vertex. Neglecting for an instant the significance of the colors, Fig. 5(c) shows five full and eight partial tiling squares in the section of lattice drawn. Returning to the CSP in which the task was to assign the variables from their respective domains (ground-state sets) so that factor constraints are satisfied, we now observe the equivalence to a tile-matching puzzle. This is comprised of a board (the space embedding ), segmented into octagons (the cells in the Voronoi tessellation), each of which is endowed with a repertoire of tiles (the ground-state set associated with the location). The tile faces could be in one of two colors (the spin values, not to be confused with the -coloring introduced to construct the CSP). The puzzle task is to select the tiles so that no two adjacent tile faces have a different color. If orange (light) and purple (dark) represent spin values, say, and , respectively, then Fig. 5 (c) shows a hypothetical tiling assignment. Here a violation occurs where the tile colors on both sides of a factor differ.

In three dimensions, the puzzle is analogous, but with a few twists. In a generalization of the two-dimensional case, the subproblem lattice will, when respecting Euclidean coordinates, form a pattern of cubes of length , each with a vertex at its centroid. Such a formation is known as the body-centered-cubic lattice. The corresponding Voronoi tessellation is the space-filling bitruncated cubic honeycomb whose base unit is the truncated octahedron Conway et al. (2016); Torquato (2013). Each of these polyhedral cells, which now define the three-dimensional tiling locations, consists of eight hexagonal facets, through which the the edges of our CSP factor graph pass, and six square facets, defining the boundaries between vertices at distance , but through which no CSP factors pass. The ground-state sets associated with each vertex of thus map to coloring configurations on the hexagonal facets of the enclosing Voronoi cell, defining a set of tiles, and again the task is to place the tiles while meeting all coloring constraints. The construction is illustrated in Fig. 6.

While technically fitting the definition of a tile-matching puzzle, the specific examples proposed here differ in a key way from conventional puzzles such as Eternity II. In the latter, each tile from the given set may be placed at more than one board location (in an allowable orientation) while presently the ground-state set associated with each site constrains the choices. This fact, along with the polynomial-time solvability Barahona (1982) of the original ground-state problem in two dimensions, could lead one to suspect that the puzzles we have contrived may not be overly interesting. In three space dimensions, however, solution by graph matching is no longer applicable, and as we will see, the problems exhibit a rich range of behavior intimately tied to the known properties of generic constraint satisfaction problems.

## V The phases of satisfaction

Hardness phase transitions in combinatorial problems have been an active area of study for several decades now Hogg et al. (1996). Of particular interest is the emergence of difficulty divergencies in random SAT problems as parameters guiding instance generation are varied Mézard et al. (2002). While the tiling CSPs we have proposed here can certainly be expressed as SAT problems com (b), the resultant formulas do not, however, follow the typical assumptions of random SAT. Most notably, the variables are not free to appear in any clause with some probability, but follow the highly structured topological constraints of the puzzle. Of course, the subproblem and disorder impose stochasticity over the factor graph variable domains , but this is a localized randomness, violating key assumptions in the analysis of random SAT. Nonetheless, we expect that commonly known facts and intuitions regarding the relation between how constrained a problem is and its empirical hardness will apply to our situation.

Because the tiling puzzle color factors are invariant to the distribution over cell problems, constrainedness depends exclusively on the features of the ground-state sets. As we alluded to earlier, a dominant correlate of problem hardness is the average size of the ground-state sets . A highly constrained regime is one in which the subproblems tend to have relatively small degeneracy and vice versa for a loosely constrained one. A subproblem demonstrating a maximal constraint level is any member of class (with a unique ground state up to ), while one constrained minimally, which was not considered in the experiments, is the zero subproblem, in which all states are ground states. Within the class of subproblems introduced in Sec. III, including members of and (with degeneracy ) tends to bias the final problem towards the maximal and minimal extremes, respectively.

By construction, all our problems are satisfiable, i.e., there is no unsatisfiable (or overconstrained) regime as usually understood. Clearly, the minimally constrained full problem is trivial as each state is a solution to the original problem. In the tiling puzzle representation, the factor constraints can be met independently because the resultant states are always in the ground-state sets. At the other extreme, highly constrained problems are also typically easy: Because there are few local choices from the ground-state sets, a tree search method will encounter relatively small branching factors when traversing the state space, while for local search, the energy landscape will be such that the algorithm can reliably locate the ground state using simple strategies to escape from local minima. These extremes strongly suggest that for some intermediate level of mean ground-state set size, a peak in problem hardness will occur. While it is not possible to increase past for unit cells without introducing more complex Ising interactions (or without trivially removing all ), an interesting question, probed experimentally in the following, is, given the proposed classes of unit-cell subproblems (, , and ) and all possible mixture distributions over the classes, did the hypothesized peak occur for some interior point of the distribution set? An affirmative answer would suggest that the corresponding CSP class is among the hardest of all tile-matching puzzles of the type we have introduced, i.e., for which the locations’ tile sets may be chosen arbitrarily, not by constraining them to map to subproblem ground states. Conversely, if the hardness maximum within the set of mixtures occurs at the set boundary, it would point to the existence of more difficult tiling puzzles not representable within our framework.

## Vi Experiments

We now numerically study the typical complexity of certain subproblem classes to illustrate their tunability. This section focuses on three-dimensional cubic lattices, followed by a demonstration that the approach works well also in two space dimensions (Sec. VII.2). We conclude with a discussion of generalizations to other graph types.

While we would ideally have performed a rigorous numerical study for a discrete “grid” representing all subproblem types we have proposed, this is computationally prohibitive. Consequently, we focus on two illustrative regimes demonstrating interesting variations in problem difficulty over three algorithms, and the role played by subproblem degeneracy. The simulation results show that highly complex problems with known ground state, and more difficult than conventional spin glasses with both Gaussian and random bimodal couplings, are accessible via our methodology, but they also strongly suggest that the hardness maximum over the distribution set occurs at a boundary, namely, when the distribution is a point mass concentrated on class . Perhaps unsurprisingly then, we conclude that tiling puzzles equivalent to Ising problems constructed using the classes are likely not the hardest among those where the tile sets are arbitrarily specifiable.

All problems are defined on three-dimensional lattices of size with periodic boundary conditions. We consider two classes of problems corresponding to one-dimensional slices of the problem mixture parameter space. The first class, gallus_26, allows subproblems to belong solely to classes and , while in the second, gallus_46, they are constrained to and . Both instance classes are parametrized by , the probability of selecting class instead of the alternative com (c). The classes each contained instances generated at uniformly spaced values of , for a total of instances per class. All subproblems are subject to uniform disorder. The selected range of values is interesting as within it, problems are more difficult than conventional spin glasses with Gaussian or random bimodal couplings. We note that the problems continue to become easier as decreases. Using the values of shown in Fig. 2, the expected ground-state set sizes in terms of are for gallus_26 and for gallus_46.

Problem difficulty is assessed through performance of three different algorithms designed for problems with rough energy landscapes. The measures of difficulty are in strong agreement across the methods, providing corroboration that the observed difficulty trends should persist robustly across a fairly large class of heuristic algorithms, including backtrack-based search Hogg et al. (1996).

Simulated annealing (SA) Kirkpatrick et al. (1983) is the most basic algorithm considered. We use the optimized implementation developed by Isakov et al. Isakov et al. (2015) with and . While we would rather have used some form of optimized time to solution (TTS) measure, which considers the best tradeoff between simulation length and number of simulations, SA runs were too time consuming on the hard problems to generate the requisite runtime histograms reliably. Consequently, we selected a fixed run length of sweeps for all problems with a single sweep per temperature. Each problem is simulated times independently. The SA TTS is defined as the computational time required to find a ground state with at least % probability, i.e., for each instance, , where is the fraction of successful runs out of the repetitions.

Furthermore, we have used a highly optimized implementation of parallel tempering Monte Carlo with isoenergetic cluster moves (ICMs) Zhu et al. (2015), an adaptive hybrid parallel tempering (PT) Hukushima and Nemoto (1996); Geyer (1991); Moreno et al. (2003) cluster algorithm Houdayer (2001). Because the ICM is considerably more efficient than SA, it allowed us to gather runtime statistics for each instance, allowing optimized s to be computed. In contrast to SA, total ICM computational time is not merely a function of overall Monte Carlo sweeps, but includes the additional cost of constructing random-sized clusters. Consequently, to track aggregated ICM computational effort we record run times in seconds on hardware dedicated entirely to the simulations. If is the empirically observed probability of finding the ground state in time or less, the optimized time to solution is defined as . We note, however, that despite the efficiency of the ICM, it commonly failed to find the solution for the harder problems within the maximum allowed total Monte Carlo sweeps, requiring min of real time, when computing the runtime histograms. For difficult classes, these “timeouts” occurred at least half of the time within the ICM repetitions used to tally the histograms. Fortunately, we were nonetheless able to infer an optimized from the conditional histogram in which solutions were found. We used temperatures spaced within and ; see Ref. Zhu et al. (2015) for further details of the ICM.

The final set of tests have used the sequential Monte Carlo Del Moral et al. (2006) method known as population annealing Monte Carlo (PAMC) Hukushima and Iba (2003); Machta (2010). This algorithm is related to SA but differs crucially in its usage of weight-based resampling, which multiplies or eliminates members of a population according to the ratios of their Boltzmann factors at adjacent temperatures, maintaining thermal equilibrium at each step. Our simulations used temperatures with and sweeps per temperature, and a population size of replicas. In Ref. Wang et al. (2015a), a PAMC-derived index of landscape roughness called the entropic family size was proposed. If is the fraction of replicas in the final population descended from initial replica , then , where . An energy landscape is thus deemed rough if is small, that is, if relatively few initial replicas survive to the final distribution, yielding a large value of . Note that converges quickly in and is easily estimated with finite populations. Thermal equilibration is ensured by requiring ; when this is not satisfied for an instance, it is rerun with a larger population size. The entropic family size is known to covary strongly with the PT autocorrelation time Wang et al. (2015a) and, as can be seen from Fig. 7, does so for the -based hardness metrics considered here as well. Figure 7 shows, respectively, the dependence of both SA and ICM measures on PAMC for the subclasses of gallus_26 [Fig. 7(a)] and gallus_46 [Fig. 7(b)] studied, where the averages at each point representing a subclass are computed over its generated instances. The plots clearly show a distinct near-linear dependence of both time-based measures on the PAMC-defined value. A larger on average implies a longer time to solution with respect to both SA and ICM algorithms, corroborating the former’s power as an objective measure of landscape roughness.

Results of the PAMC simulations are shown in Fig. 8. We observe a trend of increasing (and hence complexity) for both problem classes as the fraction of subproblems from increases towards unity. For comparison, we display the mean value of two prototypical problems with rough landscapes on the same lattice, the random and Gaussian [ a normal distribution with zero mean and variance one] spin glass, computed using and instances of each type, respectively. It is clear that problems in most of the examined subclasses of gallus_26 and gallus_46 are more difficult than both of these widely studied problem classes. Indeed, for subclasses of gallus_26 corresponding to , problems are orders of magnitude harder than bimodal spin glasses. Perhaps more surprisingly, they are also orders more difficult than Gaussian spin glasses, despite the latter possessing continuous-valued couplings while our instances restrict couplers to , which are believed to generally be easier to minimize.

Setting yields the most complex problems of those considered. In fact, we conjecture, based on less comprehensive simulations, that these instances are the hardest among all those constructed with subproblem classes in . We plan to perform a comprehensive analysis in the future. For this hard class, the unit-cell subproblems, deriving exclusively from , have eight ground states each. Because this hardness peak occurs at the boundary of the problem parameter space, it seems plausible that one could instantiate yet more complex three-dimensional tiling puzzles, where the locations were allowed more than eight (times two) tile possibilities but not so many that the problem becomes underconstrained and easy to solve.

At first, the greater pointwise difficulty shown in Fig. 8 of an mixture over one made of appears to contradict intuition that greater frustration implies higher difficulty. After all, with four bounding frustrated facets, a member of can be interpreted as more frustrated than one of , which has two. The story is somewhat more subtle though: While frustration certainly plays a role in tuning hardness in our problems, it appears to do so through its effect on constraint level, namely, on the sizes of , with inducing higher complexity than because its ground-state set size is larger. This fact is displayed in Fig. 8(b), where for the two classes is plotted against instead of . The graph shows that for a given value of , gallus_26 and gallus_46 have a very similar roughness index, implying that mean subproblem degeneracy accurately predicts difficulty regardless of the underlying subproblem mixture.

For completeness, analogous plots displaying similar difficulty trends for SAs and ICMs are shown in Fig. 9, where they are again plotted against .

In contrast to our results so far, which have considered instance-averaged difficulty measures, Fig. 10 shows histograms of optimized ICM values for three subclasses of gallus_26 [Fig. 10(a)] and of gallus_46 [Fig. 10(b)] indexed by . The leftward shift in histogram support shows clearly that problems tend to become easier with decreasing . Given the histogram shapes, one may naturally suspect that is log-normally distributed. As is unlikely to be precisely normal across the entire data range, we visualize correspondence with a Gaussian via normal probability plots Chambers et al. (1983), which relate the sample order statistics, obtained by sorting the data, with the theoretical means of the corresponding normal order statistics. Deviations from a linear relation signal lack of Gaussianity. Figure 11 shows the probability plots for the three subclasses used to generate the preceding histograms for the gallus_26 [Fig. 11(a)] and gallus_46 [Fig. 11(b)] classes, respectively. The relation is close to linear over the majority of the histogram support, implying in turn that the is approximately log-normally distributed. This is consistent with quantities such as the parallel tempering Monte Carlo autocorrelation time and other roughness measures Katzgraber et al. (2006); Yucesoy et al. (2013) having the same property.

While we have argued that is a good predictor of mean difficulty, the histogram results show that this value by no means provides a complete picture. Indeed, when (i.e., ) there is no variance in as all subproblems have eight-fold degeneracy, yet the distributions in Fig. 10 nonetheless show considerable intraclass spread in difficulty. Therefore, there are certainly other factors at play in predicating the hardness.

## Vii Discussion

In this section we discuss generalizations of the planting approach using lattice animals. Furthermore, we present a case study on how the approach generalizes to other non hypercubic lattices. Finally, we discuss the use of the planted problems for fundamental studies of spin glasses and related statistical systems.

### vii.1 Generalization via lattice animals

The proposed unit-cell planting technique shows encouraging results and properties and one may inquire how it may be generalized. In this section, we present a natural extension of the idea, still assuming lattice-structured problems, where instead of defining subproblems on the unit cells of , they are specified on subgraphs consisting of their unions. One can verify that such subgraphs, called lattice animals or polyominoes, also partition the lattice into edge-disjoint subgraphs, meaning that subproblem couplers are still not added.

A two-dimensional example of decomposition into lattice animals is shown in Fig. 12. Generalization to three-dimensional polyominoes obtained by grouping cells from the decomposition shown in Fig. 1 is straightforward. Shown in non-gray colors are six lattice animals. Only the pink (top right) one is made of a single cell. The key difference from the basic method is that unit cells of a given color are no longer constrained to have their individual ground states agree; only the complete animal ground state is relevant. This extension thus considerably extends the types of subproblems that can be employed and also introduces an additional mechanism of solution hiding, namely, via randomization of the employed lattice animals, which would of course be unknown to the would-be adversary. While the tiling puzzle and CSP interpretation of the problem, suitably modified, continues to hold in this generalization, under lattice animal randomization, the adversary would in essence no longer be certain what puzzle they are even solving.

In this work, we have considered carefully chosen families of subproblems on three-dimensional unit cells. An exploration of extensions to more general lattice animals is outside the scope of the present work. We note, however, that if subproblem couplers are sampled from a given distribution, then provided the subgraph tree widths are small, their ground states can be computed exactly Koller and Friedman (2009) and gauged to the desired overall ground state.

Finally, we note that this lattice animal generalization is also key when attempting to reduce degeneracy in the planted problems. For example, by selecting the coupler values from a Sidon set Sidon (1932); Katzgraber et al. (2015); Karimi et al. (2017) of the form with, e.g., , the degeneracy is drastically reduced. Similarly, one could select the couplers from a distribution of the form

(8) |

where is a uniform random number and close to Katzgraber (2010). However, instead of using basic tiles, more complex lattice animals must be used to construct the problems.

### vii.2 Generalization to arbitrary graphs

The need to benchmark both classical and quantum optimization heuristics has hastened the development of advanced planting techniques for solutions of spin-glass-like optimization problems. We have so far focused our attention on hypercubic lattice problems. However, quantum annealing machines typically have hardwired hardware graphs that are close to planar. A prototypical example is the chimera graph Bunyk et al. (2014) used by current versions of the D-wave quantum annealing machines. When presented with such a situation, one has two possible options to apply our planting framework. The first is to impose a lattice structure onto the available graph and the second is to specify subproblems on altogether different “unit cells” than cubes, which of course, must be graph dependent. We now discuss the first option applied to the special case of chimera.

Although in principle three-dimensional lattices can be implemented on the chimera graph Harris et al. (2017), the required overhead may limit linear (planted) problem sizes that can be practically studied. On the other hand, relatively large two-dimensional lattices with periodic boundaries can be straightforwardly embedded on chimera, with a relatively modest constant ratio of two chimera variables per lattice spin. More precisely, a chimera graph of bipartite unit cells, comprised of variables in total, can produce a toroidal two-dimensional lattice of size . Rather than formally describe the rather natural procedure, we illustrate it in Fig. 13, where a periodic lattice is created from a cell chimera graph.

So far, we have focused primarily on planting subproblems on three-dimensional unit cells, in part because planar problems without a field are solvable in polynomial time Barahona (1982), but as illustrated in Figs. 5 and 12, the two-dimensional analog is readily obtained. The analytical tractability of the planar lattice enables a deep exploration of physical and computational properties, a work which will be reported elsewhere Perera et al. (2018). Here we outline the idea and demonstrate that it does indeed perform well on planar lattices and nonplanar quasi-two-dimensional chimera graphs using population annealing simulations. The two-dimensional subproblems, i.e., the analogs of the cells shown in Fig. 2, are partitioned into classes for , within which cells have minimizing ground-state configurations each. To achieve the construction, first define two magnitudes with ; presently we take and . A cell in class is constructed by setting a random edge to be antiferromagnetic with magnitude , of the remaining edges to be ferromagnetic with magnitude , and the leftover edges to be ferromagnetic with magnitude com (d). It is readily verified that the subproblems do indeed have the specified number of local ground states which always include the ferromagnetic state.

As was done in three space dimensions, we consider two classes of problems each comprised of mixtures of two subproblem classes. Problem class C_24, consists of mixtures of and cells, while in C_23 the cells may belong to or . Both problem classes are parametrized by , the probability of choosing a cell from . The results on periodic lattices are shown in Fig. 14. Again, a wide range of landscape roughness is seen to be attainable, including a regime in which the problems are more difficult than random bimodal ( couplers) spin glasses. An interesting distinction from three-dimensional results is that the most difficult problems are not those in which cells exclusively belong to the maximally degenerate class , but rather the moderately constrained class. In fact, class is a highly underconstrained regime for this topology and gives rise to very easy problems. The fascinating connections between dimension, complexity, and phase behavior are the subject of ongoing study that goes beyond the scope of this paper.

If one wished to move beyond a regular lattice structure, the natural objective is a decomposition of the problem graph into edge-disjoint subgraphs satisfying some constraint allowing tractable minimization. An example of a constrained decomposition is into subgraphs with given minimum vertex degree Yuster (2013). More directly applicable to the planting context would be a constraint on maximum subgraph tree width thereby allowing determination of the planted subproblem ground states. The edge-disjointness property is the key common aspect with the ideas presented in this paper, as it continues to circumvent the need for adding subproblem couplers. This idea is clearly a generalization of the lattice animal methodology discussed previously. While we have not presently considered planting using such generic subgraphs, one can readily envision a heuristic decomposition algorithm that greedily grows partitioning subgraphs until their tree widths exceed some criterion.

### vii.3 Spin-glass physics

We hope researchers in the field embrace these planted problems to study physical properties of glassy systems beyond the benchmarking of optimization heuristics.

Having arbitrarily large planted solutions for hypercubic lattices allows one to address different problems in the physics of spin glasses. For example, the computation of defect energies, intimately related to fundamental properties of these paradigmatic disordered systems, strongly depends on the knowledge of ground states (see, for example, Refs. Hartmann (1997, 1999); Palassini and Young (1999, 2001); Hartmann (2001); Katzgraber et al. (2001)). Being able to plant problems would drastically reduce the computational effort in answering these fundamental questions.

Furthermore, by carefully tuning the different instance classes, problems with different disorder and frustration can be generated. A systematic study of the interplay between disorder and frustration is therefore possible for nontrivial lattices beyond hierarchical ones. Similarly, being able to tune the fraction of frustrated plaquettes allows one to carefully study the emergence of chaotic effects in spin glasses (see, for example, Refs. Katzgraber and Krzakala (2007); Thomas et al. (2011); Wang et al. (2015b) and references therein).

### vii.4 Application-based benchmark problems

It is well established that random benchmark problems Rønnow et al. (2014); Katzgraber et al. (2014) for classical and quantum solvers using spin-glass-like problems are computationally typically easy. Furthermore, the control over the hardness of the benchmark problems has been rather limited either because (post) processing is expensive Katzgraber et al. (2015); Marshall et al. (2016) or because the benchmark generation approach does not give the user enough control over the problems to match, e.g., hardware restrictions Hen et al. (2015).

Because application-based problems from industrial settings are highly structured, they pose additional challenges for the vast pool of optimization techniques designed, in general, for random unstructured problems. This has sparked the use of problems from industry to generate hard (and sometimes tunable) benchmark problems. Most notably, the use of circuit fault diagnosis Perdomo-Ortiz et al. (2015); Perdomo-Ortiz et al. () has produced superbly hard benchmarks with small number of variables. However, the use of application-based problems for benchmarking is in its infancy and while circuit fault diagnosis shows most promise Perdomo-Ortiz et al. (), many applications have produced benchmark problems that lack the richness needed to perform systematic studies; see, for example, Refs. Perdomo-Ortiz et al. (2012); Santra et al. (2014); Rieffel et al. (2015); Azinovic et al. (2017).

The problems presented in this work are highly tunable and computationally easy to generate. Furthermore, they can be embedded in more complex graphs, as is, for example, commonly done on the D-wave hardware for application benchmarks. Thus, having this tunability not only should allow for the generation of problems that might elucidate quantum speedup in analog annealers, but might also help gain a deeper understanding into quantum annealing for spin glasses in general. In parallel, having these tunable problems might elucidate the application scope of specific optimization techniques, both classical and quantum.

## Viii Conclusions

We have presented an approach for generating Ising Hamiltonians with planted ground-state solutions and a tunable complexity based on a decomposition of the model graph into edge-disjoint subgraphs. Although we have performed the construction for three-dimensional cubic lattices and illustrated the approach with the two-dimensional pendant, the approach can be generalized to other lattice structures, as shown in Sec. VII.2 for the chimera lattice. The construction allows for a wide range in computational complexity depending on the mix of the elementary building blocks used. We corroborated these results with experiments using different optimization heuristics. Subsequent studies should discuss constructions with controllable ground-state degeneracy, as well as the mapping of the complete complexity phase space.

## Acknowledgements

We acknowledge helpful discussions with Jon Machta, Bill Macready, Catherine McGeogh, Humberto Munoz-Bauza, and Martin Weigel and are grateful to Fiona Hanington for thorough suggestions to the manuscript. F. H. would like to thank the Santa Fe Institute for its hospitality and invitation to a stimulating Working Group on solution planting, and is indebted to those who introduced him to the original Eternity puzzle. H. G. K. acknowledges support from the NSF (Grant No. DMR-1151387). The Texas A&M team’s research is based upon work supported by the Office of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA), via Interagency Umbrella Agreement No. IA1-1198. The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of the ODNI, IARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon. We thank Texas A&M University and the Texas Advanced Computing Center at University of Texas at Austin for providing HPC resources.

## References

- Bach (1983) E. Bach, in Proceedings of the fifteenth annual ACM Symposium on Theory of Computing (ACM, 1983).
- Pilcher and Rardin (1987) M. G. Pilcher and R. L. Rardin, A random cut generator for symmetric traveling salesman problems with known optimal solutions, Institute for Interdisciplinary Engineering Studies, Purdue University Lafayette, IN, Tech. Rep. CC-87-4 (1987).
- Pilcher and Rardin (1992) M. G. Pilcher and R. L. Rardin, Partial polyhedral description and generation of discrete optimization problems with known optima, Naval Research Logistics (NRL) 39, 839 (1992).
- Barthel et al. (2002) W. Barthel, A. K. Hartmann, M. Leone, F. Ricci-Tersenghi, M. Weigt, and R. Zecchina, Hiding solutions in random satisfiability problems: A statistical mechanics approach, Phys. Rev. Lett. 88, 188701 (2002).
- Monasson and Zecchina (1996) R. Monasson and R. Zecchina, Entropy of the K-satisfiability problem, Phys. Rev. Lett. 76, 3881 (1996).
- Johnson et al. (2011) M. W. Johnson, M. H. S. Amin, S. Gildert, T. Lanting, F. Hamze, N. Dickson, R. Harris, A. J. Berkley, J. Johansson, P. Bunyk, et al., Quantum annealing with manufactured spins, Nature 473, 194 (2011).
- Kadowaki and Nishimori (1998) T. Kadowaki and H. Nishimori, Quantum annealing in the transverse Ising model, Phys. Rev. E 58, 5355 (1998).
- Wang et al. (2013) Z. Wang, A. Marandi, K. Wen, R. L. Byer, and Y. Yamamoto, Coherent Ising machine based on degenerate optical parametric oscillators, Phys. Rev. A 88, 063853 (2013).
- Choi (2008) V. Choi, Minor-embedding in adiabatic quantum computation. I: The parameter setting problem., Quantum Inf. Process. 7, 193 (2008).
- Demaine and Demaine (2007) E. D. Demaine and M. L. Demaine, Jigsaw puzzles, edge matching, and polyomino packing: Connections and complexity, Graphs and Combinatorics 23, 195 (2007).
- (11) Eternity II puzzle, accessed 19 May 2017, URL https://goo.gl/n3VvsX.
- Hen et al. (2015) I. Hen, J. Job, T. Albash, T. F. Rønnow, M. Troyer, and D. A. Lidar, Probing for quantum speedup in spin-glass problems with planted solutions, Phys. Rev. A 92, 042325 (2015).
- King et al. (2015) A. D. King, T. Lanting, and R. Harris, Performance of a quantum annealer on range-limited constraint satisfaction problems (2015), arXiv:1502.02098.
- Marshall et al. (2016) J. Marshall, V. Martin-Mayor, and I. Hen, Practical engineering of hard spin-glass instances, Phys. Rev. A 94, 012320 (2016).
- Yucesoy et al. (2013) B. Yucesoy, J. Machta, and H. G. Katzgraber, Correlations between the dynamics of parallel tempering and the free-energy landscape in spin glasses, Phys. Rev. E 87, 012104 (2013).
- Wang et al. (2017) W. Wang, S. Mandrà, and H. G. Katzgraber, Patch-planting spin-glass solution for benchmarking, Phys. Rev. E 96, 023312 (2017).
- Villain (1977) J. Villain, Spin glass with non-random interactions, J. Phys. C 10, 1717 (1977).
- com (a) Referring to the subproblem classes though, a few moments’ reflection reveals that interestingly, free spin degeneracy is impossible in the planted state.
- Kschischang et al. (2001) F. R. Kschischang, B. J. Frey, and H.-A. Loeliger, Factor graphs and the sum-product algorithm, IEEE Transactions on information theory 47, 498 (2001).
- Torquato (2013) S. Torquato, Random heterogeneous materials: microstructure and macroscopic properties (Springer Science & Business Media, Berlin, 2013).
- Conway et al. (2016) J. H. Conway, H. Burgiel, and C. Goodman-Strauss, The symmetries of things (CRC Press, London, 2016).
- Barahona (1982) F. Barahona, On the computational complexity of Ising spin glass models, J. Phys. A 15, 3241 (1982).
- Hogg et al. (1996) T. Hogg, B. A. Huberman, and C. P. Williams, Phase transitions and the search problem (Elsevier, Amsterdam, 1996).
- Mézard et al. (2002) M. Mézard, G. Parisi, and R. Zecchina, Analytic and Algorithmic Solution of Random Satisfiability Problems, Science 297, 812 (2002).
- com (b) If factor graph domain values and are mapped to Boolean values, say, and , respectively, the color constraints are clearly xor functions and their conjunction must be part of the Boolean formula. To obtain a conjunctive normal form (CNF) expression, a naive way to incorporate the ground sets is to straightforwardly represent each variable’s in disjunctive normal form, convert it to CNF, and include it in the formula by conjunction. The expression size is likely highly suboptimal but, for our case, nonetheless polynomially bounded and illustrates the point.
- com (c) Hence, when , the classes are identical.
- Kirkpatrick et al. (1983) S. Kirkpatrick, C. D. Gelatt, Jr., and M. P. Vecchi, Optimization by simulated annealing, Science 220, 671 (1983).
- Isakov et al. (2015) S. V. Isakov, I. N. Zintchenko, T. F. Rønnow, and M. Troyer, Optimized simulated annealing for Ising spin glasses, Comput. Phys. Commun. 192, 265 (2015), (see also ancillary material to arxiv:cond-mat/1401.1084).
- Zhu et al. (2015) Z. Zhu, A. J. Ochoa, and H. G. Katzgraber, Efficient Cluster Algorithm for Spin Glasses in Any Space Dimension, Phys. Rev. Lett. 115, 077201 (2015).
- Hukushima and Nemoto (1996) K. Hukushima and K. Nemoto, Exchange Monte Carlo method and application to spin glass simulations, J. Phys. Soc. Jpn. 65, 1604 (1996).
- Geyer (1991) C. Geyer, in 23rd Symposium on the Interface, edited by E. M. Keramidas (Interface Foundation, Fairfax Station, VA, 1991), p. 156.
- Moreno et al. (2003) J. J. Moreno, H. G. Katzgraber, and A. K. Hartmann, Finding low-temperature states with parallel tempering, simulated annealing and simple Monte Carlo, Int. J. Mod. Phys. C 14, 285 (2003).
- Houdayer (2001) J. Houdayer, A cluster Monte Carlo algorithm for 2-dimensional spin glasses, Eur. Phys. J. B. 22, 479 (2001).
- Del Moral et al. (2006) P. Del Moral, A. Doucet, and A. Jasra, Sequential monte carlo samplers, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68, 411 (2006).
- Hukushima and Iba (2003) K. Hukushima and Y. Iba, in The Monte Carlo method in the physical sciences: celebrating the 50th anniversary of the Metropolis algorithm, edited by J. E. Gubernatis (AIP, 2003), vol. 690, p. 200.
- Machta (2010) J. Machta, Population annealing with weighted averages: A Monte Carlo method for rough free-energy landscapes, Phys. Rev. E 82, 026704 (2010).
- Wang et al. (2015a) W. Wang, J. Machta, and H. G. Katzgraber, Population annealing: Theory and application in spin glasses, Phys. Rev. E 92, 063307 (2015a).
- Chambers et al. (1983) J. M. Chambers, W. S. Cleveland, B. Kleiner, and P. A. Tukey, Graphical methods for data analysis, vol. 5 (Wadsworth, Belmont, CA, 1983).
- Katzgraber et al. (2006) H. G. Katzgraber, S. Trebst, D. A. Huse, and M. Troyer, Feedback-optimized parallel tempering Monte Carlo, J. Stat. Mech. P03018 (2006).
- Koller and Friedman (2009) D. Koller and N. Friedman, Probabilistic graphical models: principles and techniques (MIT press, Cambridge, MA, 2009).
- Sidon (1932) S. Sidon, Ein Satz über trigonometrische Polynome und seine Anwendung in der Theorie der Fourier-Reihen, Mathematische Annalen 106, 536 (1932).
- Katzgraber et al. (2015) H. G. Katzgraber, F. Hamze, Z. Zhu, A. J. Ochoa, and H. Munoz-Bauza, Seeking Quantum Speedup Through Spin Glasses: The Good, the Bad, and the Ugly, Phys. Rev. X 5, 031026 (2015).
- Karimi et al. (2017) H. Karimi, G. Rosenberg, and H. G. Katzgraber, Effective optimization using sample persistence: A case study on quantum annealers and various Monte Carlo optimization methods, Phys. Rev. E 96, 043312 (2017).
- Katzgraber (2010) H. G. Katzgraber, Random Numbers in Scientific Computing: An Introduction (2010), (arXiv:1005.4117).
- Bunyk et al. (2014) P. Bunyk, E. Hoskinson, M. W. Johnson, E. Tolkacheva, F. Altomare, A. J. Berkley, R. Harris, J. P. Hilton, T. Lanting, and J. Whittaker, Architectural Considerations in the Design of a Superconducting Quantum Annealing Processor, IEEE Trans. Appl. Supercond. 24, 1 (2014).
- Harris et al. (2017) R. Harris et al., in preparation (2017).
- Perera et al. (2018) D. Perera, F. Hamze, and H. G. Katzgraber, in preparation (2018).
- com (d) Hence two bits of precision are required for these classes.
- Yuster (2013) R. Yuster, Edge-disjoint induced subgraphs with given minimum degree, The Electronic Journal of Combinatorics 20, P53 (2013).
- Hartmann (1997) A. Hartmann, Evidence for nontrivial ground-state structure of 3d spin glasses, Europhys. Lett. 40, 429 (1997).
- Hartmann (1999) A. K. Hartmann, Scaling of stiffness energy for three-dimensional Ising spin glasses, Phys. Rev. E 59, 84 (1999).
- Palassini and Young (1999) M. Palassini and A. P. Young, Triviality of the ground state structure in Ising spin glasses, Phys. Rev. Lett. 83, 5126 (1999).
- Palassini and Young (2001) M. Palassini and A. P. Young, The spin glass: Effects of ground state degeneracy, Phys. Rev. B 63, 140408(R) (2001).
- Hartmann (2001) A. K. Hartmann, Ground-state clusters of two-, three-, and four-dimensional Ising spin glasses, Phys. Rev. E 63 (2001).
- Katzgraber et al. (2001) H. G. Katzgraber, M. Palassini, and A. P. Young, Monte Carlo simulations of spin glasses at low temperatures, Phys. Rev. B 63, 184422 (2001).
- Katzgraber and Krzakala (2007) H. G. Katzgraber and F. Krzakala, Temperature and Disorder Chaos in Three-Dimensional Ising Spin Glasses, Phys. Rev. Lett. 98, 017201 (2007).
- Thomas et al. (2011) C. K. Thomas, D. A. Huse, and A. A. Middleton, Chaos and universality in two-dimensional Ising spin glasses, Phys. Rev. Lett. 107, 047203 (2011).
- Wang et al. (2015b) W. Wang, J. Machta, and H. G. Katzgraber, Chaos in spin glasses revealed through thermal boundary conditions, Phys. Rev. B 92, 094410 (2015b).
- Rønnow et al. (2014) T. F. Rønnow, Z. Wang, J. Job, S. Boixo, S. V. Isakov, D. Wecker, J. M. Martinis, D. A. Lidar, and M. Troyer, Defining and detecting quantum speedup, Science 345, 420 (2014).
- Katzgraber et al. (2014) H. G. Katzgraber, F. Hamze, and R. S. Andrist, Glassy Chimeras Could Be Blind to Quantum Speedup: Designing Better Benchmarks for Quantum Annealing Machines, Phys. Rev. X 4, 021008 (2014).
- Perdomo-Ortiz et al. (2015) A. Perdomo-Ortiz, B. O’Gorman, J. Fluegemann, R. Biswas, and V. N. Smelyanskiy, Determination and correction of persistent biases in quantum annealers (2015), (arXiv:quant-phys/1503.05679).
- (62) A. Perdomo-Ortiz, A. Feldman, A. Ozaeta, S. V. Isakov, Z. Zhu, B. O’Gorman, H. G. Katzgraber, A. Diedrich, H. Neven, J. de Kleer, et al., On the readiness of quantum optimization machines for industrial applications, (unpublished).
- Perdomo-Ortiz et al. (2012) A. Perdomo-Ortiz, N. Dickson, M. Drew-Brook, G. Rose, and A. Aspuru-Guzik, Finding low-energy conformations of lattice protein models by quantum annealing, Sci. Rep. 2, 571 (2012).
- Santra et al. (2014) S. Santra, G. Quiroz, G. Ver Steeg, and D. A. Lidar, Max 2-SAT with up to 108 qubits, New J. Phys. 16, 045006 (2014).
- Rieffel et al. (2015) E. G. Rieffel, D. Venturelli, B. O’Gorman, M. B. Do, E. M. Prystay, and V. N. Smelyanskiy, A case study in programming a quantum annealer for hard operational planning problems, Quant. Inf. Proc. 14, 1 (2015).
- Azinovic et al. (2017) M. Azinovic, D. Herr, B. Heim, E. Brown, and M. Troyer, Assessment of Quantum Annealing for the Construction of Satisfiability Filters, SciPost Phys. 2, 013 (2017).