Patch-planting spin-glass solution for benchmarking

TL;DR

This paper presents a method to generate large spin-glass instances with known solutions by stitching together smaller solved patches, enabling benchmarking of algorithms with controllable complexity.

Contribution

The authors introduce a novel patch-planting algorithm for creating large, structured spin-glass problems with known solutions, facilitating benchmarking and complexity analysis.

Findings

Patch-planting effectively creates large instances with known solutions.

The complexity of patch-planted problems can be tuned via patch size and number.

Patch-planted instances show different scaling behavior compared to random problems.

Abstract

We introduce an algorithm to generate (not solve) spin-glass instances with planted solutions of arbitrary size and structure. First, a set of small problem patches with open boundaries is solved either exactly or with a heuristic, and then the individual patches are stitched together to create a large problem with a known planted solution. Because in these problems frustration is typically smaller than in random problems, we first assess the typical computational complexity of the individual patches using population annealing Monte Carlo, and introduce an approach that allows one to fine-tune the typical computational complexity of the patch-planted system. The scaling of the typical computational complexity of these planted instances with various numbers of patches and patch sizes is investigated and compared to random instances.

Tables

Table 1. Table 1: Simulation parameters for the three-dimensional EA model experiments using population annealing Monte Carlo. Here, L 0 subscript 𝐿 0 L_{0} is the patch size, L 𝐿 L is the linear system size, R 𝑅 R is the number of replicas used in the simulation, T 0 = 1 / β 0 subscript 𝑇 0 1 subscript 𝛽 0 T_{0}=1/\beta_{0} is the lowest temperature simulated, N T subscript 𝑁 𝑇 N_{T} is the number of temperature steps (evenly spaced in β 𝛽 \beta ) in the annealing schedule, BC is the type of boundary condition [either periodic boundary condition (PBC) or free boundary condition (FBC)], and N sa subscript 𝑁 sa N_{\rm sa} is the number of disorder realizations studied. For each replica, N S = 10 subscript 𝑁 𝑆 10 N_{S}=10 Monte Carlo sweeps are performed at each temperature during the anneal. Data for PBC with L = 8 𝐿 8 L=8 are taken from Ref.  Wang et al. ( 2015a ) .

$L_0$	$L$	$R$	$T_0$	$N_T$	BC	$N_{\mathrm{sa}}$
4	4	$4\times {10}^3$	$0.2$	$101$	FBC	$345600$
4	8	${10}^4$	$0.2$	$101$	FBC	$5000$
4	12	$5\times {10}^4$	$0.2$	$201$	FBC	$5120$
4	16	$2\times {10}^5$	$0.2$	$301$	FBC	$1877$
4	20	${10}^6$	$0.2$	$401$	FBC	$194$
5	5	${10}^4$	$0.2$	$101$	FBC	$345600$
5	10	${10}^5$	$0.2$	$201$	FBC	$5000$
6	6	$2\times {10}^4$	$0.2$	$101$	FBC	$41472$
6	12	${10}^5$	$0.2$	$201$	FBC	$1752$
8	8	$5\times {10}^4$	$0.2$	$201$	FBC	$23358$
8	16	${10}^6$	$0.2$	$301$	FBC	$624$
10	10	${10}^6$	$0.2$	$301$	FBC	$8000$
10	20	$2\times {10}^6$	$0.2$	$401$	FBC	$260$
8	8	${10}^5$	$0.2$	$101$	PBC	$5099$
12	12	${10}^6$	$0.2$	$201$	PBC	$3812$

Table 2. Table 2: Statistics of the D-Wave 2X quantum annealer success probability p succ subscript 𝑝 succ p_{\rm{succ}} in Fig.  11 for random instances and patched instances with different number of patches M 𝑀 M . There are 10 3 superscript 10 3 10^{3} instances each, and the patched instances were chosen from the hardest ones out of 10 4 superscript 10 4 10^{4} instances in each class. p min subscript 𝑝 min p_{\rm{min}} , p max subscript 𝑝 max p_{\rm{max}} , and p ave subscript 𝑝 ave p_{\rm{ave}} are the minimum, maximum and average values of p succ subscript 𝑝 succ p_{\rm{succ}} , respectively and f 𝑓 f is the fraction of instances with p succ = 0 subscript 𝑝 succ 0 p_{\rm{succ}}=0 .

	Random	$M=2$	$M=3$	$M=4$
$p_{\min }$	$0$	$0$	$0$	$0$
$p_{\max }$	$0.00240\left(32\right)$	$0.00137\left(59\right)$	$0.0326\left(32\right)$	$0.173\left(9\right)$
$p_{\mathrm{ave}}$	$0.0000204\left(46\right)$	$0.0000129\left(24\right)$	$0.00115\left(10\right)$	$0.0127\left(7\right)$
$f$	$0.831\left(12\right)$	$0.775\left(13\right)$	$0.185\left(12\right)$	$0.011\left(3\right)$

Equations

$$H=-\sum _{ij}{J}_{ij}{S}_i{S}_j-\sum _i{h}_i{S}_i,$$

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

$$\text{reflectbox}R={\log }_{10}({\rho }_s).$$

$$\text{reflectbox}R(M,L_0)\le M\text{reflectbox}R(L_0),$$

$$\text{reflectbox}R(M,L_0)=\text{reflectbox}R(L_0)M^{\alpha },$$

Full text

Many optimization problems belong to the NP-hard complexity class, for which it is believed that no algorithms exist to solve them in polynomial time. Spin-glass problems without biases and on nonplanar topologies, such as the Edward-Anderson (EA) model Edwards and Anderson (1975), represent a sub-class of the NP-hard class. Because spin glasses are the simplest models with both disorder and frustration that fall into the NP-hard class, they represent the ideal model systems to benchmark algorithms, as well as novel computing architectures. A number of heuristics, as well as exhaustive search methods, have been designed and developed to minimize spin-glass Hamiltonians as efficiently as possible. These method also include simulated annealing Kirkpatrick et al. (1983), parallel tempering Monte Carlo Swendsen and Wang (1986); Geyer (1991); Hukushima and Nemoto (1996); Moreno et al. (2003), population annealing Monte Carlo Hukushima and Iba (2003); Zhou and Chen (2010); Wang et al. (2015a, b), and genetic algorithms Pal (1996); Hartmann and Rieger (2001), as well as branch-and-cut De Simone et al. (1995) algorithms, to name a few. Many of these optimization algorithms use only local updates during the minimization procedure. However, in many cases, the use of cluster algorithms with nonlocal updates can greatly enhance the searching process when the energy landscape has many metastable states with small overlap Houdayer (2001); Zhu et al. (2015, 2016a). In the last two decades, quantum heuristics have been proposed as an alternative to classical heuristics, due to their potential to exploit quantum superposition and quantum tunneling effects. Among quantum approaches, adiabatic quantum optimization (AQO) is widely used Kadowaki and Nishimori (1998); Nishimori (2001); Farhi et al. (2001); Santoro et al. (2002); Roland and Cerf (2002); Boixo et al. (2013, 2014); Rønnow et al. (2014); Boixo et al. (2016); Pudenz et al. (2015); Perdomo-Ortiz et al. (2015); Vinci et al. (2015); Mandrà et al. (2015); Venturelli et al. (2015) and likely the method most amenable to hardware implementations Johnson et al. (2011). Current state-of-the-art AQO hardware is manufactured by D-Wave System Inc., whose latest chip allows for the quantum optimization of problems of approximately up to $2000$ variables. However, whether AQO can be more efficient than classical algorithms for certain problems is still controversial Heim et al. (2015); Mandrà et al. (2016, 2017).

Given the importance of comparing optimization techniques across disciplines, it is necessary to have benchmark problems that are (1) representative of the hardness of a typical NP-hard problem, (2) scalable for large systems, and for which (3) the ground state is known a priori. While it is easy to fulfill criteria (1) and (2), it is challenging to have large problems with known solutions.

There have been previous approaches to plant solutions for benchmarking purposes. For example, Ref. Hen et al. (2015) used an approach based on constraint satisfaction problems. Although these problems are tunable in hardness, there is little control when selecting the coupler values between the individual variables. For analog machines with finite precision, such as the D-Wave quantum annealers, this could be an unnecessary restriction. Other approaches Marshall et al. (2016) start from a random coupler configuration and then stochastically update the values of the couplers with a penalty that directly correlates to the time-to-solution of a given solver. However, this approach has two shortcomings: Fist, it assumes that the typical computational hardness com (a) of a problem for a given algorithm will carry over to other optimization techniques. Second, for extremely large problems, the stochastic approach will take sizable resources to thermalize and thus will not be practical.

The method we propose here and which we call “patch planting” (see Fig. 1), where we solve small problems (patches) with open boundaries and then stitch these together to plant an arbitrarily large solution to an instance, does not suffer from these shortcomings: First, arbitrarily large problems can be generated. Second, by assessing the typical complexity using the entropic family size of population annealing Monte Carlo, a metric that characterizes the landscape of the problem and not the algorithmic complexity, we do not depend on the behavior of a particular algorithm when assessing the typical time to solution for a particular instance. Finally, the method poses no restrictions to coupler values, biases (field terms), or lattice topologies.

The paper is structured as follows. In Sec. I we introduce the benchmark problem, as well the patch planting algorithm. In Sec. II we use simulated annealing, population annealing Monte Carlo, as well as experiments on the D-Wave 2X quantum annealer to illustrate how patch planting can produce computational problems that are typically hard. Concluding remarks are presented in Sec. III.

I Patch planting

The patch planting heuristic can be described via the following steps:

• (i)

Find the ground state of patches using free boundary conditions.

• (ii)

For each patch, choose an arbitrary ground-state configuration.

• (iii)

Connect the patches with couplers between the free boundary spins ensuring that all couplings are satisfied.

Note that the patches can be chosen arbitrarily, as long as they can be glued together to form the desired problem or topology with the edges to be stitched together having free boundaries. In addition, the individual patches can be solved with any available optimization technique. As demonstrated below, it is important to solve as large patches as possible, because this will result in problems of comparable computational complexity to purely random problems. In some cases, the breakup of a problem might result in a patch that can be solved exactly, i.e., in polynomial time. Finally, when stitching the patches together, as shown in Fig. 1, it is important to “satisfy” the interaction between two spins of different patches. This means that the coupler has to be chosen as to minimize the dimer’s energy. Knowing the minimizing configuration of the individual patches and assigning the stitching couplers as to satisfy the interactions between spins of neighboring patches then results in a larger planted solution com (b).

As described in Sec. II in more detail, the typical computational complexity of the patched problem can be tuned by either changing the patch size (the larger, the harder) or using hard patches (the harder the patch, the harder the compound problem) e.g., by measuring the entropic family size via population annealing Monte Carlo. This metric can be measured with little numerical effort, and gives a good representation of the typical computational complexity of a problem. Therefore, by post-selecting individual patches, problems of different typical computational complexity can be generated.

Note that in the description of the patch planting procedure no details of the problem to be studied have been mentioned, because the approach is agnostic to the choice of couplers and topologies. We thus emphasize that the patch planting approach can be used for problems of arbitrary topology and for an arbitrary set of coupler values and biases. As such, solutions for arbitrary problems can be planted. This is of much importance when attempting to generate problem sets with particular features, such as synthetic application problems that are known to have a specific nonrandom structure, or problems where the minimum energy gap is fixed (and large) to mitigate the effects of noise on analog optimization machines Katzgraber et al. (2015); Zhu et al. (2016b).

II Experiments

II.1 Benchmark problem

To test the properties of patch planting, we use the Edward-Anderson (EA) Ising spin-glass model Edwards and Anderson (1975), initially in three space dimensions. Later, we perform experiments on the D-Wave 2X quantum annealer using the native topology of the machine Bunyk et al. (2014). The EA spin glass is defined by the following Hamiltonian

[TABLE]

where $S_{i}\in\{\pm 1\}$ are Ising spins and the first sum is over spin-spin interactions. For a three-dimensional lattice, the sum is over nearest neighbors on a cubic lattice. For simplicity, all the local magnetic fields are set to zero, i.e., $h_{i}=0$. We do emphasize, however, that patch planting also works with external biases. The spin-spin interactions $J_{ij}$ are chosen from a normal distribution with zero mean and unit variance. A set of the couplings defines an “instance.”

Given the hardware limitations of the D-Wave quantum chips, instances for the D-Wave 2X have been created by planting and patching together K44 unit cells following the two-dimensional logical structure of the Chimera graph. The couplers are randomly drawn from the Sidon set $\{\pm 5,\pm 6,\pm 7\}$ Zhu et al. (2016b). In both cases, we use free boundary conditions (FBCs) for the patches to plant larger instances. We also compare our patched instances with free boundary conditions to random instances with periodic boundary conditions (PBCs).

II.2 Simulation details

We use the entropic family size of population annealing Monte Carlo $\rho_{s}$ Wang et al. (2015b) to characterize the hardness of the instances. All simulation parameters for the three-dimensional Edwards-Anderson model are listed in Table 1. For the Chimera graph studies on the D-Wave 2X machine, we find the ground state of the patches using $R=2\times 10^{5}$ population members, $N_{T}=301$ temperature steps, $N_{S}=10$ Monte Carlo sweeps, and $T_{0}=0.1$ the lowest temperature of the anneal. The simulation for random problems are done with the same parameters, except $R=10^{6}$.

Experiments on the D-Wave 2X quantum annealer have been performed using a chip with $N=1097$ working qubits. For the Chimera graph, we used all available qubits and patched the system using either two, three, or four patches, respectively. That is, if the system has $12\times 12$ K44 cells of eight qubits each, we divide the whole lattice into two patches of $6\times 12$ K44 cells, three patches of $4\times 12$ K44 cells, or four patches of $3\times 12$ K44 cells. For the experiments, we used an annealing time of $20\ \mu$s, $100$ gauges, and $1000$ readouts for each gauge.

II.3 Correlation between typical hardness and the entropic

family size

The first crucial step in investigating the hardness of instances is to find a good metric that reliably characterizes the typical problem complexity, yet is easy to measure with little computational cost. One approach would be to use the success probability of simulated annealing as a proxy. However, even for medium-size systems, this metric becomes unreliable and computationally expensive. Another possibility consists in using specialized classical algorithms Marshall et al. (2016), such as the Hamze–de Freitas–Selby heuristic Hamze and de Freitas (2004); Selby (2014). However, in this case the typical computational complexity depends on a chosen algorithm and not on the intrinsic properties of the problem’s energy landscape. The latter can be mapped out well for random problems using parallel tempering Monte Carlo Katzgraber et al. (2015), however, at sizable computational cost for large patches. Therefore, in this work we infer the typical hardness of instances through the entropic family size $\rho_{s}$ of population annealing Monte Carlo.

Population annealing (PA) Monte Carlo Hukushima and Iba (2003); Machta (2010) is closely related to simulated annealing (SA), except that it uses a population of $R$ replicas and the population is resampled at each temperature anneal step to maintain thermal equilibrium. At each simulation step, replicas are duplicated accordingly to the ratio between the Boltzmann factors computed after and before updating the temperature. This means that replicas with lower (higher) energy tend to be duplicated (eliminated), ensuring the correct representation of the Boltzmann distribution. Therefore, PA improves the probability to find the lowest energy state over SA by more efficiently sampling phase space. We choose to normalize our replicas so that the population size stays approximately the same. Similar to SA, Metropolis sweeps are applied to each replica at the new temperature. At low temperatures, most of the original population is eliminated in the resampling steps and the final population is a descendant from a small subset of the initial population. Let $n_{i}$ be the fraction of the population from family $i$ in the initial population, then

[TABLE]

Here, $\rho_{s}$ represents the characteristic survival family size. The larger $\rho_{s}$ is, the less surviving families, i.e., the more rugged the energy landscape. Moreover, $\rho_{s}$ correlates strongly with the integrated autocorrelation time of parallel tempering, which is also a proxy towards the roughness of the energy landscape Wang et al. (2015b). Note that $\rho_{s}$ converges quickly in population size and is easily estimated with simulations. See Ref. Wang et al. (2015b) for more details on population annealing. Because $\rho_{s}$ is approximately log-normal distributed (see Fig. 2), let us define the logarithm of $\rho_{s}$ as

[TABLE]

Figure 3 shows the correlation between the probability to find the ground state for SA, $p_{\rm SA}$ at inverse temperature $\beta=1/T=5$ and

R

(data taken from Refs. Wang et al. (2014, 2015a)). As expected, the probability of success decreases by increasing

R

. Indeed, SA struggles more to find the ground state when the energy landscape is more rugged. Therefore,

R

represents a good metric to estimate the typical hardness of optimization problems. In this work, because we study large patched system sizes in three dimensions, we have used

R

at $\beta=3$, which is still a low temperature compared to the spin-glass transition temperature for this model Katzgraber et al. (2006). For the Chimera graph, where there is no phase transition, we have used

R

at a considerably lower temperature $\beta=10$.

II.4 Results in three space dimensions

We first focus on the scaling properties of

R

for patch-planted instances by either varying the patch sizes $L_{0}$ or the system size $L$. In addition, we also demonstrate that harder patches can be used to patch harder instances.

Let $M=(L/L_{0})^{3}$ the number of patches of size $L_{0}$. For random instances, $\rho_{s}$ grows exponentially with $L$ Wang et al. (2015b). Because one would expect that $\rho_{s}$ for a problem of size $L$ by patching $M$ patches of size $L_{0}$ cannot be larger than the product of $\rho_{s}$ of the individual patches, the patched instance complexity is bounded,

[TABLE]

where $\reflectbox{R}(M,L_{0})$ is

R

of the patched instance of $M$ patches of size $L_{0}$ and $\reflectbox{R}(1,L_{0})\equiv\reflectbox{R}(L_{0})$ is

R

of a patch. In Fig. 4, we show the scaling of

R

by varying the number of patches $M$ and a power-law fit of the form

[TABLE]

where $0<\alpha<1$.

R

scales sub-linearly with $M$ with an exponent $\alpha=0.31(3)$. This proves that the patch planted instances become harder by increasing the system size via the number of patches. Therefore, it is guaranteed that, for a sufficiently large number of patches, patch planted instances can become arbitrarily hard in the thermodynamic limit. Figure 5 shows the scaling of the exponent $\alpha$ by increasing the size of the patches $L_{0}$, while keeping the number of patches fixed to $M=8$. As one can see from the figure, $\alpha$ remains roughly constant for a wide range of $L_{0}$ values implying that $\alpha$ is a characteristic constant for patch planted problems. It is interesting to compare the scaling with random instances by defining an effective number of blocks as $M=(L/L_{0})^{3}$, also shown in Fig. 4. We find that both random and patch planted instances have a similar scaling form, although the random class has a larger exponent $\alpha^{\prime}=0.37(4)$, as expected. Therefore, $\rho_{s}$ for patched instances also approximately scales exponentially with system size $L$, as is the case for random instances. Note that $\alpha$ and $\alpha^{\prime}$ likely depend on the characteristics of the problem to be studied.

One may also expect to have some benefit by using either larger or harder patches. Indeed, in both cases, this results in having a larger value of

R

. In Fig. 6 we show the effects of having larger patches by analyzing the distribution of

R

at fixed size of the system, $L=16$, using two different patch sizes, $L_{0}=4$ [Fig. 6(a)] and $L_{0}=8$ [Fig. 6(b)]. As one can see, patched instances are consistently harder by using larger patches for a fixed system size. Similarly, in Fig. 7 we show the distribution of

R

by patching instances with $M=8$ patches of size $L_{0}=6$ by either using easy [Fig. 7(a)] or hard [Fig. 7(b)] patches. We defined easy patches as the $8000$ patches with the smallest

R

and hard patches as the $8000$ patches with the largest

R

from the $41472$ patches randomly generated. From these, $1000$ easy and $1000$ hard instances are then generated.

Patch-planted instances generated using hard patches are consistently harder than patch planted instances assembled from easy patches with the mean value of

R

for both cases being $3.214(5)$ and $2.957(4)$, respectively. We note that the approach pioneered in Ref. Marshall et al. (2016) applied to the production of patches could be combined with patch planting to generate unusually hard planted problems. It is also interesting to compare the complexity of the patched instances with random instances. We show the distribution of

R

for $L=8$ [Figs. 8(a) and 8(b)] and $L=12$ [Figs. 8(c) – 8(e)] with different patch sizes and random instances. One can see that while the patched instances are generally easier than the random instances, they are not necessarily trivial. There is clear overlap between the distributions, i.e., by mining the data one can obtain problems of comparable typical complexity. Note also, that the typical complexity grows with increasing patch size for a fixed system size.

Finally, we comment on the performance of parallel tempering (PT) on patched instances. Because population annealing and parallel tempering have a similar performance in both thermal sampling and optimization, and given that the entropic family size correlates strongly with the integrated autocorrelation time (characteristic measure of hardness of parallel tempering) Wang et al. (2015b), it is natural to expect the proposed patch-planted instances to be hard also for PT. To this end, it is noteworthy to mention recent results that analyze the performance of PT with isoenergetic cluster moves (ICM) (see Ref. Zhu et al. (2016a)) in solving patch-planted instances Karimi et al. (2017). PT combined with ICM has been found to be one of the best classical heuristics in solving hard optimization problems Mandrà et al. (2016). However, Ref. Karimi et al. (2017) clearly show that PT is not able to efficiently solve patch-planted instances (see Figs. 8 – 10).

II.5 Experiments on the D-Wave quantum annealer

We complement the numerical studies on three-dimensional spin glasses by experiments on the D-Wave 2X quantum annealer. For this purpose, we patch plant problems on the native topology of the machine and measure the probabilities to find the ground state $p_{\rm{succ}}$ over multiple runs. In addition, we compare to random problems and show correlation plots between the success probabilities and

R

.

The topology of the machine with $N=1097$ working qubits is cut into two, three, and four patches; see Fig. 9 for a graphical representation. For each experiment, we study $10^{3}$ instances. For the patch planted instances, we first generate $10^{4}$ patch planted problems from random patches and then use

R

to select the $10^{3}$ hardest ones. The distributions of

R

for the problems studied is shown in Fig. 10. One can see that for an increasing number of patches $M$ the problems become computationally easier. However, again by mining the data as done above results in hard problems. Figure 11 shows the sorted success probabilities for the $1000$ problems studied for different number of patches. One can see that problems with $M=4$ are computationally much easier. It remains to be tested if changing the shape of the patch could make the problems harder. For example, the four patches could be chosen to be comprised of $6\times 6$ K44 cells. Finally, Fig. 12 shows a correlation plot between success probabilities and

R

. As can be seen, there is a good correlation between these two quantities, especially for larger patches. Experiments (not shown) suggest that the correlation becomes more pronounced for larger system sizes. With some data mining and only a $10$-fold overhead, instances with two patches $M=2$ have approximately the same complexity as the random ones, which are harder than instances with three patches and four patches. Statistics of the success probabilities are shown in Table 2.

III Summary

We have introduced the concept of patch planting to create planted solutions to Ising-type optimization problems for arbitrarily large systems. The method does not restrict the values of the couplers and works for any topology that can be decomposed into patches. We studied in detail the scaling of the typical complexity of the patched instances and compared this to random instances using population annealing Monte Carlo and the D-Wave 2X machine. From our results it is clear that one should use as large patches as possible to more faithfully reproduce the hardness of random problems. Patch planting is easy to implement and could be used to generate benchmark instances for future generations of quantum devices, as well as classical algorithms and any other novel hardware. The approach is generic in that solutions could also be planted for other paradigmatic optimization problems (e.g., the traveling salesman problem) with only minor modifications.

Acknowledgements.

We thank Firas Hamze, Jon Machta, and Ethan Brown for helpful discussions. H.G.K. would like to thank I. P. Freely and I. M. A. Wiener for inspiration at the early stages of the project. W.W. and H.G.K. acknowledge support from the National Science Foundation (Grant No. DMR-1151387). The research of H.G.K. 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 for access to their Ada and Curie clusters.