Computational hardness of spin-glass problems with tile-planted solutions

TL;DR

This paper studies the computational difficulty of spin-glass problems on a lattice with planted solutions, revealing tunable hardness and phase transitions linked to magnetic ordering, using advanced Monte Carlo methods.

Contribution

It introduces a scalable method for planting solutions in spin-glass problems and analyzes the resulting hardness landscape with multiple phase transitions.

Findings

Problems exhibit a wide range of tunable hardness.

Multiple hardness phase transitions are observed.

Harder samples show magnetic ordering transitions.

Abstract

We investigate the computational hardness of spin-glass instances on a square lattice, generated via a recently introduced tunable and scalable approach for planting solutions. The method relies on partitioning the problem graph into edge-disjoint subgraphs, and planting frustrated, elementary subproblems that share a common local ground state, which guarantees that the ground state of the entire problem is known a priori. Using population annealing Monte Carlo, we compare the typical hardness of problem classes over a large region of the multi-dimensional tuning parameter space. Our results show that the problems have a wide range of tunable hardness. Moreover, we observe multiple transitions in the hardness phase space, which we further corroborate using simulated annealing and simulated quantum annealing. By investigating thermodynamic properties of these planted systems, we demonstrate that the harder samples undergo magnetic ordering transitions which are also ultimately responsible for the observed hardness transitions on changing the sample composition.

Tables

Table 1. Table 1: Parameters for the population annealing Monte Carlo simulations. L 𝐿 L is the linear system size, R 𝑅 R is the number of replicas, T 0 subscript 𝑇 0 T_{0} is the lowest temperature, and N T subscript 𝑁 𝑇 N_{T} is the number of temperature steps. At each temperature step, N S = 10 subscript 𝑁 𝑆 10 N_{S}=10 Monte Carlo sweeps per temperature step are performed on each replica. The annealing schedule is linear in the inverse temperature.

$L$	$R$	$T_0$	$N_T$
16	$1\times {10}^5$	0.2	301
24	$2\times {10}^5$	0.2	301
32	$1\times {10}^6$	0.2	301

Table 2. Table 2: A comparison of ⟨ log 10 ⁡ ρ s ⟩ delimited-⟨⟩ subscript 10 subscript 𝜌 𝑠 \langle\log_{10}\rho_{s}\rangle , SA median TTS, and SQA median TTS for a selected set of instance classes. The results are for the system size L = 32 𝐿 32 L=32 .

Instance class	$⟨{\log }_{10}{\rho }_s⟩$	SA TTS ($\mu$s)	SQA TTS ($\mu$s)
$C_1$	2.520(3)	$4.040\left(8\right)\times {10}^3$	$3.11\left(3\right)\times {10}^6$
$C_2$	4.70(2)	$1.9\left(2\right)\times {10}^7$	$4.9\left(6\right)\times {10}^8$
$C_3$	1.6447(2)	$1.45\left(2\right)\times {10}^3$	$2.86\left(6\right)\times {10}^7$
$C_4$	1.44822(7)	$2.51\left(3\right)\times {10}^0$	$2.00\left(3\right)\times {10}^5$
(0.3, 0, 0.7)	3.62(2)	$3.8\left(3\right)\times {10}^6$	$6.1\left(8\right)\times {10}^8$
(0.35, 0, 0.65)	3.75(3)	$3.7\left(2\right)\times {10}^6$	$4.6\left(5\right)\times {10}^8$
(0.4, 0, 0.6)	3.61(3)	$2.18\left(6\right)\times {10}^6$	$2.8\left(2\right)\times {10}^8$
(0.4, 0, 0)	2.39(2)	$2.6\left(3\right)\times {10}^5$	$4.4\left(5\right)\times {10}^8$
(0.45, 0, 0)	2.67(2)	$6.1\left(7\right)\times {10}^5$	$4.2\left(3\right)\times {10}^8$
(0.5, 0, 0)	2.72(2)	$5.3\left(4\right)\times {10}^5$	$2.8\left(3\right)\times {10}^8$

Equations

$${\mathcal{H}}_C=-\sum _{(ij)\in E[C]}{J}_{ij}{s}_i{s}_j{s}_i\in \{\pm 1\},$$

$$\mathcal{H}=\sum _{C\in \mathcal{C}}{\mathcal{H}}_C.$$

$$\exp(s_{0})=\left\{\begin{array}[]{ll}1.33(5)&\mbox{for}\;C_{2}\\ 2.285(2)&\mbox{for}\;C_{3}\\ 4.018(2)&\mbox{for}\;C_{4}\;.\\ \end{array}\right.$$

$${\rho }_s=\underset{R\to \infty }{\lim }R{e}^{{\sum }_i{n}_i\log n_i},$$

$$\text{TTS}(t_{\text{A}})=t_{\text{A}}\frac{\log (1-p_d)}{\log [1-p_s(t_{\text{A}})]},$$

$$g_{\mathcal{O}}=\frac{1}{2}\left(3-[⟨{\mathcal{O}}^4⟩{]}_{\text{av}}/[⟨{\mathcal{O}}^2⟩{]}_{\text{av}}^2\right),$$

$$g_{\mathcal{O}}\sim G_{\mathcal{O}}[L^{1/{\nu }_{\mathcal{O}}}(T-T_c)],$$

$$C_v=1/(k_B{T}^{2}N)[⟨{\mathcal{H}}^2⟩-⟨\mathcal{H}{⟩}^2{]}_{\text{av}},$$

Full text

I Introduction

Optimization problems with many minima occur in a multitude of fields, including finance, engineering, materials sciences, and machine learning. While a range of problems such as shortest path, maximum flow and minimum spanning tree can be elegantly solved with algorithms of a run time growing polynomially with the system size Gibbons (1985), many of the most interesting problems are known to be NP hard, and, as a consequence, polynomial-time algorithms are very unlikely to exist. Among discrete optimization problems many systems can be mapped onto the Ising spin glass Lucas (2014); Stein and Newman (2013) as a central object of study in statistical physics. Its features of frustration and random disorder that are believed to be fundamental to the existence of a spin-glass phase Binder and Young (1986) lead to a multitude of local minima separated by barriers — the complex (free) energy landscape — that are also at the heart of hard optimization problems more generally Stein and Newman (2013); Barahona (1982). Besides the overlap in model systems, statistical physics methods have also proven particularly valuable in elucidating the structure of the solution space of optimization problems Krzakała et al. (2007) and the occurrence of phase transitions in sample hardness as suitable control parameters are tuned Hartmann and Weigt (2005). Moreover, many heuristic techniques for optimization such as simulated annealing (SA) Kirkpatrick et al. (1983), parallel tempering Monte Carlo Swendsen and Wang (1986); Geyer (1991); Marinari and Parisi (1992); Hukushima and Nemoto (1996); Katzgraber et al. (2006a), population annealing Monte Carlo (PAMC) Hukushima and Iba (2003); Machta (2010); Wang et al. (2015a, b), and simulated quantum annealing (SQA) using quantum Monte Carlo Finnila et al. (1994); Kadowaki and Nishimori (1998); Santoro et al. (2002); Heim et al. (2015) have been derived from concepts in statistical and quantum physics.

Recent years have witnessed the advent of special-purpose devices for discrete optimization. Most noteworthy are the commercially available analog D-Wave quantum annealing devices Johnson et al. (2011); Bunyk et al. (2014), which strive to minimize Ising Hamiltonians by exploiting quantum tunneling and superposition. The latest machine, the D-Wave 2000Q, allows for the optimization of instances with up to $2000$ Ising variables, although sparse connectivity of the native graph and the noise due to control errors Pudenz et al. (2014, 2015); Zhu et al. (2016) pose limitations for practical applications. As of now, no experimental evidence of a quantum advantage for generic optimization applications has been discovered, while there are some early indications of a quantum advantage when using the machine as a physical simulator. However, a number of recent studies have shown a speedup for the device over selected classical algorithms for specific classes of synthetic problems Mandrá and Katzgraber (2018); Albash and Lidar (2018). In addition to D-Wave devices, a number of further experimental hardware-based solvers have been introduced, for example, the Coherent Ising machine (optical) Wang et al. (2013); Hamerly et al. (2019), and the complementary metal-oxide semiconductor (CMOS) based Fujitsu Digital Annealer Matsubara et al. (2018); Tsukamoto et al. (2017) studied recently by Aramon et al. Aramon et al. (2019).

With the ensuing recent growth in interest in hardware- and software-based Ising solvers, there is an increasing demand for hard tunable benchmark problems for performance comparisons com (a). Synthetic benchmark problems the ground states of which are known a priori by construction, commonly known as samples with planted solutions Barthel et al. (2002), are particularly advantageous in this regard. Ideally, the hardness of such problems should be readily tunable, and the construction procedure should be scalable to larger system sizes with reasonable computational effort. A number of solution planting schemes have been proposed for the case of short-range Ising spin glasses. For example, Ref. Hen et al. (2015) presented an approach in which constraint satisfaction problems in the form of cycles of couplers (i.e., frustrated loops) are embedded into an arbitrary graph. Although the problem hardness is tunable, the range of coupler strength is uncontrolled and increases with the system size, which leads to amplified control errors in analog devices with finite precision such as the D-Wave machines. Reference King et al. (2015) introduced a variation of the method with limited range of coupling strength that mitigates the precision limitations. However, both methods produce problems that are easier than other standard benchmarks for a wide range of software and hardware-based solvers Coffrin et al. (2019).

Reference Marshall et al. (2016) proposed an adaptive algorithm for generating hard problems by iteratively updating the coupler values based on the time to solution (TTS) of a chosen solver. The computational resources required by the procedure, however, severely restrict the feasible size of the problems. Moreover, the method relies on the assumption that the problem hardness as predicted by the chosen solver may also extend to other solvers. In the “patch planting” method introduced in Ref. Wang et al. (2017), the problem graphs are constructed by stitching together smaller subgraphs via couplers that ensure consistency among boundary spins. Although the method allows for the construction of arbitrarily large problems with no restrictions on coupler values or the graph topology, the resultant problems are less difficult than entirely random problems of comparable size. In addition, as the ground states of the constituent subgraphs need to be determined in advance, the problem construction may require considerable computational effort depending on the size and connectivity of the subgraphs.

Recently, some of the present authors proposed a planted solutions scheme Hamze et al. (2018) that circumvents several weaknesses of the aforementioned methods. Henceforth, let us refer to this approach as “tile planting.” The method is based on decomposing the underlying problem graph into edge-disjoint subgraphs, and embedding a selected set of subproblems in the resulting subgraphs. The subproblems are chosen such that they share a common (local) ground state, which ensures that the ground state of the entire problem is known a priori. Due to the myriad ways in which subproblems can be chosen, one can construct problems with a wide range of computational complexities. The problem construction demands very low computational effort, and facilitates the generation of arbitrarily large problems on the underlying topology.

For completeness, we mention two alternate planting approaches Hamze et al. (2019); Hen (2019) that start from a set of linear equations and then use these to plant a solution for a quadratic Ising system. In particular, the Wishart planted ensemble Hamze et al. (2019) is amenable to theoretical exploration, thus allowing one to make analytical predictions for the typical computational hardness of the planted problems with a given bit precision. While this is a huge bonus, the drawback is that the planted ensemble requires a complete graph and is thus not best suited to benchmark current quantum annealing hardware where an all-to-all connectivity is extremely hard to build Bunyk et al. (2014). However, the ensemble is ideal for the benchmarking of (classical) digital annealer hardware Matsubara et al. (2018); Tsukamoto et al. (2017), as well as optical systems Wang et al. (2013); Hamerly et al. (2019).

In the present paper, we explore the computational hardness of planted problems on square lattices constructed using tile planting. By focusing on the square lattice we are able to exhaustively investigate a large region of the tuning parameter phase space, a task that is impractical for cubic lattices due to the comparatively larger phase space. Although the planarity of the underlying graph renders them exactly solvable in polynomial time using planar graph solvers Schraudolph and Kamenetsky (2009); Galluccio et al. (2001, 2000); Saul and Kardar (1994), the free-energy landscape generally remains challenging for the majority of practical heuristic optimization methods that are agnostic to this special graph structure, and the resulting class of problems has many common features with nonplanar variants of the tile-planting approach Hamze et al. (2018) in, e.g., three space dimensions or nonplanar graphs. It should also be noted that one can break planar methods by superficial modifications such as adding fields, but the purpose here is not to obfuscate our results through such changes.

The remainder of this paper is organized as follows. In Sec. II we outline our planted solutions scheme for square lattices, and discuss the relationship of the constructed problems to constraint satisfaction problems. In Sec. III we investigate the variations in problem difficulty across different instance classes via the population annealing Monte Carlo method. Particular emphasis is given to two specific regions in the parameter space where hardness transitions are observed, which we further investigate using simulated annealing and simulated quantum annealing. Furthermore, we investigate the thermodynamic properties of a selected set of instance classes and discuss the connection between hardness and phase behavior. Section IV presents our concluding remarks.

II Tile planting

We start by describing our planting scheme for square lattices, which will be the focus of this paper. For a demonstration of the method for cubic lattices see Ref. Hamze et al. (2018). Note that the approach is adaptable to arbitrary graphs, with a carefully chosen set of subproblem classes suitable for the underlying topology.

Let us consider a square lattice of linear size $L$ with periodic boundary conditions. We denote the underlying graph as $G=(V,E)$, with vertex set $V$ and edge set $E$. As illustrated in Fig. 1, by partitioning the lattice into a checkerboard pattern, we decompose the graph $G$ into a set of vertex-sharing, unit-cell subgraphs (represented by shaded cells). Note that the subgraphs do not share any edges, and with periodic boundary conditions, each vertex in $V$ appears exactly in two subgraphs. Let $C$ denote the vertices of a given unit-cell subgraph, and let $\mathcal{C}$ represent the set of all such subgraph vertex sets constituting $G$. For each subgraph $C$, we define the Ising Hamiltonian

[TABLE]

which will henceforth be referred to as a subproblem Hamiltonian. We then form the complete Ising Hamiltonian $\mathcal{H}$ over graph $G$ by adding all subproblem Hamiltonians as

[TABLE]

As a straightforward result of constructing $\mathcal{H}$ as a direct sum of subproblem Hamiltonians, if all subproblem Hamiltonians share a common ground-state configuration, this state will also be a ground state of the overall Hamiltonian $\mathcal{H}$. Thus, the ground state energy of $\mathcal{H}$ is also known a priori. By carefully choosing the couplers for the subproblem Hamiltonians from different classes with varying levels of frustration, one can control the computational hardness of the entire problem to a great extent Hamze et al. (2018).

We now proceed to define a set of subproblem classes for the square lattice topology. Let us choose the ferromagnetic ground state, i.e., all spins either up or down as the intended shared ground state of the subproblems. Note that once the problem is constructed one can transform the planted solution from the ferromagnetic ground state to an arbitrary state via a simple gauge transformation. Here we propose four subproblem classes denoted by $\{C_{i}\}$, $i\in\{1,2,3,4\}$, corresponding to bond configurations on square plaquettes (see the illustration in Fig. 2). For the exchange couplings we use two magnitudes $J_{l}$ and $J_{s}$ with $|J_{l}|>|J_{s}|$. Here we arbitrarily choose $|J_{l}|=2$ and $|J_{s}|=1$. The use of two different coupling strengths allows us to tune the degeneracy of the overall optimization problems and hence achieve a variation in hardness of the instances created.

An instance of the subproblem class $C_{i}$ is constructed as follows. First, a chosen edge of the subproblem unit cell is set to the antiferromagnetic coupler value of $-J_{s}$. Then we assign the ferromagnetic value $+J_{s}$ to a randomly selected set of $i-1$ remaining edges, while the rest of the edges are set to $+J_{l}$. This construction procedure leads to three equivalent subproblem representations for the classes $C_{2}$ and $C_{3}$ and unique representations for $C_{1}$ and $C_{4}$, up to overall rotations in multiples of $90^{\circ}$ of the plaquettes (see Fig. 2). From this construction it is clear that subproblems of class $C_{i}$ have $i$ different ground states (apart from overall spin reversal), and the ferromagnetic configurations are always ground states.

When generating planted instances, we first choose a subproblem class for each unit-cell subproblem (i.e., shaded square plaquettes in Fig. 1) based on a chosen probability distribution over classes $\{C_{i}\}$. Then for each subproblem, we randomly choose one of the equivalent class representatives shown in Fig. 2, followed by a random rotation of the plaquette to introduce more disorder. We define each instance class based on the chosen probability distribution over the subproblem classes $\{C_{i}\}$. Let $p_{i}$ denote the probability of choosing subproblems from class $C_{i}$. We can then uniquely identify each instance class in terms of three probability parameters, for instance, $\{p_{1},p_{2},p_{3}\}$, where $p_{1}+p_{2}+p_{3}\leq 1$. As we demonstrate below, one can induce tremendous variations in problem hardness by tuning the probability parameters $\{p_{i}\}$. Moreover, due to the simplistic nature of the construction procedure, hundreds of problem instances with moderate system sizes can be generated in seconds with very low computational cost.

We note that all of the subproblems $C_{i}$ correspond to frustrated squares Toulouse (1977), whereas the remaining plaquettes (white cells in Fig. 1) are randomly frustrated with probability $15/32$, leading to an overall concentration of about 73% frustrated plaquettes in the resulting instances com (b). To study the degeneracy of the full samples resulting from subproblem classes $C_{i}$, we use exact enumeration based on the evaluation of Pfaffians which can be achieved in polynomial time for the case of planar graphs Galluccio et al. (2000). On its own, subproblem $C_{i}$ has $i$ ground states and hence $i^{\cal N}$ is an upper bound to the number of ground states of samples composed of subproblems of class $i$, where ${\cal N}$ is the number of plaquettes (shaded cells in Fig. 1). As a consequence, instances entirely composed of $C_{1}$ have a unique ferromagnetic ground state up to overall spin reversal. For the other classes the connection of subproblems provides additional constraints reducing the degeneracy. For example, for $L=24$ we find the following ground-state entropies $s_{0}=S_{0}/{\cal N}$ per shaded plaquette

[TABLE]

In comparison, the standard $\pm J$ spin-glass model in two space dimensions has $\exp(s_{0})\approx 1.38$, so it turns out to be closest to $C_{2}$ in this respect, which as we will see later on provides the hardest samples among the pure classes $C_{i}$.

As demonstrated in Ref. Hamze et al. (2018), the planted problem construction can be translated to a specific type of constraint satisfaction problem called a tiling puzzle Demaine and Demaine (2007), a well-known example of which is Eternity II. Figure 3 illustrates the steps involved in the derivation of a tiling puzzle from the unit-cell planting formulation. As shown in Fig. 3(a), we first construct a graph with respect to $\mathcal{C}$ (which we call the subproblem lattice $\tilde{G}$) by placing a vertex at the center of each planted cell. Each pair of vertices in $\tilde{G}$ is considered adjacent if and only if the corresponding subgraphs (shaded cells) share a corner vertex. Let us associate each vertex in $\tilde{G}$ with a subproblem lattice variable $\mathbf{s}_{C}$ over the four Ising variables of the corresponding subproblem. The domain of $\mathbf{s}_{C}$ is the set of ground states of the subproblem. However, viable solutions to $\mathbf{s}_{C}$ are constrained by the equality conditions enforced on the shared Ising variables. That is, for two neighboring cells $C$ and $C^{\prime}$ with a shared corner vertex $i$, $s_{C}^{i}=s_{C^{\prime}}^{i}$. A solution to the overall problem should simultaneously satisfy all such equality constraints. By interpreting each equality constraint as a factor over the subproblem lattice variables $\mathbf{s}_{C}$ and $\mathbf{s}_{C^{\prime}}$, one can construct a factor graph Kschischang et al. (2001) associated with $\tilde{G}$ as shown in Fig. 3(b), where the factors are represented by squares. We then perform a Voronoi tessellation Okabe et al. (2000) of $\tilde{G}$ [Fig. 3(c)] which partitions the surrounding space into a set of tilted square regions centered at each vertex, the boundaries of which are shown in thick lines. We can interpret each region as a tiling location, with an allowable set of tiles and orientations as determined by the domain of $\mathbf{s}_{C}$ (the ground-state set of the subproblem). In Fig. 3(c), a hypothetical tile is placed at each location, where blue (dark) and yellow (light) tile faces represent the two Ising spin values ($+1$ and $-1$). The optimal solution to the overall problem can now be interpreted as finding a tile arrangement in which no two adjacent tile faces have a different color; i.e., all factor constraints are satisfied.

III Results

III.1 Simulation details

To compare the typical hardness of planted instances across different problem classes, we primarily use population annealing Monte Carlo Hukushima and Iba (2003); Machta (2010); Wang et al. (2015a, b); Barash et al. (2017); Barzegar et al. (2018). For certain regimes with interesting variations in hardness, we also use SA Kirkpatrick et al. (1983) and SQA Finnila et al. (1994); Kadowaki and Nishimori (1998); Santoro et al. (2002); Heim et al. (2015).

PAMC is a sequential Monte Carlo algorithm somewhat similar in spirit to SA. In contrast to the latter, however, PAMC was designed with equilibrium simulations in mind, and not as an optimization method. It considers an ensemble of copies of the same system and disorder realization that is started at a high temperature and successively cooled down to the lowest temperatures of interest. In addition to equilibrating update steps applied to each copy that are typically implemented as Markov chain Monte Carlo, the population is further relaxed on lowering the temperature by a resampling step that selects and replicates (or eliminates) copies proportional to their relative Boltzmann weight at the new temperature. Large fluctuations in the weights that determine resampling are a signature of the difficulty in equilibrating a given system Machta (2010) and lead to the descendants of only a few original copies dominating the population. Hence, a measure of the effective number of surviving replica families at the lowest temperature allows one to distinguish between hard and easy problems within the instance classes. If $n_{i}$ is the fraction of the population descended from replica $i$ in the initial population (i.e., the size of the $i$th family), the entropic family size $\rho_{s}$ Wang et al. (2015a) is given by

[TABLE]

where $R$ is the initial population size. For a given set of simulation parameters, the larger the value of $\rho_{s}$, the smaller the number of surviving families, and hence the more rugged the problem’s energy landscape. Thus, $\rho_{s}$ provides a measure of hardness for algorithms that are based on local search in the classical energy landscape. As was shown previously, $\rho_{s}$ is highly correlated with other well established hardness metrics such as the integrated autocorrelation time in Markov-chain Monte Carlo and, specifically, in parallel tempering (PT) Monte Carlo Wang et al. (2015a). A unique advantage of $\rho_{s}$ is that it can be measured with relatively low computational cost in comparison to PT autocorrelation times, as well as metrics based on optimized time to solution com (c). This allows us to rigorously investigate a large region of the instance class parameter space that results from random selection of subproblem types.

Table 1 shows the PAMC simulation parameters for the three linear system sizes used in our paper to characterize hardness, $L=16$, $24$, and $32$. $\rho_{s}$ is known to converge with increasing populating size $R$, and $R$ values for different system sizes were chosen accordingly to ensure that this convergence criterion is met by all simulations.

In addition to PAMC simulations, we also perform SA and SQA calculations on a selected subset of problem classes. SA calculations are performed with the optimized multispin code by Isakov et al. Isakov et al. (2015), using a linear annealing schedule with an initial temperature of $T_{\mathrm{i}}=10$ and a final temperature of $T_{\mathrm{f}}=0.2$. For the SQA calculations, we use $64$ Trotter slices and a temperature of $T=0.03125$, while linearly decreasing the transverse field magnitude from $\Gamma=2.5$ to [math]. For both SA and SQA we estimate the time required to find the optimal solution at least once with probability $p_{d}$, i.e., the time to solution Rønnow et al. (2014); Boixo et al. (2014):

[TABLE]

where $t_{\text{A}}$ is the annealing time and $p_{s}$ is the success probability, while $p_{d}=0.99$. We measure the TTS for individual problem instances, and then calculate the median of the distribution of TTS values for a given instance class and system size. We minimize the median TTS in annealing time $t_{\text{A}}$ to obtain the optimal median TTS which takes into account the best tradeoff between the annealing time and the number of independent runs Rønnow et al. (2014). For simplicity, hereafter we use the term time to solution to refer to the aforementioned optimal median TTS obtained for each instance class and system size.

Finally, to gain insight into the physical mechanisms of hardness, we investigate thermodynamic properties of a selected subset of instance classes using PAMC, as well as using the Pfaffian enumeration method of Ref. Galluccio et al. (2000) to determine the exact density of states and hence exact thermodynamic averages for energetic quantities for each sample. Via PAMC simulations, we measure the Binder ratio $g_{\mathcal{O}}$ via

[TABLE]

where the order parameter $\mathcal{O}$ represents either the magnetization $m=(1/N)\sum_{i}s_{i}$ (ferromagnetic order), or the spin overlap $q=(1/N)\sum_{i}s_{i}^{\alpha}s_{i}^{\beta}$ (spin-glass order), with $\alpha$ and $\beta$ representing two independent replicas of the system with the same disorder. Note that $\langle\cdots\rangle$ and $[\cdots]_{\text{av}}$ represent thermal average and disorder average, respectively. As a dimensionless quantity, the Binder ratios for different system sizes are expected to cross close to the point of a continuous phase transition. From the expected scaling behavior Privman (1990),

[TABLE]

it is straightforward to extract the critical temperature $T_{c}$ and the critical exponent $\nu_{\mathcal{O}}$ from a scaling collapse Houdayer and Hartmann (2004); Katzgraber et al. (2006b). In addition to the Binder ratio, we also investigate the specific heat,

[TABLE]

using the planar graph solver, which provides us with a measure free of statistical errors.

III.2 Base instance classes

We first focus on problems constructed using subproblems drawn from a single subproblem class. Let us refer to the resulting four instance classes as base classes. Figure 4 shows the distributions of $\log_{10}\rho_{s}$ for the four base classes, along with the results for random bimodal spin glasses, i.e., $J_{ij}\in\{\pm 1\}$. For the bimodal instances, we use the same PAMC simulation parameters as used for the base classes. As Fig. 4 shows, $\rho_{s}$ has a very broad distribution such that the distribution of $\log_{10}\rho_{s}$ is at least similar to a normal shape, and we hence use $\log_{10}\rho_{s}$ as our preferred form of the hardness metric.

Figure 4 shows that both the $C_{2}$ base class and the bimodal class have very broad distributions, with the bimodal distribution being strongly skewed towards hard problem instances. In contrast, the $C_{3}$ and $C_{4}$ base classes have quite narrow distributions which appear to have vanishing overlap with the bimodal case, clearly indicating that the problems belonging to these classes are easier than the bimodal spin glasses. $C_{4}$ has the leftmost distribution among all classes, and hence produces the easiest instances, at least according to the $\rho_{s}$ metric. On the other hand, the distribution for $C_{2}$ lies the farthest to the right, only overlapping with the tail of the bimodal distribution. Thus, we conclude that $C_{2}$ instances are at least as hard as the hardest instances belonging to the bimodal class. Finally, based on the $\rho_{s}$ hardness metric, we can list the four base classes in sequence according to increasing order of hardness as $C_{4}$, $C_{3}$, $C_{1}$, and $C_{2}$. We will later demonstrate that by mixing the four subproblem types one can construct problems with varying degrees of complexity, with the level of hardness as measured by $\rho_{s}$ lying between those of the $C_{4}$ and $C_{2}$ base classes.

Before we proceed to problems constructed using mixtures of subproblem types, let us investigate thermodynamic properties of the base classes to uncover any noteworthy behavior such as the presence of phase transitions. We first present results for the $C_{3}$ and $C_{4}$ classes due to their trivial thermodynamic behavior in comparison to other base classes. Figure 5(a) shows the specific heat $C_{v}$ and Fig. 5(b) shows the Binder ratio $g_{q}$ as a function of temperature for the $C_{3}$ base class. In Fig. 6 we present the same quantities for the $C_{4}$ base class. Here, the specific heat is determined using the planar graph solver Galluccio et al. (2000) whereas the magnetization and the Binder ratio are measured using PAMC simulations. For the PAMC simulations, the results for system sizes $L=16$, $24$, and $32$ were generated using the parameters presented in Table 1. For $L=10$, we use $R=1\times 10^{5}$ and $T=201$, and for $L=48$ we use $R=2\times 10^{5}$ and $T=401$, with $T_{0}=0.2$ and $N_{S}=10$ for both cases. To ensure thermal equilibration, we require the condition $R/\rho_{s}>100$ Wang et al. (2015a) to be met, and when it is not satisfied by an instance it is rerun by increasing the number of temperature steps $N_{T}$. As is demonstrated in Figs. 5 and 6, the specific heat for both $C_{3}$ and $C_{4}$ classes does not show any system-size dependence. Moreover, the Binder cumulants $g_{q}$ for different system sizes do not show any sign of an intersection, indicating the absence of ferromagnetic or spin-glass transitions in these two base classes, which hence show properties characteristic of disordered spin systems in dimensions below the lower critical dimension. The implied absence of critical slowing down and any significant barriers in the free-energy landscape results in the relative ease of equilibrating the samples that is reflected in the small values of $\rho_{s}$ shown in Fig. 4.

In Fig. 7, we present results for the specific heat $C_{v}$, the average magnetization $[\langle|m|\rangle]_{\text{av}}$, and the Binder ratio $g_{m}$ for the $C_{1}$ base class. The specific heat and the average magnetization resemble that of the ferromagnetic Ising model, with a logarithmic divergence of $C_{v}$ and the magnetization scaling to zero at the transition point $\sim(T-T_{c})^{1/8}$ McCoy and Wu (1973). Moreover, the Binder parameter values $g_{m}$ for different system sizes clearly intersect at a common temperature $T_{c}$. Using a scaling collapse Katzgraber et al. (2006b), we determine the critical temperature to be $T_{c}=1.834(2)$, with $\nu_{m}=1.001(5)$. Our results indicate that despite being a disordered and frustrated system the $C_{1}$ base class has predominantly ferromagnetic features. This ferromagnetic ordering is a necessary consequence of the fact that these instances only have a single ground state (excluding spin-flip symmetry), which is the planted ferromagnetic solution.

Finally, we explore properties of the $C_{2}$ base class. In Fig. 8, we show the specific heat $C_{v}$, magnetization $[\langle|m|\rangle]_{\text{av}}$, and Binder parameter $g_{m}$ for the $C_{2}$ class. Here, system sizes beyond $L=32$ were not simulated due to difficulties in thermalizing $C_{2}$ instances. Figure 8(c) shows that $g_{m}$ for different system sizes intersect in the temperature interval $T\in\left(0.4,0.55\right)$, which indicates that some ferromagnetic ordering occurs. The fact that the curves do not intersect at a common point implies that corrections to scaling are significant for these system sizes. Hence, we do not attempt to extract the critical temperature via a scaling analysis. The average magnetization [Fig. 8(b)] also shows features reminiscent of a ferromagnetic transition, although, as $T\rightarrow 0$, $[\langle|m|\rangle]_{\text{av}}$ appears to approach a value close to $0.8$, rather than $1.0$ as in the case of the ferromagnetic Ising model. The specific heat [Fig. 8(a)] does not show a discernible system-size dependence, at least for the system sizes considered.

This peculiar ordering behavior can be understood by considering the clustering behavior of strong-ferromagnetic ($+2$) bonds in $C_{2}$ instances. In $C_{2}$ instances, $50\%$ of the couplers are $+2$ bonds. By investigating the percolation properties of $+2$ bonds, we find that the spanning probability asymptotically approaches unity in the thermodynamic limit (see Fig. 9). This indicates that the lattice is in the percolating phase, and that in the thermodynamic limit there is a nonvanishing fraction of sites in the largest cluster of $+2$ bonds. We find that for the lattice sizes considered approximately 80% of the spins belong to the largest bond cluster. Note that this is in contrast to an uncorrelated percolation problem in its ordered phase, where close to 100% of lattice sites would belong to the percolating cluster. The remaining spins either belong to smaller $+2$-bond clusters scattered around the lattice or they are isolated sites with no $+2$ bonds attached (see Fig. 10 for an illustration). Due to intrinsic properties of the $C_{2}$ lattices, $+2$-bond clusters and “isolated” spins have equal number of weak-ferromagnetic ($+1$) and weak-antiferromagnetic ($-1$) bonds attached. As a result, the flipping of an isolated spin, or the simultaneous flipping of all the spins in a $+2$-bond cluster, leaves the energy unchanged. By examining the ground-state spin configurations, we find that all $+2$ bonds are satisfied in the ground state, that is, spins within each cluster are aligned with one another. Thus, the ground-state degeneracy in $C_{2}$ instances arises from the energy-free flipping of entire clusters or isolated spins. As the critical temperature is approached from above, each cluster would undergo a transition from a randomly oriented (paramagnetic) state to a ferromagnetically aligned state. This is demonstrated in Fig. 11, where we show the average magnetization $\langle|m|\rangle$ as a function of temperature for the three largest clusters of a single $C_{2}$ instance. As $T\rightarrow 0$, we observe that $\langle|m|\rangle\rightarrow 1$ as the spins within each cluster are ferromagnetically aligned.

A main factor leading to the observed hardness in $C_{2}$ instances is the critical slowing down that occurs as the system undergoes a ferromagnetic transition at a considerably low temperature $T\in\left(0.4,0.55\right)$. We point out that cluster algorithms Swendsen and Wang (1987); Wolff (1989) could be effective in reducing critical slowing down and improving the performance. However, an analysis on the performance of problem-specific algorithms is beyond the scope of this paper, and hence was not attempted. The ferromagnetic nature of $C_{1}$ and $C_{2}$ base-class instances is not surprising when considering the fraction of strong-ferromagnetic bonds in $C_{1}$ and $C_{2}$ plaquettes. Three out of four bonds in a $C_{1}$ plaquette are strong-ferromagnetic bonds, and, therefore, $C_{1}$ instances have the strongest ferromagnetic properties among all classes. With the replacement of each strong-ferromagnetic bond with a weak-ferromagnetic bond, the plaquettes progressively become less ferromagnetic and more frustrated in the order $C_{1}\rightarrow C_{2}\rightarrow C_{3}\rightarrow C_{4}$. While both $C_{1}$ and $C_{2}$ order ferromagnetically with the accompanying critical slowing down, $C_{2}$ also shows a significant ground-state degeneracy [see Eq. (3)], and the expected energy barriers between the multitude of ground states lead to additional slow dynamics at low temperatures that is absent for $C_{1}$, explaining the much weaker computational challenge provided by $C_{1}$ instances as compared to $C_{2}$ samples.

III.3 Mixtures of two subproblem types

We now investigate planted instances constructed using subproblems drawn from two plaquette types. We first focus on two regimes where interesting variations in hardness can be observed: mixtures of $C_{1}$ and $C_{3}$ plaquettes, and mixtures of $C_{1}$ and $C_{4}$ plaquettes. Rather than comparing $\log_{10}\rho_{s}$ distributions, here we consider the disorder averages $\langle\log_{10}\rho_{s}\rangle$ computed over $200$ problems drawn from each instance class per system size. Figure 12(a) shows $\langle\log_{10}\rho_{s}\rangle$ for $C_{1}$–$C_{3}$ plaquette mixtures plotted as a function of $p_{1}$, which is the probability of choosing subproblems from class $C_{1}$. Thus, $p_{1}=1$ corresponds to the instance class with all subproblems drawn from $C_{1}$ plaquette type ($C_{1}$ base class), whereas $p_{1}=0$ represents instances with all $C_{3}$ subproblems ($C_{3}$ base class). Results are shown for three system sizes, $L=16$, $24$, and $32$. Figures  12(b) and  12(c) show the optimized median TTS for SA and SQA, again estimated using $200$ problem instances. Figure 13 shows the same quantities for $C_{1}$–$C_{4}$ plaquette mixtures. Results from all three algorithms clearly indicate peaks in problem hardness at $p_{1}\approx 0.35$ for $C_{1}$–$C_{3}$ mixtures, and $p_{1}\approx 0.45$ for $C_{1}$–$C_{4}$ mixtures, i.e., an easy-hard-easy transition. This easy-hard-easy transition is akin to the phase transition one observes in Boolean satisfiability ($k$-SAT) problems at a certain clauses-to-variables ratio, below which the formula is satisfiable and above which it is unsatisfiable Bollobás et al. (2001). However, we point out that our problems are exactly solvable in polynomial time, and therefore the observed hardness transition is more analogous to the transition in 2-SAT problems than the $k>0$ case for which one observes a discontinuous transition. A similar hardness peak has been observed for the frustrated loop problems on the Chimera lattice as the loop density is varied Hen et al. (2015).

It is interesting to explore whether these hardness peaks coincide with or are driven by thermodynamic phase transitions. Based on the thermodynamic behavior of the instance classes at $p_{1}=0$ ($C_{3}$ and $C_{4}$) and $p_{1}=1$ ($C_{1}$), it is reasonable to speculate that at zero temperature, as $p_{1}$ is varied from [math] to $1$, the system may undergo a transition from a disordered phase to a ferromagnetic phase. To further explore this hypothesis, we focus on the $C_{1}$–$C_{3}$ plaquette mixtures and measure $[\langle|m|\rangle]_{\text{av}}$ and $g_{m}$ at the lowest temperature $T_{\text{min}}=0.2$ for different values of $p_{1}\in[0,1]$. The results are shown in Fig. 14 as a function of $p_{1}$, where Fig. 14(a) shows $[\langle|m|\rangle]_{\text{av}}$ and Fig. 14(b) shows the quantity $Q_{m}=-\ln(1-g_{m})$. System sizes beyond $L=32$ were not simulated due to difficulties in thermalizing instance classes in the vicinity of the hardness transition. Here, the quantity $Q_{m}=-\ln(1-g_{m})$ was chosen over $g_{m}$ as it has a smaller curvature near the transition Möbius and Rössler (2009), allowing for a more accurate estimation of the crossing point. Figure 14(a) shows that $[\langle|m|\rangle]_{\text{av}}\rightarrow 1$ as $p_{1}\rightarrow 1$, consistent with the ferromagnetic properties of the $C_{1}$ base class. As Fig. 14(b) shows, $Q_{m}$ for different system sizes intersect at a nonzero $p_{1}$ value, further confirming the presence of a ferromagnetic transition. Assuming a scaling behavior of the form $Q_{m}\sim F[L^{1/\kappa}(p_{1}-p_{1}^{c})]$, we perform a finite-size scaling analysis to determine the point of intersection $p_{1}^{c}$. From finite-size scaling we expect $\kappa=\nu$, and universality implies $\nu=1$ as we are dealing with an Ising transition in two dimensions. This leads to an acceptable fit and the estimate $p_{1}^{c}=0.352(7)$.

A visual inspection of the hardness curves in Fig. 12 shows that the peak in hardness occurs in the vicinity of this transition point. To obtain a rough estimate of the asymptotic value of the hardness peak position as $L\rightarrow\infty$, for each $L$, we fit cubic polynomials to the data points of $\langle\log_{10}\rho_{s}\rangle$ close to the peak, and estimate the peak positions $(p_{1}^{h})_{L}$. Fig. 15 shows $(p_{1}^{h})_{L}$ plotted against $1/L$. By fitting a function of the form $(p_{1}^{h})_{L}=A+\lambda/L$, we determine the asymptotic peak position to be $p_{1}^{h}=A=0.370(6)$, where we have estimated the error bar using bootstrap resampling. The close proximity of the transition point $p_{1}^{c}$ to the asymptotic hardness peak position $p_{1}^{h}$ suggests that the hardness transition is, indeed, driven by the ferromagnetic transition.

We now turn to the thermodynamic properties of $C_{1}$–$C_{4}$ mixtures. Fig. 16(a) and Fig. 16(b), respectively, show $[\langle|m|\rangle]_{\text{av}}$ and $Q_{m}$ as functions of $p_{1}$ for $C_{1}$–$C_{4}$ mixtures. Both $[\langle|m|\rangle]_{\text{av}}$ and $Q_{m}$ behave in qualitatively similar ways to the corresponding quantities for $C_{1}$–$C_{3}$ mixtures. From a finite-size scaling analysis, we estimate the intersection point in $Q_{m}$ to be $p_{1}^{c}=0.464(6)$. An inspection of the hardness curves in Fig. 13 indicates that the peak in hardness occurs in the proximity of this transition point. By investigating the scaling behavior of the peak position as a function of $1/L$, we obtain a rough estimate of its asymptotic value, $p_{1}^{h}=0.50(1)$. The consistent closeness between the intersection points in $Q_{m}$ and the hardness peak positions corroborates that in both $C_{1}$–$C_{3}$ and $C_{1}$–$C_{4}$ mixtures the hardness transitions are driven by magnetic ordering.

It is noteworthy to mention that apart from the thermodynamic properties we have also analyzed some heuristic rules based on local measures of frustration Kobe and Klotz (1995), but these quantities proved not to be predictive of the hardness transitions or the relative levels of hardness of the base classes Miyazaki (2013).

Among problems constructed using a single subproblem type, $C_{2}$ base class instances are the hardest according to the $\langle\log_{10}\rho_{s}\rangle$ metric (see Fig. 4 and Table 2). Therefore it is interesting to investigate how problem hardness is affected when mixing $C_{2}$ plaquettes with other plaquette types. Fig. 17 shows $\langle\log_{10}\rho_{s}\rangle$ [panel (a)], median SA TTS [panel(b)], and median SQA TTS [panel (c)] for the three types of binary subproblem mixtures involving $C_{2}$: $C_{1}$–$C_{2}$, $C_{2}$–$C_{3}$, and $C_{2}$–$C_{4}$. The results are for the largest system $L=32$. The results are plotted against $p_{2}$, which is the fraction of $C_{2}$ plaquettes. All three hardness metrics indicate that problem hardness monotonically increases with $p_{2}$ for all three types of binary mixtures. Interestingly, one can observe variations in the relative hardness levels of the three types of binary mixtures across different $p_{2}$ values as well as different algorithms. For example, SQA finds $C_{2}$–$C_{3}$ plaquette mixtures to be harder than both $C_{1}$–$C_{2}$ and $C_{2}$–$C_{4}$ mixtures for the entire range of $p_{2}\in[0,1)$, whereas according to SA TTS and $\langle\log_{10}\rho_{s}\rangle$ the hardness of $C_{1}$–$C_{2}$ mixtures surpasses that of $C_{2}$–$C_{3}$ mixtures as $p_{2}$ approaches zero.

In Table 2, we present a quantitative comparison of $\langle\log_{10}\rho_{s}\rangle$, SA median TTS, and SQA median TTS for a selected set of instance classes, including the four base classes, three data points in the vicinity of the hardness peak in $C_{1}$–$C_{3}$ mixtures, {(0.3, 0, 0.7), (0.35, 0, 0.65), (0.4, 0, 0.6)}, and three points close to the hardness peak in $C_{1}$–$C_{4}$ mixtures, {(0.4, 0, 0), (0.45, 0, 0), (0.5, 0, 0)}. The results are for the system size $L=32$. When comparing results for the four base classes we find that, according to all three algorithms, $C_{4}$ instances are the easiest to solve whereas $C_{2}$ instances are the hardest. A comparison of the solver performances for the $C_{1}$ and $C_{3}$ base classes reveals an interesting feature with regard to the relative hardness levels of the two classes. According to the two classical algorithms PAMC and SA, instances from the $C_{1}$ base class are harder than those from the $C_{3}$ class, as conveyed by relatively higher values of $\langle\log_{10}\rho_{s}\rangle$ and SA TTS. However, SQA TTS for the $C_{1}$ base class is lower than that for $C_{3}$, indicating that SQA finds $C_{1}$ base class instances easier to solve than $C_{3}$ instances. We also note that $\langle\log_{10}\rho_{s}\rangle$ and SA TTS for the $C_{2}$ class are noticeably higher than those for the instance classes in the vicinity of the $C_{1}$–$C_{3}$ and $C_{1}$–$C_{4}$ transitions. This indicates that for the system size considered, according to the two classical algorithms $C_{2}$ problems are harder than the instance classes near the transitions. However, when comparing SQA results, we note that the TTS for the $C_{2}$ class is within the error bars of the data points with the highest TTS values near the $C_{1}$–$C_{3}$ and $C_{1}$–$C_{4}$ transitions, i.e., (0.3, 0, 0.7) and (0.4, 0, 0).

III.4 Mixtures of three subproblem types

Finally, we explore variations in problem hardness among instance classes constructed using three subproblem types. Here we discretize the parameter space into a grid by incrementing each probability parameter $p_{i}$ by a step size of $0.1$. We perform PAMC simulations for the instance classes defined on the grid points that represent mixtures of three subproblem types. $\langle\log_{10}\rho_{s}\rangle$ is computed by averaging over $200$ problems per instance class. Figure 18 shows $\langle\log_{10}\rho_{s}\rangle$ for $L=32$ for the four types of three-plaquette combinations: $C_{1}$–$C_{2}$–$C_{3}$ [Fig. 18(a)], $C_{1}$–$C_{2}$–$C_{4}$ [Fig. 18(b)], $C_{1}$–$C_{3}$–$C_{4}$ [Fig. 18(c)], and $C_{2}$–$C_{3}$–$C_{4}$ [Fig. 18(d)]. The results are represented as heat maps, with $\langle\log_{10}\rho_{s}\rangle$ for the points between the grid points estimated via interpolation. $x$ and $y$ axes in each plot represent two of the free probability parameters that uniquely define instance classes in the mixture. The color grade from dark red to light yellow represents low (red/dark) to high (yellow/light) values of $\langle\log_{10}\rho_{s}\rangle$.

First focusing on $C_{1}$–$C_{2}$–$C_{3}$ [Fig. 18(a)] and $C_{1}$–$C_{2}$–$C_{4}$ [Fig. 18(b)] mixtures, we note that $p_{2}=0$ in Fig. 18(a) and Fig. 18(b), respectively, correspond to $C_{1}$–$C_{3}$ and $C_{1}$–$C_{4}$ plaquette mixtures investigated in Sec. III.3. For any fixed value of $p_{2}$, as $p_{1}$ approaches from [math] to $1$, one observes a change in the color grade of the form dark $\rightarrow$ light $\rightarrow$ dark, which corresponds to an easy-hard-easy transition. For further clarification, we show in Fig. 19 the measured values of $\langle\log_{10}\rho_{s}\rangle$ against $p_{1}$ at fixed $p_{2}$ values for $C_{1}$–$C_{2}$–$C_{3}$ [panel (a)] and $C_{1}$–$C_{2}$–$C_{4}$ [panel (b)] mixtures. The results clearly show that the hardness transitions discussed in Sec. III.3 for the case of $p_{2}=0$ are also prevalent for nonzero $p_{2}$ in the multidimensional parameter space of three-plaquette combinations. Generalizing our conclusions on the $C_{1}$–$C_{3}$ and $C_{1}$–$C_{4}$ phase transitions, we claim that all easy-hard-easy transitions in the three-dimensional parameter space are driven by magnetic ordering transitions. As discussed in Sec. III B, the number of strong-ferromagnetic (+2) bonds in each plaquette controls the level of frustration, and the plaquettes progressively become less frustrated and more ferromagnetic in the order $C_{4}\rightarrow C_{3}\rightarrow C_{2}\rightarrow C_{1}$. As such, regions in the parameter space with a high concentration of $C_{1}$ plaquettes are predominantly ferromagnetic, whereas regions with a high concentration of $C_{4}$ and/or $C_{3}$ plaquettes are predominantly disordered. The light-colored (hard) regions in the heat maps separating the ferromagnetic and disordered problem spaces are characterized by either a high concentration of moderately ferromagnetic $C_{2}$ plaquettes or a balanced mixture of $C_{1}$ and $C_{3}$/$C_{4}$ plaquettes which results in an overall moderate concentration of +2 bonds. A transition of the form easy $\rightarrow$ hard $\rightarrow$ easy occurs as one varies the subproblem composition such that the overall concentration of +2 bonds changes from low (disordered) to high (ferromagnetic).

Fig. 19 further shows that for both $C_{1}$–$C_{2}$–$C_{3}$ and $C_{1}$–$C_{2}$–$C_{4}$ mixtures the peak in the hardness shifts to lower $p_{1}$ values as $p_{2}$ is increased. As the fraction of $C_{2}$ plaquettes increases, the concentration of +2 bonds also increases, which causes the ferromagnetic transition to occur at lower values of $p_{1}$. Consequently, for a fixed value of $p_{1}$, increasing the fraction of $C_{2}$ plaquettes leads to a transition from a disordered phase to a ferromagnetic phase, given that the fraction of $C_{1}$ plaquettes is small enough so that the initial phase is disordered. This manifests as an easy-hard-easy transition at fixed $p_{1}$ values as can be observed in Fig. 18(a) and (b). For clarity, in Fig. 20 we show projected curves of $\langle\log_{10}\rho_{s}\rangle$ at fixed $p_{1}$ values, in which we can observe the aforementioned transitions as $p_{2}$ is varied.

Another hardness transition can be observed in the parameter space of $C_{1}$–$C_{3}$–$C_{4}$ plaquette mixtures [Fig. 18(c)], as either $p_{3}$ or $p_{4}$ approaches from [math] to $1$ (or, equivalently, when decreasing $p_{1}$ from $1$ to [math]). This is a result of the decreasing concentration of $C_{1}$ plaquettes, which transforms the problem characteristics from predominantly ferromagnetic to disordered. $\langle\log_{10}\rho_{s}\rangle$ for $C_{2}$–$C_{3}$–$C_{4}$ mixtures [Fig. 18(d)], however, does not show evidence for a easy-hard-easy transition in the corresponding parameter space. This is expected since the parameter space of $C_{2}$–$C_{3}$–$C_{4}$ mixtures does not contain a region with predominantly ferromagnetic properties due to the complete lack of $C_{1}$ plaquettes. The hardness monotonically increases as either $p_{3}$ or $p_{4}$ is decreased or, in other words, as the fraction of $C_{2}$ plaquettes is increased.

Based on the results of all four types of three-plaquette mixtures, we identify that for the system sizes considered $\langle\log_{10}\rho_{s}\rangle$ for all instance classes are bounded by those of $C_{4}$ (lowest value) and $C_{2}$ (highest value) base classes. As a final note, we point out that the quantitative comparison of hardness levels across different instance classes as characterized by $\langle\log_{10}\rho_{s}\rangle$ may not exactly carry over to other solvers. Exploring the entire parameter space with other solvers using a TTS based metric is not practical due to the sheer number of instance classes involved, and thus was not attempted. However, we believe that the general trends in hardness we observe through $\langle\log_{10}\rho_{s}\rangle$, such as the hardness transitions, are common among the solvers considered here and, more general, quantum-inspired optimization techniques. This is corroborated by the comparison of multiple solver performances presented in Sec. III.3 for two sub-regions of the parameter space, where all three algorithms confirm the presence of hardness transitions.

Given the planarity of the underlying graph, it is clear from general results that specially crafted algorithms can find the ground state as well as the partition function in polynomial time Barahona (1982). This does not imply, however, that more general techniques such as Monte Carlo simulations that are based on (more or less) local explorations of the energy landscape will show polynomial scaling. Instead, due to the complex structure of states with many local minima separated by energy barriers one expects an asymptotically exponential scaling of the computational effort of such approaches on the systems considered here. Fig. 21 shows the scaling of the optimal median TTS as a function of system size $L$ for SA and SQA for two representative instance classes: the $C_{2}$ base class [panel (a)] and $C_{1}$–$C_{3}$ binary mixtures with $p_{1}=0.5$ [panel (b)]. We find that our results are best described by polynomial fits of the form $\text{TTS}=aL^{b}$. This suggests that for the range of system sizes considered the scaling of the TTS has not yet reached its asymptotically exponential behavior. For the $C_{2}$ base class, the fits yield scaling coefficients $b=9.43(4)$ for SA, and $b=7.3(1)$ for SQA. For the $C_{1}$–$C_{3}$ mixtures, the scaling coefficients are $b=6.44(5)$ for SA, and $b=5.9(1)$ for SQA. A comparison of the scaling coefficients shows that for the two instance classes considered quantum-inspired SQA outperforms classical SA in terms of scaling. However, one should bear in mind that there are superior classical algorithms for planar graphs, for example, parallel tempering with isoenergetic cluster moves (PT+ICM) Zhu et al. (2015). In fact, a scaling analysis using PT+ICM for the $C_{1}$–$C_{3}$ mixtures with $p_{1}=0.5$ has shown that PT+ICM scales better than SQA Mandrà et al. (unpublished).

IV Summary

We have investigated the computational hardness of planted spin-glass problems on a square lattice topology generated by the method of “tile planting.” Our planting scheme is based on partitioning the underlying problem graph into edge-disjoint subgraphs, and embedding subproblems chosen from predefined classes over the subgraphs. Using the PAMC method, we have mapped out a large region of the hardness phase space. For a selected subset of problem classes, we have also performed SA and SQA calculations and measured the time-to-solution performance metric.

As the building blocks for problem construction, we have considered four subproblem types (i.e. frustrated plaquettes) with different levels of frustration: $C_{1}$, $C_{2}$, $C_{3}$ and $C_{4}$. When comparing instances constructed solely using a single subproblem type, we find instances belonging to the $C_{4}$ base class to be the easiest, while we find $C_{2}$ base class instances to be the hardest. $C_{1}$ and $C_{3}$ base class instances have moderate levels of hardness. By investigating thermodynamic properties, we have shown that $C_{3}$ and $C_{4}$ base classes have properties characteristic of disordered spin systems, whereas the $C_{1}$ class has strong ferromagnetic properties. The $C_{2}$ base class exhibits properties characteristic of both ferromagnetic and disordered systems, as indicated by a nonzero (but relatively low) critical temperature and a significant ground-state degeneracy. We conjecture that the hardness in $C_{2}$ problems arises as a combined effect of the critical slowing down that occurs at the low critical temperature and the presence of local minima. We have shown that by mixing different subproblem types one can achieve tremendous variations in problem hardness. In particular, we observe several hardness transitions in the phase space, which can be identified in Fig. 18 as changes in the color grade of the form dark $\rightarrow$ light $\rightarrow$ dark. The origin of these transitions can be attributed to underlying ferromagnetic phase transitions that occur as the concentration of strong-ferromagnetic (+2) bonds is varied.

Due to the highly tunable hardness, scalability, and ease of implementation, we believe that our method could be very useful for generating benchmark problems for novel optimization methods implemented both in hardware and software. For the reader interested in constructing benchmark problems, we point out that Fig. 18 can be used as a visual guide for selecting easy versus hard problems, in which easy and hard regimes can be identified as dark and light regions, respectively. In addition, Table 2 presents a quantitative comparison of the different performance metrics for the four base classes and data points close to the hardness transitions in $C_{1}$–$C_{3}$ and $C_{1}$–$C_{4}$ binary mixtures. While $C_{2}$ instances are the hardest according to the classical algorithms PAMC and SA, according to SQA results, the instance classes in the vicinity of the hardness peaks can be as hard as $C_{2}$ instances within the error bars. A convenient way of systematically tuning hardness over a wide range is to consider binary subproblem mixtures involving $C_{2}$ (see Fig. 17), for which we expect the hardness to increase monotonically as the fraction of $C_{2}$ plaquettes is increased. $C_{2}$–$C_{4}$ mixtures can be particularly useful in this regard, as one can tune the hardness from the easiest regime $C_{4}$ to the hardest regime $C_{2}$, which corresponds to an increase in the time to solution over several orders of magnitude in the respective scales for both SA and SQA.

Acknowledgements.

We thank Mario Könz for providing his multispin simulated quantum annealing code to perform the calculations presented in this paper. D.P. would like to acknowledge Amin Barzegar and Chris Pattison for stimulating discussions. J.R. thanks Ryoji Miyazaki for informative discussions. F.H. thanks Cathy McGeogh, Andrew King, Paul Bunyk, and Fiona Hanington for helpful remarks. This research is based upon work supported in part by the Office of the Director of National Intelligence (ODNI), Intelligence Advanced Research Projects Activity (IARPA), via MIT Lincoln Laboratory Air Force Contract No. FA8721-05-C-0002. 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 ODNI, IARPA, or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purpose 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 high performance computing resources. M.W. acknowledges support by the European Commission through the IRSES network DIONICOS under Contract No. PIRSES-GA-2013-612707.